diff options
Diffstat (limited to 'geom_matching/wasserstein')
19 files changed, 4175 insertions, 0 deletions
diff --git a/geom_matching/wasserstein/example/wasserstein_dist.cpp b/geom_matching/wasserstein/example/wasserstein_dist.cpp new file mode 100644 index 0000000..e92ab54 --- /dev/null +++ b/geom_matching/wasserstein/example/wasserstein_dist.cpp @@ -0,0 +1,89 @@ +/* + +Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#include <iostream> +#include <iomanip> +#include <fstream> +#include <vector> +#include <algorithm> +#include <limits> +#include <random> + +#include "wasserstein.h" + +// any container of pairs of doubles can be used, +// we use vector in this example. + +typedef std::vector<std::pair<double, double>> PairVector; + +int main(int argc, char* argv[]) +{ + PairVector diagramA, diagramB; + + if (argc < 3 ) { + std::cerr << "Usage: " << argv[0] << " file1 file2 [wasserstein_degree] [relative_error] [internal norm]. By default power is 1.0, relative error is 0.01, internal norm is l_infinity." << std::endl; + return 1; + } + + if (!geom_ws::readDiagramPointSet(argv[1], diagramA)) { + std::exit(1); + } + + if (!geom_ws::readDiagramPointSet(argv[2], diagramB)) { + std::exit(1); + } + + double wasserPower = (4 <= argc) ? atof(argv[3]) : 1.0; + if (wasserPower < 1.0) { + std::cerr << "The third argument (wasserstein_degree) was \"" << argv[3] << "\", must be a number >= 1.0. Cannot proceed. " << std::endl; + std::exit(1); + } + + //default relative error: 1% + double delta = (5 <= argc) ? atof(argv[4]) : 0.01; + if ( delta <= 0.0) { + std::cerr << "The 4th argument (relative error) was \"" << argv[4] << "\", must be a number > 0.0. Cannot proceed. " << std::endl; + std::exit(1); + } + + // default for internal metric is l_infinity + double internal_p = ( 6 <= argc ) ? atof(argv[5]) : std::numeric_limits<double>::infinity(); + if (internal_p < 1.0) { + std::cerr << "The 5th argument (internal norm) was \"" << argv[5] << "\", must be a number >= 1.0. Cannot proceed. " << std::endl; + std::exit(1); + } + + // if you want to specify initial value for epsilon and the factor + // for epsilon-scaling + double initialEpsilon= ( 7 <= argc ) ? atof(argv[6]) : 0.0 ; + double epsFactor = ( 8 <= argc ) ? atof(argv[7]) : 0.0 ; + + double res = geom_ws::wassersteinDist(diagramA, diagramB, wasserPower, delta, internal_p, initialEpsilon, epsFactor); + std::cout << std::setprecision(15) << res << std::endl; + return 0; +} diff --git a/geom_matching/wasserstein/include/auction_oracle.h b/geom_matching/wasserstein/include/auction_oracle.h new file mode 100644 index 0000000..e803218 --- /dev/null +++ b/geom_matching/wasserstein/include/auction_oracle.h @@ -0,0 +1,305 @@ +/* + +Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#ifndef AUCTION_ORACLE_H +#define AUCTION_ORACLE_H + + +#define USE_BOOST_HEAP + +#include <map> +#include <memory> +#include <set> +#include <list> + +#ifdef USE_BOOST_HEAP +#include <boost/heap/d_ary_heap.hpp> +#endif + +#include "basic_defs_ws.h" +#include "dnn/geometry/euclidean-fixed.h" +#include "dnn/local/kd-tree.h" + +namespace geom_ws { + +struct CompPairsBySecondStruct { + bool operator()(const IdxValPair& a, const IdxValPair& b) const + { + return a.second < b.second; + } +}; + +// +struct CompPairsBySecondGreaterStruct { + bool operator()(const IdxValPair& a, const IdxValPair& b) const + { + return a.second > b.second; + } +}; + +struct CompPairsBySecondLexStruct { + bool operator()(const IdxValPair& a, const IdxValPair& b) const + { + return a.second < b.second or (a.second == b.second and a.first > b.first); + } +}; + +struct CompPairsBySecondLexGreaterStruct { + bool operator()(const IdxValPair& a, const IdxValPair& b) const + { + return a.second > b.second or (a.second == b.second and a.first > b.first); + } +}; + +using ItemsTimePair = std::pair<IdxType, int>; + +using UpdateList = std::list<ItemsTimePair>; +using UpdateListIter = UpdateList::iterator; + + +#ifdef USE_BOOST_HEAP +using LossesHeap = boost::heap::d_ary_heap<IdxValPair, boost::heap::arity<2>, boost::heap::mutable_<true>, boost::heap::compare<CompPairsBySecondGreaterStruct>>; +#else +template<class ComparisonStruct> +class IdxValHeap { +public: + using InternalKeeper = std::set<IdxValPair, ComparisonStruct>; + using handle_type = typename InternalKeeper::iterator; + // methods + handle_type push(const IdxValPair& val) + { + auto resPair = _heap.insert(val); + assert(resPair.second); + assert(resPair.first != _heap.end()); + return resPair.first; + } + + void decrease(handle_type& handle, const IdxValPair& newVal) + { + _heap.erase(handle); + handle = push(newVal); + } + + size_t size() const + { + return _heap.size(); + } + + handle_type ordered_begin() + { + return _heap.begin(); + } + + handle_type ordered_end() + { + return _heap.end(); + } + + +private: + std::set<IdxValPair, ComparisonStruct> _heap; +}; + +// if we store losses, the minimal value should come first +using LossesHeap = IdxValHeap<CompPairsBySecondLexStruct>; +#endif + +struct DebugOptimalBid { + DebugOptimalBid() : bestItemIdx(-1), bestItemValue(-666.666), secondBestItemIdx(-1), secondBestItemValue(-666.666) {}; + IdxType bestItemIdx; + double bestItemValue; + IdxType secondBestItemIdx; + double secondBestItemValue; +}; + +struct AuctionOracleAbstract { + AuctionOracleAbstract(const std::vector<DiagramPoint>& _bidders, const std::vector<DiagramPoint>& _items, const double _wassersteinPower, const double _internal_p = std::numeric_limits<double>::infinity()); + ~AuctionOracleAbstract() {} + virtual IdxValPair getOptimalBid(const IdxType bidderIdx) = 0; + virtual void setPrice(const IdxType itemsIdx, const double newPrice) = 0; + virtual void adjustPrices(void) = 0; + double getEpsilon() { return epsilon; }; + virtual void setEpsilon(double newEpsilon) { assert(newEpsilon >= 0.0); epsilon = newEpsilon; }; + std::vector<double> getPrices() { return prices; } +protected: + const std::vector<DiagramPoint>& bidders; + const std::vector<DiagramPoint>& items; + std::vector<double> prices; + double wassersteinPower; + double epsilon; + double internal_p; + double getValueForBidder(size_t bidderIdx, size_t itemsIdx); +}; + +struct AuctionOracleLazyHeap final : AuctionOracleAbstract { + AuctionOracleLazyHeap(const std::vector<DiagramPoint>& bidders, const std::vector<DiagramPoint>& items, const double wassersteinPower, const double _internal_p = std::numeric_limits<double>::infinity()); + ~AuctionOracleLazyHeap(); + // data members + // temporarily make everything public + std::vector<std::vector<double>> weightMatrix; + //double weightAdjConst; + double maxVal; + // vector of heaps to find the best items + std::vector<LossesHeap*> lossesHeap; + std::vector<std::vector<LossesHeap::handle_type>> lossesHeapHandles; + // methods + void fillInLossesHeap(void); + void setPrice(const IdxType itemsIdx, const double newPrice) override final; + IdxValPair getOptimalBid(const IdxType bidderIdx) override final; + double getMatchingWeight(const std::vector<IdxType>& biddersToItems) const; + void adjustPrices(void) override final; + // to update the queue in lazy fashion + std::vector<UpdateListIter> itemsIterators; + UpdateList updateList; + std::vector<int> biddersUpdateMoments; + int updateCounter; + void updateQueueForBidder(const IdxType bidderIdx); + // debug + DebugOptimalBid getOptimalBidDebug(const IdxType bidderIdx); +}; + +struct AuctionOracleLazyHeapRestricted final : AuctionOracleAbstract { + AuctionOracleLazyHeapRestricted(const std::vector<DiagramPoint>& bidders, const std::vector<DiagramPoint>& items, const double wassersteinPower, const double _internal_p = std::numeric_limits<double>::infinity()); + ~AuctionOracleLazyHeapRestricted(); + // data members + // temporarily make everything public + std::vector<std::vector<double>> weightMatrix; + //double weightAdjConst; + double maxVal; + // vector of heaps to find the best items + std::vector<LossesHeap*> lossesHeap; + std::vector<std::vector<size_t>> itemsIndicesForHeapHandles; + std::vector<std::vector<LossesHeap::handle_type>> lossesHeapHandles; + // methods + void fillInLossesHeap(void); + void setPrice(const IdxType itemsIdx, const double newPrice) override final; + IdxValPair getOptimalBid(const IdxType bidderIdx) override final; + double getMatchingWeight(const std::vector<IdxType>& biddersToItems) const; + void adjustPrices(void) override final; + // to update the queue in lazy fashion + std::vector<UpdateListIter> itemsIterators; + UpdateList updateList; + std::vector<int> biddersUpdateMoments; + int updateCounter; + void updateQueueForBidder(const IdxType bidderIdx); + LossesHeap diagItemsHeap; + std::vector<LossesHeap::handle_type> diagHeapHandles; + std::vector<size_t> heapHandlesIndices; + // debug + + DebugOptimalBid getOptimalBidDebug(const IdxType bidderIdx); + + // for diagonal points + bool bestDiagonalItemsComputed; + size_t bestDiagonalItemIdx; + double bestDiagonalItemValue; + size_t secondBestDiagonalItemIdx; + double secondBestDiagonalItemValue; +}; + +struct AuctionOracleKDTree final : AuctionOracleAbstract { + typedef dnn::Point<2, double> DnnPoint; + typedef dnn::PointTraits<DnnPoint> DnnTraits; + + AuctionOracleKDTree(const std::vector<DiagramPoint>& bidders, const std::vector<DiagramPoint>& items, const double wassersteinPower, const double _internal_p = std::numeric_limits<double>::infinity()); + ~AuctionOracleKDTree(); + // data members + // temporarily make everything public + double maxVal; + double weightAdjConst; + dnn::KDTree<DnnTraits>* kdtree; + std::vector<DnnPoint> dnnPoints; + std::vector<DnnPoint*> dnnPointHandles; + dnn::KDTree<DnnTraits>* kdtreeAll; + std::vector<DnnPoint> dnnPointsAll; + std::vector<DnnPoint*> dnnPointHandlesAll; + LossesHeap diagItemsHeap; + std::vector<LossesHeap::handle_type> diagHeapHandles; + std::vector<size_t> heapHandlesIndices; + std::vector<size_t> kdtreeItems; + // vector of heaps to find the best items + void setPrice(const IdxType itemsIdx, const double newPrice) override final; + IdxValPair getOptimalBid(const IdxType bidderIdx) override final; + void adjustPrices(void) override final; + // debug routines + DebugOptimalBid getOptimalBidDebug(IdxType bidderIdx); + void setEpsilon(double newVal) override final; +}; + +struct AuctionOracleKDTreeRestricted final : AuctionOracleAbstract { + typedef dnn::Point<2, double> DnnPoint; + typedef dnn::PointTraits<DnnPoint> DnnTraits; + + AuctionOracleKDTreeRestricted(const std::vector<DiagramPoint>& bidders, const std::vector<DiagramPoint>& items, const double wassersteinPower, const double _internal_p = std::numeric_limits<double>::infinity()); + ~AuctionOracleKDTreeRestricted(); + // data members + // temporarily make everything public + double maxVal; + double weightAdjConst; + dnn::KDTree<DnnTraits>* kdtree; + std::vector<DnnPoint> dnnPoints; + std::vector<DnnPoint*> dnnPointHandles; + std::vector<DnnPoint> dnnPointsAll; + std::vector<DnnPoint*> dnnPointHandlesAll; + LossesHeap diagItemsHeap; + std::vector<LossesHeap::handle_type> diagHeapHandles; + std::vector<size_t> heapHandlesIndices; + std::vector<size_t> kdtreeItems; + // vector of heaps to find the best items + void setPrice(const IdxType itemsIdx, const double newPrice) override final; + IdxValPair getOptimalBid(const IdxType bidderIdx) override final; + void adjustPrices(void) override final; + // debug routines + DebugOptimalBid getOptimalBidDebug(IdxType bidderIdx); + void setEpsilon(double newVal) override final; + + + bool bestDiagonalItemsComputed; + size_t bestDiagonalItemIdx; + double bestDiagonalItemValue; + size_t secondBestDiagonalItemIdx; + double secondBestDiagonalItemValue; +}; + +struct AuctionOracleRestricted final : AuctionOracleAbstract { + AuctionOracleRestricted(const std::vector<DiagramPoint>& bidders, const std::vector<DiagramPoint>& items, const double wassersteinPower, const double _internal_p = std::numeric_limits<double>::infinity()); + IdxValPair getOptimalBid(const IdxType bidderIdx) override; + void setPrice(const IdxType itemsIdx, const double newPrice) override; + void adjustPrices(void) override {}; + void setEpsilon(double newEpsilon) override { assert(newEpsilon >= 0.0); epsilon = newEpsilon; }; + // data + std::vector<std::vector<double>> weightMatrix; + double maxVal; + constexpr static bool isRestricted = true; +}; + +std::ostream& operator<< (std::ostream& output, const DebugOptimalBid& db); + +} // end of namespace geom_ws + +#endif diff --git a/geom_matching/wasserstein/include/auction_runner_gs.h b/geom_matching/wasserstein/include/auction_runner_gs.h new file mode 100644 index 0000000..34a91e8 --- /dev/null +++ b/geom_matching/wasserstein/include/auction_runner_gs.h @@ -0,0 +1,122 @@ +/* + +Copyright (c) 2016, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#ifndef AUCTION_RUNNER_GS_H +#define AUCTION_RUNNER_GS_H + +#include <unordered_set> + +#include "auction_oracle.h" + +//#define KEEP_UNASSIGNED_ORDERED +// if this symbol is defined, +// unassigned bidders are processed in a lexicographic order. +// See LexicogrCompDiagramPoint comparator. + + +namespace geom_ws { + +//using AuctionOracle = AuctionOracleLazyHeapRestricted; +using AuctionOracle = AuctionOracleKDTreeRestricted; + +#ifdef KEEP_UNASSIGNED_ORDERED +using IdxPointPair = std::pair<size_t, DiagramPoint>; + +struct LexicogrCompDiagramPoint { + bool operator ()(const IdxPointPair& a, const IdxPointPair& b) { + const auto& p1 = a.second; + const auto& p2 = b.second; + + return ( (not p1.isDiagonal() and p2.isDiagonal()) or + ( p1.isDiagonal() == p2.isDiagonal() and p1.getRealX() < p2.getRealX() ) or + ( p1.isDiagonal() == p2.isDiagonal() and p1.getRealX() == p2.getRealX() and p1.getRealY() < p2.getRealY() ) or + ( p1.isDiagonal() == p2.isDiagonal() and p1.getRealX() == p2.getRealX() and p1.getRealY() == p2.getRealY() and a.first < b.first ) ); + } +}; + +using OrderedUnassignedKeeper = std::set<IdxPointPair, LexicogrCompDiagramPoint>; +#endif + +// the two parameters that you can tweak in auction algorithm are: +// 1. epsilonCommonRatio +// 2. maxIterNum + +class AuctionRunnerGS { +public: + AuctionRunnerGS(const std::vector<DiagramPoint>& A, + const std::vector<DiagramPoint>& B, + const double q, + const double _delta, + const double _internal_p, + const double _initialEpsilon, + const double _epsFactor); + void setEpsilon(double newVal) { assert(epsilon > 0.0); epsilon = newVal; }; + double getEpsilon(void) const { return epsilon; } + double getWassersteinDistance(void); + static constexpr int maxIterNum { 25 }; // maximal number of iterations of epsilon-scaling +private: + // private data + std::vector<DiagramPoint> bidders, items; + const size_t numBidders; + const size_t numItems; + std::vector<IdxType> itemsToBidders; + std::vector<IdxType> biddersToItems; + double wassersteinPower; + double epsilon; + double delta; + double internal_p; + double initialEpsilon; + double epsilonCommonRatio; // next epsilon = current epsilon / epsilonCommonRatio + double weightAdjConst; + double wassersteinDistance; + // to get the 2 best items + std::unique_ptr<AuctionOracle> oracle; +#ifdef KEEP_UNASSIGNED_ORDERED + OrderedUnassignedKeeper unassignedBidders; +#else + std::unordered_set<size_t> unassignedBidders; +#endif + // private methods + void assignItemToBidder(const IdxType bidderIdx, const IdxType itemsIdx); + void clearBidTable(void); + void runAuction(void); + void runAuctionPhase(void); + void flushAssignment(void); + + // for debug only + void sanityCheck(void); + void printDebug(void); + int countUnhappy(void); + void printMatching(void); + double getDistanceToQthPowerInternal(void); + int numRounds { 0 }; +}; + +} // end of namespace geom_ws + +#endif diff --git a/geom_matching/wasserstein/include/auction_runner_jac.h b/geom_matching/wasserstein/include/auction_runner_jac.h new file mode 100644 index 0000000..ae0cb56 --- /dev/null +++ b/geom_matching/wasserstein/include/auction_runner_jac.h @@ -0,0 +1,97 @@ +/* + +Copyright (c) 2016, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#ifndef AUCTION_RUNNER_JAK_H +#define AUCTION_RUNNER_JAK_H + +#include <unordered_set> + +#include "auction_oracle.h" + +namespace geom_ws { + +using AuctionOracle = AuctionOracleKDTreeRestricted; + +// the two parameters that you can tweak in auction algorithm are: +// 1. epsilonCommonRatio +// 2. maxIterNum + +// the two parameters that you can tweak in auction algorithm are: +// 1. epsilonCommonRatio +// 2. maxIterNum + +class AuctionRunnerJak { +public: + AuctionRunnerJak(const std::vector<DiagramPoint>& A, const std::vector<DiagramPoint>& B, const double q, const double _delta, const double _internal_p); + void setEpsilon(double newVal) { assert(epsilon > 0.0); epsilon = newVal; }; + double getEpsilon(void) const { return epsilon; } + double getWassersteinDistance(void); + static constexpr double epsilonCommonRatio { 5 }; // next epsilon = current epsilon / epsilonCommonRatio + static constexpr int maxIterNum { 25 }; // maximal number of iterations of epsilon-scaling +private: + // private data + std::vector<DiagramPoint> bidders, items; + const size_t numBidders; + const size_t numItems; + std::vector<IdxType> itemsToBidders; + std::vector<IdxType> biddersToItems; + double wassersteinPower; + double epsilon; + double delta; + double internal_p; + double weightAdjConst; + double wassersteinDistance; + std::vector<IdxValPair> bidTable; + // to get the 2 best items + std::unique_ptr<AuctionOracle> oracle; + std::list<size_t> unassignedBidders; + std::vector< std::list<size_t>::iterator > unassignedBiddersIterators; + std::vector< short > itemReceivedBidVec; + std::list<size_t> itemsWithBids; + // private methods + void assignGoodToBidder(const IdxType bidderIdx, const IdxType itemsIdx); + void assignToBestBidder(const IdxType itemsIdx); + void clearBidTable(void); + void runAuction(void); + void runAuctionPhase(void); + void submitBid(IdxType bidderIdx, const IdxValPair& itemsBidValuePair); + void flushAssignment(void); + + // for debug only + void sanityCheck(void); + void printDebug(void); + int countUnhappy(void); + void printMatching(void); + double getDistanceToQthPowerInternal(void); +}; + + + +} // end of namespace geom_ws + +#endif diff --git a/geom_matching/wasserstein/include/basic_defs_ws.h b/geom_matching/wasserstein/include/basic_defs_ws.h new file mode 100644 index 0000000..20cd2a0 --- /dev/null +++ b/geom_matching/wasserstein/include/basic_defs_ws.h @@ -0,0 +1,114 @@ +/* + +Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#ifndef BASIC_DEFS_WS_H +#define BASIC_DEFS_WS_H + +#include <vector> +#include <math.h> +#include <cstddef> +#include <unordered_map> +#include <unordered_set> +#include <iostream> +#include <string> +#include <assert.h> + +#ifdef _WIN32 +#include <ciso646> +#endif + + +#include "def_debug.h" + +#define MIN_VALID_ID 10 + +namespace geom_ws { + +using IdxType = int; +using IdxValPair = std::pair<IdxType, double>; + + +struct Point { + double x, y; + bool operator==(const Point& other) const; + bool operator!=(const Point& other) const; + Point(double ax, double ay) : x(ax), y(ay) {} + Point() : x(0.0), y(0.0) {} + friend std::ostream& operator<<(std::ostream& output, const Point p); +}; + +struct DiagramPoint +{ + // data members + // Points above the diagonal have type NORMAL + // Projections onto the diagonal have type DIAG + // for DIAG points only x-coordinate is relevant + enum Type { NORMAL, DIAG}; + double x, y; + Type type; + // methods + DiagramPoint(double xx, double yy, Type ttype); + bool isDiagonal(void) const { return type == DIAG; } + bool isNormal(void) const { return type == NORMAL; } + double getRealX() const; // return the x-coord + double getRealY() const; // return the y-coord + friend std::ostream& operator<<(std::ostream& output, const DiagramPoint p); + + struct LexicographicCmp + { + bool operator()(const DiagramPoint& p1, const DiagramPoint& p2) const + { return p1.type < p2.type || (p1.type == p2.type && (p1.x < p2.x || (p1.x == p2.x && p1.y < p2.y))); } + }; +}; + +double sqrDist(const Point& a, const Point& b); +double dist(const Point& a, const Point& b); +double distLInf(const DiagramPoint& a, const DiagramPoint& b); +double distLp(const DiagramPoint& a, const DiagramPoint& b, const double p); + +template<typename DiagPointContainer> +double getFurthestDistance3Approx(DiagPointContainer& A, DiagPointContainer& B) +{ + double result { 0.0 }; + DiagramPoint begA = *(A.begin()); + DiagramPoint optB = *(B.begin()); + for(const auto& pointB : B) { + if (distLInf(begA, pointB) > result) { + result = distLInf(begA, pointB); + optB = pointB; + } + } + for(const auto& pointA : A) { + if (distLInf(pointA, optB) > result) { + result = distLInf(pointA, optB); + } + } + return result; +} + +} // end of namespace geom_ws +#endif diff --git a/geom_matching/wasserstein/include/def_debug.h b/geom_matching/wasserstein/include/def_debug.h new file mode 100644 index 0000000..7323c18 --- /dev/null +++ b/geom_matching/wasserstein/include/def_debug.h @@ -0,0 +1,36 @@ +/* + +Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#ifndef DEF_DEBUG_H +#define DEF_DEBUG_H + +//#define DEBUG_BOUND_MATCH +//#define DEBUG_NEIGHBOUR_ORACLE +//#define DEBUG_MATCHING +//#define DEBUG_AUCTION + +#endif diff --git a/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h b/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h new file mode 100644 index 0000000..a6ccef7 --- /dev/null +++ b/geom_matching/wasserstein/include/dnn/geometry/euclidean-fixed.h @@ -0,0 +1,190 @@ +#ifndef DNN_GEOMETRY_EUCLIDEAN_FIXED_H +#define DNN_GEOMETRY_EUCLIDEAN_FIXED_H + +#include <boost/operators.hpp> +#include <boost/array.hpp> +#include <boost/range/value_type.hpp> +#include <boost/serialization/access.hpp> +#include <boost/serialization/base_object.hpp> + +#include <iostream> +#include <fstream> +#include <string> +#include <sstream> +#include <cmath> + +#include "../parallel/tbb.h" // for dnn::vector<...> + +namespace dnn +{ + // TODO: wrap in another namespace (e.g., euclidean) + + template<size_t D, typename Real = double> + struct Point: + boost::addable< Point<D,Real>, + boost::subtractable< Point<D,Real>, + boost::dividable2< Point<D, Real>, Real, + boost::multipliable2< Point<D, Real>, Real > > > >, + public boost::array<Real, D> + { + public: + typedef Real Coordinate; + typedef Real DistanceType; + + + public: + Point(size_t id = 0): id_(id) {} + template<size_t DD> + Point(const Point<DD,Real>& p, size_t id = 0): + id_(id) { *this = p; } + + static size_t dimension() { return D; } + + // Assign a point of different dimension + template<size_t DD> + Point& operator=(const Point<DD,Real>& p) { for (size_t i = 0; i < (D < DD ? D : DD); ++i) (*this)[i] = p[i]; if (DD < D) for (size_t i = DD; i < D; ++i) (*this)[i] = 0; return *this; } + + Point& operator+=(const Point& p) { for (size_t i = 0; i < D; ++i) (*this)[i] += p[i]; return *this; } + Point& operator-=(const Point& p) { for (size_t i = 0; i < D; ++i) (*this)[i] -= p[i]; return *this; } + Point& operator/=(Real r) { for (size_t i = 0; i < D; ++i) (*this)[i] /= r; return *this; } + Point& operator*=(Real r) { for (size_t i = 0; i < D; ++i) (*this)[i] *= r; return *this; } + + Real norm2() const { Real n = 0; for (size_t i = 0; i < D; ++i) n += (*this)[i] * (*this)[i]; return n; } + Real max_norm() const + { + Real res = std::fabs((*this)[0]); + for (size_t i = 1; i < D; ++i) + if (std::fabs((*this)[i]) > res) + res = std::fabs((*this)[i]); + return res; + } + + Real l1_norm() const + { + Real res = std::fabs((*this)[0]); + for (size_t i = 1; i < D; ++i) + res += std::fabs((*this)[i]); + return res; + } + + Real lp_norm(const Real p) const + { + assert( !std::isinf(p) ); + if ( p == 1.0 ) + return l1_norm(); + Real res = std::pow(std::fabs((*this)[0]), p); + for (size_t i = 1; i < D; ++i) + res += std::pow(std::fabs((*this)[i]), p); + return std::pow(res, 1.0 / p); + } + + // quick and dirty for now; make generic later + //DistanceType distance(const Point& other) const { return sqrt(sq_distance(other)); } + //DistanceType sq_distance(const Point& other) const { return (other - *this).norm2(); } + + DistanceType distance(const Point& other) const { return (other - *this).max_norm(); } + DistanceType p_distance(const Point& other, const double p) const { return (other - *this).lp_norm(p); } + + size_t id() const { return id_; } + size_t& id() { return id_; } + + private: + friend class boost::serialization::access; + + template<class Archive> + void serialize(Archive& ar, const unsigned int version) { ar & boost::serialization::base_object< boost::array<Real,D> >(*this) & id_; } + + private: + size_t id_; + }; + + template<size_t D, typename Real> + std::ostream& + operator<<(std::ostream& out, const Point<D,Real>& p) + { out << p[0]; for (size_t i = 1; i < D; ++i) out << " " << p[i]; return out; } + + + template<class Point> + struct PointTraits; // intentionally undefined; should be specialized for each type + + + template<size_t D, typename Real> + struct PointTraits< Point<D, Real> > // specialization for dnn::Point + { + typedef Point<D,Real> PointType; + typedef const PointType* PointHandle; + typedef std::vector<PointType> PointContainer; + + typedef typename PointType::Coordinate Coordinate; + typedef typename PointType::DistanceType DistanceType; + + + static DistanceType + distance(const PointType& p1, const PointType& p2) { if (std::isinf(internal_p)) return p1.distance(p2); else return p1.p_distance(p2, internal_p); } + + static DistanceType + distance(PointHandle p1, PointHandle p2) { return distance(*p1,*p2); } + + static size_t dimension() { return D; } + static Real coordinate(const PointType& p, size_t i) { return p[i]; } + static Real& coordinate(PointType& p, size_t i) { return p[i]; } + static Real coordinate(PointHandle p, size_t i) { return coordinate(*p,i); } + + static size_t id(const PointType& p) { return p.id(); } + static size_t& id(PointType& p) { return p.id(); } + static size_t id(PointHandle p) { return id(*p); } + + static PointHandle + handle(const PointType& p) { return &p; } + static const PointType& + point(PointHandle ph) { return *ph; } + + void swap(PointType& p1, PointType& p2) const { return std::swap(p1, p2); } + + static PointContainer + container(size_t n = 0, const PointType& p = PointType()) { return PointContainer(n, p); } + static typename PointContainer::iterator + iterator(PointContainer& c, PointHandle ph) { return c.begin() + (ph - &c[0]); } + static typename PointContainer::const_iterator + iterator(const PointContainer& c, PointHandle ph) { return c.begin() + (ph - &c[0]); } + + // Internal_p determines which norm will be used in Wasserstein metric (not to + // be confused with wassersteinPower parameter: + // we raise \| p - q \|_{internal_p} to wassersteinPower. + static Real internal_p; + + private: + + friend class boost::serialization::access; + + template<class Archive> + void serialize(Archive& ar, const unsigned int version) {} + + }; + + template<size_t D, typename Real> + Real PointTraits< Point<D, Real> >::internal_p = std::numeric_limits<Real>::infinity(); + + + template<class PointContainer> + void read_points(const std::string& filename, PointContainer& points) + { + typedef typename boost::range_value<PointContainer>::type Point; + typedef typename PointTraits<Point>::Coordinate Coordinate; + + std::ifstream in(filename.c_str()); + std::string line; + while(std::getline(in, line)) + { + if (line[0] == '#') continue; // comment line in the file + std::stringstream linestream(line); + Coordinate x; + points.push_back(Point()); + size_t i = 0; + while (linestream >> x) + points.back()[i++] = x; + } + } +} + +#endif diff --git a/geom_matching/wasserstein/include/dnn/local/kd-tree.h b/geom_matching/wasserstein/include/dnn/local/kd-tree.h new file mode 100644 index 0000000..7e01072 --- /dev/null +++ b/geom_matching/wasserstein/include/dnn/local/kd-tree.h @@ -0,0 +1,90 @@ +#ifndef DNN_LOCAL_KD_TREE_H +#define DNN_LOCAL_KD_TREE_H + +#include "../utils.h" +#include "search-functors.h" + +#include <unordered_map> + +#include <boost/tuple/tuple.hpp> +#include <boost/shared_ptr.hpp> +#include <boost/range/value_type.hpp> + +#include <boost/static_assert.hpp> +#include <boost/type_traits.hpp> + +namespace dnn +{ + // Weighted KDTree + // Traits_ provides Coordinate, DistanceType, PointType, dimension(), distance(p1,p2), coordinate(p,i) + template< class Traits_ > + class KDTree + { + public: + typedef Traits_ Traits; + typedef dnn::HandleDistance<KDTree> HandleDistance; + + typedef typename Traits::PointType Point; + typedef typename Traits::PointHandle PointHandle; + typedef typename Traits::Coordinate Coordinate; + typedef typename Traits::DistanceType DistanceType; + typedef std::vector<PointHandle> HandleContainer; + typedef std::vector<HandleDistance> HDContainer; // TODO: use tbb::scalable_allocator + typedef HDContainer Result; + typedef std::vector<DistanceType> DistanceContainer; + typedef std::unordered_map<PointHandle, size_t> HandleMap; + + BOOST_STATIC_ASSERT_MSG(has_coordinates<Traits, PointHandle, int>::value, "KDTree requires coordinates"); + + public: + KDTree(const Traits& traits): + traits_(traits) {} + + KDTree(const Traits& traits, HandleContainer&& handles, double _wassersteinPower = 1.0); + + template<class Range> + KDTree(const Traits& traits, const Range& range, double _wassersteinPower = 1.0); + + template<class Range> + void init(const Range& range); + + DistanceType weight(PointHandle p) { return weights_[indices_[p]]; } + void increase_weight(PointHandle p, DistanceType w); + + HandleDistance find(PointHandle q) const; + Result findR(PointHandle q, DistanceType r) const; // all neighbors within r + Result findK(PointHandle q, size_t k) const; // k nearest neighbors + + HandleDistance find(const Point& q) const { return find(traits().handle(q)); } + Result findR(const Point& q, DistanceType r) const { return findR(traits().handle(q), r); } + Result findK(const Point& q, size_t k) const { return findK(traits().handle(q), k); } + + template<class ResultsFunctor> + void search(PointHandle q, ResultsFunctor& rf) const; + + const Traits& traits() const { return traits_; } + + void printWeights(void); + + private: + void init(); + + typedef typename HandleContainer::iterator HCIterator; + typedef std::tuple<HCIterator, HCIterator, size_t> KDTreeNode; + + struct CoordinateComparison; + struct OrderTree; + + private: + Traits traits_; + HandleContainer tree_; + DistanceContainer weights_; // point weight + DistanceContainer subtree_weights_; // min weight in the subtree + HandleMap indices_; + double wassersteinPower; + }; +} + +#include "kd-tree.hpp" + +#endif diff --git a/geom_matching/wasserstein/include/dnn/local/kd-tree.hpp b/geom_matching/wasserstein/include/dnn/local/kd-tree.hpp new file mode 100644 index 0000000..151a4ad --- /dev/null +++ b/geom_matching/wasserstein/include/dnn/local/kd-tree.hpp @@ -0,0 +1,303 @@ +#include <boost/range/counting_range.hpp> +#include <boost/range/algorithm_ext/push_back.hpp> +#include <boost/range.hpp> + +#include <queue> +#include <stack> + +#include "../parallel/tbb.h" + +template<class T> +dnn::KDTree<T>:: +KDTree(const Traits& traits, HandleContainer&& handles, double _wassersteinPower): + traits_(traits), tree_(std::move(handles)), wassersteinPower(_wassersteinPower) +{ assert(wassersteinPower >= 1.0); init(); } + +template<class T> +template<class Range> +dnn::KDTree<T>:: +KDTree(const Traits& traits, const Range& range, double _wassersteinPower): + traits_(traits), wassersteinPower(_wassersteinPower) +{ + assert( wassersteinPower >= 1.0); + init(range); +} + +template<class T> +template<class Range> +void +dnn::KDTree<T>:: +init(const Range& range) +{ + size_t sz = std::distance(std::begin(range), std::end(range)); + tree_.reserve(sz); + weights_.resize(sz, 0); + subtree_weights_.resize(sz, 0); + for (PointHandle h : range) + tree_.push_back(h); + init(); +} + +template<class T> +void +dnn::KDTree<T>:: +init() +{ + if (tree_.empty()) + return; + +#if defined(TBB) + task_group g; + g.run(OrderTree(tree_.begin(), tree_.end(), 0, traits())); + g.wait(); +#else + OrderTree(tree_.begin(), tree_.end(), 0, traits()).serial(); +#endif + + for (size_t i = 0; i < tree_.size(); ++i) + indices_[tree_[i]] = i; +} + +template<class T> +struct +dnn::KDTree<T>::OrderTree +{ + OrderTree(HCIterator b_, HCIterator e_, size_t i_, const Traits& traits_): + b(b_), e(e_), i(i_), traits(traits_) {} + + void operator()() const + { + if (e - b < 1000) + { + serial(); + return; + } + + HCIterator m = b + (e - b)/2; + CoordinateComparison cmp(i, traits); + std::nth_element(b,m,e, cmp); + size_t next_i = (i + 1) % traits.dimension(); + + task_group g; + if (b < m - 1) g.run(OrderTree(b, m, next_i, traits)); + if (e > m + 2) g.run(OrderTree(m+1, e, next_i, traits)); + g.wait(); + } + + void serial() const + { + std::queue<KDTreeNode> q; + q.push(KDTreeNode(b,e,i)); + while (!q.empty()) + { + HCIterator b, e; size_t i; + std::tie(b,e,i) = q.front(); + q.pop(); + HCIterator m = b + (e - b)/2; + + CoordinateComparison cmp(i, traits); + std::nth_element(b,m,e, cmp); + size_t next_i = (i + 1) % traits.dimension(); + + // Replace with a size condition instead? + if (b < m - 1) q.push(KDTreeNode(b, m, next_i)); + if (e - m > 2) q.push(KDTreeNode(m+1, e, next_i)); + } + } + + HCIterator b, e; + size_t i; + const Traits& traits; +}; + +template<class T> +template<class ResultsFunctor> +void +dnn::KDTree<T>:: +search(PointHandle q, ResultsFunctor& rf) const +{ + typedef typename HandleContainer::const_iterator HCIterator; + typedef std::tuple<HCIterator, HCIterator, size_t> KDTreeNode; + + if (tree_.empty()) + return; + + DistanceType D = std::numeric_limits<DistanceType>::infinity(); + + // TODO: use tbb::scalable_allocator for the queue + std::queue<KDTreeNode> nodes; + + + + nodes.push(KDTreeNode(tree_.begin(), tree_.end(), 0)); + + + //std::cout << "started kdtree::search" << std::endl; + + while (!nodes.empty()) + { + HCIterator b, e; size_t i; + std::tie(b,e,i) = nodes.front(); + nodes.pop(); + + CoordinateComparison cmp(i, traits()); + i = (i + 1) % traits().dimension(); + + HCIterator m = b + (e - b)/2; + DistanceType dist = pow(traits().distance(q, *m), wassersteinPower) + weights_[m - tree_.begin()]; + + + D = rf(*m, dist); + + // we are really searching w.r.t L_\infty ball; could prune better with an L_2 ball + Coordinate diff = cmp.diff(q, *m); // diff returns signed distance + DistanceType diffToWasserPower = (diff > 0 ? 1.0 : -1.0) * pow(fabs(diff), wassersteinPower); + + size_t lm = m + 1 + (e - (m+1))/2 - tree_.begin(); + if (e > m + 1 && diffToWasserPower - subtree_weights_[lm] >= -D) { + nodes.push(KDTreeNode(m+1, e, i)); + } + + size_t rm = b + (m - b) / 2 - tree_.begin(); + if (b < m && diffToWasserPower + subtree_weights_[rm] <= D) { + nodes.push(KDTreeNode(b, m, i)); + } + } + //std::cout << "exited kdtree::search" << std::endl; +} + +template<class T> +void +dnn::KDTree<T>:: +increase_weight(PointHandle p, DistanceType w) +{ + size_t idx = indices_[p]; + // weight should only increase + assert( weights_[idx] <= w ); + weights_[idx] = w; + + typedef std::tuple<HCIterator, HCIterator> KDTreeNode; + + // find the path down the tree to this node + // not an ideal strategy, but // it's not clear how to move up from the node in general + std::stack<KDTreeNode> s; + s.push(KDTreeNode(tree_.begin(),tree_.end())); + + do + { + HCIterator b,e; + std::tie(b,e) = s.top(); + + size_t im = b + (e - b)/2 - tree_.begin(); + + if (idx == im) + break; + else if (idx < im) + s.push(KDTreeNode(b, tree_.begin() + im)); + else // idx > im + s.push(KDTreeNode(tree_.begin() + im + 1, e)); + } while(1); + + // update subtree_weights_ on the path to the root + DistanceType min_w = w; + while (!s.empty()) + { + HCIterator b,e; + std::tie(b,e) = s.top(); + HCIterator m = b + (e - b)/2; + size_t im = m - tree_.begin(); + s.pop(); + + + // left and right children + if (b < m) + { + size_t lm = b + (m - b)/2 - tree_.begin(); + if (subtree_weights_[lm] < min_w) + min_w = subtree_weights_[lm]; + } + + if (e > m + 1) + { + size_t rm = m + 1 + (e - (m+1))/2 - tree_.begin(); + if (subtree_weights_[rm] < min_w) + min_w = subtree_weights_[rm]; + } + + if (weights_[im] < min_w) { + min_w = weights_[im]; + } + + if (subtree_weights_[im] < min_w ) // increase weight + subtree_weights_[im] = min_w; + else + break; + } +} + +template<class T> +typename dnn::KDTree<T>::HandleDistance +dnn::KDTree<T>:: +find(PointHandle q) const +{ + dnn::NNRecord<HandleDistance> nn; + search(q, nn); + return nn.result; +} + +template<class T> +typename dnn::KDTree<T>::Result +dnn::KDTree<T>:: +findR(PointHandle q, DistanceType r) const +{ + dnn::rNNRecord<HandleDistance> rnn(r); + search(q, rnn); + std::sort(rnn.result.begin(), rnn.result.end()); + return rnn.result; +} + +template<class T> +typename dnn::KDTree<T>::Result +dnn::KDTree<T>:: +findK(PointHandle q, size_t k) const +{ + dnn::kNNRecord<HandleDistance> knn(k); + search(q, knn); + std::sort(knn.result.begin(), knn.result.end()); + return knn.result; +} + + +template<class T> +struct dnn::KDTree<T>::CoordinateComparison +{ + CoordinateComparison(size_t i, const Traits& traits): + i_(i), traits_(traits) {} + + bool operator()(PointHandle p1, PointHandle p2) const { return coordinate(p1) < coordinate(p2); } + Coordinate diff(PointHandle p1, PointHandle p2) const { return coordinate(p1) - coordinate(p2); } + + Coordinate coordinate(PointHandle p) const { return traits_.coordinate(p, i_); } + size_t axis() const { return i_; } + + private: + size_t i_; + const Traits& traits_; +}; + +template<class T> +void +dnn::KDTree<T>:: +printWeights(void) +{ + std::cout << "weights_:" << std::endl; + for(const auto ph : indices_) { + std::cout << "idx = " << ph.second << ": (" << (ph.first)->at(0) << ", " << (ph.first)->at(1) << ") weight = " << weights_[ph.second] << std::endl; + } + std::cout << "subtree_weights_:" << std::endl; + for(size_t idx = 0; idx < subtree_weights_.size(); ++idx) { + std::cout << idx << " : " << subtree_weights_[idx] << std::endl; + } +} + + diff --git a/geom_matching/wasserstein/include/dnn/local/search-functors.h b/geom_matching/wasserstein/include/dnn/local/search-functors.h new file mode 100644 index 0000000..f257d0c --- /dev/null +++ b/geom_matching/wasserstein/include/dnn/local/search-functors.h @@ -0,0 +1,89 @@ +#ifndef DNN_LOCAL_SEARCH_FUNCTORS_H +#define DNN_LOCAL_SEARCH_FUNCTORS_H + +#include <boost/range/algorithm/heap_algorithm.hpp> + +namespace dnn +{ + +template<class NN> +struct HandleDistance +{ + typedef typename NN::PointHandle PointHandle; + typedef typename NN::DistanceType DistanceType; + typedef typename NN::HDContainer HDContainer; + + HandleDistance() {} + HandleDistance(PointHandle pp, DistanceType dd): + p(pp), d(dd) {} + bool operator<(const HandleDistance& other) const { return d < other.d; } + + PointHandle p; + DistanceType d; +}; + +template<class HandleDistance> +struct NNRecord +{ + typedef typename HandleDistance::PointHandle PointHandle; + typedef typename HandleDistance::DistanceType DistanceType; + + NNRecord() { result.d = std::numeric_limits<DistanceType>::infinity(); } + DistanceType operator()(PointHandle p, DistanceType d) { if (d < result.d) { result.p = p; result.d = d; } return result.d; } + HandleDistance result; +}; + +template<class HandleDistance> +struct rNNRecord +{ + typedef typename HandleDistance::PointHandle PointHandle; + typedef typename HandleDistance::DistanceType DistanceType; + typedef typename HandleDistance::HDContainer HDContainer; + + rNNRecord(DistanceType r_): r(r_) {} + DistanceType operator()(PointHandle p, DistanceType d) + { + if (d <= r) + result.push_back(HandleDistance(p,d)); + return r; + } + + DistanceType r; + HDContainer result; +}; + +template<class HandleDistance> +struct kNNRecord +{ + typedef typename HandleDistance::PointHandle PointHandle; + typedef typename HandleDistance::DistanceType DistanceType; + typedef typename HandleDistance::HDContainer HDContainer; + + kNNRecord(unsigned k_): k(k_) {} + DistanceType operator()(PointHandle p, DistanceType d) + { + if (result.size() < k) + { + result.push_back(HandleDistance(p,d)); + boost::push_heap(result); + if (result.size() < k) + return std::numeric_limits<DistanceType>::infinity(); + } else if (d < result[0].d) + { + boost::pop_heap(result); + result.back() = HandleDistance(p,d); + boost::push_heap(result); + } + if ( result.size() > 1 ) { + assert( result[0].d >= result[1].d ); + } + return result[0].d; + } + + unsigned k; + HDContainer result; +}; + +} + +#endif // DNN_LOCAL_SEARCH_FUNCTORS_H diff --git a/geom_matching/wasserstein/include/dnn/parallel/tbb.h b/geom_matching/wasserstein/include/dnn/parallel/tbb.h new file mode 100644 index 0000000..4aa6805 --- /dev/null +++ b/geom_matching/wasserstein/include/dnn/parallel/tbb.h @@ -0,0 +1,220 @@ +#ifndef PARALLEL_H +#define PARALLEL_H + +#include <iostream> +#include <vector> + +#include <boost/range.hpp> +#include <boost/bind.hpp> +#include <boost/foreach.hpp> + +#ifdef TBB + +#include <tbb/tbb.h> +#include <tbb/concurrent_hash_map.h> +#include <tbb/scalable_allocator.h> + +#include <boost/serialization/split_free.hpp> +#include <boost/serialization/collections_load_imp.hpp> +#include <boost/serialization/collections_save_imp.hpp> + +namespace dnn +{ + using tbb::mutex; + using tbb::task_scheduler_init; + using tbb::task_group; + using tbb::task; + + template<class T> + struct vector + { + typedef tbb::concurrent_vector<T> type; + }; + + template<class T> + struct atomic + { + typedef tbb::atomic<T> type; + static T compare_and_swap(type& v, T n, T o) { return v.compare_and_swap(n,o); } + }; + + template<class Iterator, class F> + void do_foreach(Iterator begin, Iterator end, const F& f) { tbb::parallel_do(begin, end, f); } + + template<class Range, class F> + void for_each_range_(const Range& r, const F& f) + { + for (typename Range::iterator cur = r.begin(); cur != r.end(); ++cur) + f(*cur); + } + + template<class F> + void for_each_range(size_t from, size_t to, const F& f) + { + //static tbb::affinity_partitioner ap; + //tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::range_type, F>, _1, f), ap); + tbb::parallel_for(from, to, f); + } + + template<class Container, class F> + void for_each_range(const Container& c, const F& f) + { + //static tbb::affinity_partitioner ap; + //tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::range_type, F>, _1, f), ap); + tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::const_range_type, F>, _1, f)); + } + + template<class Container, class F> + void for_each_range(Container& c, const F& f) + { + //static tbb::affinity_partitioner ap; + //tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::range_type, F>, _1, f), ap); + tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::range_type, F>, _1, f)); + } + + template<class ID, class NodePointer, class IDTraits, class Allocator> + struct map_traits + { + typedef tbb::concurrent_hash_map<ID, NodePointer, IDTraits, Allocator> type; + typedef typename type::range_type range; + }; + + struct progress_timer + { + progress_timer(): start(tbb::tick_count::now()) {} + ~progress_timer() + { std::cout << (tbb::tick_count::now() - start).seconds() << " s" << std::endl; } + + tbb::tick_count start; + }; +} + +// Serialization for tbb::concurrent_vector<...> +namespace boost +{ + namespace serialization + { + template<class Archive, class T, class A> + void save(Archive& ar, const tbb::concurrent_vector<T,A>& v, const unsigned int file_version) + { stl::save_collection(ar, v); } + + template<class Archive, class T, class A> + void load(Archive& ar, tbb::concurrent_vector<T,A>& v, const unsigned int file_version) + { + stl::load_collection<Archive, + tbb::concurrent_vector<T,A>, + stl::archive_input_seq< Archive, tbb::concurrent_vector<T,A> >, + stl::reserve_imp< tbb::concurrent_vector<T,A> > + >(ar, v); + } + + template<class Archive, class T, class A> + void serialize(Archive& ar, tbb::concurrent_vector<T,A>& v, const unsigned int file_version) + { split_free(ar, v, file_version); } + + template<class Archive, class T> + void save(Archive& ar, const tbb::atomic<T>& v, const unsigned int file_version) + { T v_ = v; ar << v_; } + + template<class Archive, class T> + void load(Archive& ar, tbb::atomic<T>& v, const unsigned int file_version) + { T v_; ar >> v_; v = v_; } + + template<class Archive, class T> + void serialize(Archive& ar, tbb::atomic<T>& v, const unsigned int file_version) + { split_free(ar, v, file_version); } + } +} + +#else + +#include <algorithm> +#include <map> +#include <boost/progress.hpp> + +namespace dnn +{ + template<class T> + struct vector + { + typedef ::std::vector<T> type; + }; + + template<class T> + struct atomic + { + typedef T type; + static T compare_and_swap(type& v, T n, T o) { if (v != o) return v; v = n; return o; } + }; + + template<class Iterator, class F> + void do_foreach(Iterator begin, Iterator end, const F& f) { std::for_each(begin, end, f); } + + template<class F> + void for_each_range(size_t from, size_t to, const F& f) + { + for (size_t i = from; i < to; ++i) + f(i); + } + + template<class Container, class F> + void for_each_range(Container& c, const F& f) + { + BOOST_FOREACH(const typename Container::value_type& i, c) + f(i); + } + + template<class Container, class F> + void for_each_range(const Container& c, const F& f) + { + BOOST_FOREACH(const typename Container::value_type& i, c) + f(i); + } + + struct mutex + { + struct scoped_lock + { + scoped_lock() {} + scoped_lock(mutex& ) {} + void acquire(mutex& ) const {} + void release() const {} + }; + }; + + struct task_scheduler_init + { + task_scheduler_init(unsigned) {} + void initialize(unsigned) {} + static const unsigned automatic = 0; + static const unsigned deferred = 0; + }; + + struct task_group + { + template<class Functor> + void run(const Functor& f) const { f(); } + void wait() const {} + }; + + template<class ID, class NodePointer, class IDTraits, class Allocator> + struct map_traits + { + typedef std::map<ID, NodePointer, + typename IDTraits::Comparison, + Allocator> type; + typedef type range; + }; + + using boost::progress_timer; +} + +#endif // TBB + +namespace dnn +{ + template<class Range, class F> + void do_foreach(const Range& range, const F& f) { do_foreach(boost::begin(range), boost::end(range), f); } +} + +#endif diff --git a/geom_matching/wasserstein/include/dnn/parallel/utils.h b/geom_matching/wasserstein/include/dnn/parallel/utils.h new file mode 100644 index 0000000..ba73814 --- /dev/null +++ b/geom_matching/wasserstein/include/dnn/parallel/utils.h @@ -0,0 +1,94 @@ +#ifndef PARALLEL_UTILS_H +#define PARALLEL_UTILS_H + +#include "../utils.h" + +namespace dnn +{ + // Assumes rng is synchronized across ranks + template<class DataVector, class RNGType, class SwapFunctor> + void shuffle(mpi::communicator& world, DataVector& data, RNGType& rng, const SwapFunctor& swap, DataVector empty = DataVector()); + + template<class DataVector, class RNGType> + void shuffle(mpi::communicator& world, DataVector& data, RNGType& rng) + { + typedef decltype(data[0]) T; + shuffle(world, data, rng, [](T& x, T& y) { std::swap(x,y); }); + } +} + +template<class DataVector, class RNGType, class SwapFunctor> +void +dnn::shuffle(mpi::communicator& world, DataVector& data, RNGType& rng, const SwapFunctor& swap, DataVector empty) +{ + // This is not a perfect shuffle: it dishes out data in chunks of 1/size. + // (It can be interpreted as generating a bistochastic matrix by taking the + // sum of size random permutation matrices.) Hopefully, it works for our purposes. + + typedef typename RNGType::result_type RNGResult; + + int size = world.size(); + int rank = world.rank(); + + // Generate local seeds + boost::uniform_int<RNGResult> uniform; + RNGResult seed; + for (size_t i = 0; i < size; ++i) + { + RNGResult v = uniform(rng); + if (i == rank) + seed = v; + } + RNGType local_rng(seed); + + // Shuffle local data + dnn::random_shuffle(data.begin(), data.end(), local_rng, swap); + + // Decide how much of our data goes to i-th processor + std::vector<size_t> out_counts(size); + std::vector<int> ranks(boost::counting_iterator<int>(0), + boost::counting_iterator<int>(size)); + for (size_t i = 0; i < size; ++i) + { + dnn::random_shuffle(ranks.begin(), ranks.end(), rng); + ++out_counts[ranks[rank]]; + } + + // Fill the outgoing array + size_t total = 0; + std::vector< DataVector > outgoing(size, empty); + for (size_t i = 0; i < size; ++i) + { + size_t count = data.size()*out_counts[i]/size; + if (total + count > data.size()) + count = data.size() - total; + + outgoing[i].reserve(count); + for (size_t j = total; j < total + count; ++j) + outgoing[i].push_back(data[j]); + + total += count; + } + + boost::uniform_int<size_t> uniform_outgoing(0,size-1); // in range [0,size-1] + while(total < data.size()) // send leftover to random processes + { + outgoing[uniform_outgoing(local_rng)].push_back(data[total]); + ++total; + } + data.clear(); + + // Exchange the data + std::vector< DataVector > incoming(size, empty); + mpi::all_to_all(world, outgoing, incoming); + outgoing.clear(); + + // Assemble our data + for(const DataVector& vec : incoming) + for (size_t i = 0; i < vec.size(); ++i) + data.push_back(vec[i]); + dnn::random_shuffle(data.begin(), data.end(), local_rng, swap); + // XXX: the final shuffle is irrelevant for our purposes. But it's also cheap. +} + +#endif diff --git a/geom_matching/wasserstein/include/dnn/utils.h b/geom_matching/wasserstein/include/dnn/utils.h new file mode 100644 index 0000000..83c2865 --- /dev/null +++ b/geom_matching/wasserstein/include/dnn/utils.h @@ -0,0 +1,41 @@ +#ifndef DNN_UTILS_H +#define DNN_UTILS_H + +#include <boost/random/uniform_int.hpp> +#include <boost/foreach.hpp> +#include <boost/typeof/typeof.hpp> + +namespace dnn +{ + +template <typename T, typename... Args> +struct has_coordinates +{ + template <typename C, typename = decltype( std::declval<C>().coordinate(std::declval<Args>()...) )> + static std::true_type test(int); + + template <typename C> + static std::false_type test(...); + + static constexpr bool value = decltype(test<T>(0))::value; +}; + +template<class RandomIt, class UniformRandomNumberGenerator, class SwapFunctor> +void random_shuffle(RandomIt first, RandomIt last, UniformRandomNumberGenerator& g, const SwapFunctor& swap) +{ + size_t n = last - first; + boost::uniform_int<size_t> uniform(0,n); + for (size_t i = n-1; i > 0; --i) + swap(first[i], first[uniform(g,i+1)]); // picks a random number in [0,i] range +} + +template<class RandomIt, class UniformRandomNumberGenerator> +void random_shuffle(RandomIt first, RandomIt last, UniformRandomNumberGenerator& g) +{ + typedef decltype(*first) T; + random_shuffle(first, last, g, [](T& x, T& y) { std::swap(x,y); }); +} + +} + +#endif diff --git a/geom_matching/wasserstein/include/wasserstein.h b/geom_matching/wasserstein/include/wasserstein.h new file mode 100644 index 0000000..38ac6bd --- /dev/null +++ b/geom_matching/wasserstein/include/wasserstein.h @@ -0,0 +1,110 @@ +/* + +Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#ifndef WASSERSTEIN_H +#define WASSERSTEIN_H + +#include <vector> +#include <map> +#include <math.h> + +#include "basic_defs_ws.h" + +// use Gauss-Seidel version; comment out to switch to Jacobi (not recommended) +#define GAUSS_SEIDEL_AUCTION + +namespace geom_ws { + +// get Wasserstein distance between two persistence diagrams +double wassersteinDistVec(const std::vector<DiagramPoint>& A, + const std::vector<DiagramPoint>& B, + const double q, + const double delta, + const double _internal_p = std::numeric_limits<double>::infinity(), + const double _initialEpsilon = 0.0, + const double _epsFactor = 0.0); + + +// compare as multisets +template<class PairContainer> +bool areEqual(PairContainer& dgm1, PairContainer& dgm2) +{ + if (dgm1.size() != dgm2.size()) { + return false; + } + + std::map<std::pair<double, double>, int> m1, m2; + + for(const auto& pair1 : dgm1) { + m1[pair1]++; + } + + for(const auto& pair2 : dgm2) { + m2[pair2]++; + } + + return m1 == m2; +} + +template<class PairContainer> +double wassersteinDist(PairContainer& A, PairContainer& B, const double q, const double delta, const double _internal_p = std::numeric_limits<double>::infinity(), const double _initialEpsilon = 0.0, const double _epsFactor = 0.0) +{ + if (areEqual(A, B)) { + return 0.0; + } + + std::vector<DiagramPoint> dgmA, dgmB; + // loop over A, add projections of A-points to corresponding positions + // in B-vector + for(auto& pairA : A) { + double x = pairA.first; + double y = pairA.second; + dgmA.push_back(DiagramPoint(x, y, DiagramPoint::NORMAL)); + dgmB.push_back(DiagramPoint(x, y, DiagramPoint::DIAG)); + } + // the same for B + for(auto& pairB : B) { + double x = pairB.first; + double y = pairB.second; + dgmA.push_back(DiagramPoint(x, y, DiagramPoint::DIAG)); + dgmB.push_back(DiagramPoint(x, y, DiagramPoint::NORMAL)); + } + + return wassersteinDistVec(dgmA, dgmB, q, delta, _internal_p, _initialEpsilon, _epsFactor); +} + + +// fill in result with points from file fname +// return false if file can't be opened +// or error occurred while reading +bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double>>& result); +bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<double, double>>& result); + +} // end of namespace geom_ws + +#endif diff --git a/geom_matching/wasserstein/src/auction_oracle.cpp b/geom_matching/wasserstein/src/auction_oracle.cpp new file mode 100644 index 0000000..f906d26 --- /dev/null +++ b/geom_matching/wasserstein/src/auction_oracle.cpp @@ -0,0 +1,1310 @@ +/* + +Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#include <assert.h> +#include <algorithm> +#include <functional> +#include <iterator> + +#include "def_debug.h" +#include "auction_oracle.h" + +namespace geom_ws { + +AuctionOracleAbstract::AuctionOracleAbstract(const std::vector<DiagramPoint>& _bidders, const std::vector<DiagramPoint>& _items, const double _wassersteinPower, const double _internal_p) : + bidders(_bidders), + items(_items), + prices(items.size(), 0.0), + wassersteinPower(_wassersteinPower), + internal_p(_internal_p) +{ +} + +double AuctionOracleAbstract::getValueForBidder(size_t bidderIdx, size_t itemIdx) +{ + return pow(distLp(bidders[bidderIdx], items[itemIdx], internal_p), wassersteinPower) + prices[itemIdx]; +} + +// ***************************** +// AuctionOracleLazyHeap +// ***************************** + +AuctionOracleLazyHeap::AuctionOracleLazyHeap(const std::vector<DiagramPoint>& b, + const std::vector<DiagramPoint>& g, + double _wassersteinPower, + const double internal_p) : + AuctionOracleAbstract(b, g, _wassersteinPower, internal_p), + maxVal(std::numeric_limits<double>::min()), + biddersUpdateMoments(b.size(), 0), + updateCounter(0) +{ + assert(b.size() == g.size() ); + assert(b.size() > 1); + + weightMatrix.reserve(b.size()); + //const double maxDistUpperBound = 3 * getFurthestDistance3Approx(b, g); + //weightAdjConst = pow(maxDistUpperBound, wassersteinPower); + //std::cout << "3getFurthestDistance3Approx = " << getFurthestDistance3Approx(b, g) << std::endl; + //std::cout << "in AuctionOracleLazyHeap weightAdjConst = " << weightAdjConst << std::endl; + // init weight matrix + for(const auto& pointA : bidders) { + std::vector<double> weightVec; + weightVec.clear(); + weightVec.reserve(b.size()); + for(const auto& pointB : items) { + double val = pow(distLp(pointA, pointB, internal_p), wassersteinPower); + if (val > maxVal) { + maxVal = val; + } + weightVec.push_back( val ); + } + weightMatrix.push_back(weightVec); + } + fillInLossesHeap(); + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + updateList.push_back(std::make_pair(static_cast<IdxType>(itemIdx), 0)); + } + for(auto updateListIter = updateList.begin(); updateListIter != updateList.end(); ++updateListIter) { + itemsIterators.push_back(updateListIter); + } +} + +void AuctionOracleLazyHeap::updateQueueForBidder(IdxType bidderIdx) +{ + assert(0 <= bidderIdx and bidderIdx < static_cast<int>(bidders.size())); + assert(bidderIdx < static_cast<int>(biddersUpdateMoments.size())); + + int bidderLastUpdateTime = biddersUpdateMoments[bidderIdx]; + auto iter = updateList.begin(); + while (iter != updateList.end() and iter->second >= bidderLastUpdateTime) { + IdxType itemIdx = iter->first; + IdxValPair newVal { itemIdx, weightMatrix[bidderIdx][itemIdx] + prices[itemIdx]}; + // to-do: change indexing of lossesHeapHandles + lossesHeap[bidderIdx]->decrease(lossesHeapHandles[bidderIdx][itemIdx], newVal); + iter++; + } + biddersUpdateMoments[bidderIdx] = updateCounter; +} + +void AuctionOracleLazyHeap::fillInLossesHeap(void) +{ + for(size_t bidderIdx = 0; bidderIdx < bidders.size(); ++bidderIdx) { + lossesHeap.push_back( new LossesHeap() ); + std::vector<LossesHeap::handle_type> handlesVec; + lossesHeapHandles.push_back(handlesVec); + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + IdxValPair vp { itemIdx, weightMatrix[bidderIdx][itemIdx] + prices[itemIdx] }; + lossesHeapHandles[bidderIdx].push_back( lossesHeap[bidderIdx]->push(vp) ); + } + } +} + +AuctionOracleLazyHeap::~AuctionOracleLazyHeap() +{ + for(auto h : lossesHeap) { + delete h; + } +} + +void AuctionOracleLazyHeap::setPrice(IdxType itemIdx, double newPrice) +{ + assert( prices.at(itemIdx) < newPrice ); +#ifdef DEBUG_AUCTION + std::cout << "price incremented by " << prices.at(itemIdx) - newPrice << std::endl; +#endif + prices[itemIdx] = newPrice; + // lazy: record the moment we updated the price of the items, + // do not update queues. + // 1. move the items with updated price to the front of the updateList, + updateList.splice(updateList.begin(), updateList, itemsIterators[itemIdx]); + // 2. record the moment we updated the price and increase the counter + updateList.front().second = updateCounter++; +} + +// subtract min. price from all prices +void AuctionOracleLazyHeap::adjustPrices(void) +{ + double minPrice = *(std::min_element(prices.begin(), prices.end())); + std::transform(prices.begin(), prices.end(), prices.begin(), [minPrice](double a) { return a - minPrice; }); +} + +DebugOptimalBid AuctionOracleLazyHeap::getOptimalBidDebug(IdxType bidderIdx) +{ + assert(bidderIdx >=0 and bidderIdx < static_cast<IdxType>(bidders.size()) ); + assert(lossesHeap.at(bidderIdx) != nullptr); + assert(lossesHeap[bidderIdx]->size() >= 2); + + updateQueueForBidder(bidderIdx); + DebugOptimalBid result; + + auto pHeap = lossesHeap[bidderIdx]; + auto topIter = pHeap->ordered_begin(); + result.bestItemIdx = topIter->first; + result.bestItemValue = topIter->second; + ++topIter; // now points to the second-best items + result.secondBestItemValue = topIter->second; + result.secondBestItemIdx = topIter->first; + + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemsDist= " << (weightAdjConst - bestItemValue) << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestDist= " << (weightAdjConst - secondBestItemValue) << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + + return result; +} + +IdxValPair AuctionOracleLazyHeap::getOptimalBid(const IdxType bidderIdx) +{ + assert(bidderIdx >=0 and bidderIdx < static_cast<IdxType>(bidders.size()) ); + assert(lossesHeap.at(bidderIdx) != nullptr); + assert(lossesHeap[bidderIdx]->size() >= 2); + + updateQueueForBidder(bidderIdx); + + auto pHeap = lossesHeap[bidderIdx]; + auto topIter = pHeap->ordered_begin(); + IdxType bestItemIdx = topIter->first; + double bestItemValue = topIter->second; + ++topIter; // now points to the second-best items + double secondBestItemValue = topIter->second; + + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemsDist= " << (weightAdjConst - bestItemValue) << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestDist= " << (weightAdjConst - secondBestItemValue) << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + + // bid value: price + value difference + epsilon + return std::make_pair(bestItemIdx, + prices[bestItemIdx] + + ( secondBestItemValue - bestItemValue ) + + epsilon ); +} + +// ***************************** +// AuctionOracleLazyHeapRestricted +// ***************************** + +AuctionOracleLazyHeapRestricted::AuctionOracleLazyHeapRestricted(const std::vector<DiagramPoint>& b, + const std::vector<DiagramPoint>& g, + double _wassersteinPower, + double internal_p) : + AuctionOracleAbstract(b, g, _wassersteinPower), + maxVal(std::numeric_limits<double>::min()), + biddersUpdateMoments(b.size(), 0), + updateCounter(0), + heapHandlesIndices(items.size(), std::numeric_limits<size_t>::max()), + bestDiagonalItemsComputed(false) +{ + assert(b.size() == g.size() ); + assert(b.size() > 1); + + weightMatrix.reserve(b.size()); + //const double maxDistUpperBound = 3 * getFurthestDistance3Approx(b, g); + //weightAdjConst = pow(maxDistUpperBound, wassersteinPower); + //std::cout << "3getFurthestDistance3Approx = " << getFurthestDistance3Approx(b, g) << std::endl; + //std::cout << "in AuctionOracleLazyHeapRestricted weightAdjConst = " << weightAdjConst << std::endl; + // init weight matrix + for(const auto& pointA : bidders) { + std::vector<double> weightVec; + weightVec.clear(); + weightVec.reserve(b.size()); + for(const auto& pointB : items) { + double val = pow(distLp(pointA, pointB, internal_p), wassersteinPower); + weightVec.push_back( val ); + if ( val > maxVal ) + maxVal = val; + } + weightMatrix.push_back(weightVec); + } + fillInLossesHeap(); + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + updateList.push_back(std::make_pair(static_cast<IdxType>(itemIdx), 0)); + } + for(auto updateListIter = updateList.begin(); updateListIter != updateList.end(); ++updateListIter) { + itemsIterators.push_back(updateListIter); + } + + size_t handleIdx {0}; + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + if (items[itemIdx].isDiagonal() ) { + heapHandlesIndices[itemIdx] = handleIdx++; + diagHeapHandles.push_back(diagItemsHeap.push(std::make_pair(itemIdx, 0))); + } + } + // todo: this must be done in readFiles procedure +} + +void AuctionOracleLazyHeapRestricted::updateQueueForBidder(IdxType bidderIdx) +{ + assert(0 <= bidderIdx and bidderIdx < static_cast<int>(bidders.size())); + assert(bidderIdx < static_cast<int>(biddersUpdateMoments.size())); + assert(lossesHeap[bidderIdx] != nullptr ); + + int bidderLastUpdateTime = biddersUpdateMoments[bidderIdx]; + auto iter = updateList.begin(); + while (iter != updateList.end() and iter->second >= bidderLastUpdateTime) { + IdxType itemIdx = iter->first; + size_t handleIdx = itemsIndicesForHeapHandles[bidderIdx][itemIdx]; + if (handleIdx < items.size() ) { + IdxValPair newVal { itemIdx, weightMatrix[bidderIdx][itemIdx] + prices[itemIdx]}; + // to-do: change indexing of lossesHeapHandles + lossesHeap[bidderIdx]->decrease(lossesHeapHandles[bidderIdx][handleIdx], newVal); + } + iter++; + } + biddersUpdateMoments[bidderIdx] = updateCounter; +} + +void AuctionOracleLazyHeapRestricted::fillInLossesHeap(void) +{ + for(size_t bidderIdx = 0; bidderIdx < bidders.size(); ++bidderIdx) { + DiagramPoint bidder { bidders[bidderIdx] }; + // no heap for diagonal bidders + if ( bidder.isDiagonal() ) { + lossesHeap.push_back( nullptr ); + lossesHeapHandles.push_back(std::vector<LossesHeap::handle_type>()); + itemsIndicesForHeapHandles.push_back( std::vector<size_t>() ); + continue; + } else { + lossesHeap.push_back( new LossesHeap() ); + assert( lossesHeap.at(bidderIdx) != nullptr ); + itemsIndicesForHeapHandles.push_back( std::vector<size_t>(items.size(), std::numeric_limits<size_t>::max() ) ); + + std::vector<LossesHeap::handle_type> handlesVec; + lossesHeapHandles.push_back(handlesVec); + size_t handleIdx { 0 }; + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + assert( itemsIndicesForHeapHandles.at(bidderIdx).at(itemIdx) > 0 ); + DiagramPoint item { items[itemIdx] }; + if ( item.isNormal() ) { + // item can be assigned to bidder, store in heap + IdxValPair vp { itemIdx, weightMatrix[bidderIdx][itemIdx] + prices[itemIdx] }; + lossesHeapHandles[bidderIdx].push_back( lossesHeap[bidderIdx]->push(vp) ); + // keep corresponding index in itemsIndicesForHeapHandles + itemsIndicesForHeapHandles[bidderIdx][itemIdx] = handleIdx++; + } + } + } + } +} + +AuctionOracleLazyHeapRestricted::~AuctionOracleLazyHeapRestricted() +{ + for(auto h : lossesHeap) { + delete h; + } +} + +void AuctionOracleLazyHeapRestricted::setPrice(IdxType itemIdx, double newPrice) +{ + assert( prices.at(itemIdx) < newPrice ); +#ifdef DEBUG_AUCTION + std::cout << "price incremented by " << prices.at(itemIdx) - newPrice << std::endl; +#endif + prices[itemIdx] = newPrice; + if (items[itemIdx].isNormal() ) { + // lazy: record the moment we updated the price of the items, + // do not update queues. + // 1. move the items with updated price to the front of the updateList, + updateList.splice(updateList.begin(), updateList, itemsIterators[itemIdx]); + // 2. record the moment we updated the price and increase the counter + updateList.front().second = updateCounter++; + } else { + // diagonal items are stored in one heap + diagItemsHeap.decrease(diagHeapHandles[heapHandlesIndices[itemIdx]], std::make_pair(itemIdx, newPrice)); + bestDiagonalItemsComputed = false; + } +} + +// subtract min. price from all prices +void AuctionOracleLazyHeapRestricted::adjustPrices(void) +{ +} + +DebugOptimalBid AuctionOracleLazyHeapRestricted::getOptimalBidDebug(IdxType bidderIdx) +{ + DebugOptimalBid result; + assert(bidderIdx >=0 and bidderIdx < static_cast<IdxType>(bidders.size()) ); + + DiagramPoint bidder = bidders[bidderIdx]; + std::vector<IdxValPair> candItems; + // corresponding point is always considered as a candidate + + size_t projItemIdx = bidderIdx; + assert( 0 <= projItemIdx and projItemIdx < items.size() ); + DiagramPoint projItem = items[projItemIdx]; + assert(projItem.type != bidder.type); + //assert(projItem.projId == bidder.id); + //assert(projItem.id == bidder.projId); + // todo: store precomputed distance? + double projItemValue = pow(distLp(bidder, projItem, internal_p), wassersteinPower) + prices[projItemIdx]; + candItems.push_back( std::make_pair(projItemIdx, projItemValue) ); + + if (bidder.isNormal()) { + assert(lossesHeap.at(bidderIdx) != nullptr); + assert(lossesHeap[bidderIdx]->size() >= 2); + updateQueueForBidder(bidderIdx); + auto pHeap = lossesHeap[bidderIdx]; + assert( pHeap != nullptr ); + auto topIter = pHeap->ordered_begin(); + candItems.push_back( *topIter ); + ++topIter; // now points to the second-best items + candItems.push_back( *topIter ); + std::sort(candItems.begin(), candItems.end(), CompPairsBySecondStruct()); + assert(candItems[1].second >= candItems[0].second); + } else { + // for diagonal bidder the only normal point has already been added + // the other 2 candidates are diagonal items only, get from the heap + // with prices + assert(diagItemsHeap.size() > 1); + auto topDiagIter = diagItemsHeap.ordered_begin(); + auto topDiag1 = *topDiagIter++; + auto topDiag2 = *topDiagIter; + candItems.push_back(topDiag1); + candItems.push_back(topDiag2); + std::sort(candItems.begin(), candItems.end(), CompPairsBySecondStruct()); + assert(candItems.size() == 3); + assert(candItems[2].second >= candItems[1].second); + assert(candItems[1].second >= candItems[0].second); + } + + result.bestItemIdx = candItems[0].first; + result.secondBestItemIdx = candItems[1].first; + result.bestItemValue = candItems[0].second; + result.secondBestItemValue = candItems[1].second; + + // checking code + + //DebugOptimalBid debugMyResult(result); + //DebugOptimalBid debugNaiveResult; + //debugNaiveResult.bestItemValue = 1e20; + //debugNaiveResult.secondBestItemValue = 1e20; + //double currItemValue; + //for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + //if ( bidders[bidderIdx].type != items[itemIdx].type and + //bidders[bidderIdx].projId != items[itemIdx].id) + //continue; + + //currItemValue = pow(distLp(bidders[bidderIdx], items[itemIdx]), wassersteinPower) + prices[itemIdx]; + //if (currItemValue < debugNaiveResult.bestItemValue) { + //debugNaiveResult.bestItemValue = currItemValue; + //debugNaiveResult.bestItemIdx = itemIdx; + //} + //} + + //for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + //if (itemIdx == debugNaiveResult.bestItemIdx) { + //continue; + //} + //if ( bidders[bidderIdx].type != items[itemIdx].type and + //bidders[bidderIdx].projId != items[itemIdx].id) + //continue; + + //currItemValue = pow(distLp(bidders[bidderIdx], items[itemIdx]), wassersteinPower) + prices[itemIdx]; + //if (currItemValue < debugNaiveResult.secondBestItemValue) { + //debugNaiveResult.secondBestItemValue = currItemValue; + //debugNaiveResult.secondBestItemIdx = itemIdx; + //} + //} + ////std::cout << "got naive result" << std::endl; + + //if ( fabs( debugMyResult.bestItemValue - debugNaiveResult.bestItemValue ) > 1e-6 or + //fabs( debugNaiveResult.secondBestItemValue - debugMyResult.secondBestItemValue) > 1e-6 ) { + //std::cerr << "bidderIdx = " << bidderIdx << "; "; + //std::cerr << bidders[bidderIdx] << std::endl; + //for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + //std::cout << itemIdx << ": " << items[itemIdx] << "; price = " << prices[itemIdx] << std::endl; + //} + //std::cerr << "debugMyResult: " << debugMyResult << std::endl; + //std::cerr << "debugNaiveResult: " << debugNaiveResult << std::endl; + //auto pHeap = lossesHeap[bidderIdx]; + //assert( pHeap != nullptr ); + //for(auto topIter = pHeap->ordered_begin(); topIter != pHeap->ordered_end(); ++topIter) { + //std::cerr << "in heap: " << topIter->first << ": " << topIter->second << "; real value = " << distLp(bidder, items[topIter->first]) + prices[topIter->first] << std::endl; + //} + //for(auto ci : candItems) { + //std::cout << "ci.idx = " << ci.first << ", value = " << ci.second << std::endl; + //} + + ////std::cerr << "twoBestItems: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; + //assert(false); + //} + + + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemsDist= " << (weightAdjConst - bestItemValue) << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestDist= " << (weightAdjConst - secondBestItemValue) << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + + return result; +} + +IdxValPair AuctionOracleLazyHeapRestricted::getOptimalBid(const IdxType bidderIdx) +{ + IdxType bestItemIdx; + //IdxType secondBestItemIdx; + double bestItemValue; + double secondBestItemValue; + + auto& bidder = bidders[bidderIdx]; + IdxType projItemIdx = bidderIdx; + assert( 0 <= projItemIdx and projItemIdx < items.size() ); + DiagramPoint projItem = items[projItemIdx]; + assert(projItem.type != bidder.type); + //assert(projItem.projId == bidder.id); + //assert(projItem.id == bidder.projId); + // todo: store precomputed distance? + double projItemValue = pow(distLp(bidder, projItem, internal_p), wassersteinPower) + prices[projItemIdx]; + + if (bidder.isDiagonal()) { + // for diagonal bidder the only normal point has already been added + // the other 2 candidates are diagonal items only, get from the heap + // with prices + assert(diagItemsHeap.size() > 1); + if (!bestDiagonalItemsComputed) { + auto topDiagIter = diagItemsHeap.ordered_begin(); + bestDiagonalItemIdx = topDiagIter->first; + bestDiagonalItemValue = topDiagIter->second; + topDiagIter++; + secondBestDiagonalItemIdx = topDiagIter->first; + secondBestDiagonalItemValue = topDiagIter->second; + bestDiagonalItemsComputed = true; + } + + if ( projItemValue < bestDiagonalItemValue) { + bestItemIdx = projItemIdx; + bestItemValue = projItemValue; + secondBestItemValue = bestDiagonalItemValue; + //secondBestItemIdx = bestDiagonalItemIdx; + } else if (projItemValue < secondBestDiagonalItemValue) { + bestItemIdx = bestDiagonalItemIdx; + bestItemValue = bestDiagonalItemValue; + secondBestItemValue = projItemValue; + //secondBestItemIdx = projItemIdx; + } else { + bestItemIdx = bestDiagonalItemIdx; + bestItemValue = bestDiagonalItemValue; + secondBestItemValue = secondBestDiagonalItemValue; + //secondBestItemIdx = secondBestDiagonalItemIdx; + } + } else { + // for normal bidder get 2 best items among non-diagonal (=normal) points + // from the corresponding heap + assert(diagItemsHeap.size() > 1); + updateQueueForBidder(bidderIdx); + auto topNormIter = lossesHeap[bidderIdx]->ordered_begin(); + IdxType bestNormalItemIdx { topNormIter->first }; + double bestNormalItemValue { topNormIter->second }; + topNormIter++; + double secondBestNormalItemValue { topNormIter->second }; + //IdxType secondBestNormalItemIdx { topNormIter->first }; + + if ( projItemValue < bestNormalItemValue) { + bestItemIdx = projItemIdx; + bestItemValue = projItemValue; + secondBestItemValue = bestNormalItemValue; + //secondBestItemIdx = bestNormalItemIdx; + } else if (projItemValue < secondBestNormalItemValue) { + bestItemIdx = bestNormalItemIdx; + bestItemValue = bestNormalItemValue; + secondBestItemValue = projItemValue; + //secondBestItemIdx = projItemIdx; + } else { + bestItemIdx = bestNormalItemIdx; + bestItemValue = bestNormalItemValue; + secondBestItemValue = secondBestNormalItemValue; + //secondBestItemIdx = secondBestNormalItemIdx; + } + } + + IdxValPair result; + + assert( secondBestItemValue >= bestItemValue ); + + result.first = bestItemIdx; + result.second = ( secondBestItemValue - bestItemValue ) + prices[bestItemIdx] + epsilon; + + + // checking code + + //DebugOptimalBid debugMyResult; + //debugMyResult.bestItemIdx = bestItemIdx; + //debugMyResult.bestItemValue = bestItemValue; + //debugMyResult.secondBestItemIdx = secondBestItemIdx; + //debugMyResult.secondBestItemValue = secondBestItemValue; + //DebugOptimalBid debugNaiveResult; + //debugNaiveResult.bestItemValue = 1e20; + //debugNaiveResult.secondBestItemValue = 1e20; + //double currItemValue; + //for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + //if ( bidders[bidderIdx].type != items[itemIdx].type and + //bidders[bidderIdx].projId != items[itemIdx].id) + //continue; + + //currItemValue = pow(distLp(bidders[bidderIdx], items[itemIdx]), wassersteinPower) + prices[itemIdx]; + //if (currItemValue < debugNaiveResult.bestItemValue) { + //debugNaiveResult.bestItemValue = currItemValue; + //debugNaiveResult.bestItemIdx = itemIdx; + //} + //} + + //for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + //if (itemIdx == debugNaiveResult.bestItemIdx) { + //continue; + //} + //if ( bidders[bidderIdx].type != items[itemIdx].type and + //bidders[bidderIdx].projId != items[itemIdx].id) + //continue; + + //currItemValue = pow(distLp(bidders[bidderIdx], items[itemIdx]), wassersteinPower) + prices[itemIdx]; + //if (currItemValue < debugNaiveResult.secondBestItemValue) { + //debugNaiveResult.secondBestItemValue = currItemValue; + //debugNaiveResult.secondBestItemIdx = itemIdx; + //} + //} + ////std::cout << "got naive result" << std::endl; + + //if ( fabs( debugMyResult.bestItemValue - debugNaiveResult.bestItemValue ) > 1e-6 or + //fabs( debugNaiveResult.secondBestItemValue - debugMyResult.secondBestItemValue) > 1e-6 ) { + //std::cerr << "bidderIdx = " << bidderIdx << "; "; + //std::cerr << bidders[bidderIdx] << std::endl; + //for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + //std::cout << itemIdx << ": " << items[itemIdx] << "; price = " << prices[itemIdx] << std::endl; + //} + //std::cerr << "debugMyResult: " << debugMyResult << std::endl; + //std::cerr << "debugNaiveResult: " << debugNaiveResult << std::endl; + //auto pHeap = lossesHeap[bidderIdx]; + //if ( pHeap != nullptr ) { + //for(auto topIter = pHeap->ordered_begin(); topIter != pHeap->ordered_end(); ++topIter) { + //std::cerr << "in heap: " << topIter->first << ": " << topIter->second << "; real value = " << distLp(bidder, items[topIter->first]) + prices[topIter->first] << std::endl; + //} + //} + ////for(auto ci : candItems) { + ////std::cout << "ci.idx = " << ci.first << ", value = " << ci.second << std::endl; + ////} + + ////std::cerr << "twoBestItems: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; + //assert(false); + // } + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemsDist= " << (weightAdjConst - bestItemValue) << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestDist= " << (weightAdjConst - secondBestItemValue) << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + + return result; +} + + +// ***************************** +// AuctionOracleKDTree +// ***************************** + +AuctionOracleKDTree::AuctionOracleKDTree(const std::vector<DiagramPoint>& _bidders, + const std::vector<DiagramPoint>& _items, + double _wassersteinPower, + double internal_p) : + AuctionOracleAbstract(_bidders, _items, _wassersteinPower, internal_p), + heapHandlesIndices(items.size(), std::numeric_limits<size_t>::max()), + kdtreeItems(items.size(), std::numeric_limits<size_t>::max()) +{ + //assert(wassersteinPower == 1.0); // temporarily, to-do: update dnn to search with any q + size_t dnnItemIdx { 0 }; + size_t trueIdx { 0 }; + dnnPoints.clear(); + // store normal items in kd-tree + for(const auto& g : items) { + if (g.isNormal()) { + kdtreeItems[trueIdx] = dnnItemIdx; + // index of items is id of dnn-point + DnnPoint p(trueIdx); + p[0] = g.getRealX(); + p[1] = g.getRealY(); + dnnPoints.push_back(p); + assert(dnnItemIdx == dnnPoints.size() - 1); + dnnItemIdx++; + } + trueIdx++; + } + + assert(dnnPoints.size() < items.size() ); + for(size_t i = 0; i < dnnPoints.size(); ++i) { + dnnPointHandles.push_back(&dnnPoints[i]); + } + DnnTraits traits; + traits.internal_p = internal_p; + //std::cout << "kdtree: " << dnnPointHandles.size() << " points" << std::endl; + kdtree = new dnn::KDTree<DnnTraits>(traits, dnnPointHandles, wassersteinPower); + + size_t dnnItemIdxAll { 0 }; + dnnPointsAll.clear(); + // store all items in kd-tree + for(const auto& g : items) { + DnnPoint p(dnnItemIdxAll++); + p[0] = g.getRealX(); + p[1] = g.getRealY(); + //std::cout << "to all tree: " << p[0] << ", " << p[1] << std::endl; + dnnPointsAll.push_back(p); + assert(dnnItemIdxAll == dnnPointsAll.size()); + } + + for(size_t i = 0; i < dnnPointsAll.size(); ++i) { + dnnPointHandlesAll.push_back(&dnnPointsAll[i]); + } + //std::cout << "kdtreeAll: " << dnnPointHandlesAll.size() << " points" << std::endl; + kdtreeAll = new dnn::KDTree<DnnTraits>(traits, dnnPointHandlesAll, wassersteinPower); + + size_t handleIdx {0}; + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + if (items[itemIdx].isDiagonal() ) { + heapHandlesIndices[itemIdx] = handleIdx++; + diagHeapHandles.push_back(diagItemsHeap.push(std::make_pair(itemIdx, 0))); + } + } + //to-do: remove maxVal from + //std::cout << "3getFurthestDistance3Approx = " << getFurthestDistance3Approx(_bidders, _items) << std::endl; + maxVal = 3*getFurthestDistance3Approx(_bidders, _items); + maxVal = pow(maxVal, wassersteinPower); + weightAdjConst = maxVal; + //std::cout << "AuctionOracleKDTree: weightAdjConst = " << weightAdjConst << std::endl; + //std::cout << "AuctionOracleKDTree constructor done" << std::endl; + // for debug +} + +DebugOptimalBid AuctionOracleKDTree::getOptimalBidDebug(IdxType bidderIdx) +{ + DebugOptimalBid result; + DiagramPoint bidder = bidders[bidderIdx]; + DnnPoint bidderDnn; + bidderDnn[0] = bidder.getRealX(); + bidderDnn[1] = bidder.getRealY(); + + //std::cout << "bidder.x = " << bidderDnn[0] << std::endl; + //std::cout << "bidder.y = " << bidderDnn[1] << std::endl; + + std::vector<IdxValPair> candItems; + + + if ( bidder.isDiagonal() ) { + // + auto twoBestItems = kdtree->findK(bidderDnn, 2); + //std::cout << "twoBestItems for non-diag: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; + candItems.push_back( std::make_pair(twoBestItems[0].p->id(), twoBestItems[0].d) ); + candItems.push_back( std::make_pair(twoBestItems[1].p->id(), twoBestItems[1].d) ); + assert(diagItemsHeap.size() > 1); + auto topDiagIter = diagItemsHeap.ordered_begin(); + auto topDiag1 = *topDiagIter++; + auto topDiag2 = *topDiagIter; + candItems.push_back(topDiag1); + candItems.push_back(topDiag2); + assert(candItems.size() == 4); + std::sort(candItems.begin(), candItems.end(), CompPairsBySecondStruct()); + assert(candItems[3].second >= candItems[2].second); + assert(candItems[2].second >= candItems[1].second); + assert(candItems[1].second >= candItems[0].second); + } else { + auto twoBestItems = kdtreeAll->findK(bidderDnn, 2); + //std::cout << "twoBestItems for all: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; + candItems.push_back( std::make_pair(twoBestItems[0].p->id(), twoBestItems[0].d) ); + candItems.push_back( std::make_pair(twoBestItems[1].p->id(), twoBestItems[1].d) ); + //size_t projItemIdx { biddersProjIndices.at(bidderIdx) }; + //assert(items[projItemIdx].projId == bidder.id); + //double projItemValue { pow(distLp(bidder, items[projItemIdx]), wassersteinPower) + prices.at(projItemIdx) }; + //candItems.push_back( std::make_pair(projItemIdx, projItemValue) ); + assert(candItems.size() == 2); + assert(candItems[1].second >= candItems[0].second); + } + + result.bestItemIdx = candItems[0].first; + result.secondBestItemIdx = candItems[1].first; + result.bestItemValue = candItems[0].second; + result.secondBestItemValue = candItems[1].second; + //std::cout << "got result: " << result << std::endl; + //double bestItemsPrice = prices[bestItemIdx]; + //if (items[result.bestItemIdx].type == DiagramPoint::DIAG) { + //double bestItemValue1 = pow( distLp(bidder, items[result.bestItemIdx]), q) + prices[result.bestItemIdx]; + //if ( fabs(result.bestItemValue - bestItemValue1) > 1e-6 ) { + //std::cerr << "XXX: " << result.bestItemValue << " vs " << bestItemValue1 << std::endl; + //result.bestItemValue = bestItemValue1; + //} + + //} + + + // checking code + /* + + DebugOptimalBid debugMyResult(result); + DebugOptimalBid debugNaiveResult; + debugNaiveResult.bestItemValue = 1e20; + debugNaiveResult.secondBestItemValue = 1e20; + double currItemValue; + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + //if ( bidders[bidderIdx].type == DiagramPoint::NORMAL and + //items[itemIdx].type == DiagramPoint::DIAG and + //bidders[bidderIdx].projId != items[itemIdx].id) + //continue; + + currItemValue = pow(distLp(bidders[bidderIdx], items[itemIdx]), wassersteinPower) + prices[itemIdx]; + if (currItemValue < debugNaiveResult.bestItemValue) { + debugNaiveResult.bestItemValue = currItemValue; + debugNaiveResult.bestItemIdx = itemIdx; + } + } + + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + if (itemIdx == debugNaiveResult.bestItemIdx) { + continue; + } + currItemValue = pow(distLp(bidders[bidderIdx], items[itemIdx]), wassersteinPower) + prices[itemIdx]; + if (currItemValue < debugNaiveResult.secondBestItemValue) { + debugNaiveResult.secondBestItemValue = currItemValue; + debugNaiveResult.secondBestItemIdx = itemIdx; + } + } + //std::cout << "got naive result" << std::endl; + + if ( fabs( debugMyResult.bestItemValue - debugNaiveResult.bestItemValue ) > 1e-6 or + fabs( debugNaiveResult.secondBestItemValue - debugMyResult.secondBestItemValue) > 1e-6 ) { + kdtreeAll->printWeights(); + std::cerr << "bidderIdx = " << bidderIdx << "; "; + std::cerr << bidders[bidderIdx] << std::endl; + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + std::cout << itemIdx << ": " << items[itemIdx] << "; price = " << prices[itemIdx] << std::endl; + } + std::cerr << "debugMyResult: " << debugMyResult << std::endl; + std::cerr << "debugNaiveResult: " << debugNaiveResult << std::endl; + //std::cerr << "twoBestItems: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; + assert(false); + } + //std::cout << "returning" << std::endl; + + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << secondBestItemIdx << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[secondBestItemIdx] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemsDist= " << (weightAdjConst - bestItemValue) << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << secondBestItemIdx << "; secondBestDist= " << (weightAdjConst - secondBestItemValue) << "; secondBestPrice = " << prices[secondBestItemIdx] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + */ + + return result; +} + +IdxValPair AuctionOracleKDTree::getOptimalBid(IdxType bidderIdx) +{ + IdxValPair result; + DebugOptimalBid debugMyResult = getOptimalBidDebug(bidderIdx); + result.first = debugMyResult.bestItemIdx; + result.second = ( debugMyResult.secondBestItemValue - debugMyResult.bestItemValue ) + prices[debugMyResult.bestItemIdx] + epsilon; + return result; +} +/* +a_{ij} = d_{ij} +value_{ij} = a_{ij} + price_j +*/ +void AuctionOracleKDTree::setPrice(IdxType itemIdx, double newPrice) +{ + assert(prices.size() == items.size()); + assert( 0 < diagHeapHandles.size() and diagHeapHandles.size() <= items.size()); + assert(newPrice > prices.at(itemIdx)); + prices[itemIdx] = newPrice; + if ( items[itemIdx].isNormal() ) { + //std::cout << "before increasing weight in kdtree " << std::endl; + //std::cout << kdtreeItems.at(itemIdx) << std::endl; + assert(0 <= itemIdx and itemIdx < kdtreeItems.size()); + assert(0 <= kdtreeItems[itemIdx] and kdtreeItems[itemIdx] < dnnPointHandles.size()); + kdtree->increase_weight( dnnPointHandles[kdtreeItems[itemIdx]], newPrice); + kdtreeAll->increase_weight( dnnPointHandlesAll[itemIdx], newPrice); + //std::cout << "after increasing weight in kdtree" << std::endl; + } else { + //std::cout << "before decreasing weight in heap" << std::endl; + //std::cout << "diagHeapHandles.size = " << diagHeapHandles.size() << std::endl; + kdtreeAll->increase_weight( dnnPointHandlesAll[itemIdx], newPrice); + //std::cout << "after increasing weight in kdtree" << std::endl; + assert(diagHeapHandles.size() > heapHandlesIndices.at(itemIdx)); + diagItemsHeap.decrease(diagHeapHandles[heapHandlesIndices[itemIdx]], std::make_pair(itemIdx, newPrice)); + } +} + +void AuctionOracleKDTree::adjustPrices(void) +{ +} + +AuctionOracleKDTree::~AuctionOracleKDTree() +{ + delete kdtree; + delete kdtreeAll; +} + +void AuctionOracleKDTree::setEpsilon(double newVal) +{ + assert(newVal >= 0.0); + epsilon = newVal; +} + +// ***************************** +// AuctionOracleRestricted +// ***************************** +AuctionOracleRestricted::AuctionOracleRestricted(const std::vector<DiagramPoint>& b, + const std::vector<DiagramPoint>& g, + double _wassersteinPower, + double internal_p) : + AuctionOracleAbstract(b, g, _wassersteinPower, internal_p), + maxVal(0.0) +{ + assert(b.size() == g.size() ); + assert(b.size() > 1); + + weightMatrix.reserve(b.size()); + for(const auto& pointA : bidders) { + std::vector<double> weightVec; + weightVec.clear(); + weightVec.reserve(b.size()); + for(const auto& pointB : items) { + double val = pow(distLp(pointA, pointB, internal_p), wassersteinPower); + if (val > maxVal) { + maxVal = val; + } + weightVec.push_back( val ); + } + weightMatrix.push_back(weightVec); + } +} + +IdxValPair AuctionOracleRestricted::getOptimalBid(const IdxType bidderIdx) +{ + assert(bidderIdx >=0 and bidderIdx < static_cast<IdxType>(bidders.size()) ); + + const auto bidder = bidders[bidderIdx]; + + IdxType bestItemIdx { -1 }; + double bestItemValue { std::numeric_limits<double>::max() }; + //IdxType secondBestItemIdx { -1 }; + double secondBestItemValue { std::numeric_limits<double>::max() }; + + // find best items idx + for(IdxType itemIdx = 0; itemIdx < static_cast<IdxType>(items.size()); ++itemIdx) { + // non-diagonal point should be matched either to another + // non-diagonal point or to its own projection + if (isRestricted and bidder.isNormal() ) { + auto item = items[itemIdx]; + if (item.isDiagonal() and itemIdx != bidderIdx) + continue; + } + auto currItemValue = weightMatrix[bidderIdx][itemIdx] + prices[itemIdx]; + if ( currItemValue < bestItemValue ) { + bestItemValue = currItemValue; + bestItemIdx = itemIdx; + } + } + + // find second best items idx and value + + for(IdxType itemIdx = 0; itemIdx < static_cast<IdxType>(items.size()); ++itemIdx) { + // non-diagonal point should be matched either to another + // non-diagonal point or to its own projection + if (isRestricted and bidder.isNormal() ) { + auto itemsItem = items[itemIdx]; + if (itemsItem.isDiagonal() and itemIdx != bidderIdx) + continue; + } + + if (static_cast<IdxType>(itemIdx) == bestItemIdx) + continue; + + auto currItemValue = weightMatrix[bidderIdx][itemIdx] + prices[itemIdx]; + if ( currItemValue < secondBestItemValue ) { + secondBestItemValue = currItemValue; + //secondBestItemIdx = itemIdx; + } + } + + assert(bestItemValue <= secondBestItemValue); + + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << secondBestItemIdx << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[secondBestItemIdx] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemsDist= " << (weightAdjConst - bestItemValue) << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << topIter->first << "; secondBestDist= " << (weightAdjConst - secondBestItemValue) << "; secondBestPrice = " << prices[topIter->first] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + + // bid value: price + value difference + epsilon + + return std::make_pair(bestItemIdx, + prices[bestItemIdx] + + ( -bestItemValue + secondBestItemValue ) + + epsilon ); +} + +void AuctionOracleRestricted::setPrice(const IdxType itemIdx, const double newPrice) +{ + assert(prices.at(itemIdx) < newPrice ); + prices[itemIdx] = newPrice; +} + +// ***************************** +// AuctionOracleKDTreeRestricted +// ***************************** + +AuctionOracleKDTreeRestricted::AuctionOracleKDTreeRestricted(const std::vector<DiagramPoint>& _bidders, + const std::vector<DiagramPoint>& _items, + const double _wassersteinPower, + const double internal_p) : + AuctionOracleAbstract(_bidders, _items, _wassersteinPower, internal_p), + heapHandlesIndices(items.size(), std::numeric_limits<size_t>::max()), + kdtreeItems(items.size(), std::numeric_limits<size_t>::max()), + bestDiagonalItemsComputed(false) +{ + size_t dnnItemIdx { 0 }; + size_t trueIdx { 0 }; + dnnPoints.clear(); + // store normal items in kd-tree + for(const auto& g : items) { + if (g.isNormal() ) { + kdtreeItems[trueIdx] = dnnItemIdx; + // index of items is id of dnn-point + DnnPoint p(trueIdx); + p[0] = g.getRealX(); + p[1] = g.getRealY(); + dnnPoints.push_back(p); + assert(dnnItemIdx == dnnPoints.size() - 1); + dnnItemIdx++; + } + trueIdx++; + } + + assert(dnnPoints.size() < items.size() ); + for(size_t i = 0; i < dnnPoints.size(); ++i) { + dnnPointHandles.push_back(&dnnPoints[i]); + } + DnnTraits traits; + //std::cout << "kdtree: " << dnnPointHandles.size() << " points" << std::endl; + kdtree = new dnn::KDTree<DnnTraits>(traits, dnnPointHandles, wassersteinPower); + + size_t handleIdx {0}; + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + if (items[itemIdx].isDiagonal()) { + heapHandlesIndices[itemIdx] = handleIdx++; + diagHeapHandles.push_back(diagItemsHeap.push(std::make_pair(itemIdx, 0))); + } + } + //to-do: remove maxVal from + //std::cout << "3getFurthestDistance3Approx = " << getFurthestDistance3Approx(_bidders, _items) << std::endl; + maxVal = 3*getFurthestDistance3Approx(_bidders, _items); + maxVal = pow(maxVal, wassersteinPower); + weightAdjConst = maxVal; + //std::cout << "AuctionOracleKDTreeRestricted: weightAdjConst = " << weightAdjConst << std::endl; + //std::cout << "AuctionOracleKDTreeRestricted constructor done" << std::endl; +} + +DebugOptimalBid AuctionOracleKDTreeRestricted::getOptimalBidDebug(IdxType bidderIdx) +{ + DebugOptimalBid result; + DiagramPoint bidder = bidders[bidderIdx]; + + //std::cout << "bidder.x = " << bidderDnn[0] << std::endl; + //std::cout << "bidder.y = " << bidderDnn[1] << std::endl; + + // corresponding point is always considered as a candidate + // if bidder is a diagonal point, projItem is a normal point, + // and vice versa. + + size_t projItemIdx = bidderIdx; + assert( 0 <= projItemIdx and projItemIdx < items.size() ); + DiagramPoint projItem = items[projItemIdx]; + assert(projItem.type != bidder.type); + //assert(projItem.projId == bidder.id); + //assert(projItem.id == bidder.projId); + // todo: store precomputed distance? + double projItemValue = pow(distLp(bidder, projItem, internal_p), wassersteinPower) + prices[projItemIdx]; + + if (bidder.isDiagonal()) { + // for diagonal bidder the only normal point has already been added + // the other 2 candidates are diagonal items only, get from the heap + // with prices + assert(diagItemsHeap.size() > 1); + if (!bestDiagonalItemsComputed) { + auto topDiagIter = diagItemsHeap.ordered_begin(); + bestDiagonalItemIdx = topDiagIter->first; + bestDiagonalItemValue = topDiagIter->second; + topDiagIter++; + secondBestDiagonalItemIdx = topDiagIter->first; + secondBestDiagonalItemValue = topDiagIter->second; + bestDiagonalItemsComputed = true; + } + + if ( projItemValue < bestDiagonalItemValue) { + result.bestItemIdx = projItemIdx; + result.bestItemValue = projItemValue; + result.secondBestItemIdx = bestDiagonalItemIdx; + result.secondBestItemValue = bestDiagonalItemValue; + } else if (projItemValue < secondBestDiagonalItemValue) { + result.bestItemIdx = bestDiagonalItemIdx; + result.bestItemValue = bestDiagonalItemValue; + result.secondBestItemIdx = projItemIdx; + result.secondBestItemValue = projItemValue; + } else { + result.bestItemIdx = bestDiagonalItemIdx; + result.bestItemValue = bestDiagonalItemValue; + result.secondBestItemIdx = secondBestDiagonalItemIdx; + result.secondBestItemValue = secondBestDiagonalItemValue; + } + } else { + // for normal bidder get 2 best items among non-diagonal points from + // kdtree + DnnPoint bidderDnn; + bidderDnn[0] = bidder.getRealX(); + bidderDnn[1] = bidder.getRealY(); + auto twoBestItems = kdtree->findK(bidderDnn, 2); + //std::cout << "twoBestItems for all: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; + size_t bestNormalItemIdx { twoBestItems[0].p->id() }; + double bestNormalItemValue { twoBestItems[0].d }; + size_t secondBestNormalItemIdx { twoBestItems[1].p->id() }; + double secondBestNormalItemValue { twoBestItems[1].d }; + + if ( projItemValue < bestNormalItemValue) { + result.bestItemIdx = projItemIdx; + result.bestItemValue = projItemValue; + result.secondBestItemIdx = bestNormalItemIdx; + result.secondBestItemValue = bestNormalItemValue; + } else if (projItemValue < secondBestNormalItemValue) { + result.bestItemIdx = bestNormalItemIdx; + result.bestItemValue = bestNormalItemValue; + result.secondBestItemIdx = projItemIdx; + result.secondBestItemValue = projItemValue; + } else { + result.bestItemIdx = bestNormalItemIdx; + result.bestItemValue = bestNormalItemValue; + result.secondBestItemIdx = secondBestNormalItemIdx; + result.secondBestItemValue = secondBestNormalItemValue; + } + } + + return result; + + //std::cout << "got result: " << result << std::endl; + //double bestItemsPrice = prices[bestItemIdx]; + //if (items[result.bestItemIdx].type == DiagramPoint::DIAG) { + //double bestItemValue1 = pow( distLp(bidder, items[result.bestItemIdx]), wassersteinPower) + prices[result.bestItemIdx]; + //if ( fabs(result.bestItemValue - bestItemValue1) > 1e-6 ) { + //std::cerr << "XXX: " << result.bestItemValue << " vs " << bestItemValue1 << std::endl; + //result.bestItemValue = bestItemValue1; + //} + + //} + + + // checking code + + /* + DebugOptimalBid debugMyResult(result); + DebugOptimalBid debugNaiveResult; + debugNaiveResult.bestItemValue = 1e20; + debugNaiveResult.secondBestItemValue = 1e20; + double currItemValue; + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + //if ( bidders[bidderIdx].type == DiagramPoint::NORMAL and + //items[itemIdx].type == DiagramPoint::DIAG and + //bidders[bidderIdx].projId != items[itemIdx].id) + //continue; + + currItemValue = pow(distLp(bidders[bidderIdx], items[itemIdx]), wassersteinPower) + prices[itemIdx]; + if (currItemValue < debugNaiveResult.bestItemValue) { + debugNaiveResult.bestItemValue = currItemValue; + debugNaiveResult.bestItemIdx = itemIdx; + } + } + + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + if (itemIdx == debugNaiveResult.bestItemIdx) { + continue; + } + currItemValue = pow(distLp(bidders[bidderIdx], items[itemIdx]), wassersteinPower) + prices[itemIdx]; + if (currItemValue < debugNaiveResult.secondBestItemValue) { + debugNaiveResult.secondBestItemValue = currItemValue; + debugNaiveResult.secondBestItemIdx = itemIdx; + } + } + //std::cout << "got naive result" << std::endl; + + if ( fabs( debugMyResult.bestItemValue - debugNaiveResult.bestItemValue ) > 1e-6 or + fabs( debugNaiveResult.secondBestItemValue - debugMyResult.secondBestItemValue) > 1e-6 ) { + kdtreeAll->printWeights(); + std::cerr << "bidderIdx = " << bidderIdx << "; "; + std::cerr << bidders[bidderIdx] << std::endl; + for(size_t itemIdx = 0; itemIdx < items.size(); ++itemIdx) { + std::cout << itemIdx << ": " << items[itemIdx] << "; price = " << prices[itemIdx] << std::endl; + } + std::cerr << "debugMyResult: " << debugMyResult << std::endl; + std::cerr << "debugNaiveResult: " << debugNaiveResult << std::endl; + //std::cerr << "twoBestItems: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; + assert(false); + } + //std::cout << "returning" << std::endl; + + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemValue = " << bestItemValue << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << secondBestItemIdx << "; secondBestValue = " << secondBestItemValue << "; secondBestPrice = " << prices[secondBestItemIdx] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + //std::cout << "getOptimalBid: bidderIdx = " << bidderIdx << "; bestItemIdx = " << bestItemIdx << "; bestItemsDist= " << (weightAdjConst - bestItemValue) << "; bestItemsPrice = " << prices[bestItemIdx] << "; secondBestItemIdx = " << secondBestItemIdx << "; secondBestDist= " << (weightAdjConst - secondBestItemValue) << "; secondBestPrice = " << prices[secondBestItemIdx] << "; bid = " << prices[bestItemIdx] + ( bestItemValue - secondBestItemValue ) + epsilon << "; epsilon = " << epsilon << std::endl; + */ + return result; +} + +IdxValPair AuctionOracleKDTreeRestricted::getOptimalBid(IdxType bidderIdx) +{ + + + DiagramPoint bidder = bidders[bidderIdx]; + + //std::cout << "bidder.x = " << bidderDnn[0] << std::endl; + //std::cout << "bidder.y = " << bidderDnn[1] << std::endl; + + // corresponding point is always considered as a candidate + // if bidder is a diagonal point, projItem is a normal point, + // and vice versa. + + size_t bestItemIdx; + double bestItemValue; + double secondBestItemValue; + + + size_t projItemIdx = bidderIdx; + assert( 0 <= projItemIdx and projItemIdx < items.size() ); + DiagramPoint projItem = items[projItemIdx]; + assert(projItem.type != bidder.type); + //assert(projItem.projId == bidder.id); + //assert(projItem.id == bidder.projId); + // todo: store precomputed distance? + double projItemValue = pow(distLp(bidder, projItem, internal_p), wassersteinPower) + prices[projItemIdx]; + + if (bidder.isDiagonal()) { + // for diagonal bidder the only normal point has already been added + // the other 2 candidates are diagonal items only, get from the heap + // with prices + + if (not bestDiagonalItemsComputed) { + auto topDiagIter = diagItemsHeap.ordered_begin(); + bestDiagonalItemIdx = topDiagIter->first; + bestDiagonalItemValue = topDiagIter->second; + if (diagItemsHeap.size() > 1) { + topDiagIter++; + secondBestDiagonalItemIdx = topDiagIter->first; + secondBestDiagonalItemValue = topDiagIter->second; + } else { + // there is only one diagonal point at all, + // ensure that second best diagonal value + // will lose to projection item + secondBestDiagonalItemValue = std::numeric_limits<double>::max(); + secondBestDiagonalItemIdx = -1; + } + bestDiagonalItemsComputed = true; + } + + if ( projItemValue < bestDiagonalItemValue) { + bestItemIdx = projItemIdx; + bestItemValue = projItemValue; + secondBestItemValue = bestDiagonalItemValue; + } else if (projItemValue < secondBestDiagonalItemValue) { + bestItemIdx = bestDiagonalItemIdx; + bestItemValue = bestDiagonalItemValue; + secondBestItemValue = projItemValue; + } else { + bestItemIdx = bestDiagonalItemIdx; + bestItemValue = bestDiagonalItemValue; + secondBestItemValue = secondBestDiagonalItemValue; + } + } else { + // for normal bidder get 2 best items among non-diagonal points from + // kdtree + DnnPoint bidderDnn; + bidderDnn[0] = bidder.getRealX(); + bidderDnn[1] = bidder.getRealY(); + auto twoBestItems = kdtree->findK(bidderDnn, 2); + //std::cout << "twoBestItems for all: " << twoBestItems[0].d << " " << twoBestItems[1].d << std::endl; + size_t bestNormalItemIdx { twoBestItems[0].p->id() }; + double bestNormalItemValue { twoBestItems[0].d }; + // if there is only one off-diagonal point in the second diagram, + // kd-tree will not return the second candidate. + // Set its value to inf, so it will always lose to the value of the projection + double secondBestNormalItemValue { twoBestItems.size() == 1 ? std::numeric_limits<double>::max() : twoBestItems[1].d }; + + if ( projItemValue < bestNormalItemValue) { + bestItemIdx = projItemIdx; + bestItemValue = projItemValue; + secondBestItemValue = bestNormalItemValue; + } else if (projItemValue < secondBestNormalItemValue) { + bestItemIdx = bestNormalItemIdx; + bestItemValue = bestNormalItemValue; + secondBestItemValue = projItemValue; + } else { + bestItemIdx = bestNormalItemIdx; + bestItemValue = bestNormalItemValue; + secondBestItemValue = secondBestNormalItemValue; + } + } + + IdxValPair result; + + assert( secondBestItemValue >= bestItemValue ); + + result.first = bestItemIdx; + result.second = ( secondBestItemValue - bestItemValue ) + prices[bestItemIdx] + epsilon; + return result; +} +/* +a_{ij} = d_{ij} +value_{ij} = a_{ij} + price_j +*/ +void AuctionOracleKDTreeRestricted::setPrice(IdxType itemIdx, double newPrice) +{ + assert(prices.size() == items.size()); + assert( 0 < diagHeapHandles.size() and diagHeapHandles.size() <= items.size()); + assert(newPrice > prices.at(itemIdx)); + prices[itemIdx] = newPrice; + if ( items[itemIdx].isNormal() ) { + //std::cout << "before increasing weight in kdtree " << std::endl; + //std::cout << kdtreeItems.at(itemIdx) << std::endl; + assert(0 <= itemIdx and itemIdx < kdtreeItems.size()); + assert(0 <= kdtreeItems[itemIdx] and kdtreeItems[itemIdx] < dnnPointHandles.size()); + kdtree->increase_weight( dnnPointHandles[kdtreeItems[itemIdx]], newPrice); + //std::cout << "after increasing weight in kdtree" << std::endl; + } else { + //std::cout << "before decreasing weight in heap" << std::endl; + //std::cout << "diagHeapHandles.size = " << diagHeapHandles.size() << std::endl; + assert(diagHeapHandles.size() > heapHandlesIndices.at(itemIdx)); + diagItemsHeap.decrease(diagHeapHandles[heapHandlesIndices[itemIdx]], std::make_pair(itemIdx, newPrice)); + bestDiagonalItemsComputed = false; + } +} + +void AuctionOracleKDTreeRestricted::adjustPrices(void) +{ +} + +AuctionOracleKDTreeRestricted::~AuctionOracleKDTreeRestricted() +{ + delete kdtree; +} + +void AuctionOracleKDTreeRestricted::setEpsilon(double newVal) +{ + assert(newVal >= 0.0); + epsilon = newVal; +} + +std::ostream& operator<< (std::ostream& output, const DebugOptimalBid& db) +{ + std::cout << "bestItemValue = " << db.bestItemValue; + std::cout << "; bestItemIdx = " << db.bestItemIdx; + std::cout << "; secondBestItemValue = " << db.secondBestItemValue; + std::cout << "; secondBestItemIdx = " << db.secondBestItemIdx; + return output; +} + +} // end of namespace geom_ws diff --git a/geom_matching/wasserstein/src/auction_runner_gs.cpp b/geom_matching/wasserstein/src/auction_runner_gs.cpp new file mode 100644 index 0000000..10d37c9 --- /dev/null +++ b/geom_matching/wasserstein/src/auction_runner_gs.cpp @@ -0,0 +1,341 @@ +/* + +Copyright (c) 2016, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + + +#include <assert.h> +#include <algorithm> +#include <functional> +#include <iterator> +#include <chrono> + +#include "def_debug.h" +#include "auction_runner_gs.h" +#include "wasserstein.h" + +//#define PRINT_DETAILED_TIMING + +namespace geom_ws { + +// ***************************** +// AuctionRunnerGS +// ***************************** + +AuctionRunnerGS::AuctionRunnerGS(const std::vector<DiagramPoint>& A, const std::vector<DiagramPoint>& B, const double q, const double _delta, const double _internal_p, const double _initialEpsilon, const double _epsFactor) : + bidders(A), + items(B), + numBidders(A.size()), + numItems(A.size()), + itemsToBidders(A.size(), -1), + biddersToItems(A.size(), -1), + wassersteinPower(q), + delta(_delta), + internal_p(_internal_p), + initialEpsilon(_initialEpsilon), + epsilonCommonRatio(_epsFactor == 0.0 ? 5.0 : _epsFactor) +{ + assert(initialEpsilon >= 0.0 ); + assert(epsilonCommonRatio >= 0.0 ); + assert(A.size() == B.size()); + oracle = std::unique_ptr<AuctionOracle>(new AuctionOracle(bidders, items, wassersteinPower, internal_p)); +} + +void AuctionRunnerGS::assignItemToBidder(IdxType itemIdx, IdxType bidderIdx) +{ + numRounds++; + //sanityCheck(); + // only unassigned bidders should submit bids and get items + assert(biddersToItems[bidderIdx] == -1); + IdxType oldItemOwner = itemsToBidders[itemIdx]; + + // set new owner + biddersToItems[bidderIdx] = itemIdx; + itemsToBidders[itemIdx] = bidderIdx; + // remove bidder from the list of unassigned bidders +#ifdef KEEP_UNASSIGNED_ORDERED + unassignedBidders.erase(std::make_pair(bidderIdx, bidders[bidderIdx])); +#else + unassignedBidders.erase(bidderIdx); +#endif + + // old owner becomes unassigned + if (oldItemOwner != -1) { + biddersToItems[oldItemOwner] = -1; +#ifdef KEEP_UNASSIGNED_ORDERED + unassignedBidders.insert(std::make_pair(oldItemOwner, bidders[oldItemOwner])); +#else + unassignedBidders.insert(oldItemOwner); +#endif + } +} + + +void AuctionRunnerGS::flushAssignment(void) +{ + for(auto& b2i : biddersToItems) { + b2i = -1; + } + for(auto& i2b : itemsToBidders) { + i2b = -1; + } + // we must flush assignment only after we got perfect matching + assert(unassignedBidders.empty()); + // all bidders become unassigned + for(size_t bidderIdx = 0; bidderIdx < numBidders; ++bidderIdx) { +#ifdef KEEP_UNASSIGNED_ORDERED + unassignedBidders.insert(std::make_pair(bidderIdx, bidders[bidderIdx])); +#else + unassignedBidders.insert(bidderIdx); +#endif + } + assert(unassignedBidders.size() == bidders.size()); + //oracle->adjustPrices(); +} + +void AuctionRunnerGS::runAuction(void) +{ +#ifdef PRINT_DETAILED_TIMING + std::chrono::high_resolution_clock hrClock; + std::chrono::time_point<std::chrono::high_resolution_clock> startMoment; + startMoment = hrClock.now(); + std::vector<double> iterResults; + std::vector<double> iterEstRelErrors; + std::vector<std::chrono::time_point<std::chrono::high_resolution_clock>> iterTimes; +#endif + // choose some initial epsilon + if (initialEpsilon == 0.0) + oracle->setEpsilon(oracle->maxVal / 4.0); + else + oracle->setEpsilon(initialEpsilon); + assert( oracle->getEpsilon() > 0 ); + int iterNum { 0 }; + bool notDone { false }; + double currentResult; + do { + flushAssignment(); + runAuctionPhase(); + iterNum++; + //std::cout << "Iteration " << iterNum << " completed. " << std::endl; + // result is d^q + currentResult = getDistanceToQthPowerInternal(); + double denominator = currentResult - numBidders * oracle->getEpsilon(); + currentResult = pow(currentResult, 1.0 / wassersteinPower); +#ifdef PRINT_DETAILED_TIMING + iterResults.push_back(currentResult); + iterTimes.push_back(hrClock.now()); + std::cout << "Iteration " << iterNum << " finished. "; + std::cout << "Current result is " << currentResult << ", epsilon = " << oracle->getEpsilon() << std::endl; + std::cout << "Number of rounds (cumulative): " << numRounds << std::endl; +#endif + if ( denominator <= 0 ) { + //std::cout << "Epsilon is too big." << std::endl; + notDone = true; + } else { + denominator = pow(denominator, 1.0 / wassersteinPower); + double numerator = currentResult - denominator; +#ifdef PRINT_DETAILED_TIMING + std::cout << " numerator: " << numerator << " denominator: " << denominator; + std::cout << "; error bound: " << numerator / denominator << std::endl; +#endif + // if relative error is greater than delta, continue + notDone = ( numerator / denominator > delta ); + } + // decrease epsilon for the next iteration + oracle->setEpsilon( oracle->getEpsilon() / epsilonCommonRatio ); + if (iterNum > maxIterNum) { + std::cerr << "Maximum iteration number exceeded, exiting. Current result is:"; + std::cerr << wassersteinDistance << std::endl; + std::exit(1); + } + } while ( notDone ); + //printMatching(); +#ifdef PRINT_DETAILED_TIMING + for(size_t iterIdx = 0; iterIdx < iterResults.size(); ++iterIdx) { + double trueRelError = ( iterResults.at(iterIdx) - currentResult ) / currentResult; + auto iterCumulativeTime = iterTimes.at(iterIdx) - startMoment; + std::chrono::duration<double, std::milli> iterTime = ( iterIdx > 0) ? iterTimes[iterIdx] - iterTimes[iterIdx - 1] : iterTimes[iterIdx] - startMoment; + std::cout << "iteration " << iterIdx << ", true rel. error " << + trueRelError << ", elapsed time " << + std::chrono::duration<double, std::milli>(iterCumulativeTime).count() << + ", iteration time " << iterTime.count() << std::endl; + } +#endif +} + +void AuctionRunnerGS::runAuctionPhase(void) +{ + //std::cout << "Entered runAuctionPhase" << std::endl; + do { +#ifdef KEEP_UNASSIGNED_ORDERED + size_t bidderIdx = unassignedBidders.begin()->first; +#else + size_t bidderIdx = *unassignedBidders.begin(); +#endif + auto optimalBid = oracle->getOptimalBid(bidderIdx); + auto optimalItemIdx = optimalBid.first; + auto bidValue = optimalBid.second; + assignItemToBidder(optimalBid.first, bidderIdx); + oracle->setPrice(optimalItemIdx, bidValue); + //printDebug(); + } while (not unassignedBidders.empty()); + //std::cout << "runAuctionPhase finished" << std::endl; + +#ifdef DEBUG_AUCTION + for(size_t bidderIdx = 0; bidderIdx < numBidders; ++bidderIdx) { + if ( biddersToItems[bidderIdx] < 0) { + std::cerr << "After auction terminated bidder " << bidderIdx; + std::cerr << " has no items assigned" << std::endl; + throw "Auction did not give a perfect matching"; + } + } +#endif + +} + +double AuctionRunnerGS::getDistanceToQthPowerInternal(void) +{ + sanityCheck(); + double result = 0.0; + for(size_t bIdx = 0; bIdx < numBidders; ++bIdx) { + auto pA = bidders[bIdx]; + assert( 0 <= biddersToItems[bIdx] and biddersToItems[bIdx] < static_cast<int>(items.size()) ); + auto pB = items[biddersToItems[bIdx]]; + result += pow(distLp(pA, pB, internal_p), wassersteinPower); + } + wassersteinDistance = pow(result, 1.0 / wassersteinPower); + return result; +} + +double AuctionRunnerGS::getWassersteinDistance(void) +{ + runAuction(); + return wassersteinDistance; +} + +// Debug routines + +void AuctionRunnerGS::printDebug(void) +{ +#ifdef DEBUG_AUCTION + sanityCheck(); + std::cout << "**********************" << std::endl; + std::cout << "Current assignment:" << std::endl; + for(size_t idx = 0; idx < biddersToItems.size(); ++idx) { + std::cout << idx << " <--> " << biddersToItems[idx] << std::endl; + } + std::cout << "Weights: " << std::endl; + //for(size_t i = 0; i < numBidders; ++i) { + //for(size_t j = 0; j < numItems; ++j) { + //std::cout << oracle->weightMatrix[i][j] << " "; + //} + //std::cout << std::endl; + //} + std::cout << "Prices: " << std::endl; + for(const auto price : oracle->getPrices()) { + std::cout << price << std::endl; + } + std::cout << "**********************" << std::endl; +#endif +} + + +void AuctionRunnerGS::sanityCheck(void) +{ +#ifdef DEBUG_AUCTION + if (biddersToItems.size() != numBidders) { + std::cerr << "Wrong size of biddersToItems, must be " << numBidders << ", is " << biddersToItems.size() << std::endl; + throw "Wrong size of biddersToItems"; + } + + if (itemsToBidders.size() != numBidders) { + std::cerr << "Wrong size of itemsToBidders, must be " << numBidders << ", is " << itemsToBidders.size() << std::endl; + throw "Wrong size of itemsToBidders"; + } + + for(size_t bidderIdx = 0; bidderIdx < numBidders; ++bidderIdx) { + if ( biddersToItems[bidderIdx] >= 0) { + + if ( std::count(biddersToItems.begin(), + biddersToItems.end(), + biddersToItems[bidderIdx]) > 1 ) { + std::cerr << "Item " << biddersToItems[bidderIdx]; + std::cerr << " appears in biddersToItems more than once" << std::endl; + throw "Duplicate in biddersToItems"; + } + + if (itemsToBidders.at(biddersToItems[bidderIdx]) != static_cast<int>(bidderIdx)) { + std::cerr << "Inconsitency: bidderIdx = " << bidderIdx; + std::cerr << ", itemIdx in biddersToItems = "; + std::cerr << biddersToItems[bidderIdx]; + std::cerr << ", bidderIdx in itemsToBidders = "; + std::cerr << itemsToBidders[biddersToItems[bidderIdx]] << std::endl; + throw "inconsistent mapping"; + } + } + } + + for(IdxType itemIdx = 0; itemIdx < static_cast<IdxType>(numBidders); ++itemIdx) { + if ( itemsToBidders[itemIdx] >= 0) { + + // check for uniqueness + if ( std::count(itemsToBidders.begin(), + itemsToBidders.end(), + itemsToBidders[itemIdx]) > 1 ) { + std::cerr << "Bidder " << itemsToBidders[itemIdx]; + std::cerr << " appears in itemsToBidders more than once" << std::endl; + throw "Duplicate in itemsToBidders"; + } + // check for consistency + if (biddersToItems.at(itemsToBidders[itemIdx]) != static_cast<int>(itemIdx)) { + std::cerr << "Inconsitency: itemIdx = " << itemIdx; + std::cerr << ", bidderIdx in itemsToBidders = "; + std::cerr << itemsToBidders[itemIdx]; + std::cerr << ", itemIdx in biddersToItems= "; + std::cerr << biddersToItems[itemsToBidders[itemIdx]] << std::endl; + throw "inconsistent mapping"; + } + } + } +#endif +} + +void AuctionRunnerGS::printMatching(void) +{ +//#ifdef DEBUG_AUCTION + sanityCheck(); + for(size_t bIdx = 0; bIdx < biddersToItems.size(); ++bIdx) { + if (biddersToItems[bIdx] >= 0) { + auto pA = bidders[bIdx]; + auto pB = items[biddersToItems[bIdx]]; + std::cout << pA << " <-> " << pB << "+" << pow(distLp(pA, pB, internal_p), wassersteinPower) << std::endl; + } else { + assert(false); + } + } +//#endif +} + +} // end of namespace geom_ws diff --git a/geom_matching/wasserstein/src/auction_runner_jac.cpp b/geom_matching/wasserstein/src/auction_runner_jac.cpp new file mode 100644 index 0000000..dcade94 --- /dev/null +++ b/geom_matching/wasserstein/src/auction_runner_jac.cpp @@ -0,0 +1,365 @@ +/* + +Copyright (c) 2016, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#include <assert.h> +#include <algorithm> +#include <functional> +#include <iterator> + +#include "def_debug.h" +#include "auction_runner_jac.h" +#include "wasserstein.h" + +namespace geom_ws { + +// ***************************** +// AuctionRunnerJak +// ***************************** + +AuctionRunnerJak::AuctionRunnerJak(const std::vector<DiagramPoint>& A, const std::vector<DiagramPoint>& B, const double q, const double _delta, const double _internal_p) : + bidders(A), + items(B), + numBidders(A.size()), + numItems(A.size()), + itemsToBidders(A.size(), -1), + biddersToItems(A.size(), -1), + wassersteinPower(q), + delta(_delta), + internal_p(_internal_p), + bidTable(A.size(), std::make_pair(-1, std::numeric_limits<double>::lowest()) ), + itemReceivedBidVec(B.size(), 0 ) +{ + assert(A.size() == B.size()); + oracle = std::unique_ptr<AuctionOracle>(new AuctionOracle(bidders, items, wassersteinPower, internal_p)); +} + +void AuctionRunnerJak::assignGoodToBidder(IdxType itemIdx, IdxType bidderIdx) +{ + //sanityCheck(); + IdxType myOldItem = biddersToItems[bidderIdx]; + IdxType currItemOwner = itemsToBidders[itemIdx]; + + // set new owner + biddersToItems[bidderIdx] = itemIdx; + itemsToBidders[itemIdx] = bidderIdx; + + + // remove bidder from the list of unassigned bidders + unassignedBidders.erase( unassignedBiddersIterators[bidderIdx] ); + assert( 0 <= bidderIdx and bidderIdx < unassignedBiddersIterators.size() ); + unassignedBiddersIterators[bidderIdx] = unassignedBidders.end(); + + if (-1 == currItemOwner) { + // the item we want to assign does not belong to anybody, + // just free myOldItem, if necessary + // RE: this cannot be necessary. I submitted the best bid, hence I was + // an unassigned bidder. + if (myOldItem != -1) { + std::cout << "This is not happening" << std::endl; + assert(false); + itemsToBidders[myOldItem] = -1; + } + } else { + // the current owner of itemIdx gets my old item (OK if it's -1) + biddersToItems[currItemOwner] = myOldItem; + // add the old owner of bids to the list of + if ( -1 != myOldItem ) { + std::cout << "This is not happening" << std::endl; + assert(false); + // if I had something, update itemsToBidders, too + // RE: nonsense: if I had something, I am not unassigned and did not + // sumbit any bid + itemsToBidders[myOldItem] = currItemOwner; + } + unassignedBidders.push_back(currItemOwner); + assert( unassignedBiddersIterators[currItemOwner] == unassignedBidders.end() ); + unassignedBiddersIterators[currItemOwner] = std::prev( unassignedBidders.end() ); + } + //sanityCheck(); +} + + +void AuctionRunnerJak::assignToBestBidder(IdxType itemIdx) +{ + assert( itemIdx >= 0 and itemIdx < static_cast<IdxType>(numItems) ); + assert( bidTable[itemIdx].first != -1); + IdxValPair bestBid { bidTable[itemIdx] }; + assignGoodToBidder(itemIdx, bestBid.first); + //std::cout << "About to call setPrice" << std::endl; + oracle->setPrice(itemIdx, bestBid.second); + //dynamic_cast<AuctionOracleKDTree*>(oracle)->setNai +} + +void AuctionRunnerJak::clearBidTable(void) +{ + for(auto& itemWithBidIdx : itemsWithBids) { + itemReceivedBidVec[itemWithBidIdx] = 0; + bidTable[itemWithBidIdx].first = -1; + bidTable[itemWithBidIdx].second = std::numeric_limits<double>::lowest(); + } + itemsWithBids.clear(); +} + +void AuctionRunnerJak::submitBid(IdxType bidderIdx, const IdxValPair& itemsBidValuePair) +{ + IdxType itemIdx = itemsBidValuePair.first; + double bidValue = itemsBidValuePair.second; + assert( itemIdx >= 0 ); + if ( bidTable[itemIdx].second < itemsBidValuePair.second ) { + bidTable[itemIdx].first = bidderIdx; + bidTable[itemIdx].second = bidValue; + } + if (0 == itemReceivedBidVec[itemIdx]) { + itemReceivedBidVec[itemIdx] = 1; + itemsWithBids.push_back(itemIdx); + } +} + +void AuctionRunnerJak::printDebug(void) +{ +#ifdef DEBUG_AUCTION + sanityCheck(); + std::cout << "**********************" << std::endl; + std::cout << "Current assignment:" << std::endl; + for(size_t idx = 0; idx < biddersToItems.size(); ++idx) { + std::cout << idx << " <--> " << biddersToItems[idx] << std::endl; + } + std::cout << "Weights: " << std::endl; + //for(size_t i = 0; i < numBidders; ++i) { + //for(size_t j = 0; j < numItems; ++j) { + //std::cout << oracle->weightMatrix[i][j] << " "; + //} + //std::cout << std::endl; + //} + std::cout << "Prices: " << std::endl; + for(const auto price : oracle->getPrices()) { + std::cout << price << std::endl; + } + //std::cout << "Value matrix: " << std::endl; + //for(size_t i = 0; i < numBidders; ++i) { + //for(size_t j = 0; j < numItems; ++j) { + //std::cout << oracle->weightMatrix[i][j] - oracle->prices[j] << " "; + //} + //std::cout << std::endl; + //} + std::cout << "**********************" << std::endl; +#endif +} + +void AuctionRunnerJak::flushAssignment(void) +{ + for(auto& b2g : biddersToItems) { + b2g = -1; + } + for(auto& g2b : itemsToBidders) { + g2b = -1; + } + //oracle->adjustPrices(); +} + +void AuctionRunnerJak::runAuction(void) +{ + // relative error + // choose some initial epsilon + oracle->setEpsilon(oracle->maxVal / 4.0); + assert( oracle->getEpsilon() > 0 ); + int iterNum { 0 }; + bool notDone { false }; + do { + flushAssignment(); + runAuctionPhase(); + iterNum++; + //std::cout << "Iteration " << iterNum << " completed. " << std::endl; + // result is d^q + double currentResult = getDistanceToQthPowerInternal(); + double denominator = currentResult - numBidders * oracle->getEpsilon(); + currentResult = pow(currentResult, 1.0 / wassersteinPower); + //std::cout << "Current result is " << currentResult << std::endl; + if ( denominator <= 0 ) { + //std::cout << "Epsilon is too big." << std::endl; + notDone = true; + } else { + denominator = pow(denominator, 1.0 / wassersteinPower); + double numerator = currentResult - denominator; + //std::cout << " numerator: " << numerator << " denominator: " << denominator << std::endl; + //std::cout << " error bound: " << numerator / denominator << std::endl; + // if relative error is greater than delta, continue + notDone = ( numerator / denominator > delta ); + } + // decrease epsilon for the next iteration + oracle->setEpsilon( oracle->getEpsilon() / epsilonCommonRatio ); + if (iterNum > maxIterNum) { + std::cerr << "Maximum iteration number exceeded, exiting. Current result is:"; + std::cerr << wassersteinDistance << std::endl; + std::exit(1); + } + } while ( notDone ); + //printMatching(); +} + +void AuctionRunnerJak::runAuctionPhase(void) +{ + //std::cout << "Entered runAuctionPhase" << std::endl; + //int numUnassignedBidders { 0 }; + + // at the beginning of a phase all bidders are unassigned + unassignedBidders.clear(); + unassignedBiddersIterators.clear(); + for(size_t bidderIdx = 0; bidderIdx < numBidders; ++bidderIdx) { + unassignedBidders.push_back(bidderIdx); + unassignedBiddersIterators.push_back( std::prev( unassignedBidders.end() )); + } + do { + // bidding phase + clearBidTable(); + for(const auto bidderIdx : unassignedBidders) { + submitBid(bidderIdx, oracle->getOptimalBid(bidderIdx)); + } + //std::cout << "Number of unassignedBidders: " << unassignedBidders.size() << std::endl; + + // assignment phase + // todo: maintain list of items that received a bid + for(auto itemIdx : itemsWithBids ) { + assignToBestBidder(itemIdx); + } + //std::cout << "Assignment phase done" << std::endl; + //sanityCheck(); + //printDebug(); + } while (unassignedBidders.size() > 0); + //std::cout << "runAuctionPhase finished" << std::endl; + + +#ifdef DEBUG_AUCTION + for(size_t bidderIdx = 0; bidderIdx < numBidders; ++bidderIdx) { + if ( biddersToItems[bidderIdx] < 0) { + std::cerr << "After auction terminated bidder " << bidderIdx; + std::cerr << " has no items assigned" << std::endl; + throw "Auction did not give a perfect matching"; + } + } +#endif + +} + +// assertion: the matching must be perfect +double AuctionRunnerJak::getDistanceToQthPowerInternal(void) +{ + sanityCheck(); + double result = 0.0; + for(size_t bIdx = 0; bIdx < numBidders; ++bIdx) { + auto pA = bidders[bIdx]; + assert( 0 <= biddersToItems[bIdx] and biddersToItems[bIdx] < items.size() ); + auto pB = items[biddersToItems[bIdx]]; + result += pow(distLp(pA, pB, internal_p), wassersteinPower); + } + wassersteinDistance = pow(result, 1.0 / wassersteinPower); + return result; +} + +double AuctionRunnerJak::getWassersteinDistance(void) +{ + runAuction(); + return wassersteinDistance; +} + +void AuctionRunnerJak::sanityCheck(void) +{ +#ifdef DEBUG_AUCTION + if (biddersToItems.size() != numBidders) { + std::cerr << "Wrong size of biddersToItems, must be " << numBidders << ", is " << biddersToItems.size() << std::endl; + throw "Wrong size of biddersToItems"; + } + + if (itemsToBidders.size() != numBidders) { + std::cerr << "Wrong size of itemsToBidders, must be " << numBidders << ", is " << itemsToBidders.size() << std::endl; + throw "Wrong size of itemsToBidders"; + } + + for(size_t bidderIdx = 0; bidderIdx < numBidders; ++bidderIdx) { + if ( biddersToItems[bidderIdx] >= 0) { + + if ( std::count(biddersToItems.begin(), + biddersToItems.end(), + biddersToItems[bidderIdx]) > 1 ) { + std::cerr << "Good " << biddersToItems[bidderIdx]; + std::cerr << " appears in biddersToItems more than once" << std::endl; + throw "Duplicate in biddersToItems"; + } + + if (itemsToBidders.at(biddersToItems[bidderIdx]) != static_cast<int>(bidderIdx)) { + std::cerr << "Inconsitency: bidderIdx = " << bidderIdx; + std::cerr << ", itemIdx in biddersToItems = "; + std::cerr << biddersToItems[bidderIdx]; + std::cerr << ", bidderIdx in itemsToBidders = "; + std::cerr << itemsToBidders[biddersToItems[bidderIdx]] << std::endl; + throw "inconsistent mapping"; + } + } + } + + for(IdxType itemIdx = 0; itemIdx < static_cast<IdxType>(numBidders); ++itemIdx) { + if ( itemsToBidders[itemIdx] >= 0) { + + // check for uniqueness + if ( std::count(itemsToBidders.begin(), + itemsToBidders.end(), + itemsToBidders[itemIdx]) > 1 ) { + std::cerr << "Bidder " << itemsToBidders[itemIdx]; + std::cerr << " appears in itemsToBidders more than once" << std::endl; + throw "Duplicate in itemsToBidders"; + } + // check for consistency + if (biddersToItems.at(itemsToBidders[itemIdx]) != static_cast<int>(itemIdx)) { + std::cerr << "Inconsitency: itemIdx = " << itemIdx; + std::cerr << ", bidderIdx in itemsToBidders = "; + std::cerr << itemsToBidders[itemIdx]; + std::cerr << ", itemIdx in biddersToItems= "; + std::cerr << biddersToItems[itemsToBidders[itemIdx]] << std::endl; + throw "inconsistent mapping"; + } + } + } +#endif +} + +void AuctionRunnerJak::printMatching(void) +{ +//#ifdef DEBUG_AUCTION + sanityCheck(); + for(size_t bIdx = 0; bIdx < biddersToItems.size(); ++bIdx) { + if (biddersToItems[bIdx] >= 0) { + auto pA = bidders[bIdx]; + auto pB = items[biddersToItems[bIdx]]; + std::cout << pA << " <-> " << pB << "+" << pow(distLp(pA, pB, internal_p), wassersteinPower) << std::endl; + } else { + assert(false); + } + } +//#endif +} + +} // end of namespace geom_ws diff --git a/geom_matching/wasserstein/src/basic_defs.cpp b/geom_matching/wasserstein/src/basic_defs.cpp new file mode 100644 index 0000000..d228123 --- /dev/null +++ b/geom_matching/wasserstein/src/basic_defs.cpp @@ -0,0 +1,138 @@ +/* + +Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + +*/ + +#include <algorithm> +#include <cfloat> +#include <set> +#include "basic_defs_ws.h" + +namespace geom_ws { +// Point + +bool Point::operator==(const Point& other) const +{ + return ((this->x == other.x) and (this->y == other.y)); +} + +bool Point::operator!=(const Point& other) const +{ + return !(*this == other); +} + +std::ostream& operator<<(std::ostream& output, const Point p) +{ + output << "(" << p.x << ", " << p.y << ")"; + return output; +} + +double sqrDist(const Point& a, const Point& b) +{ + return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y); +} + +double dist(const Point& a, const Point& b) +{ + return sqrt(sqrDist(a, b)); +} + +// DiagramPoint + +// compute l-inf distance between two diagram points +double distLInf(const DiagramPoint& a, const DiagramPoint& b) +{ + if (a.isDiagonal() and b.isDiagonal()) { + return 0.0; + } else { + return std::max(fabs(a.getRealX() - b.getRealX()), fabs(a.getRealY() - b.getRealY())); + } +} + +double distLp(const DiagramPoint& a, const DiagramPoint& b, const double p) +{ + // infinity: special case + if ( std::isinf(p) ) + return distLInf(a, b); + + // check p + assert( p >= 1.0 ); + + // avoid calling pow function + if ( p == 1.0 ) { + if ( a.isNormal() or b.isNormal() ) { + // distance between normal points is a usual l-inf distance + return fabs(a.getRealX() - b.getRealX()) + fabs(a.getRealY() - b.getRealY()); + } else + return 0.0; + } + + if ( a.isNormal() or b.isNormal() ) { + // distance between normal points is a usual l-inf distance + return std::pow(std::pow(fabs(a.getRealX() - b.getRealX()), p) + std::pow(fabs(a.getRealY() - b.getRealY()), p), 1.0/p ); + } else + return 0.0; +} + + +std::ostream& operator<<(std::ostream& output, const DiagramPoint p) +{ + if ( p.type == DiagramPoint::DIAG ) { + output << "(" << p.x << ", " << p.y << ", " << 0.5 * (p.x + p.y) << " DIAG )"; + } else { + output << "(" << p.x << ", " << p.y << ", " << " NORMAL)"; + } + return output; +} + + +DiagramPoint::DiagramPoint(double xx, double yy, Type ttype) : + x(xx), + y(yy), + type(ttype) +{ + //if ( yy < xx ) + //throw "Point is below the diagonal"; + //if ( yy == xx and ttype != DiagramPoint::DIAG) + //throw "Point on the main diagonal must have DIAG type"; +} + +double DiagramPoint::getRealX() const +{ + if (isNormal()) + return x; + else + return 0.5 * ( x + y); +} + +double DiagramPoint::getRealY() const +{ + if (isNormal()) + return y; + else + return 0.5 * ( x + y); +} + +} // end of namespace geom_ws diff --git a/geom_matching/wasserstein/src/wasserstein.cpp b/geom_matching/wasserstein/src/wasserstein.cpp new file mode 100644 index 0000000..5761f11 --- /dev/null +++ b/geom_matching/wasserstein/src/wasserstein.cpp @@ -0,0 +1,121 @@ +/* + +Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov +All rights reserved. + +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. + +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +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 HOLDER 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. + +You are under no obligation whatsoever to provide any bug fixes, patches, or +upgrades to the features, functionality or performance of the source code +(Enhancements) to anyone; however, if you choose to make your Enhancements +available either publicly, or directly to copyright holder, +without imposing a separate written license agreement for such Enhancements, +then you hereby grant the following license: a non-exclusive, royalty-free +perpetual license to install, use, modify, prepare derivative works, incorporate +into other computer software, distribute, and sublicense such enhancements or +derivative works thereof, in binary and source code form. + + */ + +#include <assert.h> +#include <algorithm> +#include <functional> +#include <iterator> + +#include "def_debug.h" +#include "wasserstein.h" + +#ifdef GAUSS_SEIDEL_AUCTION +#include "auction_runner_gs.h" +#else +#include "auction_runner_jak.h" +#endif + +namespace geom_ws { + +double wassersteinDistVec(const std::vector<DiagramPoint>& A, + const std::vector<DiagramPoint>& B, + const double q, + const double delta, + const double _internal_p, + const double _initialEpsilon, + const double _epsFactor) +{ + if (q < 1) { + std::cerr << "Wasserstein distance not defined for q = " << q << ", must be >= 1" << std::endl; + throw "Bad q in Wasserstein"; + } + if (delta < 0.0) { + std::cerr << "Relative error " << delta << ", must be > 0" << std::endl; + throw "Bad delta in Wasserstein"; + } + if (_initialEpsilon < 0.0) { + std::cerr << "Initial epsilon = " << _initialEpsilon << ", must be non-negative" << std::endl; + throw "Bad delta in Wasserstein"; + } + if (_epsFactor < 0.0) { + std::cerr << "Epsilon factor = " << _epsFactor << ", must be non-negative" << std::endl; + throw "Bad delta in Wasserstein"; + } +#ifdef GAUSS_SEIDEL_AUCTION + AuctionRunnerGS auction(A, B, q, delta, _internal_p, _initialEpsilon, _epsFactor); +#else + AuctionRunnerJak auction(A, B, q, delta, _internal_p); +#endif + return auction.getWassersteinDistance(); +} + +bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<double, double>>& result) +{ + return readDiagramPointSet(fname.c_str(), result); +} + +bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double>>& result) +{ + size_t lineNumber { 0 }; + result.clear(); + std::ifstream f(fname); + if (!f.good()) { + std::cerr << "Cannot open file " << fname << std::endl; + return false; + } + std::string line; + while(std::getline(f, line)) { + lineNumber++; + // process comments: remove everything after hash + auto hashPos = line.find_first_of("#", 0); + if( std::string::npos != hashPos) { + line = std::string(line.begin(), line.begin() + hashPos); + } + if (line.empty()) { + continue; + } + // trim whitespaces + auto whiteSpaceFront = std::find_if_not(line.begin(),line.end(),isspace); + auto whiteSpaceBack = std::find_if_not(line.rbegin(),line.rend(),isspace).base(); + if (whiteSpaceBack <= whiteSpaceFront) { + // line consists of spaces only - move to the next line + continue; + } + line = std::string(whiteSpaceFront,whiteSpaceBack); + double x, y; + std::istringstream iss(line); + if (not(iss >> x >> y)) { + std::cerr << "Error in file " << fname << ", line number " << lineNumber << ": cannot parse \"" << line << "\"" << std::endl; + return false; + } + result.push_back(std::make_pair(x,y)); + } + f.close(); + return true; +} + +} // end of namespace geom_ws |