team18-ml-based-seeding

Makefile

@@ -4,8 +4,20 @@ CPPFLAGS=	-DHAVE_KALLOC -D__AMD_SPLIT_KERNELS__ # -Wno-unused-but-set-variable -
 CPPFLAGS+= 	$(if $(MAX_MICRO_BATCH),-DMAX_MICRO_BATCH=\($(MAX_MICRO_BATCH)\))
 INCLUDES=	-I .
 OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o \
-			lchain.o align.o hit.o seed.o map.o format.o pe.o esterr.o splitidx.o \
+			lchain.o align.o hit.o map.o format.o pe.o esterr.o splitidx.o \
 			ksw2_ll_sse.o
+
+
+ifeq ($(USE_GPU),yes)
+OBJS += gpu/rmi_seed.o 
+else
+OBJS += seed.o
+endif
+
+USE_GPU ?= no
+
+
+
 # PROG=		minimap2-zerobranch-debug
 # PROG=		minimap2-nobalance-debug
 PROG=		minimap2$(SUFFIX)

gpu/gpu_config.json

@@ -2,22 +2,22 @@
     "num_streams": 1,
     "//num_streams": "Must set to 1 for current implementation of mm2-gb",
     "min_n": 512, 
-    "//min_n": "queries with less anchors will be handled on cpu",
-    "long_seg_buffer_size": 100000000,
-    "//long_seg_buffer_size": "maximum number of anchors to fit in the aggregated long segment buffer. ",
-    "max_total_n": 500000000, 
-    "max_read": 500000,
+    "//min_n": "queries with less anchors will be handled on cpu (/100000000)",
+    "long_seg_buffer_size": 25000000, 
+    "//long_seg_buffer_size": "maximum number of anchors to fit in the aggregated long segment buffer(500000000,500000). ",
+    "max_total_n": 125000000, 
+    "max_read": 125000,
      "//max_total_n, max_read": "maximum number of anchors / reads to fit in one micro batch. Make sure this fits in the device memory.",
     "range_kernel": {
         "blockdim": 512,
         "cut_check_anchors": 10,
         "//cut_check_anchors": "Number of anchors to check to attemp a cut",
         "anchor_per_block": 32768,
-        "//anchor_per_block": "Number of anchors each block handle. Must be int * blockdim"
+        "//anchor_per_block": "Number of anchors each block handle. Must be int * blockdim(32768)"
     },
     "score_kernel": {
         "micro_batch": 4,
-        "//micro_batch": "Number of micro batches to aggregate into one long kernel. Make sure your host memory size is at least micro_batch * device mem size * 2",
+        "//micro_batch": "Number of micro batches(2) to aggregate into one long kernel. Make sure your host memory size is at least micro_batch * device mem size * 2",
         "mid_blockdim": 512,
         "short_griddim": 2688,
         "long_griddim": 144,

gpu/rmi_seed.cu

@@ -0,0 +1,193 @@
+#include "../mmpriv.h"
+#include "../kalloc.h"
+#include "../ksort.h"
+#include "rmi_seed.cuh"
+#include "hipify.cuh"
+
+
+void precompute_idx_data(const mm_idx_t *mi, const mm128_v *mv, const uint64_t ***cr_values, int **t_values, int n) {
+    *cr_values = (const uint64_t **)malloc(n * sizeof(uint64_t *));
+    if (*cr_values == NULL) {
+        fprintf(stderr, "Failed to allocate cr_values\n");
+        exit(EXIT_FAILURE);
+    }
+    *t_values = (int *)malloc(n * sizeof(int));
+    if (*t_values == NULL) {
+        fprintf(stderr, "Failed to allocate t_values\n");
+        free(*cr_values);
+        exit(EXIT_FAILURE);
+    }
+    for (int i = 0; i < n; ++i) {
+        (*cr_values)[i] = mm_idx_get(mi, mv->a[i].x >> 8, &(*t_values)[i]);
+    }
+}
+
+// GPUkernel
+__global__ void collect_seeds_kernel(const mm128_v *mv, mm_seed_t *m, int32_t *n_m, int32_t n, const uint64_t **cr_values, const int *t_values) {
+    int idx = blockIdx.x * blockDim.x + threadIdx.x;
+    if (idx >= n) return;
+
+    const uint64_t *cr = cr_values[idx];
+    int t = t_values[idx];
+    if (t == 0) return;
+
+    uint32_t q_pos = mv->a[idx].y, q_span = mv->a[idx].x & 0xff;
+
+    mm_seed_t *q = &m[idx];
+    q->q_pos = q_pos;
+    q->q_span = q_span;
+    q->cr = cr;
+    q->n = t;
+    q->seg_id = mv->a[idx].y >> 32;
+    q->is_tandem = q->flt = 0;
+
+    if (idx > 0 && mv->a[idx].x >> 8 == mv->a[idx - 1].x >> 8) q->is_tandem = 1;
+    if (idx < n - 1 && mv->a[idx].x >> 8 == mv->a[idx + 1].x >> 8) q->is_tandem = 1;
+
+    atomicAdd(n_m, 1);
+}
+
+
+mm_seed_t *mm_seed_collect_all(const mm_idx_t *mi, const mm128_v *mv, int32_t *n_m_) {
+    mm_seed_t *m_d;
+    int32_t *n_m_d, n_m = 0;
+    size_t n = mv->n;
+
+    const uint64_t **cr_values;
+    int *t_values;
+    precompute_idx_data(mi, mv, &cr_values, &t_values, n);
+
+    cudaMalloc((void**)&m_d, n * sizeof(mm_seed_t));
+    cudaMalloc((void**)&n_m_d, sizeof(int32_t));
+    cudaMemcpy(n_m_d, &n_m, sizeof(int32_t), cudaMemcpyHostToDevice);
+
+    const uint64_t **cr_values_d;
+    int *t_values_d;
+    cudaMalloc((void**)&cr_values_d, n * sizeof(uint64_t *));
+    cudaMalloc((void**)&t_values_d, n * sizeof(int));
+    cudaMemcpy(cr_values_d, cr_values, n * sizeof(uint64_t *), cudaMemcpyHostToDevice);
+    cudaMemcpy(t_values_d, t_values, n * sizeof(int), cudaMemcpyHostToDevice);
+
+    int blockSize = 256;
+    int numBlocks = (n + blockSize - 1) / blockSize;
+
+    collect_seeds_kernel<<<numBlocks, blockSize>>>(mv, m_d, n_m_d, n, cr_values_d, t_values_d);
+
+    cudaMemcpy(n_m_, n_m_d, sizeof(int32_t), cudaMemcpyDeviceToHost);
+    cudaFree(n_m_d);
+    cudaFree(cr_values_d);
+    cudaFree(t_values_d);
+
+    free(cr_values);
+    free(t_values);
+
+    return m_d;
+}
+
+
+
+
+void mm_seed_mz_flt(void *km, mm128_v *mv, int32_t q_occ_max, float q_occ_frac) {
+    mm128_t *a;
+    size_t i, j, st;
+    if (mv->n <= q_occ_max || q_occ_frac <= 0.0f || q_occ_max <= 0) return;
+    KMALLOC(km, a, mv->n);
+    for (i = 0; i < mv->n; ++i)
+        a[i].x = mv->a[i].x, a[i].y = i;
+    radix_sort_128x(a, a + mv->n);
+    for (st = 0, i = 1; i <= mv->n; ++i) {
+        if (i == mv->n || a[i].x != a[st].x) {
+            int32_t cnt = i - st;
+            if (cnt > q_occ_max && cnt > mv->n * q_occ_frac)
+                for (j = st; j < i; ++j)
+                    mv->a[a[j].y].x = 0;
+            st = i;
+        }
+    }
+    kfree(km, a);
+    for (i = j = 0; i < mv->n; ++i)
+        if (mv->a[i].x != 0)
+            mv->a[j++] = mv->a[i];
+    mv->n = j;
+}
+
+#define MAX_MAX_HIGH_OCC 128
+
+void mm_seed_select(int32_t n, mm_seed_t *a, int len, int max_occ, int max_max_occ, int dist) {
+    extern void ks_heapdown_uint64_t(size_t i, size_t n, uint64_t*);
+    extern void ks_heapmake_uint64_t(size_t n, uint64_t*);
+    int32_t i, last0, m;
+    uint64_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation
+
+    if (n == 0 || n == 1) return;
+    for (i = m = 0; i < n; ++i)
+        if (a[i].n > max_occ) ++m;
+    if (m == 0) return; // no high-frequency k-mers; do nothing
+    for (i = 0, last0 = -1; i <= n; ++i) {
+        if (i == n || a[i].n <= max_occ) {
+            if (i - last0 > 1) {
+                int32_t ps = last0 < 0? 0 : (uint32_t)a[last0].q_pos >> 1;
+                int32_t pe = i == n? len : (uint32_t)a[i].q_pos >> 1;
+                int32_t j, k, st = last0 + 1, en = i;
+                int32_t max_high_occ = (int32_t)((double)(pe - ps) / dist + .499);
+                if (max_high_occ > 0) {
+                    if (max_high_occ > MAX_MAX_HIGH_OCC)
+                        max_high_occ = MAX_MAX_HIGH_OCC;
+                    for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
+                        b[k] = (uint64_t)a[j].n << 32 | j;
+                    ks_heapmake_uint64_t(k, b); // initialize the binomial heap
+                    for (; j < en; ++j) { // if there are more, choose top max_high_occ
+                        if (a[j].n < (int32_t)(b[0] >> 32)) { // then update the heap
+                            b[0] = (uint64_t)a[j].n << 32 | j;
+                            ks_heapdown_uint64_t(0, k, b);
+                        }
+                    }
+                    for (j = 0; j < k; ++j) a[(uint32_t)b[j]].flt = 1;
+                }
+                for (j = st; j < en; ++j) a[j].flt ^= 1;
+                for (j = st; j < en; ++j)
+                    if (a[j].n > max_max_occ)
+                        a[j].flt = 1;
+            }
+            last0 = i;
+        }
+    }
+}
+
+mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int max_max_occ, int dist, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos) {
+    int rep_st = 0, rep_en = 0, n_m, n_m0;
+    size_t i;
+    mm_seed_t *m;
+    *n_mini_pos = 0;
+    *mini_pos = (uint64_t*)kmalloc(km, mv->n * sizeof(uint64_t));
+    m = mm_seed_collect_all(mi, mv, &n_m0);
+    if (dist > 0 && max_max_occ > max_occ) {
+        mm_seed_select(n_m0, m, qlen, max_occ, max_max_occ, dist);
+    } else {
+        for (i = 0; i < n_m0; ++i)
+            if (m[i].n > max_occ)
+                m[i].flt = 1;
+    }
+    for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < n_m0; ++i) {
+        mm_seed_t *q = &m[i];
+        // fprintf(stderr, "X\t%d\t%d\t%d\n", q->q_pos >> 1, q->n, q->flt);
+        if (q->flt) {
+            int en = (q->q_pos >> 1) + 1, st = en - q->q_span;
+            if (st > rep_en) {
+                *rep_len += rep_en - rep_st;
+                rep_st = st, rep_en = en;
+            } else rep_en = en;
+        } else {
+            *n_a += q->n;
+            (*mini_pos)[(*n_mini_pos)++] = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
+            m[n_m++] = *q;
+        }
+    }
+    *rep_len += rep_en - rep_st;
+    *_n_m = n_m;
+
+    // 释放在 mm_seed_collect_all 中分配的 GPU 内存
+    cudaFree(m);
+
+    return m;
+}

gpu/rmi_seed.cuh

@@ -0,0 +1,31 @@
+#ifndef RMI_SEED_H
+#define RMI_SEED_H
+
+#include "../mmpriv.h"
+#include "../ksort.h"  // 包含 ksort.h 以确保声明正确
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// Function to filter out minimizers based on their occurrence
+void mm_seed_mz_flt(void *km, mm128_v *mv, int32_t q_occ_max, float q_occ_frac);
+
+// Host function to launch the GPU kernel for collecting seeds
+mm_seed_t* mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv, int32_t *n_m_);
+
+// Function to select seeds based on their frequency
+void mm_seed_select(int32_t n, mm_seed_t *a, int len, int max_occ, int max_max_occ, int dist);
+
+// Function to collect matches from minimizers
+mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int max_max_occ, int dist, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos);
+
+// Declare ks_heapmake_uint64_t and ks_heapdown_uint64_t
+void ks_heapmake_uint64_t(size_t n, uint64_t* a);
+void ks_heapdown_uint64_t(size_t i, size_t n, uint64_t* a);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // RMI_SEED_H

results.json

@@ -0,0 +1 @@
+{"args":{"name":"COPY"},"ph":"M","pid":1,"name":"process_name","sort_index":0}