diff --git a/src/sort/common_data.h b/src/sort/common_data.h index cd22fe5..346e6b8 100644 --- a/src/sort/common_data.h +++ b/src/sort/common_data.h @@ -11,15 +11,7 @@ #include "sam_io.h" #include "sort_args.h" +#include "global_vars.h" #define PROGRAM_NAME "FastSort" -#define BAM_COMPRESS_RATIAO 5 // 大概5倍压缩比，用来粗略估计解压后需要多少内存空间 - - -namespace nsgv { - -// 全局变量 for bamsort -extern SortArg gSortArg; // 参数 -extern HeaderBuf gInHdr; // 输入文件的header -extern bool gIsBigEndian; -}; // namespace nsgv \ No newline at end of file +#define BAM_COMPRESS_RATIAO 5 // 大概5倍压缩比，用来粗略估计解压后需要多少内存空间 \ No newline at end of file diff --git a/src/sort/compress.h b/src/sort/compress.h deleted file mode 100644 index bf11f4a..0000000 --- a/src/sort/compress.h +++ /dev/null @@ -1,10 +0,0 @@ -/* - Description: 压缩相关的函数 - - Copyright : All right reserved by ICT - - Author : Zhang Zhonghai - Date : 2026/02/08 -*/ - -#pragma once \ No newline at end of file diff --git a/src/sort/decompress.h b/src/sort/decompress.h deleted file mode 100644 index dcf5379..0000000 --- a/src/sort/decompress.h +++ /dev/null @@ -1,10 +0,0 @@ -/* - Description: 解压相关的函数 - - Copyright : All right reserved by ICT - - Author : Zhang Zhonghai - Date : 2026/02/08 -*/ - -#pragma once \ No newline at end of file diff --git a/src/sort/global_vars.cpp b/src/sort/global_vars.cpp index 2bfb8dc..c02de48 100644 --- a/src/sort/global_vars.cpp +++ b/src/sort/global_vars.cpp @@ -1,5 +1,8 @@ #include "global_vars.h" +// 全局参数 namespace nsgv { - -}; \ No newline at end of file +SortArg gSortArg; // 参数 +HeaderBuf gInHdr; // 输入文件的header +bool gIsBigEndian; +}; // namespace nsgv \ No newline at end of file diff --git a/src/sort/phase_1.cpp b/src/sort/phase_1.cpp index e081359..015eb2c 100644 --- a/src/sort/phase_1.cpp +++ b/src/sort/phase_1.cpp @@ -32,44 +32,33 @@ std::string getFileNameWithoutExt(const std::string& filepath) { return filepath.substr(nameStart); } -void phase1Pipeline() { +std::string getFilePathWithoutExt(const std::string& filepath) { + // 1. 找到最后一个路径分隔符 + size_t sep = filepath.find_last_of("/\\"); + size_t nameStart = (sep == std::string::npos) ? 0 : sep + 1; -#if 1 - /* set up*/ - Phase1PipelineArg phase1Arg; - phase1Arg.numThread = nsgv::gSortArg.NUM_THREADS; - const size_t kReadBufSize = 4L * 1024 * 1024 * phase1Arg.numThread; // 平均每线程4M缓冲区，累加起来，用来读入文件（BAM/SAM）（相对解压之后的缓冲区，大小可以忽略） + // 2. 找到最后一个点（作为扩展名的分隔符） + size_t dot = filepath.find_last_of('.'); - phase1Arg.midFileNamePrefix = getFileNameWithoutExt(nsgv::gSortArg.OUTPUT_FILE) + "_sort_mid."; - phase1Arg.maxMemBytes = nsgv::gSortArg.MAX_MEM; - phase1Arg.uncompressBufBytes = nsgv::gSortArg.MAX_MEM; // 比最大内存参数小点 - for (int i = 0; i nameStart) { + return filepath.substr(0, dot); } - phase1Arg.singleThreadMemBytes = kReadBufSize * BAM_COMPRESS_RATIAO; - spdlog::info("max mem: {}, uncompress mem: {}, single thread mem: {}", nsgv::gSortArg.MAX_MEM, phase1Arg.uncompressBufBytes, phase1Arg.singleThreadMemBytes); - for (int i = 0; i < phase1Arg.READ_BUF_NUM; ++i) { - phase1Arg.readData[i].Resize(kReadBufSize); - }; + return filepath; +} - // 根据最大内存参数，计算初始化开辟的空间 - phase1Arg.uncompressData.Resize(phase1Arg.uncompressBufBytes); +void phase1Pipeline(Phase1PipelineArg& p) { PROF_G_BEG(mid_all); /* create threads */ pthread_t tidArr[3]; // 2-stage pipeline - pthread_create(&tidArr[0], NULL, phase1ReadFile, &phase1Arg); - pthread_create(&tidArr[1], NULL, phase1Uncompress, &phase1Arg); - pthread_create(&tidArr[2], NULL, phase1MemCopy, &phase1Arg); + pthread_create(&tidArr[0], NULL, phase1ReadFile, &p); + pthread_create(&tidArr[1], NULL, phase1Uncompress, &p); + pthread_create(&tidArr[2], NULL, phase1MemCopy, &p); for (int i = 0; i < 3; ++i) pthread_join(tidArr[i], NULL); - spdlog::info("all bams num: {}", phase1Arg.bamNum); + spdlog::info("all bams num: {}", p.bamNum); PROF_G_END(mid_all); - - - - -#endif } \ No newline at end of file diff --git a/src/sort/phase_1.h b/src/sort/phase_1.h index 26dfec3..c087ce7 100644 --- a/src/sort/phase_1.h +++ b/src/sort/phase_1.h @@ -202,6 +202,8 @@ struct Phase1PipelineArg { // 排序第一阶段，并行流水线执行程序 -void phase1Pipeline(); +void phase1Pipeline(Phase1PipelineArg &p); -std::string getFileNameWithoutExt(const std::string& filepath); \ No newline at end of file +std::string getFileNameWithoutExt(const std::string& filepath); + +std::string getFilePathWithoutExt(const std::string& filepath); \ No newline at end of file diff --git a/src/sort/phase_2.cpp b/src/sort/phase_2.cpp index e69de29..f753149 100644 --- a/src/sort/phase_2.cpp +++ b/src/sort/phase_2.cpp @@ -0,0 +1,14 @@ +/* + Description: 排序过程的第二阶段 + + Copyright : All right reserved by ICT + + Author : Zhang Zhonghai + Date : 2026/02/08 +*/ + +#include "phase_2.h" + +void phase2Pipeline(Phase2PipelineArg& p) { + +} \ No newline at end of file diff --git a/src/sort/phase_2.h b/src/sort/phase_2.h index 46c7f16..63fef8c 100644 --- a/src/sort/phase_2.h +++ b/src/sort/phase_2.h @@ -13,6 +13,8 @@ */ #pragma once +#include + struct Phase2File { FILE* fp; @@ -20,4 +22,23 @@ struct Phase2File { // 当前读入的buffer指针df // ReadBuffer -}; \ No newline at end of file +}; + +// 循环缓冲区 +struct CircularBuffer { + +}; + + +// 循环数组 +struct CircularArray { + +}; + +/* 第二阶段的多线程流水线参数 */ +struct Phase2PipelineArg { + + +}; + +void phase2Pipeline(Phase2PipelineArg &p); \ No newline at end of file diff --git a/src/sort/compress.cpp b/src/sort/phase_2_merge.cpp similarity index 100% rename from src/sort/compress.cpp rename to src/sort/phase_2_merge.cpp diff --git a/src/sort/decompress.cpp b/src/sort/phase_2_merge.h similarity index 100% rename from src/sort/decompress.cpp rename to src/sort/phase_2_merge.h diff --git a/src/sort/process_header.cpp b/src/sort/phase_2_read.cpp similarity index 100% rename from src/sort/process_header.cpp rename to src/sort/phase_2_read.cpp diff --git a/src/sort/phase_2_read.h b/src/sort/phase_2_read.h new file mode 100644 index 0000000..e69de29 diff --git a/src/sort/phase_2_write.cpp b/src/sort/phase_2_write.cpp new file mode 100644 index 0000000..e69de29 diff --git a/src/sort/phase_2_write.h b/src/sort/phase_2_write.h new file mode 100644 index 0000000..e69de29 diff --git a/src/sort/process_header.h b/src/sort/process_header.h deleted file mode 100644 index ac1a89e..0000000 --- a/src/sort/process_header.h +++ /dev/null @@ -1,10 +0,0 @@ -/* - Description: 处理bam/sam文件的header，可能需要解压 - - Copyright : All right reserved by ICT - - Author : Zhang Zhonghai - Date : 2026/02/08 -*/ - -#pragma once \ No newline at end of file diff --git a/src/sort/sort.cpp b/src/sort/sort.cpp index 3c735e4..93e14f9 100644 --- a/src/sort/sort.cpp +++ b/src/sort/sort.cpp @@ -16,647 +16,65 @@ #include #include +#include "common_data.h" #include "const_val.h" +#include "global_vars.h" +#include "phase_1.h" +#include "phase_1_read.h" +#include "phase_1_uncompress.h" +#include "phase_1_write.h" +#include "phase_2.h" #include "sam_io.h" #include "sort_args.h" -#include "sort_impl.h" #include "util/profiling.h" #include "util/yarn.h" -#include "process_header.h" -#include "phase_1.h" -#include "phase_2.h" -#include "compress.h" -#include "decompress.h" - using std::string; -#define BAM_BLOCK_SIZE 16L * 1024 * 1024 +#define PARAM_BLOCK_SIZE 16L * 1024 * 1024 -namespace nsgv { - -SortArg gSortArg; // 参数 -HeaderBuf gInHdr; // 输入文件的header -bool gIsBigEndian; -DataBuffer gMarginBuf; // 用于存放跨两个gz block的解压数据 - -}; // namespace nsgv - -struct BamSortData { - uint16_t flag; - uint16_t qnameLen; - uint16_t bamLen; - int32_t tid; - int64_t pos; - uint8_t *uDataPos; - char *qname; // pointer to qname -}; - - -/* 将bam文件内容读取到buf，解析buf中的gz block长度信息 */ -static size_t doFirstPipeReadFile(FirstPipeArg &p, DataBuffer &halfBlock, FILE *fpr) { - ReadBuffer &readData = p.readData[p.readOrder % p.READ_BUF_NUM]; - size_t readState = 0; - size_t curReadPos = 0; - int blockLen = 0; - int maxBlockLen = 0; - - readState = fread(readData.dataBuf, 1, readData.readBufSize, fpr); - if (readState == 0) { return 0; } - - readData.startAddrArr.clear(); - - /* 处理上一个不完整的block */ // 需要一个额外的空间来保存上一个不完整的block数据，因为readbuffer里的data和block都会用来解析，而且解析之后都会清空。 - if (halfBlock.readPos > 0) { // 上一轮有剩余 - if (halfBlock.readPos < BLOCK_HEADER_LENGTH) { // 上一轮剩余数据不满足解析block长度信息 - memcpy(&halfBlock.data[halfBlock.readPos], readData.dataBuf, BLOCK_HEADER_LENGTH - halfBlock.readPos); - halfBlock.curLen = unpackInt16(&halfBlock.data[16]) + 1; // 更新一下剩余block的真正长度 - // spdlog::info("last remain, blocklen: {}, last load: {}", halfBlock.curLen, halfBlock.readPos); - } - memcpy(readData.blockBuf, halfBlock.data, halfBlock.readPos); - curReadPos = halfBlock.curLen - halfBlock.readPos; // curlen保存上一个block的长度，readPos保存上一个block在上一次读取中的长度 - memcpy(&readData.blockBuf[halfBlock.readPos], readData.dataBuf, curReadPos); // 将不完整的block剩余数据拷贝到curBlock - readData.startAddrArr.push_back(readData.blockBuf); - } - /* 解析读入buf中的文件数据，计算包含的每个block的长度信息和起始地址 */ - while (curReadPos + BLOCK_HEADER_LENGTH <= readState) { /* 确保能解析block长度 */ - blockLen = unpackInt16(&readData.dataBuf[curReadPos + 16]) + 1; - if (blockLen > maxBlockLen) { maxBlockLen = blockLen; } - if (curReadPos + blockLen <= readState) { /* 完整的block数据在buf里 */ - readData.startAddrArr.push_back(&readData.dataBuf[curReadPos]); - curReadPos += blockLen; - } else { - break; /* 当前block数据不完整，一部分在还没读入的file数据里 */ - } - } - /* 如果buf中包含不完整的block数据，先保存一下，放到下一轮里去处理 */ - halfBlock.readPos = readState - curReadPos; - halfBlock.curLen = blockLen; - if (halfBlock.readPos > 0) { - memcpy(halfBlock.data, &readData.dataBuf[curReadPos], halfBlock.readPos); // 将不完整的block拷贝到halfBlock - } - - // spdlog::info("block num-1: {}", readData.startAddrArr.size()); - //spdlog::info("read order: {}, max block len: {}", p.readOrder, maxBlockLen); - - return readState; -} - - -/* FirstPipe step-1 读取文件线程 */ -static void *firstPipeReadFile(void *data) { - FirstPipeArg &p = *(FirstPipeArg *)data; - /* 1. set up */ - FILE *fpr = fopen(nsgv::gSortArg.INPUT_FILE.c_str(), "rb"); - parseSamHeader(fpr, nsgv::gInHdr); - size_t fileSize = 0; - DataBuffer halfBlock(SINGLE_BLOCK_SIZE); - - /* 2. do the work */ - while (true) { - // self dependency - yarn::DEPENDENCY_NOT_TO_BE(p.readSig, p.READ_BUF_NUM); - - PROF_G_BEG(read); - size_t readState = doFirstPipeReadFile(p, halfBlock, fpr); - PROF_G_END(read); - - if (readState == 0) { - yarn::SIGNAL_FINISH(p.readSig, p.readFinish); - break; - } - - // update self status - yarn::UPDATE_SIG_ORDER(p.readSig, p.readOrder); - - fileSize += readState; - } - - spdlog::info("read file order: {}, file size: {}", p.readOrder, fileSize); - /* 3. clean up */ - fclose(fpr); - return nullptr; -} - -static int parseBam(uint8_t* addr, bam1_t* b) { - bam1_core_t* c = &b->core; - int32_t block_len, ret, i; - uint32_t new_l_data; - uint8_t tmp[32], *x; - - b->l_data = 0; - memcpy(&block_len, addr, 4); - if (nsgv::gIsBigEndian) - ed_swap_4p(&block_len); - if (block_len < 32) - return -4; // block_len includes core data - x = addr + 4; - - c->tid = le_to_u32(x); - c->pos = le_to_i32(x + 4); - uint32_t x2 = le_to_u32(x + 8); - c->bin = x2 >> 16; - c->qual = x2 >> 8 & 0xff; - c->l_qname = x2 & 0xff; - c->l_extranul = (c->l_qname % 4 != 0) ? (4 - c->l_qname % 4) : 0; - uint32_t x3 = le_to_u32(x + 12); - c->flag = x3 >> 16; - c->n_cigar = x3 & 0xffff; - c->l_qseq = le_to_u32(x + 16); - c->mtid = le_to_u32(x + 20); - c->mpos = le_to_i32(x + 24); - c->isize = le_to_i32(x + 28); -#if 0 - new_l_data = block_len - 32 + c->l_extranul; - if (new_l_data > INT_MAX || c->l_qseq < 0 || c->l_qname < 1) - return -4; - if (((uint64_t)c->n_cigar << 2) + c->l_qname + c->l_extranul + (((uint64_t)c->l_qseq + 1) >> 1) + c->l_qseq > (uint64_t)new_l_data) - return -4; - if (realloc_bam_data(b, new_l_data) < 0) - return -4; - b->l_data = new_l_data; - - if (bgzf_read_small(fp, b->data, c->l_qname) != c->l_qname) - return -4; - if (b->data[c->l_qname - 1] != '\0') { // try to fix missing nul termination - if (fixup_missing_qname_nul(b) < 0) - return -4; - } - for (i = 0; i < c->l_extranul; ++i) b->data[c->l_qname + i] = '\0'; - c->l_qname += c->l_extranul; - if (b->l_data < c->l_qname || bgzf_read_small(fp, b->data + c->l_qname, b->l_data - c->l_qname) != b->l_data - c->l_qname) - return -4; - if (fp->is_be) - swap_data(c, b->l_data, b->data, 0); - if (bam_tag2cigar(b, 0, 0) < 0) - return -4; - - // TODO: consider making this conditional - if (c->n_cigar > 0) { // recompute "bin" and check CIGAR-qlen consistency - hts_pos_t rlen, qlen; - bam_cigar2rqlens(c->n_cigar, bam_get_cigar(b), &rlen, &qlen); - if ((b->core.flag & BAM_FUNMAP) || rlen == 0) - rlen = 1; - b->core.bin = hts_reg2bin(b->core.pos, b->core.pos + rlen, 14, 5); - // Sanity check for broken CIGAR alignments - if (c->l_qseq > 0 && !(c->flag & BAM_FUNMAP) && qlen != c->l_qseq) { - hts_log_error("CIGAR and query sequence lengths differ for %s", bam_get_qname(b)); - return -4; - } - } -#endif - return 4 + block_len; -} - -// multi-thread uncompress bam blocks -static void mtUncompressBlock(void *data, long idx, int tid) { - PROF_T_BEG(mem_copy); - - UncompressData& p = *(UncompressData*)data; - ReadBuffer &readData = *p.readDataPtr; - - auto &blockItemArr = p.blockItemArr[tid]; - auto &bamItemArr = p.bamItemArr[tid]; - - auto &blockItem = blockItemArr.add(); - uint8_t *block = readData.startAddrArr[idx]; - - size_t dlen = SINGLE_BLOCK_SIZE; // 65535 - int block_length = unpackInt16(&block[16]) + 1; - uint32_t crc = le_to_u32(block + block_length - 8); - int ret = bgzfUncompress(blockItem.data, &dlen, (Bytef *)block + BLOCK_HEADER_LENGTH, - block_length - BLOCK_HEADER_LENGTH, crc); - if (ret != 0) { - spdlog::error("uncompress error, block id: {}, len: {}, ret: {}", idx, block_length, ret); - exit(0); - } - blockItem.blockId = idx; - blockItem.blockLen = dlen; - uint32_t nextBamStart = 0; - uint32_t bamLen = 0; - uint32_t bamNum = 0; - - /* 解析每个bam */ - while (nextBamStart + 4 <= blockItem.blockLen) { - OneBam& bam = bamItemArr.add(); - - //bam.blockThread = &blockItemArr; - //bam.blockIdx = blockItemArr.curIdx - 1; - bam.offset = nextBamStart; - uint8_t *curAddr = &blockItem.data[nextBamStart]; -#if 0 - nextBamStart += parseBam(curAddr, &bam.b); - ++bamNum; -#else - memcpy(&bamLen, curAddr, 4); - curAddr += 4; - if (nsgv::gIsBigEndian) ed_swap_4p(&bamLen); - //bam.bamLen = bamLen; - bam.tid = le_to_u32(curAddr); - bam.pos = le_to_i32(curAddr + 4); - uint32_t x2 = le_to_u32(curAddr + 8); - bam.qnameLen = x2 & 0xff; - - // spdlog::info("bam len: {}, bam struct len: {}", bamLen, sizeof(OneBam)); - // spdlog::info("bam name: {}, bam name len: {}", string((char*)(blockItem.data + bam.offset + OneBam::QnameOffset), bam.qnameLen), bam.qnameLen); - nextBamStart += 4 + bamLen; - ++bamNum; -#endif - } - // spdlog::info("bam num: {}", bamNum); - if (nextBamStart != blockItem.blockLen) { - spdlog::error("Block content does not contain integer number of bam records!"); - exit(0); - } - blockItem.bamNum = bamNum; - blockItemArr.bamNum += bamNum; - - // 判断是否超过内存阈值，只要有一个超过阈值，就转入下一阶段，进行排序和合并，但是也得把这一轮数据解压处理完 - if (blockItemArr.curIdx * SINGLE_BLOCK_SIZE * BLOCK_MEM_FACTOR >= p.father->singleThreadBytes) { - p.father->uncompressBufFull = 1; - } else { - // spdlog::info("block arr: {}, single thread: {}", blockItemArr.blockArr.size() * SINGLE_BLOCK_SIZE, p.father->singleThreadBytes); - } - PROF_T_END(tid, mem_copy); -} - -// multi-thread uncompress bam blocks -static void mtUncompressBlockBatch(void* data, long idx, int tid) { - PROF_T_BEG(mem_copy); - UncompressData& p = *(UncompressData*)data; - ReadBuffer& readData = *p.readDataPtr; - - int startIdx = START_IDX(idx, nsgv::gSortArg.NUM_THREADS, readData.startAddrArr.size()); - int stopIdx = STOP_IDX(idx, nsgv::gSortArg.NUM_THREADS, readData.startAddrArr.size()); - - auto& blockItemArr = p.blockItemArr[tid]; - auto& bamItemArr = p.bamItemArr[tid]; - - blockItemArr.add(stopIdx - startIdx); - - for (int i = startIdx; i < stopIdx; ++i) { - auto& blockItem = blockItemArr.blockArr[blockItemArr.curIdx++]; - - uint8_t* block = readData.startAddrArr[i]; - - size_t dlen = SINGLE_BLOCK_SIZE; // 65535 - int block_length = unpackInt16(&block[16]) + 1; - uint32_t crc = le_to_u32(block + block_length - 8); - int ret = bgzfUncompress(blockItem.data, &dlen, (Bytef*)block + BLOCK_HEADER_LENGTH, block_length - BLOCK_HEADER_LENGTH, crc); - if (ret != 0) { - spdlog::error("uncompress error, block id: {}, len: {}, ret: {}", idx, block_length, ret); - exit(0); - } - blockItem.blockId = i; - blockItem.blockLen = dlen; - uint32_t nextBamStart = 0; - uint32_t bamLen = 0; - uint32_t bamNum = 0; - -#if 1 - /* 解析每个bam */ - while (nextBamStart + 4 <= blockItem.blockLen) { - bamItemArr.bamArr.push_back(OneBam()); - OneBam& bam = bamItemArr.bamArr.back(); - - //bam.addr = &blockItem.data[nextBamStart]; - uint8_t* curAddr = &blockItem.data[nextBamStart]; - memcpy(&bamLen, curAddr, 4); - curAddr += 4; - if (nsgv::gIsBigEndian) - ed_swap_4p(&bamLen); - // bam.bamLen = bamLen; - bam.tid = le_to_u32(curAddr); - bam.pos = le_to_i32(curAddr + 4); - uint32_t x2 = le_to_u32(curAddr + 8); - bam.qnameLen = x2 & 0xff; - //bam.qnameAddr = (char*)(curAddr + 32); - - // spdlog::info("bam len: {}, bam struct len: {}", bamLen, sizeof(OneBam)); - nextBamStart += 4 + bamLen; - ++bamNum; - } - if (nextBamStart != blockItem.blockLen) { - spdlog::error("Block content does not contain integer number of bam records!"); - exit(0); - } -#endif - blockItem.bamNum = bamNum; - - } - // 判断是否超过内存阈值 - if (blockItemArr.curIdx * SINGLE_BLOCK_SIZE >= p.father->singleThreadBytes) { - p.father->uncompressBufFull = 1; - } else { - // spdlog::info("block arr: {}, single thread: {}", blockItemArr.blockArr.size() * SINGLE_BLOCK_SIZE, p.father->singleThreadBytes); - } - PROF_T_END(tid, mem_copy); -} - -// 线程内排序 -static void mtInThreadSort(void* data, long idx, int tid) { - UncompressData& p = *(UncompressData*)data; - auto& bamItemArr = p.bamItemArr[tid]; - auto& arr = bamItemArr.bamArr; - // 先按照坐标排序 - -#if 0 - if (tid == 0) - spdlog::info("Before: {} {} {} {} {} {} {} {} {} {}", arr[0].pos, arr[1].pos, arr[2].pos, arr[3].pos, arr[4].pos, arr[5].pos, arr[6].pos, - arr[7].pos, arr[8].pos, arr[9].pos); -#endif - - if (nsgv::gSortArg.SORT_COORIDINATE) { - std::sort(arr.begin(), arr.begin() + bamItemArr.curIdx, [](const OneBam& b1, const OneBam& b2) { return b1.pos < b2.pos; }); - } else if (nsgv::gSortArg.QUERY_NAME_TYPE == nsmd::QueryNameType::PICARD) { - std::sort(arr.begin(), arr.begin() + bamItemArr.curIdx, [](const OneBam& b1, const OneBam& b2) { - int cmp = 0; -// strncmp((char*)(b1.blockThread->blockArr[b1.blockIdx].data + b1.offset + OneBam::QnameOffset), -// (char*)(b2.blockThread->blockArr[b2.blockIdx].data + b2.offset + OneBam::QnameOffset), std::min(b1.qnameLen, b2.qnameLen)); - if (cmp == 0) - return b1.qnameLen < b2.qnameLen; - return cmp < 0; - }); - } - -#if 0 - if (tid == 0) - spdlog::info("After: {} {} {} {} {} {} {} {} {} {}", arr[0].pos, arr[1].pos, arr[2].pos, arr[3].pos, arr[4].pos, arr[5].pos, arr[6].pos, - arr[7].pos, arr[8].pos, arr[9].pos); -#endif - -} - -// 多线程压缩 -static void mtCompressBlock(void* data, long idx, int tid) { - FirstPipeArg &p = *(FirstPipeArg *)data; - auto& barr = p.sortedBamArr; - auto& t = p.taskArr[idx]; - auto& buf = p.threadBuf[tid]; - - uint8_t block[SINGLE_BLOCK_SIZE] = {0}; - uint8_t compressBlock[SINGLE_BLOCK_SIZE] = {0}; - int curAddr = 0; - for (int i = t.idx; i < t.idx + t.num; ++i) { - // copy data - //memccpy(&block[curAddr], barr[i]->blockThread->blockArr[barr[i]->blockIdx].data + barr[i]->offset, 1, barr[i]->bamLen); - //curAddr += barr[i]->bamLen; - } - - size_t dlen = BGZF_MAX_BLOCK_SIZE; - // spdlog::info("len: {}", curAddr); - bgzfCompress(block, &dlen, compressBlock, curAddr, -1); - // 先放到自己线程内部的buf里 - - // 然后在压缩完成之后，检测当前全局压缩序号，如果当前序号的block在自己线程内，则复制过去 - -} - -/* 将gz block进行解压，并进行线程内排序 */ -static void doFirstPipeUncompress(FirstPipeArg &p, int finish = 0) { - // return; - PROF_G_BEG(uncompress); - ReadBuffer &readData = p.readData[p.uncompressOrder % p.READ_BUF_NUM]; - UncompressData &uncompressData = p.uncompressData[p.uncompressOrder % p.UNCOMPRESS_BUF_NUM]; - uncompressData.readDataPtr = &readData; - -#if 1 - kt_for(p.numThread, mtUncompressBlock, &uncompressData, readData.startAddrArr.size()); -#else - kt_for(p.numThread, mtUncompressBlockBatch, &uncompressData, nsgv::gSortArg.NUM_THREADS); -#endif - - PROF_G_END(uncompress); - - // 判断是否超过内存阈值，只要有一个超过阈值，就转入下一阶段，进行排序和合并 - if (p.uncompressBufFull == 1 || finish) { - - spdlog::info("buf full - {}", p.mergeOrder); - - // spdlog::info("max mem: {}, single thread mem: {}", nsgv::gSortArg.MAX_MEM, p.singleThreadBytes); - - // sort，排序的时候应该线程利用率很高，此时不需要跟其他操作如压缩等进行重叠了 - PROF_G_BEG(sort); - kt_for(p.numThread, mtInThreadSort, &uncompressData, nsgv::gSortArg.NUM_THREADS); - PROF_G_END(sort); - // merge -#if 0 - auto& bamArr = p.sortedBamArr; - BamHeap heap; - heap.Init(&uncompressData.bamItemArr); - bamArr.resize(heap.Size()); - auto& taskArr = p.taskArr; - taskArr.clear(); - - const OneBam* bam = nullptr; - uint64_t posAll = 0; - int i = 0; - ArrayInterval intv{0, 0}; - uint64_t blockBytes = 0; - PROF_G_BEG(merge); - while ((bam = heap.Pop()) != nullptr) { - // posAll += bam->pos; - //spdlog::info("pos: {}", bam->pos); - bamArr[i++] = bam; - blockBytes += bam->bamLen; - if (blockBytes >= SINGLE_BLOCK_SIZE) { - intv.num = i - 1 - intv.idx; - taskArr.push_back(intv); - intv.idx = i - 1; - blockBytes = bam->bamLen; - } - } - PROF_G_END(merge); - spdlog::info("task size: {} {}", taskArr.size(), bamArr.size()); - if (i - 1 > intv.idx) { - intv.num = i - 1 - intv.idx; - taskArr.push_back(intv); - } - - // spdlog::info("pos all: {}", posAll); -#endif - // compress，此时并行压缩的多线程利用率应该很高了 - PROF_G_BEG(compress); -#if 0 - kt_for(p.numThread, mtCompressBlock, &p, taskArr.size()); -#endif - - for (auto &blockItemArr : uncompressData.blockItemArr) { - p.bamNum += blockItemArr.bamNum; - } - - uncompressData.ResetBlockArr(); - uncompressData.ResetBamArr(); - p.uncompressBufFull = 0; - p.mergeOrder++; - for(auto &buf : p.threadBuf) { - buf.curLen = 0; - } - PROF_G_END(compress); - - // 压缩完之后应该写入中间文件，这个应该可以overlap - } - // 等处理完这些数据之后，进入新一轮的解压和归并排序压缩输出中间文件。 - - // auto &bam = uncompressData.bamItemArr[0].bamArr.back(); - // string qname(bam.qnameAddr, bam.qnameLen); - // spdlog::info("bam name:{}", qname); -} - -/* FirstPipe step-2 并行解压gz blocks*/ -static void *firstPipeUncompress(void *data) { - FirstPipeArg &p = *(FirstPipeArg *)data; - - /* 2. do the work */ - while (true) { - // previous dependency - yarn::DEPENDENCY_NOT_TO_BE(p.readSig, 0); - // self dependency - yarn::DEPENDENCY_NOT_TO_BE(p.uncompressSig, p.UNCOMPRESS_BUF_NUM); - - if (p.readFinish) { - while (p.uncompressOrder < p.readOrder) { - yarn::DEPENDENCY_NOT_TO_BE(p.uncompressSig, p.UNCOMPRESS_BUF_NUM); - doFirstPipeUncompress(p, 1); - yarn::UPDATE_SIG_ORDER(p.uncompressSig, p.uncompressOrder); - } - yarn::SIGNAL_FINISH(p.uncompressSig, p.uncompressFinish); - break; - } - - doFirstPipeUncompress(p); - - // update status - yarn::CONSUME_SIGNAL(p.readSig); - yarn::UPDATE_SIG_ORDER(p.uncompressSig, p.uncompressOrder); - } - - spdlog::info("uncompress order: {}", p.uncompressOrder); - return nullptr; -} - -/* 将并行解压的数据放到一起，解析，到容量阈值后，进行排序并输出到中间文件 */ -void dofirstPipeWrite(FirstPipeArg &p) { - PROF_G_BEG(sort); - UncompressData &uncompressData = p.uncompressData[p.writeOrder % p.UNCOMPRESS_BUF_NUM]; - - size_t blockNum = 0; - for (int i = 0; i < p.numThread; ++i) { - blockNum += uncompressData.blockItemArr[i].curIdx; - } - size_t bamNum = blockNum * 256; - - size_t curBlockBamMemSize = blockNum * SINGLE_BLOCK_SIZE + bamNum * 32; - /* 内存使用量达到阈值后，进行排序 */ - if (curBlockBamMemSize > nsgv::gSortArg.MAX_MEM / 2) { - //uncompressData.ResetBlockArr(); - //uncompressData.ResetBamArr(); - } - - //spdlog::info("block num: {}, bam num: {}, block size: {}, bam size: {}", blockNum, bamNum, blockBytes, bamNum * 32); - // spdlog::info("block num: {}, bam num: {}, block size: {}, bam size: {}", blockNum, bamNum, blockNum * SINGLE_BLOCK_SIZE, bamNum * 32); - PROF_G_END(sort); -} - -/* FirstPipe step-3 串行写入中间文件（已经压缩好） */ -static void *firstPipeWrite(void *data) { - FirstPipeArg &p = *(FirstPipeArg *)data; - - while (true) { - yarn::DEPENDENCY_NOT_TO_BE(p.uncompressSig, 0); - - if (p.uncompressFinish) { - spdlog::info("uncompress finish, cur sort order: {}", p.writeOrder); - while (p.writeOrder < p.uncompressOrder) { - dofirstPipeWrite(p); - p.writeOrder += 1; - } - /* 需要检测一下buf中是否还有数据，如果还有，则需要进行排序，可以不输出到中间文件，直接进行second-pipe的归并排序 */ - break; - } - - dofirstPipeWrite(p); - - p.writeOrder += 1; - - // update status - yarn::CONSUME_SIGNAL(p.uncompressSig); - } - - spdlog::info("merge sort order: {}", p.writeOrder); - return nullptr; -} - -/* 对bam文件进行排序第一阶段，线程内排序，线程间merge，输入到中间bam文件 */ -static void bamSortFirstPipe() { - /* set up*/ - FirstPipeArg firstPipeArg; - firstPipeArg.numThread = nsgv::gSortArg.NUM_THREADS; - firstPipeArg.singleThreadBytes = nsgv::gSortArg.MAX_MEM / firstPipeArg.numThread; - spdlog::info("max mem: {}, single thread mem: {}", nsgv::gSortArg.MAX_MEM, firstPipeArg.singleThreadBytes); - const size_t kReadBufSize = 1L * 1024 * 1024 * firstPipeArg.numThread; // 平均每线程1M缓冲区，累加起来，用来读入文件（BAM/SAM）（相对解压之后的缓冲区，大小可以忽略） - for (int i = 0; i arr.size()) { arr.resize(curIdx + num); @@ -122,19 +120,8 @@ struct FastVector { void Clear() { curIdx = 0; } }; -/* */ - - - -class ThreadBlockArr; /* for step-2 parallel uncompress gz blocks，还有缓存满了之后的线程内排序，用来排序的*/ struct OneBam { - // uint16_t bamLen = 0; - // uint16_t qnameLen = 0; - // uint32_t tid = 0; - // char *qnameAddr = 0; // qname的地址 - // uint64_t pos = 0; // mapping 位置 - // uint8_t *addr = 0; // 地址 static constexpr int QnameOffset = 36; // 距离该sam记录开头地址的偏移量，包含了4字节bam长度，32字节的bam core信息，后面紧跟着qname字符串 uint32_t qnameLen; // 序列名字长度 uint32_t wholeBamLen = 0; // 包含4字节bam长度的所有bam内容长度，uint16只适合二代 @@ -143,308 +130,5 @@ struct OneBam { uint64_t offset = 0; // 距离首地址的偏移量 }; - - -struct ThreadBamArr { - vector bamArr; - int curIdx = 0; // - OneBam &add() { - if (curIdx < bamArr.size()) - return bamArr[curIdx++]; - else { -#if 0 - bamArr.resize((bamArr.size() + 1) << 1); - return bamArr[curIdx++]; - -#else - bamArr.push_back(OneBam()); - curIdx++; - return bamArr.back(); -#endif - } - } - void clear() { curIdx = 0; } -}; - -// 解压后的一个block数据 -struct OneBlock { - uint8_t data[SINGLE_BLOCK_SIZE]; - uint32_t blockLen = 0; // 解压后的数据长度 - uint64_t blockId = 0; // 按照顺序排列的block ID - uint64_t bamNum = 0; // 解压后的bam数量 -}; - -struct BlockIdIdx { - uint64_t blockId = 0; - int blockArrIdx = 0; // 在block数组里的索引 - // 方法一：重载 operator> 供 std::greater 使用 - bool operator>(const BlockIdIdx& other) const { - return blockId > other.blockId; // blockId 小的优先级高（最小堆） - } -}; -// 每个线程一个解压block数组 -struct ThreadBlockArr { - vector blockArr; // 解压后的数据 - int curIdx = 0; // 当前解压数据对应的vector的索引 - uint64_t bamNum = 0; // 解压后的bam数量 - - OneBlock& add() { - if (curIdx < blockArr.size()) - return blockArr[curIdx++]; - else { -#if 0 - blockArr.resize((blockArr.size() + 1) << 1); - return blockArr[curIdx++]; - -#else - blockArr.push_back(OneBlock()); - curIdx++; - return blockArr.back(); -#endif - } - } - void add(int num) { - if (curIdx + num > blockArr.size()) { - blockArr.resize(curIdx + num); - } - } - - void clear() { - curIdx = 0; - bamNum = 0; - } -}; - -class FirstPipeArg; - -// 第一阶段的解压、排序、归并、压缩 -struct UncompressData { - vector blockItemArr; // 每个thread一个，用来保存解压后的block数据 - vector bamItemArr; // 每个thread一个，用来保存解压后的bam数据 - - ReadBuffer *readDataPtr = nullptr; // 读取数据的指针 - FirstPipeArg* father = nullptr; // 保存父参数的指针，可能更新一些状态 - - UncompressData() { } - UncompressData(int numThread) { Resize(numThread); } - UncompressData(int numThread, int vecInitSize) { Resize(numThread, vecInitSize); } - void Resize(int numThread) { - Resize(numThread, 128); - } - void Resize(int numThread, int vecInitSize) { - blockItemArr.resize(numThread); - bamItemArr.resize(numThread); - for (int i = 0; i < numThread; ++i) { - blockItemArr[i].blockArr.reserve(vecInitSize); - bamItemArr[i].bamArr.reserve(vecInitSize * 192); // 192是每个block平均包含的bam数量，单线程解压一个block的数据，平均会得到192条bam记录，应该作为一个参数或者宏定义 - } - } - void ResetBlockArr() { - for (int i = 0; i < blockItemArr.size(); ++i) { - blockItemArr[i].clear(); - } - } - void ResetBamArr() { - for (int i = 0; i < bamItemArr.size(); ++i) { - bamItemArr[i].clear(); - } - } -}; -/* block 排序堆 */ -struct BlockArrIdIdx { - int arrId = 0; - size_t arrIdx = 0; // 下一个待读入数据的idx - const OneBlock *block = nullptr; -}; - -struct BlockGreaterThan{ - bool operator()(const BlockArrIdIdx &a, const BlockArrIdIdx &b) const { - return a.block->blockId > b.block->blockId; - } -}; - -/* 用来排序 */ -struct BlockHeap { - vector *arr2d; - priority_queue, BlockGreaterThan> minHeap; - size_t popNum = 0; - - int Init(vector *_arr2d) { - arr2d = _arr2d; - if (arr2d == nullptr) { - return -1; - } - for (int i = 0; i < arr2d->size(); ++i) { - auto &v = (*arr2d)[i]; - if (v.curIdx > 0) { - minHeap.push({i, 1, &v.blockArr[0]}); - } - } - return 0; - } - - const OneBlock *Pop() { - const OneBlock *ret = nullptr; - if (!minHeap.empty()) { - auto minVal = minHeap.top(); - minHeap.pop(); - ++popNum; - ret = minVal.block; - auto &v = (*arr2d)[minVal.arrId]; - if (v.curIdx > minVal.arrIdx) { - minHeap.push({minVal.arrId, minVal.arrIdx + 1, &v.blockArr[minVal.arrIdx]}); - } - } - return ret; - } - - size_t AllBlockBytes() { - size_t bytes = 0; - if (arr2d != nullptr) { - for (auto &v : *arr2d) { - for (int i = 0; i < v.curIdx; ++i) { - bytes += v.blockArr[i].blockLen; - } - } - } - return bytes; - } - - size_t Size() { - size_t len = 0; - if (arr2d != nullptr) { - for (auto &v : *arr2d) { - len += v.curIdx; - } - } - return len - popNum; - } -}; - -/* bam 排序堆 */ -struct BamArrIdIdx { - int arrId = 0; - size_t arrIdx = 0; // 下一个待读入数据的idx - const OneBam* bam = nullptr; -}; - -struct BamGreaterThan { - bool operator()(const BamArrIdIdx& a, const BamArrIdIdx& b) const { return a.bam->pos > b.bam->pos; } -}; - -/* 用来排序 bam*/ -template -struct BamHeap { - vector* arr2d; - priority_queue, GreaterThan> minHeap; - size_t popNum = 0; - - int Init(vector* _arr2d) { - arr2d = _arr2d; - if (arr2d == nullptr) { - return -1; - } - for (int i = 0; i < arr2d->size(); ++i) { - auto& v = (*arr2d)[i]; - if (v.curIdx > 0) { - minHeap.push({i, 1, &v.bamArr[0]}); - } - } - return 0; - } - - const OneBam* Pop() { - const OneBam* ret = nullptr; - if (!minHeap.empty()) { - auto minVal = minHeap.top(); - minHeap.pop(); - ++popNum; - ret = minVal.bam; - auto& v = (*arr2d)[minVal.arrId]; - if (v.curIdx > minVal.arrIdx) { - minHeap.push({minVal.arrId, minVal.arrIdx + 1, &v.bamArr[minVal.arrIdx]}); - } - } - return ret; - } - - size_t Size() { - size_t len = 0; - if (arr2d != nullptr) { - for (auto& v : *arr2d) { - len += v.curIdx; - } - } - return len - popNum; - } -}; - -struct ArrayInterval { - int idx; - int num; -}; - -/* for step-3 serial merge blocks and sort them */ - -struct MergeSortData { - DataBuffer bamData; // 用来保存解压后的数据 - // BamPtrArr bamPtrArr; // 每个bam对应的解压数据，起始地址和长度 - MergeSortData() { - // bamPtrArr.bamArr.reserve(128); - } -}; - typedef FastVector BamArr; -typedef FastVector BamPtrArr; -typedef FastVector BlockArr; - - -/* 第一阶段的多线程流水线参数 */ -struct FirstPipeArg { - static const int READ_BUF_NUM = 2; // 读入的buf数量 - static const int UNCOMPRESS_BUF_NUM = 1; // 解压的buf数量 - static const int COMPRESS_BUF_NUM = 2; // merge之后，每个线程缓冲区的数量 - - int numThread = 0; // 线程数 - uint64_t singleThreadBytes = 0; // 单线程开辟的内存字节上限 - - // for test - uint64_t bamNum = 0; // 解压后的bam数量 - - uint64_t readOrder = 0; // 读取文件 - uint64_t uncompressOrder = 0; // 并行解压gz block, 包含排序（缓冲区满之后排序），以及合并之后的压缩 - uint64_t writeOrder = 0; // 串行合并解压后的blocks，并解析每个bam的长度，达到内存阈值后，并行排序 - uint64_t mergeOrder = 0; // 用来给中间文件编号 - - volatile int readFinish = 0; - volatile int uncompressFinish = 0; - volatile int uncompressBufFull = 0; - - yarn::lock_t *readSig; - yarn::lock_t *uncompressSig; - - ReadBuffer readData[READ_BUF_NUM]; // 用来读如数据，双缓冲 - UncompressData uncompressData[UNCOMPRESS_BUF_NUM]; // 每个线程内保留一些自己的数据，比如解压缩的数据 - UncompressBlockBuffer unCompblockDataBuf; // 所有线程共用一个，串行往这里添加解压后的block数据 - - vector sortedBamArr; - vector taskArr; - vector threadBuf; - - MergeSortData mergeSortData; - - FirstPipeArg() - { - readSig = yarn::NEW_LOCK(0); - uncompressSig = yarn::NEW_LOCK(0); - for (int i = 0; i < UNCOMPRESS_BUF_NUM; ++i) { - uncompressData[i].father = this; - } - } -}; - -/* 第二阶段的多线程流水线参数 */ -struct PipeSecondArg -{ - /* data */ -}; \ No newline at end of file +typedef FastVector BamPtrArr; \ No newline at end of file diff --git a/src/sort/sort_impl.cpp b/src/sort/sort_impl.cpp deleted file mode 100644 index 884cf8c..0000000 --- a/src/sort/sort_impl.cpp +++ /dev/null @@ -1,41 +0,0 @@ -#include -#include -#include -#include - -#include "sort_impl.h" - -int ks_radixsort(size_t n, OneBam* buf, const sam_hdr_t* h, int tid) { - if (tid == 0) spdlog::info("sort in thread, {} {}", tid, n); - - int curr = 0, ret = -1; - ssize_t i; - OneBam *buf_ar2[2], *bam_a, *bam_b; - uint64_t max_pos = 1; - uint32_t max_tid = 1, tid_bytes = 0, pos_bytes = 0, byte = 0; - uint32_t tid_shift_l, tid_shift_r; - int nref = sam_hdr_nref(h); - - //if (tid == 0) - //spdlog::info("Before: {} {} {} {} {} {} {} {} {} {}", buf[0].pos, buf[1].pos, buf[2].pos, buf[3].pos, buf[4].pos, buf[5].pos, buf[6].pos, - //buf[7].pos, buf[8].pos, buf[9].pos); - // spdlog::info("name: {}", std::string(buf[0].qnameAddr, buf[0].qnameLen)); - // spdlog::info("name addr: {}", (uint64_t)buf[0].qnameLen); - // std::sort(buf, buf + n, [](const OneBam& b1, const OneBam& b2) { return b1.pos < b2.pos; }); - std::sort(buf, buf + n, [](const OneBam& b1, const OneBam& b2) { - int cmp = 0; -// strncmp((char*)(b1.blockThread->blockArr[b1.blockIdx].data + b1.offset + OneBam::QnameOffset), -// (char*)(b2.blockThread->blockArr[b2.blockIdx].data + b2.offset + OneBam::QnameOffset), std::min(b1.qnameLen, b2.qnameLen)); - if (cmp == 0) - return b1.qnameLen < b2.qnameLen; - return cmp < 0; - }); - if (tid == 0) { - //std::string s1 = (char*)(buf[0].blockThread->blockArr[buf[0].blockIdx].data + buf[0].offset + OneBam::QnameOffset); - //std::string s2 = (char*)(buf[1].blockThread->blockArr[buf[1].blockIdx].data + buf[1].offset + OneBam::QnameOffset); - //spdlog::info("After: {} {} {} {} {} {} {} {} {} {}", s1, s2, buf[2].pos, buf[3].pos, buf[4].pos, buf[5].pos, buf[6].pos, buf[7].pos, - //buf[8].pos, buf[9].pos); - } - - return 0; -} \ No newline at end of file diff --git a/src/sort/sort_impl.h b/src/sort/sort_impl.h deleted file mode 100644 index c9da2d6..0000000 --- a/src/sort/sort_impl.h +++ /dev/null @@ -1,36 +0,0 @@ -#pragma once - -#include -#include - -#include "sort.h" - -// 好像不需要用到这些了？ - -// Struct which contains the sorting key for TemplateCoordinate sort. -struct TemplateCoordinateKey { - int tid1; - int tid2; - hts_pos_t pos1; - hts_pos_t pos2; - bool neg1; - bool neg2; - const char *library; - char *mid; - char *name; - bool is_upper_of_pair; -}; - -// Struct which contains the a record, and the pointer to the sort tag (if any) or -// a combined ref / position / strand. -// Used to speed up sorts (coordinate, by-tag, and template-coordinate). -struct BamSortTag { - bam1_t *bam_record; - union { - const uint8_t *tag; - uint8_t pos_tid[12]; - TemplateCoordinateKey *key; - } u; -}; - -int ks_radixsort(size_t n, OneBam* buf, const sam_hdr_t* h, int tid); \ No newline at end of file