From f8c32f68e894fda104cd16c8a140a7758b5984e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Florian=20B=C3=A4uerlein?= Date: Wed, 20 Nov 2013 15:10:18 +0100 Subject: [PATCH] initial commit --- AuctionAlgorithm/Auction.h | 644 ++++++++++++++++++++++++++++ AuctionAlgorithm/AuctionCommon.h | 528 +++++++++++++++++++++++ AuctionAlgorithm/AuctionMT.h | 478 +++++++++++++++++++++ LSAP.cpp | 46 ++ Matrix/SparseMatrix.h | 712 +++++++++++++++++++++++++++++++ MatrixHelper.h | 88 ++++ ThreadBarrier.h | 80 ++++ internal.h | 74 ++++ 8 files changed, 2650 insertions(+) create mode 100644 AuctionAlgorithm/Auction.h create mode 100644 AuctionAlgorithm/AuctionCommon.h create mode 100644 AuctionAlgorithm/AuctionMT.h create mode 100644 LSAP.cpp create mode 100644 Matrix/SparseMatrix.h create mode 100644 MatrixHelper.h create mode 100644 ThreadBarrier.h create mode 100644 internal.h diff --git a/AuctionAlgorithm/Auction.h b/AuctionAlgorithm/Auction.h new file mode 100644 index 0000000..3eb202e --- /dev/null +++ b/AuctionAlgorithm/Auction.h @@ -0,0 +1,644 @@ +/* + * Auction.h + * + * Created on: May 21, 2013 + * Author: fb + */ + +#ifndef AUCTION_H_ +#define AUCTION_H_ + +#include "../internal.h" +#include "AuctionCommon.h" + +//#define DEBUG_AUCTION + +namespace LSAP +{ + +/** + * + */ +template > +class Auction +{ +private: + + Auction() + { + } + + virtual ~Auction() + { + } +public: + typedef typename AuctionCommon::Scalars Scalars; + typedef typename AuctionCommon::Locks Locks; + typedef typename AuctionCommon::Indizes Indizes; + typedef typename AuctionCommon::Edges Edges; + + /** + * solve the assignment problem with the auction algorithm + * use real-types as coefficients, otherwise scaling will not work properly! + * @param a nxm weight matrix of type Scalar + * @return vector of Edges which represent the assignments + */ + static const Edges solve(MatrixType & a) + { +// const Scalar e = __AUCTION_EPSILON_MULTIPLIER; + const size_t rows = a.rows(); + const size_t cols = a.cols(); + + Locks lockedRows(a.rows(), false); + Locks lockedCols(a.cols(), false); + Edges E; + + Scalar lambda = .0; + Scalar epsilon = __AUCTION_EPSILON_MULTIPLIER / a.cols(); + + // condition 3: initially set p_j >= lambda + Scalars prices(cols, .0), profits(rows, 1.); // p-Vector (1 to j) = p_j + + do + { + // Step 1 (forward auction cycle): + // Execute iterations of the forward auction algorithm until at least one + // more person becomes assigned. If there is an unassigned person left, go + // to step 2; else go to step 3. + while (forward(a, E, prices, profits, lockedRows, lockedCols, + lambda, epsilon)) + ; + + if (!AuctionCommon::allPersonsAssigned(lockedRows)) + { + + // Step 2 (reverse auction cycle): + // Execute several iterations of the reverse auction algorithm until at least + // one more object becomes assigned or until we have p_j <= lambda for all + // unassigned objects. If there is an unassigned person left, go to step 1 + // else go to step 3 + while (!reverse(a, E, prices, profits, lockedRows, lockedCols, + lambda, epsilon) + || !AuctionCommon::unassignedObjectsLTlambda(lockedCols, prices, + lambda)) + ; // reverse auction + } + + if (AuctionCommon::allPersonsAssigned(lockedRows)) + { + // Step 3 (reverse auction): + // Execute successive iterations of the reverse auction algorithm until the + // algorithm terminates with p_j <= lambda for all unassigned objects j + while (true) + { + reverse(a, E, prices, profits, lockedRows, lockedCols, + lambda, epsilon); + if (AuctionCommon::unassignedObjectsLTlambda(lockedCols, prices, lambda)) + break; + } + break; + } + + } while (true); + +#ifdef DEBUG_AUCTION + for (auto i : E ) std::cout << i.x << ", " << i.y << " => " << prices[i.y] + profits[i.x] + << " a = " << a(i.x, i.y) << std::endl; + PRINTVECTOR(prices) + PRINTVECTOR(profits) +#endif + return E; + } +private: + + /** + * template specific implementation for finding the best entry in row + * @param a input matrix + * @param prices prices + * @param i row + * @param v_i best value + * @param w_i second best value + * @param a_i_ji value of entry + * @param j_i best entry + * @return true if assignment was found, otherwise false + */ + inline static const bool findForward( Eigen::Matrix & a, const Scalars & prices, const size_t i, + Scalar & v_i, Scalar & w_i, Scalar & a_i_ji, size_t & j_i) + { + const size_t cols = a.cols(); + + bool assignmentFound = false; + + for (size_t j = 0; j < cols; j++) // for the j-th column + { + const Scalar aij = a(i,j); + if ( aij == 0 ) continue; + const Scalar diff = aij - prices[j]; + if (diff > v_i) + { + // if there already was an entry found, this is the second best + if (assignmentFound) + w_i = v_i; + + v_i = diff; + j_i = j; + a_i_ji = aij; + assignmentFound = true; + } + if (diff > w_i && j_i != j) + w_i = diff; + // if no entry is bigger than v_i, check if there's still a bigger second best entry + } + return assignmentFound; + } + + /** + * template specific implementation for finding the best entry in row + * @param a input matrix + * @param prices prices + * @param i row + * @param v_i best value + * @param w_i second best value + * @param a_i_ji value of entry + * @param j_i best entry + * @return true if assignment was found, otherwise false + */ + inline static const bool findForward( SparseMatrix & a, const Scalars & prices, const size_t i, + Scalar & v_i, Scalar & w_i, Scalar & a_i_ji, size_t & j_i) + { + bool assignmentFound = false; + + const size_t outerStart = a.outerIndizesRM()[i]; + const size_t outerEnd = a.outerIndizesRM()[i + 1]; + + for (size_t innerIdx = outerStart; innerIdx < outerEnd; innerIdx++) // for the j-th column + { + size_t j = a.innerIndizesRM()[innerIdx]; + const Scalar m_i_j = a.valuesRM()[innerIdx]; + const Scalar diff = m_i_j - prices[j]; + if (diff > v_i) + { + // if there already was an entry found, this is the second best + if (assignmentFound) + w_i = v_i; + + v_i = diff; + j_i = j; + a_i_ji = m_i_j; + assignmentFound = true; + } + if (diff > w_i && j_i != j) + w_i = diff; + // if no entry is bigger than v_i, check if there's still a bigger second best entry + } + return assignmentFound; + } + + /** + * template specific implementation for finding the best entry in row + * @param a input matrix + * @param prices prices + * @param i row + * @param v_i best value + * @param w_i second best value + * @param a_i_ji value of entry + * @param j_i best entry + * @return true if assignment was found, otherwise false + */ + inline static const bool findForward( SparseMatrixRM & a, const Scalars & prices, const size_t i, + Scalar & v_i, Scalar & w_i, Scalar & a_i_ji, size_t & j_i) + { + bool assignmentFound = false; + + const size_t outerStart = a.outerIndizes()[i]; + const size_t outerEnd = a.outerIndizes()[i + 1]; + + for (size_t innerIdx = outerStart; innerIdx < outerEnd; innerIdx++) // for the j-th column + { + size_t j = a.innerIndizes()[innerIdx]; + const Scalar m_i_j = a.values()[innerIdx]; + const Scalar diff = m_i_j - prices[j]; + if (diff > v_i) + { + // if there already was an entry found, this is the second best + if (assignmentFound) + w_i = v_i; + + v_i = diff; + j_i = j; + a_i_ji = m_i_j; + assignmentFound = true; + } + if (diff > w_i && j_i != j) + w_i = diff; + // if no entry is bigger than v_i, check if there's still a bigger second best entry + } + return assignmentFound; + } + + /** + * template specific implementation for finding the best entry in row + * @param a input matrix + * @param prices prices + * @param i row + * @param v_i best value + * @param w_i second best value + * @param a_i_ji value of entry + * @param j_i best entry + * @return true if assignment was found, otherwise false + */ + inline static const bool findForward( SparseMatrixCM & a, const Scalars & prices, const size_t i, + Scalar & v_i, Scalar & w_i, Scalar & a_i_ji, size_t & j_i) + { + const size_t cols = a.cols(); + + bool assignmentFound = false; + + for (size_t j = 0; j < cols; j++) // for the j-th column + { + const Scalar aij = a(i,j); + if ( aij == 0 ) continue; + const Scalar diff = aij - prices[j]; + if (diff > v_i) + { + // if there already was an entry found, this is the second best + if (assignmentFound) + w_i = v_i; + + v_i = diff; + j_i = j; + a_i_ji = aij; + assignmentFound = true; + } + if (diff > w_i && j_i != j) + w_i = diff; + // if no entry is bigger than v_i, check if there's still a bigger second best entry + } + return assignmentFound; + } + + /** + * forward cycle of auction algorithm + * @param a weight matrix (nxm) + * @param S assignment matrix (nxm) + * @param prices prices per object (m) + * @param profits profits per person (n) + * @param lambda bidding threshold lambda + * @param epsilon bidding increment + * @return true if assignment was made, false otherwise + */ + static bool forward(MatrixType & a, Edges & E, + Scalars & prices, Scalars & profits, Locks & lockedRows, + Locks & lockedCols, Scalar & lambda, Scalar & epsilon) + { +#ifdef DEBUG_AUCTION + std::cout << "forwarding ..." << std::endl; +#endif + const size_t rows = a.rows(); + bool assignmentFound = false; + + for (size_t i = 0; i < rows; i++) // for the i-th row/person + { + bool assignmentInThisIterationFound = false; + + // person already assigned? + if (lockedRows[i]) + continue; + + // find an unassigned person i, its best object j_i + // j_i = argmax {a_ij - p_j} for j in A(i) ( A(i) are the set of edges of the i-th row ) + // if a(i,j) = 0. it is not a valid edge + size_t j_i = 0; + + // v_i = max { a_ij - p_j} for j in A(i) // maximum profit for person i + // v_i was already found = v_i + // w_i = max { a_ij - p_j} for j in A(i) and j != j_i // second best profit + // if j_i is the only entry in A(i), w_i = - inf // there's no second best profit + Scalar w_i = -__AUCTION_INF, v_i = -__AUCTION_INF, a_i_ji = 0.; // = max { a_ij - p_j} + + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + // no possible assignment found? + if (!findForward(a, prices, i, v_i, w_i, a_i_ji, j_i)) + continue; +#ifdef DEBUG_AUCTION + std::cout << "i = " << i << " - j_i = " << j_i << " - v_i = " << v_i + << " - w_i = " << w_i << " - a_ij = " << a_i_ji << std::endl; +#endif + assignmentInThisIterationFound = false; + +// std::cout << "assignment found .." << std::endl; + const Scalar bid = a_i_ji - w_i + epsilon; + + // P_i = w_i - E + profits[i] = w_i - epsilon; // set new profit for person + + // prices(j_i) = max(lambda, a(i,j_i) - w(i) + epsilon) + // if lambda <= a_ij - w_i + E, add (i, j_i) to S + if (lambda <= bid) + { + prices[j_i] = bid; + // assignment was made, so lock row and col + lockedRows[i] = true; + lockedCols[j_i] = true; + + bool newEdge = true; + + // if j_i was assigned to different i' to begin, remove (i', j_i) from S + for (auto & e : E) + if (e.y == j_i) // change edge + { + lockedRows[e.x] = false; // unlock row i' + newEdge = false; + e.x = i; + e.v = a_i_ji; +// std::cout << " edges: " << E.size() +// << " changing edge "; + break; + } + if (newEdge) + { + Edge e; + e.x = i; + e.y = j_i; + e.v = a_i_ji; + E.push_back(e); +// std::cout << "edges: " << E.size() << " added edge "; + } + assignmentInThisIterationFound = true; +// std::cout << "assignmentInThisIterationFound = true" +// << std::endl; + } + else + { + prices[j_i] = lambda; + assignmentInThisIterationFound = false; +// std::cout << "assignmentInThisIterationFound = false" +// << std::endl; + } + if (assignmentInThisIterationFound) + assignmentFound = true; + } + +// std::cout << "returning " << ((assignmentFound) ? "true" : "false") +// << std::endl; + return assignmentFound; + + } + + inline static const bool findReverse(const Eigen::Matrix & a, const Scalars & profits, const size_t j, + Scalar & b_j, Scalar & g_j, Scalar & a_ij_j, size_t & i_j) + { + const size_t rows = a.rows(); + bool assignmentFound = false; + + for (size_t i = 0; i < rows; i++) // for the j-th column + { + const Scalar aij = a(i, j); + if ( aij == 0 ) continue; + const Scalar diff = aij - profits[i]; + if (diff > b_j) + { + // if there already was an entry found, this is the second best + if (assignmentFound) + g_j = b_j; + + b_j = diff; + i_j = i; + a_ij_j = aij; + assignmentFound = true; + } + if (diff > g_j && i_j != i) + g_j = diff; + } + return assignmentFound; + } + + inline static const bool findReverse(SparseMatrix & a, const Scalars & profits, const size_t j, + Scalar & b_j, Scalar & g_j, Scalar & a_ij_j, size_t & i_j) + { + bool assignmentFound = false; + + const size_t outerStart = a.outerIndizesCM()[j]; + const size_t outerEnd = a.outerIndizesCM()[j + 1]; + + for (size_t innerIdx = outerStart; innerIdx < outerEnd; ++innerIdx) + { + const size_t i = a.innerIndizesCM()[innerIdx]; + const Scalar m_i_j = a.valuesCM()[innerIdx]; + const Scalar diff = m_i_j - profits[i]; + if (diff > b_j) + { + // if there already was an entry found, this is the second best + if (assignmentFound) + g_j = b_j; + + b_j = diff; + i_j = i; + a_ij_j = m_i_j; + assignmentFound = true; + } + if (diff > g_j && i_j != i) + g_j = diff; + } + return assignmentFound; + } + + inline static const bool findReverse(SparseMatrixRM & a, const Scalars & profits, const size_t j, + Scalar & b_j, Scalar & g_j, Scalar & a_ij_j, size_t & i_j) + { + const size_t rows = a.rows(); + bool assignmentFound = false; + + for (size_t i = 0; i < rows; i++) // for the j-th column + { + const Scalar aij = a(i, j); + if ( aij == 0 ) continue; + const Scalar diff = aij - profits[i]; + if (diff > b_j) + { + // if there already was an entry found, this is the second best + if (assignmentFound) + g_j = b_j; + + b_j = diff; + i_j = i; + a_ij_j = aij; + assignmentFound = true; + } + if (diff > g_j && i_j != i) + g_j = diff; + } + return assignmentFound; + } + + inline static const bool findReverse(SparseMatrixCM & a, const Scalars & profits, const size_t j, + Scalar & b_j, Scalar & g_j, Scalar & a_ij_j, size_t & i_j) + { + bool assignmentFound = false; + + const size_t outerStart = a.outerIndizes()[j]; + const size_t outerEnd = a.outerIndizes()[j + 1]; + + for (size_t innerIdx = outerStart; innerIdx < outerEnd; ++innerIdx) + { + const size_t i = a.innerIndizes()[innerIdx]; + const Scalar m_i_j = a.values()[innerIdx]; + const Scalar diff = m_i_j - profits[i]; + if (diff > b_j) + { + // if there already was an entry found, this is the second best + if (assignmentFound) + g_j = b_j; + + b_j = diff; + i_j = i; + a_ij_j = m_i_j; + assignmentFound = true; + } + if (diff > g_j && i_j != i) + g_j = diff; + } + return assignmentFound; + } + + /** + * reverse cycle of auction algorithm + * @param a weight matrix (nxm) + * @param S assignment matrix (nxm) + * @param prices prices per object (m) + * @param profits profits per person (n) + * @param lambda bidding threshold lambda + * @param epsilon bidding increment + * @return true if assignment was made, false otherwise + */ + static bool reverse(MatrixType & a, Edges & E, + Scalars & prices, Scalars & profits, Locks & lockedRows, + Locks & lockedCols, Scalar & lambda, const Scalar & epsilon) + { + const size_t rows = a.rows(); + const size_t cols = a.cols(); + + bool assignmentFound = false; + + for (size_t j = 0; j < cols; j++) // for the j-th column (objects) + { + bool assignmentInThisIterationFound = false; + + // object already assigned, p_j > lambda ? + if (lockedCols[j]) + continue; + + if (!(prices[j] > lambda)) + continue; + + // Find an unassigned object j with p_j > lambda, its best person i_j + // i_j = argmax {a_ij - profits[i]) für i aus B(j) (PI !!!) + size_t i_j = 0; + + //g_j = max {a_ij - P_i} for i in B(j) and i != i_j + // if j_i is the only entry in B(j), g_j = - inf ( g_j < b_j) + //b_j = max {a_ij - P_i} for i in B(j) + Scalar b_j = -__AUCTION_INF, g_j = -__AUCTION_INF, a_ij_j = 0.; + + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + // no assignment found + if (!findReverse(a, profits, j, b_j, g_j, a_ij_j, i_j)) + continue; + +#ifdef DEBUG_AUCTION + std::cout << "j = " << j << " i_j = " << i_j << " b_j = " << b_j << " g_j = " << g_j + << " a_ij_j = " << a_ij_j + << " p_j = " << prices[j] << " P_i = " << profits[i_j]<< std::endl; +#endif + assignmentInThisIterationFound = false; + + //if b_j >= L + E, case 1: + if (b_j >= (lambda + epsilon)) + { +#ifdef DEBUG_AUCTION + std::cout << " b_j >= lambda + epsilon" << std::endl; +#endif + const Scalar diff = g_j - epsilon; // G_j - E + + const Scalar max = lambda > diff ? lambda : diff; // max { L, G_j - E} + + // p_j = max { L, G_j - E} + prices[j] = max; + + // P_i_j = a_i_jj - max {L, G_j - E} + profits[i_j] = a_ij_j - max; + + lockedRows[i_j] = true; + lockedCols[j] = true; + + bool newEdge = true; + + // if j_i was assigned to different i' to begin, remove (i', j_i) from S + for (auto & e : E) + if (e.x == i_j) // change edge + { + lockedCols[e.y] = false; // unlock row i' + newEdge = false; + e.y = j; + e.v = a_ij_j; +#ifdef DEBUG_AUCTION + std::cout << " edges: " << E.size() + << " changing edge "; +#endif + break; + + } + if (newEdge) + { + Edge e; + e.x = i_j; + e.y = j; + e.v = a_ij_j; + E.push_back(e); +#ifdef DEBUG_AUCTION + std::cout << " added edge " << E.size() << " "; +#endif + } + assignmentInThisIterationFound = true; + } + else // if B_j < L + E, case 2 + { + // p_j = B_j - E + prices[j] = b_j - epsilon; +#ifdef DEBUG_AUCTION + std::cout << " b_j < lambda + epsilon " << std::endl; +#endif + /** standard lambda scaling **/ + size_t lowerThanLambda = 0; + Scalar newLambda = lambda; + + // if the number of objectes k with p_k < lambda is bigger than (rows - cols) + for (size_t k = 0; k < cols; k++) + { + if (prices[k] < lambda) // p_k < lambda + { + lowerThanLambda++; + if (prices[k] < newLambda) + newLambda = prices[k]; + } + } + // set new lambda +#ifdef DEBUG_AUCTION + std::cerr << " changing lambda from " << lambda << " to " << newLambda << std::endl; +#endif + if (lowerThanLambda >= (cols - rows)) + lambda = newLambda; + assignmentInThisIterationFound = false; + } + if (assignmentInThisIterationFound) + assignmentFound = true; + } + return assignmentFound; + } + +private: + +}; + +} /* namespace LSAP */ +#endif /* AUCTION_H_ */ diff --git a/AuctionAlgorithm/AuctionCommon.h b/AuctionAlgorithm/AuctionCommon.h new file mode 100644 index 0000000..1fbd017 --- /dev/null +++ b/AuctionAlgorithm/AuctionCommon.h @@ -0,0 +1,528 @@ +/* + * AuctionCommon.h + * + * Created on: Jul 2, 2013 + * Author: fb + */ + +#ifndef AUCTIONCOMMON_H_ +#define AUCTIONCOMMON_H_ + +#include "../internal.h" + +//#define __AUCTION_DEBUG +#define PRINTVECTOR(VEC) \ +{\ + std::cout << #VEC << " = ";\ + for (auto V : VEC)\ + std::cout << V << " ";\ + std::cout << std::endl;} + +namespace LSAP +{ + +template +class AuctionCommon +{ +public: + typedef typename std::vector Scalars; + typedef typename std::vector Locks; + typedef typename std::vector Indizes; + typedef typename std::vector Edges; + typedef typename Eigen::Matrix WeightMatrix; + + class BidResult + { + + public: + BidResult() : + index(0), bestIndex(0), bestBid(0), secondBestBid(0), bestMatrixValue( + 0), assignmentFound(false) + { + } + + BidResult(const size_t idx, const size_t bestIdx, const Scalar bestBid, + const Scalar secondBestBid, const Scalar bestMatrixValue) : + index(idx), bestIndex(bestIdx), bestBid(bestBid), secondBestBid( + secondBestBid), bestMatrixValue(bestMatrixValue), assignmentFound( + true) + { + } + + virtual ~BidResult() + { + } + + // set by thread + size_t index; // index of examined row or column + size_t bestIndex; // index of best entry + Scalar bestBid, secondBestBid; // corresponding bids (w_i, v_i) + Scalar bestMatrixValue; // a_ij + + bool assignmentFound; // assignment found ? + }; + + typedef typename std::vector BidResults; + + /** + * gauss-seidel forward cycle for dense matrix + * @param m + * @param threadId + * @param i + * @param fromColIdx + * @param toColIdx + * @param prices + * @param results + */ + static void forwardGS(const WeightMatrix & m, const size_t threadId, + const size_t & i, const size_t & fromColIdx, const size_t & toColIdx, + const Scalars & prices, BidResults & results) + { + results[threadId].assignmentFound = false; + size_t j_i = 0; + Scalar v_i = -__AUCTION_INF, a_i_ji = 0., w_i = -__AUCTION_INF; + + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + // iterate over columns + for (size_t j = fromColIdx; j < toColIdx; j++) // for the j-th column + { + const Scalar m_i_j = m(i, j); + if (m_i_j == 0) + continue; + + const Scalar diff = m_i_j - prices[j]; + + if (diff > v_i) + { + // if there already was an entry found, this is the second best + if (results[threadId].assignmentFound) + w_i = v_i; + + a_i_ji = m_i_j; + v_i = diff; + j_i = j; + results[threadId].assignmentFound = true; + } + if (diff > w_i && j_i != j) + w_i = diff; + // if no entry is bigger than v_i, check if there's still a bigger second best entry + } + + // best object for person i + results[threadId].bestIndex = j_i; + // v_i = max { a_ij - p_j} for j in A(i) + results[threadId].bestBid = v_i; + // w_i = max { a_ij - p_j} for j in A(i) and j != j_i + results[threadId].secondBestBid = w_i; + // a_i_ji = a(i, j_i), already found above + results[threadId].bestMatrixValue = a_i_ji; + return; + } + + /** + * gauss-seidel forward cycle for sparse matrix + * @param m + * @param threadId + * @param i + * @param fromColIdx + * @param toColIdx + * @param prices + * @param results + */ + static void forwardGS(const SparseMatrix & m, const size_t threadId, + const size_t & i, const size_t & fromColIdx, const size_t & toColIdx, + const Scalars & prices, BidResults & results) + { + results[threadId].assignmentFound = false; + +// std::cerr << "forward thread " << threadId << " - fromColIdx = " << fromColIdx << " - toColIdx = " << toColIdx << std::endl; + + size_t j_i = 0; + Scalar v_i = -__AUCTION_INF, a_i_ji = 0., w_i = -__AUCTION_INF; + +// std::cout << "threadForward: thread " << threadId << " iterating from " << fromColIdx << " to " << toColIdx << std::endl; + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + // iterate over columns + for (size_t innerIndex = fromColIdx; innerIndex < toColIdx; + innerIndex++) // for the j-th column + { + + const size_t j = m.innerIndexRM(innerIndex); + const Scalar m_i_j = m.valueRM(innerIndex); + const Scalar diff = m_i_j - prices[j]; + + if (diff > v_i) + { + // if there already was an entry found, this is the second best + if (results[threadId].assignmentFound) + w_i = v_i; + + a_i_ji = m_i_j; + v_i = diff; + j_i = j; + results[threadId].assignmentFound = true; + } + if (diff > w_i && j_i != j) + w_i = diff; + // if no entry is bigger than v_i, check if there's still a bigger second best entry + } + + // best object for person i + results[threadId].bestIndex = j_i; + // v_i = max { a_ij - p_j} for j in A(i) + results[threadId].bestBid = v_i; + // w_i = max { a_ij - p_j} for j in A(i) and j != j_i + results[threadId].secondBestBid = w_i; + // a_i_ji = a(i, j_i), already found above + results[threadId].bestMatrixValue = a_i_ji; + } + + /** + * gauss-seidel forward cycle for sparse rowmajor matrix + * @param m + * @param threadId + * @param i + * @param fromColIdx + * @param toColIdx + * @param prices + * @param results + */ + static void forwardGS(const SparseMatrixRM & m, const size_t threadId, + const size_t & i, const size_t & fromColIdx, const size_t & toColIdx, + const Scalars & prices, BidResults & results) + { + results[threadId].assignmentFound = false; + +// std::cerr << "forward thread " << threadId << " - fromColIdx = " << fromColIdx << " - toColIdx = " << toColIdx << std::endl; + + size_t j_i = 0; + Scalar v_i = -__AUCTION_INF, a_i_ji = 0., w_i = -__AUCTION_INF; + +// std::cout << "threadForward: thread " << threadId << " iterating from " << fromColIdx << " to " << toColIdx << std::endl; + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + // iterate over columns + for (size_t innerIndex = fromColIdx; innerIndex < toColIdx; + innerIndex++) // for the j-th column + { + + const size_t j = m.innerIndex(innerIndex); + const Scalar m_i_j = m.value(innerIndex); + const Scalar diff = m_i_j - prices[j]; + + if (diff > v_i) + { + // if there already was an entry found, this is the second best + if (results[threadId].assignmentFound) + w_i = v_i; + + a_i_ji = m_i_j; + v_i = diff; + j_i = j; + results[threadId].assignmentFound = true; + } + if (diff > w_i && j_i != j) + w_i = diff; + // if no entry is bigger than v_i, check if there's still a bigger second best entry + } + + // best object for person i + results[threadId].bestIndex = j_i; + // v_i = max { a_ij - p_j} for j in A(i) + results[threadId].bestBid = v_i; + // w_i = max { a_ij - p_j} for j in A(i) and j != j_i + results[threadId].secondBestBid = w_i; + // a_i_ji = a(i, j_i), already found above + results[threadId].bestMatrixValue = a_i_ji; + } + + /** + * gaus-seidel reverse cycle for sparse row major matrix + * @param threadId + * @param j + * @param fromRowIdx + * @param toRowIdx + * @param m + * @param profits + * @param results + */ + static void reverseGS(const SparseMatrixRM & m, const size_t threadId, + const size_t & j, const size_t & fromRowIdx, const size_t & toRowIdx, + const Scalars & profits, BidResults & results) + { + + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + // iterate over columns + + size_t i_j = 0; + Scalar b_j = -__AUCTION_INF, a_ij_j = 0., g_j = -__AUCTION_INF; // = max { a_ij - P_i} + + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + for (size_t i = fromRowIdx; i < toRowIdx; i++) // for the i-th row + { + const Scalar m_i_j = m(i, j); + if (m_i_j == 0) + continue; + + const Scalar diff = m_i_j - profits[i]; + + if (diff > b_j) + { + // if there already was an entry found, this is the second best + if (results[threadId].assignmentFound) + g_j = b_j; + + b_j = diff; + i_j = i; + a_ij_j = m_i_j; + + results[threadId].assignmentFound = true; + } + if (diff > g_j && i_j != i) + g_j = diff; + } + + // best person for object j + results[threadId].bestIndex = i_j; + // b_i = max { a_ij - P_i} for i in A(j) + results[threadId].bestBid = b_j; + // g_i = max { a_ij - P_i} for j in A(j) and i != i_j + results[threadId].secondBestBid = g_j; + // a_ij_j = a(i_j, j), already found above + results[threadId].bestMatrixValue = a_ij_j; + + return; + } + + /** + * gaus-seidel reverse cycle for dense matrix + * @param threadId + * @param j + * @param fromRowIdx + * @param toRowIdx + * @param m + * @param profits + * @param results + */ + static void reverseGS(const WeightMatrix & m, const size_t threadId, + const size_t & j, const size_t & fromRowIdx, const size_t & toRowIdx, + const Scalars & profits, BidResults & results) + { + + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + // iterate over columns + + size_t i_j = 0; + Scalar b_j = -__AUCTION_INF, a_ij_j = 0., g_j = -__AUCTION_INF; // = max { a_ij - P_i} + + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + for (size_t i = fromRowIdx; i < toRowIdx; i++) // for the i-th row + { + const Scalar m_i_j = m(i, j); + if (m_i_j == 0) + continue; + + const Scalar diff = m_i_j - profits[i]; + + if (diff > b_j) + { + // if there already was an entry found, this is the second best + if (results[threadId].assignmentFound) + g_j = b_j; + + b_j = diff; + i_j = i; + a_ij_j = m_i_j; + + results[threadId].assignmentFound = true; + } + if (diff > g_j && i_j != i) + g_j = diff; + } + + // best person for object j + results[threadId].bestIndex = i_j; + // b_i = max { a_ij - P_i} for i in A(j) + results[threadId].bestBid = b_j; + // g_i = max { a_ij - P_i} for j in A(j) and i != i_j + results[threadId].secondBestBid = g_j; + // a_ij_j = a(i_j, j), already found above + results[threadId].bestMatrixValue = a_ij_j; + + return; + } + + /** + * gauss-seidel reverse-cycle for sparse matrices + * @param m + * @param threadId + * @param j + * @param fromRowIdx + * @param toRowIdx + * @param profits + * @param results + */ + static void reverseGS(const SparseMatrix & m, const size_t threadId, + const size_t & j, const size_t & fromRowIdx, const size_t & toRowIdx, + const Scalars & profits, BidResults & results) + { + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + // iterate over columns +// std::cerr << "reverse thread " << threadId << " - fromRowIdx = " << fromRowIdx << " - toRowIdx = " << toRowIdx << std::endl; + + size_t i_j = 0; + Scalar b_j = -__AUCTION_INF, a_ij_j = 0., g_j = -__AUCTION_INF; // = max { a_ij - P_i} + + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + for (size_t innerIndex = fromRowIdx; innerIndex < toRowIdx; + innerIndex++) // for the j-th column + { + const size_t i = m.innerIndexCM(innerIndex); + + const Scalar m_i_j = m.valueCM(innerIndex); + const Scalar diff = m_i_j - profits[i]; + + if (diff > b_j) + { + // if there already was an entry found, this is the second best + if (results[threadId].assignmentFound) + g_j = b_j; + + b_j = diff; + i_j = i; + a_ij_j = m_i_j; + + results[threadId].assignmentFound = true; + } + if (diff > g_j && i_j != i) + g_j = diff; + } + + // best person for object j + results[threadId].bestIndex = i_j; + // b_i = max { a_ij - P_i} for i in A(j) + results[threadId].bestBid = b_j; + // g_i = max { a_ij - P_i} for j in A(j) and i != i_j + results[threadId].secondBestBid = g_j; + // a_ij_j = a(i_j, j), already found above + results[threadId].bestMatrixValue = a_ij_j; + } + + /** + * build prices and profits for auction from matrix + * @param a weightmatrix + * @param prices prices = p_j + * @param profits profits = p_i + * @param epsilon epsilon used for auction + */ + static void buildPricesAndProfits(SparseMatrix & a, + Scalars & prices, Scalars & profits, const Scalar epsilon, + const size_t from_row = 0, const size_t to_row = 0) + { + const size_t rowFirst = from_row; + const size_t rowLast = (to_row == 0) ? a.rows() : to_row; + + for (auto & p : profits) + p = 0.; + + // condition 2a/b: pi_i + p_j = a_ij - epsilon => pi_i = a_ij - epsilon - p_j + // if (r, c) is a valid edge ( weight > 0.) + // find maximum profit i.e. j_i = arg max { a_ij - p_j} and second best + + for (size_t r = rowFirst; r < rowLast; ++r) + { + for (size_t c = 0; c < a.cols(); ++c) + { + const Scalar newProfit = a(r, c) - prices[c] - epsilon; + if (newProfit > profits[r]) + profits[r] = newProfit; + } + } + } + + static Indizes splitToIntervals(const size_t parts, const size_t width, + const size_t offset = 0) + { + Indizes intervals(parts, offset); + const size_t base = width / parts; + size_t reminder = width - (base * parts); + size_t last = 0; + for (size_t i = 0; i < parts; ++i) + { + intervals[i] = last + base; + if (reminder > 0) + { + reminder--; + intervals[i]++; + } + last = intervals[i]; + } + return intervals; + } + + static Indizes splitToIntervalsByMatrixType(const SparseMatrix & m, + const size_t parts, const size_t index, bool major, size_t & offset) + { + const size_t entries = + (major) ? + m.outerIndexRM(index + 1) - m.outerIndexRM(index) : + m.outerIndexCM(index + 1) - m.outerIndexCM(index); + + offset = (major) ? m.outerIndexRM(index) : m.outerIndexCM(index); + + return AuctionCommon::splitToIntervals(parts, entries, 0); + } + + static Indizes splitToIntervalsByMatrixType(const SparseMatrixRM & m, + const size_t parts, const size_t index, bool major, size_t & offset) + { + const size_t entries = + m.outerIndex(index + 1) - m.outerIndex(index); + + offset = m.outerIndex(index); + + return AuctionCommon::splitToIntervals(parts, entries, 0); + } + + static Indizes splitToIntervalsByMatrixType(const WeightMatrix & m, + const size_t parts, const size_t index, bool major, size_t & offset) + { + const size_t entries = (major) ? m.rows() : m.cols(); + offset = 0; + return AuctionCommon::splitToIntervals(parts, entries); + } + /** + * returns true if p_j <= lambda for all unassigned objects. + * + * @param c locked columns + * @param prices prices of objects + * @param lambda bidding threshold + * @return true if all prices of unassigned objects are below lambda, otherwise false + */ + static const bool unassignedObjectsLTlambda(const Locks & c, + const Scalars & prices, const Scalar lambda) + { + for (size_t j = 0; j < c.size(); ++j) + if (!c[j] && prices[j] > lambda) + return false; + return true; + } + + /** + * check if all persons are assigned + * @return true if all persons are assigned, otherwise false + */ + static const bool allPersonsAssigned(const Locks & r) + { + for (size_t i = 0; i < r.size(); ++i) + if (!r[i]) + return false; + return true; + } +private: + AuctionCommon(); + virtual ~AuctionCommon(); +} +; + +} /* namespace LSAP */ +#endif /* AUCTIONCOMMON_H_ */ diff --git a/AuctionAlgorithm/AuctionMT.h b/AuctionAlgorithm/AuctionMT.h new file mode 100644 index 0000000..6654751 --- /dev/null +++ b/AuctionAlgorithm/AuctionMT.h @@ -0,0 +1,478 @@ +/* + * Auction.h + * + * Created on: May 21, 2013 + * Author: fb + */ + +#ifndef AUCTIONMT_H_ +#define AUCTIONMT_H_ + +#include "../internal.h" +#include "../ThreadBarrier.h" +#include "AuctionCommon.h" +#include + +namespace LSAP +{ + +//#define __AUCTION_DEBUG + +/** + * + */ +template > +class AuctionMT +{ +public: + + typedef typename AuctionCommon::Scalars Scalars; + typedef typename AuctionCommon::Locks Locks; + typedef typename AuctionCommon::Indizes Indizes; + typedef typename AuctionCommon::Edges Edges; + typedef typename AuctionCommon::BidResults BidResults; + + AuctionMT() + { + } + + virtual ~AuctionMT() + { + } + + /** + * solve the assignment with the auction algorithm + * use real-types as coefficients, otherwise scaling will not work properly! + * @param a nxm weight matrix of type Scalar + * @return assignment matrix which has maximum value of objective function + */ + static Edges solve(MatrixType & m, const size_t noOfThreads = 4) //, const Scalar e = __AUCTION_EPSILON_MULTIPLIER) + { + Edges edges; + Scalars prices(m.cols(), 0.); + Scalars profits(m.rows(), 1.); + + Scalar lambda = 0.; + Scalar epsilon = (Scalar) __AUCTION_EPSILON_MULTIPLIER / m.cols(); + Locks lockedRows(m.rows(), false); + Locks lockedCols(m.cols(), false); + + std::vector threadsForward, threadsReverse; + + // barriers for synchronization + ThreadBarrier barrierParamsForward(noOfThreads); + ThreadBarrier barrierResultsForward(noOfThreads); + ThreadBarrier barrierParamsReverse(noOfThreads); + ThreadBarrier barrierResultsReverse(noOfThreads); + + Indizes iterationIntervalsForward(noOfThreads); + + Indizes iterationIntervalsReverse (noOfThreads); + + // if matrix is dense eigen-type => pre-fill the iteration intervals for threads + // otherwise do this within the forward/reverse iteration, because intervals depend on + // column/row + if ( std::is_same>::value ) + { + iterationIntervalsForward = AuctionCommon::splitToIntervals(noOfThreads, m.cols()); + iterationIntervalsReverse = AuctionCommon::splitToIntervals(noOfThreads, m.rows()); + } + + BidResults results(noOfThreads); + + size_t index = 0; // representing row or column which is searched + + bool threadsActive = true; + + //bind threads to references and run them + for (size_t t = 1; t < noOfThreads; t++) + { + std::function < void() > forwardThread = std::bind(threadForward, std::cref(m), t, + std::cref(index), std::cref(iterationIntervalsForward[t-1]), + std::cref(iterationIntervalsForward[t]), + std::cref(prices), std::ref(results), + std::ref(barrierParamsForward), + std::ref(barrierResultsForward), std::ref(threadsActive)); + + threadsForward.emplace_back(std::thread(forwardThread)); + + std::function < void() > reverseThread = std::bind(threadReverse, std::cref(m), t, + std::cref(index), std::cref(iterationIntervalsReverse[t-1]), + std::cref(iterationIntervalsReverse[t]), + std::cref(profits), std::ref(results), + std::ref(barrierParamsReverse), + std::ref(barrierResultsReverse), std::ref(threadsActive)); + + threadsReverse.emplace_back(std::thread(reverseThread)); + } + + do + { + // Step 1 (forward auction cycle): + // Execute iterations of the forward auction algorithm until at least one + // more person becomes assigned. If there is an unassigned person left, go + // to step 2; else go to step 3. + while (parallelForward(m, prices, profits, lockedRows, lockedCols, + edges, lambda, epsilon, results, iterationIntervalsForward, + noOfThreads, index, barrierParamsForward, + barrierResultsForward)) + ; + + if (!AuctionCommon::allPersonsAssigned(lockedRows)) + { + // Step 2 (reverse auction cycle): + // Execute several iterations of the reverse auction algorithm until at least + // one more object becomes assigned or until we have p_j <= lambda for all + // unassigned objects. If there is an unassigned person left, go to step 1 + // else go to step 3 + while (!parallelReverse(m, prices, profits, lockedRows, + lockedCols, edges, lambda, epsilon, results, + iterationIntervalsReverse, noOfThreads, index, + barrierParamsReverse, barrierResultsReverse) + || !AuctionCommon::unassignedObjectsLTlambda( + lockedCols, prices, lambda)) + ; // reverse auction + } + if (AuctionCommon::allPersonsAssigned(lockedRows)) + { + // Step 3 (reverse auction): + // Execute successive iterations of the reverse auction algorithm until the + // algorithm terminates with p_j <= lambda for all unassigned objects j + while (true) + { + parallelReverse(m, prices, profits, lockedRows, lockedCols, + edges, lambda, epsilon, results, + iterationIntervalsReverse, noOfThreads, index, + barrierParamsReverse, barrierResultsReverse); + if (AuctionCommon::unassignedObjectsLTlambda( + lockedCols, prices, lambda)) + break; + } + break; + } + } while (true); + + // kill threads hanging in barrier + threadsActive = false; + barrierParamsForward.kill(); + barrierParamsReverse.kill(); + barrierResultsForward.kill(); + barrierResultsReverse.kill(); + + for (size_t t = 0; t < threadsForward.size(); t++) + threadsForward[t].join(); + + for (size_t t = 0; t < threadsReverse.size(); t++) + threadsReverse[t].join(); + + return edges; + } + +private: + + static void threadForward(const MatrixType & m, const size_t threadId, const size_t & i, + const size_t & fromColIdx, const size_t & toColIdx, const Scalars & prices, + BidResults & results, ThreadBarrier & barrierParams, + ThreadBarrier & barrierResult, bool & threadActive) + { + while (threadActive) + { + if (barrierParams() ) return; + AuctionCommon::forwardGS(m, threadId, i, fromColIdx, toColIdx, prices, results); + if ( barrierResult() ) return; + } + } + + static void threadReverse(const MatrixType & m, const size_t threadId, const size_t & j, + const size_t & fromRowIdx, const size_t & toRowIdx, const Scalars & profits, + BidResults & results, ThreadBarrier & barrierParams, + ThreadBarrier & barrierResult, bool & threadActive) + { + + while (threadActive) + { + if (barrierParams()) return; + AuctionCommon::reverseGS(m, threadId, j, fromRowIdx, toRowIdx, profits, results); + if (barrierResult()) return; + } + } + + static bool parallelReverse(const MatrixType & m, Scalars & prices, + Scalars & profits, Locks & lockedRows, Locks & lockedCols, + Edges & edges, Scalar & lambda, const Scalar epsilon, + BidResults & results, Indizes & iterationIntervals, + const size_t noOfThreads, size_t & index, + ThreadBarrier & barrierParams, ThreadBarrier & barrierResult) + { + + bool assignmentFound = false; + + // wake up reverse threads + barrierParams.wakeup(); + barrierResult.wakeup(); + + // for all cols + for (size_t j = 0; j < m.cols(); ++j) + { + // column already locked ? p_j > lambda ? + if (lockedCols[j]) + continue; + + if (!(prices[j] > lambda)) + continue; + + index = j; // global index for threads + + // if sparse-matrix is used: split row/column dependent intervals for threads + // offset is used for inner index offset + size_t offset = 0; + if (std::is_same>::value) + { + Indizes it = + AuctionCommon::splitToIntervalsByMatrixType(m, noOfThreads, j, false, offset); + for ( size_t i = 0; i < it.size(); ++i ) iterationIntervals[i] = offset + it[i]; +// std::cout << "reverse: offset = " << offset << " "; PRINTVECTOR(iterationIntervals) + } + + // sync and run local function for id 0 + barrierParams(); + AuctionCommon::reverseGS(m, 0, j, offset, iterationIntervals[0], profits, results); + barrierResult(); + + // values and indizes for the best assignment + size_t i_j = 0; + Scalar b_j = -__AUCTION_INF, g_j = -__AUCTION_INF; + Scalar a_ij_j = 0.; + bool foundBest = false; // found assignment? + + // find best and second best assignment in this row + for (size_t t = 0; t < noOfThreads; t++) + { + // was an possible assignment found by the thread? + if (results[t].assignmentFound) + { + bool found_b_j = false; + + // find best b_j and corresponding assignment + if (results[t].bestBid > b_j) + { + g_j = b_j; + a_ij_j = results[t].bestMatrixValue; + i_j = results[t].bestIndex; + b_j = results[t].bestBid; + foundBest = found_b_j = true; + } + if (g_j < results[t].bestBid && !found_b_j) + g_j = results[t].bestBid; + if (g_j < results[t].secondBestBid) + g_j = results[t].secondBestBid; + } + } + + if (foundBest) // best assignment found in row? + { + if (b_j >= lambda + epsilon) + { + const Scalar diff = g_j - epsilon; // G_j - E + + const Scalar max = lambda > diff ? lambda : diff; // max { L, G_j - E} + + // p_j = max { L, G_j - E} + prices[j] = max; + + // P_i_j = a_ij_j - max {L, G_j - E} + profits[i_j] = a_ij_j - max; + + lockedRows[i_j] = true; + lockedCols[j] = true; + + bool newEdgeFound = true; + + // if j_i was assigned to different i' to begin, remove (i', j_i) from S + for (auto & e : edges) + if (e.x == i_j) // change edge + { + lockedCols[e.y] = false; // unlock row i' + newEdgeFound = false; + e.y = j; + e.v = a_ij_j; + break; + } + if (newEdgeFound) + { + Edge newEdge; + newEdge.x = i_j; + newEdge.y = j; + newEdge.v = a_ij_j; + edges.push_back(newEdge); + } + + assignmentFound = true; + } + else // if B_j < L + E, case 2 + { + // p_j = B_j - E + prices[j] = b_j - epsilon; + + /** standard lambda scaling **/ + unsigned int lowerThanLambda = 0; + Scalar newLambda = lambda; + + // if the number of objectes k with p_k < lambda is bigger than (rows - cols) + for (size_t k = 0; k < m.cols(); k++) + { + if (prices[k] < lambda) + { + lowerThanLambda++; + if (prices[k] < newLambda) + newLambda = prices[k]; + } + } + if (lowerThanLambda >= (m.cols() - m.rows())) + lambda = newLambda; + } + } + } + + barrierParams.sleep(); + barrierResult.sleep(); + + return assignmentFound; + + } + + static bool parallelForward(const MatrixType & m, Scalars & prices, + Scalars & profits, Locks & lockedRows, Locks & lockedCols, + Edges & edges, const Scalar lambda, const Scalar epsilon, + BidResults & results, Indizes & iterationIntervals, + const size_t noOfThreads, size_t & index, + ThreadBarrier & barrierParams, ThreadBarrier & barrierResult) + { + bool assignmentFound = false; + // for all rows + + barrierParams.wakeup(); + barrierResult.wakeup(); + + for (size_t i = 0; i < m.rows(); ++i) + { + // if row is locked, there's already an assignment + if (lockedRows[i]) + continue; + + index = i; + + // if sparse-matrix is used: split row/column dependent intervals for threads + // offset is used for inner index offset + size_t offset = 0; + + if (std::is_same>::value || std::is_same>::value ) + { + Indizes it = + AuctionCommon::splitToIntervalsByMatrixType(m, noOfThreads, i, true, offset ); + for ( size_t i = 0; i < it.size(); ++i ) iterationIntervals[i] = offset + it[i]; + } +// std::cout << "offset = " << offset << " "; +// PRINTVECTOR(iterationIntervals) + // BARRIER + barrierParams(); + AuctionCommon::forwardGS(m, 0, i, offset, iterationIntervals[0], prices, results); + barrierResult(); + +// for (size_t t = 0; t < noOfThreads; t++) +// pthread_join(threads[t], NULL); + + // values and indizes for the best assignment + size_t j_i = 0; + Scalar v_i = -__AUCTION_INF, w_i = -__AUCTION_INF; + Scalar a_i_ji = 0.; + bool foundBest = false; // found assignment? + + // find best and second best assignment in this row + for (size_t t = 0; t < results.size(); t++) + { + // was an possible assignment found by the thread? + if (results[t].assignmentFound) + { + bool found_v_i = false; + // find best v_i and corresponding assignment + if (results[t].bestBid > v_i) + { + w_i = v_i; + a_i_ji = results[t].bestMatrixValue; + j_i = results[t].bestIndex; + v_i = results[t].bestBid; + foundBest = found_v_i = true; + } + + if (w_i < results[t].bestBid && !found_v_i) + w_i = results[t].bestBid; + if (w_i < results[t].secondBestBid) + w_i = results[t].secondBestBid; + } + } + + if (foundBest) // best assignment found in row? + { +// std::cout << "parallelForward: found best assignment (" << i << ", " +// << j_i << ") with v_i = " << v_i +// << " & w_i = " << w_i +// << " a_i_ji = " << a_i_ji +// << std::endl; + + // P_i = w_i - epsilon + profits[i] = w_i - epsilon; + + // diff = a(i,j_i) - w(i) + epsilon + // prices(j_i) = max(lambda, a(i,j_i) - w(i) + epsilon) + const Scalar diff = a_i_ji - w_i + epsilon; + + if (lambda <= diff) + { + prices[j_i] = diff; + assignmentFound = true; // new assignment was found + + // assignment was made, so lock row and col + lockedRows[i] = true; + lockedCols[j_i] = true; + + bool newEdgeFound = true; + + // if j_i was assigned to different i' to begin, remove (i', j_i) from S + for (auto & e : edges) + if (e.y == j_i) // change edge to new assignment + { +// std::cout << "changing edge (" << e.x << ", " << e.y << ") " << "-> (" << i << ", " << j_i << ")" << std::endl; + + newEdgeFound = false; + lockedRows[e.x] = false; // unlock row + e.v = a_i_ji; + e.x = i; + break; + } + + if (newEdgeFound) // add new edge to list + { + Edge newEdge; + newEdge.x = i; + newEdge.y = j_i; + newEdge.v = a_i_ji; + edges.push_back(newEdge); +// std::cout << "adding edge (" << i << ", " << j_i << ")" << std::endl; + } + } + else + prices[j_i] = lambda; + } + } + + barrierParams.sleep(); + barrierResult.sleep(); + + return assignmentFound; + } + +}; + +} /* namespace LSAP */ +#endif /* AUCTIONMT_H_ */ diff --git a/LSAP.cpp b/LSAP.cpp new file mode 100644 index 0000000..a95e73d --- /dev/null +++ b/LSAP.cpp @@ -0,0 +1,46 @@ + +#include "AuctionAlgorithm/Auction.h" +#include "AuctionAlgorithm/AuctionMT.h" +#include "Matrix/SparseMatrix.h" +#include +#include +#include +using namespace LSAP; + +typedef double Scalar; + +typedef Eigen::Matrix WeightMatrix; + +int main(int argc, char **argv) +{ + const size_t rows = 2000, cols = 5000; + + // assert that rows <= cols and coefficients are between 0 and 1! + Eigen::MatrixXd m = Eigen::MatrixXd::Random(rows, cols); + + // shift coefficients to get positive values + for ( size_t i = 0; i < rows; ++i) + for ( size_t j = 0; j < cols; ++j ) + m(i, j) += 1.; + + m /= m.maxCoeff(); // normalize to 0..1 + + // create sparse matrix from dense + // this could take some time! + // currently there's an own sparse class, might be changed to eigen sparse type + // this matrix type stores the matrix in rowMajor AND colMajor format, so it uses + // a lot of space! (this was due to testing purposes) + SparseMatrix s(m); + + // result type + Edges solution; + + // single threaded computation and some time measurement + MEASURE_DURATION_SINGLE((solution = Auction::solve(m))); + MEASURE_DURATION_SINGLE((solution = Auction >::solve(s))); + + // multi threaded computation (2 threads) and time measurement + MEASURE_DURATION_SINGLE((solution = AuctionMT::solve(m, 2))); + MEASURE_DURATION_SINGLE((solution = AuctionMT >::solve(s, 2))); + +} diff --git a/Matrix/SparseMatrix.h b/Matrix/SparseMatrix.h new file mode 100644 index 0000000..a31092b --- /dev/null +++ b/Matrix/SparseMatrix.h @@ -0,0 +1,712 @@ +/* + * CRSMatrix.h + * + * Created on: May 8, 2013 + * Author: fb + */ + +#ifndef SPARSEMATRIX_H_ +#define SPARSEMATRIX_H_ + +#include "../internal.h" + +#define __SPARSE_MATRIX_EPSILON 0 +#define LOG_SPARSEMATRIX __PROJECT_LOG << "[SparseMatrix] " +//#define DEBUG + +namespace LSAP +{ + +template +class SparseMatrix +{ +public: + + enum CompressionType + { + rowMajor = true, colMajor = false + }; + + typedef Eigen::Matrix weightMatrix; + + typedef typename std::vector Scalars; + typedef std::vector Indizes; + + SparseMatrix() : + _rmOuterSize(0), _rmInnerSize(0), _cmOuterSize(0), _cmInnerSize(0), _NumberOfValues( + 0), _rows(0), _cols(0), _zero(0.) + { + _CompressionType = rowMajor; + } + + SparseMatrix(const weightMatrix & m) : + _rmOuterSize(0), _rmInnerSize(0), _cmOuterSize(0), _cmInnerSize(0), _NumberOfValues( + 0), _rows(0), _cols(0), _zero(0.) + { + _CompressionType = rowMajor; // compatibility + + constructFromEigenMatrix(m); + +#ifdef DEBUG + LOG_SPARSEMATRIX << "input matrix: " << std::endl << m << std::endl; + + LOG_SPARSEMATRIX << "constructing from Eigen matrix ..." << std::endl; + + LOG_SPARSEMATRIX << "using " << (_CompressionType ? "row" : "col") << " major representation! " << std::endl; + + LOG_SPARSEMATRIX << "values = "; + for ( size_t i = 0; i < _NumberOfValues; ++i) __PROJECT_LOG << _rmValues[i] << " "; + __PROJECT_LOG << "\n"; + + LOG_SPARSEMATRIX << "InnerIndizes = "; + for ( size_t i = 0; i < _NumberOfValues; ++i) __PROJECT_LOG << _rmInnerIndizes[i] << " "; + __PROJECT_LOG << "\n"; + + LOG_SPARSEMATRIX << "OuterIndizes = "; + for ( size_t i = 0; i < _rmOuterSize; ++i) __PROJECT_LOG << _rmOuterIndizes[i] << " "; + __PROJECT_LOG << "\n"; + +#endif + } + + virtual ~SparseMatrix() + { + + } + + /** + * eigen like operator for accessing the matrix-values + * @param r row + * @param c column + * @return value of (r,c) + */ + const Scalar & operator()(const size_t r, const size_t c) const + { + assert(r >= 0 && r < _rows && c >= 0 && c < _cols); + + // determine outer and inner index values + const size_t outerIdx = (_CompressionType) ? r : c; + const size_t innerIdx = (_CompressionType) ? c : r; + + // get end-index of outer indizes + for (size_t i = _rmOuterIndizes[outerIdx]; + i < _rmOuterIndizes[outerIdx + 1]; ++i) + { + if (_rmInnerIndizes[i] > innerIdx) + return _zero; + if (_rmInnerIndizes[i] == innerIdx) + return _rmValues[i]; + } + + return _zero; + } + +private: + + void constructFromEigenMatrix(const weightMatrix & m) + { + _rows = m.rows(); + _cols = m.cols(); + + // count number of non-zero values + for (size_t r = 0; r < _rows; ++r) + for (size_t c = 0; c < _cols; ++c) + if (abs(m(r, c) > __SPARSE_MATRIX_EPSILON)) + _NumberOfValues++; + + // determine array sizes dependend on compression type + _rmOuterSize = _rows + 1; + _rmInnerSize = _cols; + + // determine array sizes dependend on compression type + _cmOuterSize = _cols + 1; + _cmInnerSize = _rows; + + // allocate values and indizes + _rmOuterIndizes = Indizes(_rmOuterSize); + _rmInnerIndizes = Indizes(_NumberOfValues); + _rmValues = Scalars(_NumberOfValues); + + _cmOuterIndizes = Indizes(_cmOuterSize); + _cmInnerIndizes = Indizes(_NumberOfValues); + _cmValues = Scalars(_NumberOfValues); + + // set allocated index-array to 0 + for (auto & i : _rmOuterIndizes) + i = 0; + for (auto & i : _cmOuterIndizes) + i = 0; + + // last position holds the size of values/innerindizes + _rmOuterIndizes[_rmOuterSize - 1] = _NumberOfValues; + _cmOuterIndizes[_cmOuterSize - 1] = _NumberOfValues; + + // read matrix data for row-major representation + size_t valuesSoFar = 0; + + for (size_t o = 0; o < _rmOuterSize - 1; ++o) + { + size_t valuesPerOuter = 0; + for (size_t i = 0; i < _rmInnerSize; ++i) + { + const Scalar mVal = m(o, i); + if (mVal != 0) + { + _rmValues[valuesSoFar] = mVal; + _rmInnerIndizes[valuesSoFar] = i; + valuesPerOuter++; + valuesSoFar++; + } + } + _rmOuterIndizes[o + 1] = _rmOuterIndizes[o] + valuesPerOuter; + } + + // read matrix data for col-major representation + valuesSoFar = 0; + + for (size_t o = 0; o < _cmOuterSize - 1; ++o) + { + size_t valuesPerOuter = 0; + for (size_t i = 0; i < _cmInnerSize; ++i) + { + const Scalar mVal = m(i, o); + if (mVal != 0) + { + _cmValues[valuesSoFar] = mVal; + _cmInnerIndizes[valuesSoFar] = i; + valuesPerOuter++; + valuesSoFar++; + } + } + _cmOuterIndizes[o + 1] = _cmOuterIndizes[o] + valuesPerOuter; + } + } + +public: + + inline const Scalar maxCoeff() const + { + Scalar c = 0.; + for (size_t i = 0; i < _NumberOfValues; ++i) + c = (_rmValues[i] > c) ? _rmValues[i] : c; + return c; + } + + inline const Scalar minNonZeroCoeff() const + { + assert(_rmValues.size() > 0); + Scalar c = _rmValues[0]; + for (size_t i = 0; i < _NumberOfValues; ++i) + c = (_rmValues[i] < c) ? _rmValues[i] : c; + return c; + } + + inline Scalars & valuesRM() + { + return _rmValues; + } + inline Indizes & outerIndizesRM() + { + return _rmOuterIndizes; + } + inline Indizes & innerIndizesRM() + { + return _rmInnerIndizes; + } + + inline const size_t & outerIndexRM(const size_t o) const + { + return _rmOuterIndizes[o]; + } + inline const size_t & innerIndexRM(const size_t i) const + { + return _rmInnerIndizes[i]; + } + inline const Scalar & valueRM(const size_t i) const + { + return _rmValues[i]; + } + + inline const size_t & outerIndexCM(const size_t o) const + { + return _cmOuterIndizes[o]; + } + inline const size_t & innerIndexCM(const size_t i) const + { + return _cmInnerIndizes[i]; + } + inline const Scalar & valueCM(const size_t i) const + { + return _cmValues[i]; + } + + inline Scalars & valuesCM() + { + return _cmValues; + } + inline Indizes & outerIndizesCM() + { + return _cmOuterIndizes; + } + inline Indizes & innerIndizesCM() + { + return _cmInnerIndizes; + } + + inline const size_t & rows() const + { + return _rows; + } + inline const size_t & cols() const + { + return _cols; + } + + inline const size_t size() const + { + return _NumberOfValues; + } + + inline const bool compressionType() const + { + return _CompressionType; + } + + // compatibility to old version + inline Scalars & values() + { + return _rmValues; + } + inline Indizes & outerIndizes() + { + return _rmOuterIndizes; + } + inline Indizes & innerIndizes() + { + return _rmInnerIndizes; + } + +private: + // row-major representation + Scalars _rmValues; + Indizes _rmOuterIndizes; // outer index (holds rowIdx if CRS, colIdx if CCS) + Indizes _rmInnerIndizes; // inner index + size_t _rmOuterSize; + size_t _rmInnerSize; + + // col-major representation + Scalars _cmValues; + Indizes _cmOuterIndizes; // outer index (holds rowIdx if CRS, colIdx if CCS) + Indizes _cmInnerIndizes; // inner index + size_t _cmOuterSize; + size_t _cmInnerSize; + + // common + size_t _NumberOfValues; + + size_t _rows, _cols; + + CompressionType _CompressionType; + + Scalar _zero; // for returning reference to values + +}; + +template +class SparseMatrixRM +{ +public: + + typedef Eigen::Matrix weightMatrix; + + typedef typename std::vector Scalars; + typedef std::vector Indizes; + + SparseMatrixRM() : + _rmOuterSize(0), _rmInnerSize(0), _NumberOfValues(0), _rows(0), _cols(0), _zero(0.) + {} + + SparseMatrixRM(const weightMatrix & m) : + _rmOuterSize(0), _rmInnerSize(0), _NumberOfValues( + 0), _rows(0), _cols(0), _zero(0.) + { + constructFromEigenMatrix(m); + +#ifdef DEBUG + LOG_SPARSEMATRIX << "input matrix: " << std::endl << m << std::endl; + + LOG_SPARSEMATRIX << "constructing from Eigen matrix ..." << std::endl; + + LOG_SPARSEMATRIX << "using " << (_CompressionType ? "row" : "col") << " major representation! " << std::endl; + + LOG_SPARSEMATRIX << "values = "; + for ( size_t i = 0; i < _NumberOfValues; ++i) __PROJECT_LOG << _rmValues[i] << " "; + __PROJECT_LOG << "\n"; + + LOG_SPARSEMATRIX << "InnerIndizes = "; + for ( size_t i = 0; i < _NumberOfValues; ++i) __PROJECT_LOG << _rmInnerIndizes[i] << " "; + __PROJECT_LOG << "\n"; + + LOG_SPARSEMATRIX << "OuterIndizes = "; + for ( size_t i = 0; i < _rmOuterSize; ++i) __PROJECT_LOG << _rmOuterIndizes[i] << " "; + __PROJECT_LOG << "\n"; + +#endif + } + + virtual ~SparseMatrixRM() + { + + } + + /** + * eigen like operator for accessing the matrix-values + * @param r row + * @param c column + * @return value of (r,c) + */ + const Scalar & operator()(const size_t r, const size_t c) const + { + assert(r >= 0 && r < _rows && c >= 0 && c < _cols); + + // get end-index of outer indizes + for (size_t i = _rmOuterIndizes[r]; + i < _rmOuterIndizes[r + 1]; ++i) + { + if (_rmInnerIndizes[i] < c) continue; + if (_rmInnerIndizes[i] == c) + return _rmValues[i]; + if (_rmInnerIndizes[i] > c) + return _zero; + } + + return _zero; + } + +private: + + void constructFromEigenMatrix(const weightMatrix & m) + { + _rows = m.rows(); + _cols = m.cols(); + + // count number of non-zero values + for (size_t r = 0; r < _rows; ++r) + for (size_t c = 0; c < _cols; ++c) + if (abs(m(r, c) > __SPARSE_MATRIX_EPSILON)) + _NumberOfValues++; + + // determine array sizes dependend on compression type + _rmOuterSize = _rows + 1; + _rmInnerSize = _cols; + + // allocate values and indizes + _rmOuterIndizes = Indizes(_rmOuterSize); + _rmInnerIndizes = Indizes(_NumberOfValues); + _rmValues = Scalars(_NumberOfValues); + + // set allocated index-array to 0 + for (auto & i : _rmOuterIndizes) + i = 0; + + // last position holds the size of values/innerindizes + _rmOuterIndizes[_rmOuterSize - 1] = _NumberOfValues; + + // read matrix data for row-major representation + size_t valuesSoFar = 0; + + for (size_t o = 0; o < _rmOuterSize - 1; ++o) + { + size_t valuesPerOuter = 0; + for (size_t i = 0; i < _rmInnerSize; ++i) + { + const Scalar mVal = m(o, i); + if (mVal != 0) + { + _rmValues[valuesSoFar] = mVal; + _rmInnerIndizes[valuesSoFar] = i; + valuesPerOuter++; + valuesSoFar++; + } + } + _rmOuterIndizes[o + 1] = _rmOuterIndizes[o] + valuesPerOuter; + } + + } + +public: + + inline const Scalar maxCoeff() const + { + Scalar c = 0.; + for (size_t i = 0; i < _NumberOfValues; ++i) + c = (_rmValues[i] > c) ? _rmValues[i] : c; + return c; + } + + inline const Scalar minNonZeroCoeff() const + { + assert(_rmValues.size() > 0); + Scalar c = _rmValues[0]; + for (size_t i = 0; i < _NumberOfValues; ++i) + c = (_rmValues[i] < c) ? _rmValues[i] : c; + return c; + } + + inline Scalars & values() + { + return _rmValues; + } + inline Indizes & outerIndizes() + { + return _rmOuterIndizes; + } + inline Indizes & innerIndizes() + { + return _rmInnerIndizes; + } + + inline const size_t & outerIndex(const size_t o) const + { + return _rmOuterIndizes[o]; + } + inline const size_t & innerIndex(const size_t i) const + { + return _rmInnerIndizes[i]; + } + inline const Scalar & value(const size_t i) const + { + return _rmValues[i]; + } + + inline const size_t & rows() const + { + return _rows; + } + inline const size_t & cols() const + { + return _cols; + } + + inline const size_t size() const + { + return _NumberOfValues; + } + +private: + // row-major representation + Scalars _rmValues; + Indizes _rmOuterIndizes; // outer index (holds rowIdx if CRS, colIdx if CCS) + Indizes _rmInnerIndizes; // inner index + size_t _rmOuterSize; + size_t _rmInnerSize; + + // common + size_t _NumberOfValues; + + size_t _rows, _cols; + + Scalar _zero; // for returning reference to values + +}; + +template +class SparseMatrixCM +{ +public: + + typedef Eigen::Matrix weightMatrix; + + typedef typename std::vector Scalars; + typedef std::vector Indizes; + + SparseMatrixCM() : + _cmOuterSize(0), _cmInnerSize(0), _NumberOfValues(0), _rows(0), _cols(0), _zero(0.) + { + } + + SparseMatrixCM(const weightMatrix & m) : + _cmOuterSize(0), _cmInnerSize(0), _NumberOfValues(0), _rows(0), _cols(0), _zero(0.) + { + constructFromEigenMatrix(m); + +#ifdef DEBUG + LOG_SPARSEMATRIX << "input matrix: " << std::endl << m << std::endl; + + LOG_SPARSEMATRIX << "constructing from Eigen matrix ..." << std::endl; + + LOG_SPARSEMATRIX << "using " << (_CompressionType ? "row" : "col") << " major representation! " << std::endl; + + LOG_SPARSEMATRIX << "values = "; + for ( size_t i = 0; i < _NumberOfValues; ++i) __PROJECT_LOG << _rmValues[i] << " "; + __PROJECT_LOG << "\n"; + + LOG_SPARSEMATRIX << "InnerIndizes = "; + for ( size_t i = 0; i < _NumberOfValues; ++i) __PROJECT_LOG << _rmInnerIndizes[i] << " "; + __PROJECT_LOG << "\n"; + + LOG_SPARSEMATRIX << "OuterIndizes = "; + for ( size_t i = 0; i < _rmOuterSize; ++i) __PROJECT_LOG << _rmOuterIndizes[i] << " "; + __PROJECT_LOG << "\n"; + +#endif + } + + virtual ~SparseMatrixCM() + { + + } + + /** + * eigen like operator for accessing the matrix-values + * @param r row + * @param c column + * @return value of (r,c) + */ + const Scalar & operator()(const size_t r, const size_t c) const + { + assert(r >= 0 && r < _rows && c >= 0 && c < _cols); + + // get end-index of outer indizes + for (size_t i = _cmOuterIndizes[c]; + i < _cmOuterIndizes[c + 1]; ++i) + { + if (_cmInnerIndizes[i] < r) + continue; + if (_cmInnerIndizes[i] == r) + return _cmValues[i]; + if (_cmInnerIndizes[i] > r) + return _zero; + } + + return _zero; + } + +private: + + void constructFromEigenMatrix(const weightMatrix & m) + { + _rows = m.rows(); + _cols = m.cols(); + + // count number of non-zero values + for (size_t r = 0; r < _rows; ++r) + for (size_t c = 0; c < _cols; ++c) + if (abs(m(r, c) > __SPARSE_MATRIX_EPSILON)) + _NumberOfValues++; + + // determine array sizes dependend on compression type + _cmOuterSize = _cols + 1; + _cmInnerSize = _rows; + + // allocate values and indizes + _cmOuterIndizes = Indizes(_cmOuterSize); + _cmInnerIndizes = Indizes(_NumberOfValues); + _cmValues = Scalars(_NumberOfValues); + + // set allocated index-array to 0 + for (auto & i : _cmOuterIndizes) + i = 0; + + // last position holds the size of values/innerindizes + _cmOuterIndizes[_cmOuterSize - 1] = _NumberOfValues; + + // read matrix data for row-major representation + size_t valuesSoFar = 0; + + for (size_t o = 0; o < _cmOuterSize - 1; ++o) + { + size_t valuesPerOuter = 0; + for (size_t i = 0; i < _cmInnerSize; ++i) + { + const Scalar mVal = m(i, o); + if (mVal != 0) + { + _cmValues[valuesSoFar] = mVal; + _cmInnerIndizes[valuesSoFar] = i; + valuesPerOuter++; + valuesSoFar++; + } + } + _cmOuterIndizes[o + 1] = _cmOuterIndizes[o] + valuesPerOuter; + } + } + +public: + + inline const Scalar maxCoeff() const + { + Scalar c = 0.; + for (size_t i = 0; i < _NumberOfValues; ++i) + c = (_cmValues[i] > c) ? _cmValues[i] : c; + return c; + } + + inline const Scalar minNonZeroCoeff() const + { + assert(_cmValues.size() > 0); + Scalar c = _cmValues[0]; + for (size_t i = 0; i < _NumberOfValues; ++i) + c = (_cmValues[i] < c) ? _cmValues[i] : c; + return c; + } + + inline const size_t & outerIndex(const size_t o) const + { + return _cmOuterIndizes[o]; + } + inline const size_t & innerIndex(const size_t i) const + { + return _cmInnerIndizes[i]; + } + inline const Scalar & value(const size_t i) const + { + return _cmValues[i]; + } + + inline Scalars & values() + { + return _cmValues; + } + inline Indizes & outerIndizes() + { + return _cmOuterIndizes; + } + inline Indizes & innerIndizes() + { + return _cmInnerIndizes; + } + + inline const size_t & rows() const + { + return _rows; + } + inline const size_t & cols() const + { + return _cols; + } + + inline const size_t size() const + { + return _NumberOfValues; + } + +private: + // col-major representation + Scalars _cmValues; + Indizes _cmOuterIndizes; // outer index (holds rowIdx if CRS, colIdx if CCS) + Indizes _cmInnerIndizes; // inner index + size_t _cmOuterSize; + size_t _cmInnerSize; + + // common + size_t _NumberOfValues; + + size_t _rows, _cols; + + Scalar _zero; // for returning reference to values + +}; + +} +#endif /* SPARSEMATRIX_H_ */ diff --git a/MatrixHelper.h b/MatrixHelper.h new file mode 100644 index 0000000..99ee3d8 --- /dev/null +++ b/MatrixHelper.h @@ -0,0 +1,88 @@ +/* + * MatrixHelper.h + * + * Created on: Jun 30, 2013 + * Author: fb + */ + +#ifndef MATRIXHELPER_H_ +#define MATRIXHELPER_H_ + +#include "internal.h" +namespace LSAP +{ + +template +class MatrixHelper +{ +private: + MatrixHelper() {} + virtual ~MatrixHelper() {} +public: + + typedef typename Eigen::Matrix WeightMatrix; + + /** + * calculate objective function value (sum over x_ij * c_ij) + * @param w weight matrix + * @param s assignment matrix + * @return value of objective function + */ + static const Scalar objFuncValue(const WeightMatrix & w, const AssignmentMatrix & s) { + Scalar ret = 0.; + for (unsigned int i = 0; i < w.rows(); ++i) + for (unsigned int j = 0; j < w.cols(); ++j) + ret += w(i, j) * s(i, j); + return ret; + } + + static const bool validAssignment(const WeightMatrix & w, const AssignmentMatrix & s) + { + for (unsigned int i = 0; i < w.rows(); ++i) + for (unsigned int j = 0; j < w.cols(); ++j) + if ( w(i, j) == 0. && s(i, j) == 1 ) return false; + return true; + } + + /** + * calculate objective function value (sum over x_ij * c_ij) + * @param w weight matrix + * @param e array of edges (assumed to be of capacity rows) + * @return value of objective function + */ + static const Scalar objFuncValue(const WeightMatrix & w, Edge * e) { + Scalar ret = 0.; + for (unsigned int i = 0; i < w.rows(); ++i) + ret += w(e[i].x, e[i].y); + return ret; + } + + /** + * calculate objective function value (sum over x_ij * c_ij) + * @param w weight matrix + * @param e vector of edges + * @return value of objective function + */ + static const Scalar objFuncValue(const WeightMatrix & w, const Edges & e) { + Scalar ret = 0.; + for (unsigned int i = 0; i < e.size(); ++i) + ret += w(e[i].x, e[i].y); + + return ret; + } + + /** + * calculate objective function value (sum over x_ij * c_ij) + * @param e vector of edges with contained values + * @return value of objective function + */ + static const Scalar objFuncValue(const Edges & edges) { + Scalar ret = 0.; + for (auto & e : edges) ret += e.v; + + return ret; + } +}; + +} /* namespace LSAP */ +#endif /* MATRIXHELPER_H_ */ diff --git a/ThreadBarrier.h b/ThreadBarrier.h new file mode 100644 index 0000000..eb89475 --- /dev/null +++ b/ThreadBarrier.h @@ -0,0 +1,80 @@ +/* + * ThreadBarrier.h + * + * Created on: Jul 3, 2013 + * Author: fb + */ + +#ifndef THREADBARRIER_H_ +#define THREADBARRIER_H_ +#include + +namespace LSAP +{ +/** + * realizes a barrier with thread kill for std:threads + */ +class ThreadBarrier +{ +public: + /** + * creates a barrier which waits for n threads for continuation + * @param n number of threads + */ + ThreadBarrier(const size_t n, const bool & purespin = false, const size_t sleeptime = 5) : + _n(n), _nwait(0), _step(0), _kill(false), _spin(purespin), _sleep(true), _sleeptime(sleeptime) + { + } + virtual ~ThreadBarrier() + { + } + + void sleep() { if ( !_spin) _sleep = true; } + + void wakeup() { if ( !_spin) _sleep = false; } + + const bool operator()() + { + // return true if kill-signal + if (_kill.load()) + return true; + + unsigned int step = _step.load(); + + if (_nwait.fetch_add(1) == _n - 1) + { + _nwait.store(0); + _step.fetch_add(1); + } + else + { + // waiting + while (_step.load() == step) + { + if (_sleep) usleep(_sleeptime); +// return true if kill-signal + if (_kill.load()) + return true; + } + } + return false; + } + + /** + * set kill flag for return value + */ + void kill() + { + _kill.store(true); + } +protected: + size_t _n; + std::atomic _nwait; + std::atomic _step; + std::atomic _kill; + bool _spin; + bool _sleep; + size_t _sleeptime; +}; +} /* namespace LSAP */ +#endif /* THREADBARRIER_H_ */ diff --git a/internal.h b/internal.h new file mode 100644 index 0000000..3dc8e6a --- /dev/null +++ b/internal.h @@ -0,0 +1,74 @@ +/* + * internal.h + * + * Created on: Apr 25, 2013 + * Author: fb + */ + +#ifndef INTERNAL_H_ +#define INTERNAL_H_ + +#define __PROJECT_LOG std::cerr +#include +#include +#include +#include +#include "Matrix/SparseMatrix.h" +#define NO_CL +#define __AUCTION_EPSILON_MULTIPLIER 1e-5 // epsilon multiplier +#define __AUCTION_INF 1e6 // infinity for setting second best match + +typedef Eigen::Matrix AssignmentMatrix; + +// which type of coefficient in weight matrix +typedef double Scalar; + +// represents edges/assignments +typedef std::pair Assignment; +typedef std::vector Assignments; + +struct Edge +{ + size_t x; + size_t y; + Scalar v; +}; + +typedef std::vector Edges; + +// locking of rows/cols +typedef std::vector LockVector; + +enum APSolvingMode { MAXIMIZATION, MINIMIZATION }; + + +#define MEASURE_DURATION_SINGLE(__MD_COMMAND__) \ +{\ +double __md_accumulator = 0;\ +std::cout << "command <" << #__MD_COMMAND__ << "> took "; \ +boost::chrono::high_resolution_clock::time_point __md_start = boost::chrono::high_resolution_clock::now();\ +__MD_COMMAND__ ;\ +boost::chrono::nanoseconds __md_sec = boost::chrono::high_resolution_clock::now() - __md_start;\ +__md_accumulator += __md_sec.count();\ +std::cout << ( __md_accumulator / 1000. ) << " us " << std::endl;} \ + +#define MEASURE_DURATION_SINGLE_WITHOUT_PRINTING_COMMAND(__MD_COMMAND__) \ +{\ +double __md_accumulator = 0;\ +boost::chrono::high_resolution_clock::time_point __md_start = boost::chrono::high_resolution_clock::now();\ +__MD_COMMAND__ ;\ +boost::chrono::nanoseconds __md_sec = boost::chrono::high_resolution_clock::now() - __md_start;\ +__md_accumulator += __md_sec.count();\ +std::cout << ( __md_accumulator / 1000. ) << " us " << std::endl;} \ + + +#define MEASURE_DURATION_SINGLE_STORED(__MD_COMMAND__, __MD_STORE__) \ +{\ +double __md_accumulator = 0;\ +boost::chrono::high_resolution_clock::time_point __md_start = boost::chrono::high_resolution_clock::now();\ +__MD_COMMAND__ ;\ +boost::chrono::nanoseconds __md_sec = boost::chrono::high_resolution_clock::now() - __md_start;\ +__md_accumulator += __md_sec.count();\ +__MD_STORE__ = (double) ( __md_accumulator / 1000. );} \ + +#endif /* INTERNAL_H_ */