/* * nvbio * Copyright (c) 2011-2014, NVIDIA CORPORATION. All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: * * Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * * Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * Neither the name of the NVIDIA CORPORATION nor the * names of its contributors may be used to endorse or promote products * derived from this software without specific prior written permission. * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE * DISCLAIMED. IN NO EVENT SHALL NVIDIA CORPORATION BE LIABLE FOR ANY * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ #include #include #include #include namespace nvbio { namespace bowtie2 { namespace cuda { namespace { std::string mate_file_name(const std::string& file_name, const uint32 anchor) { const size_t dot = file_name.find('.'); if (dot == std::string::npos) return file_name + (anchor ? ".2" : ".1"); else { std::string output = file_name; output.insert( dot, (anchor ? ".2" : ".1") ); return output; } } } // anonymous namespace // clear the persisting files // void persist_clear(const std::string& file_name) { // clear mate 0 { const std::string mate_file = mate_file_name( file_name, 0u ); FILE* file = fopen( mate_file.c_str(), "w" ); fclose( file ); } // clear mate 1 { const std::string mate_file = mate_file_name( file_name, 0u ); FILE* file = fopen( mate_file.c_str(), "w" ); fclose( file ); } } // compute a set of hits // void persist_hits( const std::string& file_name, const char* name, const uint32 anchor, const uint32 count, const SeedHitDequeArray& hit_deques) { // check whether we need to persist anything if (file_name == "") return; const std::string mate_file = mate_file_name( file_name, anchor ); FILE* file = fopen( mate_file.c_str(), "a+" ); SeedHitDequeArrayHostStorage hit_deques_h = static_cast( hit_deques ); fprintf( file, "\n---------------------------------------------------------------------------------\n" ); fprintf( file, "%s\n", name ); fprintf( file, "---------------------------------------------------------------------------------\n" ); for (uint32 i = 0; i < count; ++i) fprintf( file, "cnt[%u] = %u\n", i, uint32(hit_deques_h.m_counts[i]) ); fprintf( file, "\n" ); for (uint32 i = 0; i < count; ++i) { const uint32 n_ranges = hit_deques_h.m_counts[i]; const SeedHit* hits = &hit_deques_h.m_hits[0] + hit_deques_h.m_index[i]; fprintf( file, "hits[%u] = [%u]{\n", i, n_ranges ); for (uint32 j = 0; j < n_ranges; ++j) { const SeedHit hit = hits[j]; fprintf( file, " range[%u] = { (%u, %u - %u), dir[%u], pos[%u], rc[%u] }\n", j, hit.get_range().x, hit.get_range().y, hit.get_range().y - hit.get_range().x, hit.get_indexdir(), hit.get_posinread(), hit.get_readtype() ); } fprintf( file, "}\n" ); } fprintf( file, "\n" ); fclose( file ); } // persist a set of reads // void persist_reads( const std::string& file_name, const char* name, const uint32 anchor, const uint32 count, const thrust::device_vector::iterator iterator) { // check whether we need to persist anything if (file_name == "") return; const std::string mate_file = mate_file_name( file_name, anchor ); FILE* file = fopen( mate_file.c_str(), "a+" ); thrust::host_vector hvec( count ); thrust::copy( iterator, iterator + count, hvec.begin() ); fprintf( file, "\n---------------------------------------------------------------------------------\n" ); fprintf( file, "%s\n", name ); fprintf( file, "---------------------------------------------------------------------------------\n" ); for (uint32 i = 0; i < count; ++i) fprintf( file, "read[%u] = %u\n", i, hvec[i] ); fprintf( file, "\n" ); fclose( file ); } // persist a set of selected hits // void persist_selection( const std::string& file_name, const char* name, const uint32 anchor, const uint32 read_count, const packed_read* read_infos_dptr, const uint32 n_multi, const uint32 hits_queue_size, const ReadHitsIndex& hits_index, const HitQueues& hits_queue) { // check whether we need to persist anything if (file_name == "") return; const std::string mate_file = mate_file_name( file_name, anchor ); FILE* file = fopen( mate_file.c_str(), "a+" ); fprintf( file, "\n---------------------------------------------------------------------------------\n" ); fprintf( file, "%s\n", name ); fprintf( file, "---------------------------------------------------------------------------------\n" ); thrust::host_vector loc_vec( hits_queue.loc ); thrust::host_vector seed_vec( hits_queue.seed ); if (n_multi > 1) { const uint32 link_stride = hits_index.m_stride; thrust::host_vector link_hvec( hits_index.m_links ); thrust::host_vector read_infos_hvec( read_count ); thrust::host_vector idx_hvec( read_count ); uint32* link_hptr = thrust::raw_pointer_cast( &link_hvec.front() ); uint32* read_infos_hptr = (uint32*)thrust::raw_pointer_cast( &read_infos_hvec.front() ); cudaDeviceSynchronize(); cudaMemcpy( read_infos_hptr, read_infos_dptr, read_count * sizeof(uint32), cudaMemcpyDeviceToHost ); cudaDeviceSynchronize(); // sort the reads so as to show everything in the same order all the times thrust::copy( thrust::make_counting_iterator(0u), thrust::make_counting_iterator(0u) + read_count, idx_hvec.begin() ); thrust::sort_by_key( read_infos_hvec.begin(), read_infos_hvec.end(), idx_hvec.begin() ); fprintf( file, "selection = [reads: %u - batch-size: %u - link-stride: %u]\n", read_count, n_multi, link_stride ); for (uint32 i = 0; i < read_count; ++i) { const packed_read *read_info = (packed_read *)&read_infos_hvec[i]; const uint32 idx = idx_hvec[i]; const uint32 n = link_hvec[idx]; fprintf( file, "read[%06u:%06u:%u] = [%02u]{", i, read_info->read_id, read_info->top_flag, n ); if (n) { strided_iterator link_vec( link_hptr + idx + link_stride, link_stride ); for (uint32 j = 0; j < n; ++j) { const packed_seed seed = seed_vec[ link_vec[j] ]; const uint32 loc = loc_vec[ link_vec[j] ]; fprintf( file, " seed[pos:%u,dir:%u,rc:%u,top:%u,loc:%u]\n", (uint32)seed.pos_in_read, (uint32)seed.index_dir, (uint32)seed.rc, (uint32)seed.top_flag, loc ); } } fprintf( file, "}\n" ); } } else { thrust::host_vector read_infos_hvec( read_count ); thrust::host_vector idx_hvec( read_count ); uint32* read_infos_hptr = (uint32*)thrust::raw_pointer_cast( &read_infos_hvec.front() ); cudaDeviceSynchronize(); cudaMemcpy( read_infos_hptr, read_infos_dptr, read_count * sizeof(uint32), cudaMemcpyDeviceToHost ); cudaDeviceSynchronize(); // sort the reads so as to show everything in the same order all the times thrust::copy( thrust::make_counting_iterator(0u), thrust::make_counting_iterator(0u) + read_count, idx_hvec.begin() ); thrust::sort_by_key( read_infos_hvec.begin(), read_infos_hvec.end(), idx_hvec.begin() ); fprintf( file, "selection = [reads: %u - batch-size: 1]\n", read_count ); for (uint32 i = 0; i < read_count; ++i) { const packed_read *read_info = (packed_read *)&read_infos_hvec[i]; const uint32 idx = idx_hvec[i]; fprintf( file, "read[%06u:%06u:%u] = [%02u]{", i, read_info->read_id, read_info->top_flag, 1u ); const packed_seed seed = seed_vec[ idx ]; const uint32 loc = loc_vec[ idx ]; fprintf( file, " seed[pos:%u,dir:%u,rc:%u,top:%u,loc:%u]\n", (uint32)seed.pos_in_read, (uint32)seed.index_dir, (uint32)seed.rc, (uint32)seed.top_flag, loc ); fprintf( file, "}\n" ); } } fclose( file ); } // persist a set of scores // void persist_scores( const std::string& file_name, const char* name, const uint32 anchor, const uint32 read_count, const uint32 n_multi, const uint32 hits_queue_size, const ScoringQueues& scoring_queues) { // check whether we need to persist anything if (file_name == "") return; const std::string mate_file = mate_file_name( file_name, anchor ); FILE* file = fopen( mate_file.c_str(), "a+" ); fprintf( file, "\n---------------------------------------------------------------------------------\n" ); fprintf( file, "%s\n", name ); fprintf( file, "---------------------------------------------------------------------------------\n" ); thrust::host_vector read_infos_hvec( scoring_queues.active_reads.in_queue ); thrust::host_vector score_hvec( scoring_queues.hits.score ); uint32* read_infos_hptr = (uint32*)thrust::raw_pointer_cast( &read_infos_hvec.front() ); int32* score_hptr = thrust::raw_pointer_cast( &score_hvec.front() ); if (n_multi > 1) { const uint32 link_stride = scoring_queues.hits_index.m_stride; thrust::host_vector link_hvec( scoring_queues.hits_index.m_links ); thrust::host_vector idx_hvec( read_count ); uint32* link_hptr = thrust::raw_pointer_cast( &link_hvec.front() ); // sort the reads so as to show everything in the same order all the times thrust::copy( thrust::make_counting_iterator(0u), thrust::make_counting_iterator(0u) + read_count, idx_hvec.begin() ); thrust::sort_by_key( read_infos_hptr, read_infos_hptr + read_count, idx_hvec.begin() ); fprintf( file, "scores = [reads: %u - batch-size: %u - link-stride: %u]\n", read_count, n_multi, link_stride ); for (uint32 i = 0; i < read_count; ++i) { const packed_read *read_info = (packed_read *)&read_infos_hvec[i]; const uint32 idx = idx_hvec[i]; const uint32 n = link_hvec[idx]; if (n) { strided_iterator link_vec( link_hptr + idx + link_stride, link_stride ); fprintf( file, "read[%06u:%06u:%u] = {\n", i, read_info->read_id, read_info->top_flag ); for (uint32 j = 0; j < n; ++j) { const uint32 link = link_vec[j]; fprintf( file, " score[%04u] = %d\n", j, score_hvec[link] ); } fprintf( file, "}\n" ); } } fprintf( file, "\n" ); // sort thrust::sort( score_hvec.begin(), score_hvec.end() ); // write sorted score list fprintf( file, "\n" ); for (uint32 i = 0; i < hits_queue_size; ++i) fprintf( file, "score[%u] = %d\n", i, score_hvec[i] ); fprintf( file, "\n" ); int64 sum = 0; for (uint32 i = 0; i < hits_queue_size; ++i) sum += score_hvec[i]; fprintf( file, "sum = %lld\n", sum ); // compute a crc const char* ptr = (const char*)score_hptr; const uint64 crc = crcCalc( ptr, sizeof(int32)*hits_queue_size ); fprintf( file, "crc = %llu\n", crc ); } else { thrust::host_vector idx_hvec( read_count ); // sort the reads so as to show everything in the same order all the times thrust::copy( thrust::make_counting_iterator(0u), thrust::make_counting_iterator(0u) + read_count, idx_hvec.begin() ); thrust::sort_by_key( read_infos_hptr, read_infos_hptr + read_count, idx_hvec.begin() ); for (uint32 i = 0; i < read_count; ++i) { const packed_read *read_info = (packed_read *)&read_infos_hvec[i]; const uint32 idx = idx_hvec[i]; const uint32 score = score_hvec[idx]; fprintf( file, "read[%06u:%06u:%u] : score[%d]\n", i, read_info->read_id, read_info->top_flag, score ); } fprintf( file, "\n" ); } fclose( file ); } } // namespace cuda } // namespace bowtie2 } // namespace nvbio