<commit_before>//(c) 2016-2017 by Authors
//This file is a part of ABruijn program.
//Released under the BSD license (see LICENSE file)

#include "contig_extender.h"
#include "../repeat_graph/output_generator.h"
#include <cmath>

void ContigExtender::generateUnbranchingPaths()
{
	GraphProcessor proc(_graph, _asmSeqs);
	_unbranchingPaths = proc.getUnbranchingPaths();

	_edgeToPath.clear();
	for (auto& path : _unbranchingPaths)
	{
		if (path.id.strand())
		{
			Logger::get().debug() << "UPath " << path.id.signedId() 
							<< ": " << path.edgesStr();
		}

		for (auto& edge : path.path)
		{
			_edgeToPath[edge] = &path;
		}
	}

	Logger::get().debug() << "Final graph contiain " 
		<< _unbranchingPaths.size() / 2 << " egdes";
}


void ContigExtender::generateContigs()
{
	Logger::get().debug() << "Extending contigs into repeats";

	bool graphContinue = (bool)Config::get("extend_contigs_with_repeats");

	OutputGenerator outGen(_graph, _aligner, _readSeqs);
	auto coreSeqs = outGen.generatePathSequences(_unbranchingPaths);
	std::unordered_map<UnbranchingPath*, FastaRecord*> upathsSeqs;
	for (size_t i = 0; i < _unbranchingPaths.size(); ++i)
	{
		upathsSeqs[&_unbranchingPaths[i]] = &coreSeqs[i];
	}

	std::unordered_map<GraphEdge*, 
					   std::vector<const GraphAlignment*>> alnIndex;
	for (auto& aln : _aligner.getAlignments())
	{
		if (aln.size() > 1)
		{
			for (auto edge : aln)
			{
				alnIndex[edge.edge].push_back(&aln);
			}
		}
	}

	std::unordered_set<GraphEdge*> coveredRepeats;
	std::unordered_map<const GraphEdge*, bool> repeatDirections;
	auto canTraverse = [&repeatDirections] (const GraphEdge* edge)
	{
		//if (edge->isLooped() && edge->selfComplement) return false;
		return !repeatDirections.count(edge) || 
			   repeatDirections.at(edge);
	};

	typedef std::pair<GraphPath, std::string> PathAndSeq;
	auto extendPathRight =
		[this, &coveredRepeats, &repeatDirections, &upathsSeqs, 
		 &canTraverse, &alnIndex, graphContinue] 
	(UnbranchingPath& upath)
	{

		bool extendFwd = !upath.path.back()->nodeRight->outEdges.empty();
		if (!extendFwd) return PathAndSeq();

		//first, choose the longest aligned read from this edge
		int32_t maxExtension = 0;
		GraphAlignment bestAlignment;
		for (auto pathPtr : alnIndex[upath.path.back()])
		{
			const GraphAlignment& path = *pathPtr;
			for (size_t i = 0; i < path.size(); ++i)
			{
				if (path[i].edge == upath.path.back() &&
					i < path.size() - 1)
				{
					size_t j = i + 1;
					while (j < path.size() && path[j].edge->repetitive &&
						   canTraverse(path[j].edge)) ++j;
					if (j == i + 1) break;

					int32_t alnLen = path[j - 1].overlap.curEnd - 
									 path[i + 1].overlap.curBegin;
					if (alnLen > maxExtension)
					{
						maxExtension = alnLen;
						bestAlignment.clear();
						std::copy(path.begin() + i + 1, path.begin() + j,
								  std::back_inserter(bestAlignment));
					}
					break;
				}
			}
		}
		if (maxExtension == 0) return PathAndSeq();

		auto upathAln = this->asUpathAlignment(bestAlignment);
		auto lastUpath = upathAln.back().upath;
		int32_t overhang = upathsSeqs[lastUpath]->sequence.length() - 
						   upathAln.back().aln.back().overlap.curEnd + 
						   upathAln.back().aln.front().overlap.curBegin;
		bool lastIncomplete = overhang > (int)Config::get("max_separation");
		//Logger::get().debug() << "Ctg " << upath.id.signedId() <<
		//	" overhang " << overhang << " upath " << lastUpath->id.signedId();

		for (size_t i = 0; i < upathAln.size(); ++i)
		{
			//in case we don't extend using graph structure,
			//dont mark the last upath as covered if it is
			//incomplete
			if (i == upathAln.size() - 1 && 
				lastIncomplete && !graphContinue) break;

			for (auto& aln : upathAln[i].aln)
			{
				repeatDirections[aln.edge] = true;
				repeatDirections[_graph.complementEdge(aln.edge)] = false;
				coveredRepeats.insert(aln.edge);
				coveredRepeats.insert(_graph.complementEdge(aln.edge));
			}
		}

		//generate extension sequence
		std::string extendedSeq;
		if (lastIncomplete && graphContinue && !lastUpath->isLooped())
		{
			upathAln.pop_back();
		}
		if (!upathAln.empty())
		{
			FastaRecord::Id readId = bestAlignment.front().overlap.curId;
			int32_t readStart = upathAln.front().aln.front().overlap.curBegin;
			int32_t readEnd = upathAln.back().aln.back().overlap.curEnd;
			extendedSeq = _readSeqs.getSeq(readId)
				.substr(readStart, readEnd - readStart).str();
		}
		if (lastIncomplete && graphContinue && !lastUpath->isLooped())
		{
			extendedSeq += upathsSeqs[lastUpath]->sequence.str();
		}
		
		GraphPath extendedPath;
		for (auto& ualn : upathAln)
		{
			for (auto& edgeAln : ualn.aln) extendedPath.push_back(edgeAln.edge);
		}
		if (lastIncomplete && graphContinue && !lastUpath->isLooped())
		{
			for (auto& edge : lastUpath->path) extendedPath.push_back(edge);
		}<commit_msg>allow contig extension into looped graph edges<commit_after>//(c) 2016-2017 by Authors
//This file is a part of ABruijn program.
//Released under the BSD license (see LICENSE file)

#include "contig_extender.h"
#include "../repeat_graph/output_generator.h"
#include <cmath>

void ContigExtender::generateUnbranchingPaths()
{
	GraphProcessor proc(_graph, _asmSeqs);
	_unbranchingPaths = proc.getUnbranchingPaths();

	_edgeToPath.clear();
	for (auto& path : _unbranchingPaths)
	{
		if (path.id.strand())
		{
			Logger::get().debug() << "UPath " << path.id.signedId() 
							<< ": " << path.edgesStr();
		}

		for (auto& edge : path.path)
		{
			_edgeToPath[edge] = &path;
		}
	}

	Logger::get().debug() << "Final graph contiain " 
		<< _unbranchingPaths.size() / 2 << " egdes";
}


void ContigExtender::generateContigs()
{
	Logger::get().debug() << "Extending contigs into repeats";

	bool graphContinue = (bool)Config::get("extend_contigs_with_repeats");

	OutputGenerator outGen(_graph, _aligner, _readSeqs);
	auto coreSeqs = outGen.generatePathSequences(_unbranchingPaths);
	std::unordered_map<UnbranchingPath*, FastaRecord*> upathsSeqs;
	for (size_t i = 0; i < _unbranchingPaths.size(); ++i)
	{
		upathsSeqs[&_unbranchingPaths[i]] = &coreSeqs[i];
	}

	std::unordered_map<GraphEdge*, 
					   std::vector<const GraphAlignment*>> alnIndex;
	for (auto& aln : _aligner.getAlignments())
	{
		if (aln.size() > 1)
		{
			for (auto edge : aln)
			{
				alnIndex[edge.edge].push_back(&aln);
			}
		}
	}

	std::unordered_set<GraphEdge*> coveredRepeats;
	std::unordered_map<const GraphEdge*, bool> repeatDirections;
	auto canTraverse = [&repeatDirections] (const GraphEdge* edge)
	{
		//if (edge->isLooped() && edge->selfComplement) return false;
		return !repeatDirections.count(edge) || 
			   repeatDirections.at(edge);
	};

	typedef std::pair<GraphPath, std::string> PathAndSeq;
	auto extendPathRight =
		[this, &coveredRepeats, &repeatDirections, &upathsSeqs, 
		 &canTraverse, &alnIndex, graphContinue] 
	(UnbranchingPath& upath)
	{

		bool extendFwd = !upath.path.back()->nodeRight->outEdges.empty();
		if (!extendFwd) return PathAndSeq();

		//first, choose the longest aligned read from this edge
		int32_t maxExtension = 0;
		GraphAlignment bestAlignment;
		for (auto pathPtr : alnIndex[upath.path.back()])
		{
			const GraphAlignment& path = *pathPtr;
			for (size_t i = 0; i < path.size(); ++i)
			{
				if (path[i].edge == upath.path.back() &&
					i < path.size() - 1)
				{
					size_t j = i + 1;
					while (j < path.size() && path[j].edge->repetitive &&
						   canTraverse(path[j].edge)) ++j;
					if (j == i + 1) break;

					int32_t alnLen = path[j - 1].overlap.curEnd - 
									 path[i + 1].overlap.curBegin;
					if (alnLen > maxExtension)
					{
						maxExtension = alnLen;
						bestAlignment.clear();
						std::copy(path.begin() + i + 1, path.begin() + j,
								  std::back_inserter(bestAlignment));
					}
					break;
				}
			}
		}
		if (maxExtension == 0) return PathAndSeq();

		auto upathAln = this->asUpathAlignment(bestAlignment);
		auto lastUpath = upathAln.back().upath;
		int32_t overhang = upathsSeqs[lastUpath]->sequence.length() - 
						   upathAln.back().aln.back().overlap.curEnd + 
						   upathAln.back().aln.front().overlap.curBegin;
		bool lastIncomplete = overhang > (int)Config::get("max_separation");
		//Logger::get().debug() << "Ctg " << upath.id.signedId() <<
		//	" overhang " << overhang << " upath " << lastUpath->id.signedId();

		for (size_t i = 0; i < upathAln.size(); ++i)
		{
			//in case we don't extend using graph structure,
			//dont mark the last upath as covered if it is
			//incomplete
			if (i == upathAln.size() - 1 && 
				lastIncomplete && !graphContinue) break;

			for (auto& aln : upathAln[i].aln)
			{
				repeatDirections[aln.edge] = true;
				repeatDirections[_graph.complementEdge(aln.edge)] = false;
				coveredRepeats.insert(aln.edge);
				coveredRepeats.insert(_graph.complementEdge(aln.edge));
			}
		}

		//generate extension sequence
		std::string extendedSeq;
		if (lastIncomplete && graphContinue)
		{
			upathAln.pop_back();
		}
		if (!upathAln.empty())
		{
			FastaRecord::Id readId = bestAlignment.front().overlap.curId;
			int32_t readStart = upathAln.front().aln.front().overlap.curBegin;
			int32_t readEnd = upathAln.back().aln.back().overlap.curEnd;
			extendedSeq = _readSeqs.getSeq(readId)
				.substr(readStart, readEnd - readStart).str();
		}
		if (lastIncomplete && graphContinue)
		{
			extendedSeq += upathsSeqs[lastUpath]->sequence.str();
		}
		
		GraphPath extendedPath;
		for (auto& ualn : upathAln)
		{
			for (auto& edgeAln : ualn.aln) extendedPath.push_back(edgeAln.edge);
		}
		if (lastIncomplete && graphContinue)
		{
			for (auto& edge : lastUpath->path) extendedPath.push_back(edge);
		}
{'fraction_non_alphanumeric': 0.15610413857483663, 'fraction_numerical': 0.006638315527434913, 'mean_word_length': 3.7172211350293543, 'pattern_counts': {'":': 2, '<': 83, '<?xml version=': 0, '>': 46, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4648954', 'n_tokens_mistral': 4090, 'n_tokens_neox': 3481, 'n_words': 673}