Skip to content

Commit

Permalink
Debug info
Browse files Browse the repository at this point in the history
  • Loading branch information
@eaasna
eaasna committed Sep 9, 2024
1 parent b7ea28d commit 1efd8d7
Show file tree
Hide file tree
Showing 5 changed files with 144 additions and 19 deletions.
6 changes: 6 additions & 0 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,12 @@ set (CMAKE_ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib")
set (CMAKE_LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib")
set (CMAKE_RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin")

# For debugging only
#set (CMAKE_CXX_FLAGS "-ftemplate-backtrace-limit=0")
#set (CMAKE_CXX_FLAGS "-fsanitize=address -g -O0")
set (CMAKE_CXX_FLAGS "-g -O0 -Wno-unused-parameter -Wno-unused-value -Wno-unused-but-set-variable -Wno-unused-variable -Wno-unused-local-typedefs")


# Messages
string (ASCII 27 Esc)
set (FontBold "${Esc}[1m")
Expand Down
65 changes: 59 additions & 6 deletions include/dream_stellar/verification/detail/all_or_best_local.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@
#include <dream_stellar/stellar_types.hpp>
#include <dream_stellar/utils/stellar_kernel_runtime.hpp>

#include <utilities/alphabet_wrapper/seqan/segment.hpp>

namespace dream_stellar
{

Expand All @@ -24,17 +26,30 @@ _appendNegativeSegment(TAlign const & align,
while (pos < len) {
if (isGap(row(align, 0), pos) || isGap(row(align, 1), pos)) {
score += scoreGap(scoreMatrix);
} else if (value(row(align, 0), pos) != value(row(align, 1), pos)) {
}
else{
std::cerr << "pos\t" << pos << '\n';
std::cerr << "seqan2::length(row(align, 0))\t" << seqan2::length(row(align, 0)) << '\n';
std::cerr << "seqan2::length(row(align, 1))\t" << seqan2::length(row(align, 1)) << '\n';

assert(seqan2::length(row(align, 0)) > pos);
assert(seqan2::length(row(align, 1)) > pos);
if (value(scoreMatrix, source(row(align, 0)), pos) != value(scoreMatrix, source(row(align, 1)), pos)) {
score += scoreMismatch(scoreMatrix);
} else {
break;
}
}

++pos;
}
if (pos == len) appendValue(queue, TMerger(beginPos, pos, MinValue<TScoreValue>::VALUE + 1));
else appendValue(queue, TMerger(beginPos, pos, score));
}

template <typename T>
class empty;

///////////////////////////////////////////////////////////////////////////////
// Appends a segment of only matching positions from align to queue.
template<typename TAlign, typename TPos, typename TScoreValue>
Expand All @@ -47,13 +62,51 @@ _appendPositiveSegment(TAlign const & align,
typedef Triple<TPos, TPos, TScoreValue> TMerger;
TPos beginPos = pos;

//!TODO: remove debugging
//std::cout << seqan3::to_char(value(row(align, 0), pos));
//std::cout << seqan3::to_char(value(row(align, 1), pos));

// type of row(align, 0)
// const seqan2::Gaps<seqan2::Segment<seqan2::Segment<const std::span<const seqan2::alphabet_adaptor<seqan3::dna4> >, seqan2::InfixSegment>, seqan2::InfixSegment>, seqan2::Tag<seqan2::ArrayGaps_> >&

// type of value(row(align, 0), 0)
// seqan2::Proxy<seqan2::IteratorProxy<seqan2::Iter<const seqan2::Gaps<TNestedSegment>, seqan2::Tag<seqan2::ArrayGaps_> >, seqan2::GapsIterator<seqan2::Tag<seqan2::ArrayGaps_> > > > >

// type of row(align, 0)._source
// seqan2::Holder<seqan2::Segment<seqan2::Segment<const std::span<const seqan2::alphabet_adaptor<seqan3::dna4> >, seqan2::InfixSegment>, seqan2::InfixSegment>, seqan2::Tag<seqan2::Tristate_> >

// type of row(align, 0)._source.data_value
// seqan2::Segment<seqan2::Segment<const std::span<const seqan2::alphabet_adaptor<seqan3::dna4> >, seqan2::InfixSegment>, seqan2::InfixSegment>*

// type of value(row(align, 0)._source.data_value, pos)
// seqan2::Segment<seqan2::Segment<const std::span<const seqan2::alphabet_adaptor<seqan3::dna4> >, seqan2::InfixSegment>, seqan2::InfixSegment>&

// type of source(row(align, 0))
// seqan2::Segment<seqan2::Segment<const std::span<const seqan2::alphabet_adaptor<seqan3::dna4> >, seqan2::InfixSegment>, seqan2::InfixSegment>&

// type of value(source(row(align, 0)), pos)
// const std::span<const seqan2::alphabet_adaptor<seqan3::dna4> >&

//empty<decltype(value(scoreMatrix, source(row(align, 0)), pos))> E{};
//if (value(source(row(align, 0)), pos) == value(source(row(align, 1)), pos))
// std::cout << "issame\n";

TScoreValue score = 0;
while ((pos < len) &&
(!isGap(row(align, 0), pos) &&
!isGap(row(align, 1), pos) &&
(value(row(align, 0), pos) == value(row(align, 1), pos)))) {
(value(scoreMatrix, source(row(align, 0)), pos) == value(scoreMatrix, source(row(align, 1)), pos)))) {
score += scoreMatch(scoreMatrix);
++pos;
if (pos < len)
{
std::cerr << "pos2\t" << pos << '\n';
std::cerr << "seqan2::length(row(align, 0))\t" << seqan2::length(row(align, 0)) << '\n';
std::cerr << "seqan2::length(row(align, 1))\t" << seqan2::length(row(align, 1)) << '\n';

assert(seqan2::length(row(align, 0)) > pos);
assert(seqan2::length(row(align, 1)) > pos);
}
}
appendValue(queue, TMerger(beginPos, pos, score));
}
Expand Down Expand Up @@ -127,6 +180,7 @@ _splitAtXDrops(TAlign const & align,
toViewPosition(row(align, 1), beginPosition(row(align, 1))));
appendValue(queue, TMerger(pos, pos, MinValue<TScoreValue1>::VALUE + 1));

//!TODO: limit aliLength do not go out of range over the underlying sequence
TPos aliLength = _max(toViewPosition(row(align, 0), endPosition(row(align, 0))),
toViewPosition(row(align, 1), endPosition(row(align, 1))));
TPos len;
Expand Down Expand Up @@ -247,18 +301,16 @@ allOrBestLocal(Segment<Segment<TSequence const, InfixSegment>, InfixSegment> con
if (!has_next)
break;

// while (localAlignment(localAlign, finder, scoreMatrix, minScore, lowerDiag, upperDiag, BandedWatermanEggert())) {

// std::cerr << "localAlign == \n" << localAlign << "\n";

// split local alignments containing an X-drop
String<Align<TSegment> > seedAlignments;
verification_runtime.split_at_x_drops_time.measure_time([&]()
{
_splitAtXDrops(localAlign, scoreMatrix, scoreDropOff, minScore, seedAlignments);
});

/*
typename Iterator<String<Align<TSegment> > >::Type aliIt = begin(seedAlignments);
while (aliIt != end(seedAlignments)) {
// std::cerr << "aliIt == \n" << row(*aliIt, 0) << "\n" << row(*aliIt, 1) << "\n";
// create alignment object for the complete sequences
Expand Down Expand Up @@ -292,6 +344,7 @@ allOrBestLocal(Segment<Segment<TSequence const, InfixSegment>, InfixSegment> con
++aliIt;
}
*/
if (bestLocalMethod) break;
}
}
Expand Down
17 changes: 15 additions & 2 deletions include/dream_stellar/verification/swift_hit_verifier.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ struct SwiftHitVerifier
{
static_assert(std::is_unsigned<TDelta>::value, "TDelta must be unsigned integral.");

/* try to convert std::span back to String
using TSegment = Segment<Segment<const std::span<TAlphabet const>, InfixSegment>, InfixSegment>;
TSegment finderSegment = databaseSegment.asFinderSegment();
TSegment patternSegment = querySegment.asPatternSegment();
Expand All @@ -41,10 +42,22 @@ struct SwiftHitVerifier
appendValue(host_infV, n);
Segment<const TString, InfixSegment> infV_seg{host_infV, 0, length(host_infV)};
Segment<Segment<const TString, InfixSegment>, InfixSegment> infV(infV_seg, seqan2::beginPosition(patternSegment), seqan2::endPosition(patternSegment));
*/

/*
std::cout << "Swift hit verifier\n";
for (auto n : databaseSegment.as_span())
std::cout << seqan3::to_char(n._symbol);
std::cout << '\n';
for (auto n : querySegment.as_span())
std::cout << seqan3::to_char(n._symbol);
std::cout << '\n';
*/

verifySwiftHit(
infH,
infV,
databaseSegment.asFinderSegment(),
querySegment.asPatternSegment(),
(double)eps_match_options.epsilon,
eps_match_options.minLength,
verifier_options.xDrop,
Expand Down
54 changes: 54 additions & 0 deletions include/utilities/alphabet_wrapper/seqan/segment.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
#pragma once

#include <seqan/sequence/segment_base.h>

namespace seqan2
{

template <typename TAlphabet, typename TPos>
TAlphabet value(seqan2::Segment<seqan2::Segment<std::span<TAlphabet const> const, InfixSegment>, InfixSegment> const & me,
TPos pos)
{
assert(pos >= 0);
auto const & host_span = *me.data_host.data_host;
if (pos + seqan2::beginPosition(me) >= host_span.size())
{
std::cerr << "pos\t" << pos << '\n';
std::cerr << "seqan2::beginPosition(me)" << seqan2::beginPosition(me) << '\n';
std::cerr << "host_span.size()\t" << host_span.size() << '\n';
}

assert(me.data_host.data_host);
assert(pos + seqan2::beginPosition(me) < host_span.size());
return host_span[seqan2::beginPosition(me) + pos];
}

template <typename TScore, typename TAlphabet, typename TPos>
TAlphabet value(TScore const & /* Score matrix to avoid ambiguity in overloading */,
seqan2::Segment<seqan2::Segment<std::span<TAlphabet const> const, InfixSegment>, InfixSegment> const & me,
TPos pos)
{
assert(pos >= 0);
auto const & host_span = *me.data_host.data_host;
if (pos + seqan2::beginPosition(me) >= host_span.size())
{
std::cerr << "pos\t" << pos << '\n';
std::cerr << "seqan2::beginPosition(me)" << seqan2::beginPosition(me) << '\n';
std::cerr << "host_span.size()\t" << host_span.size() << '\n';
}

assert(me.data_host.data_host);
//!TODO: it seems that the span goes out of scope at some point
assert(pos + seqan2::beginPosition(me) < host_span.size());
return host_span[seqan2::beginPosition(me) + pos];
}

template <typename TScore, typename TAlphabet, typename TPosition>
inline TAlphabet sequenceEntryForScore(TScore const & /*scoringScheme*/,
Segment<Segment<std::span<TAlphabet const> const, InfixSegment>, InfixSegment> const & seq,
TPosition pos)
{
return value(seq, pos);
}

} // namespace seqan2
21 changes: 10 additions & 11 deletions include/valik/search/search_local.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,11 +252,6 @@ bool search_local(search_arguments & arguments, search_time_statistics & time_st
dream_stellar::query_id_map<alphabet_t> query_dict{records};
stellarThreadTime.input_queries_time.manual_timing(cart_input_queries_time);

/* Debug
for (auto & query : queries)
seqan3::debug_stream << "Query sequence\n" << query << '\n';
*/

dream_stellar::_writeMoreCalculatedParams(threadOptions, threadOptions.referenceLength, queries, thread_meta.text_out);

auto swift_index_time = stellarThreadTime.swift_index_construction_time.now();
Expand All @@ -283,12 +278,10 @@ bool search_local(search_arguments & arguments, search_time_statistics & time_st
threadOptions.segmentBegin,
threadOptions.segmentEnd);

/* it works :)
using TInfixSegment = seqan2::Segment<sequence_reference_t const, seqan2::InfixSegment>;
TInfixSegment seqan2_segment(database, threadOptions.segmentBegin, threadOptions.segmentEnd);
sequence_reference_t host_sequence = host(seqan2_segment);
*/
std::cout << "Forward database\n";
for (auto n : database_segment.as_span())
std::cout << seqan3::to_char(n._symbol);
std::cout << '\n';

stellarThreadTime.forward_strand_stellar_time.measure_time([&]()
{
Expand Down Expand Up @@ -353,6 +346,12 @@ bool search_local(search_arguments & arguments, search_time_statistics & time_st
reverse_database.size() - threadOptions.segmentEnd /*begin*/,
reverse_database.size() - threadOptions.segmentBegin /*end*/);

std::cout << "Reverse database\n";
for (auto n : database_segment.as_span())
std::cout << seqan3::to_char(n._symbol);
std::cout << '\n';


stellarThreadTime.reverse_complement_database_time.manual_timing(reverse_database_time);

stellarThreadTime.reverse_strand_stellar_time.measure_time([&]()
Expand Down

0 comments on commit 1efd8d7

Please sign in to comment.