Merge branch 'v2'

.github/workflows/c-cpp.yml

@@ -0,0 +1,52 @@
+name: C/C++ CI
+
+on:
+  push:
+    branches: [ master, v2, feature/tests ]
+  pull_request:
+    branches: [ master, v2, feature/tests ]
+
+jobs:
+
+  linux-test:
+    name: Linux test
+    runs-on: ubuntu-latest
+    env: 
+      INPUT_DIR: ./test/virus
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: make
+      run: make
+
+    - name: build + all2all (default k)
+      run: |
+        ./kmer-db build -multisample-fasta ${INPUT_DIR}/sequences.list k18.db
+        ./kmer-db all2all k18.db k18.csv
+        cmp k18.csv ${INPUT_DIR}/k18.csv
+    
+    - name: build + all2all (k=24)
+      run: |
+        ./kmer-db build -k 24 -multisample-fasta ${INPUT_DIR}/sequences.list k24.db
+        ./kmer-db all2all k24.db k24.csv
+        cmp k24.csv ${INPUT_DIR}/k24.csv
+
+    - name: distance (default k)
+      run: | 
+        ./kmer-db distance jaccard min max cosine mash k18.csv
+        cmp k18.csv.jaccard ${INPUT_DIR}/k18.csv.jaccard
+        cmp k18.csv.min ${INPUT_DIR}/k18.csv.min
+        cmp k18.csv.max ${INPUT_DIR}/k18.csv.max
+        cmp k18.csv.cosine ${INPUT_DIR}/k18.csv.cosine
+        cmp k18.csv.mash ${INPUT_DIR}/k18.csv.mash
+      
+      
+  macos-build:
+    name: macOS test
+    runs-on: macOS-latest
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: make
+      run: make
+

.gitignore

@@ -14,3 +14,4 @@
 *.sqlite
 *.json
 *.suo
+/src/kmer_db.vcxproj.user

makefile

@@ -15,10 +15,20 @@ else
 	EXTRA_LIBS_DIR = ""
 endif
 
+OS_NAME := $(shell uname -s)
+ifeq ($(OS_NAME), Darwin)
+	OMP_FLAGS = -Xpreprocessor -fopenmp 
+	ABI_FLAGS = 
+else
+	OMP_FLAGS = -fopenmp
+	ABI_FLAGS = -fabi-version=6
+endif
+
+
 CC = g++
-CFLAGS	= -Wall -O3 -m64 -std=c++11 -fopenmp -pthread -mavx -I $(KMER_DB_LIBS_DIR) -I $(EXTRA_LIBS_DIR)
-CFLAGS_AVX2	= -Wall -O3 -m64 -std=c++11 -fopenmp -pthread -mavx2 -I $(KMER_DB_LIBS_DIR) -I $(EXTRA_LIBS_DIR)
-CLINK	= -lm -O3 -std=c++11 -lpthread -fopenmp -mavx -fabi-version=6 
+CFLAGS	= -Wall -O3 -m64 -std=c++11 $(OMP_FLAGS) -pthread -mavx -I $(KMER_DB_LIBS_DIR) -I $(EXTRA_LIBS_DIR)
+CFLAGS_AVX2	= -Wall -O3 -m64 -std=c++11 $(OMP_FLAGS) -pthread -mavx2 -I $(KMER_DB_LIBS_DIR) -I $(EXTRA_LIBS_DIR)
+CLINK	= -lm -O3 -std=c++11 -lpthread $(OMP_FLAGS) -mavx $(ABI_FLAGS) 
 
 OBJS := $(KMER_DB_MAIN_DIR)/analyzer.o \
 	$(KMER_DB_MAIN_DIR)/console.o \
@@ -36,7 +46,7 @@ OBJS := $(KMER_DB_MAIN_DIR)/analyzer.o \
 	$(KMER_DB_LIBS_DIR)/mmer.o 
 
 $(KMER_DB_MAIN_DIR)/parallel_sorter.o: $(KMER_DB_MAIN_DIR)/parallel_sorter.cpp
-	$(CC) -O3 -mavx -m64 -std=c++11 -pthread -fopenmp -c $< -o $@
+	$(CC) -O3 -mavx -m64 -std=c++11 -pthread $(OMP_FLAGS) -c $< -o $@
 
 ifeq ($(NO_AVX2),true)
 ## no avx2 support

src/console.cpp

@@ -355,7 +355,7 @@ int Console::runAllVsAll(const std::string& dbFilename, const std::string& simil
 	cout << "Calculating matrix of common k-mers...";
 	start = std::chrono::high_resolution_clock::now();
 	LowerTriangularMatrix<uint32_t> matrix;
-	calculator(*db, matrix);
+	calculator.all2all(*db, matrix);
 	dt = std::chrono::high_resolution_clock::now() - start;
 	cout << "OK (" << dt.count() << " seconds)" << endl;
 
@@ -446,7 +446,7 @@ int Console::runOneVsAll(const std::string& dbFilename, const std::string& singl
 	cout << "Calculating similarity vector...";
 	start = std::chrono::high_resolution_clock::now();
 	std::vector<uint32_t> sims;
-	calculator(db, queryKmers, queryKmersCount, sims);
+	calculator.one2all(db, queryKmers, queryKmersCount, sims);
 	dt = std::chrono::high_resolution_clock::now() - start;
 	cout << "OK (" << dt.count() << " seconds)" << endl;
 
@@ -528,7 +528,7 @@ int Console::runNewVsAll(const std::string& dbFilename, const std::string& multi
 			auto start = std::chrono::high_resolution_clock::now();
 
 			sims.clear();
-			calculator(db, task->kmers, task->kmersCount, sims);
+			calculator.one2all<true>(db, task->kmers, task->kmersCount, sims);
 			ofs << endl << task->sampleName << "," << task->kmersCount << ",";
 			std::copy(sims.begin(), sims.end(), ostream_iterator<uint32_t>(ofs, ","));
 
@@ -614,12 +614,16 @@ int Console::runDistanceCalculation(const std::string& similarityFilename, const
 		if ((i + 1) % 10 == 0) {
 			cout << "\r" << i + 1 << "/" << kmersCount.size() << "...";
 		}
+		// find first comma 
+		auto pos = in.find(',');
+		string queryName(in.begin(), in.begin() + pos);
 
 		std::replace(in.begin(), in.end(), ',', ' ');
 		istringstream iss(in);
+		iss.seekg(pos + 1); // move right after first comma separator
 		uint64_t queryKmersCount = 0;
-		string queryName;
-		iss >> queryName >> queryKmersCount;
+
+		iss >> queryKmersCount;
 
 		auto newEnd = std::copy(std::istream_iterator<size_t>(iss), std::istream_iterator<size_t>(), intersections.begin());
 		size_t numVals = newEnd - intersections.begin();

src/kmc_api/kmer_defs.h

@@ -33,7 +33,7 @@
 
 
 	#include <stdio.h>
-	#include <ext/algorithm>
+	#include <algorithm>
 	#include <iostream>
 	using namespace std;
 

src/parallel_sorter.cpp

@@ -11,6 +11,8 @@ Authors: Sebastian Deorowicz, Adam Gudys, Maciej Dlugosz, Marek Kokot, Agnieszka
 
 #ifdef WIN32
 #include <ppl.h>
+#elif defined __APPLE__
+#include <algorithm>
 #else
 #include <parallel/algorithm>
 #endif
@@ -24,6 +26,8 @@ void ParallelSort(kmer_t *arr, size_t arr_size)
 #ifdef WIN32
 	concurrency::parallel_sort(arr, arr + arr_size);
 	//std::stable_sort(samplePatterns.begin(), samplePatterns.end(), pid_comparer);
+#elif defined __APPLE__
+	std:: stable_sort(arr, arr + arr_size);
 #else
 	__gnu_parallel::sort(arr, arr + arr_size);
 #endif
@@ -39,6 +43,8 @@ void ParallelSort(pair<pattern_id_t, pattern_id_t*> *arr, size_t arr_size, pair<
 #ifdef WIN32
 	concurrency::parallel_sort(arr, arr + arr_size, pid_comparer);
 	//std::stable_sort(samplePatterns.begin(), samplePatterns.end(), pid_comparer);
+#elif defined __APPLE__
+	std:: stable_sort(arr, arr + arr_size, pid_comparer);
 #else
 	__gnu_parallel::sort(arr, arr + arr_size, pid_comparer);
 #endif
@@ -55,6 +61,8 @@ void ParallelSort(pair<kmer_or_pattern_t, pattern_id_t*> *arr, size_t arr_size,
 #ifdef WIN32
 	concurrency::parallel_sort(arr, arr + arr_size, pid_comparer);
 	//std::stable_sort(samplePatterns.begin(), samplePatterns.end(), pid_comparer);
+#elif defined __APPLE__
+	std:: stable_sort(arr, arr + arr_size, pid_comparer);
 #else
 	__gnu_parallel::sort(arr, arr + arr_size, pid_comparer);
 #endif

src/prefix_kmer_db.cpp

@@ -2,12 +2,7 @@
 
 #include "log.h"
 
-#ifdef WIN32
-#include <ppl.h>
-#else
-#include <parallel/algorithm>
-#endif
-#include <omp.h>
+#include <algorithm>
 
 #include <numeric>
 #include <cassert>

src/similarity_calculator.cpp

@@ -18,7 +18,7 @@ SimilarityCalculator::SimilarityCalculator(int _num_threads, size_t cacheBufferM
 
 // *****************************************************************************************
 //
-void SimilarityCalculator::operator()(PrefixKmerDb& db, LowerTriangularMatrix<uint32_t>& matrix) const 
+void SimilarityCalculator::all2all(PrefixKmerDb& db, LowerTriangularMatrix<uint32_t>& matrix) const 
 {
 	// get stuff from database
 	auto& patterns = db.getPatterns();
@@ -394,7 +394,8 @@ void SimilarityCalculator::operator()(PrefixKmerDb& db, LowerTriangularMatrix<ui
 
 // *****************************************************************************************
 //
-void  SimilarityCalculator::operator()(const PrefixKmerDb& db, const kmer_t* kmers, size_t kmersCount, std::vector<uint32_t>& similarities) const {
+template <>
+void  SimilarityCalculator::one2all<true>(const PrefixKmerDb& db, const kmer_t* kmers, size_t kmersCount, std::vector<uint32_t>& similarities) const {
 
 	// get stuff from database
 	const auto& patterns = db.getPatterns();
@@ -528,6 +529,121 @@ void  SimilarityCalculator::operator()(const PrefixKmerDb& db, const kmer_t* kme
 		});
 	}
 
+	dt = std::chrono::high_resolution_clock::now() - start;
+	LOG_VERBOSE << "Pattern unpacking time: " << dt.count() << endl;
+}
+
+// *****************************************************************************************
+//
+template <>
+void SimilarityCalculator::one2all<false>(const PrefixKmerDb& db, const kmer_t* kmers, size_t kmersCount, std::vector<uint32_t>& similarities) const
+{
+	// get stuff from database
+	const auto& patterns = db.getPatterns();
+	int samples_count = db.getSamplesCount();
+	const auto& hashtables = db.getHashtables();
+
+	similarities.resize(samples_count, 0);
+
+	std::unordered_map<pattern_id_t, int32_t> patterns2count;
+
+	std::chrono::duration<double> dt;
+	auto start = std::chrono::high_resolution_clock::now();
+
+	// iterate over kmers in analyzed sample
+	for (size_t i = 0; i < kmersCount; ++i) {
+
+		if (i + PREFETCH_DIST < kmersCount) {
+			kmer_t prefetch_kmer = kmers[i + PREFETCH_DIST];
+
+			kmer_t prefix = GET_PREFIX_SHIFTED(prefetch_kmer);
+			suffix_t suffix = GET_SUFFIX(prefetch_kmer);
+
+			hashtables[prefix].prefetch(suffix);
+		}
+
+		// check if kmer exists in a database
+		kmer_t kmer = kmers[i];
+		kmer_t prefix = GET_PREFIX_SHIFTED(kmer);
+		suffix_t suffix = GET_SUFFIX(kmer);
+
+		auto entry = hashtables[prefix].find(suffix);
+
+		if (entry != nullptr) {
+			auto pid = *entry;
+			const auto& pattern = patterns[pid];
+
+			if (pattern.get_num_kmers() == 0)
+				continue;
+
+			++patterns2count[pid];
+		}
+	}
+
+	std::vector<std::pair<pattern_id_t, int32_t>> patterns2countVector(patterns2count.size());
+	int i = 0;
+
+	uint64_t sum_pattern_lengths = 0;
+	for (auto &entry : patterns2count)
+	{
+		sum_pattern_lengths += patterns[entry.first].get_num_samples();
+		patterns2countVector[i++] = entry;	
+	}
+
+	patterns2count.clear();
+
+	dt = std::chrono::high_resolution_clock::now() - start;
+	LOG_VERBOSE << "Pattern listing time: " << dt.count() << endl;
+
+	start = std::chrono::high_resolution_clock::now();
+
+	std::vector<uint32_t> samples(samples_count);
+
+	size_t lo = 0;
+	size_t hi = patterns2countVector.size();
+
+	for (size_t id = lo; id < hi; ++id) {
+		if (id + 1 < hi)
+			_mm_prefetch((const char*)(patterns.data() + patterns2countVector[id + 1].first), _MM_HINT_T0);
+
+		auto pid = patterns2countVector[id].first;
+		const auto& pattern = patterns[pid];
+		int num_samples = pattern.get_num_samples();
+		int to_add = patterns2countVector[id].second;
+
+		uint32_t* out = samples.data() + pattern.get_num_samples(); // start from the end
+
+		int64_t current_id = pid;
+		while (current_id >= 0) {
+			const auto& cur = patterns[current_id];
+
+			out -= cur.get_num_local_samples();
+			cur.decodeSamples(out);
+
+			current_id = cur.get_parent_id();
+		}
+
+		auto *p = samples.data();
+
+		int i;
+		for (i = 0; i + 4 <= num_samples; i += 4)
+		{
+			similarities[*p++] += to_add;
+			similarities[*p++] += to_add;
+			similarities[*p++] += to_add;
+			similarities[*p++] += to_add;
+		}
+		num_samples -= i;
+
+		switch (num_samples)
+		{
+		case 3: similarities[*p++] += to_add;
+		case 2: similarities[*p++] += to_add;
+		case 1: similarities[*p++] += to_add;
+		}
+
+	}
+
 	dt = std::chrono::high_resolution_clock::now() - start;
 	LOG_VERBOSE << "Pattern unpacking time: " << dt.count() << endl;
 }

src/similarity_calculator.h

@@ -5,9 +5,10 @@ class SimilarityCalculator {
 public:
 	SimilarityCalculator(int _num_threads, size_t cacheBufferMb);
 
-	virtual void operator()(PrefixKmerDb& db, LowerTriangularMatrix<uint32_t>& matrix) const;
+	void all2all(PrefixKmerDb& db, LowerTriangularMatrix<uint32_t>& matrix) const;
 
-	virtual void operator()(const PrefixKmerDb& db, const kmer_t* kmers, size_t kmersCount, std::vector<uint32_t>& vector) const;
+	template <bool parallel = true>
+	void one2all(const PrefixKmerDb& db, const kmer_t* kmers, size_t kmersCount, std::vector<uint32_t>& vector) const;
 
 
 protected:

src/version.h

@@ -1,13 +1,16 @@
 #pragma once
 
-#define VERSION "1.7.5"
-#define DATE "13.02.2020"
+#define VERSION "1.7.6"
+#define DATE "30.03.2020"
 
 
 /* 
 
 Version history
 
+1.7.6 (31.03.2020)
+- Fixed bug in distance mode when sequence id contained spaces.
+
 1.7.5 (13.02.2020)
 - Some compilation warnings removed.
 - Fixed crash on samples with small k-mers count or very small filter values.