Option to print all scores including the ones below the threshold

Author: hani-girgis

src/Aligner.cpp

@@ -62,7 +62,7 @@ Aligner::Aligner(DataGenerator *d, ITransformer *t, Block *a, string dlmIn,
 	}
 	featList.clear();
 
-	isLengthFilter = filter;
+	canReportAll = filter;
 
 	ssPtr = new stringstream();
 }
@@ -136,7 +136,7 @@ void Aligner::processBlockHelper() {
 			init = j + 1;
 		}
 		auto p1 = blockA->at(j);
-		alignSeqVsBlock(p1.first, p1.second, blockB, kTable, monoTable,init);
+		alignSeqVsBlock(p1.first, p1.second, blockB, kTable, monoTable, init);
 	}
 
 	// Pop the block and free its memory
@@ -147,7 +147,7 @@ void Aligner::processBlockHelper() {
 template<class V>
 void Aligner::alignSeqVsBlock(std::string *info1, std::string *seq1,
 		Block *blockB, KmerHistogram<uint64_t, V> &kTable,
-		KmerHistogram<uint64_t, uint64_t> &monoTable,int init) {
+		KmerHistogram<uint64_t, uint64_t> &monoTable, int init) {
 
 	V *h1 = kTable.build(seq1);
 	uint64_t *mono1 = monoTable.build(seq1);
@@ -160,21 +160,23 @@ void Aligner::alignSeqVsBlock(std::string *info1, std::string *seq1,
 		auto p2 = blockB->at(hani);
 		string *seq2 = p2.second;
 		int l2 = seq2->size();
-		if (isLengthFilter
+
+		if (!canReportAll
 				&& (std::min(l1, l2) / (double) std::max(l1, l2) < threshold)) {
 			continue;
 		}
 
 		V *h2 = kTable.build(seq2);
 		uint64_t *mono2 = monoTable.build(seq2);
+
 		Statistician<V> s(histogramSize, k, h1, h2, mono1, mono2,
 				compositionList, keyList);
 
 		s.calculate(funIndexArray, singleFeatNum, data);
 
 		double res = predictor.calculateIdentity(data);
 
-		if (res >= threshold - error) {
+		if (canReportAll || res >= threshold - error) {
 			canWrite = true;
 
 			if (res > 1.0) {

src/Aligner.h

@@ -59,7 +59,7 @@ class Aligner {
 
 	// If enabled aligner does not align two sequences if they
 	// can achieve the minimum identity score.
-	bool isLengthFilter;
+	bool canReportAll;
 
 	double threshold;
 
@@ -83,7 +83,6 @@ class Aligner {
 			KmerHistogram<uint64_t, V> &kTable,
 			KmerHistogram<uint64_t, uint64_t> &monoTable, /*uint8_t *keyList,*/
 			int init);
-
 public:
 	Aligner(DataGenerator*, ITransformer*, Block*, string, bool, double,
 			double e, uint8_t*);

src/AlignerParallel.cpp

@@ -31,7 +31,7 @@ AlignerParallel<V>::AlignerParallel(int kmer, int hSize, double t,
 	threshold = t;
 	// error = e;
 	relaxThreshold = t - e;
-	isLengthFilter = f;
+	canReportAll = f;
 	compositionList = compList;
 	dlm = d;
 	threadNum = tNum;
@@ -128,16 +128,18 @@ void AlignerParallel<V>::setBlockA(Block *block, bool isAllVsAll) {
 					((sizeA - i - 1) / threadNum) + 1);
 
 #pragma omp parallel for schedule(static) num_threads(threadNum)
-			for (int j = i + 1; j < sizeA - 1; j++) {
-				double minimum = lenList[i];
-				double maximum = lenList[j];
-				if (maximum < minimum) {
-					minimum = lenList[j];
-					maximum = lenList[i];
-				}
-
-				if (isLengthFilter && (minimum / maximum < threshold)) {
-					continue;
+			for (int j = i + 1; j < sizeA; j++) {
+
+				if (!canReportAll) {
+					double minimum = lenList[i];
+					double maximum = lenList[j];
+					if (maximum < minimum) {
+						minimum = lenList[j];
+						maximum = lenList[i];
+					}
+					if ((minimum / maximum < threshold)) {
+						continue;
+					}
 				}
 
 				Statistician < V
@@ -148,7 +150,7 @@ void AlignerParallel<V>::setBlockA(Block *block, bool isAllVsAll) {
 				s.calculate(funIndexArray, singleFeatNum, data);
 				double res = predictor.calculateIdentity(data);
 
-				if (res >= relaxThreshold) {
+				if (canReportAll || res >= relaxThreshold) {
 					printList->at(omp_get_thread_num())->push_back(
 							std::make_pair(infoList[j], res));
 				}
@@ -239,7 +241,7 @@ void AlignerParallel<V>::output(
 				res = 0.0;
 			}
 
-			out << *info1 << dlm << *p.first << dlm << std::setprecision(4)
+			out << *info1 << dlm << *p.first << dlm << std::setprecision(8)
 					<< res << std::endl;
 		}
 		delete v;
@@ -266,15 +268,16 @@ void AlignerParallel<V>::processBlockB(Block *block) {
 #pragma omp parallel for schedule(static) num_threads(threadNum)
 		for (int h = 0; h < sizeB; h++) {
 
-			double minimum = lenList[i];
-			double maximum = lenListB[h];
-			if (maximum < minimum) {
-				minimum = lenListB[h];
-				maximum = lenList[i];
-			}
-
-			if (isLengthFilter && (minimum / maximum < threshold)) {
-				continue;
+			if (!canReportAll) {
+				double minimum = lenList[i];
+				double maximum = lenListB[h];
+				if (maximum < minimum) {
+					minimum = lenListB[h];
+					maximum = lenList[i];
+				}
+				if ((minimum / maximum < threshold)) {
+					continue;
+				}
 			}
 
 			Statistician < V
@@ -285,7 +288,7 @@ void AlignerParallel<V>::processBlockB(Block *block) {
 			s.calculate(funIndexArray, singleFeatNum, data);
 			double res = predictor.calculateIdentity(data);
 
-			if (res >= relaxThreshold) {
+			if (canReportAll || res >= relaxThreshold) {
 				printList->at(omp_get_thread_num())->push_back(
 						std::make_pair(infoListB[h], res));
 			}
@@ -309,7 +312,7 @@ void AlignerParallel<V>::processBlockB(Block *block) {
 		}
 	}
 
-	// Wait for the last print case
+// Wait for the last print case
 	if (printTask.valid()) {
 		printTask.get();
 		threadNum++;

src/AlignerParallel.h

@@ -59,7 +59,7 @@ class AlignerParallel {
 	int featNum;
 	int singleFeatNum;
 
-	bool isLengthFilter;
+	bool canReportAll;
 	double threshold;
 	double relaxThreshold;
 

src/BestFirst.cpp

@@ -1,17 +1,17 @@
 /*
-	Identity calculates DNA sequence identity scores rapidly without alignment.
+ Identity calculates DNA sequence identity scores rapidly without alignment.
 
-	Copyright (C) 2020 Hani Z. Girgis, PhD
+ Copyright (C) 2020 Hani Z. Girgis, PhD
 
-	Academic use: Affero General Public License version 1.
+ Academic use: Affero General Public License version 1.
 
-	Any restrictions to use for-profit or non-academics: Alternative commercial license is needed.
+ Any restrictions to use for-profit or non-academics: Alternative commercial license is needed.
 
-	This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
-	without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
-	Please contact Dr. Hani Z. Girgis (hzgirgis@buffalo.edu) if you need more information.
-*/
+ Please contact Dr. Hani Z. Girgis (hzgirgis@buffalo.edu) if you need more information.
+ */
 
 /*
  * BestFirst.cpp
@@ -63,13 +63,22 @@ BestFirst<T>::BestFirst(const Matrix &f, const Matrix &l,
 			best = v.first;
 			maximum = v.second;
 			eCount = 0;
+
+			// Print new features
+			std::cout << "Better performance of: " << maximum << std::endl;
+			const int *list = best.getList();
+			int s = best.getSize();
+			for (int a = 0; a < s; a++) {
+				std::cout << "\t" << fList.at(list[a])->getName() << std::endl;
+			}
 		}
 		// Step 5
 		std::vector < Node > eList = v.first.expand(fNum);
 		eCount++;
 
 		// Step 6
 		int childNum = eList.size();
+		// double resultList[childNum] { 0.0 };
 		std::vector<double> resultList(childNum, 0.0);
 #pragma omp parallel for schedule(static) num_threads(threadNum)
 		for (int i = 0; i < childNum; i++) {

src/DataGenerator.cpp

@@ -25,27 +25,15 @@
 /**
  * This constructor should be used for all vs. all
  */
-DataGenerator::DataGenerator(std::string fileName) {
-
-	FastaReader reader(fileName, Parameters::getBlockSize());
+DataGenerator::DataGenerator(std::string fileName, int blockSize) {
+	FastaReader reader(fileName, blockSize);
 	block = reader.read();
 
-	// Calculate average length
-	double sum = 0.0;
-	for (auto p : *block) {
-		sum += p.second->size();
-	}
-	int seqNum = block->size();
-	average = round(sum / seqNum);
-	cout << "Average: " << average << endl;
-
 	calculateK();
-
-	// Calculate histogram size
-	histogramSize = pow(4, k);
-	cout << "Histogram size: " << histogramSize << endl;
+	calculateHistSize();
 
 	// Find maximum length
+	int seqNum = block->size();
 	for (int h = 0; h < seqNum; h++) {
 		uint64_t len = block->at(h).second->size();
 		if (len > maxLength) {
@@ -128,34 +116,37 @@ DataGenerator::DataGenerator(std::string dbName, std::string qryFile,
 		maxLength = 2 * maxLength;
 	}
 
-	// Calculate average length
-	double sum = 0.0;
-	for (auto p : *block) {
-		sum += p.second->size();
-	}
-
-	average = round(sum / block->size());
-	cout << "Average: " << average << endl;
-
 	calculateK();
-
-	// Calculate histogram size
-	histogramSize = pow(4, k);
-	cout << "Histogram size: " << histogramSize << endl;
+	calculateHistSize();
 
 	// Free up memory
 	delete qryBlock;
 }
 
 void DataGenerator::calculateK() {
+	// Calculate average length
+	double sum = 0.0;
+	for (auto p : *block) {
+		sum += p.second->size();
+	}
+
+	average = round(sum / block->size());
+	std::cout << "Average: " << average << std::endl;
+	// Calculate k
 	k = ceil(log(average) / log(Parameters::getAlphabetSize()))
 			- Parameters::getKRelax();
 	if (k < 2) {
 		std::cerr << "DataGenerator warning: K is too small. ";
 		std::cerr << "K is set to 2." << std::endl;
 		k = 2;
 	}
-	cout << "K: " << k << endl;
+	std::cout << "K: " << k << std::endl;
+}
+
+void DataGenerator::calculateHistSize() {
+	// Calculate histogram size
+	histogramSize = pow(4, k);
+	std::cout << "Histogram size: " << histogramSize << std::endl;
 }
 
 DataGenerator::~DataGenerator() {

src/DataGenerator.h

@@ -1,17 +1,17 @@
 /*
-	Identity calculates DNA sequence identity scores rapidly without alignment.
+ Identity calculates DNA sequence identity scores rapidly without alignment.
 
-	Copyright (C) 2020 Hani Z. Girgis, PhD
+ Copyright (C) 2020 Hani Z. Girgis, PhD
 
-	Academic use: Affero General Public License version 1.
+ Academic use: Affero General Public License version 1.
 
-	Any restrictions to use for-profit or non-academics: Alternative commercial license is needed.
+ Any restrictions to use for-profit or non-academics: Alternative commercial license is needed.
 
-	This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
-	without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
-	Please contact Dr. Hani Z. Girgis (hzgirgis@buffalo.edu) if you need more information.
-*/
+ Please contact Dr. Hani Z. Girgis (hzgirgis@buffalo.edu) if you need more information.
+ */
 
 /*
  * DataGenerator.h
@@ -32,6 +32,7 @@
 class DataGenerator {
 private:
 	void calculateK();
+	void calculateHistSize();
 
 protected:
 	Block *block;
@@ -45,12 +46,12 @@ class DataGenerator {
 	double *compositionList;
 
 	// Model of all data (single statistics only)
-	Matrix *fTable;
+	Matrix *fTable = nullptr;
 	// Model of associated labels (identity scores)
-	Matrix *lTable;
+	Matrix *lTable = nullptr;
 
 public:
-	DataGenerator(std::string);
+	DataGenerator(std::string, int blockSize = Parameters::getBlockSize());
 	DataGenerator(std::string, std::string, double);
 	virtual ~DataGenerator();
 

src/FastaReader.cpp

@@ -70,6 +70,8 @@ FastaReader::FastaReader(std::string fileNameIn, int blockSizeIn,
 	codeMap['B'] = 'T';
 	codeMap['V'] = 'A';
 	codeMap['D'] = 'T';
+	// Added on 11/19/2020 to enable reading multi-sequence alignments
+	codeMap['-'] = '-';
 }
 
 FastaReader::~FastaReader() {