diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..92399ef --- /dev/null +++ b/.clang-format @@ -0,0 +1,66 @@ +# Generated from CLion C/C++ Code Style settings +BasedOnStyle: LLVM +AccessModifierOffset: -4 +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: None +AlignOperands: Align +AllowAllArgumentsOnNextLine: false +AllowAllConstructorInitializersOnNextLine: false +AllowAllParametersOfDeclarationOnNextLine: false +AllowShortBlocksOnASingleLine: Always +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: Always +AllowShortLambdasOnASingleLine: All +AllowShortLoopsOnASingleLine: true +AlwaysBreakAfterReturnType: None +AlwaysBreakTemplateDeclarations: Yes +BreakBeforeBraces: Custom +BraceWrapping: + AfterCaseLabel: false + AfterClass: false + AfterControlStatement: Never + AfterEnum: false + AfterFunction: false + AfterNamespace: false + AfterUnion: false + BeforeCatch: false + BeforeElse: false + IndentBraces: false + SplitEmptyFunction: false + SplitEmptyRecord: true +BreakBeforeBinaryOperators: None +BreakBeforeTernaryOperators: true +BreakConstructorInitializers: BeforeColon +BreakInheritanceList: BeforeColon +ColumnLimit: 0 +CompactNamespaces: false +ContinuationIndentWidth: 8 +IndentCaseLabels: true +IndentPPDirectives: None +IndentWidth: 4 +KeepEmptyLinesAtTheStartOfBlocks: true +MaxEmptyLinesToKeep: 2 +NamespaceIndentation: All +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PointerAlignment: Right +ReflowComments: false +SpaceAfterCStyleCast: true +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: false +SpaceBeforeAssignmentOperators: true +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeParens: ControlStatements +SpaceBeforeRangeBasedForLoopColon: true +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 0 +SpacesInAngles: false +SpacesInCStyleCastParentheses: false +SpacesInContainerLiterals: false +SpacesInParentheses: false +SpacesInSquareBrackets: false +TabWidth: 4 +UseTab: Never diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..1dc4554 --- /dev/null +++ b/.gitignore @@ -0,0 +1,9 @@ +.idea +cmake-build* +build +.DS_Store +data +data_bak +benchmark +__pycache__ + diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..0d24cef --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,81 @@ +cmake_minimum_required(VERSION 3.12) +project(GeoGraph) + +# default to Release building +if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + message("Defaulting CMAKE_BUILD_TYPE to Release") + set(CMAKE_BUILD_TYPE "RelWithDebInfo" CACHE STRING "Build type") +endif () + +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/") +include_directories("${PROJECT_SOURCE_DIR}/include/") + +# set default for MARCH +if (NOT MARCH) + set(MARCH native CACHE STRING "ARCH to use") +endif () + +set(CMAKE_CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSTIONS ON) + +set(CMAKE_C_STANDARD 11) +set(CMAKE_C_STANDARD_REQUIRED ON) +set(CMAKE_C_EXTENSTIONS ON) + +# setup modern c++ flags +string(REPLACE "-O2" "-O3" CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO}") +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=${MARCH} -pedantic -Wall -Wextra") +set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -march=${MARCH} -pedantic -Wall -Wextra") + +add_executable(YaoGraph YaoGraph.cpp) +add_executable(Generator Generator.cpp) + +option(WITH_CGAL "Build with CGAL for kernels and algorithms for comparision" OFF) +if (WITH_CGAL) + add_compile_definitions(WITH_CGAL) + find_package(CGAL) + include_directories(SYSTEM ${CGAL_INCLUDE_DIRS} ${INCLUDE_DIRECTORIES}) + target_link_libraries(YaoGraph ${CGAL_LIBRARIES}) +endif () + +option(WITH_CAIRO "Build with Cairo for drawing" OFF) +if (WITH_CAIRO) + add_compile_definitions(WITH_CAIRO) + find_package(Cairomm REQUIRED) + include_directories(SYSTEM "/usr/local/include" ${INCLUDE_DIRECTORIES}) + include_directories(SYSTEM ${Cairomm_INCLUDE_DIRS} ${INCLUDE_DIRECTORIES}) + target_link_libraries(YaoGraph ${Cairomm_LIBRARIES}) +endif () + +if (CMAKE_BUILD_TYPE MATCHES Debug) + message("debug mode - enable logging") + add_compile_definitions(LOG_DEBUG) +endif () + +option(WITH_VTUNE "Build for Profiling" OFF) +if (WITH_VTUNE) + add_compile_definitions(VTUNE) +endif () + +option(WITH_STATS "Collect execution statistics" OFF) +if (WITH_STATS) + add_compile_definitions(WITH_STATS) +endif () + +option(WITH_TESTS "Compile Tets" OFF) +if (WITH_TESTS) + find_package(GTest) + if (GTEST_FOUND) + add_subdirectory(test) + endif () +endif () + +set(ROAD_DIR "${PROJECT_SOURCE_DIR}/data" CACHE PATH "Path to road networks") +add_compile_definitions(ROAD_DATA="${ROAD_DIR}") + +set(DATA_DIR "${PROJECT_SOURCE_DIR}/data" CACHE PATH "Path to graph data") +add_compile_definitions(DATA_DIR="${DATA_DIR}") + +set(STAR_DIR "${PROJECT_SOURCE_DIR}/data" CACHE PATH "Path to gaia data") +add_compile_definitions(STAR_DATA="${STAR_DIR}") \ No newline at end of file diff --git a/Generator.cpp b/Generator.cpp new file mode 100644 index 0000000..49cf31d --- /dev/null +++ b/Generator.cpp @@ -0,0 +1,116 @@ +#include +#include +#include + +#include "contrib/popl.hpp" + +#include "Utils/ASSERT.hpp" + +#include "Generators/Generators.hpp" +#include "Types.hpp" + +// constants +constexpr tBox BOUNDS{{0, 0}, + {1, 1}}; +constexpr tIndex minN = 1e3; +constexpr tIndex maxN = 1e7; +constexpr tDim Cones = 6; +constexpr tIndex cellOcc = 1e2; +constexpr tDim RepsPerI = 3; +constexpr tDim RepsPerN = 3; + +const tIndex Seeds[] = {8158, 14030, 18545, 20099, 24065, 35700, 37197, 38132, 59135, 60315}; +const char Dists[] = {'u', 'g', 'd', 'r', 's', 'c', 'b'}; + +void generate(GeneratorBase &gen, const tIndex &n, const tBox &iBounds) { + + std::cout << "Generating " << gen.name() << " with target " << n << " points and seed " << gen.seed() << ": " << std::flush; + + auto [points, oBounds] = gen.generate(n, iBounds); + ASSERT(points.size() >= n); + + std::stringstream dir; + dir << DATA_DIR << "/" << gen.name(); + + std::filesystem::path p(dir.str()); + std::filesystem::create_directories(p); + + std::ofstream file(p.append("points_" + gen.name() + "_" + std::to_string(n) + "_" + std::to_string(gen.seed()) + ".csv"), std::ios::out | std::ios::trunc); + + // header + file << "# n " << points.size() << std::endl; + file << "# b " << oBounds.low[X] << " " << oBounds.low[Y] << " " << oBounds.high[X] << " " << oBounds.high[Y] << std::endl; + + for (auto &p : points) { + file << p[X] << " " << p[Y] << std::endl; + } + + std::cout << points.size() << " actually generated" << std::endl; +} + +int main(int argc, char **argv) { + + popl::OptionParser op("Point Generator"); + auto sHelp = op.add("h", "help", "produce help message"); + + // generate points + auto sAllDists = op.add("a", "all", "generate all available distributions"); + auto oN = op.add>("n", "n", "number of points to generate"); + auto oMinN = op.add>("", "minN", "minimum number of points to generate", minN); + auto oMaxN = op.add>("", "maxN", "maxium number of points to generate", maxN); + auto oSeed = op.add>("s", "seed", "seed for RNG", Seeds[0]); + auto oInst = op.add>("i", "inst", "number of instances"); + + op.parse(argc, argv); + + if (sHelp->is_set()) { + std::cout << op << "\n"; + std::cout << "list of distributions to generate [_u_ni, _g_aussian, gri_d_, _r_oad, _s_tar, _c_ircle, _b_ubbles]" << std::endl; + + return 0; + } + + std::vector seeds; + if (oInst->is_set()) { + // if oInst is set, oSeed is ignored + for (uint i = 0; i < oInst->value(); ++i) { + seeds.push_back(Seeds[i]); + } + } else { + seeds.push_back(oSeed->value()); + } + + std::vector dists; + if (sAllDists->is_set()) { + for (auto d : Dists) { + dists.push_back(d); + } + } else if (op.non_option_args().size()) { + for (uint i = 0; i < op.non_option_args().size(); ++i) { + dists.push_back(op.non_option_args()[i][0]); + } + } else { + std::cout << "Distribution must be specified" << std::endl; + } + + if (oN->is_set()) { + oMinN->set_value(oN->value()); + oMaxN->set_value(oN->value()); + } + + for (auto g : dists) { + auto gen = getGen(g, 0); + + if (!gen) { + std::cout << "Invalid distribution specified: " << g << std::endl; + continue; + } + + for (auto s : seeds) { + for (tIndex nPoints = oMinN->value(); nPoints <= oMaxN->value(); nPoints += 3 * pow(10, floor(log10(nPoints)))) { + gen->setSeed(s); + generate(*gen, nPoints, BOUNDS); + } + } + } +} diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..280050d --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2023 Daniel Funke + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..0e827d0 --- /dev/null +++ b/README.md @@ -0,0 +1,143 @@ +# Yao Graph Generator + +This tool generates the Yao-Graph of a given point set. +The Yao graph connects each point to its nearest neighbor in each of _k_ cones drawn around it. +This implementation works for two-dimensional inputs. + +## Compilation + +YaoGraph requires a modern C++ compiler with C++20 support. +To use exact predicates and constructions the CGAL library is required. + +```shell +mkdir build && cd build +cmake ../ +make +``` + +To use CGAL's kernels for predicates and constructions turn on the `WITH_CGAL` CMake option + +```shell +mkdir build && cd build +cmake -DWITH_CMAKE=ON ../ +make +``` + +Further CMake options are: +- `WITH_CAIRO` to draw images of the generated Yao graphs +- `WITH_STATS` to collect statistics during the sweepline execution (performance penalty) +- `WITH_TESTS` compile gtest-Suite +- `ROAD_DIR` directory containing DIMACS road networks (see below) +- `STAR_DIR` directory containing Gaia star catalogue (see below) + +## Usage + +### YaoGraph + +YaoGraph can either randomly generate points or use an input file with points + +Usage +```shell +./YaoGraph --help + +Yao Graph generator: + -h, --help produce help message + -b, --benchmark arg (=1) run algorithm i times + -d, --dist arg (=u) point distribution [_u_ni, _g_aussian, gri_d_, _r_oad, _s_tar, _c_ircle, _b_ubbles] + -n, --n arg number of points to generate + -s, --seed arg (=8158) seed for RNG + -f, --infile arg file with points + -a, --alg arg (=s) algorithm to use [_s_weepline, _g_rid, _n_aive, _c_gal] + -c, --kernel arg (=i) kernel to use [_i_nexact, CGALExact_p_redicatesInexactConstructions, CGALExactpredicatesInexact_c_onstructions] + -k, --cones arg (=6) number of cones, default 6 + -g, --cellOcc arg (=100) number of points per cell (grid algorithm) + -o, --outfile arg file for graph output + -p, --stdout write graph to stdout + -v, --verify verify graph (-v grid-based, -v -v naive yao (slow)) + --stats collect execution statistics (requires WITH_STATS compile flag) + +``` + +Input files have the following format +```shell +# n 1000 +# b 0 0 1 1 +0.227937 0.682662 +0.396135 0.946997 +0.844108 0.681706 +0.728227 0.537687 +``` + +The first line specifies the number of points in the file, +the second line defines the bounding box of the contained points. + +### Generator + +The generator application can be used to randomly generate input points and store the resulting points in files. +Especially the road and star distribution require some time to generate (see Data) + +```shell +./Generator --help + +Point Generator: + -h, --help produce help message + -a, --all generate all available distributions + -n, --n arg number of points to generate + --minN arg (=1000) minimum number of points to generate + --maxN arg (=10000000) maxium number of points to generate + -s, --seed arg (=8158) seed for RNG + -i, --inst arg number of instances + +list of distributions to generate [_u_ni, _g_aussian, gri_d_, _r_oad, _s_tar, _c_ircle, _b_ubbles] +``` + +## Benchmarks + +### Data + +All synthetic distributions can be generated without further inputs. +- For road networks you need to download the desired network from http://www.diag.uniroma1.it//challenge9/download.shtml +- The Gaia 2 star data can be downloaded from https://www.cosmos.esa.int/web/gaia/data-release-2 + +Benchmarks can then be executed with the Python script `scripts/runBenchmarks.py`: + +```shell +usage: runBenchmarks.py [-h] [-d {u,g,d,r,s,c,b} [{u,g,d,r,s,c,b} ...]] [--nMin NMIN] [--nMax NMAX] + [-s SEED [SEED ...]] [-i ITS] [-a {s,g,n,c} [{s,g,n,c} ...]] [-c {i,c,p} [{i,c,p} ...]] + [-k CONES] [--conesMin CONESMIN] [--conesMax CONESMAX] [-g CELLOCC] [-t TIMELIMIT] + [--stats] + +options: + -h, --help show this help message and exit + -d {u,g,d,r,s,c,b} [{u,g,d,r,s,c,b} ...], --dist {u,g,d,r,s,c,b} [{u,g,d,r,s,c,b} ...] + point distribution(s) to benchmark[_u_ni, _g_aussian, gri_d_, _r_oad, _s_tar, _c_ircle, + _b_ubbles] + --nMin NMIN minimum number of points + --nMax NMAX maximum number of points + -s SEED [SEED ...], --seed SEED [SEED ...] + seed(s) for RNG + -i ITS, --its ITS iterations per seed + -a {s,g,n,c} [{s,g,n,c} ...], --alg {s,g,n,c} [{s,g,n,c} ...] + algorithm(s) to use [_s_weepline, _g_rid, _n_aive, _c_gal] + -c {i,c,p} [{i,c,p} ...], --kernel {i,c,p} [{i,c,p} ...] + kernel(s) to use [_i_nexact, CGALExact_p_redicatesInexactConstructions, + CGALExactpredicatesInexact_c_onstructions] + -k CONES, --cones CONES + number of cones, default 6 + --conesMin CONESMIN minimum number of cones + --conesMax CONESMAX maximum number of cones + -g CELLOCC, --cellOcc CELLOCC + number of points per cell (grid algorithm) + -t TIMELIMIT, --timelimit TIMELIMIT + time limit for algorithm + --stats collect statistics + +``` + +### Plotting + +Plots can be generated with `benchmark/plot.py`. + + + + diff --git a/YaoGraph.cpp b/YaoGraph.cpp new file mode 100644 index 0000000..f6db97a --- /dev/null +++ b/YaoGraph.cpp @@ -0,0 +1,360 @@ +#include +#include + +#include "contrib/popl.hpp" + +#include "Utils/ASSERT.hpp" + +#include "Generators/Generators.hpp" +#include "Types.hpp" +#include "Utils/Timer.hpp" + +#include "Algorithms/GridYao.hpp" +#include "Algorithms/NaiveYao.hpp" +#include "Algorithms/SweepLine.hpp" + +#include "Kernels/InexactKernel.hpp" + +#ifdef WITH_CGAL +#include "Algorithms/CGAL_Yao.hpp" +#include "Kernels/CGALKernel.hpp" +#endif + +#ifdef WITH_CAIRO +#include "Painter.hpp" +#endif + + +// constants +constexpr tBox BOUNDS{{0, 0}, + {1, 1}}; +constexpr tIndex cellOcc = 1e2; +tDim Reps = 1; + +const tIndex Seeds[] = {8158, 14030, 18545, 20099, 24065, 35700, 37197, 38132, 59135, 60315}; + +std::tuple readPointsFile(const std::string &fileName) { + std::ifstream file(fileName, std::ios::in); + std::string line; + + tPoints points; + tIFloatVector minPoint = {std::numeric_limits::max(), std::numeric_limits::max()}; + tIFloatVector maxPoint = {std::numeric_limits::min(), std::numeric_limits::min()}; + + while (std::getline(file, line)) { + + if (line.starts_with("# n")) { + // number of points line + // # n XXX + auto lastSpace = line.find_last_of(' '); + tIndex pointsInFile = std::stoul(line.substr(lastSpace + 1)); + points.reserve(pointsInFile); + + continue; + } + + if (line.starts_with("# b")) { + // bounding box line + // # b minX minY maxX maxY + tIFloat minX, minY, maxX, maxY; + std::stringstream iss(line.substr(3)); + + if (!(iss >> minX >> minY >> maxX >> maxY)) { continue; }// error in file + + minPoint[X] = minX; + minPoint[Y] = minY; + maxPoint[X] = maxX; + maxPoint[Y] = maxY; + + continue; + } + + // vertex line + // X Y + std::istringstream iss(line); + tIFloat x, y; + if (!(iss >> x >> y)) { continue; }// error in file + + points.push_back({x, y}); + + if (x < minPoint[X]) minPoint[X] = x; + if (y < minPoint[Y]) minPoint[Y] = y; + if (x > maxPoint[X]) maxPoint[X] = x; + if (y > maxPoint[Y]) maxPoint[Y] = y; + } + + return std::make_tuple(points, tBox{minPoint, maxPoint}); +} + +bool gBenchmarkMode = false; + +template typename Algorithm, typename Kernel, typename... Args> +auto runAlg(const tDim &K, const tPoints &points, const tBox &bounds, Args... args) { + + tYaoGraph graph(points.size(), K); + + for (tDim rpi = 0; rpi < Reps; ++rpi) { + + if (!gBenchmarkMode) { + std::cout << "Generating Yao graph of " << points.size() << " points with " + << Algorithm::name() << " with " << KernelName::name() << " kernel" << std::endl; + } + + auto result = Timer>::time(K, points, bounds, args...); + + if (gBenchmarkMode) { + std::cout << "RESULT" + << " algorithm=" << Algorithm::name() + << " kernel=" << KernelName::name() + << " k=" << K + << " dist=" << points.getDistName() + << " n=" << points.size() + << " seed=" << points.getSeed() + << " rep=" << rpi + << " t=" << std::get<0>(result) + << std::endl; + } else { + std::cout << "Time: " << std::get<0>(result) << "ms" << std::endl; + } + + graph = std::get<1>(result); + } + + return graph; +} + +int main(int argc, char **argv) { + popl::OptionParser op("Yao Graph generator"); + auto sHelp = op.add("h", "help", "produce help message"); + + // benchmark + auto oBenchmark = op.add>("b", "benchmark", "run algorithm i times", 1); + + // generate points + auto oDist = op.add>("d", "dist", "point distribution [_u_ni, _g_aussian, gri_d_, _r_oad, _s_tar, _c_ircle, _b_ubbles]", 'u'); + auto oN = op.add>("n", "n", "number of points to generate"); + auto oSeed = op.add>("s", "seed", "seed for RNG", Seeds[0]); + + // points file + auto oPointsFile = op.add>("f", "infile", "file with points"); + + // algorithm to use + auto oAlg = op.add>("a", "alg", "algorithm to use [_s_weepline, _g_rid, _n_aive, _c_gal]", 's'); + auto oKern = op.add>("c", "kernel", "kernel to use [_i_nexact, CGALExact_p_redicatesInexactConstructions, CGALExactpredicatesInexact_c_onstructions]", 'i'); + auto oK = op.add>("k", "cones", "number of cones, default 6", 6); + auto oCellOcc = op.add>("g", "cellOcc", "number of points per cell (grid algorithm)", cellOcc); + + // output + auto oOutFile = op.add>("o", "outfile", "file for graph output"); + auto sStdOut = op.add("p", "stdout", "write graph to stdout"); +#ifdef WITH_CAIRO + auto oImageOut = op.add>("i", "image", "write image of Yao graph to specified file"); +#endif + + //verify + auto sVerify = op.add("v", "verify", "verify graph (-v grid-based, -v -v naive yao (slow))"); + auto sStats = op.add("", "stats", "collect execution statistics (requires WITH_STATS compile flag)"); + + op.parse(argc, argv); + + if (sHelp->is_set()) { + std::cout << op << "\n"; + + return 0; + } + + if (oBenchmark->is_set()) { + gBenchmarkMode = true; + Reps = oBenchmark->value(); + } + + if (sStats->is_set()) { +#ifndef WITH_STATS + std::cout << "Please compile with WITH_STATS option enabled" << std::endl; + return 4; +#endif + } + + // get points + tPoints points; + tBox bounds; + if (oN->is_set()) { + auto gen = getGen(oDist->value(), oSeed->value()); + std::tie(points, bounds) = gen->generate(oN->value(), BOUNDS); + } else if (oPointsFile->is_set()) { + std::tie(points, bounds) = readPointsFile(oPointsFile->value()); + + if (oDist->is_set()) { + points.setDistName(getGenName(oDist->value())); + } + + if (oSeed->is_set()) { + points.setSeed(oSeed->value()); + } + } else { + std::cout << "Either specify input point file or point generation option" << std::endl; + return 0; + } + +#ifdef WITH_CAIRO + if (oImageOut->is_set()) { + points = GeneratorBase::rescalePoints(points, bounds, BOUNDS); + bounds = BOUNDS; + } +#endif + + // run algorithm + tYaoGraph graph(points.size(), oK->value()); + switch (oAlg->value()) { + case 's': + switch (oKern->value()) { + case 'i': + graph = runAlg(oK->value(), points, bounds); + break; +#ifdef WITH_CGAL + case 'p': + graph = runAlg>(oK->value(), points, bounds); + break; + case 'c': + graph = runAlg>(oK->value(), points, bounds); + break; +#endif + } + + break; + case 'g': + switch (oKern->value()) { + case 'i': + graph = runAlg(oK->value(), points, bounds, oCellOcc->value()); + break; +#ifdef WITH_CGAL + case 'p': + graph = runAlg>(oK->value(), points, bounds, oCellOcc->value()); + break; + case 'c': + graph = runAlg>(oK->value(), points, bounds, oCellOcc->value()); + break; +#endif + } + + break; + case 'n': + switch (oKern->value()) { + case 'i': + graph = runAlg(oK->value(), points, bounds); + break; +#ifdef WITH_CGAL + case 'p': + graph = runAlg>(oK->value(), points, bounds); + break; + case 'c': + graph = runAlg>(oK->value(), points, bounds); + break; +#endif + } + break; +#ifdef WITH_CGAL + case 'c': + switch (oKern->value()) { + case 'p': + runAlg(oK->value(), points, bounds); + break; + case 'c': + runAlg(oK->value(), points, bounds); + break; + } +#endif + } + + if (oOutFile->is_set()) { + std::ofstream file(oOutFile->value(), std::ios::out | std::ios::trunc); + file << graph << std::endl; + } + + if (sStdOut->is_set()) { + std::cout << graph << std::endl; + } + +#ifdef WITH_CAIRO + if (oImageOut->is_set()) { + Painter painter(bounds); + painter.draw(points, false); + painter.save(oImageOut->value() + ".points"); + painter.draw(graph, points); + painter.save(oImageOut->value() + ".graph"); + } +#endif + + if (sVerify->is_set()) { + tYaoGraph exp(points.size(), oK->value()); + + if (sVerify->count() == 1) { + switch (oKern->value()) { + case 'i': + exp = runAlg(oK->value(), points, bounds, oCellOcc->value()); + break; +#ifdef WITH_CGAL + case 'p': + exp = runAlg>(oK->value(), points, bounds, oCellOcc->value()); + break; + case 'c': + exp = runAlg>(oK->value(), points, bounds, oCellOcc->value()); + break; +#endif + } + } else if (sVerify->count() == 2) { + switch (oKern->value()) { + case 'i': + exp = runAlg(oK->value(), points, bounds); + break; +#ifdef WITH_CGAL + case 'p': + exp = runAlg>(oK->value(), points, bounds); + break; + case 'c': + exp = runAlg>(oK->value(), points, bounds); + break; +#endif + } + } + + auto [valid, invalidVertices] = checkGraph(graph, exp); + +#ifdef WITH_CAIRO + if (!valid) { + Painter basePainter(bounds); + basePainter.draw(points, true); + basePainter.draw(exp, points); + + for (auto idx : invalidVertices) { + + Painter painter(basePainter); + + for (tDim k = 0; k < oK->value(); ++k) { + if (exp[idx].neighbor[k] != graph[idx].neighbor[k]) { + + if (exp[idx].neighbor[k] != INF_IDX) { + painter.setColor(0, 1, 0); + painter.draw(points[exp[idx].neighbor[k]], exp[idx].neighbor[k], true); + } + + if (graph[idx].neighbor[k] != INF_IDX) { + painter.setColor(1, 0, 0); + painter.draw(points[graph[idx].neighbor[k]], graph[idx].neighbor[k], true); + } + } + } + + painter.setColor(1, 0, 0); + painter.draw(idx, graph[idx], points); + painter.setColor(0, 0, 1); + painter.drawCones(points[idx], oK->value()); + + painter.save("invalidVertices_" + std::to_string(idx)); + } + } +#endif// ifdef WITH_CAIRO + + return !valid; + } +} diff --git a/benchmark/plot.py b/benchmark/plot.py new file mode 100644 index 0000000..6b0da58 --- /dev/null +++ b/benchmark/plot.py @@ -0,0 +1,531 @@ +import os +import re + +import numpy as np +import pandas as pd +import matplotlib as mpl +import matplotlib.pylab as plt +import seaborn as sns + +reRuntimeFilename = re.compile('benchmark_([a-zA-Z0-9]+)_([a-zA-Z]+).csv') +reConesFilename = re.compile('benchmark_cones_([a-zA-Z0-9]+)_([a-zA-Z]+).csv') +isRunningInPyCharm = 'PYCHARM_HOSTED' in os.environ + +BASE_DIR = '/home/funke/devel/geograph/benchmark/' +DATA_DIR = os.path.join(BASE_DIR, 'data') +PLOT_DIR = os.path.join(BASE_DIR, 'plots') +EXT = 'png' + + +def saveFig(f: str): + plt.tight_layout(pad=0.3) + plt.savefig(os.path.join(PLOT_DIR, f)) + if isRunningInPyCharm: + plt.show() + plt.close() + + +BASE_SIZE = 20 +LARGE_SIZE = 34 + +# sns.set(style='whitegrid', font_scale=1, rc={'text.usetex': True, +# 'text.latex.preamble': r'\usepackage{nicefrac} \usepackage{amsmath} \usepackage{lmodern}', +# "font.size": BASE_SIZE, +# "font.family": 'serif' +# "axes.titlesize": BASE_SIZE, +# "axes.labelsize": LARGE_SIZE, +# "xtick.labelsize" : LARGE_SIZE, "ytick.labelsize" : LARGE_SIZE, +# "legend.fontsize" : BASE_SIZE, +# "legend.title_fontsize" : BASE_SIZE}) + +sns.set(style='whitegrid', font_scale=2, rc={'text.usetex': True, + 'text.latex.preamble': r'\usepackage{nicefrac} \usepackage{amsmath} \usepackage{lmodern}', + "font.family": 'serif', + "lines.markersize": 10}) + +LegendLabels = { + 'gaussian': r'Gaussian', 'uni': r'Uniform', 'grid': r'Grid', 'road': r'Road', 'stars': r'Stars', + 'circle': r'Circle', 'bubbles': r'Bubbles', + 'NaiveYao': r'Naive Yao', 'Sweepline': r'Sweepline', 'GridYao': r'Grid Yao', 'CGALYao': r'CGAL Yao', + 'InexactKernel': r'Inexact Kernel', 'CGALExactPredInexactCon': r'CGAL EPIC', 'CGALExactPredExactCon': r'CGAL EPEC', + 'kernel': 'Kernel', 'algorithm': r'Algorithm', 'dist': r'Distribution', +} + + +def setLegend(o): + if isinstance(o, mpl.axes.Axes): + L = o.get_legend() + + if isinstance(o, mpl.legend.Legend): + L = o + + if (L.get_title() and L.get_title().get_text() != ''): + L.set_title(LegendLabels[L.get_title().get_text()]) + L.set_title(None) + + for i in range(len(L.get_texts())): + L.get_texts()[i].set_text(LegendLabels[L.get_texts()[i].get_text()]) + + +def unsetLegend(ax: mpl.axes.Axes): + ax.get_legend().remove() + + +###################################################################################################### +# %% Runtime plots - optics + +snsPalette = sns.color_palette(n_colors=3) +KernelPalette = {'InexactKernel': snsPalette[0], 'CGALExactPredInexactCon': snsPalette[1], + 'CGALExactPredExactCon': snsPalette[2]} +KernelDash = {'InexactKernel': (0, 5, 10), 'CGALExactPredInexactCon': (None, None), 'CGALExactPredExactCon': (0, 1, 1)} +KernelMarkers = {'InexactKernel': 's', 'CGALExactPredInexactCon': 'D', 'CGALExactPredExactCon': 'o'} + +snsPalette = sns.color_palette(n_colors=7) +DistPalette = {'gaussian': snsPalette[0], 'uni': snsPalette[1], 'grid': snsPalette[2], 'road': snsPalette[3], + 'stars': snsPalette[4], 'circle': snsPalette[5], 'bubbles': snsPalette[6]} + +snsPalette = sns.color_palette(n_colors=4) +AlgorithmDash = {'NaiveYao': (0, 5, 10), 'Sweepline': (None, None), 'GridYao': (0, 1, 1), 'CGALYao': (0, 3, 5, 1, 5)} +AlgorithmMarkers = {'NaiveYao': 's', 'Sweepline': 'D', 'GridYao': 'o', 'CGALYao': '^'} +AlgorithmPalette = {'NaiveYao': snsPalette[0], 'Sweepline': snsPalette[1], 'GridYao': snsPalette[2], + 'CGALYao': snsPalette[3]} + +figTextX = .94 +figTextY = .91 + +###################################################################################################### +# %% Priority Queue plots +pltGroup = 10 +pqDataFile = os.path.join(DATA_DIR, 'pq.csv') +if os.path.exists(pqDataFile): + with open(pqDataFile, 'r') as file: + header = file.readline() + + header = header[1:].strip().split() + + gPQ = pd.read_csv(pqDataFile, sep=' ', comment='#', names=header) + lPQ = gPQ[gPQ['k'] == 0] + + plt.stackplot(lPQ['step'], lPQ['ipPro'], lPQ['isPro'], lPQ['delPro'], + labels=['input points', 'intersections', 'deletions']) + plt.legend(loc='upper left') + plt.xlabel('algorithm step') + plt.ylabel('events processed') + saveFig('%i_pq_events_processed.%s' % (pltGroup, EXT)) + + plt.stackplot(lPQ['step'], lPQ['ipQ'], lPQ['isQ'], lPQ['delQ'], + labels=['input points', 'intersections', 'deletions']) + # setLegend() + plt.xlabel('algorithm step') + plt.ylabel('events in PQ') + saveFig('%i_pq_events_inQueue.%s' % (pltGroup, EXT)) + + plt.stackplot(lPQ['step'], lPQ['isQ'], lPQ['delQ'], labels=['intersections', 'deletions']) + # setLegend() + plt.xlabel('algorithm step') + plt.ylabel('events in PQ') + saveFig('%i_pq_events_inQueue_noInput.%s' % (pltGroup, EXT)) + + plt.stackplot(lPQ['step'], lPQ['slSize'], labels=['rays']) + # setLegend() + plt.xlabel('algorithm step') + plt.ylabel('rays in SL') + saveFig('%i_rays_inSL.%s' % (pltGroup, EXT)) + +###################################################################################################### +# %% Priority Queue aggregate plots +pltGroup = 20 +statsDataFile = os.path.join(DATA_DIR, 'stats.csv') +if os.path.exists(statsDataFile): + with open(statsDataFile, 'r') as file: + header = file.readline() + + header = header[1:].strip().split() + + gStats = pd.read_csv(statsDataFile, sep=' ', comment='#', names=header) + gStats['maxIsDelQ'] = gStats['maxIsQ'] + gStats['maxDelQ'] + gStats['stepsPN'] = gStats['steps'] / gStats['n'] + gStats['maxIsDelQSqrt'] = gStats['maxIsDelQ'] / np.sqrt(gStats['n']) + gStats['maxIsDelQPN'] = gStats['maxIsDelQ'] / gStats['n'] + gStats['maxSlSizeSqrt'] = gStats['maxSlSize'] / np.sqrt(gStats['n']) + gStats['maxSlSizePN'] = gStats['maxSlSize'] / (gStats['n'] * np.log(gStats['n'])) + + Ns = pd.Series(gStats[gStats['dist'] == 'uni']['n'].unique()) + + + def custom_round(x): + return Ns[Ns.sub(x).abs().idxmin()] + + + gStats['dN'] = gStats['n'].apply(lambda x: custom_round(x)) + + # lStats = gStats[(gStats['dist'] != 'circle') & (gStats['dist'] != 'road') & (gStats['dist'] != 'stars')] + lStats = gStats[gStats['dist'] != 'bubbles'] + + ax = sns.lineplot(data=lStats, x='dN', y='stepsPN', hue='dist', style='dist', markers=True, + palette=DistPalette) + # plt.legend(loc='lower center', ncol=2, bbox_to_anchor=(0.5, 0.8)) + unsetLegend(ax) + plt.xscale('log') + plt.xlabel(r'$n$') + plt.ylabel(r'$\nicefrac{\text{events}}{n}$') + saveFig('%i_pq_events_processed.%s' % (pltGroup, EXT)) + + ax = sns.lineplot(data=lStats, x='dN', y='maxIsDelQSqrt', hue='dist', style='dist', markers=True, + palette=DistPalette) + # plt.legend(loc='lower center', ncol=2, bbox_to_anchor=(0.5, 0.8)) + plt.legend(ncol=2) + setLegend(ax) + plt.xscale('log') + plt.xlabel(r'$n$') + plt.ylabel(r'$\nicefrac{\text{events in PQ}}{\sqrt{n}}$') + saveFig('%i_pq_events_inQueue_noInput.%s' % (pltGroup, EXT)) + + ax = sns.lineplot(data=lStats, x='dN', y='maxSlSizePN', hue='dist', style='dist', markers=True, + palette=DistPalette) + # plt.legend(loc='upper left', ncol=2, bbox_to_anchor=(0, 1.15)) + # plt.legend(ncol=2) + unsetLegend(ax) + plt.xscale('log') + # plt.ylim((3, 10)) + plt.yscale('log', base=10) + plt.xlabel('$n$') + plt.ylabel(r'$\nicefrac{\text{rays in SL}}{n}$') + saveFig('%i_rays_inSL.%s' % (pltGroup, EXT)) + +###################################################################################################### +# %% Sweepline tree update operations plots +pltGroup = 30 +upDataFile = os.path.join(DATA_DIR, 'update.csv') +if os.path.exists(upDataFile): + with open(upDataFile, 'r') as file: + header = file.readline() + + header = header[1:].strip().split() + + gUp = pd.read_csv(upDataFile, sep=' ', comment='#', names=header) + sns.displot(gUp, x='change', col='op', discrete=True, stat='probability', common_norm=False, shrink=.8) + plt.xticks([0, 1]) + + saveFig('%i_tree_updates.%s' % (pltGroup, EXT)) + +###################################################################################################### +# %% Runtime plots - read data + +llData = [] +for f in os.listdir(DATA_DIR): + match = reRuntimeFilename.match(f) + if match: + # get header from first commented line + with open(os.path.join(DATA_DIR, f), 'r') as file: + header = file.readline() + + header = header[1:].strip().split() + + lData = pd.read_csv(os.path.join(DATA_DIR, f), sep=' ', comment='#', names=header) + lData['algorithm'] = match.group(1) + lData['kernel'] = match.group(2) + llData.append(lData) + +gData = pd.concat(llData) +gData.set_index(['algorithm', 'kernel', 'dist', 'n', 'seed', 'rep'], inplace=True) +gData.reset_index(inplace=True) + +TargetNs = pd.Series(gData[gData['dist'] == 'uni']['n'].unique()) + + +def custom_round(x): + return TargetNs[TargetNs.sub(x).abs().idxmin()] + + +gData['dN'] = gData['n'].apply(lambda x: custom_round(x)) + +gData['tpn'] = gData['t'] / gData['n'] +gData['speedup'] = 1 + +# calculate speedup over naive +for dist in gData['dist'].unique(): + for alg in gData['algorithm'].unique(): + for kernel in gData[gData['algorithm'] == alg]['kernel'].unique(): + maskNav = (gData['dist'] == dist) & (gData['algorithm'] == 'NaiveYao') & (gData['kernel'] == kernel) + + Ns = gData.loc[maskNav, 'n'].unique() + + maskAlg = (gData['dist'] == dist) & (gData['algorithm'] == alg) & (gData['kernel'] == kernel) & ( + gData['n'].isin(Ns)) + + # gData.loc[maskAlg, 'speedup'] = gData.loc[maskNav, 't'].array / gData.loc[maskAlg, 't'] + +TIMELIMIT = 1800000 # ms (30 min) +Ns = gData['n'].unique() +TLpN = TIMELIMIT / Ns + +###################################################################################################### +# %% Runtime plots - overview plot + +pltGroup = 40 +ax = sns.lineplot(data=gData, x='dN', y='tpn', style='algorithm', hue='kernel', dashes=AlgorithmDash, + markers=AlgorithmMarkers, + palette=KernelPalette) +w, h = plt.gcf().get_size_inches() +plt.gcf().set_size_inches(1.5 * w, h) +plt.legend(loc='upper left', bbox_to_anchor=(1.1, 1.1)) +setLegend(ax) +plt.xscale('log') +plt.yscale('log', base=2) +plt.xlabel(r'$n$') +plt.ylabel(r'$\nicefrac{t}{n}$ [ms]') +saveFig('%i_runtime.%s' % (pltGroup, EXT)) + +ax = sns.lineplot(data=gData, x='dN', y='speedup', style='kernel', hue='algorithm', dashes=KernelDash, + markers=KernelMarkers, + palette=AlgorithmPalette) +w, h = plt.gcf().get_size_inches() +plt.gcf().set_size_inches(1.5 * w, h) +plt.legend(loc='upper left', bbox_to_anchor=(1.1, 1.1)) +setLegend(ax) +plt.xscale('log') +plt.yscale('log', base=2) +plt.xlabel(r'$n$') +plt.ylabel(r'speedup') +saveFig('%i_speedup.%s' % (pltGroup, EXT)) + +###################################################################################################### +# %% Runtime plots - overview plot (no Naive/CGAL) +pltGroup = 45 +fData = gData[(gData['algorithm'] == 'Sweepline') ^ (gData['algorithm'] == 'GridYao')] +ax = sns.lineplot(data=fData, x='dN', y='tpn', style='algorithm', hue='kernel', dashes=AlgorithmDash, + markers=AlgorithmMarkers, + palette=KernelPalette) +unsetLegend(ax) +plt.xscale('log') +plt.yscale('log', base=2) +plt.xlabel(r'$n$') +plt.ylabel(r'$\nicefrac{t}{n}$ [ms]') +saveFig('%i_runtime.%s' % (pltGroup, EXT)) + +ax = sns.lineplot(data=fData, x='dN', y='speedup', style='kernel', hue='algorithm', dashes=KernelDash, + markers=KernelMarkers, + palette=AlgorithmPalette) +unsetLegend(ax) +plt.xscale('log') +plt.yscale('log', base=2) +plt.xlabel(r'$n$') +plt.ylabel(r'speedup') +saveFig('%i_speedup.%s' % (pltGroup, EXT)) + +###################################################################################################### +# %% Runtime plots - per kernel plots +pltGroup = 50 +for kernel in gData['kernel'].unique(): + # filter data by dist + fData = gData[(gData['kernel'] == kernel) & (gData['dist'] != 'bubbles')] + + plt.close() + + fDataNoCG = fData[~((fData['dist'] == 'circle') & (fData['algorithm'] == 'GridYao'))] + ax = sns.lineplot(data=fDataNoCG, x='dN', y='tpn', style='algorithm', hue='algorithm', + dashes=AlgorithmDash, + markers=AlgorithmMarkers, + palette=AlgorithmPalette, + markersize=8) + + fDataCG = fData[((fData['dist'] == 'circle') & (fData['algorithm'] == 'GridYao'))] + ax.scatter(fDataCG['dN'], fDataCG['tpn'], color='g', marker=7, + s=6 ** 2, zorder=10) + + ax.vlines(x=fDataCG['dN'].unique(), + ymin=fDataNoCG[fDataNoCG['algorithm'] == 'GridYao'].groupby('dN').agg('mean')['tpn'], + ymax=fDataCG.groupby('dN').agg('mean')['tpn'], + colors='g', linestyles='--') + + ax.legend(loc='upper left') + setLegend(ax) + plt.xscale('log') + plt.yscale('log', base=10) + + # plt.xlim(1e3, 1e5) + + l, h = plt.ylim() + plt.plot(Ns, TLpN, c='gray') + plt.ylim(l, h) + + plt.xlabel(r'$n$') + plt.ylabel(r'$\nicefrac{t}{n}$ [ms]') + # plt.text(figTextX, figTextY, 'Kernel: %s' % (LegendLabels[kernel]), + # horizontalalignment='right', + # verticalalignment='bottom', + # transform=plt.gcf().transFigure) + saveFig('%i_%s_runtime.%s' % (pltGroup, kernel, EXT)) + + fig, axs = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=True) + + # DISTS = sorted(fData['dist'].unique()) + DISTS = ['uni', 'gaussian', 'grid', 'road', 'stars', 'circle'] + + for i, dist in enumerate(DISTS): + sns.lineplot(data=fData[fData['dist'] == dist], x='dN', y='tpn', style='algorithm', hue='dist', dashes=AlgorithmDash, + markers=AlgorithmMarkers, + palette=DistPalette, + markersize=8, + ax=fig.axes[i]) + + fig.axes[i].set_title(LegendLabels[dist]) + fig.axes[i].set_xlabel(r'$n$') + fig.axes[i].set_ylabel(r'$\nicefrac{t}{n}$ [ms]') + fig.axes[i].set_xscale('log') + fig.axes[i].set_yscale('log', base=10) + + l, h = fig.axes[i].get_ylim() + fig.axes[i].plot(Ns, TLpN, c='gray') + fig.axes[i].set_ylim(l, h) + + unsetLegend(fig.axes[i]) + + w, h = fig.get_size_inches() + fig.set_size_inches(3 * w, 2.2 * h) + + h, l = fig.axes[0].get_legend_handles_labels() + + h = h[2:] + l = l[2:] + + L = plt.figlegend(h, l, loc='center right') + setLegend(L) + + + # plt.ylabel(r'$\nicefrac{t}{n}$ [ms]') + # plt.text(figTextX, figTextY, 'Kernel: %s' % (LegendLabels[kernel]), + # horizontalalignment='right', + # verticalalignment='bottom', + # transform=plt.gcf().transFigure) + saveFig('%i_%s_runtime_dist.%s' % (pltGroup, kernel, EXT)) + +###################################################################################################### +# %% Runtime plots - per kernel plots (no Naive/CGAL) +pltGroup = 55 + +fData = gData[(gData['algorithm'] == 'Sweepline') ^ (gData['algorithm'] == 'GridYao')] +for kernel in gData['kernel'].unique(): + # filter data by dist + ffData = fData[fData['kernel'] == kernel] + + sns.lineplot(data=ffData, x='dN', y='tpn', style='algorithm', hue='algorithm', dashes=AlgorithmDash, + markers=AlgorithmMarkers, + palette=AlgorithmPalette, + markersize=8) + plt.xscale('log') + plt.yscale('log', base=2) + plt.text(figTextX, figTextY, 'Kernel: %s' % (kernel), + horizontalalignment='right', + verticalalignment='bottom', + transform=plt.gcf().transFigure) + saveFig('%i_%s_runtime.%s' % (pltGroup, kernel, EXT)) + + sns.lineplot(data=ffData, x='dN', y='tpn', style='algorithm', hue='dist', dashes=AlgorithmDash, + markers=AlgorithmMarkers, + palette=DistPalette, + markersize=8) + plt.xscale('log') + plt.yscale('log', base=2) + plt.text(figTextX, figTextY, 'Kernel: %s' % (kernel), + horizontalalignment='right', + verticalalignment='bottom', + transform=plt.gcf().transFigure) + saveFig('%i_%s_runtime_dist.%s' % (pltGroup, kernel, EXT)) + +###################################################################################################### +# %% Runtime plots - per dist plots +pltGroup = 60 +for dist in gData['dist'].unique(): + + # filter data by dist + fData = gData[gData['dist'] == dist] + + sns.lineplot(data=fData, x='dN', y='tpn', style='algorithm', hue='kernel', dashes=AlgorithmDash, + markers=AlgorithmMarkers, + palette=KernelPalette) + plt.text(figTextX, figTextY, 'Distribution: %s' % (dist), + horizontalalignment='right', + verticalalignment='bottom', + transform=plt.gcf().transFigure) + saveFig('%i_%s_runtime.%s' % (pltGroup, dist, EXT)) + + for alg in fData['algorithm'].unique(): + sns.lineplot(data=fData[fData['algorithm'] == alg], x='dN', y='tpn', style='algorithm', hue='kernel', + dashes=AlgorithmDash, + markers=AlgorithmMarkers, palette=KernelPalette) + plt.text(figTextX, figTextY, 'Distribution: %s' % (dist), + horizontalalignment='right', + verticalalignment='bottom', + transform=plt.gcf().transFigure) + saveFig('%i_%s_runtime_alg_%s.%s' % (pltGroup, dist, alg, EXT)) + + for kernel in fData['kernel'].unique(): + sns.lineplot(data=fData[fData['kernel'] == kernel], x='dN', y='tpn', style='algorithm', hue='kernel', + dashes=AlgorithmDash, + markers=AlgorithmMarkers, palette=KernelPalette) + plt.text(figTextX, figTextY, 'Distribution: %s' % (dist), + horizontalalignment='right', + verticalalignment='bottom', + transform=plt.gcf().transFigure) + saveFig('%i_%s_runtime_kernel_%s.%s' % (pltGroup, dist, kernel, EXT)) + + sns.lineplot(data=fData[~fData['algorithm'].isin(['NaiveYao', 'CGALYao'])], x='dN', y='tpn', style='algorithm', + hue='kernel', dashes=AlgorithmDash, markers=AlgorithmMarkers, palette=KernelPalette) + plt.text(figTextX, figTextY, 'Distribution: %s' % (dist), + horizontalalignment='right', + verticalalignment='bottom', + transform=plt.gcf().transFigure) + saveFig('%i_%s_runtime_sweepline_grid.%s' % (pltGroup, dist, EXT)) + + sns.lineplot( + data=fData[ + ~fData['algorithm'].isin(['NaiveYao', 'CGALYao']) & ~fData['kernel'].isin(['CGALExactPredExactCon'])], + x='dN', y='tpn', style='algorithm', hue='kernel', dashes=AlgorithmDash, markers=AlgorithmMarkers, + palette=KernelPalette) + plt.text(figTextX, figTextY, 'Distribution: %s' % (dist), + horizontalalignment='right', + verticalalignment='bottom', + transform=plt.gcf().transFigure) + saveFig('%i_%s_runtime_sweepline_grid_noExactCon.%s' % (pltGroup, dist, EXT)) + +###################################################################################################### +# %% Runtime over cones plots - read data + +llData = [] +for f in os.listdir(DATA_DIR): + match = reConesFilename.match(f) + if match: + # get header from first commented line + with open(os.path.join(DATA_DIR, f), 'r') as file: + header = file.readline() + + header = header[1:].strip().split() + + lData = pd.read_csv(os.path.join(DATA_DIR, f), sep=' ', comment='#', names=header) + lData['algorithm'] = match.group(1) + lData['kernel'] = match.group(2) + llData.append(lData) + +gData = pd.concat(llData) +gData.set_index(['algorithm', 'kernel', 'dist', 'n', 'seed', 'rep'], inplace=True) +gData.reset_index(inplace=True) + +gData['tpn'] = gData['t'] / gData['n'] +gData['tpk'] = gData['t'] / gData['k'] + +###################################################################################################### +# %% Runtime over cones plots - overview plot + +pltGroup = 70 +ax = sns.lineplot(data=gData, x='k', y='tpk', style='algorithm', hue='kernel', dashes=AlgorithmDash, + markers=AlgorithmMarkers, + palette=KernelPalette) +plt.legend(loc='lower center', ncol=2, bbox_to_anchor=(.5, .8)) +setLegend(ax) +# plt.xscale('log') +plt.yscale('log', base=10) +plt.xlabel(r'$k$') +plt.ylabel(r'$\nicefrac{t}{k}$ [ms]') +saveFig('%i_runtime_cones.%s' % (pltGroup, EXT)) diff --git a/cmake/FindCairo.cmake b/cmake/FindCairo.cmake new file mode 100644 index 0000000..acde22e --- /dev/null +++ b/cmake/FindCairo.cmake @@ -0,0 +1,34 @@ +# - Try to find Cairo +# Once done, this will define +# +# Cairo_FOUND - system has Cairo +# Cairo_INCLUDE_DIRS - the Cairo include directories +# Cairo_LIBRARIES - link these to use Cairo + +include(LibFindMacros) + +# Dependencies +libfind_package(Cairo Freetype) + +# Use pkg-config to get hints about paths +libfind_pkg_check_modules(Cairo_PKGCONF cairo) + +# Include dir +find_path(Cairo_INCLUDE_DIR + NAMES cairo.h + PATHS ${Cairo_PKGCONF_INCLUDE_DIRS} + PATH_SUFFIXES cairo +) + +# Finally the library itself +find_library(Cairo_LIBRARY + NAMES cairo + PATHS ${Cairo_PKGCONF_LIBRARY_DIRS} +) + +# Set the include dir variables and the libraries and let libfind_process do the rest. +# NOTE: Singular variables for this library, plural for libraries this this lib depends on. +set(Cairo_PROCESS_INCLUDES Cairo_INCLUDE_DIR Freetype_INCLUDE_DIRS) +set(Cairo_PROCESS_LIBS Cairo_LIBRARY Freetype_LIBRARIES) +libfind_process(Cairo) + diff --git a/cmake/FindCairomm.cmake b/cmake/FindCairomm.cmake new file mode 100644 index 0000000..d6c9fa1 --- /dev/null +++ b/cmake/FindCairomm.cmake @@ -0,0 +1,53 @@ +# - Try to find Cairomm 1.16 +# Once done, this will define +# +# Cairomm_FOUND - system has Cairomm +# Cairomm_INCLUDE_DIRS - the Cairomm include directories +# Cairomm_LIBRARIES - link these to use Cairomm + +include(LibFindMacros) + +# Dependencies +libfind_package(Cairomm Cairo) +libfind_package(Cairomm SigC++) + +foreach (dir ${CMAKE_SYSTEM_PREFIX_PATH}) + # message(INFO ${dir}) + file(GLOB Cairomm_VERSION_PATH "${dir}/include/cairomm*") + # message(INFO ${Cairomm_VERSION_PATH}) + if (NOT ${Cairomm_VERSION_PATH} STREQUAL "") + break() + endif () +endforeach () + +#message(INFO ${Cairomm_VERSION_PATH}) +get_filename_component(Cairomm_VERSION_DIR ${Cairomm_VERSION_PATH} NAME) +#message(INFO ${Cairomm_VERSION_DIR}) +string(REPLACE "cairomm-" "" Cairomm_VERSION ${Cairomm_VERSION_DIR}) +#message(INFO ${Cairomm_VERSION}) + + +# Use pkg-config to get hints about paths +libfind_pkg_check_modules(Cairomm_PKGCONF ${Cairomm_VERSION_DIR}) + +# Main include dir +find_path(Cairomm_INCLUDE_DIR + NAMES cairomm/cairomm.h + PATHS ${Cairomm_PKGCONF_INCLUDE_DIRS} + PATH_SUFFIXES ${Cairomm_VERSION_DIR} + ) + +find_path(Cairommconfig_INCLUDE_DIR + NAMES cairommconfig.h + PATHS ${Cairomm_PKGCONF_INCLUDE_DIRS} + PATH_SUFFIXES ${Cairomm_VERSION_DIR} + ) + +libfind_library(Cairomm cairomm ${Cairomm_VERSION}) + +# Set the include dir variables and the libraries and let libfind_process do the rest. +# NOTE: Singular variables for this library, plural for libraries this this lib depends on. +set(Cairomm_PROCESS_INCLUDES Cairomm_INCLUDE_DIR Cairommconfig_INCLUDE_DIR SigC++_INCLUDE_DIRS Cairo_INCLUDE_DIRS) +set(Cairomm_PROCESS_LIBS Cairomm_LIBRARY Cairo_LIBRARIES SigC++_LIBRARIES) +libfind_process(Cairomm) + diff --git a/cmake/FindEJDB.cmake b/cmake/FindEJDB.cmake new file mode 100644 index 0000000..27e8fb7 --- /dev/null +++ b/cmake/FindEJDB.cmake @@ -0,0 +1,31 @@ +# - Try to find EJDB +# Once done, this will define +# +# EJDB_FOUND - system has SigC++ +# --EJDB_INCLUDE_DIRS - the SigC++ include directories +# EJDB_LIBRARIES - link these to use SigC++ + +include(LibFindMacros) + +# Use pkg-config to get hints about paths +libfind_pkg_check_modules(EJDB_PKGCONF libejdb) + +# Include dir +find_path(EJDB_INCLUDE_DIR + NAMES ejdb.h + PATHS ${EJDB_PKGCONF_INCLUDE_DIRS} + PATH_SUFFIXES ejdb +) + +# Finally the library itself +find_library(EJDB_LIBRARY + NAMES ejdb + PATHS ${EJDB_PKGCONF_LIBRARY_DIRS} +) + +# Set the include dir variables and the libraries and let libfind_process do the rest. +# NOTE: Singular variables for this library, plural for libraries this this lib depends on. +set(EJDB_PROCESS_INCLUDES EJDB_INCLUDE_DIR) +set(EJDB_PROCESS_LIBS EJDB_LIBRARY) +libfind_process(EJDB) + diff --git a/cmake/FindFreetype.cmake b/cmake/FindFreetype.cmake new file mode 100644 index 0000000..b4c1685 --- /dev/null +++ b/cmake/FindFreetype.cmake @@ -0,0 +1,31 @@ +# - Try to find Freetype2 +# Once done, this will define +# +# Freetype_FOUND - system has Freetype +# Freetype_INCLUDE_DIRS - the Freetype include directories +# Freetype_LIBRARIES - link these to use Freetype + +include(LibFindMacros) + +# Use pkg-config to get hints about paths +libfind_pkg_check_modules(Freetype_PKGCONF freetype2) + +# Include dir +find_path(Freetype_INCLUDE_DIR + NAMES freetype/freetype.h + PATHS ${Freetype_PKGCONF_INCLUDE_DIRS} + PATH_SUFFIXES freetype2 +) + +# Finally the library itself +find_library(Freetype_LIBRARY + NAMES freetype + PATHS ${Freetype_PKGCONF_LIBRARY_DIRS} +) + +# Set the include dir variables and the libraries and let libfind_process do the rest. +# NOTE: Singular variables for this library, plural for libraries this this lib depends on. +set(Freetype_PROCESS_INCLUDES Freetype_INCLUDE_DIR) +set(Freetype_PROCESS_LIBS Freetype_LIBRARY) +libfind_process(Freetype) + diff --git a/cmake/FindKaHIP.cmake b/cmake/FindKaHIP.cmake new file mode 100644 index 0000000..08fc625 --- /dev/null +++ b/cmake/FindKaHIP.cmake @@ -0,0 +1,30 @@ +# - Try to find KaHIP +# Once done, this will define +# +# KaHIP_FOUND - system has KaHIP +# --KaHIP_INCLUDE_DIRS - the KaHIP include directories +# EJDB_LIBRARIES - link these to use KaHIP + +include(LibFindMacros) + +# Use pkg-config to get hints about paths +libfind_pkg_check_modules(KaHIP_PKGCONF libkahip) + +# Include dir +find_path(KaHIP_INCLUDE_DIR + NAMES kaHIP_interface.h + PATHS ${KaHIP_PKGCONF_INCLUDE_DIRS} +) + +# Finally the library itself +find_library(KaHIP_LIBRARY + NAMES kahip + PATHS ${KaHIP_PKGCONF_LIBRARY_DIRS} +) + +# Set the include dir variables and the libraries and let libfind_process do the rest. +# NOTE: Singular variables for this library, plural for libraries this this lib depends on. +set(KaHIP_PROCESS_INCLUDES KaHIP_INCLUDE_DIR) +set(KaHIP_PROCESS_LIBS KaHIP_LIBRARY) +libfind_process(KaHIP) + diff --git a/cmake/FindLibKDTree++.cmake b/cmake/FindLibKDTree++.cmake new file mode 100644 index 0000000..65848fe --- /dev/null +++ b/cmake/FindLibKDTree++.cmake @@ -0,0 +1,19 @@ +# - Try to find libkdtree++ +# Once done, this will define +# +# LibKDTree_FOUND - system has KaHIP +# --KaHIP_INCLUDE_DIRS - the KaHIP include directories +# EJDB_LIBRARIES - link these to use KaHIP + +# Include dir +find_path(LibKDTree_INCLUDE_DIR + NAMES kdtree++/kdtree.hpp +) + +include(FindPackageHandleStandardArgs) +# handle the QUIETLY and REQUIRED arguments and set LIBXML2_FOUND to TRUE +# if all listed variables are TRUE +find_package_handle_standard_args(LibKDTree DEFAULT_MSG + LibKDTree_INCLUDE_DIR) + +set(LibKDTree_INCLUDE_DIRS LibKDTree_INCLUDE_DIR) diff --git a/cmake/FindNanoflann.cmake b/cmake/FindNanoflann.cmake new file mode 100644 index 0000000..1dfcb87 --- /dev/null +++ b/cmake/FindNanoflann.cmake @@ -0,0 +1,18 @@ +# - Try to find nanoflann +# Once done, this will define +# +# Nanoflann_FOUND - system has Nanoflann +# Nanoflann_INCLUDE_DIRS - the Nanoflann include directories + +# Include dir +find_path(Nanoflann_INCLUDE_DIR + NAMES nanoflann.hpp +) + +include(FindPackageHandleStandardArgs) +# handle the QUIETLY and REQUIRED arguments and set LIBXML2_FOUND to TRUE +# if all listed variables are TRUE +find_package_handle_standard_args(Nanoflann DEFAULT_MSG + Nanoflann_INCLUDE_DIR) + +set(Nanoflann_INCLUDE_DIRS Nanoflann_INCLUDE_DIR) diff --git a/cmake/FindNuma.cmake b/cmake/FindNuma.cmake new file mode 100644 index 0000000..157b2b7 --- /dev/null +++ b/cmake/FindNuma.cmake @@ -0,0 +1,20 @@ +# Find the numa policy library. +# Output variables: +# NUMA_INCLUDE_DIR : e.g., /usr/include/. +# NUMA_LIBRARY : Library path of numa library +# NUMA_FOUND : True if found. +FIND_PATH(NUMA_INCLUDE_DIR NAME numa.h + HINTS $ENV{HOME}/local/include /opt/local/include /usr/local/include /usr/include) + +FIND_LIBRARY(NUMA_LIBRARY NAME numa + HINTS $ENV{HOME}/local/lib64 $ENV{HOME}/local/lib /usr/local/lib64 /usr/local/lib /opt/local/lib64 /opt/local/lib /usr/lib64 /usr/lib +) + +IF (NUMA_INCLUDE_DIR AND NUMA_LIBRARY) + SET(NUMA_FOUND TRUE) + MESSAGE(STATUS "Found numa library: inc=${NUMA_INCLUDE_DIR}, lib=${NUMA_LIBRARY}") +ELSE () + SET(NUMA_FOUND FALSE) + MESSAGE(STATUS "WARNING: Numa library not found.") + MESSAGE(STATUS "Try: 'sudo yum install numactl numactl-devel' (or sudo apt-get install libnuma libnuma-dev)") +ENDIF () diff --git a/cmake/FindSigC++.cmake b/cmake/FindSigC++.cmake new file mode 100644 index 0000000..966013a --- /dev/null +++ b/cmake/FindSigC++.cmake @@ -0,0 +1,49 @@ +# - Try to find SigC++-3.0 +# Once done, this will define +# +# SigC++_FOUND - system has SigC++ +# SigC++_INCLUDE_DIRS - the SigC++ include directories +# SigC++_LIBRARIES - link these to use SigC++ + +include(LibFindMacros) + +foreach (dir ${CMAKE_SYSTEM_PREFIX_PATH}) + # message(INFO ${dir}) + file(GLOB SigC++_VERSION_PATH "${dir}/include/sigc++*") + # message(INFO ${SigC++_VERSION_PATH}) + if (NOT ${SigC++_VERSION_PATH} STREQUAL "") + break() + endif () +endforeach () + +#message(INFO ${SigC++_VERSION_PATH}) +get_filename_component(SigC++_VERSION_DIR ${SigC++_VERSION_PATH} NAME) +#message(INFO ${SigC++_VERSION_DIR}) +string(REPLACE "sigc++-" "" SigC++_VERSION ${SigC++_VERSION_DIR}) +#message(INFO ${SigC++_VERSION}) + +# Use pkg-config to get hints about paths +libfind_pkg_check_modules(SigC++_PKGCONF ${SigC++_VERSION_DIR}) + +# Main include dir +find_path(SigC++_INCLUDE_DIR + NAMES sigc++/sigc++.h + PATHS ${SigC++_PKGCONF_INCLUDE_DIRS} + PATH_SUFFIXES ${SigC++_VERSION_DIR} + ) + +# Glib-related libraries also use a separate config header, which is in lib dir +find_path(SigC++Config_INCLUDE_DIR + NAMES sigc++config.h + PATHS ${SigC++_PKGCONF_INCLUDE_DIRS} /usr + PATH_SUFFIXES lib/${SigC++_VERSION_DIR}/include + ) + +libfind_library(SigC++ sigc ${SigC++_VERSION}) + +# Set the include dir variables and the libraries and let libfind_process do the rest. +# NOTE: Singular variables for this library, plural for libraries this this lib depends on. +set(SigC++_PROCESS_INCLUDES SigC++_INCLUDE_DIR SigC++Config_INCLUDE_DIR) +set(SigC++_PROCESS_LIBS SigC++_LIBRARY) +libfind_process(SigC++) + diff --git a/cmake/FindTBB.cmake b/cmake/FindTBB.cmake new file mode 100644 index 0000000..ccad261 --- /dev/null +++ b/cmake/FindTBB.cmake @@ -0,0 +1,291 @@ +# - Find ThreadingBuildingBlocks include dirs and libraries +# Use this module by invoking find_package with the form: +# find_package(TBB +# [REQUIRED] # Fail with error if TBB is not found +# ) # +# Once done, this will define +# +# TBB_FOUND - system has TBB +# TBB_INCLUDE_DIRS - the TBB include directories +# TBB_LIBRARIES - TBB libraries to be lined, doesn't include malloc or +# malloc proxy +# +# TBB_VERSION_MAJOR - Major Product Version Number +# TBB_VERSION_MINOR - Minor Product Version Number +# TBB_INTERFACE_VERSION - Engineering Focused Version Number +# TBB_COMPATIBLE_INTERFACE_VERSION - The oldest major interface version +# still supported. This uses the engineering +# focused interface version numbers. +# +# TBB_MALLOC_FOUND - system has TBB malloc library +# TBB_MALLOC_INCLUDE_DIRS - the TBB malloc include directories +# TBB_MALLOC_LIBRARIES - The TBB malloc libraries to be lined +# +# TBB_MALLOC_PROXY_FOUND - system has TBB malloc proxy library +# TBB_MALLOC_PROXY_INCLUDE_DIRS = the TBB malloc proxy include directories +# TBB_MALLOC_PROXY_LIBRARIES - The TBB malloc proxy libraries to be lined +# +# +# This module reads hints about search locations from variables: +# ENV TBB_ARCH_PLATFORM for windows only +# ENV TBB_ROOT or just TBB_ROOT +# +# +# +# Modified by Robert Maynard from the original OGRE source +# +#------------------------------------------------------------------- +# This file is part of the CMake build system for OGRE +# (Object-oriented Graphics Rendering Engine) +# For the latest info, see http://www.ogre3d.org/ +# +# The contents of this file are placed in the public domain. Feel +# free to make use of it in any way you like. +#------------------------------------------------------------------- +# +#============================================================================= +# Copyright 2010-2012 Kitware, Inc. +# Copyright 2012 Rolf Eike Beer +# +# Distributed under the OSI-approved BSD License (the "License"); +# see accompanying file Copyright.txt for details. +# +# This software is distributed WITHOUT ANY WARRANTY; without even the +# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +# See the License for more information. +#============================================================================= +# (To distribute this file outside of CMake, substitute the full +# License text for the above reference.) + + +#============================================================================= +# FindTBB helper functions and macros +# + +#=============================================== +# Create search paths based on prefix path +#=============================================== +macro(create_search_paths PREFIX) + foreach(dir ${${PREFIX}_PREFIX_PATH}) + set(${PREFIX}_INC_SEARCH_PATH ${${PREFIX}_INC_SEARCH_PATH} + ${dir}/include ${dir}/Include ${dir}/include/${PREFIX}) + set(${PREFIX}_LIB_SEARCH_PATH ${${PREFIX}_LIB_SEARCH_PATH} + ${dir}/lib ${dir}/Lib ${dir}/lib/${PREFIX} ${dir}/Libs) + endforeach(dir) +endmacro(create_search_paths) + +#=============================================== +# Do the final processing for the package find. +#=============================================== +macro(findpkg_finish PREFIX) + # skip if already processed during this run + if (NOT ${PREFIX}_FOUND) + if (${PREFIX}_INCLUDE_DIR AND ${PREFIX}_LIBRARY) + set(${PREFIX}_FOUND TRUE) + set (${PREFIX}_INCLUDE_DIRS ${${PREFIX}_INCLUDE_DIR}) + set (${PREFIX}_LIBRARIES ${${PREFIX}_LIBRARY}) + else () + if (${PREFIX}_FIND_REQUIRED AND NOT ${PREFIX}_FIND_QUIETLY) + message(FATAL_ERROR "Required library ${PREFIX} not found.") + endif () + endif () + + #mark the following variables as internal variables + mark_as_advanced(${PREFIX}_INCLUDE_DIR + ${PREFIX}_LIBRARY + ${PREFIX}_LIBRARY_DEBUG + ${PREFIX}_LIBRARY_RELEASE) + endif () +endmacro(findpkg_finish) + +#=============================================== +# Generate debug names from given RELEASEease names +#=============================================== +macro(get_debug_names PREFIX) + foreach(i ${${PREFIX}}) + set(${PREFIX}_DEBUG ${${PREFIX}_DEBUG} ${i}d ${i}D ${i}_d ${i}_D ${i}_debug ${i}) + endforeach(i) +endmacro(get_debug_names) + +#=============================================== +# See if we have env vars to help us find tbb +#=============================================== +macro(getenv_path VAR) + set(ENV_${VAR} $ENV{${VAR}}) + # replace won't work if var is blank + if (ENV_${VAR}) + string( REGEX REPLACE "\\\\" "/" ENV_${VAR} ${ENV_${VAR}} ) + endif () +endmacro(getenv_path) + +#=============================================== +# Couple a set of RELEASEease AND debug libraries +#=============================================== +macro(make_library_set PREFIX) + if (${PREFIX}_RELEASE AND ${PREFIX}_DEBUG) + set(${PREFIX} optimized ${${PREFIX}_RELEASE} debug ${${PREFIX}_DEBUG}) + elseif (${PREFIX}_RELEASE) + set(${PREFIX} ${${PREFIX}_RELEASE}) + elseif (${PREFIX}_DEBUG) + set(${PREFIX} ${${PREFIX}_DEBUG}) + endif () +endmacro(make_library_set) + + +#============================================================================= +# Now to actually find TBB +# + +# Get path, convert backslashes as ${ENV_${var}} +getenv_path(TBB_ROOT) +# construct search paths +set(TBB_PREFIX_PATH ${TBB_ROOT} ${ENV_TBB_ROOT}) +create_search_paths(TBB) + +# get the arch, only used by windows +if($ENV{TBB_ARCH_PLATFORM}) + set(TBB_ARCH_PLATFORM $ENV{TBB_ARCH_PLATFORM}) +endif() + +# For Windows, let's assume that the user might be using the precompiled +# TBB packages from the main website. These use a rather awkward directory +# structure (at least for automatically finding the right files) depending +# on platform and compiler, but we'll do our best to accommodate it. +# Not adding the same effort for the precompiled linux builds, though. Those +# have different versions for CC compiler versions and linux kernels which +# will never adequately match the user's setup, so there is no feasible way +# to detect the "best" version to use. The user will have to manually +# select the right files. (Chances are the distributions are shipping their +# custom version of tbb, anyway, so the problem is probably nonexistant.) +if (WIN32 AND MSVC) + set(COMPILER_PREFIX "vc7.1") + if (MSVC_VERSION EQUAL 1400) + set(COMPILER_PREFIX "vc8") + elseif(MSVC_VERSION EQUAL 1500) + set(COMPILER_PREFIX "vc9") + elseif(MSVC_VERSION EQUAL 1600) + set(COMPILER_PREFIX "vc10") + elseif(MSVC_VERSION EQUAL 1700) + set(COMPILER_PREFIX "vc11") + elseif(MSVC_VERSION EQUAL 1800) + set(COMPILER_PREFIX "vc12") + endif () + + # for each prefix path, add ia32/64\${COMPILER_PREFIX}\lib to the lib search path + foreach (dir ${TBB_PREFIX_PATH}) + if (CMAKE_CL_64) + list(APPEND TBB_LIB_SEARCH_PATH ${dir}/ia64/${COMPILER_PREFIX}/lib) + list(APPEND TBB_LIB_SEARCH_PATH ${dir}/lib/ia64/${COMPILER_PREFIX}) + list(APPEND TBB_LIB_SEARCH_PATH ${dir}/intel64/${COMPILER_PREFIX}/lib) + list(APPEND TBB_LIB_SEARCH_PATH ${dir}/lib/intel64/${COMPILER_PREFIX}) + else () + list(APPEND TBB_LIB_SEARCH_PATH ${dir}/ia32/${COMPILER_PREFIX}/lib) + list(APPEND TBB_LIB_SEARCH_PATH ${dir}/lib/ia32/${COMPILER_PREFIX}) + endif () + endforeach () +endif () + +foreach (dir ${TBB_PREFIX_PATH}) + list(APPEND TBB_LIB_SEARCH_PATH ${dir}/${TBB_ARCH_PLATFORM}/lib) + list(APPEND TBB_LIB_SEARCH_PATH ${dir}/lib/${TBB_ARCH_PLATFORM}) + list(APPEND TBB_LIB_SEARCH_PATH ${dir}/lib) +endforeach () + + +set(TBB_LIBRARY_NAMES tbb) +get_debug_names(TBB_LIBRARY_NAMES) + + +find_path(TBB_INCLUDE_DIR + NAMES tbb/tbb.h + PATHS ${TBB_INC_SEARCH_PATH}) + +find_library(TBB_LIBRARY_RELEASE + NAMES ${TBB_LIBRARY_NAMES} + PATHS ${TBB_LIB_SEARCH_PATH}) +find_library(TBB_LIBRARY_DEBUG + NAMES ${TBB_LIBRARY_NAMES_DEBUG} + PATHS ${TBB_LIB_SEARCH_PATH}) +make_library_set(TBB_LIBRARY) + +findpkg_finish(TBB) + +#on unix we need to also link to rt +if(UNIX AND NOT APPLE) + list(APPEND TBB_LIBRARIES rt) +endif() + +#if we haven't found TBB no point on going any further +if (NOT TBB_FOUND) + return() +endif () + +#============================================================================= +# Look for TBB's malloc package +set(TBB_MALLOC_LIBRARY_NAMES tbbmalloc) +get_debug_names(TBB_MALLOC_LIBRARY_NAMES) + +find_path(TBB_MALLOC_INCLUDE_DIR + NAMES tbb/tbb.h + PATHS ${TBB_INC_SEARCH_PATH}) + +find_library(TBB_MALLOC_LIBRARY_RELEASE + NAMES ${TBB_MALLOC_LIBRARY_NAMES} + PATHS ${TBB_LIB_SEARCH_PATH}) +find_library(TBB_MALLOC_LIBRARY_DEBUG + NAMES ${TBB_MALLOC_LIBRARY_NAMES_DEBUG} + PATHS ${TBB_LIB_SEARCH_PATH}) +make_library_set(TBB_MALLOC_LIBRARY) + +findpkg_finish(TBB_MALLOC) + +#============================================================================= +# Look for TBB's malloc proxy package +set(TBB_MALLOC_PROXY_LIBRARY_NAMES tbbmalloc_proxy) +get_debug_names(TBB_MALLOC_PROXY_LIBRARY_NAMES) + +find_path(TBB_MALLOC_PROXY_INCLUDE_DIR + NAMES tbb/tbbmalloc_proxy.h + PATHS ${TBB_INC_SEARCH_PATH}) + +find_library(TBB_MALLOC_PROXY_LIBRARY_RELEASE + NAMES ${TBB_MALLOC_PROXY_LIBRARY_NAMES} + PATHS ${TBB_LIB_SEARCH_PATH}) +find_library(TBB_MALLOC_PROXY_LIBRARY_DEBUG + NAMES ${TBB_MALLOC_PROXY_LIBRARY_NAMES_DEBUG} + PATHS ${TBB_LIB_SEARCH_PATH}) +make_library_set(TBB_MALLOC_PROXY_LIBRARY) + +findpkg_finish(TBB_MALLOC_PROXY) + +#----------------------------------------------------------------------------- +# setup timing libs we need to link too + + +#============================================================================= +#parse all the version numbers from tbb +if(NOT TBB_VERSION) + + #only read the start of the file + file(READ + "${TBB_INCLUDE_DIR}/tbb/tbb_stddef.h" + TBB_VERSION_CONTENTS + LIMIT 2048) + + string(REGEX REPLACE + ".*#define TBB_VERSION_MAJOR ([0-9]+).*" "\\1" + TBB_VERSION_MAJOR "${TBB_VERSION_CONTENTS}") + + string(REGEX REPLACE + ".*#define TBB_VERSION_MINOR ([0-9]+).*" "\\1" + TBB_VERSION_MINOR "${TBB_VERSION_CONTENTS}") + + string(REGEX REPLACE + ".*#define TBB_INTERFACE_VERSION ([0-9]+).*" "\\1" + TBB_INTERFACE_VERSION "${TBB_VERSION_CONTENTS}") + + string(REGEX REPLACE + ".*#define TBB_COMPATIBLE_INTERFACE_VERSION ([0-9]+).*" "\\1" + TBB_COMPATIBLE_INTERFACE_VERSION "${TBB_VERSION_CONTENTS}") + +endif() diff --git a/cmake/FindVTune.cmake b/cmake/FindVTune.cmake new file mode 100644 index 0000000..c786e31 --- /dev/null +++ b/cmake/FindVTune.cmake @@ -0,0 +1,45 @@ + +# Copyright (c) 2013-2014 Stefan.Eilemann@epfl.ch +# finds the ittnotify API + +# VTUNE_FOUND - system has VTUNE +# VTUNE_INCLUDE_DIRS - the VTUNE include directories +# VTUNE_LIBRARIES - link these to use VTUNE + +include(LibFindMacros) + +find_program(VTUNE_EXECUTABLE amplxe-cl + HINTS $ENV{VTUNE_AMPLIFIER_XE_2015_DIR} + PATH_SUFFIXES bin64 +) + +if(NOT VTUNE_EXECUTABLE) + return() +endif() + +get_filename_component(VTUNE_DIR ${VTUNE_EXECUTABLE} PATH) +set(VTUNE_DIR "${VTUNE_DIR}/..") + +# Include dir +find_path(VTUNE_INCLUDE_DIR + NAMES ittnotify.h + PATHS ${VTUNE_DIR} + PATH_SUFFIXES include +) + +# Finally the library itself +find_library(VTUNE_LIBRARY + NAMES ittnotify + PATHS ${VTUNE_DIR} + PATH_SUFFIXES lib64 +) + +if(UNIX) + list(APPEND VTUNE_LIBRARY dl) +endif() + +# Set the include dir variables and the libraries and let libfind_process do the rest. +# NOTE: Singular variables for this library, plural for libraries this this lib depends on. +set(VTUNE_PROCESS_INCLUDES VTUNE_INCLUDE_DIR) +set(VTUNE_PROCESS_LIBS VTUNE_LIBRARY) +libfind_process(VTUNE) \ No newline at end of file diff --git a/cmake/GetGitRevisionDescription.cmake b/cmake/GetGitRevisionDescription.cmake new file mode 100644 index 0000000..c8d27f2 --- /dev/null +++ b/cmake/GetGitRevisionDescription.cmake @@ -0,0 +1,130 @@ +# - Returns a version string from Git +# +# These functions force a re-configure on each git commit so that you can +# trust the values of the variables in your build system. +# +# get_git_head_revision( [ ...]) +# +# Returns the refspec and sha hash of the current head revision +# +# git_describe( [ ...]) +# +# Returns the results of git describe on the source tree, and adjusting +# the output so that it tests false if an error occurs. +# +# git_get_exact_tag( [ ...]) +# +# Returns the results of git describe --exact-match on the source tree, +# and adjusting the output so that it tests false if there was no exact +# matching tag. +# +# Requires CMake 2.6 or newer (uses the 'function' command) +# +# Original Author: +# 2009-2010 Ryan Pavlik +# http://academic.cleardefinition.com +# Iowa State University HCI Graduate Program/VRAC +# +# Copyright Iowa State University 2009-2010. +# Distributed under the Boost Software License, Version 1.0. +# (See accompanying file LICENSE_1_0.txt or copy at +# http://www.boost.org/LICENSE_1_0.txt) + +if(__get_git_revision_description) + return() +endif() +set(__get_git_revision_description YES) + +# We must run the following at "include" time, not at function call time, +# to find the path to this module rather than the path to a calling list file +get_filename_component(_gitdescmoddir ${CMAKE_CURRENT_LIST_FILE} PATH) + +function(get_git_head_revision _refspecvar _hashvar) + set(GIT_PARENT_DIR "${CMAKE_CURRENT_SOURCE_DIR}") + set(GIT_DIR "${GIT_PARENT_DIR}/.git") + while(NOT EXISTS "${GIT_DIR}") # .git dir not found, search parent directories + set(GIT_PREVIOUS_PARENT "${GIT_PARENT_DIR}") + get_filename_component(GIT_PARENT_DIR ${GIT_PARENT_DIR} PATH) + if(GIT_PARENT_DIR STREQUAL GIT_PREVIOUS_PARENT) + # We have reached the root directory, we are not in git + set(${_refspecvar} "GITDIR-NOTFOUND" PARENT_SCOPE) + set(${_hashvar} "GITDIR-NOTFOUND" PARENT_SCOPE) + return() + endif() + set(GIT_DIR "${GIT_PARENT_DIR}/.git") + endwhile() + # check if this is a submodule + if(NOT IS_DIRECTORY ${GIT_DIR}) + file(READ ${GIT_DIR} submodule) + string(REGEX REPLACE "gitdir: (.*)\n$" "\\1" GIT_DIR_RELATIVE ${submodule}) + get_filename_component(SUBMODULE_DIR ${GIT_DIR} PATH) + get_filename_component(GIT_DIR ${SUBMODULE_DIR}/${GIT_DIR_RELATIVE} ABSOLUTE) + endif() + set(GIT_DATA "${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/git-data") + if(NOT EXISTS "${GIT_DATA}") + file(MAKE_DIRECTORY "${GIT_DATA}") + endif() + + if(NOT EXISTS "${GIT_DIR}/HEAD") + return() + endif() + set(HEAD_FILE "${GIT_DATA}/HEAD") + configure_file("${GIT_DIR}/HEAD" "${HEAD_FILE}" COPYONLY) + + configure_file("${_gitdescmoddir}/GetGitRevisionDescription.cmake.in" + "${GIT_DATA}/grabRef.cmake" + @ONLY) + include("${GIT_DATA}/grabRef.cmake") + + set(${_refspecvar} "${HEAD_REF}" PARENT_SCOPE) + set(${_hashvar} "${HEAD_HASH}" PARENT_SCOPE) +endfunction() + +function(git_describe _var) + if(NOT GIT_FOUND) + find_package(Git QUIET) + endif() + get_git_head_revision(refspec hash) + if(NOT GIT_FOUND) + set(${_var} "GIT-NOTFOUND" PARENT_SCOPE) + return() + endif() + if(NOT hash) + set(${_var} "HEAD-HASH-NOTFOUND" PARENT_SCOPE) + return() + endif() + + # TODO sanitize + #if((${ARGN}" MATCHES "&&") OR + # (ARGN MATCHES "||") OR + # (ARGN MATCHES "\\;")) + # message("Please report the following error to the project!") + # message(FATAL_ERROR "Looks like someone's doing something nefarious with git_describe! Passed arguments ${ARGN}") + #endif() + + #message(STATUS "Arguments to execute_process: ${ARGN}") + + execute_process(COMMAND + "${GIT_EXECUTABLE}" + describe + ${hash} + ${ARGN} + WORKING_DIRECTORY + "${CMAKE_SOURCE_DIR}" + RESULT_VARIABLE + res + OUTPUT_VARIABLE + out + ERROR_QUIET + OUTPUT_STRIP_TRAILING_WHITESPACE) + if(NOT res EQUAL 0) + set(out "${out}-${res}-NOTFOUND") + endif() + + set(${_var} "${out}" PARENT_SCOPE) +endfunction() + +function(git_get_exact_tag _var) + git_describe(out --exact-match ${ARGN}) + set(${_var} "${out}" PARENT_SCOPE) +endfunction() diff --git a/cmake/GetGitRevisionDescription.cmake.in b/cmake/GetGitRevisionDescription.cmake.in new file mode 100644 index 0000000..888ce13 --- /dev/null +++ b/cmake/GetGitRevisionDescription.cmake.in @@ -0,0 +1,38 @@ +# +# Internal file for GetGitRevisionDescription.cmake +# +# Requires CMake 2.6 or newer (uses the 'function' command) +# +# Original Author: +# 2009-2010 Ryan Pavlik +# http://academic.cleardefinition.com +# Iowa State University HCI Graduate Program/VRAC +# +# Copyright Iowa State University 2009-2010. +# Distributed under the Boost Software License, Version 1.0. +# (See accompanying file LICENSE_1_0.txt or copy at +# http://www.boost.org/LICENSE_1_0.txt) + +set(HEAD_HASH) + +file(READ "@HEAD_FILE@" HEAD_CONTENTS LIMIT 1024) + +string(STRIP "${HEAD_CONTENTS}" HEAD_CONTENTS) +if(HEAD_CONTENTS MATCHES "ref") + # named branch + string(REPLACE "ref: " "" HEAD_REF "${HEAD_CONTENTS}") + if(EXISTS "@GIT_DIR@/${HEAD_REF}") + configure_file("@GIT_DIR@/${HEAD_REF}" "@GIT_DATA@/head-ref" COPYONLY) + elseif(EXISTS "@GIT_DIR@/logs/${HEAD_REF}") + configure_file("@GIT_DIR@/logs/${HEAD_REF}" "@GIT_DATA@/head-ref" COPYONLY) + set(HEAD_HASH "${HEAD_REF}") + endif() +else() + # detached HEAD + configure_file("@GIT_DIR@/HEAD" "@GIT_DATA@/head-ref" COPYONLY) +endif() + +if(NOT HEAD_HASH) + file(READ "@GIT_DATA@/head-ref" HEAD_HASH LIMIT 1024) + string(STRIP "${HEAD_HASH}" HEAD_HASH) +endif() diff --git a/cmake/LibFindMacros.cmake b/cmake/LibFindMacros.cmake new file mode 100644 index 0000000..49899a8 --- /dev/null +++ b/cmake/LibFindMacros.cmake @@ -0,0 +1,101 @@ +# Works the same as find_package, but forwards the "REQUIRED" and "QUIET" arguments +# used for the current package. For this to work, the first parameter must be the +# prefix of the current package, then the prefix of the new package etc, which are +# passed to find_package. +macro (libfind_package PREFIX) + set (LIBFIND_PACKAGE_ARGS ${ARGN}) + if (${PREFIX}_FIND_QUIETLY) + set (LIBFIND_PACKAGE_ARGS ${LIBFIND_PACKAGE_ARGS} QUIET) + endif (${PREFIX}_FIND_QUIETLY) + if (${PREFIX}_FIND_REQUIRED) + set (LIBFIND_PACKAGE_ARGS ${LIBFIND_PACKAGE_ARGS} REQUIRED) + endif (${PREFIX}_FIND_REQUIRED) + find_package(${LIBFIND_PACKAGE_ARGS}) +endmacro (libfind_package) + +# CMake developers made the UsePkgConfig system deprecated in the same release (2.6) +# where they added pkg_check_modules. Consequently I need to support both in my scripts +# to avoid those deprecated warnings. Here's a helper that does just that. +# Works identically to pkg_check_modules, except that no checks are needed prior to use. +macro (libfind_pkg_check_modules PREFIX PKGNAME) + if (${CMAKE_MAJOR_VERSION} EQUAL 2 AND ${CMAKE_MINOR_VERSION} EQUAL 4) + include(UsePkgConfig) + pkgconfig(${PKGNAME} ${PREFIX}_INCLUDE_DIRS ${PREFIX}_LIBRARY_DIRS ${PREFIX}_LDFLAGS ${PREFIX}_CFLAGS) + else (${CMAKE_MAJOR_VERSION} EQUAL 2 AND ${CMAKE_MINOR_VERSION} EQUAL 4) + find_package(PkgConfig) + if (PKG_CONFIG_FOUND) + pkg_check_modules(${PREFIX} ${PKGNAME}) + endif (PKG_CONFIG_FOUND) + endif (${CMAKE_MAJOR_VERSION} EQUAL 2 AND ${CMAKE_MINOR_VERSION} EQUAL 4) +endmacro (libfind_pkg_check_modules) + +# Do the final processing once the paths have been detected. +# If include dirs are needed, ${PREFIX}_PROCESS_INCLUDES should be set to contain +# all the variables, each of which contain one include directory. +# Ditto for ${PREFIX}_PROCESS_LIBS and library files. +# Will set ${PREFIX}_FOUND, ${PREFIX}_INCLUDE_DIRS and ${PREFIX}_LIBRARIES. +# Also handles errors in case library detection was required, etc. +macro (libfind_process PREFIX) + # Skip processing if already processed during this run + if (NOT ${PREFIX}_FOUND) + # Start with the assumption that the library was found + set (${PREFIX}_FOUND TRUE) + + # Process all includes and set _FOUND to false if any are missing + foreach (i ${${PREFIX}_PROCESS_INCLUDES}) + if (${i}) + set (${PREFIX}_INCLUDE_DIRS ${${PREFIX}_INCLUDE_DIRS} ${${i}}) + mark_as_advanced(${i}) + else (${i}) + message(STATUS "Missing ${i}") + set (${PREFIX}_FOUND FALSE) + endif (${i}) + endforeach (i) + + # Process all libraries and set _FOUND to false if any are missing + foreach (i ${${PREFIX}_PROCESS_LIBS}) + if (${i}) + set (${PREFIX}_LIBRARIES ${${PREFIX}_LIBRARIES} ${${i}}) + mark_as_advanced(${i}) + else (${i}) + message(STATUS "Missing ${i}") + set (${PREFIX}_FOUND FALSE) + endif (${i}) + endforeach (i) + + # Print message and/or exit on fatal error + if (${PREFIX}_FOUND) + if (NOT ${PREFIX}_FIND_QUIETLY) + message (STATUS "Found ${PREFIX} ${${PREFIX}_VERSION}") + endif (NOT ${PREFIX}_FIND_QUIETLY) + else (${PREFIX}_FOUND) + if (${PREFIX}_FIND_REQUIRED) + foreach (i ${${PREFIX}_PROCESS_INCLUDES} ${${PREFIX}_PROCESS_LIBS}) + message("${i}=${${i}}") + endforeach (i) + message (FATAL_ERROR "Required library ${PREFIX} NOT FOUND.\nInstall the library (dev version) and try again. If the library is already installed, use ccmake to set the missing variables manually.") + endif (${PREFIX}_FIND_REQUIRED) + endif (${PREFIX}_FOUND) + endif (NOT ${PREFIX}_FOUND) +endmacro (libfind_process) + +macro(libfind_library PREFIX basename) + set(TMP "") + if(MSVC80) + set(TMP -vc80) + endif(MSVC80) + if(MSVC90) + set(TMP -vc90) + endif(MSVC90) + set(${PREFIX}_LIBNAMES ${basename}${TMP}) + if(${ARGC} GREATER 2) + set(${PREFIX}_LIBNAMES ${basename}${TMP}-${ARGV2}) + string(REGEX REPLACE "\\." "_" TMP ${${PREFIX}_LIBNAMES}) + set(${PREFIX}_LIBNAMES ${${PREFIX}_LIBNAMES} ${TMP}) + endif(${ARGC} GREATER 2) + find_library(${PREFIX}_LIBRARY + NAMES ${${PREFIX}_LIBNAMES} + PATHS ${${PREFIX}_PKGCONF_LIBRARY_DIRS} + ) +endmacro(libfind_library) + diff --git a/include/Algorithms/CGAL_Delaunay.hpp b/include/Algorithms/CGAL_Delaunay.hpp new file mode 100644 index 0000000..5cc5fd7 --- /dev/null +++ b/include/Algorithms/CGAL_Delaunay.hpp @@ -0,0 +1,37 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include +#include + +#include "Types.hpp" + +class CGAL_Delaunay2D { + +public: + using K = CGAL::Exact_predicates_inexact_constructions_kernel; + using DT = CGAL::Delaunay_triangulation_2; + using Point = DT::Point; + + static std::string name() { + return "CGALDelaunay"; + } + +public: + auto operator()(const tPoints &points) { + + DT dt; + + auto transform = [](const auto &p) -> auto { + return Point(p[0], p[1]); + }; + + dt.insert(boost::make_transform_iterator(points.begin(), transform), + boost::make_transform_iterator(points.end(), transform)); + + return dt; + } +}; \ No newline at end of file diff --git a/include/Algorithms/CGAL_Theta.hpp b/include/Algorithms/CGAL_Theta.hpp new file mode 100644 index 0000000..62eb298 --- /dev/null +++ b/include/Algorithms/CGAL_Theta.hpp @@ -0,0 +1,45 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include + +#include "Kernels/CGALKernel.hpp" +#include "Types.hpp" + +template +class CGAL_Theta2D { + +public: + using K = Kernel; + using Point = typename K::Point_2; + using Direction = typename K::Direction_2; + using Graph = boost::adjacency_list; + using Theta = CGAL::Construct_theta_graph_2; + + static std::string name() { + return "CGALTheta"; + } + +public: + auto operator()(const tDim &K, const tPoints &points, [[maybe_unused]] const tBox &bounds) { + + Theta theta(K); + Graph g; + + auto transform = [](const auto &p) -> auto{ + return Point(p[0], p[1]); + }; + + theta(boost::make_transform_iterator(points.begin(), transform), + boost::make_transform_iterator(points.end(), transform), + g); + + return g; + } +}; \ No newline at end of file diff --git a/include/Algorithms/CGAL_Yao.hpp b/include/Algorithms/CGAL_Yao.hpp new file mode 100644 index 0000000..5caf2a0 --- /dev/null +++ b/include/Algorithms/CGAL_Yao.hpp @@ -0,0 +1,46 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include + +#include "Kernels/CGALKernel.hpp" +#include "Types.hpp" + +template +class CGAL_Yao2D { + +public: + using K = Kernel; + using Point = typename K::Point_2; + using Direction = typename K::Direction_2; + using Graph = boost::adjacency_list; + using Yao = CGAL::Construct_yao_graph_2; + + static std::string name() { + return "CGALYao"; + } + +public: + auto operator()(const tDim &K, const tPoints &points, [[maybe_unused]] const tBox &bounds) { + + Yao yao(K); + Graph g; + + auto transform = [](const auto &p) -> auto{ + return Point(p[0], p[1]); + }; + + yao(boost::make_transform_iterator(points.begin(), transform), + boost::make_transform_iterator(points.end(), transform), + g); + + return tYaoGraph(points.size(), K); + //return g; + } +}; \ No newline at end of file diff --git a/include/Algorithms/GridYao.hpp b/include/Algorithms/GridYao.hpp new file mode 100644 index 0000000..3da6bcc --- /dev/null +++ b/include/Algorithms/GridYao.hpp @@ -0,0 +1,185 @@ +// +// Created by Daniel Funke on 18.03.21. +// + +#pragma once + +#include "Predicates.hpp" +#include "Types.hpp" + +template +class GridYao { + +private: + using tEFloat = typename Kernel::Float; + using tEFloatVector = std::array; + + using tKPoint = typename Kernel::Point; + using tKPoints = std::vector; + + class Grid { + + private: + friend class GridYao; + + using tGridCell = std::vector; + using tGrid = std::vector; + + public: + Grid(const tBox &_bounds, const tIndex _nCells) : bounds(_bounds) { + numCellsPerDim = std::ceil(std::sqrt(_nCells)); + numCells = numCellsPerDim * numCellsPerDim; + + for (tDim d = 0; d < bounds.high.size(); ++d) { + cellSize[d] = (bounds.high[d] - bounds.low[d]) / numCellsPerDim; + } + minCellSize = *std::min_element(cellSize.begin(), cellSize.end()); + + cells.resize(numCells); + } + + void buildGrid(const tKPoints &points) { + for (tIndex i = 0; i < points.size(); ++i) { + tIndex idx = getIndex(points[i]); + cells[idx].push_back(i); + } + } + + tIndexVector getIndexVector(const tKPoint &p) const { + tIndexVector idx; + for (tDim d = 0; d < idx.size(); ++d) { + idx[d] = std::floor(Kernel::to_float_exact((p[d] - bounds.low[d]) / cellSize[d]));//TODO exact? + + if (idx[d] == numCellsPerDim) [[unlikely]] {// if point is exactly on boundary + idx[d]--; + } + } + return idx; + } + + tIndex getIndex(const tIndexVector &idx) const { + tIndex i = 0; + for (tDim d = idx.size() - 1; d < idx.size(); --d) { + i += std::pow(numCellsPerDim, d) * idx[d]; + } + return i; + } + + tIndex getIndex(const tKPoint &p) const { + return getIndex(getIndexVector(p)); + } + + const tGridCell &operator[](const tIndex &i) const { + return cells[i]; + } + + tGridCell &operator[](const tIndex &i) { + return cells[i]; + } + + const tGridCell &operator[](const tIndexVector &i) const { + return cells[getIndex(i)]; + } + + tGridCell &operator[](const tIndexVector &i) { + return cells[getIndex(i)]; + } + + private: + tIndex numCellsPerDim; + tIndex numCells; + tEFloatVector cellSize; + tEFloat minCellSize; + tBox bounds; + tGrid cells; + }; + +public: + static std::string name() { + return "GridYao"; + } + + auto operator()(const tDim &K, const tPoints &iPoints, const tBox &bounds, const tIndex &cellOcc) const { + tYaoGraph graph(iPoints.size(), K); + + auto rays = Kernel::computeCones(K); + auto kPoints = Kernel::mkPoints(iPoints); + + Grid grid(bounds, std::ceil(kPoints.size() / static_cast(cellOcc))); + grid.buildGrid(kPoints); + + for (tIndex i = 0; i < kPoints.size(); ++i) { + + auto cLines = Kernel::computePointCones(kPoints[i], rays); + + for (tIndex radius = 0; + !isFinalized(graph[i], radius, grid.minCellSize) && radius < grid.numCellsPerDim; ++radius) { + searchRadius(i, radius, graph, grid, kPoints, cLines); + } + } + + return graph; + } + +private: + bool isFinalized(const tYaoVertex &v, const tIndex &radius, const tEFloat &minCellSize) const { + + auto d = static_cast(radius - 1) * minCellSize;//TODO better type handling, exact computation; + auto d2 = d * d; + + for (tDim k = 0; k < v.neighbor.size(); ++k) { + if (v.neighbor[k] == INF_IDX) { + // no neighbor found so far + return false; + } else { + // we found a neighbor for given search radius, can we find a closer one still + if (v.distance[k] > d2) { + // found neighbor is further away than any point for next cell radius could be + return false; + } + } + } + + return true; + } + + void searchRadius(const tIndex &point, const int &radius, + tYaoGraph &graph, const Grid &grid, const tKPoints &points, const auto &cLines) const { + + auto idx = grid.getIndexVector(points[point]); + + for (int dx = -radius; dx <= radius; ++dx) { + for (int dy = -radius; dy <= radius; ++dy) { + + if (std::max(std::abs(dx), std::abs(dy)) < radius) { + // skip inner cells, already visited in prior calls; + continue; + } + + auto cell = idx + std::array({dx, dy}); + if (std::all_of(cell.begin(), cell.end(), + [&grid](const auto &x) { return x < grid.numCellsPerDim; })) { + visitGridCell(point, cell, graph, grid, points, cLines); + } + } + } + } + + void visitGridCell(const tIndex &point, const tIndexVector &cell, + tYaoGraph &graph, const Grid &grid, const tKPoints &points, const auto &cLines) const { + + auto &vertex = graph[point]; + + for (const auto &p : grid[cell]) { + + if (p == point) continue; + + auto sec = Kernel::getCone(points[point], points[p], cLines); + + if (vertex.neighbor[sec] == INF_IDX || Kernel::compareDistance(points[point], points[p], points[vertex.neighbor[sec]], vertex.distance[sec])) { + vertex.neighbor[sec] = p; + vertex.distance[sec] = Kernel::to_float(Kernel::distance2(points[point], points[p])); + } + } + } +}; \ No newline at end of file diff --git a/include/Algorithms/NaiveYao.hpp b/include/Algorithms/NaiveYao.hpp new file mode 100644 index 0000000..7fd8089 --- /dev/null +++ b/include/Algorithms/NaiveYao.hpp @@ -0,0 +1,42 @@ +// +// Created by Daniel Funke on 18.03.21. +// + +#pragma once + +#include "Types.hpp" + +template +class NaiveYao { + +public: + static std::string name() { + return "NaiveYao"; + } + + auto operator()(const tDim &K, const tPoints &iPoints, [[maybe_unused]] const tBox &bounds) const { + tYaoGraph g(iPoints.size(), K); + + auto rays = Kernel::computeCones(K); + auto kPoints = Kernel::mkPoints(iPoints); + + for (tIndex i = 0; i < kPoints.size(); ++i) { + + auto cLines = Kernel::computePointCones(kPoints[i], rays); + + for (tIndex j = 0; j < kPoints.size(); ++j) { + + if (i == j) continue; + + auto sec = Kernel::getCone(kPoints[i], kPoints[j], cLines); + + if (g[i].neighbor[sec] == INF_IDX || Kernel::compareDistance(kPoints[i], kPoints[j], kPoints[g[i].neighbor[sec]], g[i].distance[sec])) { + g[i].neighbor[sec] = j; + g[i].distance[sec] = Kernel::to_float(Kernel::distance2(kPoints[i], kPoints[j])); + } + } + } + + return g; + } +}; \ No newline at end of file diff --git a/include/Algorithms/SweepLine.hpp b/include/Algorithms/SweepLine.hpp new file mode 100644 index 0000000..17de5e9 --- /dev/null +++ b/include/Algorithms/SweepLine.hpp @@ -0,0 +1,824 @@ +// +// Created by Daniel Funke on 18.03.21. +// + +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include "Predicates.hpp" +#include "Types.hpp" +#include "Utils/ListIndexTree.hpp" +#include "Utils/PriorityQueue.hpp" + +#include "Utils/Logging.hpp" + +#include "contrib/robin_hood.h" + +#ifdef WITH_CAIRO +#include "Painter.hpp" +//#define PAINT_STEPS +#endif + +template +class SweepLine { + + using tEFloat = typename Kernel::Float;// possibly exact float + using tDirection = typename Kernel::Direction; + using tPoint = typename Kernel::Point; + + using tRay = typename Kernel::Ray; + + struct SweeplineDS { + + using tRays = SearchTree; + using tRayHandle = typename tRays::Iterator; + + tDirection slDirection; + tRays slRays; + +#ifdef WITH_STATS + tIndex m_size = 0; +#endif + + SweeplineDS(const tDirection slDir_) + : slDirection(slDir_) {} + + + tRayHandle begin() { return slRays.begin(); } + tRayHandle end() { return slRays.end(); } + + auto insert(tRayHandle &pos, const tRay &ray) { +#ifdef WITH_STATS + m_size++; +#endif + return slRays.insert(pos, ray); + } + + auto insert_pair(tRayHandle &pos, const tRay &left, const tRay &right) { +#ifdef WITH_STATS + m_size += 2; +#endif + return slRays.insert_pair(pos, left, right); + } + + auto erase(tRayHandle &pos) { +#ifdef WITH_STATS + m_size--; +#endif + return slRays.erase(pos); + } + + auto replace(tRayHandle &pos, const tRay &ray) { + return slRays.replace(pos, ray); + } + +#ifdef WITH_STATS + auto size() const { + return m_size; + } +#endif + + tEFloat prj(const tPoint &p) { + return slDirection.prj(p); + } + + tRayHandle linearFind(const tPoint &p) {// TODO: make const? + auto exp = begin(); + while (exp != end() && exp->orientedSide(p) == tOrientedSide::LEFT) { + ++exp; + } + + return exp; + } + + // return iterator to ray immediatly right of point p + tRayHandle find(const tPoint &p) {// TODO: make const? + + auto cmp = [](const tPoint &p, const tRay &r) { + return r.orientedSide(p) != tOrientedSide::LEFT; + }; + + auto it = slRays.find(p, cmp); + + ASSERT(it == end() || it->orientedSide(p) != tOrientedSide::LEFT); + ASSERT(std::prev(it) == end() || std::prev(it)->orientedSide(p) != tOrientedSide::RIGHT); + ASSERT(it == linearFind(p)); + + return it; + } + +#ifdef WITH_CAIRO + + void draw(const tPoint &pos, Painter &painter) { + + // draw sweepline + painter.setDash(); + painter.drawLine(pos, {pos[X] + std::cos(slDirection.angle() + tIFloat(M_PI_2)), + pos[Y] + std::sin(slDirection.angle() + tIFloat(M_PI_2))}); + painter.drawLine(pos, {pos[X] + std::cos(slDirection.angle() - tIFloat(M_PI_2)), + pos[Y] + std::sin(slDirection.angle() - tIFloat(M_PI_2))}); + painter.unsetDash(); + + // draw direction + painter.setDash(); + //painter.drawLine(pos, {pos[X] + slDirVector[X], pos[Y] + slDirVector[Y]}); + painter.unsetDash(); + + // draw rays + for (auto &r : slRays) { + r.draw(painter); + } + } + +#endif + }; + + struct Event { + + enum Type { Input, + Intersection, + Deletion }; + + using tRayHandle = typename SweeplineDS::tRayHandle; + + Event(const tPoint &pos_) + : type(Type::Deletion), pos(pos_), idx(INF_IDX) {} + + Event(const tPoint &pos_, const tIndex &idx_) + : type(Type::Input), pos(pos_), idx(idx_) {} + + Event(const tPoint &pos_, const tRayHandle &left) + : type(Type::Deletion), pos(pos_), idx(INF_IDX), leftRay(left) {} + + Event(const tPoint &pos_, const tRayHandle &left, + const tRayHandle &right) + : type(Type::Intersection), pos(pos_), idx(INF_IDX), leftRay(left), + rightRay(right) {} + + Type type; + tPoint pos; + + tIndex idx; + + tRayHandle leftRay; + tRayHandle rightRay; + }; + +public: + static std::string name() { + return "Sweepline"; + } + + auto operator()(const tDim &K, const tPoints &points, const tBox &bounds) const { + tYaoGraph g(points.size(), K); + + auto rays = Kernel::computeCones(K); + + for (tDim k = 0; k < K; ++k) { + sweepline(points, k, K, rays, g, bounds); + } + + return g; + } + +private: + // Only for pairs of std::hash-able types for simplicity. + // You can of course template this struct to allow other hash functions + struct pair_hash { + template + std::size_t operator()(const std::pair &p) const { + auto h1 = std::hash{}(p.first.addr()); + auto h2 = std::hash{}(p.second.addr()); + + // from Boost HashCombine + return h1 ^ (h2 + 0x9e3779b9 + (h1 << 6) + (h2 >> 2)); + } + }; + + struct it_hash { + template + std::size_t operator()(const IT &p) const { + return std::hash{}(p.addr()); + } + }; + + void sweepline(const tPoints &iPoints, tDim k, const tDim &K, const auto &rays, tYaoGraph &graph, + const tBox &iBounds) const { + + using pqItem = std::pair; + auto pqCmp = [](const pqItem &a, const pqItem &b) { + return a.first > b.first; + }; + + using PQ = PriQueueAdapter; + using isKey = std::pair; + robin_hood::unordered_map isMap; + robin_hood::unordered_map extMap; + + auto bounds = Kernel::mkBBox(iBounds); + + tDirection lRay(-rays[k]); // range [0..2PI] + tDirection rRay(-rays[(k + 1) % K]);// range [0..2PI] + //ASSERT(lRay.angle() <= rRay.angle()); + + tDirection slDir(Kernel::Bisector(lRay, rRay));// range [0..2PI] + //ASSERT(lRay.angle() <= slDir.angle() && slDir.angle() <= rRay.angle()); + + SweeplineDS sl(slDir);// range [0..2PI] + + typename PQ::arrType pqArr; + pqArr.reserve(iPoints.size()); + for (tIndex i = 0; i < iPoints.size(); ++i) { + auto p = Kernel::mkPoint(iPoints[i]); + pqArr.push_back({sl.prj(p), Event(p, i)}); + } + + PQ pq(std::move(pqArr), std::sqrt(iPoints.size())); + +#ifdef WITH_CAIRO + + [[maybe_unused]] constexpr std::array BLACK = {0, 0, 0}; + [[maybe_unused]] constexpr std::array BLUE = {0, 0, 1}; + [[maybe_unused]] constexpr std::array RED = {1, 0, 0}; + [[maybe_unused]] constexpr std::array GREEN = {0, 1, 0}; + [[maybe_unused]] constexpr std::array CYAN = {0, 1, 1}; + [[maybe_unused]] constexpr std::array YELLOW = {1, 1, 0}; + [[maybe_unused]] constexpr std::array ORANGE = {.8, .4, 0}; + [[maybe_unused]] constexpr std::array PINK = {.7, .3, .7}; + + Painter basePainter(iBounds); + basePainter.setColor(BLACK); + basePainter.draw(iPoints, true); +#endif + + LOG("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl); + LOG("Processing cone " << k << " (" << lRay.angle() << ", " << rRay.angle() << ") - sweepline " << slDir.angle() << std::endl); + + tIndex idx = 0; +#ifdef WITH_STATS + tIndex ipProcessed = 0; + tIndex isProcessed = 0; + tIndex delProcessed = 0; + + tIndex maxIpInQueue = iPoints.size(); + tIndex maxIsInQueue = 0; + tIndex maxDelInQueue = 0; + tIndex maxSlSize = 0; +#endif + + while (!pq.empty()) { + +#ifdef WITH_STATS + if (maxIsInQueue < isMap.size()) maxIsInQueue = isMap.size(); + if (maxDelInQueue < extMap.size()) maxDelInQueue = extMap.size(); + if (maxSlSize < sl.size()) maxSlSize = sl.size(); +#endif + + auto cKey = pq.top().first; + auto cPoint = pq.top().second; + // pq.pop(); // we need the pq for drawing + + LOG(idx << ": " + << cPoint.idx << " (" << cPoint.pos << "): key: " << cKey << std::endl); + +#if defined(WITH_CAIRO) && defined(PAINT_STEPS) + //basePainter.draw(cPoint.pos, idx); + + Painter stepPainter(basePainter); + + stepPainter.setColor(RED); + if (cPoint.type == Event::Type::Input) { + stepPainter.draw(cPoint.pos, idx, false); + } else if (cPoint.type == Event::Type::Intersection) { + stepPainter.drawSquare(cPoint.pos); + } else if (cPoint.type == Event::Type::Deletion) { + stepPainter.drawTri(cPoint.pos); + } + + stepPainter.setColor(BLUE); + sl.draw(cPoint.pos, stepPainter); + + stepPainter.setColor(ORANGE); + for (const auto &is : isMap) { + stepPainter.drawSquare(pq.get_key(is.second).second.pos); + } + + stepPainter.setColor(PINK); + for (const auto &ext : extMap) { + stepPainter.drawTri(pq.get_key(ext.second).second.pos); + } + + std::ostringstream stepFilename; + stepFilename << "img_k" << k << "_s" << std::setfill('0') + << std::setw(static_cast(std::ceil(std::log10(iPoints.size()) + 2))) << idx; + stepPainter.save(stepFilename.str()); +#endif + // now pop PQ after drawing + pq.pop(); + + auto handleRayIntersection = [&isMap, &extMap, &pq, &bounds, &sl, &idx, &cKey](auto &leftRay, auto &rightRay, bool left) { + auto is = leftRay->intersection(*rightRay, bounds); + if (is) { + const tPoint *p = std::get_if(&*is); + if (p && sl.prj(*p) > cKey) { + isMap[std::make_pair(leftRay, rightRay)] = pq.push({sl.prj(*p), Event({*p, leftRay, rightRay})}); + LOG(idx << ": " + << " added " << (left ? "left" : "right") + << " intersection point at (" << *p << ") key: " << sl.prj(*p) << std::endl); + } + tRay *r = std::get_if(&*is); + if (r) { + if (r->direction() != leftRay->direction()) { + r->invert(); + } + + LOG(idx << ": " + << "RI left: " << *leftRay << std::endl); + LOG(idx << ": " + << "RI right: " << *rightRay << std::endl); + LOG(idx << ": " + << "RI IS Ray: " << *r << std::endl); + + ASSERT(leftRay->rightRegion == rightRay->leftRegion); + ASSERT(r->direction() == leftRay->direction() && r->direction() == rightRay->direction()); + + if (left) { + if (rightRay->isExtended()) { + auto pqR = extMap.find(rightRay); + ASSERT(pqR != extMap.end()); + pq.remove(pqR->second); + LOG(idx << ": " + << " RI deleted Deletion event " << rightRay->extOrigin() << " key: " << sl.prj(rightRay->extOrigin()) << std::endl); + extMap.erase(pqR); + } + + // remove rightRay, as leftRay replaces it + // point iterator to right ray to left ray + + leftRay->leftRegion = leftRay->leftRegion; + leftRay->rightRegion = rightRay->rightRegion; + sl.erase(rightRay); + rightRay = leftRay; + } else { + if (leftRay->isExtended()) { + auto pqL = extMap.find(leftRay); + ASSERT(pqL != extMap.end()); + pq.remove(pqL->second); + LOG(idx << ": " + << " RI deleted Deletion event " << leftRay->extOrigin() << " key: " << sl.prj(leftRay->extOrigin()) << std::endl); + extMap.erase(pqL); + } + + // remove leftRay, as rightRay replaces it + // point iterator to left ray to right ray + + rightRay->leftRegion = leftRay->leftRegion; + rightRay->rightRegion = rightRay->rightRegion; + sl.erase(leftRay); + leftRay = rightRay; + } + + LOG(idx << ": " + << "new Ray: " << (left ? *leftRay : *rightRay) << std::endl); + } + } + }; + + switch (cPoint.type) { + case Event::Type::Input: { + + LOG(idx << ": " + << " type: input point" << std::endl); + + auto itBr = sl.find(cPoint.pos);// right ray + + if (itBr != sl.end() && itBr->orientedSide(cPoint.pos) == tOrientedSide::LINE) {// check whether point is ON right boundary + if (itBr->direction() == lRay) { + LOG(idx << ": point ON right boundary that is a left ray" << std::endl); + LOG(idx << ": old Ray " << *itBr << std::endl); + + itBr = std::next(itBr); + + LOG(idx << ": new Ray " << *itBr << std::endl); + } + } + + auto itBl = (itBr == sl.begin() ? sl.end() : std::prev(itBr));// left ray + ASSERT(itBl == sl.end() || itBl->orientedSide(cPoint.pos) != tOrientedSide::RIGHT); + + LOG(idx << ": " + << " left ray: " << (itBl != sl.end() ? to_string(*itBl) : "NULL") << std::endl); + LOG(idx << ": " + << " right ray: " << (itBr != sl.end() ? to_string(*itBr) : "NULL") << std::endl); + +#if defined(WITH_CAIRO) && defined(PAINT_STEPS) + stepPainter.setColor(CYAN); + if (itBl != sl.end()) itBl->draw(stepPainter); + stepPainter.setColor(GREEN); + if (itBr != sl.end()) itBr->draw(stepPainter); + stepPainter.save(stepFilename.str()); +#endif + + // add graph edge + // if p lies on right boundary, it does not belong to this cone + if (itBr != sl.end() && itBr->leftRegion != INF_IDX) { + ASSERT(itBl->rightRegion == itBr->leftRegion); + + auto r = itBr->orientedSide(cPoint.pos) != tOrientedSide::LINE ? itBr->leftRegion : itBr->rightRegion; + if (r != INF_IDX) { + graph[cPoint.idx].neighbor[k] = r; + graph[cPoint.idx].distance[k] = Kernel::to_float(Kernel::distance2(cPoint.pos, Kernel::mkPoint(iPoints[r]))); + } + + LOG(idx << ": " + << " edge added: (" << cPoint.idx << ", " << itBr->leftRegion << ") w: " << Kernel::distance2(cPoint.pos, Kernel::mkPoint(iPoints[itBr->leftRegion])) << std::endl); + +#ifdef WITH_CAIRO + basePainter.setColor(BLACK); + basePainter.drawLine(Kernel::mkPoint(iPoints[itBr->leftRegion]), cPoint.pos); +#endif + } + + // check if Bl and Br intersect, check only required if they don't originate from same point + if (itBl != sl.end() && itBr != sl.end()) { + auto itIs = isMap.find(std::make_pair(itBl, itBr)); + if (itIs != isMap.end()) { + LOG(idx << ": " + << " deleted intersection point between " << *itIs->first.first << " and " << *itIs->first.second << std::endl); + pq.remove(itIs->second); + isMap.erase(itIs); + } + } + + // create new rays and insert them into SL + auto itBln = sl.insert_pair( + itBr, + tRay({cPoint.pos, lRay, + itBl != sl.end() ? itBl->rightRegion : INF_IDX, cPoint.idx}), + tRay({cPoint.pos, rRay, cPoint.idx, + itBr != sl.end() ? itBr->leftRegion : INF_IDX})); + auto itBrn = std::next(itBln); + + ASSERT(itBln == std::prev(itBrn)); + ASSERT(itBrn == std::next(itBln)); + + LOG(idx << ": " + << " left ray " << *itBln << std::endl); + LOG(idx << ": " + << " right ray " << *itBrn << std::endl); + + // insert intersection points into PQ + if (itBl != sl.end()) { + handleRayIntersection(itBl, itBln, true); + } + + if (itBr != sl.end()) { + handleRayIntersection(itBrn, itBr, false); + } + +#ifdef WITH_STATS + ipProcessed++; +#endif + break; + } + + case Event::Type::Intersection: { + + LOG(idx << ": " + << " type: intersection point" << std::endl); + + auto itBl = cPoint.leftRay; + auto itBr = cPoint.rightRay; + +#if defined(WITH_CAIRO) && defined(PAINT_STEPS) + stepPainter.setColor(CYAN); + itBl->draw(stepPainter); + itBr->draw(stepPainter); + stepPainter.save(stepFilename.str()); +#endif + + ASSERT(std::next(itBl) == itBr); + ASSERT(itBl == std::prev(itBr)); + + LOG(idx << ": " + << " left ray: " << (itBl != sl.end() ? to_string(*itBl) : "NULL") << std::endl); + LOG(idx << ": " + << " right ray: " << (itBr != sl.end() ? to_string(*itBr) : "NULL") << std::endl); + + // delete intersection point from hash map + [[maybe_unused]] auto chk = isMap.erase(std::make_pair(itBl, itBr)); + ASSERT(chk); + + // delete intersection points from PQ + if (itBl != sl.begin()) { + auto itBll = std::prev(itBl); + ASSERT(itBll != sl.end()); + + auto itIs = isMap.find(std::make_pair(itBll, itBl)); + if (itIs != isMap.end()) { + LOG(idx << ": " + << " deleted intersection point between " << *itIs->first.first << " and " << *itIs->first.second << std::endl); + pq.remove(itIs->second); + isMap.erase(itIs); + } + } + + if (itBr != sl.end()) { + auto itBrr = std::next(itBr); + + if (itBrr != sl.end()) { + auto itIs = isMap.find(std::make_pair(itBr, itBrr)); + if (itIs != isMap.end()) { + LOG(idx << ": " + << " deleted intersection point between " << *itIs->first.first << " and " << *itIs->first.second << std::endl); + pq.remove(itIs->second); + isMap.erase(itIs); + } + } + } + + tRay rL({cPoint.pos, lRay, itBl->leftRegion, itBr->rightRegion}); + tRay rR({cPoint.pos, rRay, itBl->leftRegion, itBr->rightRegion}); + + // bisector line between A and B + auto pL = Kernel::mkPoint(iPoints[itBl->leftRegion]); + auto pR = Kernel::mkPoint(iPoints[itBr->rightRegion]); + + auto pMid = Kernel::Midpoint(pL, pR); + auto aBs = Kernel::Bisector(pL, pR, sl.slDirection); + + // ASSERT(Kernel::approxEQ(pMid, Kernel::Midpoint(pR, pL))); + // ASSERT(Kernel::approxEQ(aBs.angle(), Kernel::Bisector(pR, pL, sl.slDirection).angle())); + + tRay Bs(pMid, aBs, itBl->leftRegion, itBr->rightRegion); + +#if defined(WITH_CAIRO) && defined(PAINT_STEPS) + stepPainter.setColor(GREEN); + rL.draw(stepPainter); + rR.draw(stepPainter); + + stepPainter.setColor(YELLOW); + stepPainter.drawLine(pL, pR); + stepPainter.drawSquare(pMid); + Bs.draw(stepPainter); +#endif + + LOG(idx << ": " + << " left ray: " << rL << std::endl); + LOG(idx << ": " + << " right ray: " << rR << std::endl); + LOG(idx << ": " + << " bisector: " << Bs << std::endl); + + // check for intersections of Bln, Blr and BS + auto BsL = Bs.intersection(rL, bounds); + const tPoint *pBsL = nullptr; + if (BsL) { + ASSERT(std::holds_alternative(*BsL)); + pBsL = std::get_if(&*BsL); + } + + auto BsR = Bs.intersection(rR, bounds); + const tPoint *pBsR = nullptr; + if (BsR) { + ASSERT(std::holds_alternative(*BsR)); + pBsR = std::get_if(&*BsR); + } + + if (pBsL && pBsR && Kernel::approxEQ(*pBsL, *pBsR) && Kernel::approxEQ(*pBsL, cPoint.pos)) {//TODO better way to test? + ASSERT(Kernel::approxEQ(*pBsR, cPoint.pos)); + // we don't check for sweepline projection to avoid floating point problems + } else { + if (pBsL && Kernel::approxLT(sl.prj(*pBsL), cKey)) {// check whether IS is before SL + pBsL = nullptr; + } + if (pBsR && Kernel::approxLT(sl.prj(*pBsR), cKey)) {// check whether IS is before SL + pBsR = nullptr; + } + } + +#if defined(WITH_CAIRO) && defined(PAINT_STEPS) + if (pBsL) { + stepPainter.setColor(RED); + stepPainter.drawSquare(*pBsL); + } + + if (pBsR) { + stepPainter.setColor(RED); + stepPainter.drawSquare(*pBsR); + } +#endif + + // check whether rays have extension before deletion: delete delete event + if (itBl->isExtended()) { + auto pqBl = extMap.find(itBl); + ASSERT(pqBl != extMap.end()); + pq.remove(pqBl->second); + LOG(idx << ": " + << " deleted Deletion event " << itBl->extOrigin() << " key: " << sl.prj(itBl->extOrigin()) << std::endl); + extMap.erase(pqBl); + } + if (itBr->isExtended()) { + auto pqBr = extMap.find(itBr); + ASSERT(pqBr != extMap.end()); + pq.remove(pqBr->second); + LOG(idx << ": " + << " deleted Deletion event " << itBr->extOrigin() << " key: " << sl.prj(itBr->extOrigin()) << std::endl); + extMap.erase(pqBr); + } + +#ifdef WITH_CAIRO + basePainter.setColor(BLACK); + basePainter.setDash(); + basePainter.drawSquare(cPoint.pos); + basePainter.drawLine(itBl->origin(), cPoint.pos); + basePainter.drawLine(itBr->origin(), cPoint.pos); + basePainter.unsetDash(); +#endif + + // remove old Rays from list + sl.erase(itBl); + //auto itInsertPos = sl.erase(itBr); + auto itBn = sl.end(); + + if (!pBsL && !pBsR) { + LOG(idx << ": " + << " case a) no intersection" << std::endl); + // bisector intersects no ray from P + if (sl.prj(pR) < sl.prj(pL)) { + // right point is further from v -> right ray becomes border + itBn = sl.replace(itBr, rR); + } else { + // left point is further from v -> left ray becomes border + itBn = sl.replace(itBr, rL); + } + } else { + if (pBsL && pBsR) { + // bisector intersects both rays - > must be in same point v + ASSERT(Kernel::approxEQ(*pBsL, *pBsR)); + ASSERT(Kernel::approxEQ(*pBsL, cPoint.pos)); + ASSERT(Kernel::approxEQ(*pBsR, cPoint.pos)); + LOG(idx << ": " + << " case c) both intersect" << std::endl); + // boundary originates at v with bisector angle + Bs.setOrigin(cPoint.pos); + itBn = sl.replace(itBr, Bs); + } else { + LOG(idx << ": " + << " case b) one intersection" << std::endl); + if (pBsL) { + // bisector intersects left ray + // boundary to intersection point, then ray with bisector angle + Bs.setOrigin(*pBsL); + if (!Kernel::approxEQ(rL.origin(), Bs.origin())) { + itBn = sl.replace(itBr, tRay({rL, Bs})); + extMap[itBn] = pq.push({sl.prj(Bs.origin()), Event(Bs.origin(), itBn)}); + LOG(idx << ": " + << " added deletion point at (" << Bs.origin() << ") key: " << sl.prj(Bs.origin()) + << " for " << *itBn << std::endl); + } else {// if they are almost equal immediately use bisector ray + itBn = sl.replace(itBr, Bs); + } + } else { + // bisector intersects right ray + // boundary to intersection point, then ray with bisector angle + Bs.setOrigin(*pBsR); + if (!Kernel::approxEQ(rR.origin(), Bs.origin())) { + itBn = sl.replace(itBr, tRay({rR, Bs})); + extMap[itBn] = pq.push({sl.prj(Bs.origin()), Event(Bs.origin(), itBn)}); + LOG(idx << ": " + << " added deletion point at (" << Bs.origin() << ") key: " << sl.prj(Bs.origin()) + << " for " << *itBn << std::endl); + } else {// if they are almost equal immediately use bisector ray + itBn = sl.replace(itBr, Bs); + } + } + } + } + ASSERT(itBn != sl.end());// some boundary was inserted into SL + + LOG(idx << ": " + << "found boundary: " << *itBn << std::endl); + + // insert intersection points into PQ + + auto saveItBn = itBn; + auto saveItL = itBn; + if (itBn != sl.begin()) { + auto itL = std::prev(itBn); + saveItL = itL; + handleRayIntersection(itL, itBn, true); + } + + auto itR = std::next(itBn); + if (itR != sl.end()) { + handleRayIntersection(itBn, itR, false); + } + + //check if itBn was modified and re-route intersection points + if (itBn != saveItBn && saveItL != saveItBn) { + // saveItL is different from saveItBn so there was a previous IT + + // the left neighbor should still be the same with the new itBn + ASSERT(saveItL == std::prev(itBn)); + + auto itIs = isMap.find(std::make_pair(saveItL, saveItBn)); + if (itIs != isMap.end()) { + auto oldEvent = pq.get_key(itIs->second); + pq.remove(itIs->second); + isMap.erase(itIs); + + isMap[std::make_pair(saveItL, itBn)] = pq.push({sl.prj(oldEvent.second.pos), + Event({oldEvent.second.pos, saveItL, itBn})}); + + LOG(idx << ": " + << " updated intersection point between " << *saveItL << " and " << *itBn << std::endl); + } + } + +#if defined(WITH_CAIRO) && defined(PAINT_STEPS) + stepPainter.save(stepFilename.str()); +#endif + +#ifdef WITH_STATS + isProcessed++; +#endif + break; + } + + case Event::Type::Deletion: { + + LOG(idx << ": " + << " type: deletion point" << std::endl); + + auto itB = cPoint.leftRay;// we store the ray to be deleted as left ray + ASSERT(itB != sl.end()); + ASSERT(itB->isExtended()); + //ASSERT(itB->leftRegion == itB->ext->leftRegion); + //ASSERT(itB->rightRegion == itB->ext->rightRegion); + //ASSERT(!itB->ext->ext); + + // delete extension point from hash map + [[maybe_unused]] auto chk = extMap.erase(itB); + ASSERT(chk); + +#ifdef WITH_CAIRO + basePainter.setColor(BLACK); + basePainter.setDash(); + basePainter.drawLine(itB->origin(), itB->extOrigin()); + basePainter.drawTri(cPoint.pos); + basePainter.unsetDash(); +#endif + + LOG(idx << ": " + << " old ray: " << *itB << std::endl); + + // we replace the RayUnion with its lower ray, so all intersections pointers should still be valid + // all intersections processed after this point will be with lower ray + itB->foldExtension(); + + LOG(idx << " new ray: " << *itB << std::endl); + + ASSERT(!itB->isExtended()); +#ifdef WITH_STATS + delProcessed++; +#endif + break; + } + } + + ++idx; + + LOG(std::endl); + } + +#ifdef WITH_STATS + std::cout << "STATS" + << " algorithm=" << name() + << " kernel=" << Kernel::name() + << " dist=" << iPoints.getDistName() + << " n=" << iPoints.size() + << " k=" << k + << " steps=" << idx + << " ipPro=" << ipProcessed + << " isPro=" << isProcessed + << " delPro=" << delProcessed + << " maxIpQ=" << maxIpInQueue + << " maxIsQ=" << maxIsInQueue + << " maxDelQ=" << maxDelInQueue + << " maxSlSize=" << maxSlSize + << std::endl; +#endif + +#ifdef WITH_CAIRO + basePainter.save("img_k" + std::to_string(k)); +#endif + } +}; diff --git a/include/Generators/Bubbles.hpp b/include/Generators/Bubbles.hpp new file mode 100644 index 0000000..e8776a8 --- /dev/null +++ b/include/Generators/Bubbles.hpp @@ -0,0 +1,76 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include "GeneratorBase.hpp" + +class Bubbles : public Generator { + + friend class Generator; + +public: + Bubbles(const tIndex &seed) : Generator(seed), + dist(0, std::nextafter(1, std::numeric_limits::max())) {} + + std::string name() const override { + return "bubbles"; + } + + std::tuple _generate(const tIndex n, const tBox &bounds) override { + tPoints points; + points.reserve(n); + + auto gMidpoint = bounds.low; + + for (uint d = 0; d < D; ++d) { + gMidpoint[d] += (bounds.high[d] - bounds.low[d]) / 2; + } + + std::vector bubbles; + bubbles.resize(pow(2, D)); + + for (uint i = 0; i < pow(2, D); ++i) { + for (uint d = 0; d < D; ++d) { + bubbles[i].low[d] = + i & (1 << d) ? gMidpoint[d] : bounds.low[d]; + bubbles[i].high[d] = + i & (1 << d) ? bounds.high[d] : gMidpoint[d]; + } + } + + tIndex ppB = n / bubbles.size();// points per bubble + for (uint b = 0; b < bubbles.size(); ++b) { + const auto &bubble = bubbles[b]; + + auto midpoint = bubble.low; + auto span = bubble.high; + + for (uint d = 0; d < D; ++d) { + midpoint[d] += (bubble.high[d] - bubble.low[d]) / 2; + span[d] -= bubble.low[d]; + } + + for (tIndex i = b * ppB + 1; i <= (b + 1) * ppB; ++i) { + tIFloatVector p; + + do { + p = midpoint; + + for (uint d = 0; d < D; ++d) { + p[d] += span[d] * rand(); + } + } while (!bounds.contains(p)); + + + points.push_back(p); + } + } + + return std::make_tuple(points, bounds); + } + +private: + std::uniform_real_distribution dist; +}; \ No newline at end of file diff --git a/include/Generators/Circle.hpp b/include/Generators/Circle.hpp new file mode 100644 index 0000000..f326d7b --- /dev/null +++ b/include/Generators/Circle.hpp @@ -0,0 +1,65 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include "GeneratorBase.hpp" + +class Circle : public Generator { + + friend class Generator; + +public: + Circle(const tIndex &seed) : Generator(seed), + dist(0, std::nextafter(1, std::numeric_limits::max())) {} + + std::string name() const override { + return "circle"; + } + + std::tuple _generate(const tIndex n, const tBox &bounds) override { + tPoints points; + points.reserve(n); + + tIFloatVector radius; + auto center = bounds.low; + + for (uint d = 0; d < D; ++d) { + center[d] += (bounds.high[d] - bounds.low[d]) / 2; + radius[d] = .45 * (bounds.high[d] - bounds.low[d]);//give it some space at the edge + } + + // we generate normal distributed points and map them to points on a unit circle/sphere + // these points are uniform distributed on the circle/sphere + // we then map the circle/sphere points to the ellipse/ellipsoid + // the resulting points are not uniform distributed, but close enough + + for (tIndex i = 1; i <= n; ++i) { + tIFloatVector p; + + //normal distributed point around origin + tIFloat r = 0;//radius + for (uint d = 0; d < D; ++d) { + p[d] = rand(); + r += p[d] * p[d]; + } + + r = std::sqrt(r); + + // now map them to the ellipse/ellipsoid + for (uint d = 0; d < D; ++d) { + // origin distortion point on circle/sphere + p[d] = center[d] + radius[d] * (p[d] / r); + } + + + points.push_back(p); + } + + return std::make_tuple(points, bounds); + } + +private: + std::normal_distribution dist; +}; \ No newline at end of file diff --git a/include/Generators/Gaussian.hpp b/include/Generators/Gaussian.hpp new file mode 100644 index 0000000..3282c12 --- /dev/null +++ b/include/Generators/Gaussian.hpp @@ -0,0 +1,46 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include "GeneratorBase.hpp" + +class Gaussian : public Generator { + + friend class Generator; + +public: + Gaussian(const tIndex &seed) : Generator(seed), + dist(0, .25) {} + + std::string name() const override { + return "gaussian"; + } + + std::tuple _generate(const tIndex n, const tBox &bounds) override { + tPoints points; + points.reserve(n); + + auto halfLength = .5 * (bounds.high - bounds.low); + auto midpoint = bounds.low + halfLength; + + for (tIndex i = 0; i < n; ++i) { + + tIFloatVector p; + do { + for (tDim d = 0; d < D; ++d) { + p[d] = midpoint[d] + halfLength[d] * rand(); + } + } while (!bounds.contains(p)); + + ASSERT(bounds.contains(p)); + points.emplace_back(p); + } + + return std::make_tuple(points, bounds); + } + +private: + std::normal_distribution dist; +}; \ No newline at end of file diff --git a/include/Generators/GeneratorBase.hpp b/include/Generators/GeneratorBase.hpp new file mode 100644 index 0000000..3346672 --- /dev/null +++ b/include/Generators/GeneratorBase.hpp @@ -0,0 +1,206 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "Predicates.hpp" +#include "Types.hpp" +#include "Utils/ASSERT.hpp" + +class GeneratorBase { +public: + std::tuple generate(const tIndex n, const tBox &bounds) { + auto [p, b] = _generate(n, bounds); + p.setDistName(name()); + p.setSeed(seed()); + + return std::make_tuple(p, b); + } + + virtual ~GeneratorBase() = default; + + virtual std::string name() const = 0; + virtual tIndex seed() const = 0; + virtual void setSeed(const tIndex &seed) = 0; + +protected: + virtual std::tuple _generate(const tIndex n, const tBox &bounds) = 0; + +public: + static tPoints rescalePoints(const tPoints &inPoints, const tBox &oldBounds, const tBox &newBounds) { + tPoints outPoints; + outPoints.reserve(inPoints.size()); + + for (const auto &i : inPoints) { + + tIFloatVector p; + for (tDim d = 0; d < D; ++d) { + p[d] = newBounds.low[d] + (newBounds.high[d] - newBounds.low[d]) * ((i[d] - oldBounds.low[d]) / static_cast((oldBounds.high[d] - oldBounds.low[d]))); + } + + ASSERT(newBounds.contains(p)); + outPoints.emplace_back(p); + } + + return outPoints; + } +}; + +template +class Generator : public GeneratorBase { + +public: + Generator(const tIndex &seed) : gen(seed), _seed(seed) {} + +protected: + auto rand() { + return static_cast(*this).dist(gen); + } + + // See https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle + std::vector FisherYatesShuffle(const tIndex &sample_size, const tIndex &pop_size) { + if (sample_size > pop_size) { + std::vector sample(pop_size); + std::iota(sample.begin(), sample.end(), 0); + + return sample; + } + + std::vector sample(sample_size); + + for (tIndex i = 0; i != pop_size; ++i) { + auto lDist = std::uniform_int_distribution(0, i); + tIndex j = lDist(gen); + if (j < sample.size()) { + if (i < sample.size()) { + sample[i] = sample[j]; + } + sample[j] = i; + } + } + return sample; + } + + template + std::vector sort_indices(const std::vector &v, Compare comp) { + + // initialize original index locations + std::vector idx(v.size()); + std::iota(idx.begin(), idx.end(), 0); + + // sort indexes based on comparing values in v + // using std::stable_sort instead of std::sort + // to avoid unnecessary index re-orderings + // when v contains elements of equal values + std::stable_sort(idx.begin(), idx.end(), + [&v, &comp](tIndex i1, tIndex i2) { return comp(v[i1], v[i2]); }); + + return idx; + } + + template + std::tuple extract_points(const tIndex &n, const std::vector> &inPoints, + std::array, D2> coordSort) { + ASSERT(inPoints.size() > n); + + // sample random point to start with + std::array seedPoint; + for (tDim d = 0; d < D2; ++d) { + seedPoint[d] = static_cast(std::floor((.25 + .5 * rand()) * coordSort[d].size())); + } + + tIFloat windowSize = .0; + tIFloat windowSizeInc = .01; + const tIFloat maxWindowSize = .5; + + std::vector pointIdx; + + while (pointIdx.size() < n && windowSize <= maxWindowSize) { + pointIdx.clear();//TODO: this can be done smarter + + windowSize += windowSizeInc; + + std::array, D2> windowIdx; + + for (tDim d = 0; d < D2; ++d) { + + auto begin = coordSort[d].begin(); + std::advance(begin, + static_cast(std::floor(std::max(0.0, + seedPoint[X] - (coordSort[d].size() * windowSize))))); + + auto end = coordSort[d].begin(); + std::advance(end, + static_cast(std::floor(std::min(static_cast(coordSort[d].size()), + seedPoint[X] + (coordSort[d].size() * windowSize))))); + + std::copy(begin, end, std::back_inserter(windowIdx[d])); + std::sort(windowIdx[d].begin(), windowIdx[d].end()); + + if (windowIdx[d].size() < n) { + continue; + } + } + + std::array, D2> windowIdxIS; + windowIdxIS[0] = windowIdx[0]; + + for (tDim d = 1; d < D2; ++d) { + std::set_intersection(windowIdxIS[d - 1].begin(), windowIdxIS[d - 1].end(), windowIdx[d].begin(), windowIdx[d].end(), std::back_inserter(windowIdxIS[d])); + } + + pointIdx = windowIdxIS[D2 - 1]; + + if (pointIdx.size() > 1.05 * n) { + windowSize -= windowSizeInc; + windowSizeInc /= 2; + pointIdx.clear(); + } + } + + tPoints points; + points.reserve(pointIdx.size()); + + tIFloatVector minPoint = {std::numeric_limits::max(), std::numeric_limits::max()}; + tIFloatVector maxPoint = {std::numeric_limits::min(), std::numeric_limits::min()}; + + for (const auto &i : pointIdx) { + tIFloatVector p; + for (tDim d = 0; d < D; ++d) { + p[d] = inPoints[i][d]; + } + + if (p[X] < minPoint[X]) minPoint[X] = p[X]; + if (p[Y] < minPoint[Y]) minPoint[Y] = p[Y]; + if (p[X] > maxPoint[X]) maxPoint[X] = p[X]; + if (p[Y] > maxPoint[Y]) maxPoint[Y] = p[Y]; + + points.emplace_back(p); + } + + return std::make_tuple(points, tBox{minPoint, maxPoint}); + } + +public: + tIndex seed() const override { + return _seed; + } + + void setSeed(const tIndex &seed) override { + _seed = seed; + gen.seed(seed); + } + + +private: + std::mt19937 gen; + tIndex _seed; +}; diff --git a/include/Generators/Generators.hpp b/include/Generators/Generators.hpp new file mode 100644 index 0000000..6c712a6 --- /dev/null +++ b/include/Generators/Generators.hpp @@ -0,0 +1,63 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include + +#include "Bubbles.hpp" +#include "Circle.hpp" +#include "Gaussian.hpp" +#include "Grid.hpp" +#include "Road.hpp" +#include "Stars.hpp" +#include "Uniform.hpp" + +std::unique_ptr getGen(const char &dist, const tIndex &seed) { + // [_u_ni, _g_aussian, gri_d_, _r_oad, _s_tar] + + switch (dist) { + case 'u': + return std::make_unique(seed); + case 'g': + return std::make_unique(seed); + case 'd': + return std::make_unique(seed); + case 'r': + return std::make_unique(seed); + case 's': + return std::make_unique(seed); + case 'b': + return std::make_unique(seed); + case 'c': + return std::make_unique(seed); + + default: + return nullptr; + } +} + +std::string getGenName(const char &dist) { + // [_u_ni, _g_aussian, gri_d_, _r_oad, _s_tar] + + switch (dist) { + case 'u': + return "uni"; + case 'g': + return "gaussian"; + case 'd': + return "grid"; + case 'r': + return "road"; + case 's': + return "stars"; + case 'c': + return "circle"; + case 'b': + return "bubbles"; + + default: + return nullptr; + } +} \ No newline at end of file diff --git a/include/Generators/Grid.hpp b/include/Generators/Grid.hpp new file mode 100644 index 0000000..02c4ecd --- /dev/null +++ b/include/Generators/Grid.hpp @@ -0,0 +1,44 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include "GeneratorBase.hpp" + +class Grid : public Generator { + + friend class Generator; + +public: + Grid(const tIndex &seed) : Generator(seed) {} + + std::string name() const override { + return "grid"; + } + + std::tuple _generate(const tIndex n, const tBox &bounds) override { + tPoints points; + points.reserve(n); + + tIndex pointsPerDim = std::floor(std::sqrt(n)); + + tIFloatVector pointSpacing; + for (tDim d = 0; d < bounds.high.size(); ++d) { + pointSpacing[d] = (bounds.high[d] - bounds.low[d]) / pointsPerDim; + } + + for (tIndex i = 0; i < pointsPerDim; ++i) { + for (tIndex j = 0; j < pointsPerDim; ++j) { + tIFloatVector p; + p[X] = bounds.low[X] + (pointSpacing[X] / 2) + pointSpacing[X] * i; + p[Y] = bounds.low[Y] + (pointSpacing[Y] / 2) + pointSpacing[Y] * j; + + ASSERT(bounds.contains(p)); + points.emplace_back(p); + } + } + + return std::make_tuple(points, bounds); + } +}; diff --git a/include/Generators/Road.hpp b/include/Generators/Road.hpp new file mode 100644 index 0000000..8969cc8 --- /dev/null +++ b/include/Generators/Road.hpp @@ -0,0 +1,66 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include "GeneratorBase.hpp" + +class Road : public Generator { + + friend class Generator; + +public: + Road(const tIndex &seed) : Generator(seed) { + + //std::string fileName = ROAD_DATA "/USA-road-d.USA.co"; + std::string fileName = ROAD_DATA "/USA-road-d.NY.co"; + std::ifstream file(fileName, std::ios::in); + std::string line; + + while (std::getline(file, line)) { + if (line.starts_with('c')) { + // comment line + continue; + } + if (line.starts_with('p')) { + // problem line, occurs before any input points by definition + // p aux co N + auto lastSpace = line.find_last_of(' '); + tIndex pointsInFile = std::stoul(line.substr(lastSpace + 1)); + filePoints.resize(pointsInFile); + } + if (line.starts_with('v')) { + // vertex line + // v IDX X Y + + std::istringstream iss(line); + char c; + int idx, x, y; + if (!(iss >> c >> idx >> x >> y)) { continue; }// error in file + + filePoints[idx - 1] = {x / 1e6, y / 1e6}; + } + } + + for (tDim d = 0; d < 2; ++d) { + coordSort[d] = sort_indices(filePoints, [d](const tFilePoint &a, const tFilePoint &b) { return a[d] < b[d]; }); + } + } + + std::string name() const override { + return "road"; + } + + std::tuple _generate(const tIndex n, [[maybe_unused]] const tBox &bounds) override { + + return extract_points(n, filePoints, coordSort); + } + +private: + std::uniform_real_distribution dist; + + using tFilePoint = std::array; + std::vector filePoints; + std::array, 2> coordSort; +}; diff --git a/include/Generators/Stars.hpp b/include/Generators/Stars.hpp new file mode 100644 index 0000000..18657cc --- /dev/null +++ b/include/Generators/Stars.hpp @@ -0,0 +1,52 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include "GeneratorBase.hpp" + +class Stars : public Generator { + + friend class Generator; + +public: + Stars(const tIndex &seed) : Generator(seed) { + + //std::string fileName = STAR_DATA "/gaia_1331909727.points"; + std::string fileName = STAR_DATA "/gaia_1.points"; + std::ifstream file(fileName, std::ios::in); + std::string line; + + while (std::getline(file, line)) { + // star line + // x y z + + std::istringstream iss(line); + double x, y, z; + if (!(iss >> x >> y >> z)) { continue; }// error in file + + filePoints.push_back({x, y, z}); + } + + for (tDim d = 0; d < 3; ++d) { + coordSort[d] = sort_indices(filePoints, [d](const tFilePoint &a, const tFilePoint &b) { return a[d] < b[d]; }); + } + } + + std::string name() const override { + return "stars"; + } + + std::tuple _generate(const tIndex n, [[maybe_unused]] const tBox &bounds) override { + + return extract_points(n, filePoints, coordSort); + } + +private: + std::uniform_real_distribution dist; + + using tFilePoint = std::array; + std::vector filePoints; + std::array, 3> coordSort; +}; \ No newline at end of file diff --git a/include/Generators/Uniform.hpp b/include/Generators/Uniform.hpp new file mode 100644 index 0000000..1cabb91 --- /dev/null +++ b/include/Generators/Uniform.hpp @@ -0,0 +1,41 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include "GeneratorBase.hpp" + +class Uniform : public Generator { + + friend class Generator; + +public: + Uniform(const tIndex &seed) : Generator(seed), + dist(0, std::nextafter(1, std::numeric_limits::max())) {} + + std::string name() const override { + return "uni"; + } + + std::tuple _generate(const tIndex n, const tBox &bounds) override { + tPoints points; + points.reserve(n); + + for (tIndex i = 0; i < n; ++i) { + + tIFloatVector p; + for (tDim d = 0; d < D; ++d) { + p[d] = bounds.low[d] + (bounds.high[d] - bounds.low[d]) * rand(); + } + + ASSERT(bounds.contains(p)); + points.emplace_back(p); + } + + return std::make_tuple(points, bounds); + } + +private: + std::uniform_real_distribution dist; +}; \ No newline at end of file diff --git a/include/Kernels/CGALKernel.hpp b/include/Kernels/CGALKernel.hpp new file mode 100644 index 0000000..a60aded --- /dev/null +++ b/include/Kernels/CGALKernel.hpp @@ -0,0 +1,523 @@ +// +// Created by funke on 10/1/21. +// + +#pragma once + +#include +#include +#include +#include + +#include "Predicates.hpp" +#include "Types.hpp" + +#ifdef WITH_CAIRO +#include "Painter.hpp" +#endif + +using ExactPredicatesInexactConstructions = CGAL::Exact_predicates_inexact_constructions_kernel; +using ExactPredicatesExactConstructions = CGAL::Exact_predicates_exact_constructions_kernel; + +template +struct KernelName { + static std::string name() { + return K::name(); + } +}; + +template<> +struct KernelName { + static std::string name() { + return "CGALExactPredExactCon"; + } +}; + +template<> +struct KernelName { + static std::string name() { + return "CGALExactPredInexactCon"; + } +}; + +template +class CGALKernel { + +public: + using Float = typename K::FT; + using Vector = CGAL::Vector_2; + using Point = typename K::Point_2; + using Segment = CGAL::Segment_2; + using Box = CGAL::Bbox_2; + using Line = CGAL::Line_2; + + static std::string name() { + return KernelName::name(); + } + + class Direction { + public: + Direction(const tIFloat _dir) : dir(Vector(std::cos(_dir), std::sin(_dir))) {} + Direction(const Vector &v) : dir(v) {} + Direction(const CGAL::Direction_2 &d) : dir(d) {} + + Float prj(const Point &p) { + return CGAL::scalar_product(Vector(p[X], p[Y]), dir.vector()) / CGAL::scalar_product(dir.vector(), dir.vector()); + } + + tIFloat angle() const { + return atan2P(CGAL::to_double(dir.dy()), CGAL::to_double(dir.dx())); + } + + auto &base() const { + return dir; + } + + Direction operator-() const { + return Direction(-dir); + } + + Direction perp() const { + return {dir.perpendicular(CGAL::CLOCKWISE)};// TODO: check angle orientation + } + + Direction perp(const Direction &ref) const { + auto p = perp().dir.vector(); + if (CGAL::angle(p, ref.dir.vector()) == CGAL::OBTUSE) { + p *= -1; + } + return Direction(p); + } + + bool operator==(const Direction &o) const { + return dir == o.dir; + } + + private: + CGAL::Direction_2 dir; + }; + + struct Ray { + private: + CGAL::Ray_2 iRay; + + + public: + tIndex leftRegion; + tIndex rightRegion; + + private: + // extension of ray with a segment before ray + std::unique_ptr ext;// initialized with nullptr + + public: + Ray(const Point &p_, const Direction &angle_, + const tIndex &lr, const tIndex &rr) + : iRay(p_, angle_.base()), leftRegion(lr), rightRegion(rr) {} + + Ray(const Ray &o) + : iRay(o.iRay), leftRegion(o.leftRegion), + rightRegion(o.rightRegion) { + if (o.ext) { + ext = std::make_unique(*o.ext); + } + } + + std::string str() const { + std::stringstream os; + + os << (leftRegion != INF_IDX ? std::to_string(leftRegion) : "INF_IDX") << "/" + << (rightRegion != INF_IDX ? std::to_string(rightRegion) : "INF_IDX"); + + // use fabs(atan2P) to eliminate -0 output + + if (ext) { + os << " p: " << ext->start() << " a: " << std::fabs(atan2P(CGAL::to_double(ext->direction().dy()), CGAL::to_double(ext->direction().dx()))); + os << " EXT: "; + } + os << " p: " << iRay.start() << " a: " << std::fabs(atan2P(CGAL::to_double(iRay.direction().dy()), CGAL::to_double(iRay.direction().dx()))); + + return os.str(); + } + + std::string sstr() const { + std::stringstream os; + + // use fabs(angle()) to eliminate -0 output + + os << (leftRegion != INF_IDX ? std::to_string(leftRegion) : "INF_IDX") << "/" + << (rightRegion != INF_IDX ? std::to_string(rightRegion) : "INF_IDX"); + + return os.str(); + } + + Ray &operator=(const Ray &o) { + iRay = o.iRay; + leftRegion = o.leftRegion; + rightRegion = o.rightRegion; + + ext.reset(); + if (o.ext) { + ext = std::make_unique(*o.ext); + } + + return *this; + } + + // convenience constructor for ray unions + Ray(const Ray &upper, const Ray &lower) { + *this = lower; + this->ext = std::make_unique(upper.origin(), lower.origin()); + } + + + tOrientedSide orientedSide(const Point &x) const { + // extension is before ray, as long as its active we check against it + if (ext) { + return static_cast(ext->supporting_line().oriented_side(x)); + } else { + return static_cast(iRay.supporting_line().oriented_side(x)); + } + } + + bool isExtended() const { + return (bool) ext; + } + + Point origin() const { + if (ext) { + return ext->start(); + } else { + return iRay.start(); + } + } + + Direction direction() const { + if (ext) { + return ext->direction(); + } else { + return iRay.direction(); + } + } + + void invert() { + if (ext) { + *ext = ext->opposite(); + } + iRay = iRay.opposite(); + } + + + Point extOrigin() const { + ASSERT(ext); + return iRay.start(); + } + + void setOrigin(const Point &p_) { + iRay = CGAL::Ray_2(p_, iRay.direction()); + } + + void foldExtension() { + ASSERT(ext); + + // delete the upper ray + ext.reset(); + + ASSERT(!ext); + } + + using tIntersectionRetVal = std::optional>; + + tIntersectionRetVal intersection(const Ray &b, const Box &bounds) const { + return intersection(*this, b, bounds); + } + + public: + static tIntersectionRetVal intersection(const Ray &a, const Ray &b, const Box &bounds) { + + if (a.origin() == b.origin()) { + // we ignore rays originating at he same source + return std::nullopt; + } else { + if (!a.ext && !b.ext) { + return isRR(a, b, bounds); + } else if (a.ext && b.ext) { + return isUU(a, b, bounds); + } else if (a.ext) { + return isUR(a, b, bounds); + } else {// b.ext + return isRU(a, b, bounds); + } + } + } + +#ifdef WITH_CAIRO + + void draw(Painter &painter) const { + if (!ext) { + painter.draw(iRay.start(), 0, false); + painter.drawLine(iRay.start(), iRay.point(2)); + } else { + painter.draw(ext->start(), 0, false); + painter.drawLine(ext->start(), ext->end()); + painter.drawSquare(ext->end()); + painter.drawLine(iRay.start(), iRay.point(2)); + } + } + +#endif + + private: + template + static tIntersectionRetVal CGALintersection(const A &a, const B &b, const Box &bounds) { + + auto result = CGAL::intersection(a, b); + + if (result) { + const Point *p = boost::get(&*result); + if (p) { + // intersection is a single point, check whether its within our bounds + //TODO bound check required at all? + if (CGAL::do_intersect(*p, bounds)) { + return *p; + } else { + return std::nullopt; + } + } else { + if constexpr (std::is_same_v> && std::is_same_v>) { + const CGAL::Ray_2 *r = boost::get>(&*result); + if (r) { + // one roy lies on other ray + // direction must be the same + ASSERT(a.direction() == b.direction()); + ASSERT(r->direction() == a.direction()); + return Ray(r->source(), r->direction(), INF_IDX, INF_IDX); + } + } + + const CGAL::Segment_2 *s = boost::get>(&*result); + if (s) { + // a) one ray lies on other ray with opposite direction, segment between origin's + // both inputs are Rays + // => this should not happen + // b) extension segment lies on ray + // one input must be Segment + ASSERT((std::is_same_v> || std::is_same_v>) ); + return Ray(s->source(), s->direction(), INF_IDX, INF_IDX); + } + return std::nullopt; + } + } else { + return std::nullopt; + } + } + + static tIntersectionRetVal isRR(const Ray &a, const Ray &b, const Box &bounds) { + return CGALintersection(a.iRay, b.iRay, bounds); + } + + static tIntersectionRetVal isUR(const Ray &ua, const Ray &b, const Box &bounds) { + ASSERT(ua.ext); + + // check for intersection of pre-segment of ua and b + auto result = CGALintersection(*ua.ext, b.iRay, bounds); + + if (result) { + return result; + } else { + // check intersecton of main ray of ua and b + return isRR(ua, b, bounds); + } + } + + static tIntersectionRetVal isRU(const Ray &a, const Ray &ub, const Box &bounds) { + ASSERT(ub.ext); + + // check for intersection of pre-segment of ua and b + auto result = CGALintersection(a.iRay, *ub.ext, bounds); + + if (result) { + return result; + } else { + // check intersecton of a and main ray of ub + return isRR(a, ub, bounds); + } + } + + static tIntersectionRetVal isUU(const Ray &ua, const Ray &ub, const Box &bounds) { + ASSERT(ua.ext && ub.ext); + + // check for intersection of pre-segment of ua and ub + auto result = CGALintersection(*ua.ext, *ub.ext, bounds); + if (result) { + return result; + } + + // check for intersection of main ray of ua and pre-segment of ub + result = CGALintersection(ua.iRay, *ub.ext, bounds); + if (result) { + return result; + } + + // check for intersection of pre-segment of ua and main ray of ub + result = CGALintersection(*ua.ext, ub.iRay, bounds); + if (result) { + return result; + } + + // check for intersection of main rays of ua and ub + return isRR(ua, ub, bounds); + } + }; + + static Point mkPoint(const tIFloatVector &p) { + return Point(p[X], p[Y]); + } + + static std::vector mkPoints(const std::vector &p) { + std::vector points; + points.reserve(p.size()); + for (const auto &i : p) { + points.push_back(mkPoint(i)); + } + + return points; + } + + static Box mkBBox(const tBox &box) { + return Box(box.low[X], box.low[Y], box.high[X], box.high[Y]); + } + + static Point Midpoint(const Point &a, const Point &b) { + return CGAL::midpoint(a, b); + } + + static Direction Bisector(const Direction &l, const Direction &r) { + return Direction(0.5 * (l.base().vector() + r.base().vector())); + } + + static Direction Bisector(const Point &l, const Point &r) { + return Direction(CGAL::bisector(l, r).direction()); + } + + static Direction Bisector(const Point &l, const Point &r, const Direction &ref) { + auto b = CGAL::bisector(l, r).direction().vector(); + if (CGAL::angle(b, ref.base().vector()) == CGAL::OBTUSE) { + b *= -1; + } + return Direction(b); + } + + static Float distance2(const Point &a, const Point &b) { + return CGAL::squared_distance(a, b); + } + + // static Float distance(const Point &a, const Point &b) { + // return CGAL::sqrt(distance2(a, b)); + // } + + static bool approxEQ(const Point &a, const Point &b) { + return distance2(a, b) < MaxError::value;//TODO: use more meaningful test + } + + static bool approxEQ(const Float &a, const Float &b) { + return std::fabs(CGAL::to_double(a) - CGAL::to_double(b)) < MaxError::value;//TODO: use more meaningful test + } + + static bool approxLT(const Float &a, const Float &b) { + return a - b < 0;//TODO: use more meaningful test + } + + static bool approxGT(const Float &a, const Float &b) { + return a - b > 0;//TODO: use more meaningful test + } + + static tIFloat to_float(const Float &x) { + return CGAL::to_double(x); + } + + static tIFloat to_float_exact(const Float &x) { + if constexpr (std::is_same_v) { + x.exact(); + } + return CGAL::to_double(x); + } + + static bool compareDistance(const Point &origin, const Point &newPoint, [[maybe_unused]] const Point &oldPoint, const Float &oldDist) { + if constexpr (std::is_same_v) { + return compareDistance(origin, newPoint, oldPoint); + } else { + // we have exact constructions + return compareDistance(origin, newPoint, oldDist); + } + } + + static bool compareDistance(const Point &origin, const Point &newPoint, const Point &oldPoint) { + return CGAL::compare_distance_to_point(origin, newPoint, oldPoint) == CGAL::SMALLER; + } + + static bool compareDistance(const Point &origin, const Point &newPoint, const Float &oldDist) { + return CGAL::compare_squared_distance(origin, newPoint, oldDist) == CGAL::SMALLER; + } + + static auto computeCones(const tDim &cones) { + + CGAL::Compute_cone_boundaries_2 functor; + std::vector> rays(cones); + functor(cones, CGAL::Direction_2(1, 0), rays.begin()); + + return rays; + } + + static auto computePointCones(const Point &p, const std::vector> &rays) { + std::vector cLines; + cLines.reserve(rays.size()); + for (tDim k = 0; k < rays.size(); ++k) { + cLines.emplace_back(p, rays[k]); + } + + return cLines; + } + + static tDim getCone(const Point &origin, const Point &newPoint, const std::vector &cones) { + + tDim sec, secGuess; + sec = secGuess = std::floor(atan2P(to_float(newPoint[1] - origin[1]), to_float(newPoint[0] - origin[0])) / (2 * M_PI / cones.size())); + +#ifndef NDEBUG + bool found = false; +#endif + + for (; sec < cones.size() || sec % cones.size() < secGuess; ++sec) { + auto osL = cones[sec % cones.size()].oriented_side(newPoint); + auto osR = cones[(sec + 1) % cones.size()].oriented_side(newPoint); + + if ((osL == CGAL::ON_POSITIVE_SIDE || osL == CGAL::ON_ORIENTED_BOUNDARY) && osR == CGAL::ON_NEGATIVE_SIDE) { +#ifndef NDEBUG + found = true; +#endif + break; + } + } + ASSERT(found); + + return sec % cones.size(); + } +}; + +//TODO: there has to be a better way for template deduction +std::ostream &operator<<(std::ostream &os, const typename CGALKernel::Ray &r) { + os << r.str(); + return os; +} +std::string to_string(const typename CGALKernel::Ray &r) { + return r.str(); +} + +std::ostream &operator<<(std::ostream &os, const typename CGALKernel::Ray &r) { + os << r.str(); + return os; +} +std::string to_string(const typename CGALKernel::Ray &r) { + return r.str(); +} diff --git a/include/Kernels/InexactKernel.hpp b/include/Kernels/InexactKernel.hpp new file mode 100644 index 0000000..f8e35ec --- /dev/null +++ b/include/Kernels/InexactKernel.hpp @@ -0,0 +1,578 @@ +// +// Created by funke on 10/1/21. +// + +#pragma once + +#include +#include +#include +#include + +#include "Predicates.hpp" +#include "Types.hpp" + +#ifdef WITH_CAIRO +#include "Painter.hpp" +#endif + +class InexactKernel { + +public: + using Float = tIFloat; + using Vector = tIFloatVector; + using Point = tIFloatVector; + + static std::string name() { + return "InexactKernel"; + } + + class Direction { + public: + Direction(const tIFloat _dir) : dir(_dir), tanDir(std::tan(dir)), vec({std::cos(dir), std::sin(dir)}) {} + Direction(const Vector _vec) : Direction(atan2P(_vec[Y], _vec[X])) {} + + Float prj(const Point &p) { + return (p[X] * vec[X] + p[Y] * vec[Y]) / dot(vec, vec); + } + + Float tan() const { + return tanDir; + } + + Float cos() const { + return std::cos(dir); + } + + Float sin() const { + return std::sin(dir); + } + + tIFloat angle() const { + return dir; + } + + Direction operator-() const { + return Direction(wrapAngle(dir + M_PI)); + } + + Direction perp() const { + return {wrapAngle(dir + M_PI_2)};// TODO: check angle orientation + } + + Direction perp(const Direction &ref) const { + auto p = perp(); + return angleBetween(p.dir, ref.dir) <= M_PI_2 ? p : Direction(wrapAngle(p.dir + M_PI)); + } + + bool operator==(const Direction &o) const { + return dir == o.dir; + } + + private: + Float dir; + Float tanDir; + Vector vec; + }; + + class Line { + public: + Line(const Point &_a, const Point &_b) : a(_a), b(_b) {} + + private: + Point a; + Point b; + }; + + struct Ray { + private: + Point p; + Direction dir; + + public: + tIndex leftRegion; + tIndex rightRegion; + + private: + // extension of ray with another ray + std::unique_ptr ext;// initialized with nullptr + + public: + Ray(const Point &p_, const Direction &angle_) : p(p_), dir(angle_), leftRegion(INF_IDX), rightRegion(INF_IDX) {} + + Ray(const Point &p_, const Direction &angle_, const tIndex &lr, const tIndex &rr) + : p(p_), dir(angle_), leftRegion(lr), rightRegion(rr) {} + + Ray(const Point &p_, const Direction &angle_, const tIndex &lr, + const tIndex &rr, const Ray &ext_) + : p(p_), dir(angle_), leftRegion(lr), rightRegion(rr), + ext(std::make_unique(ext_)) { + + // the starting points of upper and lower ray should be different + ASSERT(!approxEQ(p, ext->p)); + + // the starting point of the lower ray must be on the upper rayl; + ASSERT(approxEQ(ext->p, {ext->p[X], dir.tan() * (ext->p[X] - p[X]) + p[Y]})); + ASSERT(leftRegion == ext->leftRegion); + ASSERT(rightRegion == ext->rightRegion); + } + + Ray(const Ray &o) + : p(o.p), dir(o.dir), leftRegion(o.leftRegion), + rightRegion(o.rightRegion) { + if (o.ext) { + ext = std::make_unique(*o.ext); + } + } + + std::string str() const { + std::stringstream os; + + // use fabs(angle()) to eliminate -0 output + + os << (leftRegion != INF_IDX ? std::to_string(leftRegion) : "INF_IDX") << "/" + << (rightRegion != INF_IDX ? std::to_string(rightRegion) : "INF_IDX") + << " p: " << p << " a: " << std::fabs(dir.angle()); + if (ext) { + os << " EXT: " + << " p: " << ext->p << " a: " << std::fabs(ext->dir.angle()); + ; + } + + return os.str(); + } + + std::string sstr() const { + std::stringstream os; + + // use fabs(angle()) to eliminate -0 output + + os << (leftRegion != INF_IDX ? std::to_string(leftRegion) : "INF_IDX") << "/" + << (rightRegion != INF_IDX ? std::to_string(rightRegion) : "INF_IDX"); + + return os.str(); + } + + Ray &operator=(const Ray &o) { + p = o.p; + dir = o.dir; + leftRegion = o.leftRegion; + rightRegion = o.rightRegion; + + ext.reset(); + if (o.ext) { + ext = std::make_unique(*o.ext); + } + + return *this; + } + + // convenience constructor for ray unions + Ray(const Ray &upper, const Ray &lower) : Ray(upper.p, upper.dir, upper.leftRegion, upper.rightRegion, lower) {} + + bool isExtended() const { + return (bool) ext; + } + + Point origin() const { + return p; + } + + Direction direction() const { + return dir; + } + + void invert() { + if (ext) { + ext->dir = -ext->dir; + } + + dir = -dir; + } + + Point extOrigin() const { + ASSERT(ext); + return ext->p; + } + + void setOrigin(const Point &p_) { + p = p_; + } + + void foldExtension() { + ASSERT(ext); + + // we replace the RayUnion with its lower ray, so all intersections pointers should still be valid + // all intersections processed after this point will be with lower ray + p = ext->p; + dir = ext->dir; + // and delete the lower ray + ext.reset(); + + ASSERT(!ext); + } + + tOrientedSide orientedSide(const Point &x) const { + // we only consider main ray, when the starting point of lower ray is + // swept, this ray will be replaced by it + Float res = (((p[X] + dir.cos()) - p[X]) * (x[Y] - p[Y]) - ((p[Y] + dir.sin()) - p[Y]) * (x[X] - p[X])); + return res > 0 ? tOrientedSide::LEFT : res < 0 ? tOrientedSide::RIGHT + : tOrientedSide::LINE; + } + + using tIntersectionRetVal = std::optional>; + + tIntersectionRetVal intersection(const Ray &b, const tBox &bounds) const { + return intersection(*this, b, bounds); + } + + public: + static tIntersectionRetVal intersection(const Ray &a, const Ray &b, const tBox &bounds) { + + if (a.origin() == b.origin()) { + // we ignore rays originating at he same source + return std::nullopt; + } else { + if (!a.ext && !b.ext) { + return isRR(a, b, bounds); + } else if (a.ext && b.ext) { + return isUU(a, b, bounds); + } else if (a.ext) { + return isUR(a, b, bounds); + } else {// b.ext + return isRU(a, b, bounds); + } + } + } + +#ifdef WITH_CAIRO + + void draw(Painter &painter) const { + if (!ext) { + painter.draw(p, 0, false); + painter.drawLine(p, {p[X] + dir.cos(), p[Y] + dir.sin()}); + } else { + painter.draw(p, 0, false); + painter.drawLine(p, ext->p); + painter.drawSquare(ext->p); + painter.drawLine(ext->p, {ext->p[X] + ext->dir.cos(), ext->p[Y] + ext->dir.sin()}); + } + } + +#endif + + private: + static tIntersectionRetVal isRR(const Ray &a, const Ray &b, const tBox &bounds) { + if (std::fmod(std::fmod(a.dir.angle() - b.dir.angle(), M_PI) + M_PI, M_PI) == 0) { + // parallel lines + + Direction v = b.p - a.p; + if (std::fmod(std::fmod(v.angle() - b.dir.angle(), M_PI) + M_PI, M_PI) == 0) { + // lines are on each other + if (a.dir.angle() == b.dir.angle()) { + // both lines point in the same direction -> we can determine a prior and after point + if (v.angle() == a.dir.angle()) { + // vector BA is in same direction as line -> b is before a + return Ray{a.p, b.dir}; + } else { + // vector BA is in opposite direction as line -> a is before b + return Ray{b.p, a.dir}; + } + } else { + throw std::logic_error("Parallel lines in opposing directions"); + } + } else { + // lines are parallel but not on each other -> no intersecion + return std::nullopt; + } + } + + if (std::fmod(std::fmod(a.dir.angle(), M_PI) + M_PI, M_PI) == M_PI_2) { + // vertical line this's x + return Point{a.p[X], b.dir.tan() * (a.p[X] - b.p[X]) + b.p[Y]}; + } else if (std::fmod(std::fmod(b.dir.angle(), M_PI) + M_PI, M_PI) == M_PI_2) { + // vertical line at b's x + return Point{b.p[X], a.dir.tan() * (b.p[X] - a.p[X]) + a.p[Y]}; + } + + Float m0 = a.dir.tan();// Line 0: y = m0 (x - x0) + y0 + Float m1 = b.dir.tan();// Line 1: y = m1 (x - x1) + y1 + + Float x = ((m0 * a.p[X] - m1 * b.p[X]) - (a.p[Y] - b.p[Y])) / (m0 - m1); + Float y = m0 * (x - a.p[X]) + a.p[Y]; + + if (!bounds.contains({x, y})) { + return std::nullopt; + } + + return Point{x, y}; + } + + static tIntersectionRetVal isUR(const Ray &ua, const Ray &b, const tBox &bounds) { + ASSERT(ua.ext); + + // check for intersection of main ray of ua and b + auto is = isRR(ua, b, bounds); + // check whether IS is before starting point of extension ray of ur + if (is) { + const Point *p = std::get_if(&*is); + if (p && distance2(ua.p, *p) < distance2(ua.p, ua.ext->p)) { + return is; + } + + const Ray *r = std::get_if(&*is); + if (r && distance2(ua.p, r->origin()) < distance2(ua.p, ua.ext->p)) { + return is; + } + } + + // upper ray of ur does not intersect r OR intersection is below starting point of extension ray + is = isRR(*ua.ext, b, bounds); + // check whether IS is after starting point of extension ray of ur + if (is) { + const Point *p = std::get_if(&*is); + if (p && distance2(ua.p, *p) >= distance2(ua.p, ua.ext->p)) { + return is; + } + + const Ray *r = std::get_if(&*is); + if (r && distance2(ua.p, r->origin()) >= distance2(ua.p, ua.ext->p)) { + return is; + } + } + + return std::nullopt; + } + + static tIntersectionRetVal isRU(const Ray &a, const Ray &ub, const tBox &bounds) { + ASSERT(ub.ext); + + // check for intersection of a and main ray of ub + auto is = isRR(a, ub, bounds); + // check whether IS is before starting point of extension ray of ur + if (is) { + const Point *p = std::get_if(&*is); + if (p && distance2(ub.p, *p) < distance2(ub.p, ub.ext->p)) { + return is; + } + + const Ray *r = std::get_if(&*is); + if (r && distance2(ub.p, r->origin()) < distance2(ub.p, ub.ext->p)) { + return is; + } + } + + // upper ray of ur does not intersect r OR intersection is below starting point of extension ray + is = isRR(a, *ub.ext, bounds); + // check whether IS is after starting point of extension ray of ur + if (is) { + const Point *p = std::get_if(&*is); + if (p && distance2(ub.p, *p) >= distance2(ub.p, ub.ext->p)) { + return is; + } + + const Ray *r = std::get_if(&*is); + if (r && distance2(ub.p, r->origin()) >= distance2(ub.p, ub.ext->p)) { + return is; + } + } + + return std::nullopt; + } + + static tIntersectionRetVal isUU(const Ray &ua, const Ray &ub, const tBox &bounds) { + ASSERT(ua.ext && ub.ext); + + // check for intersection of main ray of ua and ub + auto is = isRR(ua, ub, bounds); + // check whether IS is before starting point of lowerRay of ua + // check whether IS is before starting point of lowerRay of ub + if (is) { + const Point *p = std::get_if(&*is); + if (p && + distance2(ua.p, *p) < distance2(ua.p, ua.ext->p) && + distance2(ub.p, *p) < distance2(ub.p, ub.ext->p)) { + + return is; + } + + const Ray *r = std::get_if(&*is); + if (r && + distance2(ua.p, r->origin()) < distance2(ua.p, ua.ext->p) && + distance2(ub.p, r->origin()) < distance2(ub.p, ub.ext->p)) { + + return is; + } + } + + // check for intersection of main ray of ua and lower ray of ub + is = isRR(ua, *ub.ext, bounds); + // check whether IS is before starting point of lowerRay of ua + // check whether IS is after starting point of lowerRay of ub + if (is) { + const Point *p = std::get_if(&*is); + if (p && + distance2(ua.p, *p) < distance2(ua.p, ua.ext->p) && + distance2(ub.p, *p) >= distance2(ub.p, ub.ext->p)) { + return is; + } + + const Ray *r = std::get_if(&*is); + if (r && + distance2(ua.p, r->origin()) < distance2(ua.p, ua.ext->p) && + distance2(ub.p, r->origin()) >= distance2(ub.p, ub.ext->p)) { + + return is; + } + } + + // check for intersection of lower ray of ua and main ray of ub + is = isRR(*ua.ext, ub, bounds); + // check whether IS is after starting point of ua's lowerRay + // check whether IS is before starting point of ub's lowerRay + if (is) { + const Point *p = std::get_if(&*is); + if (p && + distance2(ua.p, *p) >= distance2(ua.p, ua.ext->p) && + distance2(ub.p, *p) < distance2(ub.p, ub.ext->p)) { + return is; + } + + const Ray *r = std::get_if(&*is); + if (r && + distance2(ua.p, r->origin()) >= distance2(ua.p, ua.ext->p) && + distance2(ub.p, r->origin()) < distance2(ub.p, ub.ext->p)) { + + return is; + } + } + + // check for intersection of lower rays of ua and ub + is = isRR(*ua.ext, *ub.ext, bounds); + // check whether IS is after starting point of ua's lowerRay + // check whether IS is after starting point of ub's lowerRay + if (is) { + const Point *p = std::get_if(&*is); + if (p && + distance2(ua.p, *p) >= distance2(ua.p, ua.ext->p) && + distance2(ub.p, *p) >= distance2(ub.p, ub.ext->p)) { + return is; + } + + const Ray *r = std::get_if(&*is); + if (r && + distance2(ua.p, r->origin()) >= distance2(ua.p, ua.ext->p) && + distance2(ub.p, r->origin()) >= distance2(ub.p, ub.ext->p)) { + + return is; + } + } + + return std::nullopt; + } + }; + + static Point mkPoint(const tIFloatVector &p) { + return p; + } + + static std::vector mkPoints(const std::vector &p) { + return p; + } + + static tBox mkBBox(const tBox &box) { + return box; + } + + static Point Midpoint(const Point &a, const Point &b) { + return 0.5 * (a + b); + } + + static Direction Bisector(const Direction &l, const Direction &r) { + auto dir = Direction(wrapAngle(0.5 * (l.angle() + r.angle()))); + return l.angle() < dir.angle() && dir.angle() < r.angle() ? dir : -dir; + } + + static Direction Bisector(const Point &l, const Point &r) { + return Direction(r - l).perp(); + } + + static Direction Bisector(const Point &l, const Point &r, const Direction &ref) { + return Direction(r - l).perp(ref); + } + + static Float distance2(const Point &a, const Point &b) { + Float dist2 = 0; + for (tDim d = 0; d < a.size(); ++d) { + dist2 += (b[d] - a[d]) * (b[d] - a[d]); + } + return dist2; + } + + // static Float distance(const Point &a, const Point &b) { + // return std::sqrt(distance2(a, b)); + // } + + static bool approxEQ(const Point &a, const Point &b) { + return distance2(a, b) < MaxError::value;//TODO: use more meaningful test + } + + static bool approxEQ(const Float &a, const Float &b) { + return std::fabs(a - b) < MaxError::value;//TODO: use more meaningful test + } + + static bool approxLT(const Float &a, const Float &b) { + return a - b < MaxError::value;//TODO: use more meaningful test + } + + static bool approxGT(const Float &a, const Float &b) { + return a - b > MaxError::value;//TODO: use more meaningful test + } + + static tIFloat to_float(const Float &x) { + return x; + } + + static tIFloat to_float_exact(const Float &x) { + return x; + } + + static bool compareDistance(const Point &origin, const Point &newPoint, [[maybe_unused]] const Point &oldPoint, const Float &oldDist) { + return compareDistance(origin, newPoint, oldDist); + } + + static bool compareDistance(const Point &origin, const Point &newPoint, const Point &oldPoint) { + return distance2(origin, newPoint) < distance2(origin, oldPoint); + } + + static bool compareDistance(const Point &origin, const Point &newPoint, const Float &oldDist) { + return distance2(origin, newPoint) < oldDist; + } + + static auto computeCones(const tDim &cones) { + std::vector rays; + + for (tDim k = 0; k < cones; ++k) { + rays.emplace_back(k * (2 * M_PI / cones)); + } + + return rays; + } + + static auto computePointCones([[maybe_unused]] const Point &p, const std::vector &rays) { + return rays; + } + + static tDim getCone(const Point &origin, const Point &newPoint, const std::vector &cones) { + return std::floor(atan2P(to_float(newPoint[1] - origin[1]), to_float(newPoint[0] - origin[0])) / (2 * M_PI / cones.size())); + } +}; + +std::ostream &operator<<(std::ostream &os, const InexactKernel::Ray &r) { + os << r.str(); + return os; +} + +std::string to_string(const InexactKernel::Ray &r) { + return r.str(); +} diff --git a/include/Painter.hpp b/include/Painter.hpp new file mode 100644 index 0000000..e9a92d8 --- /dev/null +++ b/include/Painter.hpp @@ -0,0 +1,280 @@ +#pragma once + +#include + +#include "Predicates.hpp" +#include "Types.hpp" + +#if (CAIROMM_MAJOR_VERSION == 1 and CAIROMM_MINOR_VERSION >= 16) +#define MY_CAIRO_FORMAT Cairo::ImageSurface::Format::ARGB32 +#define MY_CAIRO_SLANT Cairo::ToyFontFace::Slant::NORMAL +#define MY_CAIRO_WEIGHT Cairo::ToyFontFace::Weight::NORMAL +#else +#define MY_CAIRO_FORMAT Cairo::FORMAT_ARGB32 +#define MY_CAIRO_SLANT Cairo::FONT_SLANT_NORMAL +#define MY_CAIRO_WEIGHT Cairo::FONT_WEIGHT_NORMAL +#endif + +#ifdef WITH_CGAL +#include +#include +#endif + +class Painter { + +public: + Painter(const tBox &_bounds, const uint & width = 1000) { + bounds = _bounds; + + for (uint d = 0; d < 2; ++d) { + img.low[d] = 0;// image starts at 0 + img.high[d] = img.low[d] + (bounds.high[d] - bounds.low[d]); + } + + bool maxD = img.high[0] < img.high[1]; + double aspectRatio = img.high[maxD] / img.high[!maxD]; + img.high[maxD] = width; + img.high[!maxD] = width * aspectRatio; + + cs = Cairo::ImageSurface::create(MY_CAIRO_FORMAT, imgDim(0), imgDim(1)); + cr = Cairo::Context::create(cs); + + // draw background white + cr->set_source_rgb(1, 1, 1); + cr->paint(); + + cr->set_line_width(1.0); + cr->set_source_rgb(0, 0, 0); + + // set font options + cr->select_font_face("serif", MY_CAIRO_SLANT, + MY_CAIRO_WEIGHT); + cr->set_font_size(FONT_SIZE); + + // draw bounds + cr->rectangle(translatePoint(bounds.low[0], 0), + translatePoint(bounds.low[1], 1), + translateLength(bounds.high[0] - bounds.low[0], 0), + translateLength(bounds.high[1] - bounds.low[1], 1)); + cr->stroke(); + } + + Painter(const Painter &a) { + bounds = a.bounds; + img = a.img; + + cs = Cairo::ImageSurface::create(MY_CAIRO_FORMAT, imgDim(0), imgDim(1)); + cr = Cairo::Context::create(cs); + + // set font options + cr->select_font_face("serif", MY_CAIRO_SLANT, + MY_CAIRO_WEIGHT); + cr->set_font_size(FONT_SIZE); + + // copy a's surface + cr->save(); + + cr->set_source(a.cs, 0, 0); + cr->paint(); + + cr->restore(); + } + + void save(const std::string &file) const { + cs->write_to_png(file + ".png"); + } + + void setColor(const std::array &c) { + setColor(c[0], c[1], c[2]); + } + + void setColor(float r, float g, float b, + float alpha = 1) { + cr->set_source_rgba(r, g, b, alpha); + } + + void draw(const tIFloatVector &point) { + cr->arc(translatePoint(point[0], 0), + translatePoint(point[1], 1), 5, 0, 2 * M_PI); + cr->fill(); + } + + void drawSquare(const tIFloatVector &point) { + cr->rectangle(translatePoint(point[0], 0) - 5, + translatePoint(point[1], 1) - 5, 10, 10); + cr->fill(); + } + + void drawTri(const tIFloatVector &point) { + cr->move_to(translatePoint(point[0], 0) - 5, translatePoint(point[1], 1) - 5); + cr->line_to(translatePoint(point[0], 0) + 5, translatePoint(point[1], 1) - 5); + cr->line_to(translatePoint(point[0], 0), translatePoint(point[1], 1) + 5); + cr->line_to(translatePoint(point[0], 0) - 5, translatePoint(point[1], 1) - 5); + cr->fill(); + } + +#ifdef WITH_CGAL + template + void drawSquare(const CGAL::Point_2 &point) { + cr->rectangle(translatePoint(CGAL::to_double(point[0]), 0) - 5, + translatePoint(CGAL::to_double(point[1]), 1) - 5, 10, 10); + cr->fill(); + } + + template + void drawTri(const CGAL::Point_2 &point) { + cr->move_to(translatePoint(CGAL::to_double(point[0]), 0) - 5, translatePoint(CGAL::to_double(point[1]), 1) - 5); + cr->line_to(translatePoint(CGAL::to_double(point[0]), 0) + 5, translatePoint(CGAL::to_double(point[1]), 1) - 5); + cr->line_to(translatePoint(CGAL::to_double(point[0]), 0), translatePoint(CGAL::to_double(point[1]), 1) + 5); + cr->line_to(translatePoint(CGAL::to_double(point[0]), 0) - 5, translatePoint(CGAL::to_double(point[1]), 1) - 5); + cr->fill(); + } +#endif + + void draw(const tIFloatVector &point, const tIndex &idx, bool label = true) { + cr->arc(translatePoint(point[0], 0), + translatePoint(point[1], 1), 5, 0, 2 * M_PI); + cr->fill(); + + if (label) { + cr->move_to(translatePoint(point[0], 0) + 7, + translatePoint(point[1], 1) + 7); + cr->set_font_size(10); + cr->show_text(std::to_string(idx)); + cr->set_font_size(FONT_SIZE); + } + } + +#ifdef WITH_CGAL + template + void draw(const CGAL::Point_2 &point, const tIndex &idx, bool label = true) { + draw({CGAL::to_double(point[0]), CGAL::to_double(point[1])}, idx, label); + } +#endif + + void draw(const tPoints &points, bool label = true) { + for (tIndex i = 0; i < points.size(); ++i) { + draw(points[i], i, label); + } + } + + void draw(const tIndexSet &pointSet, const tPoints &points) { + for (const auto &p : pointSet) { + draw(points[p], p); + } + } + + void drawCircle(const tIFloatVector ¢er, const tIFloat &radius) { + + // draw center + cr->arc(translatePoint(center[0], 0), + translatePoint(center[1], 1), 5, 0, 2 * M_PI); + cr->fill(); + + cr->arc(translatePoint(center[0], 0), + translatePoint(center[1], 1), + translateLength(radius, 0), 0, 2 * M_PI); + cr->stroke(); + } + + void drawLine(const tIFloatVector &a, const tIFloatVector &b) { + cr->move_to(translatePoint(a[0], 0), translatePoint(a[1], 1)); + cr->line_to(translatePoint(b[0], 0), translatePoint(b[1], 1)); + cr->stroke(); + } + +#ifdef WITH_CGAL + template + void drawLine(const CGAL::Point_2 &a, const CGAL::Point_2 &b) { + cr->move_to(translatePoint(CGAL::to_double(a[0]), 0), translatePoint(CGAL::to_double(a[1]), 1)); + cr->line_to(translatePoint(CGAL::to_double(b[0]), 0), translatePoint(CGAL::to_double(b[1]), 1)); + cr->stroke(); + } +#endif + + void drawCone(const tIFloatVector &apex, const tDim &k, const tDim &K, const tIFloat length = 1) { + tIFloat lowerPhi = k * 2 * M_PI / K; + tIFloat upperPhi = (k + 1) * 2 * M_PI / K; + + drawCone(apex, lowerPhi, upperPhi, length); + } + + void drawCone(const tIFloatVector &apex, const tIFloat &lowerPhi, const tIFloat &upperPhi, const tIFloat length = 1) { + drawLine(apex, apex + tIFloatVector({length * std::cos(lowerPhi), + length * std::sin(lowerPhi)})); + drawLine(apex, apex + tIFloatVector({length * std::cos(upperPhi), + length * std::sin(upperPhi)})); + } + + void drawGrid(const tIndex _nCells) { + auto numCellsPerDim = std::ceil(std::sqrt(_nCells)); + + tIFloatVector cellSize; + for (tDim d = 0; d < bounds.high.size(); ++d) { + cellSize[d] = (bounds.high[d] - bounds.low[d]) / numCellsPerDim; + } + + for (tIndex x = 0; x < numCellsPerDim; ++x) { + for (tIndex y = 0; y < numCellsPerDim; ++y) { + cr->rectangle(translatePoint(x * cellSize[0], 0), + translatePoint(y * cellSize[1], 1), + translateLength(cellSize[0], 0), + translateLength(cellSize[1], 1)); + } + } + + cr->stroke(); + } + + void drawCones(const tIFloatVector &apex, const tDim &K, const tIFloat length = 1) { + for (tDim k = 0; k < K; ++k) { + drawCone(apex, k, K, length); + } + } + + void draw(const tIndex &i, const tYaoVertex &v, const tPoints &points) { + for (tDim k = 0; k < v.neighbor.size(); ++k) { + if (v.neighbor[k] != INF_IDX) { + drawLine(points[i], points[v.neighbor[k]]); + } + } + } + + void draw(const tYaoGraph &g, const tPoints &points) { + for (tIndex i = 0; i < g.size(); ++i) { + draw(i, g[i], points); + } + } + + void setDash() { + cr->set_dash(std::vector({2, 2}), 0); + } + + void unsetDash() { + cr->unset_dash(); + } + +private: + float imgDim(const uint dim) { + return img.high[dim] - img.low[dim]; + } + + float translatePoint(float in, uint dim) { + return ((in - bounds.low[dim]) / + ((bounds.high[dim] - bounds.low[dim]))) * + imgDim(dim); + } + + float translateLength(float in, uint dim) { + return (in / ((bounds.high[dim] - bounds.low[dim]))) * imgDim(dim); + } + +private: + Cairo::RefPtr cs; + Cairo::RefPtr cr; + + tBox bounds; + tBox img; + + static constexpr uint FONT_SIZE = 15; +}; \ No newline at end of file diff --git a/include/Predicates.hpp b/include/Predicates.hpp new file mode 100644 index 0000000..05aad7d --- /dev/null +++ b/include/Predicates.hpp @@ -0,0 +1,155 @@ +#pragma once + +#include + +#include "Types.hpp" + +tIFloat dot(const tIFloatVector &a, const tIFloatVector &b) { + tIFloat dot = 0; + for (tDim d = 0; d < a.size(); ++d) { + dot += a[d] * b[d]; + } + return dot; +} + +tIFloat wrapAngle(const tIFloat &a) { + const auto TwoPi = 2 * M_PI; + return a - TwoPi * std::floor(a / TwoPi); +} + +tIFloat angleBetween(const tIFloat &a, const tIFloat &b) { + const auto TwoPi = 2 * M_PI; + auto phi = std::fmod(std::abs(b - a), TwoPi);// This is either the distance or 360 - distance + auto distance = phi > M_PI ? TwoPi - phi : phi; + return distance; +} + +tIFloat atan2P(const tIFloat &y, const tIFloat &x) { + auto rad = std::atan2(y, x); + return rad >= 0 ? rad : 2 * M_PI + rad; +} + +tIFloat atan2P(const tIFloatVector &a, const tIFloatVector &b) { + return atan2P(b[1] - a[1], b[0] - a[0]); +} + +template::value, TA>::type, + typename = typename std::enable_if::value, TB>::type> +auto operator+(const std::array &a, const std::array &b) { + auto result = a; + for (tDim d = 0; d < a.size(); ++d) { + result[d] += b[d]; + } + return result; +} + +template::value, TA>::type, + typename = typename std::enable_if::value, TB>::type> +auto operator*(const std::array &a, const std::array &b) { + auto result = a; + for (tDim d = 0; d < a.size(); ++d) { + result[d] *= b[d]; + } + return result; +} + +template::value, TA>::type, + typename = typename std::enable_if::value, TB>::type> +auto operator-(const std::array &a, const std::array &b) { + auto result = a; + for (tDim d = 0; d < a.size(); ++d) { + result[d] -= b[d]; + } + return result; +} + +template::value, TA>::type, + typename = typename std::enable_if::value, TB>::type> +auto operator*(const TA &a, const std::array &b) { + auto result = b; + for (tDim d = 0; d < b.size(); ++d) { + result[d] *= a; + } + return result; +} + +bool operator==(const tYaoVertex &a, const tYaoVertex &b) { + for (tDim d = 0; d < a.neighbor.size(); ++d) { + if (a.neighbor[d] != b.neighbor[d]) { + return false; + } + } + + return true; +} + +std::ostream &operator<<(std::ostream &os, const tYaoVertex &v) { + + char sep = 0; + + for (tDim d = 0; d < v.neighbor.size(); ++d) { + os << sep << (v.neighbor[d] != INF_IDX ? std::to_string(v.neighbor[d]) : "I"); + sep = ' '; + } + + return os; +} + +std::ostream &operator<<(std::ostream &os, const tYaoGraph &g) { + + char sep = 0; + + for (tDim i = 0; i < g.size(); ++i) { + os << sep << g[i]; + sep = '\n'; + } + + return os; +} + +std::ostream &operator<<(std::ostream &os, const tIFloatVector &v) { + // os << "(" << v[X] << ", " << v[Y] << ")"; + os << v[X] << " " << v[Y]; + return os; +} + +std::tuple checkGraph(const tYaoGraph &is, const tYaoGraph &exp) { + + bool valid = true; + tIndexSet invalidVertices; + + if (is.size() != exp.size()) { + valid = false; + std::cerr << "size error: is " << is.size() << " exp " << exp.size() << std::endl; + } + + for (tIndex i = 0; i < is.size(); ++i) { + if (is[i] != exp[i]) { + valid = false; + invalidVertices.insert(i); + + std::cerr << "vertex error " << i << ":\n\tis: " << is[i] << "\n\texp: " << exp[i] << std::endl; + + for (tIndex k = 0; k < is[i].neighbor.size(); ++k) { + if (is[i].neighbor[k] != exp[i].neighbor[k]) { + std::cerr << "\tcone " << k << ": is: " << is[i].neighbor[k] << " (" << is[i].distance[k] + << ") exp: " << exp[i].neighbor[k] << " (" << exp[i].distance[k] << ")" << std::endl; + } + } + } + } + + if (!valid) { + std::cerr << "graph not valid"; + if (!invalidVertices.empty()) { + std::cerr << " - " << invalidVertices.size() << "/" << is.size() << " vertices invalid"; + } + std::cerr << std::endl; + } + + return {valid, invalidVertices}; +} \ No newline at end of file diff --git a/include/Types.hpp b/include/Types.hpp new file mode 100644 index 0000000..fdf8aa1 --- /dev/null +++ b/include/Types.hpp @@ -0,0 +1,108 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "Utils/ASSERT.hpp" + +using tIFloat = double;// Inexact Float +using tDim = unsigned short; +using tIndex = std::size_t; + +constexpr tDim X = 0; +constexpr tDim Y = 1; +constexpr tDim Z = 2; +constexpr tIndex INF_IDX = tIndex(-1); + +constexpr tDim D = 2; + +using tIFloatVector = std::array; +using tIndexVector = std::array; + +class tPoints : public std::vector { +public: + void setDistName(const std::string &name) { + _distName = name; + } + + std::string getDistName() const { + return _distName; + } + + void setSeed(const tIndex &seed) { + _seed = seed; + } + + tIndex getSeed() const { + return _seed; + } + +private: + std::string _distName; + tIndex _seed; +}; + +using tIndexSet = std::unordered_set; + +enum tOrientedSide { + RIGHT = -1, + LINE, + LEFT +}; + +std::string to_string(const tOrientedSide &o) { + switch (o) { + case RIGHT: + return "Right"; + case LINE: + return "Line"; + case LEFT: + return "Left"; + default: + return ""; + } +} + +template +struct MaxError; + +template<> +struct MaxError { + static constexpr float value = 1e-10; +}; + +template<> +struct MaxError { + static constexpr double value = 1e-22; +}; + +struct tBox { + tIFloatVector low; + tIFloatVector high; + + bool contains(const tIFloatVector &p) const { + ASSERT(p.size() == low.size()); + for (tDim d = 0; d < p.size(); ++d) { + if (p[d] < low[d] || high[d] < p[d]) { + return false; + } + } + + return true; + } +}; + +struct tYaoVertex { + + explicit tYaoVertex(const tDim &k) : neighbor(k, INF_IDX), distance(k, std::numeric_limits::max()) {} + + std::vector neighbor; + std::vector distance; +}; + +struct tYaoGraph : public std::vector { + tYaoGraph(const tIndex &n, const tDim &k) : std::vector(n, tYaoVertex(k)) {} +}; diff --git a/include/Utils/ASSERT.hpp b/include/Utils/ASSERT.hpp new file mode 100644 index 0000000..2ba83d6 --- /dev/null +++ b/include/Utils/ASSERT.hpp @@ -0,0 +1,52 @@ +#pragma once + +#pragma once + +#include +#include + +#ifndef NDEBUG + +#include + +#endif + +class AssertionException : public std::exception { + +public: + AssertionException() { }; + + AssertionException(const std::string &expr, const std::string &file, + const uint line) + : _what(file + ":" + std::to_string(line) + ": failed assertion `" + + expr + "'") { } + + AssertionException(const std::string &expr, const std::string &file, + const uint line, const std::string & msg) + : _what(file + ":" + std::to_string(line) + ": failed assertion `" + + expr + "' - " + msg) { } + + const char *what() const noexcept { return _what.c_str(); } + +private: + std::string _what; +}; + +#ifdef NDEBUG + +#define ASSERT(e) ((void)(0)) +#define MASSERT(e, msg) ((void)(0)) +#define RASSERT(e) ((void)(0)) + +#else // NDEBUG + +#define ASSERT(e) \ + ((void)((e) ? 0 : throw AssertionException(#e, __FILE__, __LINE__))) + +#define MASSERT(e, msg) \ + ((void)((e) ? 0 : throw AssertionException(#e, __FILE__, __LINE__, msg))) + +#define RASSERT(e) \ + ((void)((e) ? 0 : raise(SIGINT))) + +#endif diff --git a/include/Utils/ListIndexTree.hpp b/include/Utils/ListIndexTree.hpp new file mode 100644 index 0000000..bc6828b --- /dev/null +++ b/include/Utils/ListIndexTree.hpp @@ -0,0 +1,1060 @@ +// +// Created by Daniel Funke on 01.07.21. +// + +#pragma once + +#include +#include +#include + +#include "Utils/ASSERT.hpp" + +template +class SearchTree { + +private: + using index_type = std::size_t; + + template + class MemoryManager { + + public: + MemoryManager(const index_type &block_size_ = 1e5) : block_size(block_size_) { + blocks.emplace_front(); + blocks.front().reserve(block_size); + } + + N *acquire() { + if (freelist.size()) { + N *p = freelist[freelist.size() - 1]; + freelist.pop_back(); + return p; + } else { + if (blocks.begin()->size() == block_size) { + + blocks.emplace_front(); + blocks.front().reserve(block_size); + } + + blocks.front().emplace_back(); + return &blocks.front()[blocks.front().size() - 1]; + } + } + + void release(N *p) { + freelist.push_back(p); + } + + private: + std::forward_list> blocks; + index_type block_size; + std::vector freelist; + }; + +private: + // forward declaration of node types + struct InternalNode; + struct Leaf; + + using height_type = unsigned short; + + // abstract node type + struct Node { + + Node(const height_type &h) : height(h) {} + + virtual ~Node() {} + + bool isRoot() const { + return parent == nullptr; + } + + virtual bool isLeaf() const = 0; + virtual bool isNode() const = 0; + + Leaf *asLeaf() { + ASSERT(isLeaf()); + return static_cast(this); + } + + InternalNode *asNode() { + ASSERT(isNode()); + return static_cast(this); + } + + // parent is always a internal node + InternalNode *parent = nullptr; + height_type height; + }; + + struct InternalNode : public Node { + + InternalNode() : Node(0) {} + + bool isLeaf() const override { return false; } + bool isNode() const override { return true; } + + int getBalance() const { + return (left ? left->height : 0) - (right ? right->height : 0); + } + + Node *left = nullptr; + Node *right = nullptr; + + Leaf *maxRep = nullptr; + }; + + struct Leaf : public Node { + + Leaf() : Node(1) {} + + bool isLeaf() const override { return true; } + bool isNode() const override { return false; } + + std::unique_ptr obj; + + Leaf *prev = nullptr; + Leaf *next = nullptr; + }; + +public: + class Iterator { + + friend class SearchTree; + + protected: + Leaf *leaf_; + + public: + using difference_type = short; + using value_type = T; + using pointer = T *; + using const_pointer = const T *; + using reference = T &; + using const_reference = const T &; + using iterator_category = std::bidirectional_iterator_tag; + + using it_pointer = Leaf *; + using const_it_pointer = const Leaf *; + + Iterator() : leaf_(nullptr) {} + Iterator(Leaf *leaf) : leaf_(leaf) {} + + public: + Iterator &operator++() { + leaf_ = leaf_->next; + return *this; + } + + Iterator &operator--() { + leaf_ = leaf_->prev; + return *this; + } + + bool operator==(Iterator b) const { return leaf_ == b.leaf_; } + bool operator!=(Iterator b) const { return leaf_ != b.leaf_; } + + reference operator*() { return *(leaf_->obj.get()); } + pointer operator->() { return leaf_->obj.get(); } + + const_reference operator*() const { return *(leaf_->obj.get()); } + const_pointer operator->() const { return leaf_->obj.get(); } + + it_pointer addr() { + return leaf_; + } + + const_it_pointer addr() const { + return leaf_; + } + }; + +public: + SearchTree() { + m_root = newInternalNode(); + + m_endLeaf = newLeaf(); + m_endLeaf->prev = m_endLeaf; + m_endLeaf->next = m_endLeaf; + + m_beginLeaf = m_endLeaf; + } + + Iterator begin() { + return Iterator(_begin()); + } + + Iterator end() { + return Iterator(_end()); + } + +private: + Leaf *_begin() { + return m_beginLeaf; + } + + Leaf *_end() { + return m_endLeaf; + } + +public: + Iterator insert(Iterator pos, const T &obj) { + ASSERT(checkInvariants()); + ASSERT(checkList()); + auto it = Iterator(insert(pos.leaf_, obj)); + ASSERT(checkInvariants()); + ASSERT(checkList()); + return it; + } + + Iterator insert_pair(Iterator pos, const T &left, const T &right) { + ASSERT(checkInvariants()); + ASSERT(checkList()); + auto it = Iterator(insert_pair(pos.leaf_, left, right)); + ASSERT(checkInvariants()); + ASSERT(checkList()); + return it; + } + +private: + Leaf *newLeaf() { + return mm_leafs.acquire(); + } + + InternalNode *newInternalNode() { + return mm_inodes.acquire(); + } + + void deleteNode(Node *p) { + if (p->isLeaf()) { + deleteLeaf(p->asLeaf()); + } else { + deleteInternalNode(p->asNode()); + } + } + + void deleteLeaf(Leaf *p) { + p->next = nullptr; + p->prev = nullptr; + p->parent = nullptr; + p->height = 1; + p->obj.reset(); + + mm_leafs.release(p); + } + + void deleteInternalNode(InternalNode *p) { + p->left = nullptr; + p->right = nullptr; + p->parent = nullptr; + p->height = 0; + p->maxRep = nullptr; + + mm_inodes.release(p); + } + +private: + Leaf *insert(Leaf *pos, const T &obj) { + + // construct new leaf + Leaf *leaf = newLeaf(); + leaf->obj = std::make_unique(obj); + + // special case empty list + if (pos == _begin() && pos == _end()) { + + ASSERT(!m_root->left && !m_root->right); + ASSERT(_begin() == _end()); + + leaf->prev = _end(); + leaf->next = _end(); + leaf->parent = m_root; + + m_beginLeaf = leaf; + m_endLeaf->prev = leaf; + + m_root->left = leaf; + updateAndRebalance(m_root); + + return m_root->left->asLeaf(); + } + + // take care of the leaf double linked list + if (pos == _end()) { + leaf->prev = m_endLeaf->prev; + leaf->next = _end(); + + leaf->prev->next = leaf; + + m_endLeaf->prev = leaf; + } else if (pos == _begin()) { + leaf->prev = _end(); + leaf->next = m_beginLeaf; + + leaf->next->prev = leaf; + + m_beginLeaf = leaf; + } else { + leaf->prev = pos->prev; + ASSERT(leaf->prev != _end()); + leaf->next = pos; + ASSERT(leaf->next == pos); + + leaf->prev->next = leaf; + leaf->next->prev = leaf; + } + + // insert leaf into tree + if (leaf->prev != _end()) { + // we have a left neighbor, join to it from right + joinFromRight(leaf->prev->parent, leaf); + } else { + ASSERT(leaf->next != _end()); + // we have no left neighbor but a right one, join to it from left + joinFromLeft(leaf->next->parent, leaf); + } + + return leaf; + } + + Leaf *insert_pair(Leaf *pos, const T &left, const T &right) { + + // construct new leafs + Leaf *leftLeaf = newLeaf(); + leftLeaf->obj = std::make_unique(left); + + Leaf *rightLeaf = newLeaf(); + rightLeaf->obj = std::make_unique(right); + + leftLeaf->next = rightLeaf; + rightLeaf->prev = leftLeaf; + + // special case empty list + if (pos == _begin() && pos == _end()) { + + ASSERT(!m_root->left && !m_root->right); + ASSERT(_begin() == _end()); + + leftLeaf->prev = _end(); + rightLeaf->next = _end(); + leftLeaf->parent = m_root; + rightLeaf->parent = m_root; + + m_beginLeaf = leftLeaf; + m_endLeaf->prev = rightLeaf; + + m_root->left = leftLeaf; + m_root->right = rightLeaf; + updateAndRebalance(m_root); + + return m_root->left->asLeaf(); + } + + // take care of the leaf double linked list + if (pos == _end()) { + leftLeaf->prev = m_endLeaf->prev; + rightLeaf->next = _end(); + + leftLeaf->prev->next = leftLeaf; + + m_endLeaf->prev = rightLeaf; + } else if (pos == _begin()) { + leftLeaf->prev = _end(); + rightLeaf->next = m_beginLeaf; + + rightLeaf->next->prev = rightLeaf; + + m_beginLeaf = leftLeaf; + } else { + leftLeaf->prev = pos->prev; + ASSERT(leftLeaf->prev != _end()); + rightLeaf->next = pos; + ASSERT(rightLeaf->next == pos); + + leftLeaf->prev->next = leftLeaf; + rightLeaf->next->prev = rightLeaf; + } + + // insert leaf into tree + if (leftLeaf->prev != _end()) { + // we have a left neighbor, join to it from right + joinFromRight(leftLeaf->prev->parent, leftLeaf, rightLeaf); + } else { + ASSERT(leftLeaf->next != _end()); + // we have no left neighbor but a right one, join to it from left + joinFromLeft(rightLeaf->next->parent, leftLeaf, rightLeaf); + } + + return leftLeaf; + } + +private: + void joinFromRight(InternalNode *parent, Leaf *leaf) { + + ASSERT(parent->left); + ASSERT(leaf->prev != _end()); + + bool prevIsLeftChild = (parent->left == leaf->prev); + ASSERT(prevIsLeftChild || (parent->right && parent->right == leaf->prev)); + + if (!parent->right) { + // parent has empty right child, prev must be left child + // insert new leaf as right child + ASSERT(prevIsLeftChild); + + leaf->parent = parent; + + parent->right = leaf; + + updateAndRebalance(parent); + + return; + } + + InternalNode *newNode = newInternalNode(); + newNode->parent = parent; + + newNode->left = prevIsLeftChild ? parent->left : parent->right; + newNode->left->parent = newNode; + + newNode->right = leaf; + newNode->right->parent = newNode; + + (prevIsLeftChild ? parent->left : parent->right) = newNode; + updateAndRebalance((prevIsLeftChild ? parent->left : parent->right)->asNode()); + } + + void joinFromRight(InternalNode *parent, Leaf *leftLeaf, Leaf *rightLeaf) { + + ASSERT(parent->left); + ASSERT(leftLeaf->prev != _end()); + + bool prevIsLeftChild = (parent->left == leftLeaf->prev); + ASSERT(prevIsLeftChild || (parent->right && parent->right == leftLeaf->prev)); + + InternalNode *pNode = newInternalNode(); + pNode->left = leftLeaf; + leftLeaf->parent = pNode; + pNode->right = rightLeaf; + rightLeaf->parent = pNode; + updateNodeInfo(pNode); + + if (!parent->right) { + // parent has empty right child, prev must be left child + // insert new leaf as right child + ASSERT(prevIsLeftChild); + + pNode->parent = parent; + + parent->right = pNode; + + updateAndRebalance(parent); + + return; + } + + InternalNode *newNode = newInternalNode(); + newNode->parent = parent; + + newNode->left = prevIsLeftChild ? parent->left : parent->right; + newNode->left->parent = newNode; + + newNode->right = pNode; + newNode->right->parent = newNode; + + (prevIsLeftChild ? parent->left : parent->right) = newNode; + updateAndRebalance((prevIsLeftChild ? parent->left : parent->right)->asNode()); + } + + void joinFromLeft(InternalNode *parent, Leaf *leaf) { + + ASSERT(parent->left); + // this is only called when object is inserted to the beginning of list + // new leaf must become leftmost leaf + ASSERT(leaf->prev == _end()); + + if (!parent->right) { + // parent has empty right child, move left to right, then insert + + leaf->parent = parent; + + parent->right = parent->left; + parent->left = leaf; + + updateAndRebalance(parent); + + return; + } + + InternalNode *newNode = newInternalNode(); + newNode->parent = parent; + newNode->left = leaf; + newNode->left->parent = newNode; + + newNode->right = parent->left; + newNode->right->parent = newNode; + + parent->left = newNode; + updateAndRebalance(parent->left->asNode()); + } + + void joinFromLeft(InternalNode *parent, Leaf *leftLeaf, Leaf *rightLeaf) { + + ASSERT(parent->left); + // this is only called when object is inserted to the beginning of list + // new leaf must become leftmost leaf + ASSERT(leftLeaf->prev == _end()); + + InternalNode *pNode = newInternalNode(); + pNode->left = leftLeaf; + leftLeaf->parent = pNode; + pNode->right = rightLeaf; + rightLeaf->parent = pNode; + updateNodeInfo(pNode); + + if (!parent->right) { + // parent has empty right child, move left to right, then insert + + pNode->parent = parent; + + parent->right = parent->left; + parent->left = pNode; + + updateAndRebalance(parent); + + return; + } + + InternalNode *newNode = newInternalNode(); + newNode->parent = parent; + newNode->left = pNode; + newNode->left->parent = newNode; + + newNode->right = parent->left; + newNode->right->parent = newNode; + + parent->left = newNode; + updateAndRebalance(parent->left->asNode()); + } + +public: + Iterator erase(Iterator pos) { + ASSERT(checkInvariants()); + ASSERT(checkList()); + auto it = Iterator(erase(pos.leaf_)); + ASSERT(checkInvariants()); + ASSERT(checkList()); + return it; + } + + Iterator replace(Iterator pos, const T &obj) { + ASSERT(checkInvariants()); + ASSERT(checkList()); + auto it = Iterator(replace(pos.leaf_, obj)); + ASSERT(checkInvariants()); + ASSERT(checkList()); + return it; + } + +private: + Leaf *erase(Leaf *pos) { + ASSERT(pos != nullptr); + ASSERT(pos != _end()); + ASSERT(pos->parent != nullptr); + + // special case singleton list + if (pos == _begin() && pos == _end()->prev) { + + ASSERT(m_root->left && !m_root->right); + ASSERT(m_root->left->isLeaf()); + + m_beginLeaf = _end(); + m_endLeaf->prev = _end(); + + deleteLeaf(m_root->left->asLeaf()); + m_root->left = nullptr; + updateAndRebalance(m_root); + + return _end(); + } + + // take care of the leaf double linked list + Leaf *retValue = nullptr;// save return value -> node after deleted one + + if (pos == _end()->prev) { + ASSERT(pos->prev != _end());// singleton case already handled above + pos->prev->next = _end(); + m_endLeaf->prev = pos->prev; + retValue = _end(); + } else if (pos == _begin()) { + ASSERT(pos->next != _end());// singleton case already handled above + pos->next->prev = _end(); + m_beginLeaf = pos->next; + retValue = _begin(); + } else { + pos->prev->next = pos->next; + pos->next->prev = pos->prev; + retValue = pos->next; + } + + // delete actual leaf + auto parent = pos->parent; + ASSERT(parent != nullptr); + ASSERT(parent->left == pos || parent->right == pos); + + if (parent->right == pos) { + ASSERT(parent->left); + // parent has two children and we are the right child, simply remove + deleteNode(parent->right); + parent->right = nullptr; + } else { + ASSERT(parent->left == pos); + // we are the left child, move right child into left, maybe nullptr + deleteNode(parent->left); + parent->left = parent->right;//maybe nullptr + parent->right = nullptr; + } + + //parent cannot have a right child anymore + ASSERT(!parent->right); + + // recursivley delete InternalNodes if they become empty or singletons + eraseRec(parent); + + return retValue; + } + + Leaf *replace(Leaf *pos, const T &obj) { + ASSERT(pos != nullptr); + ASSERT(pos != _end()); + ASSERT(pos->parent != nullptr); + ASSERT(pos->isLeaf()); + ASSERT(pos->obj); + ASSERT(pos->prev); + ASSERT(pos->next); + + *pos->obj = obj; + + return pos; + } + +private: + void eraseRec(InternalNode *node) { + ASSERT(node != nullptr); + + if (node->isRoot()) { + ASSERT(node->parent == nullptr); + updateAndRebalance(node); + return; + } + + auto parent = node->parent; + ASSERT(parent != nullptr); + + if (node->left && node->right) { + // node has both children, abort deletion recursion and just update and rebalance + updateAndRebalance(node); + } else { + // node has definitely no right child + ASSERT(!node->right); + + ASSERT(parent->left == node || parent->right == node); + + if (parent->right == node) { + parent->right = node->left;// maybe nullptr + if (parent->right) { + parent->right->parent = parent; + } + } else { + ASSERT(parent->left == node); + parent->left = node->left;// maybe nullptr + if (parent->left) { + parent->left->parent = parent; + } + } + + deleteNode(node); + eraseRec(parent); + } + } + +private: + // A utility function to right + // rotate subtree rooted with y + // See the diagram given above. + InternalNode *rightRotate(InternalNode *y) { + + // Node *yFromParent; + bool leftChild = false; + + if (y->parent == nullptr) { + ASSERT(y->isRoot() && y == m_root); + // yFromParent = m_root; + } else { + leftChild = (y == y->parent->left); + ASSERT(leftChild || y == y->parent->right); + // yFromParent = (leftChild ? y->parent->left : y->parent->right); + } + + // ASSERT(yFromParent && y == yFromParent); + + ASSERT(y->left); + Node *x = y->left; + ASSERT(x->isNode() && x->asNode()->left); + Node *T2 = x->asNode()->right; + + // perform rotation + x->parent = y->parent; + + if (T2) { + y->left = T2; + y->left->parent = y; + } else if (y->right) { + // if T2 is empty move right child over + y->left = y->right; + y->right = nullptr; + } else { + // y becomes empty, delete it + deleteInternalNode(y); + y = nullptr; + } + + // perform rotation + x->asNode()->right = y; + + if (y) { + // update info + updateNodeInfo(y); + x->asNode()->right->parent = x->asNode(); + } + + // update info + updateNodeInfo(x->asNode()); + + // update fromParentPointer + if (x->parent == nullptr) { + ASSERT(x->isRoot()); + m_root = static_cast(x); + + // Return new root + return m_root; + + } else { + Node **pPtr = (leftChild ? &(x->parent->left) : &(x->parent->right)); + (*pPtr) = x; + + // Return new root + return (*pPtr)->asNode(); + } + } + + // A utility function to left + // rotate subtree rooted with x + // See the diagram given above. + InternalNode *leftRotate(InternalNode *x) { + + // Node *xFromParent; + bool leftChild = false; + + if (x->parent == nullptr) { + ASSERT(x->isRoot() && x == m_root); + // xFromParent = m_root; + } else { + leftChild = (x == x->parent->left); + ASSERT(leftChild || x == x->parent->right); + // xFromParent = (leftChild ? x->parent->left : x->parent->right); + } + + // ASSERT(xFromParent && x == xFromParent); + + ASSERT(x->left && x->right); + Node *y = x->right; + ASSERT(y->isNode() && y->asNode()->left); + Node *T2 = y->asNode()->left; + + // perform rotation + x->right = T2; + x->right->parent = x; + + // update info + updateNodeInfo(x); + + // perform rotation + y->parent = x->parent; + y->asNode()->left = x; + y->asNode()->left->parent = y->asNode(); + + // update info + updateNodeInfo(y->asNode()); + + // update fromParentPointer + if (y->parent == nullptr) { + ASSERT(y->isRoot()); + m_root = static_cast(y); + + // Return new root + return m_root; + } else { + Node **pPtr = (leftChild ? &(y->parent->left) : &(y->parent->right)); + (*pPtr) = y; + + // Return new root + return (*pPtr)->asNode(); + } + } + + bool updateNodeInfo(InternalNode *node) { + + auto oldRep = node->maxRep; + auto oldHeight = node->height; + + if (node->right) { + ASSERT(node->left); + node->maxRep = (node->right->isNode() ? node->right->asNode()->maxRep : node->right->asLeaf()); + node->height = 1 + std::max(height_type(node->left->height), + height_type(node->right->height)); + } else if (node->left) { + ASSERT(!node->right); + node->maxRep = (node->left->isNode() ? node->left->asNode()->maxRep : node->left->asLeaf()); + node->height = 1 + (node->left->height); + } else { + ASSERT(!node->left && !node->right); + node->maxRep = nullptr; + node->height = 0; + } + + return not(oldHeight == node->height && oldRep == node->maxRep); + } + + void updateAndRebalance(InternalNode *node /*, bool rebalanceRequired = true, bool updateRequired = true*/) { + + //if (updateRequired) { + /*updateRequired = */ updateNodeInfo(node); + // } + + //if (rebalanceRequired) { + int balance = node->getBalance(); + + // If this node becomes unbalanced, + // then there are 4 cases + // if this node is unbalanced and we rebalance it, then we do _not_ need to traverese further upward + + if (balance > 1) { + int leftBalance = (node->left && node->left->isNode() ? node->left->asNode()->getBalance() : 0); + + // Left Left Case + if (balance > 1 && leftBalance >= 0) { + node = rightRotate(node); + } + + // Left Right Case + if (balance > 1 && leftBalance < 0) { + ASSERT(node->left && node->left->isNode()); + leftRotate(node->left->asNode()); + node = rightRotate(node); + } + } else if (balance < -1) { + int rightBalance = (node->right && node->right->isNode() ? node->right->asNode()->getBalance() : 0); + + // Right Right Case + if (balance < -1 && rightBalance <= 0) { + node = leftRotate(node); + } + + // Right Left Case + if (balance < -1 && rightBalance > 0) { + ASSERT(node->right && node->right->isNode()); + rightRotate(node->right->asNode()); + node = leftRotate(node); + } + } + + /* + // if we rebalanced the node, we do not need to rebalance further up the tree + rebalanceRequired = balance <= std::abs(1); + + // if we rebalanced we need to update + updateRequired = updateRequired || balance >= std::abs(1); + */ + //} + + if (node->parent != nullptr /* && (updateRequired || rebalanceRequired)*/) { + // recursively traverse up to check whether node above became unbalanced + updateAndRebalance(node->parent /*, rebalanceRequired, updateRequired*/); + } + } + +public: + template + Iterator find(const O &obj, const Compare &cmp) { + return Iterator(find(m_root, obj, cmp)); + } + +private: + template + Leaf *find(InternalNode *node, const O &obj, const Compare &cmp) { + + if (node->maxRep != nullptr && cmp(obj, *(node->maxRep->obj))) { + ASSERT(node->left); + + if (node->left->isNode()) { + if (cmp(obj, *(node->left->asNode()->maxRep->obj))) { + return find(node->left->asNode(), obj, cmp); + } else { + ASSERT(node->right); + + if (node->right->isNode()) { + return find(node->right->asNode(), obj, cmp); + } else { + return node->right->asLeaf(); + } + } + } else { + if (cmp(obj, *(node->left->asLeaf()->obj))) { + return node->left->asLeaf(); + } else { + ASSERT(node->right); + + if (node->right->isNode()) { + return find(node->right->asNode(), obj, cmp); + } else { + return node->right->asLeaf(); + } + } + } + } else { + return _end(); + } + } + +public: + bool checkInvariants() const { + return checkInvariants(m_root); + } + + bool checkList() { + bool valid = true; + auto it = begin(); + valid &= (std::prev(it) == end()); + + it = std::next(it); + for (; it != end() && valid; ++it) { + valid &= (it == std::next(std::prev(it))); + valid &= (it == std::prev(std::next(it))); + } + + return valid; + } + + std::string to_string() { + return to_string(m_root, 0); + } + +private: + bool checkInvariants(InternalNode *node) const { + if (std::abs(node->getBalance()) > 1) { + return false; + } + + if (!node->left && node->right) { + return false; + } + + if (node->maxRep) { + Node *n = node; + while (!n->isLeaf()) { + ASSERT(n && n->isNode()); + if (n->asNode()->right) { + n = n->asNode()->right; + } else { + n = n->asNode()->left; + } + } + ASSERT(n && n->isLeaf()); + + if (node->maxRep->obj != n->asLeaf()->obj) { + return false; + } + } else { + // only root my not have a rep if it is empty + if (!node->isRoot() && !node->left && !node->right) { + return false; + } + } + + bool valid = true; + if (node->left) { + if (node->left->isNode()) { + valid = checkInvariants(node->left->asNode()); + } else { + ASSERT(node->left->isLeaf()); + valid = bool(node->left->asLeaf()->obj); + } + } + + if (!valid) { return false; } + + if (node->right) { + if (node->right->isNode()) { + valid = checkInvariants(node->right->asNode()); + } else { + ASSERT(node->right->isLeaf()); + valid = bool(node->right->asLeaf()->obj); + } + } + + return valid; + } + + std::string to_string(Node *root, int space) { + // Base case + if (root == NULL) + return ""; + + // Increase distance between levels + const int SINC = 20; + space += SINC; + std::stringstream ss; + + // Process right child first + if (root->isNode()) { + ss << to_string(root->asNode()->right, space); + } + + // Print current node after space + // count + ss << '\n'; + for (int i = SINC; i < space; i++) + ss << ' '; + + if (root->isNode()) { + ss << root << "(" << root->height << ", " << root->asNode()->getBalance() << ")"; + } else { + ASSERT(root->isLeaf()); + if constexpr (std::is_class_v) { + ss << root->asLeaf()->obj->sstr(); + } else { + ss << *(root->asLeaf()->obj); + } + } + + ss << "\n"; + + // Process left child + if (root->isNode()) { + ss << to_string(root->asNode()->left, space); + } + + return ss.str(); + } + +private: + InternalNode *m_root; + + Leaf *m_beginLeaf = nullptr; + Leaf *m_endLeaf; + + MemoryManager mm_leafs; + MemoryManager mm_inodes; +}; diff --git a/include/Utils/Logging.hpp b/include/Utils/Logging.hpp new file mode 100644 index 0000000..458d5f3 --- /dev/null +++ b/include/Utils/Logging.hpp @@ -0,0 +1,9 @@ +#pragma once + +#ifdef LOG_DEBUG +#define LOG(msg) \ + std::cout << msg +#else +#define LOG(msg) \ + ((void)(0)) +#endif \ No newline at end of file diff --git a/include/Utils/PriorityQueue.hpp b/include/Utils/PriorityQueue.hpp new file mode 100644 index 0000000..97f0f63 --- /dev/null +++ b/include/Utils/PriorityQueue.hpp @@ -0,0 +1,310 @@ +#pragma once + +#include +#include + +template> +class PriQueueD { +public: + PriQueueD(size_t capacity, size_t log_degree = 3) + : _capacity(capacity), _ldeg(log_degree) { + _heap_element_vec.reserve(capacity); + _heap_hindex_vec.reserve(capacity); + _handle_vec.reserve(capacity); + _free_handle_vec.reserve(capacity); + } + ~PriQueueD() {} + +private: + // identifier object created on push, and deleted on pop + +public: + using handle = size_t; + + const T &top() const { + ASSERT(!empty()); + return _heap_element_vec[0]; + } + + bool empty() const { return _heap_element_vec.empty(); } + size_t size() const { return _heap_element_vec.size(); } + + handle push(T value) { + auto h = get_handle(); + + auto pos = _heap_element_vec.size(); + // * First Iteration of sift_up **************************************** + size_t par = parent(pos); + if (_heap_element_vec.empty() || compare(_heap_element_vec[par], value)) { + _heap_element_vec.emplace_back(std::move(value)); + _heap_hindex_vec.emplace_back(h); + _handle_vec[h] = pos; + return h; + } + _heap_element_vec.emplace_back(std::move(_heap_element_vec[par])); + _heap_hindex_vec.emplace_back(_heap_hindex_vec[par]); + _handle_vec[_heap_hindex_vec[par]] = pos; + + pos = sift_up(par, value); + _heap_element_vec[pos] = std::move(value); + _heap_hindex_vec[pos] = h; + _handle_vec[_heap_hindex_vec[pos]] = pos; + + return _heap_hindex_vec[pos]; + } + + void pop() { + ASSERT(!empty()); + + auto temp_element = std::move(_heap_element_vec[_heap_element_vec.size() - 1]); + auto temp_handle = _heap_hindex_vec[_heap_hindex_vec.size() - 1]; + + _heap_element_vec.pop_back(); + _heap_hindex_vec.pop_back(); + + _free_handle_vec.push_back(_heap_hindex_vec[0]); + if (!_heap_element_vec.size()) return; + + auto pos = sift_down(0, temp_element); + _heap_element_vec[pos] = std::move(temp_element); + _heap_hindex_vec[pos] = temp_handle; + _handle_vec[temp_handle] = pos; + } + + void remove(handle h) { + ASSERT(!empty()); + + auto pos = _handle_vec[h]; + + auto temp_element = std::move(_heap_element_vec[_heap_element_vec.size() - 1]); + auto temp_handle = _heap_hindex_vec[_heap_hindex_vec.size() - 1]; + + _heap_element_vec.pop_back(); + _heap_hindex_vec.pop_back(); + + _free_handle_vec.push_back(h); + if (pos >= _heap_element_vec.size()) return; + + if (compare(temp_element, _heap_element_vec[pos])) { + pos = sift_up(pos, temp_element); + } else { + pos = sift_down(pos, temp_element); + } + _heap_element_vec[pos] = std::move(temp_element); + _heap_hindex_vec[pos] = temp_handle; + _handle_vec[temp_handle] = pos; + } + + void change_key(handle h, T newvalue) { + + auto pos = _handle_vec[h]; + if (compare(newvalue, _heap_element_vec[pos])) { + pos = sift_up(pos, newvalue); + } else { + pos = sift_down(pos, newvalue); + } + _heap_element_vec[pos] = std::move(newvalue); + _heap_hindex_vec[pos] = h; + _handle_vec[h] = pos; + } + + const T &get_key(handle h) const { + return _heap_element_vec[_handle_vec[h]]; + } + + bool verify() const { + return verify_node(0); + } + +private: + /* member definitions *****************************************************/ + //std::priority_queue, Comp> pq; + size_t _capacity; + size_t _ldeg; + Comp _comp; + std::vector _heap_hindex_vec; + std::vector _heap_element_vec; + std::vector _handle_vec; + std::vector _free_handle_vec; + + + bool compare(const T &v1, const T &v2) const { + return _comp(v2, v1); + } + + void move(size_t from, size_t too) { + _heap_hindex_vec[too] = _heap_hindex_vec[from]; + _heap_element_vec[too] = std::move(_heap_element_vec[from]); + _handle_vec[_heap_hindex_vec[too]] = too; + } + + size_t get_handle() { + if (_free_handle_vec.size()) { + size_t b = _free_handle_vec[_free_handle_vec.size() - 1]; + _free_handle_vec.pop_back(); + return b; + } else { + if (_handle_vec.size() == _capacity) { + _capacity *= 2; + + _heap_element_vec.reserve(_capacity); + _heap_hindex_vec.reserve(_capacity); + _handle_vec.reserve(_capacity); + } + + _handle_vec.push_back(0ull); + return _handle_vec.size() - 1; + } + } + + size_t parent(size_t index) const { + return (index - 1) >> _ldeg; + } + + size_t child(size_t index) const { + return (index << _ldeg) + 1; + } + + size_t sift_up(size_t pos, const T &val) { + while (pos != 0) { + size_t par = parent(pos); + if (compare(_heap_element_vec[par], val)) return pos; + move(par, pos); + pos = par; + } + return 0; + } + + size_t sift_down(size_t pos, const T &val) { + // while not in leaf + while (child(pos) < _heap_element_vec.size()) { + size_t min_idx = child(pos); + + // find min child + for (size_t i = child(pos) + 1; + i < child(pos + 1) && i < _heap_element_vec.size(); + ++i) { + min_idx = compare(_heap_element_vec[min_idx], _heap_element_vec[i]) ? min_idx : i; + } + + if (compare(val, _heap_element_vec[min_idx])) return pos; + //_vec[pos] = std::move(_vec[min_idx]); + move(min_idx, pos); + pos = min_idx; + } + return pos; + } + + bool verify_node(size_t pos) const { + // check all children + for (size_t i = child(pos); + i < child(pos + 1) && i < _heap_element_vec.size(); + ++i) { + if (!compare(_heap_element_vec[pos], _heap_element_vec[i]) || !verify_node(i)) + return false; + } + + return true; + } +}; + + +template> +class PriQueueAdapter { + +public: + using pqType = PriQueueD; + using arrType = std::vector; + using handle = typename pqType::handle; + + PriQueueAdapter(arrType &&static_elements, size_t capacityPQ, size_t log_degree = 3) + : arr(std::move(static_elements)), pq(capacityPQ, log_degree) { + + // initially sort the static_elements according to Comp + std::sort(arr.begin(), arr.end(), comp); + arrIdx = arr.size(); + } + +public: + const T &top() const { + ASSERT(!empty()); + + if (arrIdx > 0 && !pq.empty()) { + // both are not empty, we need to compare + + if (compare(arr[arrIdx - 1], pq.top())) { + return arr[arrIdx - 1]; + } else { + return pq.top(); + } + } else { + // one is empty, can't be both per ASSERT + + if (!pq.empty()) { + ASSERT(arrIdx == 0); + return pq.top(); + } else { + ASSERT(arrIdx > 0); + return arr[arrIdx - 1]; + } + } + } + + bool empty() const { return arrIdx == 0 && pq.empty(); } + size_t size() const { return arrIdx + pq.size(); } + + handle push(T value) { + return pq.push(value); + } + + void pop() { + ASSERT(!empty()); + + if (arrIdx > 0 && !pq.empty()) { + // both are not empty, we need to compare + + if (compare(arr[arrIdx - 1], pq.top())) { + arrIdx--; + } else { + pq.pop(); + } + } else { + // one is empty, can't be both per ASSERT + + if (!pq.empty()) { + ASSERT(arrIdx == 0); + pq.pop(); + } else { + ASSERT(arrIdx > 0); + arrIdx--; + } + } + } + + void remove(handle h) { + pq.remove(h); + } + + void change_key(handle h, T newvalue) { + pq.change_key(h, newvalue); + } + + const T &get_key(handle h) const { + return pq.get_key(h); + } + + bool verify() const { + return pq.verify(); + } + +private: + bool compare(const T &v1, const T &v2) const { + return comp(v2, v1); + } + +private: + arrType arr; + pqType pq; + typename arrType::size_type arrIdx; + Comp comp; +}; \ No newline at end of file diff --git a/include/Utils/Timer.hpp b/include/Utils/Timer.hpp new file mode 100644 index 0000000..552e4e8 --- /dev/null +++ b/include/Utils/Timer.hpp @@ -0,0 +1,23 @@ +// +// Created by Daniel Funke on 17.03.21. +// + +#pragma once + +#include + +template +class Timer { + +public: + template + static auto time(Args... args) { + T obj; + + auto t1 = std::chrono::high_resolution_clock::now(); + auto res = obj(args...); + auto t2 = std::chrono::high_resolution_clock::now(); + + return std::make_tuple(std::chrono::duration_cast(t2 - t1).count(), res); + } +}; \ No newline at end of file diff --git a/include/contrib/popl.hpp b/include/contrib/popl.hpp new file mode 100644 index 0000000..8fd668b --- /dev/null +++ b/include/contrib/popl.hpp @@ -0,0 +1,1331 @@ +/*** + ____ __ ____ __ + ( _ \ / \( _ \( ) + ) __/( O )) __// (_/\ + (__) \__/(__) \____/ + version 1.3.0 + https://github.com/badaix/popl + + This file is part of popl (program options parser lib) + Copyright (C) 2015-2021 Johannes Pohl + + This software may be modified and distributed under the terms + of the MIT license. See the LICENSE file for details. +***/ + +/// checked with clang-tidy: +/// run-clang-tidy-3.8.py -header-filter='.*' +/// -checks='*,-misc-definitions-in-headers,-google-readability-braces-around-statements,-readability-braces-around-statements,-cppcoreguidelines-pro-bounds-pointer-arithmetic,-google-build-using-namespace,-google-build-using-namespace' + +#ifndef POPL_HPP +#define POPL_HPP + +#ifndef NOMINMAX +#define NOMINMAX +#endif // NOMINMAX + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef WINDOWS +#include +#endif + + +namespace popl +{ + +#define POPL_VERSION "1.3.0" + + +/// Option's argument type +/** + * Switch has "no" argument + * Value has "required" argument + * Implicit has "optional" argument + */ +enum class Argument +{ + no = 0, // option never takes an argument + required, // option always requires an argument + optional // option may take an argument +}; + + +/// Option's attribute +/** + * inactive: Option is not set and will not be parsed + * hidden: Option is active, but will not show up in the help message + * required: Option must be set on the command line. Otherwise an exception will be thrown + * optional: Option must not be set. Default attribute. + * advanced: Option is advanced and will only show up in the advanced help message + * expoert: Option is expert and will only show up in the expert help message + */ +enum class Attribute +{ + inactive = 0, + hidden = 1, + required = 2, + optional = 3, + advanced = 4, + expert = 5 +}; + + +/// Option name type. Used in invalid_option exception. +/** + * unspecified: not specified + * short_name: The option's short name + * long_name: The option's long name + */ +enum class OptionName +{ + unspecified, + short_name, + long_name +}; + + +/// Abstract Base class for Options +/** + * Base class for Options + * holds just configuration data, no runtime data. + * Option is not bound to a special type "T" + */ +class Option +{ + friend class OptionParser; + +public: + /// Construct an Option + /// @param short_name the options's short name. Must be empty or one character. + /// @param long_name the option's long name. Can be empty. + /// @param description the Option's description that will be shown in the help message + Option(const std::string& short_name, const std::string& long_name, std::string description); + + /// Destructor + virtual ~Option() = default; + + /// default copy constructor + Option(const Option&) = default; + + /// default move constructor + Option(Option&&) = default; + + /// default assignement operator + Option& operator=(const Option&) = default; + + /// default move assignement operator + Option& operator=(Option&&) = default; + + /// Get the Option's short name + /// @return character of the options's short name or 0 if no short name is defined + char short_name() const; + + /// Get the Option's long name + /// @return the long name of the Option. Empty string if no long name is defined + std::string long_name() const; + + /// Get the Option's long or short name + /// @param what_name the option's name to return + /// @param what_hyphen preced the returned name with (double-)hypen + /// @return the requested name of the Option. Empty string if not defined. + std::string name(OptionName what_name, bool with_hypen = false) const; + + /// Get the Option's description + /// @return the description + std::string description() const; + + /// Get the Option's default value + /// @param out stream to write the default value to + /// @return true if a default value is available, false if not + virtual bool get_default(std::ostream& out) const = 0; + + /// Set the Option's attribute + /// @param attribute + void set_attribute(const Attribute& attribute); + + /// Get the Option's attribute + /// @return the Options's attribute + Attribute attribute() const; + + /// Get the Option's argument type + /// @return argument type (no, required, optional) + virtual Argument argument_type() const = 0; + + /// Check how often the Option is set on command line + /// @return the Option's count on command line + virtual size_t count() const = 0; + + /// Check if the Option is set + /// @return true if set at least once + virtual bool is_set() const = 0; + + +protected: + /// Parse the command line option and fill the internal data structure + /// @param what_name short or long option name + /// @param value the value as given on command line + virtual void parse(OptionName what_name, const char* value) = 0; + + /// Clear the internal data structure + virtual void clear() = 0; + + std::string short_name_; + std::string long_name_; + std::string description_; + Attribute attribute_; +}; + + + +/// Value option with optional default value +/** + * Value option with optional default value + * If set, it requires an argument + */ +template +class Value : public Option +{ +public: + /// Construct an Value Option + /// @param short_name the option's short name. Must be empty or one character. + /// @param long_name the option's long name. Can be empty. + /// @param description the Option's description that will be shown in the help message + Value(const std::string& short_name, const std::string& long_name, const std::string& description); + + /// Construct an Value Option + /// @param short_name the option's short name. Must be empty or one character. + /// @param long_name the option's long name. Can be empty. + /// @param description the Option's description that will be shown in the help message + /// @param default_val the Option's default value + /// @param assign_to pointer to a variable to assign the parsed command line value to + Value(const std::string& short_name, const std::string& long_name, const std::string& description, const T& default_val, T* assign_to = nullptr); + + size_t count() const override; + bool is_set() const override; + + /// Assign the last parsed command line value to "var" + /// @param var pointer to the variable where is value is written to + void assign_to(T* var); + + /// Manually set the Option's value. Deletes current value(s) + /// @param value the new value of the option + void set_value(const T& value); + + /// Get the Option's value. Will throw if option at index idx is not available + /// @param idx the zero based index of the value (if set multiple times) + /// @return the Option's value at index "idx" + T value(size_t idx = 0) const; + + /// Get the Option's value, return default_value if not set. + /// @param default_value return value if value is not set + /// @param idx the zero based index of the value (if set multiple times) + /// @return the Option's value at index "idx" or the default value or default_value + T value_or(const T& default_value, size_t idx = 0) const; + + /// Set the Option's default value + /// @param value the default value if not specified on command line + void set_default(const T& value); + + /// Check if the Option has a default value + /// @return true if the Option has a default value + bool has_default() const; + + /// Get the Option's default value. Will throw if no default is set. + /// @return the Option's default value + T get_default() const; + bool get_default(std::ostream& out) const override; + + Argument argument_type() const override; + +protected: + void parse(OptionName what_name, const char* value) override; + std::unique_ptr default_; + + virtual void update_reference(); + virtual void add_value(const T& value); + void clear() override; + + T* assign_to_; + std::vector values_; +}; + + + +/// Value option with implicit default value +/** + * Value option with implicit default value + * If set, an argument is optional + * -without argument it carries the implicit default value + * -with argument it carries the explicit value + */ +template +class Implicit : public Value +{ +public: + Implicit(const std::string& short_name, const std::string& long_name, const std::string& description, const T& implicit_val, T* assign_to = nullptr); + + Argument argument_type() const override; + +protected: + void parse(OptionName what_name, const char* value) override; +}; + + + +/// Value option without value +/** + * Value option without value + * Does not require an argument + * Can be either set or not set + */ +class Switch : public Value +{ +public: + Switch(const std::string& short_name, const std::string& long_name, const std::string& description, bool* assign_to = nullptr); + + void set_default(const bool& value) = delete; + Argument argument_type() const override; + +protected: + void parse(OptionName what_name, const char* value) override; +}; + + + +using Option_ptr = std::shared_ptr