pipe lines

src/dssp.cpp

@@ -27,6 +27,7 @@
 // Calculate DSSP-like secondary structure information
 
 #include "dssp.hpp"
+#include "queue.hpp"
 
 #include <deque>
 #include <iomanip>
@@ -40,6 +41,8 @@ using helix_type = dssp::helix_type;
 using helix_position_type = dssp::helix_position_type;
 using chain_break_type = dssp::chain_break_type;
 
+using queue_type = blocking_queue<std::tuple<uint32_t,uint32_t>>;
+
 // --------------------------------------------------------------------
 
 const double
@@ -765,35 +768,21 @@ double CalculateHBondEnergy(residue &inDonor, residue &inAcceptor)
 
 // --------------------------------------------------------------------
 
-void CalculateHBondEnergies(std::vector<residue> &inResidues)
+void CalculateHBondEnergies(std::vector<residue> &inResidues, queue_type &q1)
 {
-	if (cif::VERBOSE)
-		std::cerr << "calculating hbond energies" << std::endl;
-
-	// Calculate the HBond energies
-
-	std::unique_ptr<cif::progress_bar> progress;
-	if (cif::VERBOSE == 0)
-		progress.reset(new cif::progress_bar((inResidues.size() * (inResidues.size() - 1) / 2), "calculate hbond energies"));
-
-	for (uint32_t i = 0; i + 1 < inResidues.size(); ++i)
+	for (;;)
 	{
-		auto &ri = inResidues[i];
+		const auto &[i, j] = q1.pop();
 
-		for (uint32_t j = i + 1; j < inResidues.size(); ++j)
-		{
-			auto &rj = inResidues[j];
+		if (i == 0 and j == 0)
+			break;
 
-			if (distance_sq(ri.mCAlpha, rj.mCAlpha) < kMinimalCADistance * kMinimalCADistance)
-			{
-				CalculateHBondEnergy(ri, rj);
-				if (j != i + 1)
-					CalculateHBondEnergy(rj, ri);
-			}
+		auto &ri = inResidues[i];
+		auto &rj = inResidues[j];
 
-			if (progress)
-				progress->consumed(1);
-		}
+		CalculateHBondEnergy(ri, rj);
+		if (j != i + 1)
+			CalculateHBondEnergy(rj, ri);
 	}
 }
 
@@ -869,7 +858,7 @@ bool Linked(const bridge &a, const bridge &b)
 
 // --------------------------------------------------------------------
 
-void CalculateBetaSheets(std::vector<residue> &inResidues, statistics &stats)
+void CalculateBetaSheets(std::vector<residue> &inResidues, statistics &stats, std::vector<std::tuple<uint32_t, uint32_t>> &q)
 {
 	if (cif::VERBOSE)
 		std::cerr << "calculating beta sheets" << std::endl;
@@ -877,58 +866,50 @@ void CalculateBetaSheets(std::vector<residue> &inResidues, statistics &stats)
 	// Calculate Bridges
 	std::vector<bridge> bridges;
 
-	std::unique_ptr<cif::progress_bar> progress;
-	if (cif::VERBOSE == 0)
-		progress.reset(new cif::progress_bar(((inResidues.size() - 5) * (inResidues.size() - 4) / 2), "calculate beta sheets"));
-
-	for (uint32_t i = 1; i + 4 < inResidues.size(); ++i)
+	for (const auto &[i, j] : q)
 	{
 		auto &ri = inResidues[i];
+		auto &rj = inResidues[j];
 
-		for (uint32_t j = i + 3; j + 1 < inResidues.size(); ++j)
-		{
-			auto &rj = inResidues[j];
+		bridge_type type = TestBridge(ri, rj);
+		if (type == bridge_type::None)
+			continue;
 
-			bridge_type type = TestBridge(ri, rj);
-			if (type == bridge_type::None)
+		bool found = false;
+		for (bridge &bridge : bridges)
+		{
+			if (type != bridge.type or i != bridge.i.back() + 1)
 				continue;
 
-			bool found = false;
-			for (bridge &bridge : bridges)
+			if (type == bridge_type::Parallel and bridge.j.back() + 1 == j)
 			{
-				if (type != bridge.type or i != bridge.i.back() + 1)
-					continue;
-
-				if (type == bridge_type::Parallel and bridge.j.back() + 1 == j)
-				{
-					bridge.i.push_back(i);
-					bridge.j.push_back(j);
-					found = true;
-					break;
-				}
-
-				if (type == bridge_type::AntiParallel and bridge.j.front() - 1 == j)
-				{
-					bridge.i.push_back(i);
-					bridge.j.push_front(j);
-					found = true;
-					break;
-				}
+				bridge.i.push_back(i);
+				bridge.j.push_back(j);
+				found = true;
+				break;
 			}
 
-			if (not found)
+			if (type == bridge_type::AntiParallel and bridge.j.front() - 1 == j)
 			{
-				bridge bridge = {};
-
-				bridge.type = type;
 				bridge.i.push_back(i);
-				bridge.chainI = ri.mAsymID;
-				bridge.j.push_back(j);
-				bridge.chainJ = rj.mAsymID;
-
-				bridges.push_back(bridge);
+				bridge.j.push_front(j);
+				found = true;
+				break;
 			}
 		}
+
+		if (not found)
+		{
+			bridge bridge = {};
+
+			bridge.type = type;
+			bridge.i.push_back(i);
+			bridge.chainI = ri.mAsymID;
+			bridge.j.push_back(j);
+			bridge.chainJ = rj.mAsymID;
+
+			bridges.push_back(bridge);
+		}
 	}
 
 	// extend ladders
@@ -1563,11 +1544,49 @@ void DSSP_impl::calculateSecondaryStructure()
 		mSSBonds.emplace_back(&*r1, &*r2);
 	}
 
-	CalculateHBondEnergies(mResidues);
-	CalculateBetaSheets(mResidues, mStats);
+	// Calculate the HBond energies
+
+	queue_type q1, q2;
+	std::vector<std::tuple<uint32_t,uint32_t>> near;
+
+	std::thread hbond_thread(std::bind(&CalculateHBondEnergies, std::ref(mResidues), std::ref(q1)));
+
+	std::unique_ptr<cif::progress_bar> progress;
+	if (cif::VERBOSE == 0)
+		progress.reset(new cif::progress_bar((mResidues.size() * (mResidues.size() - 1) / 2), "calculate hbond energies"));
+
+	for (uint32_t i = 0; i + 1 < mResidues.size(); ++i)
+	{
+		auto &ri = mResidues[i];
+
+		for (uint32_t j = i + 1; j < mResidues.size(); ++j)
+		{
+			auto &rj = mResidues[j];
+
+			if (distance_sq(ri.mCAlpha, rj.mCAlpha) < kMinimalCADistance * kMinimalCADistance)
+			{
+				q1.push({ i, j });
+				near.emplace_back(i, j);
+			}
+
+			if (progress)
+				progress->consumed(1);
+		}
+	}
+
+	q1.push({0, 0});
+
+	hbond_thread.join();
+
+	std::thread bsheet_thread(std::bind(&CalculateBetaSheets, std::ref(mResidues), std::ref(mStats), std::ref(near)));
+
+	// CalculateHBondEnergies(mResidues);
+	// CalculateBetaSheets(mResidues, mStats);
 	CalculateAlphaHelices(mResidues, mStats);
 	CalculatePPHelices(mResidues, mStats, m_min_poly_proline_stretch_length);
 
+	bsheet_thread.join();
+
 	if (cif::VERBOSE > 1)
 	{
 		for (auto &r : mResidues)

src/queue.hpp

+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ * 
+ * Copyright (c) 2021 NKI/AVL, Netherlands Cancer Institute
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. 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.
+ * 
+ * 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 THE COPYRIGHT OWNER OR CONTRIBUTORS 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.
+ */
+
+#pragma once
+
+#include <condition_variable>
+#include <mutex>
+#include <queue>
+
+template <typename T, size_t N = 100>
+class blocking_queue
+{
+  public:
+	void push(T const &value)
+	{
+		std::unique_lock<std::mutex> lock(m_guard);
+
+		while (m_queue.size() >= N)
+			m_full_signal.wait(lock);
+
+		m_queue.push(value);
+
+		m_empty_signal.notify_one();
+	}
+
+	T pop()
+	{
+		std::unique_lock<std::mutex> lock(m_guard);
+		while (m_queue.empty())
+			m_empty_signal.wait(lock);
+
+		auto value = m_queue.front();
+		m_queue.pop();
+
+		m_full_signal.notify_one();
+
+		return value;
+	}
+
+	template<class Rep, class Period>
+	std::tuple<bool,T> pop(const std::chrono::duration<Rep, Period>& wait_for)
+	{
+		std::unique_lock<std::mutex> lock(m_guard);
+		while (m_queue.empty())
+		{
+		    auto now = std::chrono::system_clock::now();
+			if (m_empty_signal.wait_until(lock, now + wait_for) == std::cv_status::timeout)
+				return { true , T{} };
+		}
+
+		auto value = m_queue.front();
+		m_queue.pop();
+
+		m_full_signal.notify_one();
+
+		return { false, value };
+	}
+
+	bool is_full() const
+	{
+		std::unique_lock<std::mutex> lock(m_guard);
+		return m_queue.size() >= N;
+	}
+
+  private:
+	std::queue<T> m_queue;
+	mutable std::mutex m_guard;
+	std::condition_variable m_empty_signal, m_full_signal;
+};
+
+template <typename T, size_t N = 10>
+class non_blocking_queue
+{
+  public:
+	bool push(T const &value)
+	{
+		bool result = false;
+
+		std::unique_lock<std::mutex> lock(m_guard);
+
+		if (m_queue.size() < N)
+		{
+			m_queue.push(value);
+			m_empty_signal.notify_one();
+			result = true;
+		}
+			
+		return result;
+	}
+
+	T pop()
+	{
+		std::unique_lock<std::mutex> lock(m_guard);
+		while (m_queue.empty())
+			m_empty_signal.wait(lock);
+
+		auto value = m_queue.front();
+		m_queue.pop();
+
+		return value;
+	}
+
+  private:
+	std::queue<T> m_queue;
+	mutable std::mutex m_guard;
+	std::condition_variable m_empty_signal;
+};
\ No newline at end of file