#define __STDC_LIMIT_MACROS #include #include #include "rcut.h" #include "Purge_Dups.h" #include "Correct.h" #include "ksort.h" #include "kthread.h" #include "hic.h" #define VERBOSE_CUT 0 #define mc_edge_key(e) ((e).x) KRADIX_SORT_INIT(mce, mc_edge_t, mc_edge_key, member_size(mc_edge_t, x)) #define mb_edge_key(e) ((e).x) KRADIX_SORT_INIT(mbe, mb_edge_t, mb_edge_key, member_size(mb_edge_t, x)) #define mc_generic_key(x) (x) KRADIX_SORT_INIT(mc64, uint64_t, mc_generic_key, 8) KRADIX_SORT_INIT(mc32, uint32_t, mc_generic_key, 4) #define pt_a(x, id) ((x).ma.a + ((x).idx.a[(id)]>>32)) #define pt_n(x, id) ((uint32_t)((x).idx.a[(id)])) #define ma_x(z) (((z).x>>32)) #define ma_y(z) (((uint32_t)((z).x))) #define mcb_pat(x, id) (&((x).m.a[(id)-1])) #define mcp_de(x, id, m) ((x).z[((id)<<(x).hapN)+(m)]) uint8_t bit_filed[8] = {1, 2, 4, 8, 16, 32, 64, 128}; #define is_bit_set(id, a) ((a)[(id)>>3]&bit_filed[(id)&7]) typedef struct { int32_t max_iter; int32_t n_perturb; int32_t n_b_perturb; int32_t n_s_perturb; double f_perturb; uint64_t seed; } mc_opt_t; typedef struct { t_w_t z[2]; } mc_pairsc_t; typedef struct { uint64_t x; // RNG uint32_t cc_off, cc_size; kvec_t(uint64_t) cc_edge; uint32_t *cc_node; uint32_t *bfs, *bfs_mark; mc_pairsc_t *z, *z_opt;///keep scores to nodes(1) and nodes(-1) int8_t *s, *s_opt; uint8_t *f; } mc_svaux_t; typedef struct { uint64_t x; // RNG uint32_t cc_off, cc_size; kvec_t(uint64_t) cc_edge; uint32_t *cc_node; uint32_t *bfs, *bfs_mark; mb_node_t *u, *u_opt;///keep status uint8_t *f; } mb_svaux_t; typedef struct{ uint64_t chain_id, bid, uid; }mc_bp_iter; typedef struct{ uint32_t chain_id; uint32_t f_bid; uint32_t f_uid; uint32_t l_bid; uint32_t l_uid; uint32_t id; t_w_t w; }mc_bp_res; typedef struct{ size_t n, m; uint8_t *a; }bits_p; typedef struct{ bubble_type* b_b; uint64_t *idx, idx_n, occ, n_thread; mc_bp_res *res; uint8_t *lock; mc_svaux_t *b_aux; mc_match_t *ma; bits_p *vis; }mc_bp_t; typedef struct { uint64_t x; // RNG uint32_t cc_off, cc_size; kvec_t(uint32_t) cc_node; kvec_t(uint32_t) bfs; ///share uint32_t *bfs_mark; mc_pairsc_t *z, *z_opt;///keep scores to nodes(1) and nodes(-1) int8_t *s, *s_opt; } mc_svaux_t_s; typedef struct { kvec_t(mc_svaux_t_s) bs; kvec_t(uint64_t) cc_edge; uint32_t *bfs_mark; mc_pairsc_t *z, *z_opt;///keep scores to nodes(1) and nodes(-1) int8_t *s, *s_opt; uint32_t n_t, n_g; } mc_svaux_t_mul; void mc_opt_init(mc_opt_t *opt, int32_t n_perturb, double f_perturb, uint64_t seed) { memset(opt, 0, sizeof(mc_opt_t)); opt->n_perturb = n_perturb; opt->f_perturb = f_perturb; opt->max_iter = 1000; opt->seed = seed; opt->n_s_perturb = n_perturb; opt->n_b_perturb = n_perturb*0.5; } void mc_merge_dup(mc_g_t *mg) // MUST BE sorted { uint32_t i, j, k, st; w_t w; for (st = 0, i = 1, k = 0; i <= mg->e->ma.n; ++i) { if (i == mg->e->ma.n || mg->e->ma.a[i].x != mg->e->ma.a[st].x) { if (i - st > 1) { for (j = st, w = 0; j < i; ++j) { w += mg->e->ma.a[j].w; } mg->e->ma.a[k] = mg->e->ma.a[st]; mg->e->ma.a[k++].w = w; } else mg->e->ma.a[k++] = mg->e->ma.a[st]; st = i; } } mg->e->ma.n = k; } void mb_merge_dup(mb_g_t *mbg) // MUST BE sorted { uint32_t i, j, k, st; t_w_t w[4]; for (st = 0, i = 1, k = 0; i <= mbg->e->ma.n; ++i) { if (i == mbg->e->ma.n || mbg->e->ma.a[i].x != mbg->e->ma.a[st].x) { if (i - st > 1) { w[0] = w[1] = w[2] = w[3] = 0; for (j = st; j < i; ++j) { w[0] += mbg->e->ma.a[j].w[0]; w[1] += mbg->e->ma.a[j].w[1]; w[2] += mbg->e->ma.a[j].w[2]; w[3] += mbg->e->ma.a[j].w[3]; } mbg->e->ma.a[k] = mbg->e->ma.a[st]; mbg->e->ma.a[k].w[0] = w[0]; mbg->e->ma.a[k].w[1] = w[1]; mbg->e->ma.a[k].w[2] = w[2]; mbg->e->ma.a[k].w[3] = w[3]; k++; } else mbg->e->ma.a[k++] = mbg->e->ma.a[st]; st = i; } } mbg->e->ma.n = k; } static void mc_edges_idx(mc_match_t *ma) { uint32_t st, i; kv_resize(uint64_t, ma->idx, ma->n_seq); ma->idx.n = ma->n_seq; memset(ma->idx.a, 0, ma->idx.n*sizeof(uint64_t)); for (st = 0, i = 1; i <= ma->ma.n; ++i) if (i == ma->ma.n || (ma->ma.a[i].x>>32) != (ma->ma.a[st].x>>32)) ma->idx.a[ma->ma.a[st].x>>32] = (uint64_t)st << 32 | (i - st), st = i; } static void mb_edges_idx(mb_match_t *ma) { uint32_t st, i; kv_resize(uint64_t, ma->idx, ma->n_seq); ma->idx.n = ma->n_seq; memset(ma->idx.a, 0, ma->idx.n*sizeof(uint64_t)); for (st = 0, i = 1; i <= ma->ma.n; ++i) if (i == ma->ma.n || (ma->ma.a[i].x>>32) != (ma->ma.a[st].x>>32)) ma->idx.a[ma->ma.a[st].x>>32] = (uint64_t)st << 32 | (i - st), st = i; } mc_g_t *init_mc_g_t(ma_ug_t *ug, asg_t *read_g, int8_t *s, uint32_t renew_s) { mc_g_t *p = NULL; CALLOC(p, 1); p->ug = ug; p->rg = read_g; kv_init(p->s); if(s) { p->s.a = s; p->s.n = ug->g->n_seq; p->s.m = 0; if(renew_s) memset(p->s.a, 0, p->s.n); } else { CALLOC(p->s.a, ug->g->n_seq); p->s.n = p->s.m = ug->g->n_seq; } return p; } static mc_edge_t *get_mc_edge(const mc_match_t *ma, uint32_t sid1, uint32_t sid2) { mc_edge_t *o = pt_a(*ma, sid1); uint32_t n = pt_n(*ma, sid1), k; for (k = 0; k < n; ++k) if (((uint32_t)o[k].x) == sid2) return &(o[k]); return NULL; } static mb_edge_t *get_mb_edge(const mb_match_t *ma, uint32_t sid1, uint32_t sid2) { mb_edge_t *o = pt_a(*ma, sid1); uint32_t n = pt_n(*ma, sid1), k; for (k = 0; k < n; ++k) if (((uint32_t)o[k].x) == sid2) return &(o[k]); return NULL; } uint32_t mb_edges_symm(mb_match_t *ma); uint32_t debug_mb_edges_symm(mb_match_t *ma) { uint32_t i, n = 0; mb_edge_t *t = NULL, *m = NULL; for (i = 0; i < ma->ma.n; ++i) { m = &ma->ma.a[i]; if (ma_x(*m) == ma_y(*m)) { fprintf(stderr, "ERROR-0-::%s\n", __func__); continue; } t = get_mb_edge(ma, ma_y(*m), ma_x(*m)); if(!t) { fprintf(stderr, "ERROR-1-::%s\n", __func__); continue; } if(m->w[0] != t->w[0]) { fprintf(stderr, "\nERROR-2-::%s\n", __func__); fprintf(stderr, "ma_x(*m): %lu, ma_y(*m): %u\n", ma_x(*m), ma_y(*m)); // fprintf(stderr, "m->w[0]: %ld, m->w[1]: %ld, m->w[2]: %ld, m->w[3]: %ld\n", // m->w[0], m->w[1], m->w[2], m->w[3]); // fprintf(stderr, "t->w[0]: %ld, t->w[1]: %ld, t->w[2]: %ld, t->w[3]: %ld\n", // t->w[0], t->w[1], t->w[2], t->w[3]); fprintf(stderr, "m->w[0]: %f, m->w[1]: %f, m->w[2]: %f, m->w[3]: %f\n", m->w[0], m->w[1], m->w[2], m->w[3]); fprintf(stderr, "t->w[0]: %f, t->w[1]: %f, t->w[2]: %f, t->w[3]: %f\n", t->w[0], t->w[1], t->w[2], t->w[3]); } if(m->w[3] != t->w[3]) { fprintf(stderr, "\nERROR-3-::%s\n", __func__); fprintf(stderr, "ma_x(*m): %lu, ma_y(*m): %u\n", ma_x(*m), ma_y(*m)); // fprintf(stderr, "m->w[0]: %ld, m->w[1]: %ld, m->w[2]: %ld, m->w[3]: %ld\n", // m->w[0], m->w[1], m->w[2], m->w[3]); // fprintf(stderr, "t->w[0]: %ld, t->w[1]: %ld, t->w[2]: %ld, t->w[3]: %ld\n", // t->w[0], t->w[1], t->w[2], t->w[3]); fprintf(stderr, "m->w[0]: %f, m->w[1]: %f, m->w[2]: %f, m->w[3]: %f\n", m->w[0], m->w[1], m->w[2], m->w[3]); fprintf(stderr, "t->w[0]: %f, t->w[1]: %f, t->w[2]: %f, t->w[3]: %f\n", t->w[0], t->w[1], t->w[2], t->w[3]); } if(m->w[1] != t->w[2]) { fprintf(stderr, "\nERROR-4-::%s\n", __func__); fprintf(stderr, "ma_x(*m): %lu, ma_y(*m): %u\n", ma_x(*m), ma_y(*m)); // fprintf(stderr, "m->w[0]: %ld, m->w[1]: %ld, m->w[2]: %ld, m->w[3]: %ld\n", // m->w[0], m->w[1], m->w[2], m->w[3]); // fprintf(stderr, "t->w[0]: %ld, t->w[1]: %ld, t->w[2]: %ld, t->w[3]: %ld\n", // t->w[0], t->w[1], t->w[2], t->w[3]); fprintf(stderr, "m->w[0]: %f, m->w[1]: %f, m->w[2]: %f, m->w[3]: %f\n", m->w[0], m->w[1], m->w[2], m->w[3]); fprintf(stderr, "t->w[0]: %f, t->w[1]: %f, t->w[2]: %f, t->w[3]: %f\n", t->w[0], t->w[1], t->w[2], t->w[3]); } if(m->w[2] != t->w[1]) { fprintf(stderr, "\nERROR-5-::%s\n", __func__); fprintf(stderr, "ma_x(*m): %lu, ma_y(*m): %u\n", ma_x(*m), ma_y(*m)); // fprintf(stderr, "m->w[0]: %ld, m->w[1]: %ld, m->w[2]: %ld, m->w[3]: %ld\n", // m->w[0], m->w[1], m->w[2], m->w[3]); // fprintf(stderr, "t->w[0]: %ld, t->w[1]: %ld, t->w[2]: %ld, t->w[3]: %ld\n", // t->w[0], t->w[1], t->w[2], t->w[3]); fprintf(stderr, "m->w[0]: %f, m->w[1]: %f, m->w[2]: %f, m->w[3]: %f\n", m->w[0], m->w[1], m->w[2], m->w[3]); fprintf(stderr, "t->w[0]: %f, t->w[1]: %f, t->w[2]: %f, t->w[3]: %f\n", t->w[0], t->w[1], t->w[2], t->w[3]); } } return n; } inline void decode_mb_node(mb_g_t *mbg, uint32_t id, uint32_t **a0, uint32_t *n0, mc_node_t *s0, uint32_t **a1, uint32_t *n1, mc_node_t *s1) { if(a0) (*a0) = mbg->u->bid.a + mbg->u->u.a[id].a[0]; if(n0) (*n0) = mbg->u->u.a[id].occ[0]; if(s0) (*s0) = mbg->u->u.a[id].s[0]; if(a1) (*a1) = mbg->u->bid.a + mbg->u->u.a[id].a[1]; if(n1) (*n1) = mbg->u->u.a[id].occ[1]; if(s1) (*s1) = mbg->u->u.a[id].s[1]; } uint64_t *mb_nodes_core(kv_u_trans_t *ref, mc_match_t* ma, uint32_t occ, uint32_t *fg, kvec_t_u32_warp *sk) { uint32_t i, x, y; uint64_t *group; u_trans_t *o = NULL; for (i = 0; i < occ; ++i) fg[i] = (uint32_t)-1; // connected componets for (i = 0; i < occ; ++i) { if (fg[i] != (uint32_t)-1) continue; if (pt_n(*ma, i) == 0) { fg[i] = i;///group id continue; } sk->a.n = 0; kv_push(uint32_t, sk->a, i); while (sk->a.n > 0) { uint32_t k, j, st, n, t; sk->a.n--; k = sk->a.a[sk->a.n]; fg[k] = i;///group id o = u_trans_a(*ref, k); n = u_trans_n(*ref, k); for (st = 0, j = 1; j <= n; ++j) { if(j == n || o[j].tn != o[st].tn) { t = o[st].tn; if ((pt_n(*ma, t) != 0) && (fg[t] == (uint32_t)-1)) { kv_push(uint32_t, sk->a, t); } st = j; } } } } // precalculate the size of each group CALLOC(group, occ); for (i = 0; i < occ; ++i) group[i] = (uint64_t)fg[i] << 32 | i; radix_sort_mc64(group, group + occ); for (i = 1, x = y = 0; i <= occ; ++i) { if (i == occ || group[i]>>32 != group[x]>>32) { uint32_t j; for (j = x; j < i; ++j) group[j] = (uint64_t)y << 32 | (uint32_t)group[j];///(group id)|first element in this group ++y, x = i; } } return group; } void assgin_mb_node(mb_nodes_t *x, kv_u_trans_t *ref, mc_match_t* ma, uint32_t *fg, kvec_t_u32_warp *sk, uint64_t *cc, uint32_t cc_off, uint32_t cc_size) { uint32_t i, v, n, st, j, t, pass, uid; u_trans_t *o = NULL; mb_node_t *p = NULL; for (i = 0; i < cc_size; ++i) { v = (uint32_t)cc[cc_off + i];///node id fg[v] = (uint32_t)-1; } sk->a.n = 0; v = (uint32_t)cc[cc_off]; kv_push(uint32_t, sk->a, v); fg[v] = 0; pass = 1; while (sk->a.n > 0) { sk->a.n--; v = sk->a.a[sk->a.n]; if(pt_n(*ma, v) == 0) continue; o = u_trans_a(*ref, v); n = u_trans_n(*ref, v); for (st = 0, j = 1; j <= n; ++j) { if(j == n || o[j].tn != o[st].tn) { t = o[st].tn; if (pt_n(*ma, t) == 0) { st = j; continue; } if (fg[t] == (uint32_t)-1)///uncolor { fg[t] = 1 - fg[v]; kv_push(uint32_t, sk->a, t); } else if(fg[t] == fg[v]) ///color { pass = 0; break; } st = j; } } if(pass == 0) break; } if(pass) { uid = x->u.n; kv_pushp(mb_node_t, x->u, &p); p->z[0] = p->z[1] = p->z[2] = p->z[3] = 0; p->s[0] = p->s[1] = 0; p->occ[0] = p->occ[1] = 0; for (i = 0, p->a[0] = x->bid.n; i < cc_size; ++i) { v = (uint32_t)cc[cc_off + i];///node id if(fg[v] == 0) { kv_push(uint32_t, x->bid, v); p->occ[0]++; x->idx.a[v] = (uid<<1); } } for (i = 0, p->a[1] = x->bid.n; i < cc_size; ++i) { v = (uint32_t)cc[cc_off + i];///node id if(fg[v] == 1) { kv_push(uint32_t, x->bid, v); p->occ[1]++; x->idx.a[v] = (uid<<1) + 1; } } } else { for (i = 0; i < cc_size; ++i) { v = (uint32_t)cc[cc_off + i];///node id uid = x->u.n; kv_pushp(mb_node_t, x->u, &p); p->z[0] = p->z[1] = p->z[2] = p->z[3] = 0; p->s[0] = p->s[1] = 0; p->occ[0] = p->occ[1] = 0; p->a[0] = x->bid.n; kv_push(uint32_t, x->bid, v); p->occ[0]++; x->idx.a[v] = (uid<<1); p->a[1] = x->bid.n; } } } void debug_mb_nodes(mb_nodes_t *x, kv_u_trans_t *ref, mc_match_t* ma) { uint32_t i, k, bid, ori, *a = NULL, a_n, v; for (i = 0; i < x->idx.n; i++) { if(x->idx.a[i] == (uint32_t)-1) { fprintf(stderr, "ERROR-0-::%s\n", __func__); continue; } bid = x->idx.a[i]>>1; ori = x->idx.a[i] & 1; a = x->bid.a + x->u.a[bid].a[ori]; a_n = x->u.a[bid].occ[ori]; for (k = 0; k < a_n; k++) { if(a[k] == i) break; } if(k >= a_n) fprintf(stderr, "ERROR-1-::%s\n", __func__); } uint32_t *tt[2], t_n[2], o_n, m, t, found; int8_t *vis = NULL; CALLOC(vis, x->idx.n); u_trans_t *o = NULL; for (i = 0; i < x->u.n; i++)///each block { tt[0] = x->bid.a + x->u.a[i].a[0]; t_n[0] = x->u.a[i].occ[0]; tt[1] = x->bid.a + x->u.a[i].a[1]; t_n[1] = x->u.a[i].occ[1]; if(t_n[0] == 1 && t_n[1] == 0) continue; for (k = 0; k < t_n[0]; k++) { vis[tt[0][k]] = 1; } for (k = 0; k < t_n[1]; k++) { vis[tt[1][k]] = -1; } for (k = 0; k < t_n[0]; k++) { v = tt[0][k]; o = u_trans_a(*ref, v); o_n = u_trans_n(*ref, v); found = 0; for (m = 0; m < o_n; m++) { t = o[m].tn; if (pt_n(*ma, t) == 0) continue; if (vis[t] != -1) fprintf(stderr, "ERROR-2-::%s\n", __func__); else found = 1; } if(found == 0) fprintf(stderr, "ERROR-2-*::%s\n", __func__); } for (k = 0; k < t_n[1]; k++) { v = tt[1][k]; o = u_trans_a(*ref, v); o_n = u_trans_n(*ref, v); found = 0; for (m = 0; m < o_n; m++) { t = o[m].tn; if (pt_n(*ma, t) == 0) continue; if (vis[t] != 1) fprintf(stderr, "ERROR-3-::%s\n", __func__); else found = 1; } if(found == 0) fprintf(stderr, "ERROR-3-*::%s\n", __func__); } for (k = 0; k < t_n[0]; k++) { vis[tt[0][k]] = 0; } for (k = 0; k < t_n[1]; k++) { vis[tt[1][k]] = 0; } } free(vis); } mb_nodes_t *update_mb_nodes_t(kv_u_trans_t *ref, mc_match_t* ma, uint32_t occ) { uint32_t i, st; uint64_t *cc = NULL; mb_nodes_t *x = NULL; CALLOC(x, 1); x->bid.n = x->u.n = x->idx.n = 0; kv_resize(uint32_t, x->idx, occ); x->idx.n = occ; memset(x->idx.a, -1, sizeof(uint32_t)*x->idx.n); kvec_t_u32_warp stack; kv_init(stack.a); uint32_t *flag = NULL; MALLOC(flag, occ); cc = mb_nodes_core(ref, ma, occ, flag, &stack); for (st = 0, i = 1; i <= occ; ++i) { if (i == occ || cc[st]>>32 != cc[i]>>32) { assgin_mb_node(x, ref, ma, flag, &stack, cc, st, i - st); st = i; } } free(cc); free(flag); kv_destroy(stack.a); /*******************************for debug************************************/ // debug_mb_nodes(x, ref, ma); /*******************************for debug************************************/ return x; } mb_g_t *init_mb_g_t(mc_match_t* e, kv_u_trans_t *ref, uint32_t is_sys) { mc_edge_t *o = NULL; uint32_t i, k, m, n, a_n[2], *a[2], qn, tn, qb, tb; mb_edge_t *t = NULL; mb_g_t *p = NULL; CALLOC(p, 1); p->u = update_mb_nodes_t(ref, e, e->n_seq); p->e = NULL; CALLOC(p->e, 1); p->e->n_seq = p->u->u.n; kv_init(p->e->ma); kv_init(p->e->idx); for (i = 0; i < p->u->u.n; i++)///each block { qb = i; p->u->u.a[i].s[0] = p->u->u.a[i].s[1] = 0; decode_mb_node(p, i, &(a[0]), &(a_n[0]), NULL, &(a[1]), &(a_n[1]), NULL); for (k = 0; k < a_n[0]; k++) { qn = a[0][k]; o = pt_a(*e, qn); n = pt_n(*e, qn); for (m = 0; m < n; m++) { tn = ma_y(o[m]); tb = p->u->idx.a[tn]>>1;///tn is the unitig id; tb is the block id if(tb == qb) { continue; } kv_pushp(mb_edge_t, p->e->ma, &t); t->x = (uint64_t)qb << 32 | tb; t->w[0] = t->w[1] = t->w[2] = t->w[3] = 0; t->w[p->u->idx.a[tn]&1] = o[m].w; } } for (k = 0; k < a_n[1]; k++) { qn = a[1][k]; o = pt_a(*e, qn); n = pt_n(*e, qn); for (m = 0; m < n; m++) { tn = ma_y(o[m]); tb = p->u->idx.a[tn]>>1; if(tb == qb) { continue; } kv_pushp(mb_edge_t, p->e->ma, &t); t->x = (uint64_t)qb << 32 | tb; t->w[0] = t->w[1] = t->w[2] = t->w[3] = 0; t->w[(p->u->idx.a[tn]&1)+2] = o[m].w; } } } radix_sort_mbe(p->e->ma.a, p->e->ma.a + p->e->ma.n); mb_merge_dup(p); // MUST BE sorted mb_edges_idx(p->e); /*******************************for debug************************************/ // debug_mb_edges_symm(p->e); /*******************************for debug************************************/ if(is_sys) mb_edges_symm(p->e); return p; } void destory_mc_g_t(mc_g_t **p) { if(!p || !(*p)) return; if((*p)->s.m == 0) (*p)->s.a = NULL; kv_destroy((*p)->s); if((*p)->e) { kv_destroy((*p)->e->idx); kv_destroy((*p)->e->ma); free((*p)->e->cc); free((*p)->e); } free((*p)); } void destory_mb_g_t(mb_g_t **p) { if(!p || !(*p)) return; kv_destroy((*p)->e->idx); kv_destroy((*p)->e->ma); free((*p)->e->cc); free((*p)->e); kv_destroy((*p)->u->bid); kv_destroy((*p)->u->idx); kv_destroy((*p)->u->u); free((*p)->u); free((*p)); } static void ks_shuffle_uint32_t(size_t n, uint32_t a[], uint64_t *x) { size_t i, j; for (i = n; i > 1; --i) { uint32_t tmp; j = (size_t)(kr_drand_r(x) * i); tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp; } } static void normalize_mc_edge(mc_edge_t *a, mc_edge_t *b) { if(a->w >= b->w) { b->x = (uint32_t)a->x; b->x <<= 32; b->x |= (a->x>>32); b->w = a->w; } else { a->x = (uint32_t)b->x; a->x <<= 32; a->x |= (b->x>>32); a->w = b->w; } } uint32_t mc_edges_symm(mc_match_t *ma) { uint8_t *del = NULL; uint32_t i, k, n = 0; mc_edge_t *t = NULL, *m = NULL; CALLOC(del, ma->ma.n); for (i = 0; i < ma->ma.n; ++i) { m = &ma->ma.a[i]; if (ma_x(*m) == ma_y(*m)) { del[i] = 1, ++n;///self overlap continue; } t = get_mc_edge(ma, ma_y(*m), ma_x(*m)); if(!t) { del[i] = 1, ++n;///self overlap continue; } normalize_mc_edge(m, t); } if (n > 0) { for (i = k = 0; i < ma->ma.n; ++i) if (!del[i]) ma->ma.a[k++] = ma->ma.a[i]; ma->ma.n = k; mc_edges_idx(ma); } free(del); return n; } static void normalize_mb_edge(mb_edge_t *a, mb_edge_t *b) { if(a->w >= b->w) { b->x = (uint32_t)a->x; b->x <<= 32; b->x |= (a->x>>32); b->w[0] = a->w[0]; b->w[3] = a->w[3]; b->w[1] = a->w[2]; b->w[2] = a->w[1]; } else { a->x = (uint32_t)b->x; a->x <<= 32; a->x |= (b->x>>32); a->w[0] = b->w[0]; a->w[3] = b->w[3]; a->w[1] = b->w[2]; a->w[2] = b->w[1]; } } uint32_t mb_edges_symm(mb_match_t *ma) { uint8_t *del = NULL; uint32_t i, k, n = 0; mb_edge_t *t = NULL, *m = NULL; CALLOC(del, ma->ma.n); for (i = 0; i < ma->ma.n; ++i) { m = &ma->ma.a[i]; if (ma_x(*m) == ma_y(*m)) { del[i] = 1, ++n;///self overlap continue; } t = get_mb_edge(ma, ma_y(*m), ma_x(*m)); if(!t) { del[i] = 1, ++n;///self overlap continue; } normalize_mb_edge(m, t); } if (n > 0) { for (i = k = 0; i < ma->ma.n; ++i) if (!del[i]) ma->ma.a[k++] = ma->ma.a[i]; ma->ma.n = k; mb_edges_idx(ma); } free(del); return n; } void debug_mc_interval_t(mc_interval_t *p, uint32_t p_n, uint32_t *p_idx, ma_ug_t *ug, asg_t *rg, trans_chain* t_ch) { fprintf(stderr, "0----------[M::%s]----------\n", __func__); uint32_t i, offset, v, sid, eid, spos, epos, p_status, p_uid, occ; ma_utg_t *u = NULL; mc_interval_t *a = NULL, *t = NULL; for (v = 0; v < ug->u.n; v++) { a = p + p_idx[v]; occ = p_idx[v+1] - p_idx[v]; for (i = 0; i < occ; i++) { if(a[i].uID != v) fprintf(stderr, "ERROR-s\n"); } } for (v = 0, p_status = (uint32_t)-1, p_uid = (uint32_t)-1; v < p_n; v++) { t = &(p[v]); sid = t->nS; eid = t->nE; spos = t->bS; epos = t->bE; if(p_uid == t->uID && p_status == t->hs) { fprintf(stderr, "ERROR-a\n"); } p_status = t->hs; p_uid = t->uID; u = &(ug->u.a[t->uID]); for (i = offset = 0; i < u->n; i++) { if(i == sid) { if(spos != offset) { fprintf(stderr, "ERROR-b\n"); } } if(i == eid) { if(epos != (offset+rg->seq[u->a[i]>>33].len - 1)) { fprintf(stderr, "ERROR-c, real end: %u\n", (uint32_t)(offset+rg->seq[u->a[i]>>33].len - 1)); } } offset += (uint32_t)u->a[i]; if(i >= sid && i <= eid) { if(t_ch->ir_het[u->a[i]>>33] != t->hs) { fprintf(stderr, "ERROR-d: is_r_het: %u, h_status: %u\n", t_ch->ir_het[u->a[i]>>33], t->hs); } } } } } double get_w_scale(kv_u_trans_t *ta) { uint32_t i, max_w_i, min_w_i; double max_w = 0, min_w = 0, w; max_w_i = min_w_i = (uint32_t)-1; for (i = 0; i < ta->n; ++i) { if(ta->a[i].del) continue; if(ta->a[i].nw == 0) continue; w = (ta->a[i].nw >= 0? ta->a[i].nw:-ta->a[i].nw); if(max_w_i == (uint32_t)-1 || max_w < w) max_w = w, max_w_i = i; if(min_w_i == (uint32_t)-1 || min_w > w) min_w = w, min_w_i = i; } if(max_w_i == (uint32_t)-1 || min_w_i == (uint32_t)-1) return 1; if(min_w > 1.1) return 1; double sc_max = (double)(1<<30), sc_min = 1.1; return MIN(sc_max/max_w, sc_min/min_w) + sc_min; } void update_mc_edges(mc_g_t *mg, hap_overlaps_list* ha, kv_u_trans_t *ta, trans_chain* t_ch, double f_rate, uint32_t is_sys) { uint32_t v, i, k, qn, tn, qs, qe, ts, te, occ, as, ae, l, offset, l_pos; uint64_t hetLen, homLen, oLen; mc_interval_t *a = NULL; mc_edge_t *ma = NULL; asg_t* nsg = mg->ug->g; ma_utg_t *u = NULL; mc_interval_t *t = NULL; kvec_t(mc_interval_t) p; kv_init(p); kvec_t(uint32_t) p_idx; kv_init(p_idx); if(t_ch) { kv_push(uint32_t, p_idx, 0); for (v = 0; v < nsg->n_seq; v++) { u = &(mg->ug->u.a[v]); for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k) { if (k == u->n || t_ch->ir_het[u->a[k]>>33] != t_ch->ir_het[u->a[l]>>33]) { kv_pushp(mc_interval_t, p, &t); t->uID = v; t->hs = t_ch->ir_het[u->a[l]>>33]; t->bS = l_pos; t->bE = offset + mg->rg->seq[u->a[k-1]>>33].len - 1; t->nS = l; t->nE = k - 1; l = k; l_pos = offset + (uint32_t)u->a[k-1]; } offset += (uint32_t)u->a[k-1]; } kv_push(uint32_t, p_idx, p.n); } ///debug_mc_interval_t(p.a, p.n, p_idx.a, mg->ug, mg->rg, t_ch); } if(!mg->e) { CALLOC(mg->e, 1); mg->e->n_seq = mg->ug->g->n_seq; kv_init(mg->e->idx); kv_init(mg->e->ma); } if(ha) { for (v = 0; v < ha->num; v++) { for (i = 0; i < ha->x[v].a.n; i++) { if(ha->x[v].a.a[i].score <= 0) continue; if(ha->x[v].a.a[i].xUid == ha->x[v].a.a[i].yUid) continue; if(p.n > 0 && p_idx.n > 0) { /*****************qn*****************/ qn = ha->x[v].a.a[i].xUid; qs = ha->x[v].a.a[i].x_beg_pos; qe = ha->x[v].a.a[i].x_end_pos - 1; a = p.a + p_idx.a[qn]; occ = p_idx.a[qn+1] - p_idx.a[qn]; for (k = 0, hetLen = 0, homLen = 0; k < occ; k++) { as = a[k].bS; ae = a[k].bE; oLen = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0); if(homLen + hetLen > 0 && oLen == 0) break; if(oLen == 0) continue; if(a[k].hs == N_HET) { homLen += oLen; } else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET)) { homLen += oLen; } else { hetLen += oLen; } } if(hetLen <= ((hetLen + homLen)*f_rate)) continue; /*****************qn*****************/ /*****************tn*****************/ tn = ha->x[v].a.a[i].yUid; ts = ha->x[v].a.a[i].y_beg_pos; te = ha->x[v].a.a[i].y_end_pos - 1; a = p.a + p_idx.a[tn]; occ = p_idx.a[tn+1] - p_idx.a[tn]; for (k = 0, hetLen = 0, homLen = 0; k < occ; k++) { as = a[k].bS; ae = a[k].bE; oLen = ((MIN(te, ae) >= MAX(ts, as))? MIN(te, ae) - MAX(ts, as) + 1 : 0); if(homLen + hetLen > 0 && oLen == 0) break; if(oLen == 0) continue; if(a[k].hs == N_HET) { homLen += oLen; } else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET)) { homLen += oLen; } else { hetLen += oLen; } } if(hetLen <= ((hetLen + homLen)*f_rate)) continue; /*****************tn*****************/ } kv_pushp(mc_edge_t, mg->e->ma, &ma); ma->x = (uint64_t)ha->x[v].a.a[i].xUid << 32 | ha->x[v].a.a[i].yUid; ma->w = ha->x[v].a.a[i].score; } } } if(ta) { // double sc = get_w_scale(ta); // fprintf(stderr, "sc: %f\n", sc); for (i = 0; i < ta->n; ++i) { if(ta->a[i].del) continue; if(p.n > 0 && p_idx.n > 0) { /*****************qn*****************/ qn = ta->a[i].qn; qs = ta->a[i].qs; qe = ta->a[i].qe - 1; a = p.a + p_idx.a[qn]; occ = p_idx.a[qn+1] - p_idx.a[qn]; for (k = 0, hetLen = 0, homLen = 0; k < occ; k++) { as = a[k].bS; ae = a[k].bE; oLen = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0); if(homLen + hetLen > 0 && oLen == 0) break; if(oLen == 0) continue; if(a[k].hs == N_HET) { homLen += oLen; } else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET)) { homLen += oLen; } else { hetLen += oLen; } } if(hetLen <= ((hetLen + homLen)*f_rate)) continue; /*****************qn*****************/ /*****************tn*****************/ tn = ta->a[i].tn; ts = ta->a[i].ts; te = ta->a[i].te - 1; a = p.a + p_idx.a[tn]; occ = p_idx.a[tn+1] - p_idx.a[tn]; for (k = 0, hetLen = 0, homLen = 0; k < occ; k++) { as = a[k].bS; ae = a[k].bE; oLen = ((MIN(te, ae) >= MAX(ts, as))? MIN(te, ae) - MAX(ts, as) + 1 : 0); if(homLen + hetLen > 0 && oLen == 0) break; if(oLen == 0) continue; if(a[k].hs == N_HET) { homLen += oLen; } else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET)) { homLen += oLen; } else { hetLen += oLen; } } if(hetLen <= ((hetLen + homLen)*f_rate)) continue; /*****************tn*****************/ } kv_pushp(mc_edge_t, mg->e->ma, &ma); ma->x = (uint64_t)ta->a[i].qn << 32 | ta->a[i].tn; // ma->w = w_cast((ta->a[i].nw*sc)); ma->w = w_cast((ta->a[i].nw)); } } for (i = k = 0; i < mg->e->ma.n; i++) { if(mg->e->ma.a[i].w == 0) continue; mg->e->ma.a[k] = mg->e->ma.a[i]; k++; } mg->e->ma.n = k; radix_sort_mce(mg->e->ma.a, mg->e->ma.a + mg->e->ma.n); mc_merge_dup(mg); mc_edges_idx(mg->e); if(is_sys) mc_edges_symm(mg->e); kv_destroy(p); kv_destroy(p_idx); } void debug_mc_g_t(mc_g_t *mg) { fprintf(stderr, "0----------[M::%s]----------\n", __func__); mc_edge_t *o = NULL, *s = NULL; uint32_t i, k, n, cnt; for (i = 0; i < mg->e->n_seq; ++i) { o = pt_a(*(mg->e), i); n = pt_n(*(mg->e), i); for (k = 0; k < n; ++k) { if(ma_x(o[k]) != i) fprintf(stderr, "ERROR-g\n"); s = get_mc_edge(mg->e, ma_y(o[k]), ma_x(o[k])); if(!s) fprintf(stderr, "ERROR-e\n"); if(s) { if(!(ma_x(*s) == ma_y(o[k]) && ma_y(*s) == ma_x(o[k]) && s->w == o[k].w)) { fprintf(stderr, "ERROR-f\n"); } } } for (k = cnt = 0; k < mg->e->ma.n; ++k) { if(ma_x(mg->e->ma.a[k]) == i) cnt++; } if(cnt != n) fprintf(stderr, "ERROR-h\n"); } } uint64_t *mc_g_cc_core(mc_match_t *ma) { uint32_t i, x, y, *flag; uint64_t *group; mc_edge_t *o = NULL; kvec_t(uint32_t) stack; kv_init(stack); MALLOC(flag, ma->n_seq); for (i = 0; i < ma->n_seq; ++i) flag[i] = (uint32_t)-1; // connected componets for (i = 0; i < ma->n_seq; ++i) { if (flag[i] != (uint32_t)-1) continue; stack.n = 0; kv_push(uint32_t, stack, i); while (stack.n > 0) { uint32_t k, j, n; stack.n--; k = stack.a[stack.n]; flag[k] = i;///group id // n = (uint32_t)ma->idx[k]; // s = ma->idx[k] >> 32; o = pt_a(*ma, k); n = pt_n(*ma, k); for (j = 0; j < n; ++j) { uint32_t t = ma_y(o[j]); if (flag[t] != (uint32_t)-1) continue; // if (ns == ms) PT_EXPAND(stack, ms); // stack[ns++] = t; kv_push(uint32_t, stack, t); } } } kv_destroy(stack); // precalculate the size of each group CALLOC(group, ma->n_seq); for (i = 0; i < ma->n_seq; ++i) group[i] = (uint64_t)flag[i] << 32 | i; radix_sort_mc64(group, group + ma->n_seq); for (i = 1, x = y = 0; i <= ma->n_seq; ++i) { if (i == ma->n_seq || group[i]>>32 != group[x]>>32) { uint32_t j; for (j = x; j < i; ++j) group[j] = (uint64_t)y << 32 | (uint32_t)group[j];///(group id)|first element in this group ++y, x = i; } } free(flag); return group; } void mc_g_cc(mc_match_t *ma) { ma->cc = mc_g_cc_core(ma); } mc_bp_t *mc_bp_t_init(mc_match_t *ma, mc_svaux_t *b_aux, bubble_type* bub, uint64_t n_thread) { uint32_t i, k, n; mc_bp_t *bp = NULL; ma_utg_t *u = NULL; CALLOC(bp, 1); bp->b_aux = b_aux; bp->ma = ma; bp->b_b = bub; bp->n_thread = n_thread; CALLOC(bp->lock, bub->ug->g->n_seq); CALLOC(bp->res, n_thread); CALLOC(bp->vis, n_thread); for (i = 0; i < n_thread; i++) { bp->vis[i].m = bp->vis[i].n = bub->ug->g->n_seq; CALLOC(bp->vis[i].a, bp->vis[i].n); } MALLOC(bp->idx, bub->chain_weight.n+1); bp->idx_n = bp->occ = 0; for (i = 0; i < bub->chain_weight.n; i++) { bp->idx[i] = bp->occ; bp->idx_n++; if(bub->chain_weight.a[i].del) continue; u = &(bub->b_ug->u.a[bub->chain_weight.a[i].id]);///list of bubbles for (k = 0; k < u->n; k++) { get_bubbles(bub, u->a[k]>>33, NULL, NULL, NULL, &n, NULL); bp->occ += n; } } bp->idx[i] = bp->occ; fprintf(stderr, "# nodes in chains: %lu, # chains: %lu\n", bp->occ, bp->idx_n); return bp; } void destroy_mc_bp_t(mc_bp_t **bp) { uint32_t i; for (i = 0; i < (*bp)->n_thread; i++) { free((*bp)->vis[i].a); } free((*bp)->idx); free((*bp)->res); free((*bp)->lock); free((*bp)); } mc_svaux_t *mc_svaux_init(const mc_g_t *mg, uint64_t x) { uint32_t st, i, max_cc = 0; mc_match_t *ma = mg->e; mc_svaux_t *b; CALLOC(b, 1); b->x = x; for (st = 0, i = 1; i <= ma->n_seq; ++i) if (i == ma->n_seq || ma->cc[st]>>32 != ma->cc[i]>>32) max_cc = max_cc > i - st? max_cc : i - st, st = i; kv_init(b->cc_edge); MALLOC(b->cc_node, max_cc); ///CALLOC(b->s, ma->n_seq); b->s = mg->s.a; CALLOC(b->s_opt, ma->n_seq); MALLOC(b->bfs, ma->n_seq); MALLOC(b->bfs_mark, ma->n_seq); memset(b->bfs_mark, -1, ma->n_seq*sizeof(uint32_t)); CALLOC(b->z, ma->n_seq); CALLOC(b->z_opt, ma->n_seq); CALLOC(b->f, ma->n_seq); return b; } void mc_svaux_destroy(mc_svaux_t *b) { b->s = NULL; kv_destroy(b->cc_edge); free(b->cc_node); free(b->s); free(b->s_opt); free(b->z); free(b->z_opt); free(b->bfs); free(b->bfs_mark); free(b->f); free(b); } mb_svaux_t *mb_svaux_init(const mb_g_t *mg, uint64_t x) { uint32_t st, i, max_cc = 0; mb_match_t *ma = mg->e; mb_svaux_t *b; CALLOC(b, 1); b->x = x; for (st = 0, i = 1; i <= ma->n_seq; ++i) if (i == ma->n_seq || ma->cc[st]>>32 != ma->cc[i]>>32) max_cc = max_cc > i - st? max_cc : i - st, st = i; kv_init(b->cc_edge); MALLOC(b->cc_node, max_cc); b->u = mg->u->u.a; CALLOC(b->u_opt, ma->n_seq); MALLOC(b->bfs, ma->n_seq); MALLOC(b->bfs_mark, ma->n_seq); memset(b->bfs_mark, -1, ma->n_seq*sizeof(uint32_t)); CALLOC(b->f, ma->n_seq); return b; } void mb_svaux_destroy(mb_svaux_t *b) { b->u = NULL; kv_destroy(b->cc_edge); free(b->cc_node); free(b->u); free(b->u_opt); free(b->bfs); free(b->bfs_mark); free(b->f); free(b); } mc_svaux_t_mul *init_mc_svaux_t_mul(const mc_g_t *mg, uint64_t n_threads) { uint32_t st, i; mc_match_t *ma = mg->e; mc_svaux_t_mul *b; CALLOC(b, 1); for (st = 0, i = 1, b->n_t = 0; i <= ma->n_seq; ++i) { if (i == ma->n_seq || ma->cc[st]>>32 != ma->cc[i]>>32) { b->n_t++; } } b->n_g = b->n_t; if(n_threads < b->n_t) b->n_t = n_threads; kv_init(b->cc_edge); MALLOC(b->bfs_mark, ma->n_seq); memset(b->bfs_mark, -1, ma->n_seq*sizeof(uint32_t)); CALLOC(b->s_opt, ma->n_seq); CALLOC(b->z, ma->n_seq); CALLOC(b->z_opt, ma->n_seq); b->s = mg->s.a; kv_init(b->bs); CALLOC(b->bs.a, b->n_t); b->bs.n = b->bs.m = b->n_t; for (i = 0; i < b->bs.n; i++) { kv_init(b->bs.a[i].cc_node); kv_init(b->bs.a[i].bfs); b->bs.a[i].bfs_mark = b->bfs_mark; b->bs.a[i].z = b->z; b->bs.a[i].z_opt = b->z_opt; b->bs.a[i].s = b->s; b->bs.a[i].s_opt = b->s_opt; } return b; } void destroy_mc_svaux_t_mul(mc_svaux_t_mul *b) { uint32_t i; kv_destroy(b->cc_edge); free(b->bfs_mark); free(b->s_opt); free(b->z); free(b->z_opt); b->s = NULL; for (i = 0; i < b->bs.n; i++) { kv_destroy(b->bs.a[i].cc_node); kv_destroy(b->bs.a[i].bfs); } kv_destroy(b->bs); free(b); } uint32_t mc_best(const mc_match_t *ma, mc_svaux_t *b) { uint32_t i, max_i = (uint32_t)-1; t_w_t w, max_w; for (i = 0, max_w = -1; i < b->cc_size; ++i) { uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///uid if(b->f[k] || b->s[k] == 0) continue; ///z += -((t_w_t)(b->s[k])) * (b->z[k].z[0] - b->z[k].z[1]); ///-((t_w_t)(b->s[k])) * (b->z[k].z[0] - b->z[k].z[1]) current ///((t_w_t)(b->s[k])) * (b->z[k].z[0] - b->z[k].z[1]) flipped w = ((t_w_t)(b->s[k])) * (b->z[k].z[0] - b->z[k].z[1]) * 2; if(w <= 0) continue; if(w > max_w) { w = max_w; max_i = k; } } return max_i; } t_w_t mb_score(mb_g_t *mbg, mb_svaux_t *b) { uint32_t i; t_w_t z = 0; mb_match_t *ma = mbg->e; for (i = 0; i < b->cc_size; ++i) { uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///uid ///a[0] z += -(t_w_t)(mbg->u->u.a[k].s[0]) * (mbg->u->u.a[k].z[0] - mbg->u->u.a[k].z[1]); ///a[1] z += -(t_w_t)(mbg->u->u.a[k].s[1]) * (mbg->u->u.a[k].z[2] - mbg->u->u.a[k].z[3]); } return z; } t_w_t mc_score(const mc_match_t *ma, mc_svaux_t *b) { uint32_t i; t_w_t z = 0; for (i = 0; i < b->cc_size; ++i) { uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///uid z += -((t_w_t)(b->s[k])) * (b->z[k].z[0] - b->z[k].z[1]); } return z; } t_w_t mc_score_all(const mc_match_t *ma, mc_svaux_t *b) { uint32_t k; t_w_t z = 0; for (k = 0; k < ma->n_seq; ++k) { z += -((t_w_t)(b->s[k])) * (b->z[k].z[0] - b->z[k].z[1]); } return z; } void mc_reset_z(const mc_match_t *ma, mc_svaux_t *b) { uint32_t i; for (i = 0; i < b->cc_size; ++i) { uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///uid uint32_t o = ma->idx.a[k] >> 32; uint32_t j, n = (uint32_t)ma->idx.a[k]; b->z[k].z[0] = b->z[k].z[1] = 0; for (j = 0; j < n; ++j) { const mc_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e); if (b->s[t] > 0) b->z[k].z[0] += e->w; else if (b->s[t] < 0) b->z[k].z[1] += e->w; } } } void mc_reset_z_debug(const mc_match_t *ma, mc_svaux_t *b) { uint32_t i; t_w_t z[2]; for (i = 0; i < b->cc_size; ++i) { uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///uid uint32_t o = ma->idx.a[k] >> 32; uint32_t j, n = (uint32_t)ma->idx.a[k]; z[0] = b->z[k].z[0]; z[1] = b->z[k].z[1]; b->z[k].z[0] = b->z[k].z[1] = 0; for (j = 0; j < n; ++j) { const mc_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e); if (b->s[t] > 0) b->z[k].z[0] += e->w; else if (b->s[t] < 0) b->z[k].z[1] += e->w; } if(z[0] != b->z[k].z[0]) fprintf(stderr, "ERROR1\n"); if(z[1] != b->z[k].z[1]) fprintf(stderr, "ERROR2\n"); } } t_w_t mc_init_spin(const mc_match_t *ma, mc_svaux_t *b) { uint32_t i; b->cc_edge.n = 0; if(b->cc_size <= 2)//mannually clear { for (i = 0; i < b->cc_size; ++i) {///how many nodes b->s[(uint32_t)ma->cc[b->cc_off + i]] = 0; } } for (i = 0; i < b->cc_size; ++i) {///how many nodes uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///node id b->cc_node[i] = k; if(b->s[k] == 0) break; } if(i >= b->cc_size) goto passed; for (i = 0; i < b->cc_size; ++i) {///how many nodes uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///node id uint32_t o = ma->idx.a[k] >> 32;///cc group id uint32_t n = (uint32_t)ma->idx.a[k], j; b->cc_node[i] = k; for (j = 0; j < n; ++j) { w_t w = ma->ma.a[o + j].w; w = w > 0? w : -w; kv_push(uint64_t, b->cc_edge, (uint64_t)((uint32_t)-1 - ((uint32_t)w)) << 32 | (o + j)); } } radix_sort_mc64(b->cc_edge.a, b->cc_edge.a + b->cc_edge.n); for (i = 0; i < b->cc_edge.n; ++i) { // from the strongest edge to the weakest const mc_edge_t *e = &ma->ma.a[(uint32_t)b->cc_edge.a[i]]; uint32_t n1 = ma_x(*e), n2 = ma_y(*e); if (b->s[n1] == 0 && b->s[n2] == 0) { b->x = kr_splitmix64(b->x); b->s[n1] = b->x&1? 1 : -1; b->s[n2] = e->w > 0? -b->s[n1] : b->s[n1]; }/****************************may have bugs********************************/ else if(b->s[n1] == 0) { b->s[n1] = e->w > 0? -b->s[n2] : b->s[n2]; } else if(b->s[n2] == 0) { b->s[n2] = e->w > 0? -b->s[n1] : b->s[n1]; } /****************************may have bugs********************************/ } passed: mc_reset_z(ma, b); return mc_score(ma, b); } ///k is uid static void mc_set_spin(const mc_match_t *ma, mc_svaux_t *b, uint32_t k, int8_t s) { uint32_t o, j, n; int8_t s0 = b->s[k]; if (s0 == s) return; o = ma->idx.a[k] >> 32; n = (uint32_t)ma->idx.a[k]; for (j = 0; j < n; ++j) { const mc_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e);///1->z[0]; (-1)->z[1]; if (s0 != 0) b->z[t].z[(s0 < 0)] -= e->w; if (s != 0) b->z[t].z[(s < 0)] += e->w; } b->s[k] = s; } void debug_mb_z(mb_g_t *mbg, uint32_t k) { t_w_t z[4]; z[0] = z[1] = z[2] = z[3] = 0; mb_match_t *ma = mbg->e; uint32_t o = ma->idx.a[k] >> 32; uint32_t j, n = (uint32_t)ma->idx.a[k]; for (j = 0; j < n; ++j) { const mb_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e); ///a[0]->b[0] if(mbg->u->u.a[t].s[0] > 0) z[0] += e->w[0]; else if(mbg->u->u.a[t].s[0] < 0) z[1] += e->w[0]; ///a[0]->b[1] if(mbg->u->u.a[t].s[1] > 0) z[0] += e->w[1]; else if(mbg->u->u.a[t].s[1] < 0) z[1] += e->w[1]; ///a[1]->b[0] if(mbg->u->u.a[t].s[0] > 0) z[2] += e->w[2]; else if(mbg->u->u.a[t].s[0] < 0) z[3] += e->w[2]; ///a[1]->b[1] if(mbg->u->u.a[t].s[1] > 0) z[2] += e->w[3]; else if(mbg->u->u.a[t].s[1] < 0) z[3] += e->w[3]; } if(mbg->u->u.a[k].z[0] != z[0]) fprintf(stderr, "ERROR-0-::%s, z[0]: %f, n-z[0]: %f\n", __func__, z[0], mbg->u->u.a[k].z[0]); if(mbg->u->u.a[k].z[1] != z[1]) fprintf(stderr, "ERROR-1-::%s, z[1]: %f, n-z[1]: %f\n", __func__, z[1], mbg->u->u.a[k].z[1]); if(mbg->u->u.a[k].z[2] != z[2]) fprintf(stderr, "ERROR-2-::%s, z[2]: %f, n-z[2]: %f\n", __func__, z[2], mbg->u->u.a[k].z[2]); if(mbg->u->u.a[k].z[3] != z[3]) fprintf(stderr, "ERROR-3-::%s, z[3]: %f, n-z[3]: %f\n", __func__, z[3], mbg->u->u.a[k].z[3]); } ///k is uid static void mb_flip_spin(mb_g_t *mbg, mb_svaux_t *b, uint32_t k) { if(mbg->u->u.a[k].s[0] == 0 && mbg->u->u.a[k].s[1] == 0) return; mb_match_t *ma = mbg->e; uint32_t o, j, n; o = ma->idx.a[k] >> 32; n = (uint32_t)ma->idx.a[k]; mbg->u->u.a[k].s[0] = -mbg->u->u.a[k].s[0]; mbg->u->u.a[k].s[1] = -mbg->u->u.a[k].s[1]; for (j = 0; j < n; ++j) { const mb_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e);///1->z[0]; (-1)->z[1]; if(mbg->u->u.a[k].s[0] != 0) { ///note: e is from k to t ///t[0] ---> k[0], so update z[0 + x] and e[0] mbg->u->u.a[t].z[mbg->u->u.a[k].s[0] < 0] += e->w[0]; mbg->u->u.a[t].z[mbg->u->u.a[k].s[0] > 0] -= e->w[0]; ///t[1] ---> k[0], so update z[2 + x] and e[1] mbg->u->u.a[t].z[(mbg->u->u.a[k].s[0] < 0) + 2] += e->w[1]; mbg->u->u.a[t].z[(mbg->u->u.a[k].s[0] > 0) + 2] -= e->w[1]; } if(mbg->u->u.a[k].s[1] != 0) { ///note: e is from k to t ///t[0] ---> k[1], so update z[0 + x] and e[2] mbg->u->u.a[t].z[mbg->u->u.a[k].s[1] < 0] += e->w[2]; mbg->u->u.a[t].z[mbg->u->u.a[k].s[1] > 0] -= e->w[2]; ///t[1] ---> k[1], so update z[2 + x] and e[3] mbg->u->u.a[t].z[(mbg->u->u.a[k].s[1] < 0) + 2] += e->w[3]; mbg->u->u.a[t].z[(mbg->u->u.a[k].s[1] > 0) + 2] -= e->w[3]; } /*******************************for debug************************************/ // debug_mb_z(mbg, t); /*******************************for debug************************************/ } } void mc_best_flip(const mc_match_t *ma, mc_svaux_t *b) { uint32_t idx; for (idx = 0; idx < b->cc_size; ++idx) { b->f[(uint32_t)ma->cc[b->cc_off + idx]] = 0; ///uint32_t k = (uint32_t)ma->cc[b->cc_off + idx];///uid } while (1) { idx = mc_best(ma, b); if(idx == (uint32_t)-1) break; mc_set_spin(ma, b, idx, -b->s[idx]); b->f[idx] = 1; } } static t_w_t mc_optimize_local(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_t *b, uint32_t *n_iter) { uint32_t i, n_flip = 0; int32_t n_iter_local = 0; while (n_iter_local < opt->max_iter) { ++(*n_iter); ks_shuffle_uint32_t(b->cc_size, b->cc_node, &b->x); for (i = n_flip = 0; i < b->cc_size; ++i) { uint32_t k = b->cc_node[i];///uid int8_t s; if (b->z[k].z[0] == b->z[k].z[1]) continue; s = b->z[k].z[0] > b->z[k].z[1]? -1 : 1; if (b->s[k] != s) { mc_set_spin(ma, b, k, s);///no need to change the score of k itself ++n_flip; } } ++n_iter_local; if (n_flip == 0) break; } // if(n_flip != 0) mc_best_flip(ma, b); return mc_score(ma, b); } static t_w_t mb_optimize_local(const mc_opt_t *opt, mb_g_t *mbg, mb_svaux_t *b, uint32_t *n_iter) { uint32_t i, n_flip = 0; int32_t n_iter_local = 0; mb_node_t *u = NULL; t_w_t z = 0; while (n_iter_local < opt->max_iter) { ++(*n_iter); ks_shuffle_uint32_t(b->cc_size, b->cc_node, &b->x); for (i = n_flip = 0; i < b->cc_size; ++i) { uint32_t k = b->cc_node[i];///uid u = &(mbg->u->u.a[k]); if(u->z[0] == u->z[1] && u->z[2] == u->z[3]) continue; z = 0; ///a[0] z += -(t_w_t)(u->s[0]) * (u->z[0] - u->z[1]); ///a[1] z += -(t_w_t)(u->s[1]) * (u->z[2] - u->z[3]); if(z >= 0) continue; mb_flip_spin(mbg, b, k);///no need to change the score of k itself ++n_flip; } ++n_iter_local; if (n_flip == 0) break; } // if(n_flip != 0) mc_best_flip(ma, b); return mb_score(mbg, b); } static void mc_perturb(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_t *b) { uint32_t i; for (i = 0; i < b->cc_size; ++i) { uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///node id double y; y = kr_drand_r(&b->x); if (y < opt->f_perturb) mc_set_spin(ma, b, k, -b->s[k]); } } static void mb_perturb(const mc_opt_t *opt, mb_g_t *mbg, mb_svaux_t *b) { uint32_t i; for (i = 0; i < b->cc_size; ++i) { uint32_t k = (uint32_t)mbg->e->cc[b->cc_off + i];///node id double y; y = kr_drand_r(&b->x); if (y < opt->f_perturb) mb_flip_spin(mbg, b, k); } } static uint32_t mc_bfs(const mc_match_t *ma, mc_svaux_t *b, uint32_t k0, uint32_t bfs_round, uint32_t max_size) { uint32_t i, n_bfs = 0, st, en, r; b->bfs[n_bfs++] = k0, b->bfs_mark[k0] = k0; st = 0, en = n_bfs; for (r = 0; r < bfs_round; ++r) { for (i = st; i < en; ++i) { uint32_t k = b->bfs[i]; uint32_t o = ma->idx.a[k] >> 32; uint32_t n = (uint32_t)ma->idx.a[k], j; for (j = 0; j < n; ++j) { uint32_t t = (uint32_t)ma->ma.a[o + j].x; if (b->bfs_mark[t] != k0) b->bfs[n_bfs++] = t, b->bfs_mark[t] = k0; } } st = en, en = n_bfs; if (max_size > 0 && n_bfs > max_size) break; } return n_bfs; } ///bfs_round is 3 static void mc_perturb_node(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_t *b, int32_t bfs_round) { uint32_t i, k, n_bfs = 0; k = (uint32_t)(kr_drand_r(&b->x) * b->cc_size + .499); if(k >= b->cc_size) k = b->cc_size - 1; k = (uint32_t)ma->cc[b->cc_off + k];///node id n_bfs = mc_bfs(ma, b, k, bfs_round, (int32_t)(b->cc_size * opt->f_perturb)); for (i = 0; i < n_bfs; ++i) mc_set_spin(ma, b, b->bfs[i], -b->s[b->bfs[i]]); } static uint32_t mb_bfs(const mb_match_t *ma, mb_svaux_t *b, uint32_t k0, uint32_t bfs_round, uint32_t max_size) { uint32_t i, n_bfs = 0, st, en, r; b->bfs[n_bfs++] = k0, b->bfs_mark[k0] = k0; st = 0, en = n_bfs; for (r = 0; r < bfs_round; ++r) { for (i = st; i < en; ++i) { uint32_t k = b->bfs[i]; uint32_t o = ma->idx.a[k] >> 32; uint32_t n = (uint32_t)ma->idx.a[k], j; for (j = 0; j < n; ++j) { uint32_t t = (uint32_t)ma->ma.a[o + j].x; if (b->bfs_mark[t] != k0) b->bfs[n_bfs++] = t, b->bfs_mark[t] = k0; } } st = en, en = n_bfs; if (max_size > 0 && n_bfs > max_size) break; } return n_bfs; } ///bfs_round is 3 static void mb_perturb_node(const mc_opt_t *opt, mb_g_t *mbg, mb_svaux_t *b, int32_t bfs_round) { uint32_t i, k, n_bfs = 0; k = (uint32_t)(kr_drand_r(&b->x) * b->cc_size + .499); if(k >= b->cc_size) k = b->cc_size - 1; k = (uint32_t)mbg->e->cc[b->cc_off + k];///node id n_bfs = mb_bfs(mbg->e, b, k, bfs_round, (int32_t)(b->cc_size * opt->f_perturb)); for (i = 0; i < n_bfs; ++i) mb_flip_spin(mbg, b, b->bfs[i]); } void clean_mc_bp_res(mc_bp_res *res) { res->chain_id = (uint32_t)-1; res->f_bid = res->f_uid = res->l_bid = res->l_uid = (uint32_t)-1; res->id = (uint32_t)-1; res->w = -1; } void reset_mc_bp_iter(mc_bp_t* bp, mc_bp_iter *x, uint64_t id) { ma_utg_t *u = NULL; uint32_t i, n, occ; for (i = 0; i < bp->idx_n; i++) { if(id >= bp->idx[i]) break; } id -= bp->idx[i]; x->chain_id = bp->b_b->chain_weight.a[i].id; u = &(bp->b_b->b_ug->u.a[x->chain_id]); for (i = occ = 0; i < u->n; i++) { get_bubbles(bp->b_b, u->a[i]>>33, NULL, NULL, NULL, &n, NULL); occ += n; if(id < occ) { x->bid = i; x->uid = id - (occ -n); break; } } } inline uint32_t next_uid(mc_bp_iter *iter, bubble_type* bub, uint32_t *c_bid, uint32_t *c_uid) { ma_utg_t *u = &(bub->b_ug->u.a[iter->chain_id]); uint32_t *a, n, uid; while (1) { if(iter->bid >= u->n) break; get_bubbles(bub, u->a[iter->bid]>>33, NULL, NULL, &a, &n, NULL); while (1) { if(iter->uid >= n) break; uid = a[iter->uid]>>1; if(c_bid) (*c_bid) = iter->bid; if(c_uid) (*c_uid) = iter->uid; iter->uid++; return uid; } iter->bid++, iter->uid = 0; } return (uint32_t)-1; } t_w_t incre_weight(mc_svaux_t *b_aux, mc_match_t *ma, uint8_t* vis, uint32_t uid) { mc_edge_t *o = NULL; uint32_t n, i, t; t_w_t w = ((t_w_t)(b_aux->s[uid])) * (b_aux->z[uid].z[0] - b_aux->z[uid].z[1]) * 2; t_w_t w_off = 0; o = pt_a(*ma, uid); n = pt_n(*ma, uid); for (i = 0; i < n; ++i) { t = ma_y(o[i]); if(vis[t] == 0) continue; if(t == uid) continue; w_off += (b_aux->s[uid]*b_aux->s[t]*o[i].w); } return w - (w_off*4);//2 for self; 4 for both directions } void select_min_bp(bubble_type* bub, mc_match_t *ma, mc_svaux_t *b_aux, mc_bp_t* bp, uint8_t *lock, bits_p *vis, uint64_t id, mc_bp_res* r) { uint32_t uid, val = 0, max_bid, max_uid, c_bid, c_uid, f_bid, f_uid; t_w_t w = 0, max_w = -1; mc_bp_iter i; reset_mc_bp_iter(bp, &i, id); memset(vis->a, 0, vis->n); max_bid = max_uid = (uint32_t)-1; f_bid = i.bid; f_uid = i.uid;///f_bid::bubble id, f_uid::unitig id while (1) { uid = next_uid(&i, bub, &c_bid, &c_uid); if(uid == (uint32_t)-1) break; if(vis->a[uid]) continue; ///already flip uid ///update w w += incre_weight(b_aux, ma, vis->a, uid); vis->a[uid] = 1; if(lock[uid] == 0) val = 1; if(val == 0) continue; ///update max_w if(max_w < w) { max_w = w; max_bid = c_bid; max_uid = c_uid; } } if(max_w <= 0 || max_bid == (uint32_t)-1 || max_uid == (uint32_t)-1) return; if((max_w > r->w) || (max_w == r->w && id < r->id)) { r->w = max_w; r->id = id; r->chain_id = i.chain_id; r->l_bid = max_bid; r->l_uid = max_uid; r->f_bid = f_bid; r->f_uid = f_uid; } ///if((res->min_w > i_b->weight) || (res->min_w == i_b->weight && id < res->min_idx)) } static void worker_for_min_bp(void *data, long i, int tid) // callback for kt_for() { mc_bp_t* bp = (mc_bp_t *)data; select_min_bp(bp->b_b, bp->ma, bp->b_aux, bp, bp->lock, &(bp->vis[tid]), i, &(bp->res[tid])); } uint32_t best_bp(mc_bp_t *bp, mc_bp_res *res) { uint32_t i; clean_mc_bp_res(res); for (i = 0; i < bp->n_thread; i++) { clean_mc_bp_res(&(bp->res[i])); } kt_for(bp->n_thread, worker_for_min_bp, bp, bp->occ); for (i = 0; i < bp->n_thread; i++) { if(bp->res[i].chain_id == (uint32_t)-1) continue; if(bp->res[i].w <= 0) continue; if((bp->res[i].w > res->w) || (bp->res[i].w == res->w && bp->res[i].id < res->id)) { (*res) = bp->res[i]; } } if(res->chain_id != (uint32_t)-1) return 1; return 0; } void mc_set_bp_spin(bubble_type* bub, mc_match_t *ma, mc_svaux_t *b_aux, mc_bp_t* bp, uint8_t *lock, bits_p *vis, mc_bp_res *res) { uint32_t uid, val = 0, c_bid, c_uid; mc_bp_iter i; i.chain_id = res->chain_id; i.bid = res->f_bid; i.uid = res->f_uid; memset(vis->a, 0, vis->n); while (1) { uid = next_uid(&i, bub, &c_bid, &c_uid); if(uid == (uint32_t)-1) break; if(lock[uid] == 0) { val = 1; break; } if(c_bid == res->l_bid && c_uid == res->l_uid) break; } if(val == 0) { fprintf(stderr, "ERROR-1\n"); return; } i.bid = res->f_bid; i.uid = res->f_uid; while (1) { uid = next_uid(&i, bub, &c_bid, &c_uid); if(uid == (uint32_t)-1) break; if(vis->a[uid] == 1) continue; lock[uid] = 1; vis->a[uid] = 1; mc_set_spin(ma, b_aux, uid, -b_aux->s[uid]); if(c_bid == res->l_bid && c_uid == res->l_uid) break; } } t_w_t mc_solve_bp_cc(mc_bp_t *bp) { mc_bp_res res; memset(bp->lock, 0, bp->b_b->ug->g->n_seq); while (best_bp(bp, &res)) { mc_set_bp_spin(bp->b_b, bp->ma, bp->b_aux, bp, bp->lock, &(bp->vis[0]), &res); } return mc_score_all(bp->ma, bp->b_aux); } void mc_reset_z_all(const mc_match_t *ma, mc_svaux_t *b) { // t_w_t z[2]; uint32_t k; for (k = 0; k < ma->n_seq; ++k) { uint32_t o = ma->idx.a[k] >> 32; uint32_t j, n = (uint32_t)ma->idx.a[k]; // z[0] = b->z[k].z[0]; z[1] = b->z[k].z[1]; b->z[k].z[0] = b->z[k].z[1] = 0; for (j = 0; j < n; ++j) { const mc_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e); if (b->s[t] > 0) b->z[k].z[0] += e->w; else if (b->s[t] < 0) b->z[k].z[1] += e->w; } // if(z[0] != b->z[k].z[0]) fprintf(stderr, "ERROR1-all\n"); // if(z[1] != b->z[k].z[1]) fprintf(stderr, "ERROR2-all\n"); } } void mc_solve_bp(mc_bp_t *bp) { double index_time = yak_realtime(); uint32_t r = 1; t_w_t sc_opt, sc; mc_reset_z_all(bp->ma, bp->b_aux); sc_opt = mc_score_all(bp->ma, bp->b_aux); while (1) { sc = mc_solve_bp_cc(bp); // fprintf(stderr, "[M::%s::# round: %u] sc_opt: %ld, sc: %ld\n", __func__, r, sc_opt, sc); fprintf(stderr, "[M::%s::# round: %u] sc_opt: %f, sc: %f\n", __func__, r, sc_opt, sc); if(sc <= sc_opt) break; sc_opt = sc; r++; } fprintf(stderr, "[M::%s::%.3f] ==> round %u\n", __func__, yak_realtime()-index_time, r); } t_w_t mc_score_all_advance(const mc_match_t *ma, int8_t *s) { uint32_t k; t_w_t z[2], zt = 0; for (k = 0; k < ma->n_seq; ++k) { uint32_t o = ma->idx.a[k] >> 32; uint32_t j, n = (uint32_t)ma->idx.a[k]; z[0] = z[1] = 0; for (j = 0; j < n; ++j) { const mc_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e); if (s[t] > 0) z[0] += e->w; else if (s[t] < 0) z[1] += e->w; } zt += -((t_w_t)(s[k])) * (z[0] - z[1]); } return zt; } void mc_status_all(const mc_match_t *ma, int8_t *s) { t_w_t tt_w = mc_score_all_advance(ma, s); uint64_t k, nn = 0, ne = 0; for (k = 0; k < ma->n_seq; ++k) { if((uint32_t)ma->idx.a[k]) nn++; ne += (uint32_t)ma->idx.a[k]; } fprintf(stderr, "[M::%s::] ==> # nodes: %lu, # edges: %lu, weight: %f\n", __func__, nn, ne, tt_w); } t_w_t mb_score_all_advance(const mc_match_t *ma, mb_g_t *mbg) { uint32_t k; t_w_t z[2], zt = 0; for (k = 0; k < ma->n_seq; ++k) { uint32_t o = ma->idx.a[k] >> 32; uint32_t j, n = (uint32_t)ma->idx.a[k]; z[0] = z[1] = 0; for (j = 0; j < n; ++j) { const mc_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e); if (mbg->u->u.a[mbg->u->idx.a[t]>>1].s[mbg->u->idx.a[t]&1] > 0) z[0] += e->w; else if (mbg->u->u.a[mbg->u->idx.a[t]>>1].s[mbg->u->idx.a[t]&1] < 0) z[1] += e->w; } zt += -((t_w_t)(mbg->u->u.a[mbg->u->idx.a[k]>>1].s[mbg->u->idx.a[k]&1])) * (z[0] - z[1]); } return zt; } void print_sc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, t_w_t sc_opt, uint32_t n_iter) { t_w_t w = mc_score(mg->e, b); // if(w != sc_opt) fprintf(stderr, "ERROR\n"); fprintf(stderr, "# iter: %u, sc_opt: %f, sc-local: %f, sc-global: %f\n", n_iter, sc_opt, w, mc_score_all_advance(mg->e, mg->s.a)); } void print_mc_node(const mc_match_t *ma, mc_svaux_t *b, uint32_t id) { fprintf(stderr, "[M::%s::utg%.6ul-hap%u]\n", __func__, id, b->s[id]>0?1:b->s[id]<0?2:0); w_t w[128], z[4]; uint32_t o, n, i, hn = 4; int8_t s; for (i = 0; i < hn; i++) w[i] = 0; o = ma->idx.a[id] >> 32; n = (uint32_t)ma->idx.a[id]; for (i = 0; i < n; ++i) { s = b->s[ma_y(ma->ma.a[o + i])]; w[s>0?1:s<0?2:0] += ma->ma.a[o + i].w; } z[0] = z[3] = 0; z[1] = b->z[id].z[0]; z[2] = b->z[id].z[1]; for (i = 0; i < hn; i++) fprintf(stderr, "w[%u]-%f, z[%u]-%f\n", i, w[i], i, z[i]); } uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint32_t cc_off, uint32_t cc_size) { uint32_t j, k, n_iter = 0, flush = opt->max_iter * 50; t_w_t sc_opt = -(1<<30), sc;///problem-w b->cc_off = cc_off, b->cc_size = cc_size; if (b->cc_size < 2) return 0; sc_opt = mc_init_spin(mg->e, b); // print_sc(opt, mg, b, sc_opt, n_iter); if (b->cc_size == 2) return 0; for (j = 0; j < b->cc_size; ++j) {///backup s and z in s_opt and z_opt b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; ///hap status of each unitig b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; ///z[0]: positive weight; z[1]: positive weight } // print_sc(opt, mg->e, b, sc_opt, n_iter); sc = mc_optimize_local(opt, mg->e, b, &n_iter); if (sc > sc_opt) { for (j = 0; j < b->cc_size; ++j) { b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; } sc_opt = sc; } else { for (j = 0; j < b->cc_size; ++j) { b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; } } // print_mc_node(mg->e, b, 36880); // print_sc(opt, mg, b, sc_opt, n_iter); // mc_reset_z_debug(mg->e, b); // print_sc(opt, mg->e, b, sc_opt, n_iter); // fprintf(stderr, "\ncc_size: %u, cc_off: %u\n", b->cc_size, b->cc_off); for (k = 0; k < (uint32_t)opt->n_perturb; ++k) { if (k&1) mc_perturb(opt, mg->e, b); else mc_perturb_node(opt, mg->e, b, 3); sc = mc_optimize_local(opt, mg->e, b, &n_iter); // fprintf(stderr, "(%u) sc_opt: %f, sc: %f\n", k, sc_opt, sc); if (sc > sc_opt) { for (j = 0; j < b->cc_size; ++j) { b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; } sc_opt = sc; // print_sc(opt, mg, b, sc_opt, n_iter); } else { for (j = 0; j < b->cc_size; ++j) { b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; } } if((n_iter%flush) == 0) { mc_reset_z(mg->e, b); sc = mc_score(mg->e, b); for (j = 0; j < b->cc_size; ++j) { b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; } sc_opt = sc; } } for (j = 0; j < b->cc_size; ++j) { b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; } return n_iter; } void reset_mb_g_t_z(mb_g_t *mbg); uint32_t mb_solve_cc(const mc_opt_t *opt, mb_g_t *mbg, mb_svaux_t *b, uint32_t cc_off, uint32_t cc_size) { uint32_t j, k, n_iter = 0, flush = opt->max_iter * 50; t_w_t sc_opt = -(1<<30), sc;///problem-w b->cc_off = cc_off, b->cc_size = cc_size; if (b->cc_size <= 2) return 0; for (j = 0; j < b->cc_size; ++j) {///how many nodes k = (uint32_t)mbg->e->cc[b->cc_off + j];///node id b->cc_node[j] = k; } sc_opt = mb_score(mbg, b); for (j = 0; j < b->cc_size; ++j) {///backup s and z in s_opt and z_opt memcpy(b->u_opt[b->cc_node[j]].s, b->u[b->cc_node[j]].s, 2*sizeof(mc_node_t)); memcpy(b->u_opt[b->cc_node[j]].z, b->u[b->cc_node[j]].z, 4*sizeof(t_w_t)); } // print_sc(opt, mg->e, b, sc_opt, n_iter); sc = mb_optimize_local(opt, mbg, b, &n_iter); if (sc > sc_opt) { for (j = 0; j < b->cc_size; ++j) { memcpy(b->u_opt[b->cc_node[j]].s, b->u[b->cc_node[j]].s, 2*sizeof(mc_node_t)); memcpy(b->u_opt[b->cc_node[j]].z, b->u[b->cc_node[j]].z, 4*sizeof(t_w_t)); } sc_opt = sc; } else { for (j = 0; j < b->cc_size; ++j) { memcpy(b->u[b->cc_node[j]].s, b->u_opt[b->cc_node[j]].s, 2*sizeof(mc_node_t)); memcpy(b->u[b->cc_node[j]].z, b->u_opt[b->cc_node[j]].z, 4*sizeof(t_w_t)); } } // mc_reset_z_debug(mg->e, b); // print_sc(opt, mg->e, b, sc_opt, n_iter); // fprintf(stderr, "\ncc_size: %u, cc_off: %u\n", b->cc_size, b->cc_off); for (k = 0; k < (uint32_t)opt->n_perturb; ++k) { if (k&1) mb_perturb(opt, mbg, b); else mb_perturb_node(opt, mbg, b, 3); sc = mb_optimize_local(opt, mbg, b, &n_iter); // fprintf(stderr, "(%u) sc_opt: %f, sc: %f\n", k, sc_opt, sc); if (sc > sc_opt) { for (j = 0; j < b->cc_size; ++j) { memcpy(b->u_opt[b->cc_node[j]].s, b->u[b->cc_node[j]].s, 2*sizeof(mc_node_t)); memcpy(b->u_opt[b->cc_node[j]].z, b->u[b->cc_node[j]].z, 4*sizeof(t_w_t)); } sc_opt = sc; } else { for (j = 0; j < b->cc_size; ++j) { memcpy(b->u[b->cc_node[j]].s, b->u_opt[b->cc_node[j]].s, 2*sizeof(mc_node_t)); memcpy(b->u[b->cc_node[j]].z, b->u_opt[b->cc_node[j]].z, 4*sizeof(t_w_t)); } } if((n_iter%flush) == 0) { reset_mb_g_t_z(mbg); sc = mb_score(mbg, b); for (j = 0; j < b->cc_size; ++j) { memcpy(b->u_opt[b->cc_node[j]].z, b->u[b->cc_node[j]].z, 4*sizeof(t_w_t)); } sc_opt = sc; } // print_sc(opt, mg->e, b, sc_opt, n_iter); } for (j = 0; j < b->cc_size; ++j) { memcpy(b->u[b->cc_node[j]].s, b->u_opt[b->cc_node[j]].s, 2*sizeof(mc_node_t)); memcpy(b->u[b->cc_node[j]].z, b->u_opt[b->cc_node[j]].z, 4*sizeof(t_w_t)); } return n_iter; } void reset_mb_g_t_z(mb_g_t *mbg) { uint32_t k; mb_match_t *ma = mbg->e; for (k = 0; k < ma->n_seq; ++k) ///each block { uint32_t o = ma->idx.a[k] >> 32; uint32_t j, n = (uint32_t)ma->idx.a[k]; mbg->u->u.a[k].z[0] = mbg->u->u.a[k].z[1] = 0; mbg->u->u.a[k].z[2] = mbg->u->u.a[k].z[3] = 0; for (j = 0; j < n; ++j) { const mb_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e); ///a[0]->b[0] if(mbg->u->u.a[t].s[0] > 0) mbg->u->u.a[k].z[0] += e->w[0]; else if(mbg->u->u.a[t].s[0] < 0) mbg->u->u.a[k].z[1] += e->w[0]; ///a[0]->b[1] if(mbg->u->u.a[t].s[1] > 0) mbg->u->u.a[k].z[0] += e->w[1]; else if(mbg->u->u.a[t].s[1] < 0) mbg->u->u.a[k].z[1] += e->w[1]; ///a[1]->b[0] if(mbg->u->u.a[t].s[0] > 0) mbg->u->u.a[k].z[2] += e->w[2]; else if(mbg->u->u.a[t].s[0] < 0) mbg->u->u.a[k].z[3] += e->w[2]; ///a[1]->b[1] if(mbg->u->u.a[t].s[1] > 0) mbg->u->u.a[k].z[2] += e->w[3]; else if(mbg->u->u.a[t].s[1] < 0) mbg->u->u.a[k].z[3] += e->w[3]; } } } void debug_mbg(mb_g_t *mbg, mc_g_t *mg, mc_svaux_t *b) { uint32_t i, k, a_n[2], *a[2], qn, found; int8_t s[2]; for (i = 0; i < mbg->u->u.n; i++)///each block { fprintf(stderr, "i: %u, mbg->u->u.n: %u\n", i, (uint32_t)mbg->u->u.n); decode_mb_node(mbg, i, &(a[0]), &(a_n[0]), &(s[0]), &(a[1]), &(a_n[1]), &(s[1])); found = 0; for (k = 0; k < a_n[0]; k++) { qn = a[0][k]; if(b->s[qn] == s[0]) found = 1; if(b->s[qn] == 0) continue; if(b->s[qn] != s[0]) fprintf(stderr, "ERROR-0-::%s\n", __func__); } if(found == 0) fprintf(stderr, "ERROR-0-*::%s\n", __func__); found = 0; for (k = 0; k < a_n[1]; k++) { qn = a[1][k]; if(b->s[qn] == s[1]) found = 1; if(b->s[qn] == 0) continue; if(b->s[qn] != s[1]) fprintf(stderr, "ERROR-1-::%s\n", __func__); } if(found == 0 && a_n[1] > 0) fprintf(stderr, "ERROR-1-*::%s\n", __func__); } // t_w_t mcw = mc_score_all(mg->e, b); // t_w_t mbw = 0; // mb_match_t *ma = mbg->e; // for (k = 0; k < ma->n_seq; ++k) { // ///a[0] // mbw += -(t_w_t)(mbg->u->u.a[k].s[0]) * (mbg->u->u.a[k].z[0] - mbg->u->u.a[k].z[1]); // ///a[1] // mbw += -(t_w_t)(mbg->u->u.a[k].s[1]) * (mbg->u->u.a[k].z[2] - mbg->u->u.a[k].z[3]); // } } void mc_init_spin_all(const mc_opt_t *opt, mc_g_t *mg, mb_g_t *mbg, mc_svaux_t *b) { uint32_t st, i, k, m, a_n[2], *a[2], qn; for (st = 0, i = 1; i <= mg->e->n_seq; ++i) { if (i == mg->e->n_seq || mg->e->cc[st]>>32 != mg->e->cc[i]>>32) { b->cc_off = st, b->cc_size = i - st; if (b->cc_size >= 2) { mc_init_spin(mg->e, b); } st = i; } } if(!mbg) return; // memcpy(b->s_opt, b->s, sizeof(int8_t)*mg->e->n_seq); ///adjust by block kvec_t(uint32_t) s; kv_init(s); uint32_t ws, we, n_flip; uint8_t *vis = NULL; CALLOC(vis, mg->e->n_seq); t_w_t w = 0, max_w = 0; for (i = 0; i < mbg->u->u.n; i++)///each block { decode_mb_node(mbg, i, &(a[0]), &(a_n[0]), NULL, &(a[1]), &(a_n[1]), NULL); ///if(a_n[0] == 0 || a_n[1] == 0) continue; if(a_n[0] + a_n[1] <= 1) continue; s.n = 0; for (k = 0; k < a_n[0]; k++) { qn = a[0][k]; if(b->s[qn] == 0) continue; // qn <<= 1; if(b->s[qn] < 0) qn += ((uint32_t)1<<31); kv_push(uint32_t, s, qn); } for (k = 0; k < a_n[1]; k++) { qn = a[1][k]; if(b->s[qn] == 0) continue; // qn <<= 1; qn++; if(b->s[qn] > 0) qn += ((uint32_t)1<<31); kv_push(uint32_t, s, qn); } if(s.n == 0) continue; radix_sort_mc32(s.a, s.a + s.n); ws = we = (uint32_t)-1; n_flip = 0; for (st = 0, k = 1; k <= s.n; k++) { if(k == s.n || (s.a[k]>>31) != (s.a[st]>>31)) { w = 0; for (m = st; m < k; m++) vis[(s.a[m]<<1)>>1] = 1; for (m = st; m < k; m++) w += incre_weight(b, mg->e, vis, (s.a[m]<<1)>>1); for (m = st; m < k; m++) vis[(s.a[m]<<1)>>1] = 0; if(ws == (uint32_t)-1 || max_w < w) max_w = w, ws = st, we = k; st = k; n_flip++; } } if(n_flip <= 1) continue; /*******************************for debug************************************/ // t_w_t sc_opt = mc_score_all(mg->e, b); /*******************************for debug************************************/ for (m = ws; m < we; m++) mc_set_spin(mg->e, b, (s.a[m]<<1)>>1, -b->s[(s.a[m]<<1)>>1]); /*******************************for debug************************************/ // t_w_t sc_cur = mc_score_all(mg->e, b); // if(sc_cur - sc_opt != max_w*2) // { // fprintf(stderr, "ERROR-0-::%s\n", __func__); // // fprintf(stderr, "sc_opt: %ld, sc_cur: %ld, max_w: %ld\n", sc_opt, sc_cur, max_w); // fprintf(stderr, "sc_opt: %f, sc_cur: %f, max_w: %f\n", sc_opt, sc_cur, max_w); // } /*******************************for debug************************************/ } kv_destroy(s); free(vis); for (i = 0; i < mbg->u->u.n; i++)///each block { decode_mb_node(mbg, i, &(a[0]), &(a_n[0]), NULL, &(a[1]), &(a_n[1]), NULL); mbg->u->u.a[i].s[0] = mbg->u->u.a[i].s[1] = 0; mbg->u->u.a[i].z[0] = mbg->u->u.a[i].z[1] = 0; mbg->u->u.a[i].z[2] = mbg->u->u.a[i].z[3] = 0; for (k = 0; k < a_n[0]; k++) { qn = a[0][k]; if(b->s[qn] == 0) continue; mbg->u->u.a[i].s[0] = b->s[qn]; break; } for (k = 0; k < a_n[1]; k++) { qn = a[1][k]; if(b->s[qn] == 0) continue; mbg->u->u.a[i].s[1] = b->s[qn]; break; } } reset_mb_g_t_z(mbg); /*******************************for debug************************************/ // debug_mbg(mbg, mg, b); /*******************************for debug************************************/ // memcpy(b->s, b->s_opt, sizeof(int8_t)*mg->e->n_seq); } uint64_t *mb_g_cc_core(mb_match_t *ma) { uint32_t i, x, y, *flag; uint64_t *group; mb_edge_t *o = NULL; kvec_t(uint32_t) stack; kv_init(stack); MALLOC(flag, ma->n_seq); for (i = 0; i < ma->n_seq; ++i) flag[i] = (uint32_t)-1; // connected componets for (i = 0; i < ma->n_seq; ++i) { if (flag[i] != (uint32_t)-1) continue; stack.n = 0; kv_push(uint32_t, stack, i); while (stack.n > 0) { uint32_t k, j, n; stack.n--; k = stack.a[stack.n]; flag[k] = i;///group id // n = (uint32_t)ma->idx[k]; // s = ma->idx[k] >> 32; o = pt_a(*ma, k); n = pt_n(*ma, k); for (j = 0; j < n; ++j) { uint32_t t = ma_y(o[j]); if (flag[t] != (uint32_t)-1) continue; // if (ns == ms) PT_EXPAND(stack, ms); // stack[ns++] = t; kv_push(uint32_t, stack, t); } } } kv_destroy(stack); // precalculate the size of each group CALLOC(group, ma->n_seq); for (i = 0; i < ma->n_seq; ++i) group[i] = (uint64_t)flag[i] << 32 | i; radix_sort_mc64(group, group + ma->n_seq); for (i = 1, x = y = 0; i <= ma->n_seq; ++i) { if (i == ma->n_seq || group[i]>>32 != group[x]>>32) { uint32_t j; for (j = x; j < i; ++j) group[j] = (uint64_t)y << 32 | (uint32_t)group[j];///(group id)|first element in this group ++y, x = i; } } free(flag); return group; } void mb_g_cc(mb_g_t *mbg) { mbg->e->cc = mb_g_cc_core(mbg->e); } void mc_set_by_mbg(mc_g_t *mg, mb_g_t *mbg) { // t_w_t z0 = mc_score_all_advance(mg->e, mg->s.a); // t_w_t z1 = mb_score_all_advance(mg->e, mbg); // if(z0 >= z1) return; uint32_t i, k, qn, *a[2], a_n[2]; int8_t s[2]; for (i = 0; i < mbg->u->u.n; i++) { decode_mb_node(mbg, i, &(a[0]), &(a_n[0]), &(s[0]), &(a[1]), &(a_n[1]), &(s[1])); for (k = 0; k < a_n[0]; k++) { qn = a[0][k]; if(mg->s.a[qn] == 0) continue; mg->s.a[qn] = s[0]; } for (k = 0; k < a_n[1]; k++) { qn = a[1][k]; if(mg->s.a[qn] == 0) continue; mg->s.a[qn] = s[1]; } } } void debug_mb_solve_core(mb_g_t *mbg) { uint32_t i; for (i = 0; i < mbg->u->u.n; i++) { debug_mb_z(mbg, i); } } void print_mb_g_blcok(mb_g_t *mbg) { uint32_t i, k, a_n[2], *a[2]; for (i = 0; i < mbg->u->u.n; i++)///each block { decode_mb_node(mbg, i, &(a[0]), &(a_n[0]), NULL, &(a[1]), &(a_n[1]), NULL); { if(a_n[0] + a_n[1] <= 1) continue; fprintf(stderr, "\nB0-%s\n", a_n[0]>1?"mul":"single"); for (k = 0; k < a_n[0]; k++) { fprintf(stderr,"s-utg%.6ul\n", a[0][k]+1); } fprintf(stderr, "B1-%s\n", a_n[0]>1?"mul":"single"); for (k = 0; k < a_n[1]; k++) { fprintf(stderr,"d-utg%.6ul\n", a[1][k]+1); } } } } void mb_solve_core(mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys) { if(!ref) return; double index_time = yak_realtime(); uint32_t st, i; mb_g_t *mbg = init_mb_g_t(mg->e, ref, is_sys); mb_svaux_t *bb; /**************************init**************************/ if(VERBOSE_CUT) { fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); } mc_svaux_t *b; mc_g_cc(mg->e); b = mc_svaux_init(mg, opt->seed); mc_init_spin_all(opt, mg, mbg, b); mc_svaux_destroy(b); free(mg->e->cc); mg->e->cc = NULL; /**************************init**************************/ mb_g_cc(mbg); bb = mb_svaux_init(mbg, opt->seed); if(VERBOSE_CUT) { fprintf(stderr, "*********before-[M::%s::mc_score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); fprintf(stderr, "*********before-[M::%s::mb_score->%f] ==> Partition\n", __func__, mb_score_all_advance(mg->e, mbg)); /*******************************for debug************************************/ // print_mb_g_blcok(mbg); /*******************************for debug************************************/ } opt->n_perturb = opt->n_b_perturb; for (st = 0, i = 1; i <= mbg->e->n_seq; ++i) { if (i == mbg->e->n_seq || mbg->e->cc[st]>>32 != mbg->e->cc[i]>>32) { mb_solve_cc(opt, mbg, bb, st, i - st); st = i; } } opt->n_perturb = opt->n_s_perturb - opt->n_b_perturb; mc_set_by_mbg(mg, mbg); if(VERBOSE_CUT) { /*******************************for debug************************************/ // debug_mb_solve_core(mbg); /*******************************for debug************************************/ fprintf(stderr, "##############end-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); } destory_mb_g_t(&mbg); mb_svaux_destroy(bb); fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time); } void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub) { double index_time = yak_realtime(); uint32_t st, i; mc_svaux_t *b; mc_bp_t *bp = NULL; mc_g_cc(mg->e); b = mc_svaux_init(mg, opt->seed); if(bub) bp = mc_bp_t_init(mg->e, b, bub, asm_opt.thread_num); /*******************************for debug************************************/ if(bp) mc_init_spin_all(opt, mg, NULL, b); if(bp) mc_solve_bp(bp); /*******************************for debug************************************/ if(VERBOSE_CUT) { fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); } for (st = 0, i = 1; i <= mg->e->n_seq; ++i) { if (i == mg->e->n_seq || mg->e->cc[st]>>32 != mg->e->cc[i]>>32) { mc_solve_cc(opt, mg, b, st, i - st); st = i; } } if(VERBOSE_CUT) { fprintf(stderr, "##############end-[---M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b)); mc_status_all(mg->e, mg->s.a); } if(bp) mc_solve_bp(bp); ///mc_write_info(g, b); mc_svaux_destroy(b); if(bp) destroy_mc_bp_t(&bp); fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time); } void set_p_flag(mc_g_t *mg, uint32_t uID, uint8_t* trio_flag, trans_chain* t_ch, int8_t s) { uint32_t i; ma_utg_t *u = &(mg->ug->u.a[uID]); for (i = 0; i < u->n; i++) { trio_flag[u->a[i]>>33] |= SET_TRIO; if(t_ch->ir_het[u->a[i]>>33] == N_HET) continue; if(s == 0) { if(t_ch->ir_het[u->a[i]>>33]&P_HET)//special case { trio_flag[u->a[i]>>33] |= FATHER; } continue; } trio_flag[u->a[i]>>33] |= (s > 0? FATHER:MOTHER); } } void filter_ovlp_by_mc(ma_ug_t *ug, asg_t *read_g, uint32_t uID, hap_overlaps_list* ha, mc_match_t *ma, int8_t *s) { mc_edge_t *o = pt_a(*ma, uID); uint32_t n = pt_n(*ma, uID), k, qn, tn; hap_overlaps *p = NULL; int index; for (k = 0; k < n; ++k) { qn = ma_x(o[k]); tn = ma_y(o[k]); p = NULL; if((s[qn]*s[tn])!=-1) continue; index = get_specific_hap_overlap(&(ha->x[qn]), qn, tn); if(index != -1 && ha->x[qn].a.a[index].score == (long long)o[k].w) { p = &(ha->x[qn].a.a[index]); } else { index = get_specific_hap_overlap(&(ha->x[tn]), tn, qn); if(index != -1 && ha->x[tn].a.a[index].score == (long long)o[k].w) { p = &(ha->x[tn].a.a[index]); } } if(!p) fprintf(stderr, "ERROR\n"); p->status = FLIP; } } void clean_ovlp_by_mc(mc_g_t *mg, hap_overlaps_list* ha) { uint32_t v, i, k, qn, tn, types[4]; types[X2Y] = Y2X; types[Y2X] = X2Y; types[XCY] = YCX; types[YCX] = XCY; int index; hap_overlaps *x = NULL, *y = NULL; for (i = 0; i < mg->e->n_seq; ++i) { filter_ovlp_by_mc(mg->ug, mg->rg, i, ha, mg->e, mg->s.a); } for (v = 0; v < ha->num; v++) { for (i = 0; i < ha->x[v].a.n; i++) { qn = ha->x[v].a.a[i].xUid; tn = ha->x[v].a.a[i].yUid; x = &(ha->x[v].a.a[i]); if(x->status != FLIP) continue; index = get_specific_hap_overlap(&(ha->x[tn]), tn, qn); if(index != -1) { y = &(ha->x[tn].a.a[index]); set_reverse_hap_overlap(y, x, types); y->status = FLIP; } if(index == -1) fprintf(stderr, "ERROR\n"); } } for (v = 0; v < ha->num; v++) { for (i = k = 0; i < ha->x[v].a.n; i++) { if(ha->x[v].a.a[i].status != FLIP) continue; ha->x[v].a.a[k] = ha->x[v].a.a[i]; k++; } ha->x[v].a.n = k; } } void filter_ta_by_mc(ma_ug_t *ug, asg_t *read_g, uint32_t uID, kv_u_trans_t* ta, mc_match_t *ma, int8_t *s) { mc_edge_t *o = pt_a(*ma, uID); uint32_t n = pt_n(*ma, uID), k, qn, tn; u_trans_t *p = NULL; for (k = 0; k < n; ++k) { qn = ma_x(o[k]); tn = ma_y(o[k]); p = NULL; if((s[qn]*s[tn])!=-1) continue; get_u_trans_spec(ta, qn, tn, &p, NULL); if(p && p->nw == o[k].w) { p->del = 0; } else { get_u_trans_spec(ta, tn, qn, &p, NULL); if(p && p->nw == o[k].w) p->del = 0; } if(!p) fprintf(stderr, "ERROR-ta-p\n"); } } void clean_ta_by_mc(mc_g_t *mg, kv_u_trans_t *ta) { uint32_t v, i; u_trans_t *p = NULL; for (i = 0; i < ta->n; i++) ta->a[i].del = 1; for (i = 0; i < mg->e->n_seq; ++i) filter_ta_by_mc(mg->ug, mg->rg, i, ta, mg->e, mg->s.a); for (i = v = 0; i < ta->n; i++) { if(ta->a[i].del) continue; ta->a[v++] = ta->a[i]; } ta->n = v; for (i = 0; i < v; i++) { kv_pushp(u_trans_t, *ta, &p); (*p) = ta->a[i]; p->qn = ta->a[i].tn; p->qs = ta->a[i].ts; p->qe = ta->a[i].te; p->tn = ta->a[i].qn; p->ts = ta->a[i].qs; p->te = ta->a[i].qe; } kt_u_trans_t_idx(ta, mg->ug->g->n_seq); } void p_nodes(mc_g_t *mg, trans_chain* t_ch, uint8_t* trio_flag) { uint32_t i; for (i = 0; i < mg->e->n_seq; ++i) { set_p_flag(mg, i, trio_flag, t_ch, mg->s.a[i]); } } void write_mc_g_t(mc_opt_t *opt, mc_g_t *mg, const char *name) { FILE* fp = fopen(name, "w"); fwrite(opt, sizeof(mc_opt_t), 1, fp); fwrite(&(mg->s.n), sizeof(mg->s.n), 1, fp); fwrite(mg->s.a, sizeof(mc_node_t), mg->s.n, fp); fwrite(&(mg->e->n_seq), sizeof(mg->e->n_seq), 1, fp); fwrite(&(mg->e->ma.n), sizeof(mg->e->ma.n), 1, fp); fwrite(mg->e->ma.a, sizeof(mc_edge_t), mg->e->ma.n, fp); fwrite(&(mg->e->idx.n), sizeof(mg->e->idx.n), 1, fp); fwrite(mg->e->idx.a, sizeof(uint64_t), mg->e->idx.n, fp); fclose(fp); } mc_g_t* load_mc_g_t(mc_opt_t *opt, const char *name) { FILE* fp = NULL; fp = fopen(name, "r"); if(!fp) return NULL; uint64_t flag = 0; mc_g_t *mg = NULL; CALLOC(mg, 1); kv_init(mg->s); CALLOC(mg->e, 1); flag += fread(opt, sizeof(mc_opt_t), 1, fp); flag += fread(&(mg->s.n), sizeof(mg->s.n), 1, fp); mg->s.m = mg->s.n; MALLOC(mg->s.a, mg->s.n); flag += fread(mg->s.a, sizeof(mc_node_t), mg->s.n, fp); flag += fread(&(mg->e->n_seq), sizeof(mg->e->n_seq), 1, fp); flag += fread(&(mg->e->ma.n), sizeof(mg->e->ma.n), 1, fp); mg->e->ma.m = mg->e->ma.n; MALLOC(mg->e->ma.a, mg->e->ma.n); flag += fread(mg->e->ma.a, sizeof(mc_edge_t), mg->e->ma.n, fp); flag += fread(&(mg->e->idx.n), sizeof(mg->e->idx.n), 1, fp); mg->e->idx.m = mg->e->idx.n; MALLOC(mg->e->idx.a, mg->e->idx.n); flag += fread(mg->e->idx.a, sizeof(uint64_t), mg->e->idx.n, fp); fclose(fp); return mg; } void debug_mc_g_t(const char* name) { mc_opt_t opt; mc_g_t *mg = load_mc_g_t(&opt, name); mc_solve_core(&opt, mg, NULL); destory_mc_g_t(&mg); exit(1); } void print_hap_s(int8_t *s, uint32_t sn) { fprintf(stderr, "\n[M::%s]\n", __func__); uint32_t i; for (i = 0; i < sn; i++) { fprintf(stderr,"utg%.6ul\t", i + 1); if(s[i] > 0) fprintf(stderr, "h%u\n", 1); else if(s[i] < 0) fprintf(stderr, "h%u\n", 2); else fprintf(stderr,"\n"); } } void dump_debug_phasing(const char* fn, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref); void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov, int is_dump) { if(is_dump) { dump_debug_phasing(MC_NAME, ta, ug, read_g, f_rate, renew_s, s, is_sys, bub, ref); } mc_opt_t opt; mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed); mc_g_t *mg = init_mc_g_t(ug, read_g, s, renew_s); update_mc_edges(mg, ovlp, ta, t_ch, f_rate, is_sys); fprintf(stderr, "[M::%s:: # edges: %u]\n", __func__, (uint32_t)mg->e->ma.n); mb_solve_core(&opt, mg, ref, is_sys); ///debug_mc_g_t(mg); // if(renew_s == 0) write_mc_g_t(&opt, mg, MC_NAME); mc_solve_core(&opt, mg, bub); if((asm_opt.flag & HA_F_PARTITION) && t_ch) { p_nodes(mg, t_ch, trio_flag); } if(clean_ov){ if(ovlp) clean_ovlp_by_mc(mg, ovlp); if(ta) clean_ta_by_mc(mg, ta); } // print_hap_s(s, ug->u.n); destory_mc_g_t(&mg); } void comp(int m, int N, int M, mcb_t *p, int *c) { if (m == M + 1) { int i; mcg_node_t x = 0; for (i = 0; i < M; i++) x |= ((mcg_node_t)1<<(c[i+1]-1)); kv_push(mcg_node_t, *p, x); } else { for (c[m] = c[m - 1] + 1; c[m] <= N - M + m; c[m]++) { comp(m + 1, N, M, p, c); } } } void get_mcb(uint32_t n, uint32_t m, mcb_t *p, int *c) { memset(c, 0, sizeof(int)*n+1); p->n = 0; comp(1, n, m, p, c); } mc_gg_t *init_mc_gg_t(uint32_t un, kv_gg_status *s, uint16_t hapN) { uint32_t i; int *c = NULL; CALLOC(c, hapN+1); mc_gg_t *p = NULL; CALLOC(p, 1); // p->ug = ug; p->rg = read_g; p->un = un; p->s = s; p->hN = hapN; CALLOC(p->m.a, p->hN); p->m.n = p->m.m = hapN; for (i = 0; i < p->hN; i++) get_mcb(hapN, i+1, &(p->m.a[i]), c); p->mask = (1<mask--; free(c); return p; } void destory_mc_gg_t(mc_gg_t **p) { uint32_t i; if(!p || !(*p)) return; for (i = 0; i < (*p)->m.m; i++) { free((*p)->m.a[i].a); } free((*p)->m.a); if((*p)->e) { kv_destroy((*p)->e->idx); kv_destroy((*p)->e->ma); free((*p)->e->cc); free((*p)->e); } free((*p)); } kv_gg_status *init_mc_gg_status(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint64_t t_cov, uint16_t hapN) { fprintf(stderr, "t_cov-%lu\n", t_cov); uint64_t *covs = NULL, i, k, k_i, c = t_cov/hapN, c_min, c_max; uint32_t len[2]; uint8_t *vis = NULL; CALLOC(vis, read_g->n_seq); kv_gg_status *p = NULL; CALLOC(covs, hapN); for (i = 0; i < hapN; i++) { c_min = ((i+1)*c) - (0.6*c); c_max = ((i+1)*c) + (0.6*c); if(i == 0) c_min = 0; if(i == hapN) c_max = (uint32_t)-1; covs[i] = (c_max<<32)|c_min; } CALLOC(p, 1); p->n = p->m = ug->u.n; CALLOC(p->a, p->n); for (i = 0; i < p->n; i++) { c = get_utg_cov(ug, i, read_g, coverage_cut, sources, ruIndex, vis); p->a[i].h[0] = p->a[i].h[1] = (uint16_t)-1; p->a[i].s = 0; k_i = 0; p->a[i].hw[0] = 1; p->a[i].hw[1] = 0; p->a[i].hc = 0; for (k = 0; k < hapN; k++) { c_min = (uint32_t)covs[k]; c_max = covs[k]>>32; if(c < c_min || c >= c_max) continue; p->a[i].h[k_i++] = k + 1; } if(p->a[i].h[1] == (uint16_t)-1) continue; len[0] = (c >= (p->a[i].h[0]*(t_cov/hapN))? c - (p->a[i].h[0]*(t_cov/hapN)) : (p->a[i].h[0]*(t_cov/hapN)) - c); len[1] = (c >= (p->a[i].h[1]*(t_cov/hapN))? c - (p->a[i].h[1]*(t_cov/hapN)) : (p->a[i].h[1]*(t_cov/hapN)) - c); if(len[0] > len[1]) { k_i = p->a[i].h[0]; p->a[i].h[0] = p->a[i].h[1]; p->a[i].h[1] = k_i; k_i = len[0]; len[0] = len[1]; len[1] = k_i; } p->a[i].hw[0] = (double)len[1]/(double)(len[0]+len[1]); p->a[i].hw[1] = (double)len[0]/(double)(len[0]+len[1]); } free(covs); free(vis); return p; } void update_mc_edges_general(mc_gg_t *mg, kv_u_trans_t *ta, uint16_t hapN) { uint32_t i, k; mc_edge_t *ma = NULL; CALLOC(mg->e, 1); mg->e->n_seq = mg->un; kv_init(mg->e->idx); kv_init(mg->e->ma); for (i = 0; i < ta->n; ++i) { if(ta->a[i].del) continue; if(mg->s->a[ta->a[i].qn].h[0] >= hapN && mg->s->a[ta->a[i].qn].h[1] >= hapN) continue; if(mg->s->a[ta->a[i].tn].h[0] >= hapN && mg->s->a[ta->a[i].tn].h[1] >= hapN) continue; kv_pushp(mc_edge_t, mg->e->ma, &ma); ma->x = (uint64_t)ta->a[i].qn << 32 | ta->a[i].tn; ma->w = w_cast((ta->a[i].nw)); } for (i = k = 0; i < mg->e->ma.n; i++) { if(mg->e->ma.a[i].w == 0) continue; mg->e->ma.a[k] = mg->e->ma.a[i]; k++; } mg->e->ma.n = k; radix_sort_mce(mg->e->ma.a, mg->e->ma.a + mg->e->ma.n); mc_edges_idx(mg->e); mc_edges_symm(mg->e); } typedef struct { t_w_t *z; uint32_t hapN; } mc_poy_t; typedef struct { uint64_t x; // RNG uint32_t cc_off, cc_size; kvec_t(uint64_t) cc_edge; uint32_t *cc_node; uint32_t *bfs, *bfs_mark; mc_poy_t *z, *z_opt;///keep scores to nodes(1) and nodes(-1) mc_gg_status *s, *s_opt; mcb_t *m; mcg_node_t mask; uint32_t hapN; } mcgg_svaux_t; mc_poy_t *init_mc_poy_t(uint32_t un, uint32_t hapN) { uint32_t i; mc_poy_t *z = NULL; CALLOC(z, 1); MALLOC(z->z, ((uint32_t)un<z[i] = 0; z->hapN = hapN; return z; } void destroy_mc_poy_t(mc_poy_t **p) { if(!p || !(*p)) return; free((*p)->z); free(*p); } mcgg_svaux_t *mcgg_svaux_init(const mc_gg_t *mg, uint64_t x, uint32_t hapN) { uint32_t st, i, max_cc = 0; mc_match_t *ma = mg->e; mcgg_svaux_t *b; CALLOC(b, 1); b->x = x; for (st = 0, i = 1; i <= ma->n_seq; ++i) if (i == ma->n_seq || ma->cc[st]>>32 != ma->cc[i]>>32) max_cc = max_cc > i - st? max_cc : i - st, st = i; kv_init(b->cc_edge); MALLOC(b->cc_node, max_cc); b->s = mg->s->a; CALLOC(b->s_opt, ma->n_seq); MALLOC(b->bfs, ma->n_seq); MALLOC(b->bfs_mark, ma->n_seq); memset(b->bfs_mark, -1, ma->n_seq*sizeof(uint32_t)); b->z = init_mc_poy_t(ma->n_seq, hapN); b->z_opt = init_mc_poy_t(ma->n_seq, hapN); b->m = mg->m.a; b->mask = ((mcg_node_t)1)<mask--; b->hapN = hapN; return b; } void mcgg_svaux_destroy(mcgg_svaux_t *b) { b->s = NULL; kv_destroy(b->cc_edge); free(b->cc_node); free(b->s); free(b->s_opt); destroy_mc_poy_t(&(b->z)); destroy_mc_poy_t(&(b->z_opt)); free(b->bfs); free(b->bfs_mark); free(b); } static inline mcg_node_t kr_drand_node(uint64_t id, mcgg_svaux_t *b, uint16_t *hc) { uint16_t hn = b->s[id].h[0]-1; if(hc) (*hc) = 0; b->x = kr_splitmix64(b->x); if(b->s[id].h[1] != (uint16_t)-1) { union { uint64_t i; double d; } u; u.i = 0x3FFULL << 52 | (b->x) >> 12; if((u.d - 1.0) > b->s[id].hw[0]) { hn = b->s[id].h[1]-1; if(hc) (*hc) = 1; } } return b->m[hn].a[(b->x)%b->m[hn].n]; } static inline mcg_node_t kr_drand_node_ref(uint64_t id, mcgg_svaux_t *b, uint16_t *hc, mcg_node_t ref, uint64_t rev) { uint64_t i, k, m, mm, mn, mi, cn, hn = b->s[id].h[0]; mcg_node_t t; if(hc) (*hc) = 0; b->x = kr_splitmix64(b->x); if(b->s[id].h[1] != (uint16_t)-1) { union { uint64_t i; double d; } u; u.i = 0x3FFULL << 52 | (b->x) >> 12; if((u.d - 1.0) > b->s[id].hw[0]) { hn = b->s[id].h[1]; if(hc) (*hc) = 1; } } if(rev) ref ^= (mcg_node_t)-1; ref &= b->mask; t = ref; cn = 0; while (t) { cn += t&1; t >>= 1; } if(cn == hn) return ref; if(cn>=hn) mm=cn, m=1, mn=cn-hn;///1->0 else mm=b->hapN-cn, m=0, mn=hn-cn;///0->1 for (i = 0; i < mn; i++) { mi = (b->x%mm); for (k = 0; k < b->hapN; k++) { if(((ref>>k)&1)!=m) continue; if(mi == 0) { ref ^= ((mcg_node_t)1<h[s->hc]; mcg_node_t ref = s->s; for (i = cn = 0; i < 32; i++) cn += ((ref>>i)&1); if (cn != hn) fprintf(stderr, "%s-ERROR-cn, cn-%u, hn-%u\n", cmd, cn, hn); else fprintf(stderr, "%s-pass-cn, cn-%u, hn-%u\n", cmd, cn, hn); } void mcgg_reset_z(const mc_match_t *ma, mcgg_svaux_t *b) { uint32_t i; for (i = 0; i < b->cc_size; ++i) { uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///uid uint32_t o = ma->idx.a[k] >> 32; uint32_t j, n = (1<z->hapN); for (j = 0; j < n; j++) mcp_de(*(b->z), k, j) = 0; n = (uint32_t)ma->idx.a[k]; for (j = 0; j < n; ++j) { const mc_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e); mcp_de(*(b->z), k, b->s[t].s) += e->w; } } } t_w_t mcgg_score(const mc_match_t *ma, mcgg_svaux_t *b) { uint32_t i, j, n = (1<z->hapN); t_w_t z = 0; for (i = 0; i < b->cc_size; ++i) { uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///uid for (j = 0; j < n; j++) { z += ((b->s[k].s&((mcg_node_t)j))?-mcp_de(*(b->z), k, j):mcp_de(*(b->z), k, j)); } } return z; } t_w_t mcgg_init_spin(const mc_match_t *ma, mcgg_svaux_t *b) { uint32_t i; b->cc_edge.n = 0; for (i = 0; i < b->cc_size; ++i) {///how many nodes uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///node id b->cc_node[i] = k; if(b->s[k].s == 0) break; } if(i >= b->cc_size) goto passed; for (i = 0; i < b->cc_size; ++i) {///how many nodes uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///node id uint32_t o = ma->idx.a[k] >> 32;///cc group id uint32_t n = (uint32_t)ma->idx.a[k], j; b->cc_node[i] = k; for (j = 0; j < n; ++j) { w_t w = ma->ma.a[o + j].w; w = w > 0? w : -w; kv_push(uint64_t, b->cc_edge, (uint64_t)((uint32_t)-1 - ((uint32_t)w)) << 32 | (o + j)); } } radix_sort_mc64(b->cc_edge.a, b->cc_edge.a + b->cc_edge.n); for (i = 0; i < b->cc_edge.n; ++i) { // from the strongest edge to the weakest const mc_edge_t *e = &ma->ma.a[(uint32_t)b->cc_edge.a[i]]; uint32_t n1 = ma_x(*e), n2 = ma_y(*e); if (b->s[n1].s == 0 && b->s[n2].s == 0) { b->s[n1].s = kr_drand_node(n1, b, &(b->s[n1].hc)); // debug_hapM(&(b->s[n1]), "s0"); b->s[n2].s = kr_drand_node_ref(n2, b, &(b->s[n2].hc), b->s[n1].s, e->w>0?1:0); // debug_hapM(&(b->s[n2]), "s1"); } else if(b->s[n1].s == 0) { b->s[n1].s = kr_drand_node_ref(n1, b, &(b->s[n1].hc), b->s[n2].s, e->w>0?1:0); // debug_hapM(&(b->s[n1]), "s2"); } else if(b->s[n2].s == 0) { b->s[n2].s = kr_drand_node_ref(n2, b, &(b->s[n2].hc), b->s[n1].s, e->w>0?1:0); // debug_hapM(&(b->s[n2]), "s3"); } } passed: mcgg_reset_z(ma, b); return mcgg_score(ma, b); } static mcg_node_t get_max_m(uint64_t id, mcgg_svaux_t *b) { uint32_t hn = b->s[id].h[0]-1, k, j, n = (1<z->hapN); mcg_node_t m, *p = NULL; t_w_t z, max_z = -(1<<30); for (k = 0; k < b->m[hn].n; k++) { m = b->m[hn].a[k]; z = 0; for (j = 0; j < n; j++) { z += ((m&((mcg_node_t)j))?-mcp_de(*(b->z), id, j):mcp_de(*(b->z), id, j)); } if(!p || max_z < z || (max_z == z && m == b->s[id].s)) p = &(b->m[hn].a[k]), max_z = z; } if(b->s[id].h[1] != (uint16_t)-1) { hn = b->s[id].h[1]-1; for (k = 0; k < b->m[hn].n; k++) { m = b->m[hn].a[k]; z = 0; for (j = 0; j < n; j++) { z += ((m&((mcg_node_t)j))?-mcp_de(*(b->z), id, j):mcp_de(*(b->z), id, j)); } if(!p || max_z < z || (max_z == z && m == b->s[id].s)) p = &(b->m[hn].a[k]), max_z = z; } } return (*p); } ///k is uid static void mcgg_set_spin(const mc_match_t *ma, mcgg_svaux_t *b, uint32_t k, mcg_node_t s, const char* cmd) { uint32_t o, j, n; mcg_node_t s0 = b->s[k].s; /*******************************for debug************************************/ // mcg_node_t t = s; // o = 0; // while (t) { // o += (t&1); t>>=1; // } // if(o != b->s[k].h[0] && o != b->s[k].h[1]) fprintf(stderr, "cmd-%s, ERROR-mcgg-1\n", cmd); // if(s0 == s) fprintf(stderr, "cmd-%s, ERROR-mcgg-2\n", cmd); /*******************************for debug************************************/ if (s0 == s) return; o = ma->idx.a[k] >> 32; n = (uint32_t)ma->idx.a[k]; for (j = 0; j < n; ++j) { const mc_edge_t *e = &ma->ma.a[o + j]; uint32_t t = ma_y(*e);///1->z[0]; (-1)->z[1]; mcp_de(*(b->z), t, s0) -= e->w; mcp_de(*(b->z), t, s) += e->w; } b->s[k].s = s; } static t_w_t mcgg_optimize_local(const mc_opt_t *opt, const mc_match_t *ma, mcgg_svaux_t *b, uint32_t *n_iter) { uint32_t i, n_flip = 0; int32_t n_iter_local = 0; mcg_node_t ms; while (n_iter_local < opt->max_iter) { ++(*n_iter); ks_shuffle_uint32_t(b->cc_size, b->cc_node, &b->x); for (i = n_flip = 0; i < b->cc_size; ++i) { uint32_t k = b->cc_node[i];///uid ms = get_max_m(k, b); if(ms != b->s[k].s) { mcgg_set_spin(ma, b, k, ms, __func__);///no need to change the score of k itself // debug_hapM(&(b->s[k]), "s4"); ++n_flip; } } ++n_iter_local; if (n_flip == 0) break; } return mcgg_score(ma, b); } void inline back_status(mcgg_svaux_t *b, uint32_t id, uint32_t to_opt) { if(to_opt) { b->s_opt[id] = b->s[id]; memcpy(b->z_opt->z+(id<z->hapN), b->z->z+(id<z->hapN), (1<z->hapN)*sizeof(t_w_t)); } else { b->s[id] = b->s_opt[id]; memcpy(b->z->z+(id<z->hapN), b->z_opt->z+(id<z->hapN), (1<z->hapN)*sizeof(t_w_t)); } } /** static inline mcg_node_t kr_drand_node_new(uint64_t id, mcgg_svaux_t *b) { if(id == 10070) fprintf(stderr, "id-%lu, h[0]-%u, h[1]-%u\n", id, b->s[id].h[0], b->s[id].h[1]); uint32_t hn = b->s[id].h[0]-1, is_old = 1, k; if((b->s[id].h[1] != (uint16_t)-1) && (kr_drand_r(&b->x) > b->s[id].hw[0])) { hn = b->s[id].h[1]-1; is_old = 0; } if(id == 10070) fprintf(stderr, "id-%lu, hn-%u, b->m[hn].n-%u, is_old-%u\n", id, hn, b->m[hn].n, is_old); b->x = kr_splitmix64(b->x); k = b->x%(b->m[hn].n-is_old); if(id == 10070) fprintf(stderr, "id-%lu, hn-%u, k-%u\n", id, hn, k); if(b->m[hn].a[k] == b->s[id].s) k = b->m[hn].n-1; return b->m[hn].a[k]; } **/ static inline mcg_node_t kr_drand_node_new(uint64_t id, mcgg_svaux_t *b) { uint32_t hn = b->s[id].h[0]-1, k; if((b->s[id].h[1] != (uint16_t)-1) && (kr_drand_r(&b->x) > b->s[id].hw[0])) { hn = b->s[id].h[1]-1; } b->x = kr_splitmix64(b->x); k = b->x%(b->m[hn].n); if(b->m[hn].a[k] == b->s[id].s && b->m[hn].n > 1) { k = b->x%(b->m[hn].n-1); if(b->m[hn].a[k] == b->s[id].s) k = b->m[hn].n-1; } return b->m[hn].a[k]; } static void mcgg_perturb(const mc_opt_t *opt, const mc_match_t *ma, mcgg_svaux_t *b) { uint32_t i; for (i = 0; i < b->cc_size; ++i) { uint32_t k = (uint32_t)ma->cc[b->cc_off + i];///node id double y; y = kr_drand_r(&b->x); if (y < opt->f_perturb) mcgg_set_spin(ma, b, k, kr_drand_node_new(k, b), __func__); } } static uint32_t mcgg_bfs(const mc_match_t *ma, mcgg_svaux_t *b, uint32_t k0, uint32_t bfs_round, uint32_t max_size) { uint32_t i, n_bfs = 0, st, en, r; b->bfs[n_bfs++] = k0, b->bfs_mark[k0] = k0; st = 0, en = n_bfs; for (r = 0; r < bfs_round; ++r) { for (i = st; i < en; ++i) { uint32_t k = b->bfs[i]; uint32_t o = ma->idx.a[k] >> 32; uint32_t n = (uint32_t)ma->idx.a[k], j; for (j = 0; j < n; ++j) { uint32_t t = (uint32_t)ma->ma.a[o + j].x; if (b->bfs_mark[t] != k0) b->bfs[n_bfs++] = t, b->bfs_mark[t] = k0; } } st = en, en = n_bfs; if (max_size > 0 && n_bfs > max_size) break; } return n_bfs; } ///bfs_round is 3 static void mcgg_perturb_node(const mc_opt_t *opt, const mc_match_t *ma, mcgg_svaux_t *b, int32_t bfs_round) { uint32_t i, k, n_bfs = 0; k = (uint32_t)(kr_drand_r(&b->x) * b->cc_size + .499); if(k >= b->cc_size) k = b->cc_size - 1; k = (uint32_t)ma->cc[b->cc_off + k];///node id n_bfs = mcgg_bfs(ma, b, k, bfs_round, (int32_t)(b->cc_size * opt->f_perturb)); for (i = 0; i < n_bfs; ++i) mcgg_set_spin(ma, b, b->bfs[i], kr_drand_node_new(b->bfs[i], b), __func__); } void print_mcgg_node(const mc_match_t *ma, mcgg_svaux_t *b, uint32_t id) { fprintf(stderr, "[M::%s::utg%.6ul-hap%u]\n", __func__, id, b->s[id].s); w_t w[128]; uint32_t o, n, i, hn = (1<z->hapN); for (i = 0; i < hn; i++) w[i] = 0; o = ma->idx.a[id] >> 32; n = (uint32_t)ma->idx.a[id]; for (i = 0; i < n; ++i) w[b->s[ma_y(ma->ma.a[o + i])].s] += ma->ma.a[o + i].w; for (i = 0; i < hn; i++) fprintf(stderr, "w[%u]-%f, z[%u]-%f\n", i, w[i], i, mcp_de(*(b->z), id, i)); } uint32_t mcgg_solve_cc(const mc_opt_t *opt, const mc_gg_t *mg, mcgg_svaux_t *b, uint32_t cc_off, uint32_t cc_size) { // double t0, t1, tt0, tt1; uint32_t j, k, n_iter = 0, flush = opt->max_iter * 50; t_w_t sc_opt = -(1<<30), sc;///problem-w b->cc_off = cc_off, b->cc_size = cc_size; if (b->cc_size < 2) return 0; sc_opt = mcgg_init_spin(mg->e, b); if (b->cc_size == 2) return 0; for (j = 0; j < b->cc_size; ++j) back_status(b, b->cc_node[j], 1); sc = mcgg_optimize_local(opt, mg->e, b, &n_iter); if (sc > sc_opt) { for (j = 0; j < b->cc_size; ++j) back_status(b, b->cc_node[j], 1); sc_opt = sc; } else { for (j = 0; j < b->cc_size; ++j) back_status(b, b->cc_node[j], 0); } fprintf(stderr, "\nBeg-[M::%s::score->%f]\n", __func__, mcgg_score(mg->e, b)); // print_mcgg_node(mg->e, b, 3838); // print_mcgg_node(mg->e, b, 36880); // tt0 = tt1 = 0; for (k = 0; k < (uint32_t)opt->n_perturb; ++k) { // t0 = yak_realtime(); if (k&1) mcgg_perturb(opt, mg->e, b); else mcgg_perturb_node(opt, mg->e, b, 3); // tt0 += yak_realtime()-t0; // t1 = yak_realtime(); sc = mcgg_optimize_local(opt, mg->e, b, &n_iter); // tt1 += yak_realtime()-t1; // fprintf(stderr, "++(%u) sc_after: %f\n", k, mcgg_score(mg->e, b)); if (sc > sc_opt) { for (j = 0; j < b->cc_size; ++j) back_status(b, b->cc_node[j], 1); sc_opt = sc; } else { for (j = 0; j < b->cc_size; ++j) back_status(b, b->cc_node[j], 0); } if((n_iter%flush) == 0) { mcgg_reset_z(mg->e, b); sc = mcgg_score(mg->e, b); for (j = 0; j < b->cc_size; ++j) back_status(b, b->cc_node[j], 1); sc_opt = sc; } // fprintf(stderr, "\n"); // print_mcgg_node(mg->e, b, 3838); // print_mcgg_node(mg->e, b, 36880); // if((k&31)==0)fprintf(stderr, "+++(%u) sc: %f, sc_opt: %f, tt0: %.3f, tt1: %.3f\n", k, sc, sc_opt, tt0, tt1); } for (j = 0; j < b->cc_size; ++j) back_status(b, b->cc_node[j], 0); fprintf(stderr, "End-[M::%s::score->%f]\n", __func__, mcgg_score(mg->e, b)); return n_iter; } void mc_solve_core_genral(const mc_opt_t *opt, mc_gg_t *mg, uint32_t hapN) { double index_time = yak_realtime(); uint32_t st, i; mcgg_svaux_t *b; mc_g_cc(mg->e); b = mcgg_svaux_init(mg, opt->seed, hapN); // if(VERBOSE_CUT) // { // fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); // } for (st = 0, i = 1; i <= mg->e->n_seq; ++i) { if (i == mg->e->n_seq || mg->e->cc[st]>>32 != mg->e->cc[i]>>32) { mcgg_solve_cc(opt, mg, b, st, i - st); st = i; } } // if(VERBOSE_CUT) // { // fprintf(stderr, "##############end-[---M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b)); // } ///mc_write_info(g, b); mcgg_svaux_destroy(b); fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time); } void print_mcb(mc_gg_t *mg) { uint32_t i, k, m; mcg_node_t t; mcb_t *p; for (i = 0; i < mg->m.n; i++) { p = mcb_pat(*mg, i + 1); fprintf(stderr, "# haplotypes: %u, # combination: %u\n", i+1, (uint32_t)p->n); for (k = 0; k < p->n; k++) { t = p->a[k]; for (m = 0; m < 32; m++) { if((t>>m)&1) fprintf(stderr, "%u\t", m); } fprintf(stderr, "\n"); } } } void print_hap_p(kv_gg_status *s) { fprintf(stderr, "\n[M::%s]\n", __func__); uint32_t i, h; mcg_node_t m; for (i = 0; i < s->n; i++) { fprintf(stderr,"utg%.6ul\t", i + 1); m = s->a[i].s; h = 0; while (m) { h++; if(m&1) fprintf(stderr, "h%u\t", h); m>>=1; } fprintf(stderr,"\n"); } } void write_mc_gg_dump(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s, uint16_t hapN, const char* fn) { fprintf(stderr, "\n[M::%s]\n", __func__); char *buf = (char*)calloc(strlen(fn) + 50, 1); sprintf(buf, "%s.hic.dbg.dump.bin", fn); FILE* fp = fopen(buf, "w"); fwrite(&ta->n, sizeof(ta->n), 1, fp); fwrite(ta->a, sizeof(u_trans_t), ta->n, fp); fwrite(&ta->idx.n, sizeof(ta->idx.n), 1, fp); fwrite(ta->idx.a, sizeof(uint64_t), ta->idx.n, fp); fwrite(&un, sizeof(un), 1, fp); fwrite(&(s->n), sizeof(s->n), 1, fp); fwrite(s->a, sizeof(mc_gg_status), s->n, fp); fwrite(&hapN, sizeof(hapN), 1, fp); fclose(fp); free(buf); } uint32_t load_mc_gg_dump(kv_u_trans_t **rta, uint32_t *un, kv_gg_status **rs, uint16_t *hapN, const char* fn) { fprintf(stderr, "\n[M::%s]\n", __func__); kv_u_trans_t *ta = NULL; kv_gg_status *s = NULL; uint64_t flag = 0; char *buf = (char*)calloc(strlen(fn) + 25, 1); sprintf(buf, "%s.hic.dbg.dump.bin", fn); FILE* fp = NULL; fp = fopen(buf, "r"); if(!fp) { free(buf); return 0; } CALLOC(ta, 1); flag += fread(&ta->n, sizeof(ta->n), 1, fp); ta->m = ta->n; MALLOC(ta->a, ta->n); flag += fread(ta->a, sizeof(u_trans_t), ta->n, fp); flag += fread(&ta->idx.n, sizeof(ta->idx.n), 1, fp); ta->idx.m = ta->idx.n; MALLOC(ta->idx.a, ta->idx.n); flag += fread(ta->idx.a, sizeof(uint64_t), ta->idx.n, fp); flag += fread(un, sizeof(*un), 1, fp); CALLOC(s, 1); flag += fread(&(s->n), sizeof(s->n), 1, fp); s->m = s->n; MALLOC(s->a, s->n); flag += fread(s->a, sizeof(mc_gg_status), s->n, fp); flag += fread(hapN, sizeof(*hapN), 1, fp); *rta = ta; *rs = s; fclose(fp); free(buf); return 1; } mc_g_t* to_mc_g_t(kv_u_trans_t *ta, kv_gg_status *s, uint32_t un) { fprintf(stderr, "[M::%s]\n", __func__); uint32_t i, k; mc_edge_t *ma = NULL; mc_g_t *mg = NULL; CALLOC(mg, 1); CALLOC(mg->e, 1); mg->e->n_seq = un; kv_init(mg->e->idx); kv_init(mg->e->ma); // double sc = get_w_scale(ta); // fprintf(stderr, "sc: %f\n", sc); for (i = 0; i < ta->n; ++i) { if(ta->a[i].del) continue; kv_pushp(mc_edge_t, mg->e->ma, &ma); ma->x = (uint64_t)ta->a[i].qn << 32 | ta->a[i].tn; // ma->w = w_cast((ta->a[i].nw*sc)); ma->w = w_cast((ta->a[i].nw)); } for (i = k = 0; i < mg->e->ma.n; i++) { if(mg->e->ma.a[i].w == 0) continue; mg->e->ma.a[k] = mg->e->ma.a[i]; k++; } mg->e->ma.n = k; radix_sort_mce(mg->e->ma.a, mg->e->ma.a + mg->e->ma.n); mc_merge_dup(mg); mc_edges_idx(mg->e); mc_edges_symm(mg->e); kv_init(mg->s); mg->s.m = mg->s.n = s->n; CALLOC(mg->s.a, mg->s.n); for (i = 0; i < mg->s.n; ++i) { if(s->a[i].s != 1 && s->a[i].s != 2) continue; mg->s.a[i] = (s->a[i].s == 1? 1:-1); } return mg; } void debug_mc_gg_t(const char* fn, uint32_t update_ta, uint32_t convert_mc_g_t) { kv_u_trans_t *ta = NULL; kv_gg_status *s = NULL; uint32_t un; uint16_t hapN; if(load_mc_gg_dump(&ta, &un, &s, &hapN, fn)) { if(convert_mc_g_t) { mc_opt_t opt; mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed); mc_g_t *mg = NULL; mg = to_mc_g_t(ta, s, un); char *o_file = get_outfile_name(asm_opt.output_file_name); trans_chain* t_ch = load_hc_trans(o_file); mb_solve_core(&opt, mg, &(t_ch->k_trans), 1); mc_solve_core(&opt, mg, NULL); } else { mc_solve_general(ta, un, s, hapN, update_ta, 0); } } exit(1); } void clean_solve_general_ovlp(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s) { uint32_t i; for (i = 0; i < ta->n; i++) ta->a[i].del = !!(s->a[ta->a[i].qn].s&s->a[ta->a[i].tn].s); kt_u_trans_t_simple_symm(ta, un, 0); } void mc_solve_general(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s, uint16_t hapN, uint16_t update_ta, uint16_t write_dump) { print_hap_p(s); exit(1); mc_opt_t opt; mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed); mc_gg_t *mg = init_mc_gg_t(un, s, hapN); // print_mcb(mg); update_mc_edges_general(mg, ta, hapN); mc_solve_core_genral(&opt, mg, hapN); if(update_ta) clean_solve_general_ovlp(ta, un, s); // print_hap_p(s); destory_mc_gg_t(&mg); if(write_dump) write_mc_gg_dump(ta, un, s, hapN, MC_NAME); } void dump_kv_u_trans_t(kv_u_trans_t *z, FILE *fp) { fwrite(&(z->n), sizeof(z->n), 1, fp); fwrite(z->a, sizeof((*(z->a))), z->n, fp); fwrite(&z->idx.n, sizeof(z->idx.n), 1, fp); fwrite(z->idx.a, sizeof((*(z->idx.a))), z->idx.n, fp); } void dump_asg_t(asg_t *g, FILE *fp) { uint32_t tmp, Len; tmp = g->n_arc; fwrite(&tmp, sizeof(tmp), 1, fp); tmp = g->is_srt; fwrite(&tmp, sizeof(tmp), 1, fp); tmp = g->n_seq; fwrite(&tmp, sizeof(tmp), 1, fp); tmp = g->is_symm; fwrite(&tmp, sizeof(tmp), 1, fp); tmp = g->r_seq; fwrite(&tmp, sizeof(tmp), 1, fp); // Len = g->n_seq*2; // fwrite(g->seq_vis, sizeof((*(g->seq_vis))), Len, fp); Len = g->n_seq*2; fwrite(g->idx, sizeof((*(g->idx))), Len, fp); fwrite(g->arc, sizeof((*(g->arc))), g->n_arc, fp); fwrite(g->seq, sizeof((*(g->seq))), g->n_seq, fp); } void dump_ma_ug_t(ma_ug_t *ug, FILE *fp) { ma_utg_t *u = NULL; uint32_t t, i; fwrite(&(ug->u.n), sizeof(ug->u.n), 1, fp); for (i = 0; i < ug->u.n; i++) { u = &(ug->u.a[i]); t = u->len; fwrite(&t, sizeof(t), 1, fp); t = u->circ; fwrite(&t, sizeof(t), 1, fp); fwrite(&(u->start), sizeof(u->start), 1, fp); fwrite(&(u->end), sizeof(u->end), 1, fp); fwrite(&(u->n), sizeof(u->n), 1, fp); fwrite(u->a, sizeof(uint64_t), u->n, fp); } dump_asg_t(ug->g, fp); } void dump_bubble_type(bubble_type* bub, FILE *fp) { fwrite(&(bub->chain_weight.n), sizeof(bub->chain_weight.n), 1, fp); fwrite(bub->chain_weight.a, sizeof((*(bub->chain_weight.a))), bub->chain_weight.n, fp); fwrite(&(bub->list.n), sizeof(bub->list.n), 1, fp); fwrite(bub->list.a, sizeof((*(bub->list.a))), bub->list.n, fp); fwrite(&(bub->num.n), sizeof(bub->num.n), 1, fp); fwrite(bub->num.a, sizeof((*(bub->num.a))), bub->num.n, fp); fwrite(&(bub->pathLen.n), sizeof(bub->pathLen.n), 1, fp); fwrite(bub->pathLen.a, sizeof((*(bub->pathLen.a))), bub->pathLen.n, fp); dump_ma_ug_t(bub->b_ug, fp); dump_asg_t(bub->b_g, fp); dump_ma_ug_t(bub->ug, fp); } void dump_debug_phasing(const char* fn, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref) { fprintf(stderr, "\n[M::%s]\n", __func__); char *buf = (char*)calloc(strlen(fn) + 50, 1); sprintf(buf, "%s.hic.dbg.dump.bin", fn); FILE *fp = fopen(buf, "w"); dump_kv_u_trans_t(ta, fp);///ta dump_ma_ug_t(ug, fp);///ug dump_asg_t(read_g, fp);///read_g fwrite(&f_rate, sizeof(f_rate), 1, fp);///f_rate fwrite(&renew_s, sizeof(renew_s), 1, fp);///renew_s fwrite(s, sizeof((*s)), ug->g->n_seq, fp);///s fwrite(&is_sys, sizeof(is_sys), 1, fp);///is_sys dump_bubble_type(bub, fp);///bub dump_kv_u_trans_t(ref, fp);///ta fclose(fp); free(buf); } void load_kv_u_trans_t(kv_u_trans_t *z, FILE *fp) { fread(&(z->n), sizeof(z->n), 1, fp); z->m = z->n; MALLOC(z->a, z->n); fread(z->a, sizeof((*(z->a))), z->n, fp); fread(&z->idx.n, sizeof(z->idx.n), 1, fp); z->idx.m = z->idx.n; MALLOC(z->idx.a, z->idx.n); fread(z->idx.a, sizeof((*(z->idx.a))), z->idx.n, fp); } void load_asg_t(asg_t *g, FILE *fp) { uint32_t tmp, Len; fread(&tmp, sizeof(tmp), 1, fp); g->n_arc = g->m_arc = tmp; fread(&tmp, sizeof(tmp), 1, fp); g->is_srt = tmp; fread(&tmp, sizeof(tmp), 1, fp); g->n_seq = g->m_seq = tmp; fread(&tmp, sizeof(tmp), 1, fp); g->is_symm = tmp; fread(&tmp, sizeof(tmp), 1, fp); g->r_seq = tmp; // Len = g->n_seq*2; // fwrite(g->seq_vis, sizeof((*(g->seq_vis))), Len, fp); Len = g->n_seq*2; CALLOC(g->idx, Len); fread(g->idx, sizeof((*(g->idx))), Len, fp); CALLOC(g->arc, g->n_arc); fread(g->arc, sizeof((*(g->arc))), g->n_arc, fp); CALLOC(g->seq, g->n_seq); fread(g->seq, sizeof((*(g->seq))), g->n_seq, fp); } void load_ma_ug_t(ma_ug_t *ug, FILE *fp) { ma_utg_t *u = NULL; uint32_t t, i; fread(&(ug->u.n), sizeof(ug->u.n), 1, fp); ug->u.m = ug->u.n; CALLOC(ug->u.a, ug->u.n); for (i = 0; i < ug->u.n; i++) { u = &(ug->u.a[i]); fread(&t, sizeof(t), 1, fp); u->len = t; fread(&t, sizeof(t), 1, fp); u->circ = t; fread(&(u->start), sizeof(u->start), 1, fp); fread(&(u->end), sizeof(u->end), 1, fp); fread(&(u->n), sizeof(u->n), 1, fp); u->m = u->n; MALLOC(u->a, u->n); fread(u->a, sizeof((*(u->a))), u->n, fp); } CALLOC(ug->g, 1); load_asg_t(ug->g, fp); } void load_bubble_type(bubble_type* bub, FILE *fp) { fread(&(bub->chain_weight.n), sizeof(bub->chain_weight.n), 1, fp); bub->chain_weight.m = bub->chain_weight.n; CALLOC(bub->chain_weight.a, bub->chain_weight.n); fread(bub->chain_weight.a, sizeof((*(bub->chain_weight.a))), bub->chain_weight.n, fp); fread(&(bub->list.n), sizeof(bub->list.n), 1, fp); bub->list.m = bub->list.n; CALLOC(bub->list.a, bub->list.n); fread(bub->list.a, sizeof((*(bub->list.a))), bub->list.n, fp); fread(&(bub->num.n), sizeof(bub->num.n), 1, fp); bub->num.m = bub->num.n; CALLOC(bub->num.a, bub->num.n); fread(bub->num.a, sizeof((*(bub->num.a))), bub->num.n, fp); fread(&(bub->pathLen.n), sizeof(bub->pathLen.n), 1, fp); bub->pathLen.m = bub->pathLen.n; CALLOC(bub->pathLen.a, bub->pathLen.n); fread(bub->pathLen.a, sizeof((*(bub->pathLen.a))), bub->pathLen.n, fp); CALLOC(bub->b_ug, 1); load_ma_ug_t(bub->b_ug, fp); CALLOC(bub->b_g, 1); load_asg_t(bub->b_g, fp); CALLOC(bub->ug, 1); load_ma_ug_t(bub->ug, fp); } void load_debug_phasing(const char* fn, kv_u_trans_t **ta, ma_ug_t **ug, asg_t **read_g, double *f_rate, uint32_t *renew_s, int8_t **s, uint32_t *is_sys, bubble_type **bub, kv_u_trans_t **ref) { fprintf(stderr, "\n[M::%s]\n", __func__); char *buf = (char*)calloc(strlen(fn) + 50, 1); sprintf(buf, "%s.hic.dbg.dump.bin", fn); FILE *fp = fopen(buf, "r"); CALLOC((*ta), 1); load_kv_u_trans_t(*ta, fp);///ta CALLOC((*ug), 1); load_ma_ug_t(*ug, fp);///ug CALLOC((*read_g), 1); load_asg_t((*read_g), fp);///read_g fread(f_rate, sizeof((*f_rate)), 1, fp);///f_rate fread(renew_s, sizeof((*renew_s)), 1, fp);///renew_s CALLOC((*s), (*ug)->g->n_seq); fread(*s, sizeof((*(*s))), (*ug)->g->n_seq, fp);///s fread(is_sys, sizeof((*is_sys)), 1, fp);///is_sys CALLOC((*bub), 1); load_bubble_type(*bub, fp);///bub CALLOC((*ref), 1); load_kv_u_trans_t(*ref, fp);///ta fclose(fp); free(buf); } void quick_debug_phasing(const char* fn) { kv_u_trans_t *ta; ma_ug_t *ug; asg_t *read_g; double f_rate; uint32_t renew_s; int8_t *s; uint32_t is_sys, k; bubble_type *bub; kv_u_trans_t *ref; load_debug_phasing(fn, &ta, &ug, &read_g, &f_rate, &renew_s, &s, &is_sys, &bub, &ref); mc_solve(NULL, NULL, ta, ug, read_g, f_rate, NULL, renew_s, s, is_sys, bub, ref, 0, 0); for (k = 0; k < ug->g->n_seq; k++) { fprintf(stderr, "utg%.6dl(len::%u)\n", (int32_t)(k)+1, ug->g->seq[k].len); } exit(1); }