From 488e4f2c532eda96757e593598e2acaa76f1a2e0 Mon Sep 17 00:00:00 2001 From: Victor Lin Date: Fri, 17 Jul 2026 15:06:13 +0800 Subject: [PATCH 1/3] safely handle invalid shapes, polygons, links, and graph connections build against modern MSVC and windows 11 --- CMakeLists.txt | 73 +++++++--- depthmapX/UI/AboutDlg.ui | 2 +- depthmapX/UI/licenseagreement.ui | 4 +- depthmapX/dialogs/AboutDlg.cpp | 2 +- depthmapX/version_defs.h | 29 ++++ depthmapx-master-msvc-build-fix.patch | 159 ++++++++++++++++++++++ genlib/pafmath.h | 6 +- releases/licenses.txt | 1 + salaTest/testshapemaps.cpp | 35 +++++ salalib/alllinemap.cpp | 82 +++++++---- salalib/axialminimiser.cpp | 30 ++-- salalib/axialmodules/axialintegration.cpp | 6 +- salalib/shapemap.cpp | 115 ++++++++++------ salalib/vgamodules/vgavisualglobal.cpp | 6 +- version.h | 2 +- 15 files changed, 440 insertions(+), 112 deletions(-) create mode 100644 depthmapX/version_defs.h create mode 100644 depthmapx-master-msvc-build-fix.patch diff --git a/CMakeLists.txt b/CMakeLists.txt index 6bf2539a..57c2234a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,14 +1,31 @@ +cmake_minimum_required(VERSION 3.13.0) set(CMAKE_OSX_DEPLOYMENT_TARGET "10.9" CACHE STRING "Minimum macOS deployment version") project(depthmapX) -cmake_minimum_required(VERSION 3.13.0) set(CMAKE_CXX_STANDARD 17) +option(DEPTHMAPX_ENABLE_IPO "Enable interprocedural optimisation for Release builds" OFF) + +if (DEPTHMAPX_ENABLE_IPO) + include(CheckIPOSupported) + check_ipo_supported(RESULT DEPTHMAPX_IPO_SUPPORTED OUTPUT DEPTHMAPX_IPO_ERROR) + if (DEPTHMAPX_IPO_SUPPORTED) + set(CMAKE_INTERPROCEDURAL_OPTIMIZATION_RELEASE ON) + else() + message(WARNING "IPO/LTO is not supported by this toolchain: ${DEPTHMAPX_IPO_ERROR}") + endif() +endif() + +if (MSVC) + # /MP improves build time. /Oi, /Gy and linker folding improve Release runtime/size. + add_compile_options(/MP "$<$:/Oi>" "$<$:/Gy>") + add_link_options("$<$:/OPT:REF>" "$<$:/OPT:ICF>") +endif() + if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") set(warnings "-Wall -Wextra") elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") set(warnings "/W4 /EHsc") - execute_process(COMMAND make_version_header.bat WORKING_DIRECTORY depthmapX) endif() # policy for target sources - we don't expect any old CMakes @@ -16,27 +33,39 @@ cmake_policy(SET CMP0076 NEW) include_directories(".") -# Get the current working branch -execute_process( - COMMAND git rev-parse --abbrev-ref HEAD - WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} - OUTPUT_VARIABLE APP_BRANCH - OUTPUT_STRIP_TRAILING_WHITESPACE -) - -# Get the latest abbreviated commit hash of the working branch -execute_process( - COMMAND git log -1 --format=%h - WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} - OUTPUT_VARIABLE APP_COMMIT - OUTPUT_STRIP_TRAILING_WHITESPACE -) - -# generate version_defs.h -include_directories(${CMAKE_BINARY_DIR}) -configure_file("${CMAKE_SOURCE_DIR}/version_defs.h.in" "${CMAKE_BINARY_DIR}/version_defs.h" @ONLY) +# Generate version metadata without relying on the deprecated Windows WMIC tool. +set(APP_BRANCH "unknown") +set(APP_COMMIT "unknown") +find_package(Git QUIET) +if (GIT_FOUND AND EXISTS "${CMAKE_SOURCE_DIR}/.git") + execute_process( + COMMAND "${GIT_EXECUTABLE}" rev-parse --abbrev-ref HEAD + WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}" + RESULT_VARIABLE APP_BRANCH_RESULT + OUTPUT_VARIABLE APP_BRANCH + OUTPUT_STRIP_TRAILING_WHITESPACE + ERROR_QUIET + ) + if (NOT APP_BRANCH_RESULT EQUAL 0 OR APP_BRANCH STREQUAL "") + set(APP_BRANCH "unknown") + endif() -string(TIMESTAMP APP_DATE %Y-%m-%d) + execute_process( + COMMAND "${GIT_EXECUTABLE}" log -1 --format=%h + WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}" + RESULT_VARIABLE APP_COMMIT_RESULT + OUTPUT_VARIABLE APP_COMMIT + OUTPUT_STRIP_TRAILING_WHITESPACE + ERROR_QUIET + ) + if (NOT APP_COMMIT_RESULT EQUAL 0 OR APP_COMMIT STREQUAL "") + set(APP_COMMIT "unknown") + endif() +endif() + +string(TIMESTAMP APP_DATE "%Y-%m-%d") +include_directories("${CMAKE_BINARY_DIR}") +configure_file("${CMAKE_SOURCE_DIR}/version_defs.h.in" "${CMAKE_BINARY_DIR}/version_defs.h" @ONLY) set(modules_core "" CACHE INTERNAL "modules_core" FORCE) set(MODULES_GUI FALSE) diff --git a/depthmapX/UI/AboutDlg.ui b/depthmapX/UI/AboutDlg.ui index 82f683e7..196119a1 100644 --- a/depthmapX/UI/AboutDlg.ui +++ b/depthmapX/UI/AboutDlg.ui @@ -78,7 +78,7 @@ - <html><head/><body><p><a href="https://github.com/SpaceGroupUCL/depthmapX"><span style=" text-decoration: underline; color:#0000ff;">https://github.com/varoudis/depthmapX</span></a></p></body></html> + <html><head/><body><p><a href="https://github.com/SpaceGroupUCL/depthmapX"><span style=" text-decoration: underline; color:#0000ff;">https://github.com/SpaceGroupUCL/depthmapX</span></a></p></body></html> diff --git a/depthmapX/UI/licenseagreement.ui b/depthmapX/UI/licenseagreement.ui index 4ae9ba1e..3e3ae01d 100644 --- a/depthmapX/UI/licenseagreement.ui +++ b/depthmapX/UI/licenseagreement.ui @@ -34,7 +34,7 @@ - <html><head/><body><p><span style=" font-size:24pt;">depthmapX</span></p><p><span style=" font-size:12pt;">Multi-Platform Spatial Network Analysis Software</span></p><p><a href="https://github.com/varoudis/depthmapX"><span style=" text-decoration: underline; color:#0000ff;">https://github.com/varoudis/depthmapX</span></a></p><p><span style=" font-size:10pt;">(C) 2000-2010 University College London, Alasdair Turner, Eva Friedrich<br/>(C) 2011-2014 Tasos Varoudis, UCL<br/>(C) 2017 Christian Sailer, Petros Koutsolampros</span></p><p><br/></p><p><span style=" font-size:12pt;">In memory of Alasdair Turner </span></p></body></html> + <html><head/><body><p><span style=" font-size:24pt;">depthmapX</span></p><p><span style=" font-size:12pt;">Multi-Platform Spatial Network Analysis Software</span></p><p><a href="https://github.com/SpaceGroupUCL/depthmapX"><span style=" text-decoration: underline; color:#0000ff;">https://github.com/SpaceGroupUCL/depthmapX</span></a></p><p><span style=" font-size:10pt;">(C) 2000-2010 University College London, Alasdair Turner, Eva Friedrich<br/>(C) 2011-2014 Tasos Varoudis, UCL<br/>(C) 2017 Christian Sailer, Petros Koutsolampros<br/>(C) 2026 Victor Yu-Chieh Lin</span></p><p><br/></p><p><span style=" font-size:12pt;">In memory of Alasdair Turner </span></p></body></html> @@ -47,7 +47,7 @@ <html><head><meta name="qrichtext" content="1" /><style type="text/css"> p, li { white-space: pre-wrap; } </style></head><body style=" font-family:'.SF NS Text'; font-size:13pt; font-weight:400; font-style:normal;"> -<p style=" margin-top:12px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><span style=" font-family:'Courier New,courier'; font-size:14pt; color:#000000;">(c) 2000-2010 University College London, Alasdair Turner, Eva Friedrich<br />(c) 2011-2014, Tasos Varoudis<br />(c) 2017 Christian Sailer, Petros Koutsolampros</span></p> +<p style=" margin-top:12px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><span style=" font-family:'Courier New,courier'; font-size:14pt; color:#000000;">(c) 2000-2010 University College London, Alasdair Turner, Eva Friedrich<br />(c) 2011-2014, Tasos Varoudis<br />(c) 2017 Christian Sailer, Petros Koutsolampros<br />(c) 2026 Victor Yu-Chieh Lin</span></p> <p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><span style=" font-family:'MS Shell Dlg 2'; font-size:8.25pt;"><br /></span></p> <p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><span style=" font-family:'Courier New,courier'; font-size:14pt; color:#000000;">This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.</span></p> <p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><span style=" font-family:'MS Shell Dlg 2'; font-size:8.25pt;"><br /></span></p> diff --git a/depthmapX/dialogs/AboutDlg.cpp b/depthmapX/dialogs/AboutDlg.cpp index 346b595f..06847b30 100644 --- a/depthmapX/dialogs/AboutDlg.cpp +++ b/depthmapX/dialogs/AboutDlg.cpp @@ -23,7 +23,7 @@ CAboutDlg::CAboutDlg(QWidget *parent) QString m_version_info; m_version_info = QString(tr("Version %1.%2.%3 (%4, %5)")).arg(DEPTHMAPX_MAJOR_VERSION).arg(DEPTHMAPX_MINOR_VERSION).arg(DEPTHMAPX_REVISION_VERSION).arg(APP_GIT_BRANCH).arg(APP_GIT_COMMIT); QString m_copyright; - m_copyright = QString(tr("(C) 2000-2010 University College London, Alasdair Turner, Eva Friedrich\n(C) 2011-2014 Tasos Varoudis\n(C) 2017 Christian Sailer, Petros Koutsolampros")); + m_copyright = QString(tr("(C) 2000-2010 University College London, Alasdair Turner, Eva Friedrich\n(C) 2011-2014 Tasos Varoudis\n(C) 2017 Christian Sailer, Petros Koutsolampros\n(C) 2026 Victor Yu-Chieh Lin ")); QString m_agreement; m_agreement = QString(tr("This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.\x0D\x0D\x0A\x0D\x0D\x0AThis program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.\x0D\x0D\x0A\x0D\x0D\x0AYou should have received a copy of the GNU General Public License along with this program. If not, see .")); diff --git a/depthmapX/version_defs.h b/depthmapX/version_defs.h new file mode 100644 index 00000000..043cfea2 --- /dev/null +++ b/depthmapX/version_defs.h @@ -0,0 +1,29 @@ +// Copyright (C) 2018 Christian Sailer +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. + +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +// This file is autogenerated - do not modify it directly! + +#pragma once + +#ifndef APP_DATE +#define APP_DATE "/00/00" +#endif + +#ifndef APP_GIT_BRANCH +#define APP_GIT_BRANCH "master" +#endif + +#ifndef APP_GIT_COMMIT +#define APP_GIT_COMMIT "02ceaec" +#endif diff --git a/depthmapx-master-msvc-build-fix.patch b/depthmapx-master-msvc-build-fix.patch new file mode 100644 index 00000000..504ed48e --- /dev/null +++ b/depthmapx-master-msvc-build-fix.patch @@ -0,0 +1,159 @@ +diff --git a/CMakeLists.txt b/CMakeLists.txt +index fc3d844..57c2234 100644 +--- a/CMakeLists.txt ++++ b/CMakeLists.txt +@@ -1,6 +1,6 @@ ++cmake_minimum_required(VERSION 3.13.0) + set(CMAKE_OSX_DEPLOYMENT_TARGET "10.9" CACHE STRING "Minimum macOS deployment version") + project(depthmapX) +-cmake_minimum_required(VERSION 3.13.0) + set(CMAKE_CXX_STANDARD 17) + + option(DEPTHMAPX_ENABLE_IPO "Enable interprocedural optimisation for Release builds" OFF) +@@ -26,7 +26,6 @@ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU" OR + set(warnings "-Wall -Wextra") + elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") + set(warnings "/W4 /EHsc") +- execute_process(COMMAND make_version_header.bat WORKING_DIRECTORY depthmapX) + endif() + + # policy for target sources - we don't expect any old CMakes +@@ -34,28 +33,40 @@ cmake_policy(SET CMP0076 NEW) + + include_directories(".") + +-# Get the current working branch +-execute_process( +- COMMAND git rev-parse --abbrev-ref HEAD +- WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} +- OUTPUT_VARIABLE APP_BRANCH +- OUTPUT_STRIP_TRAILING_WHITESPACE +-) ++# Generate version metadata without relying on the deprecated Windows WMIC tool. ++set(APP_BRANCH "unknown") ++set(APP_COMMIT "unknown") ++find_package(Git QUIET) ++if (GIT_FOUND AND EXISTS "${CMAKE_SOURCE_DIR}/.git") ++ execute_process( ++ COMMAND "${GIT_EXECUTABLE}" rev-parse --abbrev-ref HEAD ++ WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}" ++ RESULT_VARIABLE APP_BRANCH_RESULT ++ OUTPUT_VARIABLE APP_BRANCH ++ OUTPUT_STRIP_TRAILING_WHITESPACE ++ ERROR_QUIET ++ ) ++ if (NOT APP_BRANCH_RESULT EQUAL 0 OR APP_BRANCH STREQUAL "") ++ set(APP_BRANCH "unknown") ++ endif() + +-# Get the latest abbreviated commit hash of the working branch +-execute_process( +- COMMAND git log -1 --format=%h +- WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} +- OUTPUT_VARIABLE APP_COMMIT +- OUTPUT_STRIP_TRAILING_WHITESPACE +-) ++ execute_process( ++ COMMAND "${GIT_EXECUTABLE}" log -1 --format=%h ++ WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}" ++ RESULT_VARIABLE APP_COMMIT_RESULT ++ OUTPUT_VARIABLE APP_COMMIT ++ OUTPUT_STRIP_TRAILING_WHITESPACE ++ ERROR_QUIET ++ ) ++ if (NOT APP_COMMIT_RESULT EQUAL 0 OR APP_COMMIT STREQUAL "") ++ set(APP_COMMIT "unknown") ++ endif() ++endif() + +-# generate version_defs.h +-include_directories(${CMAKE_BINARY_DIR}) ++string(TIMESTAMP APP_DATE "%Y-%m-%d") ++include_directories("${CMAKE_BINARY_DIR}") + configure_file("${CMAKE_SOURCE_DIR}/version_defs.h.in" "${CMAKE_BINARY_DIR}/version_defs.h" @ONLY) + +-string(TIMESTAMP APP_DATE %Y-%m-%d) +- + set(modules_core "" CACHE INTERNAL "modules_core" FORCE) + set(MODULES_GUI FALSE) + set(MODULES_CLI FALSE) +diff --git a/genlib/pafmath.h b/genlib/pafmath.h +index 202df5f..20775c6 100644 +--- a/genlib/pafmath.h ++++ b/genlib/pafmath.h +@@ -58,8 +58,6 @@ inline double prandomr(int set = 0) { return double(pafrand(set)) / double(PAF_R + // note, in order to stop confusing myself I have ln defined: + #define ln(X) log(X) + +-inline double log2(double a) { return (ln(a) * M_1_LN2); } +- + // Hillier Hanson dvalue + /* + inline double dvalue(double k) +@@ -69,10 +67,10 @@ inline double dvalue(double k) + */ + + // Hillier Hanson dvalue (from Kruger 1989 -- see Teklenburg et al) +-inline double dvalue(double k) { return 2.0 * (k * (log2((k + 2.0) / 3.0) - 1.0) + 1.0) / ((k - 1.0) * (k - 2.0)); } ++inline double dvalue(double k) { return 2.0 * (k * (std::log2((k + 2.0) / 3.0) - 1.0) + 1.0) / ((k - 1.0) * (k - 2.0)); } + + // Hillier Hanson pvalue +-inline double pvalue(double k) { return 2.0 * (k - log2(k) - 1.0) / ((k - 1.0) * (k - 2.0)); } ++inline double pvalue(double k) { return 2.0 * (k - std::log2(k) - 1.0) / ((k - 1.0) * (k - 2.0)); } + + // Teklenburg integration (correction 31.01.11 due to Ulrich Thaler + inline double teklinteg(double nodecount, double totaldepth) { +diff --git a/salalib/axialmodules/axialintegration.cpp b/salalib/axialmodules/axialintegration.cpp +index ee6e4ac..c647bfd 100644 +--- a/salalib/axialmodules/axialintegration.cpp ++++ b/salalib/axialmodules/axialintegration.cpp +@@ -21,6 +21,8 @@ + #include "genlib/pflipper.h" + #include "genlib/stringutils.h" + ++#include ++ + bool AxialIntegration::run(Communicator *comm, ShapeGraph &map, bool simple_version) { + // note, from 10.0, Depthmap no longer includes *self* connections on axial lines + // self connections are stripped out on loading graph files, as well as no longer made +@@ -426,11 +428,11 @@ bool AxialIntegration::run(Communicator *comm, ShapeGraph &map, bool simple_vers + // some debate over whether or not this should be node count - 1 + // (i.e., including or not including the node itself) + double prob = double(depthcounts[k]) / double(node_count); +- entropy -= prob * log2(prob); ++ entropy -= prob * std::log2(prob); + // Formula from Turner 2001, "Depthmap" + factorial *= double(k + 1); + double q = (pow(mean_depth, double(k)) / double(factorial)) * exp(-mean_depth); +- rel_entropy += (double)prob * log2(prob / q); ++ rel_entropy += (double)prob * std::log2(prob / q); + // + harmonic += 1.0 / double(depthcounts[k]); + } +diff --git a/salalib/vgamodules/vgavisualglobal.cpp b/salalib/vgamodules/vgavisualglobal.cpp +index ed29f4d..734f0f8 100644 +--- a/salalib/vgamodules/vgavisualglobal.cpp ++++ b/salalib/vgamodules/vgavisualglobal.cpp +@@ -20,6 +20,8 @@ + + #include "genlib/stringutils.h" + ++#include ++ + bool VGAVisualGlobal::run(Communicator *comm, PointMap &map, bool simple_version) { + time_t atime = 0; + if (comm) { +@@ -173,11 +175,11 @@ bool VGAVisualGlobal::run(Communicator *comm, PointMap &map, bool simple_version + for (size_t k = 1; k < distribution.size(); k++) { + if (distribution[k] > 0) { + double prob = double(distribution[k]) / double(total_nodes - 1); +- entropy -= prob * log2(prob); ++ entropy -= prob * std::log2(prob); + // Formula from Turner 2001, "Depthmap" + factorial *= double(k + 1); + double q = (pow(mean_depth, double(k)) / double(factorial)) * exp(-mean_depth); +- rel_entropy += (float)prob * log2(prob / q); ++ rel_entropy += (float)prob * std::log2(prob / q); + } + } + if (!simple_version) { diff --git a/genlib/pafmath.h b/genlib/pafmath.h index 202df5f3..20775c6c 100644 --- a/genlib/pafmath.h +++ b/genlib/pafmath.h @@ -58,8 +58,6 @@ inline double prandomr(int set = 0) { return double(pafrand(set)) / double(PAF_R // note, in order to stop confusing myself I have ln defined: #define ln(X) log(X) -inline double log2(double a) { return (ln(a) * M_1_LN2); } - // Hillier Hanson dvalue /* inline double dvalue(double k) @@ -69,10 +67,10 @@ inline double dvalue(double k) */ // Hillier Hanson dvalue (from Kruger 1989 -- see Teklenburg et al) -inline double dvalue(double k) { return 2.0 * (k * (log2((k + 2.0) / 3.0) - 1.0) + 1.0) / ((k - 1.0) * (k - 2.0)); } +inline double dvalue(double k) { return 2.0 * (k * (std::log2((k + 2.0) / 3.0) - 1.0) + 1.0) / ((k - 1.0) * (k - 2.0)); } // Hillier Hanson pvalue -inline double pvalue(double k) { return 2.0 * (k - log2(k) - 1.0) / ((k - 1.0) * (k - 2.0)); } +inline double pvalue(double k) { return 2.0 * (k - std::log2(k) - 1.0) / ((k - 1.0) * (k - 2.0)); } // Teklenburg integration (correction 31.01.11 due to Ulrich Thaler inline double teklinteg(double nodecount, double totaldepth) { diff --git a/releases/licenses.txt b/releases/licenses.txt index b9bbe994..d5335cbe 100644 --- a/releases/licenses.txt +++ b/releases/licenses.txt @@ -2,6 +2,7 @@ depthmapX C 2000-2010 University College London, Alasdair Turner, Eva Friedrich C 2011-2014 Tasos Varoudis, UCL C 2017 Christian Sailer, Petros Koutsolampros +C 2026 Victor Yu-Chieh Lin Released under Gnu Public License (GPL v3) Please see gplv3.txt for license details diff --git a/salaTest/testshapemaps.cpp b/salaTest/testshapemaps.cpp index 987a980b..c50224cd 100644 --- a/salaTest/testshapemaps.cpp +++ b/salaTest/testshapemaps.cpp @@ -20,6 +20,7 @@ #include "catch.hpp" #include #include +#include TEST_CASE("Testing ShapeMap::getAllShapes variants") { @@ -121,3 +122,37 @@ TEST_CASE("Testing ShapeMap::getAllShapes variants") REQUIRE(colour.bluef() == Approx(0.2f).epsilon(EPSILON)); } } + +TEST_CASE("SalaShape handles empty and degenerate polylines") { + SECTION("Empty shape") { + SalaShape shape(SalaShape::SHAPE_POLY); + shape.setCentroidAreaPerim(); + REQUIRE(shape.getArea() == Approx(0.0)); + REQUIRE(shape.getPerimeter() == Approx(0.0)); + REQUIRE(std::isfinite(shape.getCentroid().x)); + REQUIRE(std::isfinite(shape.getCentroid().y)); + REQUIRE(shape.getAngDev() == Approx(0.0)); + } + + SECTION("Collinear points") { + SalaShape shape(SalaShape::SHAPE_POLY); + shape.m_points = {Point2f(0.0, 0.0), Point2f(2.0, 0.0), Point2f(4.0, 0.0)}; + shape.setCentroidAreaPerim(); + REQUIRE(shape.getArea() == Approx(0.0)); + REQUIRE(shape.getCentroid().x == Approx(2.0)); + REQUIRE(shape.getCentroid().y == Approx(0.0)); + } +} + +TEST_CASE("ShapeMap rejects stale polygon edit references") { + ShapeMap shapeMap("Test ShapeMap"); + const int shapeRef = shapeMap.polyBegin(Line(Point2f(0.0, 0.0), Point2f(1.0, 0.0))); + + REQUIRE_FALSE(shapeMap.polyAppend(shapeRef + 1000, Point2f(2.0, 0.0))); + REQUIRE_FALSE(shapeMap.polyClose(shapeRef + 1000)); + REQUIRE_FALSE(shapeMap.polyCancel(shapeRef + 1000)); + REQUIRE(shapeMap.getAllShapes().size() == 1); + + shapeMap.removeShape(shapeRef + 1000, false); + REQUIRE(shapeMap.getAllShapes().size() == 1); +} diff --git a/salalib/alllinemap.cpp b/salalib/alllinemap.cpp index 73830474..69eabd5e 100644 --- a/salalib/alllinemap.cpp +++ b/salalib/alllinemap.cpp @@ -4,6 +4,7 @@ #include "genlib/exceptions.h" #include #include +#include AllLineMap::AllLineMap(Communicator *comm, std::vector &drawingLayers, @@ -104,22 +105,38 @@ AllLineMap::AllLineMap(Communicator *comm, comm->CommPostMessage( Communicator::CURRENT_RECORD, 0 ); } - // cut out duplicates: - int removed = 0; // for testing purposes + // Mark duplicates first and compact once. Repeated vector erases added an + // O(n) shift to every duplicate and size() - 1 underflowed when empty. + std::vector duplicate(axiallines.size(), false); + const double maxdim = __max(region.width(), region.height()); for (size_t j = 0; j < axiallines.size(); j++) { - for (size_t k = axiallines.size() - 1; k > j; k--) { - double maxdim = __max(region.width(),region.height()); - if (approxeq(axiallines[j].start(), axiallines[k].start(), maxdim * TOLERANCE_B) && approxeq(axiallines[j].end(), axiallines[k].end(), maxdim * TOLERANCE_B)) { - for (int preaxiali: preaxialdata[k]) { - preaxialdata[j].insert(preaxiali); - } - preaxialdata.erase(preaxialdata.begin() + int(k)); - axiallines.erase(axiallines.begin() + int(k)); - removed++; + if (duplicate[j]) { + continue; + } + for (size_t k = j + 1; k < axiallines.size(); k++) { + if (!duplicate[k] && + approxeq(axiallines[j].start(), axiallines[k].start(), maxdim * TOLERANCE_B) && + approxeq(axiallines[j].end(), axiallines[k].end(), maxdim * TOLERANCE_B)) { + preaxialdata[j].insert(preaxialdata[k].begin(), preaxialdata[k].end()); + duplicate[k] = true; } } } + size_t writeIndex = 0; + for (size_t readIndex = 0; readIndex < axiallines.size(); readIndex++) { + if (duplicate[readIndex]) { + continue; + } + if (writeIndex != readIndex) { + axiallines[writeIndex] = std::move(axiallines[readIndex]); + preaxialdata[writeIndex] = std::move(preaxialdata[readIndex]); + } + ++writeIndex; + } + axiallines.resize(writeIndex); + preaxialdata.resize(writeIndex); + region.grow(0.99); // <- this paired with crop code below to prevent error init(axiallines.size(), m_polygons.getRegion()); // used to be double density here initialiseAttributesAxial(); @@ -296,23 +313,41 @@ void AllLineMap::makeDivisions(const std::vector& polyconnections comm->CommPostMessage( Communicator::NUM_RECORDS, polyconnections.size() ); } + // std::distance on a std::map is linear. Cache its stable map-order index + // once rather than walking from begin() for every polygon connector. + std::map radialIndices; + size_t radialIndex = 0; + for (const auto &entry : radialdivisions) { + radialIndices.emplace(entry.first, radialIndex++); + } + for (size_t i = 0; i < polyconnections.size(); i++) { PixelRefVector pixels = pixelateLine(polyconnections[i].line); - std::vector testedshapes; + std::unordered_set testedshapes; auto connIter = radialdivisions.find(polyconnections[i].key); - size_t connindex = std::distance(radialdivisions.begin(), connIter); + auto connIndexIter = radialIndices.find(polyconnections[i].key); + if (connIter == radialdivisions.end() || connIndexIter == radialIndices.end()) { + continue; + } + const size_t connindex = connIndexIter->second; + if (connindex >= radiallines.size()) { + continue; + } double tolerance = sqrt(TOLERANCE_A);// * polyconnections[i].line.length(); for (size_t j = 0; j < pixels.size(); j++) { PixelRef pix = pixels[j]; auto& shapes = m_pixel_shapes(static_cast(pix.y), static_cast(pix.x)); for (const ShapeRef& shape: shapes) { - auto iter = depthmapX::findBinary( testedshapes, shape.m_shape_ref ); - if (iter != testedshapes.end()) { + if (!testedshapes.insert(shape.m_shape_ref).second) { continue; } - testedshapes.insert(iter, int(shape.m_shape_ref)); - const Line& line = m_shapes.find(shape.m_shape_ref)->second.getLine(); + auto shapeIter = m_shapes.find(shape.m_shape_ref); + auto axialDividerIter = axialdividers.find(static_cast(shape.m_shape_ref)); + if (shapeIter == m_shapes.end() || axialDividerIter == axialdividers.end()) { + continue; + } + const Line& line = shapeIter->second.getLine(); // if (intersect_region(line, polyconnections[i].line, tolerance * line.length()) ) { switch ( intersect_line_distinguish(line, polyconnections[i].line, tolerance * line.length()) ) { @@ -320,24 +355,15 @@ void AllLineMap::makeDivisions(const std::vector& polyconnections break; case 2: { - size_t index = depthmapX::findIndexFromKey(axialdividers, (int) shape.m_shape_ref); - if (index != shape.m_shape_ref) { - throw 1; // for the code to work later this can't be true! - } - axialdividers[index].insert(connindex); + axialDividerIter->second.insert(static_cast(connindex)); connIter->second.insert(shape.m_shape_ref); } break; case 1: { - size_t index = depthmapX::findIndexFromKey(axialdividers, (int) shape.m_shape_ref); - if (index != shape.m_shape_ref) { - throw 1; // for the code to work later this can't be true! - } - // // this makes sure actually crosses between the line and the openspace properly if (radiallines[connindex].cuts(line)) { - axialdividers[index].insert(connindex); + axialDividerIter->second.insert(static_cast(connindex)); connIter->second.insert(shape.m_shape_ref); } } diff --git a/salalib/axialminimiser.cpp b/salalib/axialminimiser.cpp index f82665b6..1d91ab11 100644 --- a/salalib/axialminimiser.cpp +++ b/salalib/axialminimiser.cpp @@ -91,6 +91,9 @@ void AxialMinimiser::removeSubsets(std::map >& axsegcuts, std bool subset = false; for (size_t j = 0; j < axa.m_connections.size(); j++) { int indextob = axa.m_connections[j]; + if (indextob < 0 || static_cast(indextob) >= m_axialconns.size()) { + continue; + } if (indextob == ii || m_removed[indextob]) { // <- removed[indextob] should never happen as it should have been removed below continue; } @@ -101,22 +104,27 @@ void AxialMinimiser::removeSubsets(std::map >& axsegcuts, std // first check it's a connection subset // note that changes in 10.08 mean that lines no longer connect to themselves // this means that the subset 1 connects {2,3} and 2 connects {1,3} are equivalent - for (size_t axai = 0, axbi = 0; axai < axa.m_connections.size() && axbi < axb.m_connections.size(); axai++, axbi++) { - // extra 10.08 -> step over connection to b + size_t axai = 0; + size_t axbi = 0; + while (axai < axa.m_connections.size() && axbi < axb.m_connections.size()) { + // Step over the implicit a<->b connection on each side before comparing. if (axa.m_connections[axai] == indextob) { - axai++; + ++axai; + continue; } - // extra 10.08 add axb.m_connections[axbi] == ii -> step over connection to a - while (axbi < axb.m_connections.size() && (axb.m_connections[axbi] == ii || axa.m_connections[axai] > axb.m_connections[axbi])) { - axbi++; + if (axb.m_connections[axbi] == static_cast(ii)) { + ++axbi; + continue; } - if (axbi >= axb.m_connections.size()) { - break; + if (axa.m_connections[axai] == axb.m_connections[axbi]) { + ++coconnecting; + ++axai; + ++axbi; } - else if (axa.m_connections[axai] == axb.m_connections[axbi]) { - coconnecting++; + else if (axa.m_connections[axai] > axb.m_connections[axbi]) { + ++axbi; } - else if (axa.m_connections[axai] < axb.m_connections[axbi]) { + else { break; } } diff --git a/salalib/axialmodules/axialintegration.cpp b/salalib/axialmodules/axialintegration.cpp index ee6e4acc..c647bfdb 100644 --- a/salalib/axialmodules/axialintegration.cpp +++ b/salalib/axialmodules/axialintegration.cpp @@ -21,6 +21,8 @@ #include "genlib/pflipper.h" #include "genlib/stringutils.h" +#include + bool AxialIntegration::run(Communicator *comm, ShapeGraph &map, bool simple_version) { // note, from 10.0, Depthmap no longer includes *self* connections on axial lines // self connections are stripped out on loading graph files, as well as no longer made @@ -426,11 +428,11 @@ bool AxialIntegration::run(Communicator *comm, ShapeGraph &map, bool simple_vers // some debate over whether or not this should be node count - 1 // (i.e., including or not including the node itself) double prob = double(depthcounts[k]) / double(node_count); - entropy -= prob * log2(prob); + entropy -= prob * std::log2(prob); // Formula from Turner 2001, "Depthmap" factorial *= double(k + 1); double q = (pow(mean_depth, double(k)) / double(factorial)) * exp(-mean_depth); - rel_entropy += (double)prob * log2(prob / q); + rel_entropy += (double)prob * std::log2(prob / q); // harmonic += 1.0 / double(depthcounts[k]); } diff --git a/salalib/shapemap.cpp b/salalib/shapemap.cpp index fb486374..ed444071 100644 --- a/salalib/shapemap.cpp +++ b/salalib/shapemap.cpp @@ -79,6 +79,10 @@ void SalaShape::setCentroidAreaPerim() { m_area = 0.0; m_perimeter = 0.0; m_centroid = Point2f(0, 0); + if (m_points.empty()) { + return; + } + for (size_t i = 0; i < m_points.size(); i++) { Point2f &p1 = m_points[i]; Point2f &p2 = m_points[(i + 1) % m_points.size()]; @@ -90,12 +94,25 @@ void SalaShape::setCentroidAreaPerim() { Point2f side = p2 - p1; m_perimeter += side.length(); } + const double signedArea = m_area; m_type &= ~SHAPE_CCW; - if (sgn(m_area) == 1) { + if (sgn(signedArea) == 1) { m_type |= SHAPE_CCW; } - m_centroid.scale(2.0 / m_area); // note, *not* fabs(m_area) as it is then confused by clockwise ordered shapes - m_area = fabs(m_area); + + if (std::isfinite(signedArea) && fabs(signedArea) > TOLERANCE_A) { + // Do not use fabs here: clockwise polygons need the signed-area direction. + m_centroid.scale(2.0 / signedArea); + } else { + // Repeated or collinear points have zero area. Keep a finite centroid + // instead of introducing Inf/NaN into rendering and spatial indexing. + m_centroid = Point2f(0, 0); + for (const Point2f &point : m_points) { + m_centroid += point; + } + m_centroid /= static_cast(m_points.size()); + } + m_area = std::isfinite(signedArea) ? fabs(signedArea) : 0.0; if (isOpen()) { // take off the automatically collected final side Point2f side = m_points.back() - m_points.front(); @@ -108,6 +125,10 @@ void SalaShape::setCentroid(const Point2f &p) { m_centroid = p; } // get the angular deviation along the length of a poly line: double SalaShape::getAngDev() const { + if (m_points.size() < 3) { + return 0.0; + } + double dev = 0.0; for (size_t i = 1; i < m_points.size() - 1; i++) { double ang = angle(m_points[i - 1], m_points[i], m_points[i + 1]); @@ -223,7 +244,7 @@ void ShapeMap::copy(const ShapeMap &sourcemap, int copyflags) { } } - if ((copyflags & ShapeMap::COPY_ATTRIBUTES) == ShapeMap::COPY_GRAPH) { + if ((copyflags & ShapeMap::COPY_GRAPH) == ShapeMap::COPY_GRAPH) { if (sourcemap.m_hasgraph) { m_hasgraph = true; // straight copy: @@ -813,8 +834,11 @@ int ShapeMap::polyBegin(const Line &line) { } bool ShapeMap::polyAppend(int shape_ref, const Point2f &point) { - // don't do anything too complex: - SalaShape &shape = m_shapes.rbegin()->second; + auto shapeIter = m_shapes.find(shape_ref); + if (shapeIter == m_shapes.end()) { + return false; + } + SalaShape &shape = shapeIter->second; // check you can actually do this first if (!(shape.isLine() || shape.isPolyLine())) { @@ -855,8 +879,11 @@ bool ShapeMap::polyAppend(int shape_ref, const Point2f &point) { } bool ShapeMap::polyClose(int shape_ref) { - // don't do anything too complex: - SalaShape &shape = m_shapes.rbegin()->second; + auto shapeIter = m_shapes.find(shape_ref); + if (shapeIter == m_shapes.end()) { + return false; + } + SalaShape &shape = shapeIter->second; // check you can actually do this first if (!shape.isPolyLine()) { @@ -874,14 +901,22 @@ bool ShapeMap::polyClose(int shape_ref) { } bool ShapeMap::polyCancel(int shape_ref) { - // don't do anything too complex: - SalaShape &shape = m_shapes.rbegin()->second; + auto shapeIter = m_shapes.find(shape_ref); + if (shapeIter == m_shapes.end()) { + return false; + } + SalaShape &shape = shapeIter->second; // check you can actually do this first if (!(shape.isLine() || shape.isPolyLine())) { return false; } + if (m_undobuffer.empty() || + m_undobuffer.back().m_action != SalaEvent::SALA_CREATED || + m_undobuffer.back().m_shape_ref != shape_ref) { + return false; + } m_undobuffer.pop_back(); removeShape(shape_ref, true); @@ -913,12 +948,20 @@ bool ShapeMap::removeSelected() { } void ShapeMap::removeShape(int shaperef, bool undoing) { + auto shapeIter = m_shapes.find(shaperef); + if (shapeIter == m_shapes.end()) { + return; + } + + const size_t rowid = static_cast(std::distance(m_shapes.begin(), shapeIter)); + if (m_hasgraph && m_connectors.size() != m_shapes.size()) { + // Connector indices are positional. Avoid partially mutating a corrupt map. + return; + } + // remove shape from four keys: the pixel grid, the poly list, the attributes and the connections removePolyPixels(shaperef); // done first, as all interface references use this list - auto shapeIter = m_shapes.find(shaperef); - size_t rowid = std::distance(m_shapes.begin(), shapeIter); - if (!undoing) { // <- if not currently undoing another event, then add to the undo buffer: m_undobuffer.push_back(SalaEvent(SalaEvent::SALA_DELETED, shaperef)); m_undobuffer.back().m_geometry = shapeIter->second; @@ -933,11 +976,11 @@ void ShapeMap::removeShape(int shaperef, bool undoing) { int conn_col = m_attributes->getColumnIndex("Connectivity"); // TODO: Replace with iterators - for (size_t i = m_connectors.size() - 1; static_cast(i) != -1; i--) { + for (size_t i = m_connectors.size(); i-- > 0;) { if (i == rowid) { continue; // it's going to be removed anyway } - for (size_t j = m_connectors[i].m_connections.size() - 1; static_cast(j) != -1; j--) { + for (size_t j = m_connectors[i].m_connections.size(); j-- > 0;) { if (m_connectors[i].m_connections[j] == int(rowid)) { m_connectors[i].m_connections.erase(m_connectors[i].m_connections.begin() + int(j)); if (conn_col != -1) { @@ -953,32 +996,28 @@ void ShapeMap::removeShape(int shaperef, bool undoing) { m_connectors.erase(m_connectors.begin() + int(rowid)); - // take out explicit links and unlinks (note, undo won't restore these): - for (auto revIter = m_links.rbegin(); revIter != m_links.rend(); ++revIter) { - if (revIter->a == static_cast(rowid) || revIter->b == static_cast(rowid)) { - m_links.erase(std::next(revIter).base()); - } else { - if (revIter->a > int(rowid)) - revIter->a -= 1; - if (revIter->b > int(rowid)) - revIter->b -= 1; - } - } - for (auto revIter = m_unlinks.rbegin(); revIter != m_unlinks.rend(); ++revIter) { - if (revIter->a == static_cast(rowid) || revIter->b == static_cast(rowid)) { - m_unlinks.erase(std::next(revIter).base()); - } else { - if (revIter->a > int(rowid)) - revIter->a -= 1; - if (revIter->b > int(rowid)) - revIter->b -= 1; + // Take out explicit links and unlinks (note, undo won't restore these). + // Erasing through a reverse iterator invalidated the iterator used by the old loop. + auto removeAndReindex = [rowid](std::vector &pairs) { + for (auto iter = pairs.begin(); iter != pairs.end();) { + if (iter->a == static_cast(rowid) || iter->b == static_cast(rowid)) { + iter = pairs.erase(iter); + continue; + } + if (iter->a > static_cast(rowid)) { + --iter->a; + } + if (iter->b > static_cast(rowid)) { + --iter->b; + } + ++iter; } - } + }; + removeAndReindex(m_links); + removeAndReindex(m_unlinks); } - if (shapeIter != m_shapes.end()) { - shapeIter = m_shapes.erase(shapeIter); - } + m_shapes.erase(shapeIter); // n.b., shaperef should have been used to create the row in the first place: const AttributeKey shapeRefKey(shaperef); m_attributes->removeRow(shapeRefKey); diff --git a/salalib/vgamodules/vgavisualglobal.cpp b/salalib/vgamodules/vgavisualglobal.cpp index ed29f4dd..734f0f8b 100644 --- a/salalib/vgamodules/vgavisualglobal.cpp +++ b/salalib/vgamodules/vgavisualglobal.cpp @@ -20,6 +20,8 @@ #include "genlib/stringutils.h" +#include + bool VGAVisualGlobal::run(Communicator *comm, PointMap &map, bool simple_version) { time_t atime = 0; if (comm) { @@ -173,11 +175,11 @@ bool VGAVisualGlobal::run(Communicator *comm, PointMap &map, bool simple_version for (size_t k = 1; k < distribution.size(); k++) { if (distribution[k] > 0) { double prob = double(distribution[k]) / double(total_nodes - 1); - entropy -= prob * log2(prob); + entropy -= prob * std::log2(prob); // Formula from Turner 2001, "Depthmap" factorial *= double(k + 1); double q = (pow(mean_depth, double(k)) / double(factorial)) * exp(-mean_depth); - rel_entropy += (float)prob * log2(prob / q); + rel_entropy += (float)prob * std::log2(prob / q); } } if (!simple_version) { diff --git a/version.h b/version.h index 34843804..0297fcd6 100644 --- a/version.h +++ b/version.h @@ -22,7 +22,7 @@ // use these to define the depthmap versions #define DEPTHMAPX_MAJOR_VERSION 0 #define DEPTHMAPX_MINOR_VERSION 8 -#define DEPTHMAPX_REVISION_VERSION 0 +#define DEPTHMAPX_REVISION_VERSION 1 #define DEPTHMAP_MODULE_VERSION 10.04 From 59b6ce02cfde29c6093d2699b724da2c37a4d160 Mon Sep 17 00:00:00 2001 From: Victor Lin Date: Sat, 18 Jul 2026 18:39:38 +0800 Subject: [PATCH 2/3] Multithreaded analyses, optimizations --- CMakeLists.txt | 1 + depthmapX/GraphDoc.cpp | 2 +- depthmapX/UI/SegmentAnalysisDlg.ui | 4 +- depthmapx-master-msvc-build-fix.patch | 159 ------- genlib/comm.h | 9 +- genlib/parallel.h | 138 +++++++ salaTest/CMakeLists.txt | 1 + salaTest/testparallelanalysis.cpp | 334 +++++++++++++++ salalib/CMakeLists.txt | 1 + salalib/alllinemap.cpp | 287 +++++++++---- salalib/axialmap.cpp | 322 +++++++++------ salalib/axialmodules/axialintegration.cpp | 480 ++++++++++++---------- salalib/segmmodules/segmangular.cpp | 186 +++++---- salalib/segmmodules/segmmetric.cpp | 308 +++++++------- salalib/segmmodules/segmtopological.cpp | 312 +++++++------- salalib/shapemap.cpp | 31 +- salalib/shapemap.h | 5 +- version.h | 2 +- 18 files changed, 1590 insertions(+), 992 deletions(-) delete mode 100644 depthmapx-master-msvc-build-fix.patch create mode 100644 genlib/parallel.h create mode 100644 salaTest/testparallelanalysis.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 57c2234a..71f5b916 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,6 +37,7 @@ include_directories(".") set(APP_BRANCH "unknown") set(APP_COMMIT "unknown") find_package(Git QUIET) +find_package(Threads REQUIRED) if (GIT_FOUND AND EXISTS "${CMAKE_SOURCE_DIR}/.git") execute_process( COMMAND "${GIT_EXECUTABLE}" rev-parse --abbrev-ref HEAD diff --git a/depthmapX/GraphDoc.cpp b/depthmapX/GraphDoc.cpp index c58d54b7..978f27f6 100644 --- a/depthmapX/GraphDoc.cpp +++ b/depthmapX/GraphDoc.cpp @@ -292,7 +292,7 @@ void QGraphDoc::cancel_wait() else { m_communicator->Cancel(); } - // Don't cancel --- cancel should be handled by the thread! + // Don't cancel --- cancel should be handled by the thread! } } diff --git a/depthmapX/UI/SegmentAnalysisDlg.ui b/depthmapX/UI/SegmentAnalysisDlg.ui index cc9b50c0..6ccf5252 100644 --- a/depthmapX/UI/SegmentAnalysisDlg.ui +++ b/depthmapX/UI/SegmentAnalysisDlg.ui @@ -23,7 +23,7 @@ - Tulip Analysis (Faster) + Tulip Analysis @@ -65,7 +65,7 @@ - Full Angular (Slower) + Full Angular (Multi-Threaded) diff --git a/depthmapx-master-msvc-build-fix.patch b/depthmapx-master-msvc-build-fix.patch deleted file mode 100644 index 504ed48e..00000000 --- a/depthmapx-master-msvc-build-fix.patch +++ /dev/null @@ -1,159 +0,0 @@ -diff --git a/CMakeLists.txt b/CMakeLists.txt -index fc3d844..57c2234 100644 ---- a/CMakeLists.txt -+++ b/CMakeLists.txt -@@ -1,6 +1,6 @@ -+cmake_minimum_required(VERSION 3.13.0) - set(CMAKE_OSX_DEPLOYMENT_TARGET "10.9" CACHE STRING "Minimum macOS deployment version") - project(depthmapX) --cmake_minimum_required(VERSION 3.13.0) - set(CMAKE_CXX_STANDARD 17) - - option(DEPTHMAPX_ENABLE_IPO "Enable interprocedural optimisation for Release builds" OFF) -@@ -26,7 +26,6 @@ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU" OR - set(warnings "-Wall -Wextra") - elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") - set(warnings "/W4 /EHsc") -- execute_process(COMMAND make_version_header.bat WORKING_DIRECTORY depthmapX) - endif() - - # policy for target sources - we don't expect any old CMakes -@@ -34,28 +33,40 @@ cmake_policy(SET CMP0076 NEW) - - include_directories(".") - --# Get the current working branch --execute_process( -- COMMAND git rev-parse --abbrev-ref HEAD -- WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} -- OUTPUT_VARIABLE APP_BRANCH -- OUTPUT_STRIP_TRAILING_WHITESPACE --) -+# Generate version metadata without relying on the deprecated Windows WMIC tool. -+set(APP_BRANCH "unknown") -+set(APP_COMMIT "unknown") -+find_package(Git QUIET) -+if (GIT_FOUND AND EXISTS "${CMAKE_SOURCE_DIR}/.git") -+ execute_process( -+ COMMAND "${GIT_EXECUTABLE}" rev-parse --abbrev-ref HEAD -+ WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}" -+ RESULT_VARIABLE APP_BRANCH_RESULT -+ OUTPUT_VARIABLE APP_BRANCH -+ OUTPUT_STRIP_TRAILING_WHITESPACE -+ ERROR_QUIET -+ ) -+ if (NOT APP_BRANCH_RESULT EQUAL 0 OR APP_BRANCH STREQUAL "") -+ set(APP_BRANCH "unknown") -+ endif() - --# Get the latest abbreviated commit hash of the working branch --execute_process( -- COMMAND git log -1 --format=%h -- WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} -- OUTPUT_VARIABLE APP_COMMIT -- OUTPUT_STRIP_TRAILING_WHITESPACE --) -+ execute_process( -+ COMMAND "${GIT_EXECUTABLE}" log -1 --format=%h -+ WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}" -+ RESULT_VARIABLE APP_COMMIT_RESULT -+ OUTPUT_VARIABLE APP_COMMIT -+ OUTPUT_STRIP_TRAILING_WHITESPACE -+ ERROR_QUIET -+ ) -+ if (NOT APP_COMMIT_RESULT EQUAL 0 OR APP_COMMIT STREQUAL "") -+ set(APP_COMMIT "unknown") -+ endif() -+endif() - --# generate version_defs.h --include_directories(${CMAKE_BINARY_DIR}) -+string(TIMESTAMP APP_DATE "%Y-%m-%d") -+include_directories("${CMAKE_BINARY_DIR}") - configure_file("${CMAKE_SOURCE_DIR}/version_defs.h.in" "${CMAKE_BINARY_DIR}/version_defs.h" @ONLY) - --string(TIMESTAMP APP_DATE %Y-%m-%d) -- - set(modules_core "" CACHE INTERNAL "modules_core" FORCE) - set(MODULES_GUI FALSE) - set(MODULES_CLI FALSE) -diff --git a/genlib/pafmath.h b/genlib/pafmath.h -index 202df5f..20775c6 100644 ---- a/genlib/pafmath.h -+++ b/genlib/pafmath.h -@@ -58,8 +58,6 @@ inline double prandomr(int set = 0) { return double(pafrand(set)) / double(PAF_R - // note, in order to stop confusing myself I have ln defined: - #define ln(X) log(X) - --inline double log2(double a) { return (ln(a) * M_1_LN2); } -- - // Hillier Hanson dvalue - /* - inline double dvalue(double k) -@@ -69,10 +67,10 @@ inline double dvalue(double k) - */ - - // Hillier Hanson dvalue (from Kruger 1989 -- see Teklenburg et al) --inline double dvalue(double k) { return 2.0 * (k * (log2((k + 2.0) / 3.0) - 1.0) + 1.0) / ((k - 1.0) * (k - 2.0)); } -+inline double dvalue(double k) { return 2.0 * (k * (std::log2((k + 2.0) / 3.0) - 1.0) + 1.0) / ((k - 1.0) * (k - 2.0)); } - - // Hillier Hanson pvalue --inline double pvalue(double k) { return 2.0 * (k - log2(k) - 1.0) / ((k - 1.0) * (k - 2.0)); } -+inline double pvalue(double k) { return 2.0 * (k - std::log2(k) - 1.0) / ((k - 1.0) * (k - 2.0)); } - - // Teklenburg integration (correction 31.01.11 due to Ulrich Thaler - inline double teklinteg(double nodecount, double totaldepth) { -diff --git a/salalib/axialmodules/axialintegration.cpp b/salalib/axialmodules/axialintegration.cpp -index ee6e4ac..c647bfd 100644 ---- a/salalib/axialmodules/axialintegration.cpp -+++ b/salalib/axialmodules/axialintegration.cpp -@@ -21,6 +21,8 @@ - #include "genlib/pflipper.h" - #include "genlib/stringutils.h" - -+#include -+ - bool AxialIntegration::run(Communicator *comm, ShapeGraph &map, bool simple_version) { - // note, from 10.0, Depthmap no longer includes *self* connections on axial lines - // self connections are stripped out on loading graph files, as well as no longer made -@@ -426,11 +428,11 @@ bool AxialIntegration::run(Communicator *comm, ShapeGraph &map, bool simple_vers - // some debate over whether or not this should be node count - 1 - // (i.e., including or not including the node itself) - double prob = double(depthcounts[k]) / double(node_count); -- entropy -= prob * log2(prob); -+ entropy -= prob * std::log2(prob); - // Formula from Turner 2001, "Depthmap" - factorial *= double(k + 1); - double q = (pow(mean_depth, double(k)) / double(factorial)) * exp(-mean_depth); -- rel_entropy += (double)prob * log2(prob / q); -+ rel_entropy += (double)prob * std::log2(prob / q); - // - harmonic += 1.0 / double(depthcounts[k]); - } -diff --git a/salalib/vgamodules/vgavisualglobal.cpp b/salalib/vgamodules/vgavisualglobal.cpp -index ed29f4d..734f0f8 100644 ---- a/salalib/vgamodules/vgavisualglobal.cpp -+++ b/salalib/vgamodules/vgavisualglobal.cpp -@@ -20,6 +20,8 @@ - - #include "genlib/stringutils.h" - -+#include -+ - bool VGAVisualGlobal::run(Communicator *comm, PointMap &map, bool simple_version) { - time_t atime = 0; - if (comm) { -@@ -173,11 +175,11 @@ bool VGAVisualGlobal::run(Communicator *comm, PointMap &map, bool simple_version - for (size_t k = 1; k < distribution.size(); k++) { - if (distribution[k] > 0) { - double prob = double(distribution[k]) / double(total_nodes - 1); -- entropy -= prob * log2(prob); -+ entropy -= prob * std::log2(prob); - // Formula from Turner 2001, "Depthmap" - factorial *= double(k + 1); - double q = (pow(mean_depth, double(k)) / double(factorial)) * exp(-mean_depth); -- rel_entropy += (float)prob * log2(prob / q); -+ rel_entropy += (float)prob * std::log2(prob / q); - } - } - if (!simple_version) { diff --git a/genlib/comm.h b/genlib/comm.h index 7885c7fc..70a16a09 100644 --- a/genlib/comm.h +++ b/genlib/comm.h @@ -17,6 +17,7 @@ #pragma once //#include +#include #include #include #include @@ -66,7 +67,7 @@ class Communicator { enum { NUM_STEPS, CURRENT_STEP, NUM_RECORDS, CURRENT_RECORD }; protected: - bool m_cancelled; + std::atomic_bool m_cancelled; bool m_delete_flag; // nb. converted to Win32 UTF-16 Unicode path (AT 31.01.11) Linux, MacOS use UTF-8 (AT 29.04.11) std::string m_infilename; @@ -80,7 +81,7 @@ class Communicator { m_infile = NULL; m_infile2 = NULL; m_outfile = NULL; - m_cancelled = false; + m_cancelled.store(false, std::memory_order_relaxed); m_delete_flag = false; } // @@ -130,8 +131,8 @@ class Communicator { } void SetOutfile(const char *filename) { m_outfile = new std::ofstream(filename); } // - bool IsCancelled() const { return m_cancelled; } - void Cancel() { m_cancelled = true; } + bool IsCancelled() const { return m_cancelled.load(std::memory_order_acquire); } + void Cancel() { m_cancelled.store(true, std::memory_order_release); } // std::ifstream &getInFileStream() { return *m_infile; } std::ifstream &GetInfile2() { return *m_infile2; } diff --git a/genlib/parallel.h b/genlib/parallel.h new file mode 100644 index 00000000..17641810 --- /dev/null +++ b/genlib/parallel.h @@ -0,0 +1,138 @@ +// genlib - a component of the depthmapX - spatial network analysis platform +// Copyright (C) 2026 depthmapX contributors + +#pragma once + +#include "genlib/comm.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace depthmapX { + +inline size_t configuredAnalysisThreadCount(size_t taskCount, size_t estimatedBytesPerWorker = 0) { + size_t requested = 0; + if (const char *value = std::getenv("DEPTHMAPX_THREADS")) { + try { + requested = static_cast(std::stoul(value)); + } catch (...) { + requested = 0; + } + } + + if (requested == 0 && taskCount < 256) { + return 1; + } + + size_t hardwareThreads = static_cast(std::thread::hardware_concurrency()); + if (hardwareThreads == 0) { + hardwareThreads = 2; + } + + size_t workers = requested > 0 ? std::min(requested, hardwareThreads) + : std::min(hardwareThreads, 8); + workers = std::max(1, std::min(workers, taskCount)); + + // Keep per-analysis scratch allocations under a conservative 512 MiB budget by default. + if (estimatedBytesPerWorker > 0) { + constexpr size_t MEMORY_BUDGET = size_t(512) * 1024 * 1024; + const size_t memoryLimitedWorkers = std::max(1, MEMORY_BUDGET / estimatedBytesPerWorker); + workers = std::min(workers, memoryLimitedWorkers); + } + + return std::max(1, workers); +} + +template +void parallelFor(size_t taskCount, size_t workerCount, Communicator *comm, Function &&function) { + if (taskCount == 0) { + return; + } + + workerCount = std::max(1, std::min(workerCount, taskCount)); + if (workerCount == 1) { + for (size_t task = 0; task < taskCount; ++task) { + if (comm && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + function(task, size_t(0)); + if (comm && (((task + 1) & 0x3f) == 0 || task + 1 == taskCount)) { + comm->CommPostMessage(Communicator::CURRENT_RECORD, static_cast(task + 1)); + } + } + return; + } + + std::atomic completed{0}; + std::atomic_bool stop{false}; + std::exception_ptr workerException; + std::mutex exceptionMutex; + + std::vector workers; + workers.reserve(workerCount); + for (size_t worker = 0; worker < workerCount; ++worker) { + workers.emplace_back([&, worker]() { + // Static cyclic scheduling keeps reductions repeatable while still spreading + // neighbouring roots across the available workers. + for (size_t task = worker; task < taskCount && !stop.load(std::memory_order_relaxed); + task += workerCount) { + if (comm && comm->IsCancelled()) { + stop.store(true, std::memory_order_relaxed); + break; + } + + try { + function(task, worker); + completed.fetch_add(1, std::memory_order_release); + } catch (...) { + { + std::lock_guard lock(exceptionMutex); + if (!workerException) { + workerException = std::current_exception(); + } + } + stop.store(true, std::memory_order_relaxed); + break; + } + } + }); + } + + size_t lastReported = std::numeric_limits::max(); + while (!stop.load(std::memory_order_relaxed) && completed.load(std::memory_order_acquire) < taskCount) { + if (comm && comm->IsCancelled()) { + stop.store(true, std::memory_order_relaxed); + break; + } + const size_t current = completed.load(std::memory_order_acquire); + if (comm && current != lastReported) { + comm->CommPostMessage(Communicator::CURRENT_RECORD, static_cast(current)); + lastReported = current; + } + std::this_thread::sleep_for(std::chrono::milliseconds(25)); + } + + for (std::thread &worker : workers) { + worker.join(); + } + + if (workerException) { + std::rethrow_exception(workerException); + } + if (comm && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + if (comm) { + comm->CommPostMessage(Communicator::CURRENT_RECORD, static_cast(completed.load())); + } +} + +} // namespace depthmapX diff --git a/salaTest/CMakeLists.txt b/salaTest/CMakeLists.txt index 13f30707..57fe33c1 100644 --- a/salaTest/CMakeLists.txt +++ b/salaTest/CMakeLists.txt @@ -26,6 +26,7 @@ set(salaTest_SRCS testpointinpoly.cpp testpushvalues.cpp testisovist.cpp + testparallelanalysis.cpp ) # salaTest_SRCS include_directories("../ThirdParty/Catch" "../ThirdParty/FakeIt") diff --git a/salaTest/testparallelanalysis.cpp b/salaTest/testparallelanalysis.cpp new file mode 100644 index 00000000..f523603e --- /dev/null +++ b/salaTest/testparallelanalysis.cpp @@ -0,0 +1,334 @@ +// Copyright (C) 2026 depthmapX contributors + +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. + +#include "catch.hpp" + +#include "genlib/parallel.h" +#include "salalib/alllinemap.h" +#include "salalib/axialmodules/axialintegration.h" +#include "salalib/mapconverter.h" +#include "salalib/mgraph.h" +#include "salalib/segmmodules/segmangular.h" +#include "salalib/segmmodules/segmmetric.h" +#include "salalib/segmmodules/segmtopological.h" + +#include +#include +#include +#include +#include + +namespace { + +void setAnalysisThreads(const char *value) { +#ifdef _WIN32 + _putenv_s("DEPTHMAPX_THREADS", value ? value : ""); +#else + if (value) { + setenv("DEPTHMAPX_THREADS", value, 1); + } else { + unsetenv("DEPTHMAPX_THREADS"); + } +#endif +} + +class AnalysisThreadSetting { + public: + explicit AnalysisThreadSetting(const char *value) { + const char *current = std::getenv("DEPTHMAPX_THREADS"); + if (current) { + m_previous = current; + m_hadPrevious = true; + } + setAnalysisThreads(value); + } + + ~AnalysisThreadSetting() { setAnalysisThreads(m_hadPrevious ? m_previous.c_str() : nullptr); } + + private: + std::string m_previous; + bool m_hadPrevious = false; +}; + +std::unique_ptr makeAxialGraph() { + std::unique_ptr metaGraph(new MetaGraph("Parallel axial test")); + metaGraph->m_drawingFiles.emplace_back("Drawing"); + metaGraph->m_drawingFiles.back().m_spacePixels.emplace_back("Lines"); + ShapeMap &drawing = metaGraph->m_drawingFiles.back().m_spacePixels.back(); + + drawing.makeLineShape(Line(Point2f(0.0, 0.0), Point2f(6.0, 0.0))); + drawing.makeLineShape(Line(Point2f(1.0, -1.0), Point2f(1.0, 1.0))); + drawing.makeLineShape(Line(Point2f(2.0, -1.0), Point2f(2.0, 1.0))); + drawing.makeLineShape(Line(Point2f(3.0, -1.0), Point2f(3.0, 1.0))); + drawing.makeLineShape(Line(Point2f(4.0, -1.0), Point2f(4.0, 1.0))); + drawing.makeLineShape(Line(Point2f(5.0, -1.0), Point2f(5.0, 1.0))); + drawing.makeLineShape(Line(Point2f(1.0, 1.0), Point2f(5.0, 1.0))); + drawing.makeLineShape(Line(Point2f(1.0, -1.0), Point2f(5.0, -1.0))); + + return MapConverter::convertDrawingToAxial(nullptr, "Axial", metaGraph->m_drawingFiles); +} + + +std::unique_ptr makeAllLineMap() { + std::vector drawingLayers; + drawingLayers.emplace_back("Drawing"); + drawingLayers.back().m_spacePixels.emplace_back("Boundaries"); + ShapeMap &drawing = drawingLayers.back().m_spacePixels.back(); + + drawing.makeLineShape(Line(Point2f(0.0, 0.0), Point2f(12.0, 0.0))); + drawing.makeLineShape(Line(Point2f(12.0, 0.0), Point2f(12.0, 10.0))); + drawing.makeLineShape(Line(Point2f(12.0, 10.0), Point2f(0.0, 10.0))); + drawing.makeLineShape(Line(Point2f(0.0, 10.0), Point2f(0.0, 0.0))); + + drawing.makeLineShape(Line(Point2f(4.0, 3.0), Point2f(8.0, 3.0))); + drawing.makeLineShape(Line(Point2f(8.0, 3.0), Point2f(8.0, 7.0))); + drawing.makeLineShape(Line(Point2f(8.0, 7.0), Point2f(4.0, 7.0))); + drawing.makeLineShape(Line(Point2f(4.0, 7.0), Point2f(4.0, 3.0))); + + return std::unique_ptr( + new AllLineMap(nullptr, drawingLayers, Point2f(1.0, 1.0), "All lines")); +} + +void requireSameGraphGeometryAndConnections(const ShapeGraph &serial, const ShapeGraph ¶llel) { + REQUIRE(serial.getShapeCount() == parallel.getShapeCount()); + REQUIRE(serial.getConnections().size() == parallel.getConnections().size()); + + auto serialShape = serial.getAllShapes().begin(); + auto parallelShape = parallel.getAllShapes().begin(); + for (; serialShape != serial.getAllShapes().end(); ++serialShape, ++parallelShape) { + REQUIRE(serialShape->first == parallelShape->first); + const Line &serialLine = serialShape->second.getLine(); + const Line ¶llelLine = parallelShape->second.getLine(); + REQUIRE(parallelLine.start().x == Approx(serialLine.start().x)); + REQUIRE(parallelLine.start().y == Approx(serialLine.start().y)); + REQUIRE(parallelLine.end().x == Approx(serialLine.end().x)); + REQUIRE(parallelLine.end().y == Approx(serialLine.end().y)); + } + + for (size_t index = 0; index < serial.getConnections().size(); ++index) { + REQUIRE(serial.getConnections()[index].m_connections == + parallel.getConnections()[index].m_connections); + } +} + +std::unique_ptr makeSegmentGraph() { + std::unique_ptr metaGraph(new MetaGraph("Parallel segment test")); + metaGraph->m_drawingFiles.emplace_back("Drawing"); + metaGraph->m_drawingFiles.back().m_spacePixels.emplace_back("Lines"); + ShapeMap &drawing = metaGraph->m_drawingFiles.back().m_spacePixels.back(); + + drawing.makeLineShape(Line(Point2f(0.0, 0.0), Point2f(1.0, 0.0))); + drawing.makeLineShape(Line(Point2f(1.0, 0.0), Point2f(2.0, 0.0))); + drawing.makeLineShape(Line(Point2f(2.0, 0.0), Point2f(3.0, 1.0))); + drawing.makeLineShape(Line(Point2f(2.0, 0.0), Point2f(3.0, -1.0))); + drawing.makeLineShape(Line(Point2f(3.0, 1.0), Point2f(4.0, 1.0))); + drawing.makeLineShape(Line(Point2f(3.0, -1.0), Point2f(4.0, -1.0))); + drawing.makeLineShape(Line(Point2f(4.0, 1.0), Point2f(5.0, 0.0))); + drawing.makeLineShape(Line(Point2f(4.0, -1.0), Point2f(5.0, 0.0))); + + return MapConverter::convertDrawingToSegment(nullptr, "Segment", metaGraph->m_drawingFiles); +} + +void requireSameSegmentConnections(const ShapeGraph &serial, const ShapeGraph ¶llel) { + const std::vector &serialConnections = serial.getConnections(); + const std::vector ¶llelConnections = parallel.getConnections(); + REQUIRE(serialConnections.size() == parallelConnections.size()); + + auto compareDirectional = [](const std::map &expected, + const std::map &actual) { + REQUIRE(expected.size() == actual.size()); + auto expectedIter = expected.begin(); + auto actualIter = actual.begin(); + for (; expectedIter != expected.end(); ++expectedIter, ++actualIter) { + REQUIRE(expectedIter->first.ref == actualIter->first.ref); + REQUIRE(expectedIter->first.dir == actualIter->first.dir); + REQUIRE(actualIter->second == Approx(expectedIter->second)); + } + }; + + for (size_t index = 0; index < serialConnections.size(); ++index) { + REQUIRE(serialConnections[index].m_segment_axialref == + parallelConnections[index].m_segment_axialref); + compareDirectional(serialConnections[index].m_back_segconns, + parallelConnections[index].m_back_segconns); + compareDirectional(serialConnections[index].m_forward_segconns, + parallelConnections[index].m_forward_segconns); + } +} + +void requireSameAttributes(const ShapeGraph &serial, const ShapeGraph ¶llel) { + const AttributeTable &serialAttributes = serial.getAttributeTable(); + const AttributeTable ¶llelAttributes = parallel.getAttributeTable(); + REQUIRE(serial.getShapeCount() == parallel.getShapeCount()); + REQUIRE(serialAttributes.getNumColumns() == parallelAttributes.getNumColumns()); + + for (size_t column = 0; column < serialAttributes.getNumColumns(); ++column) { + INFO("Column: " << serialAttributes.getColumnName(column)); + REQUIRE(serialAttributes.getColumnName(column) == parallelAttributes.getColumnName(column)); + for (size_t row = 0; row < serial.getShapeCount(); ++row) { + const float serialValue = serial.getAttributeRowFromShapeIndex(row).getValue(column); + const float parallelValue = parallel.getAttributeRowFromShapeIndex(row).getValue(column); + if (std::isfinite(serialValue) && std::isfinite(parallelValue)) { + REQUIRE(parallelValue == Approx(serialValue).epsilon(1e-5)); + } else { + REQUIRE(std::isfinite(serialValue) == std::isfinite(parallelValue)); + } + } + } +} + +void requireAngularIntegration(const ShapeGraph &graph, const std::vector &radiusSuffixes) { + const AttributeTable &attributes = graph.getAttributeTable(); + for (const std::string &suffix : radiusSuffixes) { + const int countColumn = attributes.getColumnIndex(std::string("Angular Node Count") + suffix); + const int totalDepthColumn = attributes.getColumnIndex(std::string("Angular Total Depth") + suffix); + const int integrationColumn = attributes.getColumnIndex(std::string("Angular Integration") + suffix); + REQUIRE(countColumn != -1); + REQUIRE(totalDepthColumn != -1); + REQUIRE(integrationColumn != -1); + + for (size_t rowIndex = 0; rowIndex < graph.getShapeCount(); ++rowIndex) { + const AttributeRow &row = graph.getAttributeRowFromShapeIndex(rowIndex); + const double count = row.getValue(countColumn); + const double totalDepth = row.getValue(totalDepthColumn); + const double integration = row.getValue(integrationColumn); + if (totalDepth > 1e-9) { + REQUIRE(integration == Approx(count * count / totalDepth).epsilon(1e-5)); + } else { + REQUIRE(integration == Approx(-1.0)); + } + } + } +} + +} // namespace + +TEST_CASE("Parallel axial analysis matches single-thread analysis") { + std::unique_ptr serial = makeAxialGraph(); + std::unique_ptr parallel = makeAxialGraph(); + + { + AnalysisThreadSetting threads("1"); + REQUIRE(AxialIntegration({2.0, -1.0}, -1, true, true, true).run(nullptr, *serial, false)); + } + { + AnalysisThreadSetting threads("4"); + REQUIRE(AxialIntegration({2.0, -1.0}, -1, true, true, true).run(nullptr, *parallel, false)); + } + + requireSameAttributes(*serial, *parallel); +} + +TEST_CASE("Parallel segment analyses match single-thread analyses") { + std::unique_ptr serial = makeSegmentGraph(); + std::unique_ptr parallel = makeSegmentGraph(); + + { + AnalysisThreadSetting threads("1"); + REQUIRE(SegmentAngular({1.0, -1.0}).run(nullptr, *serial, false)); + REQUIRE(SegmentMetric(-1.0, false).run(nullptr, *serial, false)); + REQUIRE(SegmentTopological(-1.0, false).run(nullptr, *serial, false)); + } + { + AnalysisThreadSetting threads("4"); + REQUIRE(SegmentAngular({1.0, -1.0}).run(nullptr, *parallel, false)); + REQUIRE(SegmentMetric(-1.0, false).run(nullptr, *parallel, false)); + REQUIRE(SegmentTopological(-1.0, false).run(nullptr, *parallel, false)); + } + + requireSameAttributes(*serial, *parallel); + requireAngularIntegration(*serial, {" R1.00", ""}); + requireAngularIntegration(*parallel, {" R1.00", ""}); +} + +TEST_CASE("Parallel line connection discovery matches single-thread discovery") { + std::unique_ptr serial; + std::unique_ptr parallel; + { + AnalysisThreadSetting threads("1"); + serial = makeAxialGraph(); + } + { + AnalysisThreadSetting threads("4"); + parallel = makeAxialGraph(); + } + requireSameGraphGeometryAndConnections(*serial, *parallel); + requireSameAttributes(*serial, *parallel); +} + +TEST_CASE("Parallel All-Line Stage 1 matches single-thread map creation") { + std::unique_ptr serial; + std::unique_ptr parallel; + { + AnalysisThreadSetting threads("1"); + serial = makeAllLineMap(); + } + { + AnalysisThreadSetting threads("4"); + parallel = makeAllLineMap(); + } + requireSameGraphGeometryAndConnections(*serial, *parallel); + requireSameAttributes(*serial, *parallel); + + std::tuple, std::unique_ptr> serialFewest; + std::tuple, std::unique_ptr> parallelFewest; + { + AnalysisThreadSetting threads("1"); + serialFewest = serial->extractFewestLineMaps(nullptr); + } + { + AnalysisThreadSetting threads("4"); + parallelFewest = parallel->extractFewestLineMaps(nullptr); + } + + requireSameGraphGeometryAndConnections(*std::get<0>(serialFewest), + *std::get<0>(parallelFewest)); + requireSameGraphGeometryAndConnections(*std::get<1>(serialFewest), + *std::get<1>(parallelFewest)); +} + +TEST_CASE("Parallel direct segment-map creation matches single-thread creation") { + std::unique_ptr serial; + std::unique_ptr parallel; + { + AnalysisThreadSetting threads("1"); + serial = makeSegmentGraph(); + } + { + AnalysisThreadSetting threads("4"); + parallel = makeSegmentGraph(); + } + + requireSameGraphGeometryAndConnections(*serial, *parallel); + requireSameSegmentConnections(*serial, *parallel); + requireSameAttributes(*serial, *parallel); +} + +TEST_CASE("Parallel axial-to-segment conversion matches single-thread conversion") { + std::unique_ptr serialAxial; + std::unique_ptr parallelAxial; + { + AnalysisThreadSetting threads("1"); + serialAxial = makeAxialGraph(); + parallelAxial = makeAxialGraph(); + } + + std::unique_ptr serial; + std::unique_ptr parallel; + { + AnalysisThreadSetting threads("1"); + serial = MapConverter::convertAxialToSegment(nullptr, *serialAxial, "Segment", true, false, 0.0); + } + { + AnalysisThreadSetting threads("4"); + parallel = MapConverter::convertAxialToSegment(nullptr, *parallelAxial, "Segment", true, false, 0.0); + } + + requireSameGraphGeometryAndConnections(*serial, *parallel); + requireSameSegmentConnections(*serial, *parallel); + requireSameAttributes(*serial, *parallel); +} diff --git a/salalib/CMakeLists.txt b/salalib/CMakeLists.txt index c4c435a7..2370bd7b 100644 --- a/salalib/CMakeLists.txt +++ b/salalib/CMakeLists.txt @@ -33,6 +33,7 @@ set(salalib_SRCS add_compile_definitions(_DEPTHMAP SALALIB_LIBRARY) add_library(${salalib} STATIC ${salalib_SRCS} ${vgamodules_SRCS} ${axialmodules_SRCS} ${segmmodules_SRCS} ${parsers_SRCS}) +target_link_libraries(${salalib} PUBLIC Threads::Threads) add_subdirectory(vgamodules) add_subdirectory(axialmodules) diff --git a/salalib/alllinemap.cpp b/salalib/alllinemap.cpp index 69eabd5e..2acd13d0 100644 --- a/salalib/alllinemap.cpp +++ b/salalib/alllinemap.cpp @@ -2,10 +2,148 @@ #include "salalib/axialminimiser.h" #include "salalib/tolerances.h" #include "genlib/exceptions.h" +#include "genlib/parallel.h" #include +#include +#include +#include #include +#include +#include #include +namespace { + +bool axialLinesMatch(const Line &a, const Line &b, double tolerance) { + return approxeq(a.start(), b.start(), tolerance) && approxeq(a.end(), b.end(), tolerance); +} + +void compactDuplicateAxialLinesSerial(std::vector &axiallines, KeyVertices &preaxialdata, + double tolerance) { + std::vector duplicate(axiallines.size(), false); + for (size_t j = 0; j < axiallines.size(); ++j) { + if (duplicate[j]) { + continue; + } + for (size_t k = j + 1; k < axiallines.size(); ++k) { + if (!duplicate[k] && axialLinesMatch(axiallines[j], axiallines[k], tolerance)) { + preaxialdata[j].insert(preaxialdata[k].begin(), preaxialdata[k].end()); + duplicate[k] = true; + } + } + } + + size_t writeIndex = 0; + for (size_t readIndex = 0; readIndex < axiallines.size(); ++readIndex) { + if (duplicate[readIndex]) { + continue; + } + if (writeIndex != readIndex) { + axiallines[writeIndex] = std::move(axiallines[readIndex]); + preaxialdata[writeIndex] = std::move(preaxialdata[readIndex]); + } + ++writeIndex; + } + axiallines.resize(writeIndex); + preaxialdata.resize(writeIndex); +} + +void compactDuplicateAxialLines(std::vector &axiallines, KeyVertices &preaxialdata, + double tolerance) { + const size_t lineCount = axiallines.size(); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(lineCount); + if (workerCount == 1 || lineCount < 512 || + lineCount > static_cast(std::numeric_limits::max())) { + compactDuplicateAxialLinesSerial(axiallines, preaxialdata, tolerance); + return; + } + + using DuplicatePair = std::pair; + constexpr size_t MAX_STORED_PAIRS = size_t(4) * 1024 * 1024; + const size_t perWorkerLimit = std::max(1, MAX_STORED_PAIRS / workerCount); + std::vector> workerPairs(workerCount); + std::atomic_bool overflow{false}; + + depthmapX::parallelFor(lineCount, workerCount, nullptr, [&](size_t j, size_t worker) { + if (overflow.load(std::memory_order_relaxed)) { + return; + } + std::vector &pairs = workerPairs[worker]; + for (size_t k = j + 1; k < lineCount; ++k) { + if (axialLinesMatch(axiallines[j], axiallines[k], tolerance)) { + if (pairs.size() >= perWorkerLimit) { + overflow.store(true, std::memory_order_relaxed); + return; + } + pairs.emplace_back(static_cast(j), static_cast(k)); + } + } + }); + + if (overflow.load(std::memory_order_relaxed)) { + // Pathological duplicate-heavy inputs would require too much temporary + // memory. Fall back to the exact serial algorithm instead. + compactDuplicateAxialLinesSerial(axiallines, preaxialdata, tolerance); + return; + } + + size_t pairCount = 0; + for (const auto &pairs : workerPairs) { + pairCount += pairs.size(); + } + std::vector pairs; + pairs.reserve(pairCount); + for (auto &workerResult : workerPairs) { + pairs.insert(pairs.end(), workerResult.begin(), workerResult.end()); + } + std::sort(pairs.begin(), pairs.end()); + + // Apply the matches in the same j/k order as the original serial loop. This + // preserves which representative survives even though comparisons ran in parallel. + std::vector duplicate(lineCount, false); + size_t pairIndex = 0; + for (size_t j = 0; j < lineCount; ++j) { + while (pairIndex < pairs.size() && pairs[pairIndex].first < j) { + ++pairIndex; + } + size_t groupEnd = pairIndex; + while (groupEnd < pairs.size() && pairs[groupEnd].first == j) { + ++groupEnd; + } + if (!duplicate[j]) { + for (size_t index = pairIndex; index < groupEnd; ++index) { + const size_t k = pairs[index].second; + if (!duplicate[k]) { + preaxialdata[j].insert(preaxialdata[k].begin(), preaxialdata[k].end()); + duplicate[k] = true; + } + } + } + pairIndex = groupEnd; + } + + size_t writeIndex = 0; + for (size_t readIndex = 0; readIndex < lineCount; ++readIndex) { + if (duplicate[readIndex]) { + continue; + } + if (writeIndex != readIndex) { + axiallines[writeIndex] = std::move(axiallines[readIndex]); + preaxialdata[writeIndex] = std::move(preaxialdata[readIndex]); + } + ++writeIndex; + } + axiallines.resize(writeIndex); + preaxialdata.resize(writeIndex); +} + +struct DivisionTaskResult { + size_t radialIndex = std::numeric_limits::max(); + std::vector axialShapeRefs; +}; + +} // namespace + AllLineMap::AllLineMap(Communicator *comm, std::vector &drawingLayers, const Point2f& seed, @@ -105,37 +243,11 @@ AllLineMap::AllLineMap(Communicator *comm, comm->CommPostMessage( Communicator::CURRENT_RECORD, 0 ); } - // Mark duplicates first and compact once. Repeated vector erases added an - // O(n) shift to every duplicate and size() - 1 underflowed when empty. - std::vector duplicate(axiallines.size(), false); + // Duplicate comparisons are independent and can be distributed across + // workers. Their matches are merged in the original serial order so the + // surviving line and key-vertex data remain deterministic. const double maxdim = __max(region.width(), region.height()); - for (size_t j = 0; j < axiallines.size(); j++) { - if (duplicate[j]) { - continue; - } - for (size_t k = j + 1; k < axiallines.size(); k++) { - if (!duplicate[k] && - approxeq(axiallines[j].start(), axiallines[k].start(), maxdim * TOLERANCE_B) && - approxeq(axiallines[j].end(), axiallines[k].end(), maxdim * TOLERANCE_B)) { - preaxialdata[j].insert(preaxialdata[k].begin(), preaxialdata[k].end()); - duplicate[k] = true; - } - } - } - - size_t writeIndex = 0; - for (size_t readIndex = 0; readIndex < axiallines.size(); readIndex++) { - if (duplicate[readIndex]) { - continue; - } - if (writeIndex != readIndex) { - axiallines[writeIndex] = std::move(axiallines[readIndex]); - preaxialdata[writeIndex] = std::move(preaxialdata[readIndex]); - } - ++writeIndex; - } - axiallines.resize(writeIndex); - preaxialdata.resize(writeIndex); + compactDuplicateAxialLines(axiallines, preaxialdata, maxdim * TOLERANCE_B); region.grow(0.99); // <- this paired with crop code below to prevent error init(axiallines.size(), m_polygons.getRegion()); // used to be double density here @@ -307,10 +419,8 @@ void AllLineMap::makeDivisions(const std::vector& polyconnections std::map >& radialdivisions, std::map > &axialdividers, Communicator *comm) { - time_t atime = 0; if (comm) { - qtimer( atime, 0 ); - comm->CommPostMessage( Communicator::NUM_RECORDS, polyconnections.size() ); + comm->CommPostMessage(Communicator::NUM_RECORDS, static_cast(polyconnections.size())); } // std::distance on a std::map is linear. Cache its stable map-order index @@ -321,66 +431,75 @@ void AllLineMap::makeDivisions(const std::vector& polyconnections radialIndices.emplace(entry.first, radialIndex++); } - for (size_t i = 0; i < polyconnections.size(); i++) { - PixelRefVector pixels = pixelateLine(polyconnections[i].line); - std::unordered_set testedshapes; - auto connIter = radialdivisions.find(polyconnections[i].key); - auto connIndexIter = radialIndices.find(polyconnections[i].key); - if (connIter == radialdivisions.end() || connIndexIter == radialIndices.end()) { - continue; - } - const size_t connindex = connIndexIter->second; - if (connindex >= radiallines.size()) { - continue; + const auto &shapesByRef = m_shapes; + const auto &pixelShapes = m_pixel_shapes; + const AllLineMap &readOnlyMap = *this; + std::vector taskResults(polyconnections.size()); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(polyconnections.size()); + + depthmapX::parallelFor(polyconnections.size(), workerCount, comm, [&](size_t i, size_t) { + const auto connIndexIter = radialIndices.find(polyconnections[i].key); + if (connIndexIter == radialIndices.end() || connIndexIter->second >= radiallines.size()) { + return; } - double tolerance = sqrt(TOLERANCE_A);// * polyconnections[i].line.length(); - for (size_t j = 0; j < pixels.size(); j++) { - PixelRef pix = pixels[j]; - auto& shapes = m_pixel_shapes(static_cast(pix.y), - static_cast(pix.x)); - for (const ShapeRef& shape: shapes) { - if (!testedshapes.insert(shape.m_shape_ref).second) { + + DivisionTaskResult &result = taskResults[i]; + result.radialIndex = connIndexIter->second; + const Line &connectionLine = polyconnections[i].line; + const double tolerance = sqrt(TOLERANCE_A); + const PixelRefVector pixels = readOnlyMap.pixelateLine(connectionLine); + std::unordered_set testedShapes; + + for (const PixelRef &pixel : pixels) { + const auto &shapes = pixelShapes(static_cast(pixel.y), + static_cast(pixel.x)); + for (const ShapeRef &shape : shapes) { + if (!testedShapes.insert(shape.m_shape_ref).second) { + continue; + } + + const auto shapeIter = shapesByRef.find(static_cast(shape.m_shape_ref)); + if (shapeIter == shapesByRef.end()) { continue; } - auto shapeIter = m_shapes.find(shape.m_shape_ref); - auto axialDividerIter = axialdividers.find(static_cast(shape.m_shape_ref)); - if (shapeIter == m_shapes.end() || axialDividerIter == axialdividers.end()) { + const Line &line = shapeIter->second.getLine(); + const double lineTolerance = tolerance * line.length(); + if (!intersect_region(line, connectionLine, lineTolerance)) { continue; } - const Line& line = shapeIter->second.getLine(); - // - if (intersect_region(line, polyconnections[i].line, tolerance * line.length()) ) { - switch ( intersect_line_distinguish(line, polyconnections[i].line, tolerance * line.length()) ) { - case 0: - break; - case 2: - { - axialDividerIter->second.insert(static_cast(connindex)); - connIter->second.insert(shape.m_shape_ref); - } - break; - case 1: - { - // this makes sure actually crosses between the line and the openspace properly - if (radiallines[connindex].cuts(line)) { - axialDividerIter->second.insert(static_cast(connindex)); - connIter->second.insert(shape.m_shape_ref); - } - } - break; - default: - break; - } + + const int intersectionType = intersect_line_distinguish(line, connectionLine, lineTolerance); + if (intersectionType == 2 || + (intersectionType == 1 && radiallines[result.radialIndex].cuts(line))) { + result.axialShapeRefs.push_back(static_cast(shape.m_shape_ref)); } } } - if (comm) { - if (qtimer( atime, 500 )) { - if (comm->IsCancelled()) { - throw Communicator::CancelledException(); - } - comm->CommPostMessage( Communicator::CURRENT_RECORD, i ); + + std::sort(result.axialShapeRefs.begin(), result.axialShapeRefs.end()); + result.axialShapeRefs.erase( + std::unique(result.axialShapeRefs.begin(), result.axialShapeRefs.end()), + result.axialShapeRefs.end()); + }); + + // Commit results in polygon-connection order. Only this phase mutates the + // maps, keeping the worker phase read-only and deterministic. + for (size_t i = 0; i < taskResults.size(); ++i) { + const DivisionTaskResult &result = taskResults[i]; + if (result.radialIndex == std::numeric_limits::max()) { + continue; + } + auto radialDivisionIter = radialdivisions.find(polyconnections[i].key); + if (radialDivisionIter == radialdivisions.end()) { + continue; + } + for (int shapeRef : result.axialShapeRefs) { + auto axialDividerIter = axialdividers.find(shapeRef); + if (axialDividerIter == axialdividers.end()) { + continue; } + axialDividerIter->second.insert(static_cast(result.radialIndex)); + radialDivisionIter->second.insert(shapeRef); } } } diff --git a/salalib/axialmap.cpp b/salalib/axialmap.cpp index 9b3e114e..a02ae9b6 100644 --- a/salalib/axialmap.cpp +++ b/salalib/axialmap.cpp @@ -30,10 +30,12 @@ #include "genlib/containerutils.h" #include "genlib/readwritehelpers.h" #include "genlib/pflipper.h" +#include "genlib/parallel.h" #include #include #include +#include #ifndef _WIN32 #define _finite finite @@ -66,23 +68,45 @@ void ShapeGraph::makeConnections(const KeyVertices &keyvertices) m_keyvertices.clear(); // note, expects these to be numbered 0, 1... - int conn_col = m_attributes->getColumnIndex("Connectivity"); - int leng_col = m_attributes->getColumnIndex("Line Length"); + const int conn_col = m_attributes->getColumnIndex("Connectivity"); + const int leng_col = m_attributes->getColumnIndex("Line Length"); + const size_t shapeCount = m_shapes.size(); + const double tolerance = TOLERANCE_B * __max(m_region.height(), m_region.width()); + + std::vector shapeKeys; + std::vector lineLengths; + shapeKeys.reserve(shapeCount); + lineLengths.reserve(shapeCount); + + std::unordered_map shapeIndexLookup; + shapeIndexLookup.reserve(shapeCount); + int shapeIndex = 0; + for (const auto &shape : m_shapes) { + shapeKeys.push_back(shape.first); + lineLengths.push_back(static_cast(shape.second.getLine().length())); + shapeIndexLookup.emplace(shape.first, shapeIndex++); + } - int i = -1; - for (auto shape: m_shapes) { - i++; - int key = shape.first; - AttributeRow &row = - m_attributes->getRow(AttributeKey(key)); - // all indices should match... - m_connectors.push_back( Connector() ); - m_connectors[i].m_connections = getLineConnections( key, TOLERANCE_B*__max(m_region.height(),m_region.width())); - row.setValue(conn_col, float(m_connectors[i].m_connections.size()) ); - row.setValue(leng_col, float(shape.second.getLine().length()) ); - if (keyvertices.size()) { + // Connection discovery only reads the completed geometry and spatial index. + // Store each line's result in its own connector, then update attributes serially. + m_connectors.resize(shapeCount); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(shapeCount); + depthmapX::parallelFor(shapeCount, workerCount, nullptr, [&](size_t index, size_t) { + m_connectors[index].m_connections = + getLineConnections(shapeKeys[index], tolerance, &shapeIndexLookup); + }); + + if (!keyvertices.empty()) { + m_keyvertices.resize(shapeCount); + } + + for (size_t index = 0; index < shapeCount; ++index) { + AttributeRow &row = m_attributes->getRow(AttributeKey(shapeKeys[index])); + row.setValue(conn_col, static_cast(m_connectors[index].m_connections.size())); + row.setValue(leng_col, lineLengths[index]); + if (!keyvertices.empty()) { // note: depends on lines being recorded in same order as keyvertices... - m_keyvertices.push_back( keyvertices[i] ); + m_keyvertices[index] = keyvertices[index]; } } @@ -440,114 +464,111 @@ void ShapeGraph::unlinkFromShapeMap(const ShapeMap& shapemap) // Method 1: direct linkage of endpoints where they touch void ShapeGraph::makeNewSegMap(Communicator *comm) { - // now make a connection set from the ends of lines: - struct LineConnector { - const Line &m_line; - Connector &m_connector; - int m_index; - LineConnector(const Line &line, Connector &connector, int index) - : m_line(line), m_connector(connector), m_index(index) {} + struct EndpointConnection { + size_t m_a; + size_t m_b; + bool m_aForward; + bool m_bForward; + int m_aDirection; + int m_bDirection; + float m_angle; }; - std::vector connectionset; - for (auto &shape : m_shapes) { - if (shape.second.isLine()) { - connectionset.emplace_back(); - } - } - std::map lineConnectors; - auto connectionIter = connectionset.begin(); - int connectionIdx = 0; - for (auto &shape : m_shapes) { + std::vector segmentLines; + segmentLines.reserve(m_shapes.size()); + std::unordered_map lineIndexLookup; + lineIndexLookup.reserve(m_shapes.size()); + + for (const auto &shape : m_shapes) { if (shape.second.isLine()) { - lineConnectors.insert(std::make_pair( - shape.first, LineConnector(shape.second.getLine(), *connectionIter, connectionIdx))); - connectionIter++; - connectionIdx++; + const size_t index = segmentLines.size(); + segmentLines.push_back(&shape.second.getLine()); + lineIndexLookup.emplace(static_cast(shape.first), index); } } - time_t atime = 0; + std::vector connectionset(segmentLines.size()); + std::vector> pendingConnections(segmentLines.size()); + if (comm) { - qtimer(atime, 0); - comm->CommPostMessage(Communicator::NUM_RECORDS, lineConnectors.size()); + comm->CommPostMessage(Communicator::NUM_RECORDS, segmentLines.size()); } - double maxdim = __max(m_region.width(), m_region.height()); - - int count = 0; - for (auto &lineConnector_a : lineConnectors) { - Connector &connectionset_a = lineConnector_a.second.m_connector; - const Line &line_a = lineConnector_a.second.m_line; - int idx_a = lineConnector_a.second.m_index; - // n.b., vector() is based on t_start and t_end, so we must use t_start and t_end here and throughout - PixelRef pix1 = pixelate(line_a.t_start()); - std::vector &shapes1 = m_pixel_shapes(static_cast(pix1.y), static_cast(pix1.x)); - for (auto &shape : shapes1) { - auto lineConnector_b = lineConnectors.find(shape.m_shape_ref); - if (lineConnector_b != lineConnectors.end() && idx_a < lineConnector_b->second.m_index) { - - Connector &connectionset_b = lineConnector_b->second.m_connector; - const Line &line_b = lineConnector_b->second.m_line; - int idx_b = lineConnector_b->second.m_index; - - Point2f alpha = line_a.vector(); - Point2f beta = line_b.vector(); - alpha.normalise(); - beta.normalise(); - if (approxeq(line_a.t_start(), line_b.t_start(), (maxdim * TOLERANCE_B))) { - float x = float(2.0 * acos(__min(__max(-dot(alpha, beta), -1.0), 1.0)) / M_PI); - depthmapX::addIfNotExists(connectionset_a.m_back_segconns, SegmentRef(1, idx_b), x); - depthmapX::addIfNotExists(connectionset_b.m_back_segconns, SegmentRef(1, idx_a), x); + const double maxdim = __max(m_region.width(), m_region.height()); + const auto &pixelShapes = static_cast(m_pixel_shapes); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(segmentLines.size()); + + // Endpoint discovery is read-only. Each line records private updates, + // which are committed below in the original line and pixel order. + depthmapX::parallelFor(segmentLines.size(), workerCount, comm, [&](size_t idx_a, size_t) { + const Line &line_a = *segmentLines[idx_a]; + std::vector &updates = pendingConnections[idx_a]; + + auto inspectEndpoint = [&](const Point2f &endpoint, bool aForward) { + const PixelRef pixel = pixelate(endpoint); + const std::vector &nearbyShapes = + pixelShapes(static_cast(pixel.y), static_cast(pixel.x)); + + for (const ShapeRef &shapeRef : nearbyShapes) { + const auto lookup = lineIndexLookup.find(shapeRef.m_shape_ref); + if (lookup == lineIndexLookup.end() || idx_a >= lookup->second) { + continue; } - if (approxeq(line_a.t_start(), line_b.t_end(), (maxdim * TOLERANCE_B))) { - float x = float(2.0 * acos(__min(__max(-dot(alpha, -beta), -1.0), 1.0)) / M_PI); - depthmapX::addIfNotExists(connectionset_a.m_back_segconns, SegmentRef(-1, idx_b), x); - depthmapX::addIfNotExists(connectionset_b.m_forward_segconns, SegmentRef(1, idx_a), x); - } - } - } - - PixelRef pix2 = pixelate(line_a.t_end()); - std::vector &shapes2 = m_pixel_shapes(static_cast(pix2.y), static_cast(pix2.x)); - for (auto &shape : shapes2) { - auto lineConnector_b = lineConnectors.find(shape.m_shape_ref); - if (lineConnector_b != lineConnectors.end() && idx_a < lineConnector_b->second.m_index) { - - Connector &connectionset_b = lineConnector_b->second.m_connector; - const Line &line_b = lineConnector_b->second.m_line; - int idx_b = lineConnector_b->second.m_index; + const size_t idx_b = lookup->second; + const Line &line_b = *segmentLines[idx_b]; Point2f alpha = line_a.vector(); Point2f beta = line_b.vector(); - alpha.normalise(); beta.normalise(); - if (approxeq(line_a.t_end(), line_b.t_start(), (maxdim * TOLERANCE_B))) { - float x = float(2.0 * acos(__min(__max(-dot(-alpha, beta), -1.0), 1.0)) / M_PI); - depthmapX::addIfNotExists(connectionset_a.m_forward_segconns, SegmentRef(1, idx_b), x); - depthmapX::addIfNotExists(connectionset_b.m_back_segconns, SegmentRef(-1, idx_a), x); - } - if (approxeq(line_a.t_end(), line_b.t_end(), (maxdim * TOLERANCE_B))) { - float x = float(2.0 * acos(__min(__max(-dot(-alpha, -beta), -1.0), 1.0)) / M_PI); - depthmapX::addIfNotExists(connectionset_a.m_forward_segconns, SegmentRef(-1, idx_b), x); - depthmapX::addIfNotExists(connectionset_b.m_forward_segconns, SegmentRef(-1, idx_a), x); - } - } - } - if (comm) { - if (qtimer(atime, 500)) { - if (comm->IsCancelled()) { - throw Communicator::CancelledException(); + if (!aForward) { + if (approxeq(line_a.t_start(), line_b.t_start(), maxdim * TOLERANCE_B)) { + const float angle = + float(2.0 * acos(__min(__max(-dot(alpha, beta), -1.0), 1.0)) / M_PI); + updates.push_back({idx_a, idx_b, false, false, 1, 1, angle}); + } + if (approxeq(line_a.t_start(), line_b.t_end(), maxdim * TOLERANCE_B)) { + const float angle = + float(2.0 * acos(__min(__max(-dot(alpha, -beta), -1.0), 1.0)) / M_PI); + updates.push_back({idx_a, idx_b, false, true, -1, 1, angle}); + } + } else { + if (approxeq(line_a.t_end(), line_b.t_start(), maxdim * TOLERANCE_B)) { + const float angle = + float(2.0 * acos(__min(__max(-dot(-alpha, beta), -1.0), 1.0)) / M_PI); + updates.push_back({idx_a, idx_b, true, false, 1, -1, angle}); + } + if (approxeq(line_a.t_end(), line_b.t_end(), maxdim * TOLERANCE_B)) { + const float angle = + float(2.0 * acos(__min(__max(-dot(-alpha, -beta), -1.0), 1.0)) / M_PI); + updates.push_back({idx_a, idx_b, true, true, -1, -1, angle}); + } } - comm->CommPostMessage(Communicator::CURRENT_RECORD, count); } + }; + + inspectEndpoint(line_a.t_start(), false); + inspectEndpoint(line_a.t_end(), true); + }); + + for (size_t idx_a = 0; idx_a < pendingConnections.size(); ++idx_a) { + for (const EndpointConnection &update : pendingConnections[idx_a]) { + auto &connectionsA = update.m_aForward + ? connectionset[update.m_a].m_forward_segconns + : connectionset[update.m_a].m_back_segconns; + auto &connectionsB = update.m_bForward + ? connectionset[update.m_b].m_forward_segconns + : connectionset[update.m_b].m_back_segconns; + depthmapX::addIfNotExists(connectionsA, + SegmentRef(update.m_aDirection, static_cast(update.m_b)), + update.m_angle); + depthmapX::addIfNotExists(connectionsB, + SegmentRef(update.m_bDirection, static_cast(update.m_a)), + update.m_angle); } - count++; } - // initialise attributes now separated from making the connections makeSegmentConnections(connectionset); } @@ -565,39 +586,77 @@ void ShapeGraph::makeNewSegMap(Communicator *comm) { void ShapeGraph::makeSegmentMap(std::vector& lines, std::vector& connectors, double stubremoval) { - // the first (key) pair is the line / line intersection, second is the pair of associated segments for the first line + struct SegmentBreakData { + int m_axis = XAXIS; + int m_parity = 1; + std::vector> m_breaks; + }; + + // The first pair identifies a line intersection. The second pair stores + // the adjacent segments for the line that was processed first. std::map> segmentlist; - // this code relies on the polygon order being the same as the connections + // Shape order and connector order are positional throughout this routine. + const size_t axialCount = m_connectors.size(); + std::vector axialShapes; + std::vector axialRefs; + axialShapes.reserve(axialCount); + axialRefs.reserve(axialCount); + auto shapeIter = m_shapes.begin(); + for (size_t i = 0; i < axialCount; ++i) { + if (shapeIter == m_shapes.end()) { + throw depthmapX::RuntimeException("Axial shapes and connectors are inconsistent"); + } + axialRefs.push_back(shapeIter->first); + axialShapes.push_back(&shapeIter->second); + ++shapeIter; + } + + std::vector breakData(axialCount); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(axialCount); - auto iter = m_shapes.begin(); - for (size_t i = 0; i < m_connectors.size(); i++) { - auto shape = iter->second; - int axialRef = iter->first; - iter++; + // Intersection calculation and sorting are independent for each axial line. + // The actual segment numbering and graph mutation remain serial below. + depthmapX::parallelFor(axialCount, workerCount, nullptr, [&](size_t i, size_t) { + const SalaShape &shape = *axialShapes[i]; if (!shape.isLine()) { - continue; + return; } - const Line& line = shape.getLine(); - std::vector > breaks; // this is a vector instead of a map because the - // original code allowed for duplicate keys - int axis = line.width() >= line.height() ? XAXIS : YAXIS; - // we need the breaks ordered from start to end of the line - // this is automatic for XAXIS, but on YAXIS, need to know - // if the line is ascending or decending - int parity = (axis == XAXIS) ? 1 : line.sign(); - - auto& connections = m_connectors[i].m_connections; - for (size_t j = 0; j < connections.size(); j++) { - // find the intersection point and add... - // note: more than one break at the same place allowed - auto shapeJ = depthmapX::getMapAtIndex(m_shapes, connections[j])->second; - if (static_cast(i) != connections[j] && shapeJ.isLine()) { - breaks.push_back(std::make_pair(parity * line.intersection_point( shapeJ.getLine(), axis, TOLERANCE_A ), - connections[j])); + + const Line &line = shape.getLine(); + SegmentBreakData &data = breakData[i]; + data.m_axis = line.width() >= line.height() ? XAXIS : YAXIS; + data.m_parity = data.m_axis == XAXIS ? 1 : line.sign(); + const std::vector &connections = m_connectors[i].m_connections; + data.m_breaks.reserve(connections.size()); + + for (int connection : connections) { + if (connection < 0 || static_cast(connection) >= axialShapes.size() || + static_cast(i) == connection) { + continue; } + const SalaShape &connectedShape = *axialShapes[static_cast(connection)]; + if (connectedShape.isLine()) { + data.m_breaks.emplace_back( + data.m_parity * line.intersection_point(connectedShape.getLine(), data.m_axis, TOLERANCE_A), + connection); + } + } + std::sort(data.m_breaks.begin(), data.m_breaks.end()); + }); + + const double effectiveStubRemoval = stubremoval == 0.0 ? TOLERANCE_C : stubremoval; + + for (size_t i = 0; i < axialCount; ++i) { + const SalaShape &shape = *axialShapes[i]; + const int axialRef = axialRefs[i]; + if (!shape.isLine()) { + continue; } - std::sort(breaks.begin(), breaks.end()); + const Line &line = shape.getLine(); + const std::vector> &breaks = breakData[i].m_breaks; + const int axis = breakData[i].m_axis; + const int parity = breakData[i].m_parity; // okay, now we have a list from one end of the other of lines this line connects with Point2f lastpoint = line.start(); int seg_a = -1, seg_b = -1; @@ -606,11 +665,8 @@ void ShapeGraph::makeSegmentMap(std::vector& lines, std::vector // found in July 2006. It has been set "high" deliberately (1e-6 = a millionth of the line height / width) // in order to catch small errors made by operators or floating point errors in other systems // when drawing, for example, three axial lines intersecting - if (stubremoval == 0.0) { - // if 0, convert to tolerance - stubremoval = TOLERANCE_C; - } - neardist = (axis == XAXIS) ? (line.width() * stubremoval) : (line.height() * stubremoval); + neardist = (axis == XAXIS) ? (line.width() * effectiveStubRemoval) + : (line.height() * effectiveStubRemoval); double overlapdist = (axis == XAXIS) ? (line.width() * TOLERANCE_C) : (line.height() * TOLERANCE_C); // for (auto breaksIter = breaks.begin(); breaksIter != breaks.end();) { diff --git a/salalib/axialmodules/axialintegration.cpp b/salalib/axialmodules/axialintegration.cpp index c647bfdb..e31553ff 100644 --- a/salalib/axialmodules/axialintegration.cpp +++ b/salalib/axialmodules/axialintegration.cpp @@ -18,19 +18,22 @@ #include "salalib/axialmodules/axialintegration.h" +#include "genlib/parallel.h" #include "genlib/pflipper.h" #include "genlib/stringutils.h" +#include #include +#include +#include +#include bool AxialIntegration::run(Communicator *comm, ShapeGraph &map, bool simple_version) { // note, from 10.0, Depthmap no longer includes *self* connections on axial lines // self connections are stripped out on loading graph files, as well as no longer made - time_t atime = 0; if (comm) { - qtimer(atime, 0); - comm->CommPostMessage(Communicator::NUM_RECORDS, map.getShapeCount()); + comm->CommPostMessage(Communicator::NUM_RECORDS, static_cast(map.getShapeCount())); } AttributeTable &attributes = map.getAttributeTable(); @@ -219,300 +222,323 @@ bool AxialIntegration::run(Communicator *comm, ShapeGraph &map, bool simple_vers } } - // for choice - AnalysisInfo **audittrail; - if (m_choice) { - audittrail = new AnalysisInfo *[map.getShapeCount()]; - for (size_t i = 0; i < map.getShapeCount(); i++) { - audittrail[i] = new AnalysisInfo[radii.size()]; - } + const size_t shapeCount = map.getShapeCount(); + const size_t radiusCount = radii.size(); + + struct RadiusResult { + float nodeCount = 0.0f; + float totalWeight = -1.0f; + float meanDepth = -1.0f; + float weightedMeanDepth = -1.0f; + float integrationDv = -1.0f; + float integrationPv = -1.0f; + float integrationTekl = -1.0f; + float ra = -1.0f; + float rra = -1.0f; + float totalDepth = -1.0f; + float penn = -1.0f; + float entropy = -1.0f; + float relativeEntropy = -1.0f; + float intensity = -1.0f; + float harmonic = -1.0f; + }; + struct ChoiceValue { + double choice = 0.0; + double weightedChoice = 0.0; + }; + struct WorkerState { + std::vector covered; + std::vector previous; + std::vector choice; + + WorkerState(size_t count, size_t radiiCount, bool calculateChoice) + : covered(count), previous(count), choice(calculateChoice ? count * radiiCount : 0) {} + }; + + std::vector results(shapeCount * radiusCount); + std::vector controlResults(shapeCount, -1.0f); + std::vector controllabilityResults(shapeCount, -1.0f); + + const size_t estimatedBytesPerWorker = + shapeCount * (sizeof(std::uint8_t) + sizeof(int) + (m_choice ? radiusCount * sizeof(ChoiceValue) : 0)); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(shapeCount, estimatedBytesPerWorker); + std::vector workerStates; + workerStates.reserve(workerCount); + for (size_t worker = 0; worker < workerCount; ++worker) { + workerStates.emplace_back(shapeCount, radiusCount, m_choice); } - // n.b., for this operation we assume continuous line referencing from zero (this is silly?) - // has already failed due to this! when intro hand drawn fewest line (where user may have deleted) - // it's going to get worse... - - bool *covered = new bool[map.getShapeCount()]; - - size_t i = -1; - for (auto & iter : attributes) { - i++; - AttributeRow &row = iter.getRow(); - for (size_t j = 0; j < map.getShapeCount(); j++) { - covered[j] = false; - } + depthmapX::parallelFor(shapeCount, workerCount, comm, [&](size_t root, size_t workerIndex) { + WorkerState &state = workerStates[workerIndex]; + std::fill(state.covered.begin(), state.covered.end(), std::uint8_t(0)); if (m_choice) { - for (size_t k = 0; k < map.getShapeCount(); k++) { - audittrail[k][0].previous.ref = -1; // note, 0th member used as radius doesn't matter - // note, choice columns are not cleared, but cummulative over all shortest path pairs - } + std::fill(state.previous.begin(), state.previous.end(), -1); } if (m_local) { double control = 0.0; - const std::vector &connections = map.getConnections()[i].m_connections; + const std::vector &connections = map.getConnections()[root].m_connections; std::vector totalneighbourhood; for (int connection : connections) { - // n.b., as of Depthmap 10.0, connections[j] and i cannot coexist - // if (connections[j] != i) { + if (connection < 0 || static_cast(connection) >= shapeCount) { + continue; + } depthmapX::addIfNotExists(totalneighbourhood, connection); - int retro_size = 0; - auto &retconnectors = map.getConnections()[size_t(connection)].m_connections; - for (auto retconnector : retconnectors) { - retro_size++; + const auto &retconnectors = map.getConnections()[size_t(connection)].m_connections; + for (int retconnector : retconnectors) { depthmapX::addIfNotExists(totalneighbourhood, retconnector); } - control += 1.0 / double(retro_size); - //} + if (!retconnectors.empty()) { + control += 1.0 / double(retconnectors.size()); + } } - - if (!simple_version) { - if (connections.size() > 0) { - row.setValue(control_col, float(control)); - row.setValue(controllability_col, - float(double(connections.size()) / double(totalneighbourhood.size() - 1))); - } else { - row.setValue(control_col, -1); - row.setValue(controllability_col, -1); + if (!simple_version && !connections.empty()) { + controlResults[root] = float(control); + if (totalneighbourhood.size() > 1) { + controllabilityResults[root] = + float(double(connections.size()) / double(totalneighbourhood.size() - 1)); } } } - std::vector depthcounts; - depthcounts.push_back(0); - + std::vector depthcounts(1, 0); pflipper>> foundlist; - foundlist.a().push_back(std::pair(i, -1)); - covered[i] = true; - int total_depth = 0, depth = 1, node_count = 1, pos = -1, previous = -1; // node_count includes this 1 - double weight = 0.0, rootweight = 0.0, total_weight = 0.0, w_total_depth = 0.0; + foundlist.a().push_back(std::make_pair(static_cast(root), -1)); + state.covered[root] = 1; + + int totalDepth = 0; + int depth = 1; + int nodeCount = 1; + double rootweight = 0.0; + double totalWeight = 0.0; + double weightedTotalDepth = 0.0; if (m_weighted_measure_col != -1) { - rootweight = weights[i]; - // include this line in total weights (as per nodecount) - total_weight += rootweight; + rootweight = weights[root]; + totalWeight = rootweight; } - int index = -1; - int r = 0; - for (int radius : radii) { - while (foundlist.a().size()) { - if (!m_choice) { - index = foundlist.a().back().first; - } else { - pos = pafrand() % foundlist.a().size(); - index = foundlist.a().at(pos).first; - previous = foundlist.a().at(pos).second; - audittrail[index][0].previous.ref = - previous; // note 0th member used here: can be used individually different radius previous + + std::minstd_rand random(static_cast(0x9e3779b9U ^ static_cast(root))); + size_t processed = 0; + + for (size_t radiusIndex = 0; radiusIndex < radiusCount; ++radiusIndex) { + const int radius = radii[radiusIndex]; + while (!foundlist.a().empty()) { + if (comm && ((processed++ & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); } - Connector &line = map.getConnections()[index]; - for (size_t k = 0; k < line.m_connections.size(); k++) { - if (!covered[line.m_connections[k]]) { - covered[line.m_connections[k]] = true; - foundlist.b().push_back(std::pair(line.m_connections[k], index)); - if (m_weighted_measure_col != -1) { - // the weight is taken from the discovered node: - weight = weights[line.m_connections[k]]; - total_weight += weight; - w_total_depth += depth * weight; + + size_t position = foundlist.a().size() - 1; + if (m_choice) { + position = static_cast(random()) % foundlist.a().size(); + } + const int index = foundlist.a()[position].first; + const int previous = foundlist.a()[position].second; + if (m_choice) { + state.previous[index] = previous; + } + + const Connector &line = map.getConnections()[size_t(index)]; + for (int connection : line.m_connections) { + if (connection < 0 || static_cast(connection) >= shapeCount || state.covered[connection]) { + continue; + } + state.covered[connection] = 1; + foundlist.b().push_back(std::make_pair(connection, index)); + + double weight = 0.0; + if (m_weighted_measure_col != -1) { + weight = weights[connection]; + totalWeight += weight; + weightedTotalDepth += depth * weight; + } + + if (m_choice && previous != -1) { + int here = index; + while (here != static_cast(root) && here >= 0 && static_cast(here) < shapeCount) { + ChoiceValue &choice = state.choice[size_t(here) * radiusCount + radiusIndex]; + choice.choice += 1.0; + choice.weightedChoice += weight * rootweight; + here = state.previous[here]; } - if (m_choice && previous != -1) { - // both directional paths are now recorded for choice - // (coincidentally fixes choice problem which was completely wrong) - size_t here = index; // note: start counting from index as actually looking ahead here - while (here != i) { // not i means not the current root for the path - audittrail[here][r].choice += 1; - audittrail[here][r].weighted_choice += weight * rootweight; - here = - audittrail[here][0].previous.ref; // <- note, just using 0th position: radius for - // the previous doesn't matter in this analysis - } - if (m_weighted_measure_col != -1) { - // in weighted choice, root node and current node receive values: - audittrail[i][r].weighted_choice += (weight * rootweight) * 0.5; - audittrail[line.m_connections[k]][r].weighted_choice += (weight * rootweight) * 0.5; - } + if (m_weighted_measure_col != -1) { + state.choice[root * radiusCount + radiusIndex].weightedChoice += + (weight * rootweight) * 0.5; + state.choice[size_t(connection) * radiusCount + radiusIndex].weightedChoice += + (weight * rootweight) * 0.5; } - total_depth += depth; - node_count++; - depthcounts.back() += 1; } + totalDepth += depth; + ++nodeCount; + ++depthcounts.back(); } - if (!m_choice) + + if (m_choice) { + foundlist.a().erase(foundlist.a().begin() + static_cast(position)); + } else { foundlist.a().pop_back(); - else - foundlist.a().erase(foundlist.a().begin() + pos); - if (!foundlist.a().size()) { + } + if (foundlist.a().empty()) { foundlist.flip(); - depth++; + ++depth; depthcounts.push_back(0); if (radius != -1 && depth > radius) { break; } } } - // set the attributes for this node: - row.setValue(count_col[r], float(node_count)); + + RadiusResult &result = results[root * radiusCount + radiusIndex]; + result.nodeCount = float(nodeCount); if (m_weighted_measure_col != -1) { - row.setValue(total_weight_col[r], float(total_weight)); + result.totalWeight = float(totalWeight); } - // node count > 1 to avoid divide by zero (was > 2) - if (node_count > 1) { - // note -- node_count includes this one -- mean depth as per p.108 Social Logic of Space - double mean_depth = double(total_depth) / double(node_count - 1); - row.setValue(depth_col[r], float(mean_depth)); - if (m_weighted_measure_col != -1) { - // weighted mean depth: - row.setValue(w_depth_col[r], float(w_total_depth / total_weight)); + + if (nodeCount > 1) { + const double meanDepth = double(totalDepth) / double(nodeCount - 1); + result.meanDepth = float(meanDepth); + if (m_weighted_measure_col != -1 && totalWeight != 0.0) { + result.weightedMeanDepth = float(weightedTotalDepth / totalWeight); } - // total nodes > 2 to avoid divide by 0 (was > 3) - if (node_count > 2 && mean_depth > 1.0) { - double ra = 2.0 * (mean_depth - 1.0) / double(node_count - 2); - // d-value / p-value from Depthmap 4 manual, note: node_count includes this one - double rra_d = ra / dvalue(node_count); - double rra_p = ra / dvalue(node_count); - double integ_tk = teklinteg(node_count, total_depth); - row.setValue(integ_dv_col[r], float(1.0 / rra_d)); + if (nodeCount > 2 && meanDepth > 1.0) { + const double ra = 2.0 * (meanDepth - 1.0) / double(nodeCount - 2); + const double rraD = ra / dvalue(nodeCount); + const double rraP = ra / dvalue(nodeCount); + result.integrationDv = float(1.0 / rraD); if (!simple_version) { - row.setValue(integ_pv_col[r], float(1.0 / rra_p)); - if (total_depth - node_count + 1 > 1) { - row.setValue(integ_tk_col[r], float(integ_tk)); - } else { - row.setValue(integ_tk_col[r], -1.0f); - } + result.integrationPv = float(1.0 / rraP); + result.integrationTekl = + totalDepth - nodeCount + 1 > 1 ? float(teklinteg(nodeCount, totalDepth)) : -1.0f; } - if (m_fulloutput) { - row.setValue(ra_col[r], float(ra)); - - if (!simple_version) { - row.setValue(rra_col[r], float(rra_d)); - } - row.setValue(td_col[r], float(total_depth)); - + result.ra = float(ra); + result.totalDepth = float(totalDepth); if (!simple_version) { - // alan's palm-tree normalisation: palmtree - double dmin = node_count - 1; - double dmax = palmtree(node_count, depth - 1); + result.rra = float(rraD); + const double dmin = nodeCount - 1; + const double dmax = palmtree(nodeCount, depth - 1); if (dmax != dmin) { - row.setValue(penn_norm_col[r], float((dmax - total_depth) / (dmax - dmin))); + result.penn = float((dmax - totalDepth) / (dmax - dmin)); } } } - } else { - row.setValue(integ_dv_col[r], -1.0f); - - if (!simple_version) { - row.setValue(integ_pv_col[r], -1.0f); - row.setValue(integ_tk_col[r], -1.0f); - } - if (m_fulloutput) { - row.setValue(ra_col[r], -1.0f); - - if (!simple_version) { - row.setValue(rra_col[r], -1.0f); - } - - row.setValue(td_col[r], -1.0f); - - if (!simple_version) { - row.setValue(penn_norm_col[r], -1.0f); - } - } } if (!simple_version) { - double entropy = 0.0, intensity = 0.0, rel_entropy = 0.0, factorial = 1.0, harmonic = 0.0; - for (size_t k = 0; k < depthcounts.size(); k++) { - if (depthcounts[k] != 0) { - // some debate over whether or not this should be node count - 1 - // (i.e., including or not including the node itself) - double prob = double(depthcounts[k]) / double(node_count); - entropy -= prob * std::log2(prob); - // Formula from Turner 2001, "Depthmap" - factorial *= double(k + 1); - double q = (pow(mean_depth, double(k)) / double(factorial)) * exp(-mean_depth); - rel_entropy += (double)prob * std::log2(prob / q); - // - harmonic += 1.0 / double(depthcounts[k]); + double entropy = 0.0; + double intensity = 0.0; + double relativeEntropy = 0.0; + double factorial = 1.0; + double harmonicDenominator = 0.0; + for (size_t k = 0; k < depthcounts.size(); ++k) { + if (depthcounts[k] == 0) { + continue; } + const double probability = double(depthcounts[k]) / double(nodeCount); + entropy -= probability * std::log2(probability); + factorial *= double(k + 1); + const double q = + (std::pow(meanDepth, double(k)) / factorial) * std::exp(-meanDepth); + if (q > 0.0) { + relativeEntropy += probability * std::log2(probability / q); + } + harmonicDenominator += 1.0 / double(depthcounts[k]); } - harmonic = double(depthcounts.size()) / harmonic; - if (total_depth > node_count) { - intensity = node_count * entropy / (total_depth - node_count); - } else { - intensity = -1; + if (harmonicDenominator > 0.0) { + result.harmonic = float(double(depthcounts.size()) / harmonicDenominator); } - row.setValue(entropy_col[r], float(entropy)); - row.setValue(rel_entropy_col[r], float(rel_entropy)); - row.setValue(intensity_col[r], float(intensity)); - row.setValue(harmonic_col[r], float(harmonic)); - } - } else { - row.setValue(depth_col[r], -1.0f); - row.setValue(integ_dv_col[r], -1.0f); - - if (!simple_version) { - row.setValue(integ_pv_col[r], -1.0f); - row.setValue(integ_tk_col[r], -1.0f); - row.setValue(entropy_col[r], -1.0f); - row.setValue(rel_entropy_col[r], -1.0f); - row.setValue(harmonic_col[r], -1.0f); + intensity = totalDepth > nodeCount ? nodeCount * entropy / (totalDepth - nodeCount) : -1.0; + result.entropy = float(entropy); + result.relativeEntropy = float(relativeEntropy); + result.intensity = float(intensity); } } - ++r; } - // - if (comm) { - if (qtimer(atime, 500)) { - if (comm->IsCancelled()) { - delete[] covered; - throw Communicator::CancelledException(); + }); + + size_t root = 0; + for (auto &attributeIter : attributes) { + if (comm && ((root & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + AttributeRow &row = attributeIter.getRow(); + if (m_local && !simple_version) { + row.setValue(control_col, controlResults[root]); + row.setValue(controllability_col, controllabilityResults[root]); + } + for (size_t radiusIndex = 0; radiusIndex < radiusCount; ++radiusIndex) { + const RadiusResult &result = results[root * radiusCount + radiusIndex]; + row.setValue(count_col[radiusIndex], result.nodeCount); + row.setValue(depth_col[radiusIndex], result.meanDepth); + row.setValue(integ_dv_col[radiusIndex], result.integrationDv); + if (m_weighted_measure_col != -1) { + row.setValue(total_weight_col[radiusIndex], result.totalWeight); + row.setValue(w_depth_col[radiusIndex], result.weightedMeanDepth); + } + if (!simple_version) { + row.setValue(integ_pv_col[radiusIndex], result.integrationPv); + row.setValue(integ_tk_col[radiusIndex], result.integrationTekl); + row.setValue(entropy_col[radiusIndex], result.entropy); + row.setValue(rel_entropy_col[radiusIndex], result.relativeEntropy); + row.setValue(intensity_col[radiusIndex], result.intensity); + row.setValue(harmonic_col[radiusIndex], result.harmonic); + } + if (m_fulloutput) { + row.setValue(ra_col[radiusIndex], result.ra); + row.setValue(td_col[radiusIndex], result.totalDepth); + if (!simple_version) { + row.setValue(rra_col[radiusIndex], result.rra); + row.setValue(penn_norm_col[radiusIndex], result.penn); } - comm->CommPostMessage(Communicator::CURRENT_RECORD, i); } } + ++root; } - delete[] covered; + if (m_choice) { - i = -1; - for (auto & iter: attributes) { - i++; - AttributeRow &row = iter.getRow(); - double total_choice = 0.0, w_total_choice = 0.0; - for (size_t r = 0; r < radii.size(); r++) { - total_choice += audittrail[i][r].choice; - w_total_choice += audittrail[i][r].weighted_choice; - // n.b., normalise choice according to (n-1)(n-2)/2 (maximum possible through routes) - double node_count = row.getValue(count_col[r]); - double total_weight = 0; - if (m_weighted_measure_col != -1) { - total_weight = row.getValue(total_weight_col[r]); + root = 0; + for (auto &attributeIter : attributes) { + if (comm && ((root & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + AttributeRow &row = attributeIter.getRow(); + double cumulativeChoice = 0.0; + double cumulativeWeightedChoice = 0.0; + for (size_t radiusIndex = 0; radiusIndex < radiusCount; ++radiusIndex) { + for (size_t worker = 0; worker < workerCount; ++worker) { + const ChoiceValue &choice = + workerStates[worker].choice[root * radiusCount + radiusIndex]; + cumulativeChoice += choice.choice; + cumulativeWeightedChoice += choice.weightedChoice; } - if (node_count > 2) { - row.setValue(choice_col[r], float(total_choice)); - row.setValue(n_choice_col[r], float(2.0 * total_choice / ((node_count - 1) * (node_count - 2)))); + const RadiusResult &result = results[root * radiusCount + radiusIndex]; + const double nodeCount = result.nodeCount; + const double totalWeight = result.totalWeight; + if (nodeCount > 2.0) { + row.setValue(choice_col[radiusIndex], float(cumulativeChoice)); + row.setValue(n_choice_col[radiusIndex], + float(2.0 * cumulativeChoice / ((nodeCount - 1.0) * (nodeCount - 2.0)))); if (m_weighted_measure_col != -1) { - row.setValue(w_choice_col[r], float(w_total_choice)); - row.setValue(nw_choice_col[r], float(2.0 * w_total_choice / (total_weight * total_weight))); + row.setValue(w_choice_col[radiusIndex], float(cumulativeWeightedChoice)); + row.setValue(nw_choice_col[radiusIndex], + totalWeight != 0.0 + ? float(2.0 * cumulativeWeightedChoice / (totalWeight * totalWeight)) + : -1.0f); } } else { - row.setValue(choice_col[r], -1); - row.setValue(n_choice_col[r], -1); + row.setValue(choice_col[radiusIndex], -1.0f); + row.setValue(n_choice_col[radiusIndex], -1.0f); if (m_weighted_measure_col != -1) { - row.setValue(w_choice_col[r], -1); - row.setValue(nw_choice_col[r], -1); + row.setValue(w_choice_col[radiusIndex], -1.0f); + row.setValue(nw_choice_col[radiusIndex], -1.0f); } } } + ++root; } - for (size_t i = 0; i < map.getShapeCount(); i++) { - delete[] audittrail[i]; - } - delete[] audittrail; } - map.setDisplayedAttribute(-1); // <- override if it's already showing + map.setDisplayedAttribute(-1); map.setDisplayedAttribute(integ_dv_col.back()); return true; diff --git a/salalib/segmmodules/segmangular.cpp b/salalib/segmmodules/segmangular.cpp index 93d34e5e..6b30cf76 100644 --- a/salalib/segmmodules/segmangular.cpp +++ b/salalib/segmmodules/segmangular.cpp @@ -19,8 +19,12 @@ #include "salalib/segmmodules/segmangular.h" #include "salalib/options.h" +#include "genlib/parallel.h" #include "genlib/stringutils.h" +#include +#include + bool SegmentAngular::run(Communicator *comm, ShapeGraph &map, bool) { if (map.getMapType() != ShapeMap::SEGMENTMAP) { @@ -28,12 +32,7 @@ bool SegmentAngular::run(Communicator *comm, ShapeGraph &map, bool) { } AttributeTable &attributes = map.getAttributeTable(); - - time_t atime = 0; - if (comm) { - qtimer(atime, 0); - comm->CommPostMessage(Communicator::NUM_RECORDS, map.getConnections().size()); - } + const size_t shapeCount = map.getShapeCount(); // note: radius must be sorted lowest to highest, but if -1 occurs ("radius n") it needs to be last... // ...to ensure no mess ups, we'll re-sort here: @@ -50,121 +49,128 @@ bool SegmentAngular::run(Communicator *comm, ShapeGraph &map, bool) { radii.push_back(-1.0); } - std::vector depth_col, count_col, total_col; - // first enter table values - for (int radius : radii) { + std::vector depth_col, count_col, total_col, integration_col; + for (double radius : radii) { std::string radius_text = makeRadiusText(Options::RADIUS_ANGULAR, radius); - std::string depth_col_text = std::string("Angular Mean Depth") + radius_text; - attributes.insertOrResetColumn(depth_col_text.c_str()); - std::string count_col_text = std::string("Angular Node Count") + radius_text; - attributes.insertOrResetColumn(count_col_text.c_str()); - std::string total_col_text = std::string("Angular Total Depth") + radius_text; - attributes.insertOrResetColumn(total_col_text.c_str()); + depth_col.push_back(attributes.insertOrResetColumn(std::string("Angular Mean Depth") + radius_text)); + count_col.push_back(attributes.insertOrResetColumn(std::string("Angular Node Count") + radius_text)); + total_col.push_back(attributes.insertOrResetColumn(std::string("Angular Total Depth") + radius_text)); + integration_col.push_back(attributes.insertOrResetColumn(std::string("Angular Integration") + radius_text)); } - for (int radius : radii) { - std::string radius_text = makeRadiusText(Options::RADIUS_ANGULAR, radius); - std::string depth_col_text = std::string("Angular Mean Depth") + radius_text; - depth_col.push_back(attributes.getColumnIndex(depth_col_text.c_str())); - std::string count_col_text = std::string("Angular Node Count") + radius_text; - count_col.push_back(attributes.getColumnIndex(count_col_text.c_str())); - std::string total_col_text = std::string("Angular Total Depth") + radius_text; - total_col.push_back(attributes.getColumnIndex(total_col_text.c_str())); + if (comm) { + comm->CommPostMessage(Communicator::NUM_RECORDS, static_cast(shapeCount)); } - std::vector covered(map.getShapeCount()); - size_t i = 0; - for (auto & iter : attributes){ - for (size_t j = 0; j < map.getShapeCount(); j++) { - covered[j] = false; - } + const size_t radiusCount = radii.size(); + std::vector resultDepth(shapeCount * radiusCount, -1.0); + std::vector resultCount(shapeCount * radiusCount, 0); + + const size_t bytesPerWorker = shapeCount * sizeof(std::uint8_t); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(shapeCount, bytesPerWorker); + std::vector> covered(workerCount, std::vector(shapeCount)); + + depthmapX::parallelFor(shapeCount, workerCount, comm, [&](size_t root, size_t worker) { + std::vector &workerCovered = covered[worker]; + std::fill(workerCovered.begin(), workerCovered.end(), std::uint8_t(0)); + std::vector> anglebins; - anglebins.push_back(std::make_pair(0.0f, SegmentData(0, i, SegmentRef(), 0, 0.0, 0))); + anglebins.push_back(std::make_pair(0.0f, SegmentData(0, root, SegmentRef(), 0, 0.0, 0))); - std::vector total_depth; - std::vector node_count; - for (size_t r = 0; r < radii.size(); r++) { - total_depth.push_back(0.0); - node_count.push_back(0); - } - // node_count includes this one, but will be added in next algo: - while (anglebins.size()) { - auto iter = anglebins.begin(); - SegmentData lineindex = iter->second; - if (!covered[lineindex.ref]) { - covered[lineindex.ref] = true; - double depth_to_line = iter->first; - total_depth[lineindex.coverage] += depth_to_line; - node_count[lineindex.coverage] += 1; - anglebins.erase(iter); - Connector &line = map.getConnections()[lineindex.ref]; + std::vector totalDepth(radiusCount, 0.0); + std::vector nodeCount(radiusCount, 0); + size_t processed = 0; + + while (!anglebins.empty()) { + if (comm && ((processed++ & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + + auto angleIter = anglebins.begin(); + SegmentData lineindex = angleIter->second; + if (!workerCovered[lineindex.ref]) { + workerCovered[lineindex.ref] = 1; + const double depthToLine = angleIter->first; + totalDepth[lineindex.coverage] += depthToLine; + nodeCount[lineindex.coverage] += 1; + anglebins.erase(angleIter); + + const Connector &line = map.getConnections()[lineindex.ref]; if (lineindex.dir != -1) { - for (auto &segconn : line.m_forward_segconns) { - if (!covered[segconn.first.ref]) { - double angle = depth_to_line + segconn.second; + for (const auto &segconn : line.m_forward_segconns) { + if (!workerCovered[segconn.first.ref]) { + const double angle = depthToLine + segconn.second; size_t rbin = lineindex.coverage; - while (rbin != radii.size() && radii[rbin] != -1 && angle > radii[rbin]) { - rbin++; + while (rbin != radiusCount && radii[rbin] != -1 && angle > radii[rbin]) { + ++rbin; } - if (rbin != radii.size()) { + if (rbin != radiusCount) { depthmapX::insert_sorted( - anglebins, std::make_pair(float(angle), - SegmentData(segconn.first, SegmentRef(), 0, 0.0, rbin))); + anglebins, + std::make_pair(float(angle), SegmentData(segconn.first, SegmentRef(), 0, 0.0, rbin))); } } } } if (lineindex.dir != 1) { - for (auto &segconn : line.m_back_segconns) { - if (!covered[segconn.first.ref]) { - double angle = depth_to_line + segconn.second; + for (const auto &segconn : line.m_back_segconns) { + if (!workerCovered[segconn.first.ref]) { + const double angle = depthToLine + segconn.second; size_t rbin = lineindex.coverage; - while (rbin != radii.size() && radii[rbin] != -1 && angle > radii[rbin]) { - rbin++; + while (rbin != radiusCount && radii[rbin] != -1 && angle > radii[rbin]) { + ++rbin; } - if (rbin != radii.size()) { + if (rbin != radiusCount) { depthmapX::insert_sorted( - anglebins, std::make_pair(float(angle), - SegmentData(segconn.first, SegmentRef(), 0, 0.0, rbin))); + anglebins, + std::make_pair(float(angle), SegmentData(segconn.first, SegmentRef(), 0, 0.0, rbin))); } } } } } else { - anglebins.erase(iter); + anglebins.erase(angleIter); } } - AttributeRow &row = iter.getRow(); - // set the attributes for this node: - int curs_node_count = 0; - double curs_total_depth = 0.0; - for (size_t r = 0; r < radii.size(); r++) { - curs_node_count += node_count[r]; - curs_total_depth += total_depth[r]; - row.setValue(count_col[r], float(curs_node_count)); - if (curs_node_count > 1) { - // note -- node_count includes this one -- mean depth as per p.108 Social Logic of Space - double mean_depth = curs_total_depth / double(curs_node_count - 1); - row.setValue(depth_col[r], float(mean_depth)); - row.setValue(total_col[r], float(curs_total_depth)); - } else { - row.setValue(depth_col[r], -1); - row.setValue(total_col[r], -1); - } + + int cumulativeCount = 0; + double cumulativeDepth = 0.0; + for (size_t radiusIndex = 0; radiusIndex < radiusCount; ++radiusIndex) { + cumulativeCount += nodeCount[radiusIndex]; + cumulativeDepth += totalDepth[radiusIndex]; + const size_t resultIndex = root * radiusCount + radiusIndex; + resultCount[resultIndex] = cumulativeCount; + resultDepth[resultIndex] = cumulativeDepth; } - // - if (comm) { - if (qtimer(atime, 500)) { - if (comm->IsCancelled()) { - throw Communicator::CancelledException(); - } - comm->CommPostMessage(Communicator::CURRENT_RECORD, i); + }); + + size_t root = 0; + for (auto &attributeIter : attributes) { + if (comm && ((root & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + AttributeRow &row = attributeIter.getRow(); + for (size_t radiusIndex = 0; radiusIndex < radiusCount; ++radiusIndex) { + const size_t resultIndex = root * radiusCount + radiusIndex; + const int nodeCount = resultCount[resultIndex]; + const double totalDepth = resultDepth[resultIndex]; + row.setValue(count_col[radiusIndex], float(nodeCount)); + if (nodeCount > 1) { + row.setValue(depth_col[radiusIndex], float(totalDepth / double(nodeCount - 1))); + row.setValue(total_col[radiusIndex], float(totalDepth)); + row.setValue(integration_col[radiusIndex], + totalDepth > 1e-9 ? float(double(nodeCount) * double(nodeCount) / totalDepth) + : -1.0f); + } else { + row.setValue(depth_col[radiusIndex], -1.0f); + row.setValue(total_col[radiusIndex], -1.0f); + row.setValue(integration_col[radiusIndex], -1.0f); } } - i++; + ++root; } - map.setDisplayedAttribute(-2); // <- override if it's already showing + map.setDisplayedAttribute(-2); map.setDisplayedAttribute(depth_col.back()); return true; diff --git a/salalib/segmmodules/segmmetric.cpp b/salalib/segmmodules/segmmetric.cpp index 64b8ba93..d35be69f 100644 --- a/salalib/segmmodules/segmmetric.cpp +++ b/salalib/segmmodules/segmmetric.cpp @@ -18,188 +18,218 @@ #include "salalib/segmmodules/segmmetric.h" +#include "genlib/parallel.h" #include "genlib/stringutils.h" -bool SegmentMetric::run(Communicator *comm, ShapeGraph &map, bool) { +#include +#include +#include +#include +#include + +namespace { +struct SegmentMetricResult { + double meanDepth = -1.0; + double weightedMeanDepth = -1.0; + double totalDepth = -1.0; + double totalNodes = 0.0; + double totalLength = 0.0; +}; + +struct SegmentMetricWorkerState { + std::vector seen; + std::vector audittrail; + std::array, 512> bins; + std::vector choice; + + SegmentMetricWorkerState(size_t shapeCount, bool calculateChoice) + : seen(shapeCount), audittrail(shapeCount), choice(calculateChoice ? shapeCount : 0) {} +}; +} // namespace +bool SegmentMetric::run(Communicator *comm, ShapeGraph &map, bool) { AttributeTable &attributes = map.getAttributeTable(); + const size_t shapeCount = map.getShapeCount(); - bool retvar = true; - - time_t atime = 0; - - if (comm) { - qtimer(atime, 0); - comm->CommPostMessage(Communicator::NUM_RECORDS, - (m_sel_only ? map.getSelSet().size() : map.getConnections().size())); - } - int reccount = 0; - - // record axial line refs for topological analysis - std::vector axialrefs; - // quick through to find the longest seg length std::vector seglengths; + seglengths.reserve(shapeCount); float maxseglength = 0.0f; - for (size_t cursor = 0; cursor < map.getShapeCount(); cursor++) { - AttributeRow& row = map.getAttributeRowFromShapeIndex(cursor); - axialrefs.push_back(row.getValue(attributes.getColumnIndex("Axial Line Ref"))); - seglengths.push_back(row.getValue(attributes.getColumnIndex("Segment Length"))); - if (seglengths.back() > maxseglength) { - maxseglength = seglengths.back(); + std::vector roots; + roots.reserve(shapeCount); + for (size_t cursor = 0; cursor < shapeCount; ++cursor) { + AttributeRow &row = map.getAttributeRowFromShapeIndex(cursor); + const float length = row.getValue(attributes.getColumnIndex("Segment Length")); + seglengths.push_back(length); + maxseglength = std::max(maxseglength, length); + if (!m_sel_only || row.isSelected()) { + roots.push_back(cursor); } } - std::string prefix, suffix; - int maxbin = 512; - prefix = "Metric "; - + std::string suffix; if (m_radius != -1.0) { suffix = dXstring::formatString(m_radius, " R%.f metric"); } - std::string choicecol = prefix + "Choice" + suffix; - std::string wchoicecol = prefix + "Choice [SLW]" + suffix; - std::string meandepthcol = prefix + "Mean Depth" + suffix; - std::string wmeandepthcol = prefix + std::string("Mean Depth [SLW]") + suffix; - std::string totaldcol = prefix + "Total Depth" + suffix; - std::string totalcol = prefix + "Total Nodes" + suffix; - std::string wtotalcol = prefix + "Total Length" + suffix; - // + const std::string prefix = "Metric "; + const std::string choicecol = prefix + "Choice" + suffix; + const std::string wchoicecol = prefix + "Choice [SLW]" + suffix; + const std::string meandepthcol = prefix + "Mean Depth" + suffix; + const std::string wmeandepthcol = prefix + "Mean Depth [SLW]" + suffix; + const std::string totaldcol = prefix + "Total Depth" + suffix; + const std::string totalcol = prefix + "Total Nodes" + suffix; + const std::string wtotalcol = prefix + "Total Length" + suffix; + if (!m_sel_only) { - attributes.insertOrResetColumn(choicecol.c_str()); - attributes.insertOrResetColumn(wchoicecol.c_str()); + attributes.insertOrResetColumn(choicecol); + attributes.insertOrResetColumn(wchoicecol); } - attributes.insertOrResetColumn(meandepthcol.c_str()); - attributes.insertOrResetColumn(wmeandepthcol.c_str()); - attributes.insertOrResetColumn(totaldcol.c_str()); - attributes.insertOrResetColumn(totalcol.c_str()); - attributes.insertOrResetColumn(wtotalcol.c_str()); - // - std::vector seen(map.getShapeCount()); - std::vector audittrail(map.getShapeCount()); - std::vector choicevals(map.getShapeCount()); - for (size_t cursor = 0; cursor < map.getShapeCount(); cursor++) { - AttributeRow& row = map.getAttributeRowFromShapeIndex(cursor); - if (m_sel_only && !row.isSelected()) { - continue; - } - for (size_t i = 0; i < map.getShapeCount(); i++) { - seen[i] = 0xffffffff; + attributes.insertOrResetColumn(meandepthcol); + attributes.insertOrResetColumn(wmeandepthcol); + attributes.insertOrResetColumn(totaldcol); + attributes.insertOrResetColumn(totalcol); + attributes.insertOrResetColumn(wtotalcol); + + if (comm) { + comm->CommPostMessage(Communicator::NUM_RECORDS, static_cast(roots.size())); + } + + std::vector results(shapeCount); + const size_t estimatedBytesPerWorker = + shapeCount * (sizeof(unsigned int) + sizeof(TopoMetSegmentRef) + (!m_sel_only ? sizeof(TopoMetSegmentChoice) : 0)); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(roots.size(), estimatedBytesPerWorker); + + std::vector workerStates; + workerStates.reserve(workerCount); + for (size_t worker = 0; worker < workerCount; ++worker) { + workerStates.emplace_back(shapeCount, !m_sel_only); + } + + depthmapX::parallelFor(roots.size(), workerCount, comm, [&](size_t task, size_t workerIndex) { + const size_t cursor = roots[task]; + SegmentMetricWorkerState &state = workerStates[workerIndex]; + std::fill(state.seen.begin(), state.seen.end(), std::numeric_limits::max()); + for (auto &bin : state.bins) { + bin.clear(); } - std::vector list[512]; // 512 bins! + int bin = 0; - list[bin].push_back(cursor); - double rootseglength = seglengths[cursor]; - audittrail[cursor] = TopoMetSegmentRef(cursor, Connector::SEG_CONN_ALL, rootseglength * 0.5, -1); + state.bins[bin].push_back(static_cast(cursor)); + const double rootseglength = seglengths[cursor]; + state.audittrail[cursor] = + TopoMetSegmentRef(static_cast(cursor), Connector::SEG_CONN_ALL, rootseglength * 0.5, -1); + int open = 1; unsigned int segdepth = 0; - double total = 0.0, wtotal = 0.0, wtotaldepth = 0.0, totalsegdepth = 0.0, totalmetdepth = 0.0; + double total = 0.0; + double wtotal = 0.0; + double wtotaldepth = 0.0; + double totalmetdepth = 0.0; + size_t processed = 0; + while (open != 0) { - while (list[bin].size() == 0) { - bin++; - segdepth += 1; - if (bin == maxbin) { - bin = 0; - } + if (comm && ((processed++ & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); } - // - TopoMetSegmentRef &here = audittrail[list[bin].back()]; - list[bin].pop_back(); - open--; - // + while (state.bins[bin].empty()) { + bin = (bin + 1) % 512; + ++segdepth; + } + + TopoMetSegmentRef &here = state.audittrail[state.bins[bin].back()]; + state.bins[bin].pop_back(); + --open; if (here.done) { continue; - } else { - here.done = true; } - // - double len = seglengths[here.ref]; - totalsegdepth += segdepth; - totalmetdepth += here.dist - len * 0.5; // preloaded with length ahead + here.done = true; + + const double len = seglengths[here.ref]; + totalmetdepth += here.dist - len * 0.5; wtotal += len; wtotaldepth += len * (here.dist - len * 0.5); - total += 1; - // - Connector &axline = map.getConnections().at(here.ref); - int connected_cursor = -2; - - auto iter = axline.m_back_segconns.begin(); - bool backsegs = true; - - while (connected_cursor != -1) { - if (backsegs && iter == axline.m_back_segconns.end()) { - iter = axline.m_forward_segconns.begin(); - backsegs = false; - } - if (!backsegs && iter == axline.m_forward_segconns.end()) { - break; - } + total += 1.0; - connected_cursor = iter->first.ref; - - if (seen[connected_cursor] > segdepth && static_cast(connected_cursor) != cursor) { - bool seenalready = (seen[connected_cursor] == 0xffffffff) ? false : true; - float length = seglengths[connected_cursor]; - audittrail[connected_cursor] = + const Connector &axline = map.getConnections().at(here.ref); + auto processConnection = [&](int connected_cursor) { + if (connected_cursor < 0 || static_cast(connected_cursor) >= shapeCount) { + return; + } + if (state.seen[connected_cursor] > segdepth && static_cast(connected_cursor) != cursor) { + const bool seenalready = state.seen[connected_cursor] != std::numeric_limits::max(); + const float length = seglengths[connected_cursor]; + state.audittrail[connected_cursor] = TopoMetSegmentRef(connected_cursor, here.dir, here.dist + length, here.ref); - seen[connected_cursor] = segdepth; + state.seen[connected_cursor] = segdepth; if (m_radius == -1 || here.dist + length < m_radius) { - // puts in a suitable bin ahead of us... - open++; - // - // better to divide by 511 but have 512 bins... - list[(bin + int(floor(0.5 + 511 * length / maxseglength))) % 512].push_back(connected_cursor); + ++open; + const int offset = maxseglength > 0.0f + ? int(std::floor(0.5 + 511.0 * length / maxseglength)) + : 0; + state.bins[(bin + offset) % 512].push_back(connected_cursor); } - // not sure why this is outside the radius restriction - // (sel_only: with restricted selection set, not all lines will be labelled) - // (seenalready: need to check that we're not doing this twice, given the seen can go twice) - - // Quick mod - TV - if (!m_sel_only && connected_cursor > int(cursor) && - !seenalready) { // only one way paths, saves doing this twice + if (!m_sel_only && connected_cursor > static_cast(cursor) && !seenalready) { int subcur = connected_cursor; while (subcur != -1) { - // in this method of choice, start and end lines are included - choicevals[subcur].choice += 1; - choicevals[subcur].wchoice += (rootseglength * length); - subcur = audittrail[subcur].previous; + state.choice[subcur].choice += 1.0; + state.choice[subcur].wchoice += rootseglength * length; + subcur = state.audittrail[subcur].previous; } } } - iter++; + }; + + for (const auto &connection : axline.m_back_segconns) { + processConnection(connection.first.ref); } - } - // also put in mean depth: - // - row.setValue(meandepthcol.c_str(), totalmetdepth / (total - 1)); - row.setValue(totaldcol.c_str(), totalmetdepth); - row.setValue(wmeandepthcol.c_str(), wtotaldepth / (wtotal - rootseglength)); - row.setValue(totalcol.c_str(), total); - row.setValue(wtotalcol.c_str(), wtotal); - // - if (comm) { - if (qtimer(atime, 500)) { - if (comm->IsCancelled()) { - throw Communicator::CancelledException(); - } + for (const auto &connection : axline.m_forward_segconns) { + processConnection(connection.first.ref); } - comm->CommPostMessage(Communicator::CURRENT_RECORD, reccount); } - reccount++; - } - if (!m_sel_only) { - // note, I've stopped sel only from calculating choice values: - for (size_t cursor = 0; cursor < map.getShapeCount(); cursor++) { - AttributeRow& row = map.getAttributeRowFromShapeIndex(cursor); - row.setValue(choicecol.c_str(), choicevals[cursor].choice); - row.setValue(wchoicecol.c_str(), choicevals[cursor].wchoice); + + SegmentMetricResult &result = results[cursor]; + result.totalDepth = totalmetdepth; + result.totalNodes = total; + result.totalLength = wtotal; + if (total > 1.0) { + result.meanDepth = totalmetdepth / (total - 1.0); + } + if (wtotal > rootseglength) { + result.weightedMeanDepth = wtotaldepth / (wtotal - rootseglength); } + }); + + for (size_t cursor : roots) { + if (comm && ((cursor & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + const SegmentMetricResult &result = results[cursor]; + AttributeRow &row = map.getAttributeRowFromShapeIndex(cursor); + row.setValue(meandepthcol, float(result.meanDepth)); + row.setValue(totaldcol, float(result.totalDepth)); + row.setValue(wmeandepthcol, float(result.weightedMeanDepth)); + row.setValue(totalcol, float(result.totalNodes)); + row.setValue(wtotalcol, float(result.totalLength)); } if (!m_sel_only) { - map.setDisplayedAttribute(attributes.getColumnIndex(choicecol.c_str())); + for (size_t cursor = 0; cursor < shapeCount; ++cursor) { + if (comm && ((cursor & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + double choice = 0.0; + double weightedChoice = 0.0; + for (size_t worker = 0; worker < workerCount; ++worker) { + choice += workerStates[worker].choice[cursor].choice; + weightedChoice += workerStates[worker].choice[cursor].wchoice; + } + AttributeRow &row = map.getAttributeRowFromShapeIndex(cursor); + row.setValue(choicecol, float(choice)); + row.setValue(wchoicecol, float(weightedChoice)); + } + map.setDisplayedAttribute(attributes.getColumnIndex(choicecol)); } else { - map.setDisplayedAttribute(attributes.getColumnIndex(meandepthcol.c_str())); + map.setDisplayedAttribute(attributes.getColumnIndex(meandepthcol)); } - return retvar; + return true; } diff --git a/salalib/segmmodules/segmtopological.cpp b/salalib/segmmodules/segmtopological.cpp index 7110ffb7..a6608964 100644 --- a/salalib/segmmodules/segmtopological.cpp +++ b/salalib/segmmodules/segmtopological.cpp @@ -18,196 +18,222 @@ #include "salalib/segmmodules/segmtopological.h" +#include "genlib/parallel.h" #include "genlib/stringutils.h" -bool SegmentTopological::run(Communicator *comm, ShapeGraph &map, bool) { +#include +#include +#include +#include +#include + +namespace { +struct SegmentTopologicalResult { + double meanDepth = -1.0; + double weightedMeanDepth = -1.0; + double totalDepth = -1.0; + double totalNodes = 0.0; + double totalLength = 0.0; +}; + +struct SegmentTopologicalWorkerState { + std::vector seen; + std::vector audittrail; + std::array, 512> bins; + std::vector choice; + + SegmentTopologicalWorkerState(size_t shapeCount, bool calculateChoice) + : seen(shapeCount), audittrail(shapeCount), choice(calculateChoice ? shapeCount : 0) {} +}; +} // namespace +bool SegmentTopological::run(Communicator *comm, ShapeGraph &map, bool) { AttributeTable &attributes = map.getAttributeTable(); + const size_t shapeCount = map.getShapeCount(); - bool retvar = true; - - time_t atime = 0; - - if (comm) { - qtimer(atime, 0); - comm->CommPostMessage(Communicator::NUM_RECORDS, - (m_sel_only ? map.getSelSet().size() : map.getConnections().size())); - } - int reccount = 0; - - // record axial line refs for topological analysis std::vector axialrefs; - // quick through to find the longest seg length + axialrefs.reserve(shapeCount); std::vector seglengths; - float maxseglength = 0.0f; - for (size_t cursor = 0; cursor < map.getShapeCount(); cursor++) { - AttributeRow& row = map.getAttributeRowFromShapeIndex(cursor); - axialrefs.push_back(row.getValue(attributes.getColumnIndex("Axial Line Ref"))); - seglengths.push_back(row.getValue(attributes.getColumnIndex("Segment Length"))); - if (seglengths.back() > maxseglength) { - maxseglength = seglengths.back(); + seglengths.reserve(shapeCount); + std::vector roots; + roots.reserve(shapeCount); + for (size_t cursor = 0; cursor < shapeCount; ++cursor) { + AttributeRow &row = map.getAttributeRowFromShapeIndex(cursor); + axialrefs.push_back(static_cast(row.getValue(attributes.getColumnIndex("Axial Line Ref")))); + const float length = row.getValue(attributes.getColumnIndex("Segment Length")); + seglengths.push_back(length); + if (!m_sel_only || row.isSelected()) { + roots.push_back(cursor); } } - std::string prefix, suffix; - int maxbin; - prefix = "Topological "; - maxbin = 2; + std::string suffix; if (m_radius != -1.0) { suffix = dXstring::formatString(m_radius, " R%.f metric"); } - std::string choicecol = prefix + "Choice" + suffix; - std::string wchoicecol = prefix + "Choice [SLW]" + suffix; - std::string meandepthcol = prefix + "Mean Depth" + suffix; - std::string wmeandepthcol = prefix + std::string("Mean Depth [SLW]") + suffix; - std::string totaldcol = prefix + "Total Depth" + suffix; - std::string totalcol = prefix + "Total Nodes" + suffix; - std::string wtotalcol = prefix + "Total Length" + suffix; - // + const std::string prefix = "Topological "; + const std::string choicecol = prefix + "Choice" + suffix; + const std::string wchoicecol = prefix + "Choice [SLW]" + suffix; + const std::string meandepthcol = prefix + "Mean Depth" + suffix; + const std::string wmeandepthcol = prefix + "Mean Depth [SLW]" + suffix; + const std::string totaldcol = prefix + "Total Depth" + suffix; + const std::string totalcol = prefix + "Total Nodes" + suffix; + const std::string wtotalcol = prefix + "Total Length" + suffix; + if (!m_sel_only) { - attributes.insertOrResetColumn(choicecol.c_str()); - attributes.insertOrResetColumn(wchoicecol.c_str()); + attributes.insertOrResetColumn(choicecol); + attributes.insertOrResetColumn(wchoicecol); } - attributes.insertOrResetColumn(meandepthcol.c_str()); - attributes.insertOrResetColumn(wmeandepthcol.c_str()); - attributes.insertOrResetColumn(totaldcol.c_str()); - attributes.insertOrResetColumn(totalcol.c_str()); - attributes.insertOrResetColumn(wtotalcol.c_str()); - // - std::vector seen(map.getShapeCount()); - std::vector audittrail(map.getShapeCount()); - std::vector choicevals(map.getShapeCount()); - for (size_t cursor = 0; cursor < map.getShapeCount(); cursor++) { - AttributeRow& row = map.getAttributeRowFromShapeIndex(cursor); - if (m_sel_only && !row.isSelected()) { - continue; - } - for (size_t i = 0; i < map.getShapeCount(); i++) { - seen[i] = 0xffffffff; + attributes.insertOrResetColumn(meandepthcol); + attributes.insertOrResetColumn(wmeandepthcol); + attributes.insertOrResetColumn(totaldcol); + attributes.insertOrResetColumn(totalcol); + attributes.insertOrResetColumn(wtotalcol); + + if (comm) { + comm->CommPostMessage(Communicator::NUM_RECORDS, static_cast(roots.size())); + } + + std::vector results(shapeCount); + const size_t estimatedBytesPerWorker = + shapeCount * (sizeof(unsigned int) + sizeof(TopoMetSegmentRef) + (!m_sel_only ? sizeof(TopoMetSegmentChoice) : 0)); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(roots.size(), estimatedBytesPerWorker); + + std::vector workerStates; + workerStates.reserve(workerCount); + for (size_t worker = 0; worker < workerCount; ++worker) { + workerStates.emplace_back(shapeCount, !m_sel_only); + } + + depthmapX::parallelFor(roots.size(), workerCount, comm, [&](size_t task, size_t workerIndex) { + const size_t cursor = roots[task]; + SegmentTopologicalWorkerState &state = workerStates[workerIndex]; + std::fill(state.seen.begin(), state.seen.end(), std::numeric_limits::max()); + for (auto &bin : state.bins) { + bin.clear(); } - std::vector list[512]; // 512 bins! + int bin = 0; - list[bin].push_back(cursor); - double rootseglength = seglengths[cursor]; - audittrail[cursor] = TopoMetSegmentRef(cursor, Connector::SEG_CONN_ALL, rootseglength * 0.5, -1); + state.bins[bin].push_back(static_cast(cursor)); + const double rootseglength = seglengths[cursor]; + state.audittrail[cursor] = + TopoMetSegmentRef(static_cast(cursor), Connector::SEG_CONN_ALL, rootseglength * 0.5, -1); + int open = 1; unsigned int segdepth = 0; - double total = 0.0, wtotal = 0.0, wtotaldepth = 0.0, totalsegdepth = 0.0, totalmetdepth = 0.0; + double total = 0.0; + double wtotal = 0.0; + double wtotaldepth = 0.0; + double totalmetdepth = 0.0; + size_t processed = 0; + while (open != 0) { - while (list[bin].size() == 0) { - bin++; - segdepth += 1; - if (bin == maxbin) { - bin = 0; - } + if (comm && ((processed++ & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + while (state.bins[bin].empty()) { + bin = (bin + 1) % 2; + ++segdepth; } - // - TopoMetSegmentRef &here = audittrail[list[bin].back()]; - list[bin].pop_back(); - open--; - // + + TopoMetSegmentRef &here = state.audittrail[state.bins[bin].back()]; + state.bins[bin].pop_back(); + --open; if (here.done) { continue; - } else { - here.done = true; } - // - double len = seglengths[here.ref]; - totalsegdepth += segdepth; - totalmetdepth += here.dist - len * 0.5; // preloaded with length ahead + here.done = true; + + const double len = seglengths[here.ref]; + totalmetdepth += segdepth; wtotal += len; wtotaldepth += len * segdepth; + total += 1.0; - total += 1; - // - Connector &axline = map.getConnections().at(here.ref); - int connected_cursor = -2; - - auto iter = axline.m_back_segconns.begin(); - bool backsegs = true; - - while (connected_cursor != -1) { - if (backsegs && iter == axline.m_back_segconns.end()) { - iter = axline.m_forward_segconns.begin(); - backsegs = false; - } - if (!backsegs && iter == axline.m_forward_segconns.end()) { - break; + const Connector &axline = map.getConnections().at(here.ref); + auto processConnection = [&](int connected_cursor) { + if (connected_cursor < 0 || static_cast(connected_cursor) >= shapeCount) { + return; } - - connected_cursor = iter->first.ref; - - if (seen[connected_cursor] > segdepth && static_cast(connected_cursor) != cursor) { - bool seenalready = (seen[connected_cursor] == 0xffffffff) ? false : true; - float length = seglengths[connected_cursor]; - int axialref = axialrefs[connected_cursor]; - audittrail[connected_cursor] = + if (state.seen[connected_cursor] > segdepth && static_cast(connected_cursor) != cursor) { + const bool seenalready = state.seen[connected_cursor] != std::numeric_limits::max(); + const float length = seglengths[connected_cursor]; + state.audittrail[connected_cursor] = TopoMetSegmentRef(connected_cursor, here.dir, here.dist + length, here.ref); - seen[connected_cursor] = segdepth; + state.seen[connected_cursor] = segdepth; if (m_radius == -1 || here.dist + length < m_radius) { - // puts in a suitable bin ahead of us... - open++; - // + ++open; + const int axialref = axialrefs[connected_cursor]; if (axialrefs[here.ref] == axialref) { - list[bin].push_back(connected_cursor); + state.bins[bin].push_back(connected_cursor); } else { - list[(bin + 1) % 2].push_back(connected_cursor); - seen[connected_cursor] = - segdepth + 1; // this is so if another node is connected directly to this one but - // is found later it is still handled -- note it can result in the - // connected cursor being added twice + state.bins[(bin + 1) % 2].push_back(connected_cursor); + state.seen[connected_cursor] = segdepth + 1; } } - // not sure why this is outside the radius restriction - // (sel_only: with restricted selection set, not all lines will be labelled) - // (seenalready: need to check that we're not doing this twice, given the seen can go twice) - - // Quick mod - TV - if (!m_sel_only && connected_cursor > int(cursor) && - !seenalready) { // only one way paths, saves doing this twice + if (!m_sel_only && connected_cursor > static_cast(cursor) && !seenalready) { int subcur = connected_cursor; while (subcur != -1) { - // in this method of choice, start and end lines are included - choicevals[subcur].choice += 1; - choicevals[subcur].wchoice += (rootseglength * length); - subcur = audittrail[subcur].previous; + state.choice[subcur].choice += 1.0; + state.choice[subcur].wchoice += rootseglength * length; + subcur = state.audittrail[subcur].previous; } } } - iter++; + }; + + for (const auto &connection : axline.m_back_segconns) { + processConnection(connection.first.ref); } - } - // also put in mean depth: - row.setValue(meandepthcol.c_str(), totalsegdepth / (total - 1)); - row.setValue(totaldcol.c_str(), totalsegdepth); - row.setValue(wmeandepthcol.c_str(), wtotaldepth / (wtotal - rootseglength)); - row.setValue(totalcol.c_str(), total); - row.setValue(wtotalcol.c_str(), wtotal); - // - if (comm) { - if (qtimer(atime, 500)) { - if (comm->IsCancelled()) { - throw Communicator::CancelledException(); - } + for (const auto &connection : axline.m_forward_segconns) { + processConnection(connection.first.ref); } - comm->CommPostMessage(Communicator::CURRENT_RECORD, reccount); } - reccount++; - } - if (!m_sel_only) { - // note, I've stopped sel only from calculating choice values: - for (size_t cursor = 0; cursor < map.getShapeCount(); cursor++) { - AttributeRow& row = map.getAttributeRowFromShapeIndex(cursor); - row.setValue(choicecol.c_str(), choicevals[cursor].choice); - row.setValue(wchoicecol.c_str(), choicevals[cursor].wchoice); + + SegmentTopologicalResult &result = results[cursor]; + result.totalDepth = totalmetdepth; + result.totalNodes = total; + result.totalLength = wtotal; + if (total > 1.0) { + result.meanDepth = totalmetdepth / (total - 1.0); + } + if (wtotal > rootseglength) { + result.weightedMeanDepth = wtotaldepth / (wtotal - rootseglength); } + }); + + for (size_t cursor : roots) { + if (comm && ((cursor & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + const SegmentTopologicalResult &result = results[cursor]; + AttributeRow &row = map.getAttributeRowFromShapeIndex(cursor); + row.setValue(meandepthcol, float(result.meanDepth)); + row.setValue(totaldcol, float(result.totalDepth)); + row.setValue(wmeandepthcol, float(result.weightedMeanDepth)); + row.setValue(totalcol, float(result.totalNodes)); + row.setValue(wtotalcol, float(result.totalLength)); } if (!m_sel_only) { - map.setDisplayedAttribute(attributes.getColumnIndex(choicecol.c_str())); + for (size_t cursor = 0; cursor < shapeCount; ++cursor) { + if (comm && ((cursor & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + double choice = 0.0; + double weightedChoice = 0.0; + for (size_t worker = 0; worker < workerCount; ++worker) { + choice += workerStates[worker].choice[cursor].choice; + weightedChoice += workerStates[worker].choice[cursor].wchoice; + } + AttributeRow &row = map.getAttributeRowFromShapeIndex(cursor); + row.setValue(choicecol, float(choice)); + row.setValue(wchoicecol, float(weightedChoice)); + } + map.setDisplayedAttribute(attributes.getColumnIndex(choicecol)); } else { - map.setDisplayedAttribute(attributes.getColumnIndex(meandepthcol.c_str())); + map.setDisplayedAttribute(attributes.getColumnIndex(meandepthcol)); } - return retvar; + return true; } diff --git a/salalib/shapemap.cpp b/salalib/shapemap.cpp index ed444071..608bd4b8 100644 --- a/salalib/shapemap.cpp +++ b/salalib/shapemap.cpp @@ -2061,14 +2061,15 @@ int ShapeMap::connectIntersected(int rowid, bool linegraph) { // use the other version, getShapeConnections for arbitrary shape-shape connections // note, connections are listed by rowid in list, *not* reference number // (so they may vary: must be checked carefully when shapes are removed / added) -std::vector ShapeMap::getLineConnections(int lineref, double tolerance) { +std::vector ShapeMap::getLineConnections( + int lineref, double tolerance, const std::unordered_map *shapeIndexLookup) const { std::vector connections; - SalaShape &poly = m_shapes.find(lineref)->second; - if (!poly.isLine()) { - return std::vector(); + const auto polyIter = m_shapes.find(lineref); + if (polyIter == m_shapes.end() || !polyIter->second.isLine()) { + return connections; } - const Line &l = poly.getLine(); + const Line &l = polyIter->second.getLine(); std::unordered_set shapesToTest; @@ -2091,13 +2092,27 @@ std::vector ShapeMap::getLineConnections(int lineref, double tolerance) { } for (ShapeRef shape : shapesToTest) { if ((shape.m_tags & ShapeRef::SHAPE_OPEN) == ShapeRef::SHAPE_OPEN) { - const Line &line = m_shapes.find(int(shape.m_shape_ref))->second.getLine(); + const auto shapeIter = m_shapes.find(static_cast(shape.m_shape_ref)); + if (shapeIter == m_shapes.end()) { + continue; + } + const Line &line = shapeIter->second.getLine(); if (intersect_region(line, l, line.length() * tolerance)) { // n.b. originally this followed the logic that we must normalise intersect_line properly: tolerance * // line length one * line length two in fact, works better if it's just line.length() * tolerance... if (intersect_line(line, l, line.length() * tolerance)) { - depthmapX::insert_sorted(connections, - depthmapX::findIndexFromKey(m_shapes, int(shape.m_shape_ref))); + int shapeIndex = -1; + if (shapeIndexLookup) { + const auto indexIter = shapeIndexLookup->find(static_cast(shape.m_shape_ref)); + if (indexIter != shapeIndexLookup->end()) { + shapeIndex = indexIter->second; + } + } else { + shapeIndex = depthmapX::findIndexFromKey(m_shapes, int(shape.m_shape_ref)); + } + if (shapeIndex >= 0) { + depthmapX::insert_sorted(connections, shapeIndex); + } } } } diff --git a/salalib/shapemap.h b/salalib/shapemap.h index 01847716..ad5cef53 100644 --- a/salalib/shapemap.h +++ b/salalib/shapemap.h @@ -36,6 +36,7 @@ #include #include #include +#include #include // each pixel has various lists of information: @@ -413,7 +414,9 @@ class ShapeMap : public PixelBase { // Connect a particular shape into the graph int connectIntersected(int rowid, bool linegraph); // Get the connections for a particular line - std::vector getLineConnections(int lineref, double tolerance); + std::vector getLineConnections( + int lineref, double tolerance, + const std::unordered_map *shapeIndexLookup = nullptr) const; // Get arbitrary shape connections for a particular shape std::vector getShapeConnections(int polyref, double tolerance); // Make all connections diff --git a/version.h b/version.h index 0297fcd6..c793f474 100644 --- a/version.h +++ b/version.h @@ -22,7 +22,7 @@ // use these to define the depthmap versions #define DEPTHMAPX_MAJOR_VERSION 0 #define DEPTHMAPX_MINOR_VERSION 8 -#define DEPTHMAPX_REVISION_VERSION 1 +#define DEPTHMAPX_REVISION_VERSION 2 #define DEPTHMAP_MODULE_VERSION 10.04 From 564fce9da472916f62fd385ae5bfb80b0e9e54fe Mon Sep 17 00:00:00 2001 From: Victor Lin Date: Sun, 19 Jul 2026 13:50:00 +0800 Subject: [PATCH 3/3] Multithread Tulip Analysis --- depthmapX/UI/SegmentAnalysisDlg.ui | 4 +- salaTest/testparallelanalysis.cpp | 86 +++++++--- salalib/segmmodules/segmangular.cpp | 7 +- salalib/segmmodules/segmtulip.cpp | 254 ++++++++++++++++++++++++++++ 4 files changed, 318 insertions(+), 33 deletions(-) diff --git a/depthmapX/UI/SegmentAnalysisDlg.ui b/depthmapX/UI/SegmentAnalysisDlg.ui index 6ccf5252..5c345b91 100644 --- a/depthmapX/UI/SegmentAnalysisDlg.ui +++ b/depthmapX/UI/SegmentAnalysisDlg.ui @@ -23,7 +23,7 @@ - Tulip Analysis + Tulip Analysis (Multi-Threaded) @@ -58,7 +58,7 @@ - Include choice (betweenness) + Include choice (betweenness — Serial / Slower) diff --git a/salaTest/testparallelanalysis.cpp b/salaTest/testparallelanalysis.cpp index f523603e..867beab6 100644 --- a/salaTest/testparallelanalysis.cpp +++ b/salaTest/testparallelanalysis.cpp @@ -14,6 +14,7 @@ #include "salalib/mgraph.h" #include "salalib/segmmodules/segmangular.h" #include "salalib/segmmodules/segmmetric.h" +#include "salalib/segmmodules/segmtulip.h" #include "salalib/segmmodules/segmtopological.h" #include @@ -181,29 +182,6 @@ void requireSameAttributes(const ShapeGraph &serial, const ShapeGraph ¶llel) } } -void requireAngularIntegration(const ShapeGraph &graph, const std::vector &radiusSuffixes) { - const AttributeTable &attributes = graph.getAttributeTable(); - for (const std::string &suffix : radiusSuffixes) { - const int countColumn = attributes.getColumnIndex(std::string("Angular Node Count") + suffix); - const int totalDepthColumn = attributes.getColumnIndex(std::string("Angular Total Depth") + suffix); - const int integrationColumn = attributes.getColumnIndex(std::string("Angular Integration") + suffix); - REQUIRE(countColumn != -1); - REQUIRE(totalDepthColumn != -1); - REQUIRE(integrationColumn != -1); - - for (size_t rowIndex = 0; rowIndex < graph.getShapeCount(); ++rowIndex) { - const AttributeRow &row = graph.getAttributeRowFromShapeIndex(rowIndex); - const double count = row.getValue(countColumn); - const double totalDepth = row.getValue(totalDepthColumn); - const double integration = row.getValue(integrationColumn); - if (totalDepth > 1e-9) { - REQUIRE(integration == Approx(count * count / totalDepth).epsilon(1e-5)); - } else { - REQUIRE(integration == Approx(-1.0)); - } - } - } -} } // namespace @@ -241,8 +219,6 @@ TEST_CASE("Parallel segment analyses match single-thread analyses") { } requireSameAttributes(*serial, *parallel); - requireAngularIntegration(*serial, {" R1.00", ""}); - requireAngularIntegration(*parallel, {" R1.00", ""}); } TEST_CASE("Parallel line connection discovery matches single-thread discovery") { @@ -332,3 +308,63 @@ TEST_CASE("Parallel axial-to-segment conversion matches single-thread conversion requireSameSegmentConnections(*serial, *parallel); requireSameAttributes(*serial, *parallel); } + + +TEST_CASE("Parallel Tulip analysis without Choice matches single-thread analysis") { + std::unique_ptr serial = makeSegmentGraph(); + std::unique_ptr parallel = makeSegmentGraph(); + + const int serialWeight = serial->getAttributeTable().insertOrResetColumn("Tulip Weight"); + const int parallelWeight = parallel->getAttributeTable().insertOrResetColumn("Tulip Weight"); + const int serialRouteWeight = serial->getAttributeTable().insertOrResetColumn("Tulip Route Weight"); + const int parallelRouteWeight = parallel->getAttributeTable().insertOrResetColumn("Tulip Route Weight"); + for (size_t row = 0; row < serial->getShapeCount(); ++row) { + const float weight = static_cast(1 + (row % 5)); + const float routeWeight = static_cast(1 + (row % 7)); + serial->getAttributeRowFromShapeIndex(row).setValue(serialWeight, weight); + parallel->getAttributeRowFromShapeIndex(row).setValue(parallelWeight, weight); + serial->getAttributeRowFromShapeIndex(row).setValue(serialRouteWeight, routeWeight); + parallel->getAttributeRowFromShapeIndex(row).setValue(parallelRouteWeight, routeWeight); + } + + { + AnalysisThreadSetting threads("1"); + REQUIRE(SegmentTulip({1.0, 2.0, -1.0}, false, 1024, serialWeight, Options::RADIUS_STEPS, false) + .run(nullptr, *serial, false)); + REQUIRE(SegmentTulip({0.5, 1.0, -1.0}, false, 1024, serialWeight, Options::RADIUS_ANGULAR, false) + .run(nullptr, *serial, false)); + REQUIRE(SegmentTulip({2.0, 5.0, -1.0}, false, 1024, serialWeight, Options::RADIUS_METRIC, false, false, -1, + serialRouteWeight) + .run(nullptr, *serial, false)); + } + { + AnalysisThreadSetting threads("4"); + REQUIRE(SegmentTulip({1.0, 2.0, -1.0}, false, 1024, parallelWeight, Options::RADIUS_STEPS, false) + .run(nullptr, *parallel, false)); + REQUIRE(SegmentTulip({0.5, 1.0, -1.0}, false, 1024, parallelWeight, Options::RADIUS_ANGULAR, false) + .run(nullptr, *parallel, false)); + REQUIRE(SegmentTulip({2.0, 5.0, -1.0}, false, 1024, parallelWeight, Options::RADIUS_METRIC, false, false, -1, + parallelRouteWeight) + .run(nullptr, *parallel, false)); + } + + requireSameAttributes(*serial, *parallel); +} + +TEST_CASE("Tulip Choice remains serial and matches regardless of thread setting") { + std::unique_ptr serial = makeSegmentGraph(); + std::unique_ptr requestedParallel = makeSegmentGraph(); + + { + AnalysisThreadSetting threads("1"); + REQUIRE(SegmentTulip({1.0, -1.0}, false, 1024, -1, Options::RADIUS_STEPS, true) + .run(nullptr, *serial, false)); + } + { + AnalysisThreadSetting threads("4"); + REQUIRE(SegmentTulip({1.0, -1.0}, false, 1024, -1, Options::RADIUS_STEPS, true) + .run(nullptr, *requestedParallel, false)); + } + + requireSameAttributes(*serial, *requestedParallel); +} diff --git a/salalib/segmmodules/segmangular.cpp b/salalib/segmmodules/segmangular.cpp index 6b30cf76..9c04e3d9 100644 --- a/salalib/segmmodules/segmangular.cpp +++ b/salalib/segmmodules/segmangular.cpp @@ -49,13 +49,12 @@ bool SegmentAngular::run(Communicator *comm, ShapeGraph &map, bool) { radii.push_back(-1.0); } - std::vector depth_col, count_col, total_col, integration_col; + std::vector depth_col, count_col, total_col; for (double radius : radii) { std::string radius_text = makeRadiusText(Options::RADIUS_ANGULAR, radius); depth_col.push_back(attributes.insertOrResetColumn(std::string("Angular Mean Depth") + radius_text)); count_col.push_back(attributes.insertOrResetColumn(std::string("Angular Node Count") + radius_text)); total_col.push_back(attributes.insertOrResetColumn(std::string("Angular Total Depth") + radius_text)); - integration_col.push_back(attributes.insertOrResetColumn(std::string("Angular Integration") + radius_text)); } if (comm) { @@ -158,13 +157,9 @@ bool SegmentAngular::run(Communicator *comm, ShapeGraph &map, bool) { if (nodeCount > 1) { row.setValue(depth_col[radiusIndex], float(totalDepth / double(nodeCount - 1))); row.setValue(total_col[radiusIndex], float(totalDepth)); - row.setValue(integration_col[radiusIndex], - totalDepth > 1e-9 ? float(double(nodeCount) * double(nodeCount) / totalDepth) - : -1.0f); } else { row.setValue(depth_col[radiusIndex], -1.0f); row.setValue(total_col[radiusIndex], -1.0f); - row.setValue(integration_col[radiusIndex], -1.0f); } } ++root; diff --git a/salalib/segmmodules/segmtulip.cpp b/salalib/segmmodules/segmtulip.cpp index 9f0f6fee..81c39e78 100644 --- a/salalib/segmmodules/segmtulip.cpp +++ b/salalib/segmmodules/segmtulip.cpp @@ -18,8 +18,11 @@ #include "salalib/segmmodules/segmtulip.h" +#include "genlib/parallel.h" #include "genlib/stringutils.h" +#include + bool SegmentTulip::run(Communicator *comm, ShapeGraph &map, bool) { if (map.getMapType() != ShapeMap::SEGMENTMAP) { @@ -341,6 +344,257 @@ bool SegmentTulip::run(Communicator *comm, ShapeGraph &map, bool) { radiusmask |= (1 << i); } + // Choice accumulates data across all origins and therefore keeps the original + // serial implementation below. Without Choice, each origin is independent and can + // safely run with private traversal state on a worker thread. + if (!m_choice) { + const size_t shapeCount = map.getConnections().size(); + std::vector roots; + roots.reserve(shapeCount); + for (size_t cursor = 0; cursor < shapeCount; ++cursor) { + if (!m_sel_only || map.getAttributeRowFromShapeIndex(cursor).isSelected()) { + roots.push_back(cursor); + } + } + + struct TulipResult { + double nodeCount = 0.0; + double totalDepth = 0.0; + double totalWeight = 0.0; + double totalWeightedDepth = 0.0; + }; + + struct TulipWorkerState { + std::vector> bins; + std::vector uncovered; + std::vector depths; + + TulipWorkerState(size_t binCount, size_t shapeCount, size_t radiusCount) + : bins(binCount), uncovered(shapeCount * 2), depths(shapeCount * radiusCount * 2) {} + }; + + const size_t resultStride = static_cast(radiussize); + std::vector results(roots.size() * resultStride); + + const size_t bytesPerWorker = + shapeCount * 2 * sizeof(unsigned int) + + shapeCount * resultStride * 2 * sizeof(int) + + static_cast(tulip_bins) * sizeof(std::vector); + const size_t workerCount = depthmapX::configuredAnalysisThreadCount(roots.size(), bytesPerWorker); + + std::vector workerStates; + workerStates.reserve(workerCount); + for (size_t worker = 0; worker < workerCount; ++worker) { + workerStates.emplace_back(static_cast(tulip_bins), shapeCount, resultStride); + } + + auto depthIndex = [resultStride](size_t ref, size_t radiusIndex, int dir) { + return ((ref * resultStride + radiusIndex) * 2) + static_cast(dir); + }; + auto uncoveredIndex = [](size_t ref, int dir) { + return ref * 2 + static_cast(dir); + }; + + depthmapX::parallelFor(roots.size(), workerCount, comm, [&](size_t task, size_t workerIndex) { + const size_t cursor = roots[task]; + TulipWorkerState &state = workerStates[workerIndex]; + + for (std::vector &bin : state.bins) { + bin.clear(); + } + std::fill(state.uncovered.begin(), state.uncovered.end(), static_cast(radiusmask)); + std::fill(state.depths.begin(), state.depths.end(), 0); + + const double rootseglength = lengths[cursor]; + SegmentData rootData(0, static_cast(cursor), SegmentRef(), 0, + static_cast(0.5 * rootseglength), radiusmask); + auto rootIter = std::lower_bound(state.bins[0].begin(), state.bins[0].end(), rootData); + if (rootIter == state.bins[0].end() || rootData != *rootIter) { + state.bins[0].insert(rootIter, rootData); + } + + int depthlevel = 0; + int opencount = 1; + size_t currentbin = 0; + size_t processed = 0; + while (opencount) { + if (comm && ((processed++ & 0x3ff) == 0) && comm->IsCancelled()) { + throw Communicator::CancelledException(); + } + + while (state.bins[currentbin].empty()) { + ++depthlevel; + ++currentbin; + if (currentbin == static_cast(tulip_bins)) { + currentbin = 0; + } + } + + SegmentData lineindex = state.bins[currentbin].back(); + state.bins[currentbin].pop_back(); + --opencount; + + const int ref = lineindex.ref; + const int dir = (lineindex.dir == 1) ? 0 : 1; + const unsigned int coverage = + lineindex.coverage & state.uncovered[uncoveredIndex(static_cast(ref), dir)]; + if (coverage == 0) { + continue; + } + + int rbin = 0; + int rbinbase = 0; + if (lineindex.previous.ref != -1) { + state.uncovered[uncoveredIndex(static_cast(ref), dir)] &= ~coverage; + while (((coverage >> rbin) & 0x1U) == 0U) { + ++rbin; + } + rbinbase = rbin; + while (rbin < radiussize) { + if (((coverage >> rbin) & 0x1U) == 1U) { + state.depths[depthIndex(static_cast(ref), static_cast(rbin), dir)] = + depthlevel; + } + ++rbin; + } + } else { + state.uncovered[uncoveredIndex(static_cast(ref), 0)] &= ~coverage; + state.uncovered[uncoveredIndex(static_cast(ref), 1)] &= ~coverage; + } + + const Connector &line = map.getConnections()[static_cast(ref)]; + auto addConnections = [&](const std::map &connections, char previousDirection) { + for (const auto &segconn : connections) { + rbin = rbinbase; + const SegmentRef conn = segconn.first; + const int connDir = (conn.dir == 1) ? 0 : 1; + if ((state.uncovered[uncoveredIndex(static_cast(conn.ref), connDir)] & coverage) == 0) { + continue; + } + + const int extradepth = routeweight_col != -1 + ? static_cast(floor(segconn.second * tulip_bins * 0.5 * + routeweights[static_cast(conn.ref)])) + : static_cast(floor(segconn.second * tulip_bins * 0.5)); + const float seglength = lengths[static_cast(conn.ref)]; + + switch (m_radius_type) { + case Options::RADIUS_ANGULAR: + while (rbin != radiussize && radius[static_cast(rbin)] != -1 && + depthlevel + extradepth > static_cast(radius[static_cast(rbin)])) { + ++rbin; + } + break; + case Options::RADIUS_METRIC: + while (rbin != radiussize && radius[static_cast(rbin)] != -1 && + lineindex.metricdepth + seglength * 0.5 > radius[static_cast(rbin)]) { + ++rbin; + } + break; + case Options::RADIUS_STEPS: + if (rbin != radiussize && radius[static_cast(rbin)] != -1 && + lineindex.segdepth >= static_cast(radius[static_cast(rbin)])) { + ++rbin; + } + break; + } + + if ((coverage >> rbin) != 0U) { + SegmentData next(conn, SegmentRef(previousDirection, lineindex.ref), + lineindex.segdepth + 1, lineindex.metricdepth + seglength, + (coverage >> rbin) << rbin); + const size_t bin = + (currentbin + static_cast(tulip_bins) + static_cast(extradepth)) % + static_cast(tulip_bins); + depthmapX::insert_sorted(state.bins[bin], next); + ++opencount; + } + } + }; + + if (lineindex.dir != -1) { + addConnections(line.m_forward_segconns, 1); + } + if (lineindex.dir != 1) { + addConnections(line.m_back_segconns, -1); + } + } + + for (int radiusIndex = 0; radiusIndex < radiussize; ++radiusIndex) { + TulipResult result; + for (size_t ref = 0; ref < shapeCount; ++ref) { + const bool reached0 = + ((state.uncovered[uncoveredIndex(ref, 0)] >> radiusIndex) & 0x1U) == 0U; + const bool reached1 = + ((state.uncovered[uncoveredIndex(ref, 1)] >> radiusIndex) & 0x1U) == 0U; + if (!reached0 && !reached1) { + continue; + } + + int reachedDir; + if (reached0 && reached1) { + const int depth0 = state.depths[depthIndex(ref, static_cast(radiusIndex), 0)]; + const int depth1 = state.depths[depthIndex(ref, static_cast(radiusIndex), 1)]; + reachedDir = depth0 < depth1 ? 0 : 1; + } else { + reachedDir = reached0 ? 0 : 1; + } + + const int reachedDepth = + state.depths[depthIndex(ref, static_cast(radiusIndex), reachedDir)]; + result.nodeCount += 1.0; + result.totalDepth += reachedDepth; + result.totalWeight += weights[ref]; + result.totalWeightedDepth += reachedDepth * weights[ref]; + } + + result.totalDepth /= ((tulip_bins - 1.0f) * 0.5f); + result.totalWeightedDepth /= ((tulip_bins - 1.0f) * 0.5f); + results[task * resultStride + static_cast(radiusIndex)] = result; + } + }); + + for (size_t task = 0; task < roots.size(); ++task) { + AttributeRow &row = map.getAttributeRowFromShapeIndex(roots[task]); + for (size_t radiusIndex = 0; radiusIndex < resultStride; ++radiusIndex) { + const TulipResult &result = results[task * resultStride + radiusIndex]; + row.setValue(count_col[radiusIndex], static_cast(result.nodeCount)); + if (result.nodeCount > 1.0) { + row.setValue(td_col[radiusIndex], static_cast(result.totalDepth)); + if (m_weighted_measure_col != -1) { + row.setValue(total_weight_col[radiusIndex], static_cast(result.totalWeight)); + row.setValue(w_td_col[radiusIndex], static_cast(result.totalWeightedDepth)); + } + } else { + row.setValue(td_col[radiusIndex], -1.0f); + if (m_weighted_measure_col != -1) { + row.setValue(total_weight_col[radiusIndex], -1.0f); + row.setValue(w_td_col[radiusIndex], -1.0f); + } + } + + if (result.totalDepth > 1e-9) { + row.setValue(integ_col[radiusIndex], + static_cast(result.nodeCount * result.nodeCount / result.totalDepth)); + if (m_weighted_measure_col != -1) { + row.setValue(w_integ_col[radiusIndex], + static_cast(result.totalWeight * result.totalWeight / + result.totalWeightedDepth)); + } + } else { + row.setValue(integ_col[radiusIndex], -1.0f); + if (m_weighted_measure_col != -1) { + row.setValue(w_integ_col[radiusIndex], -1.0f); + } + } + } + } + + map.setDisplayedAttribute(-2); + map.setDisplayedAttribute(td_col.back()); + return !roots.empty(); + } + for (size_t cursor = 0; cursor < map.getConnections().size(); cursor++) { AttributeRow &row = map.getAttributeRowFromShapeIndex(cursor);