Speed up map matching up to 23x

CMakeLists.txt

@@ -1,78 +1,84 @@
-cmake_minimum_required( VERSION 3.5.1)
+cmake_minimum_required(VERSION 3.16...4.4)
 
 message(STATUS "CMAKE version ${CMAKE_VERSION}")
 
-if (POLICY CMP0074)
-    message(STATUS "Set CMP0074 state to NEW")
-    cmake_policy(SET CMP0074 NEW)
-endif()
-
-if (POLICY CMP0086)
-    message(STATUS "Set CMP0086 state to NEW")
-    cmake_policy(SET CMP0086 NEW)
-endif()
-
-if (POLICY CMP0078)
-    message(STATUS "Set CMP0078 state to NEW")
-    cmake_policy(SET CMP0078 NEW)
-endif()
-
 # Prevent in source build
 set(CMAKE_DISABLE_SOURCE_CHANGES  OFF)
 set(CMAKE_DISABLE_IN_SOURCE_BUILD ON)
 set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS ON)
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 
-project(fmm)
+project(fmm LANGUAGES C CXX)
 
-set(CMAKE_BUILD_TYPE "Release")
-set(CMAKE_CXX_FLAGS "-O3 -DSPDLOG_ACTIVE_LEVEL=SPDLOG_LEVEL_TRACE")
-set(CMAKE_CXX_STANDARD 11)
-set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/build")
+# Default to an optimised build. An unknown build type would silently drop
+# the optimisation flags, so anything but the standard types is reset too.
+set(FMM_KNOWN_BUILD_TYPES Debug Release RelWithDebInfo MinSizeRel)
+if(NOT CMAKE_BUILD_TYPE IN_LIST FMM_KNOWN_BUILD_TYPES)
+  if(CMAKE_BUILD_TYPE)
+    message(STATUS "Unknown build type '${CMAKE_BUILD_TYPE}', using Release")
+  endif()
+  set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Build type" FORCE)
+endif()
+message(STATUS "Build type ${CMAKE_BUILD_TYPE}")
+
+# Boost.Geometry >= 1.75 requires C++14 and GDAL >= 3.9 headers require
+# C++17, so C++17 is the floor for building against a current toolchain.
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+# Compile-time log level of spdlog. Messages below this level are removed
+# from the binary entirely (SPDLOG_LEVEL_TRACE keeps everything and lets the
+# CLI --log_level flag switch them on at run time).
+set(FMM_SPDLOG_ACTIVE_LEVEL "SPDLOG_LEVEL_TRACE" CACHE STRING
+  "Compile-time spdlog level (SPDLOG_LEVEL_TRACE, SPDLOG_LEVEL_DEBUG, SPDLOG_LEVEL_INFO, ...)")
+add_compile_definitions(SPDLOG_ACTIVE_LEVEL=${FMM_SPDLOG_ACTIVE_LEVEL})
+
+# Optimise for the CPU of the build machine. The resulting binaries only run
+# on CPUs supporting the same instruction set extensions, so this is OFF by
+# default.
+option(FMM_NATIVE_ARCH "Compile with -march=native" OFF)
+if(FMM_NATIVE_ARCH)
+  add_compile_options(-march=native)
+endif()
 
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake")
 
 find_package(Conda)
 if(CONDA_PREFIX)
   message(STATUS "Set CONDA_PREFIX ${CONDA_PREFIX}")
-  set(CMAKE_PREFIX_PATH ${CMAKE_PREFIX_PATH} ${CONDA_PREFIX})
+  list(APPEND CMAKE_PREFIX_PATH ${CONDA_PREFIX})
 else()
   message(STATUS "Non conda exist, search library in default path")
 endif()
 
 find_package(GDAL 2.2 REQUIRED)
-if (GDAL_FOUND)
-  message(STATUS "GDAL headers found at ${GDAL_INCLUDE_DIR}")
-  message(STATUS "GDAL library found at ${GDAL_LIBRARIES}")
-  include_directories(${GDAL_INCLUDE_DIR})
-else()
-  message(FATAL_ERROR "GDAL Not Found!")
-endif (GDAL_FOUND)
-link_libraries(${GDAL_LIBRARIES})
+message(STATUS "GDAL ${GDAL_VERSION} headers found at ${GDAL_INCLUDE_DIRS}")
+message(STATUS "GDAL library found at ${GDAL_LIBRARIES}")
 
 if (WIN32)
-  find_package(Boost 1.56.0 REQUIRED serialization exception)
+  find_package(Boost 1.70 REQUIRED COMPONENTS serialization exception)
 else ()
-  find_package(Boost 1.56.0 REQUIRED serialization)
+  find_package(Boost 1.70 REQUIRED COMPONENTS serialization)
 endif (WIN32)
-if (Boost_FOUND)
-  message(STATUS "Boost headers found at ${Boost_INCLUDE_DIR}")
-  message(STATUS "Boost library found at ${Boost_LIBRARIES}")
-  message(STATUS "Boost library version ${Boost_LIB_VERSION}")
+message(STATUS "Boost ${Boost_VERSION} headers found at ${Boost_INCLUDE_DIRS}")
+message(STATUS "Boost library found at ${Boost_LIBRARIES}")
+if(TARGET Boost::headers)
+  set(FMM_BOOST_HEADERS Boost::headers)
 else()
-  message(FATAL_ERROR "Boost Not Found!")
-endif (Boost_FOUND)
-link_libraries(${Boost_LIBRARIES})
-
+  set(FMM_BOOST_HEADERS Boost::boost)
+endif()
 
 find_package(OpenMP REQUIRED)
-if(OPENMP_FOUND)
-  set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
-  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
-  message(STATUS "OpenMP_HEADERS found at ${OpenMP_CXX_INCLUDE_DIRS}")
-  message(STATUS "OpenMP_CXX_LIBRARIES found at ${OpenMP_CXX_LIBRARIES}")
-endif()
-link_libraries(${OpenMP_CXX_LIBRARIES})
+message(STATUS "OpenMP ${OpenMP_CXX_VERSION} flags ${OpenMP_CXX_FLAGS}")
+
+# Every target below (object libraries, FMMLIB, executables and the python
+# module in the subdirectory) picks up include directories, compile flags and
+# link libraries from these imported targets.
+link_libraries(GDAL::GDAL Boost::serialization ${FMM_BOOST_HEADERS}
+  OpenMP::OpenMP_CXX)
+if (WIN32)
+  link_libraries(Boost::exception)
+endif (WIN32)
 
 ### Set RPATH properties
 

python/CMakeLists.txt

@@ -1,43 +1,31 @@
 # Make sure the swig package is loaded.
-find_package(SWIG REQUIRED)
+find_package(SWIG 4.0 REQUIRED)
 include(${SWIG_USE_FILE})
-if (SWIG_FOUND)
-  message(STATUS "Swig version is ${SWIG_VERSION}")
-else()
-  message(FATAL_ERROR "Swig not found!\n")
-endif()
+message(STATUS "Swig version is ${SWIG_VERSION}")
 
-find_package(PythonInterp REQUIRED)
-find_package(PythonLibs REQUIRED)
-if (PYTHONLIBS_FOUND)
-  message(STATUS "Python header found at ${PYTHON_INCLUDE_DIRS}")
-  message(STATUS "Python library found at ${PYTHON_LIBRARIES}")
-  include_directories(${PYTHON_INCLUDE_DIRS})
-else()
-  message(FATAL_ERROR "Python library not found!\n")
-endif()
+# Development.Module locates the headers needed to build an extension module
+# without requiring a shared libpython (which is the recommended way of
+# building CPython extensions and also works for statically linked Pythons).
+find_package(Python3 COMPONENTS Interpreter Development.Module REQUIRED)
+message(STATUS "Python ${Python3_VERSION} interpreter found at ${Python3_EXECUTABLE}")
+message(STATUS "Python header found at ${Python3_INCLUDE_DIRS}")
 
-execute_process(
-  COMMAND python -c "from distutils.sysconfig import get_python_lib;print(get_python_lib())"
-  OUTPUT_VARIABLE PYTHON_SITE_PACKAGES OUTPUT_STRIP_TRAILING_WHITESPACE)
-message(STATUS "Python packages ${PYTHON_SITE_PACKAGES}")
+set(FMM_PYTHON_INSTALL_DIR "${Python3_SITEARCH}" CACHE PATH
+  "Directory where the fmm python module is installed")
+message(STATUS "Python packages ${FMM_PYTHON_INSTALL_DIR}")
 
-set_source_files_properties(
-${PROJECT_SOURCE_DIR}/python/fmm.i PROPERTIES CPLUSPLUS ON)
+set_property(SOURCE ${PROJECT_SOURCE_DIR}/python/fmm.i
+  PROPERTY CPLUSPLUS ON)
 
-if (${CMAKE_VERSION} VERSION_LESS 3.13.0)
-  message(STATUS "Using swig add module")
-  SWIG_ADD_MODULE(fmm python ${PROJECT_SOURCE_DIR}/python/fmm.i)
-  swig_link_libraries(fmm FMMLIB ${PYTHON_LIBRARIES})
-  install(TARGETS ${SWIG_MODULE_fmm_REAL_NAME} DESTINATION ${PYTHON_SITE_PACKAGES})
-else()
-  message(STATUS "Using swig add library")
-  SWIG_ADD_LIBRARY(pyfmm
-    LANGUAGE python
-    SOURCES ${PROJECT_SOURCE_DIR}/python/fmm.i)
-  set_property(TARGET pyfmm PROPERTY OUTPUT_NAME fmm)
-  swig_link_libraries(pyfmm FMMLIB ${PYTHON_LIBRARIES})
-  install(TARGETS pyfmm DESTINATION ${PYTHON_SITE_PACKAGES})
-endif()
+swig_add_library(pyfmm
+  LANGUAGE python
+  SOURCES ${PROJECT_SOURCE_DIR}/python/fmm.i)
+set_target_properties(pyfmm PROPERTIES
+  OUTPUT_NAME fmm
+  SWIG_USE_TARGET_INCLUDE_DIRECTORIES TRUE)
+target_include_directories(pyfmm PRIVATE ${FMM_INCLUDE_PATH})
+target_link_libraries(pyfmm PRIVATE FMMLIB Python3::Module)
 
-install(FILES ${CMAKE_BINARY_DIR}/python/fmm.py DESTINATION ${PYTHON_SITE_PACKAGES})
+install(TARGETS pyfmm LIBRARY DESTINATION ${FMM_PYTHON_INSTALL_DIR})
+install(FILES $<TARGET_PROPERTY:pyfmm,SWIG_SUPPORT_FILES>
+  DESTINATION ${FMM_PYTHON_INSTALL_DIR})

src/algorithm/geom_algorithm.cpp

@@ -154,8 +154,10 @@ void FMM::ALGORITHM::boundingbox_geometry(
   int Npoints = linestring.get_num_points();
   *x1 = DBL_MAX;
   *y1 = DBL_MAX;
-  *x2 = DBL_MIN;
-  *y2 = DBL_MIN;
+  // DBL_MIN is the smallest positive double, so it is not a valid start
+  // value for a maximum when coordinates can be negative.
+  *x2 = -DBL_MAX;
+  *y2 = -DBL_MAX;
   double x, y;
   for (int i = 0; i < Npoints; ++i) {
     x = linestring.get_x(i);

src/io/gps_reader.cpp

@@ -8,9 +8,11 @@
 #include "util/debug.hpp"
 #include "util/util.hpp"
 #include "config/gps_config.hpp"
+#include <charconv>
 #include <iostream>
 #include <stdexcept>
 #include <string>
+#include <system_error>
 
 #include <boost/format.hpp>
 
@@ -19,6 +21,41 @@ using namespace FMM::CORE;
 using namespace FMM::IO;
 using namespace FMM::UTIL;
 
+namespace {
+/**
+ * Convert line[start, end) to int exactly like std::stoi does with the
+ * token. std::from_chars handles the common case without allocating, and
+ * anything it does not fully accept (leading blanks or '+', trailing
+ * characters, out of range values) is delegated to std::stoi so that the
+ * value and the exceptions are the same as before.
+ */
+int parse_int_field(const std::string &line, std::size_t start,
+                    std::size_t end) {
+  int value = 0;
+  const char *first = line.data() + start;
+  const char *last = line.data() + end;
+  auto result = std::from_chars(first, last, value);
+  if (result.ec == std::errc() && result.ptr == last) return value;
+  return std::stoi(line.substr(start, end - start));
+}
+
+/**
+ * Convert line[start, end) to float exactly like std::stof does with the
+ * token (both are correctly rounded), see parse_int_field for the fallback.
+ * The value is deliberately parsed as float to keep the precision of the
+ * previous implementation.
+ */
+float parse_float_field(const std::string &line, std::size_t start,
+                        std::size_t end) {
+  float value = 0;
+  const char *first = line.data() + start;
+  const char *last = line.data() + end;
+  auto result = std::from_chars(first, last, value);
+  if (result.ec == std::errc() && result.ptr == last) return value;
+  return std::stof(line.substr(start, end - start));
+}
+} // namespace
+
 std::vector<Trajectory> ITrajectoryReader::read_next_N_trajectories(int N) {
   std::vector<Trajectory> trajectories;
   int i = 0;
@@ -271,25 +308,32 @@ Trajectory CSVPointReader::read_next_trajectory() {
       line = prev_line;
       prev_line.clear();
     }
-    std::stringstream ss(line);
     int id = 0;
     double x = 0, y = 0;
     double timestamp = 0;
     int index = 0;
-    while (std::getline(ss, intermediate, delim)) {
+    // Split the line on delim and convert the columns of interest. This is
+    // equivalent to tokenising with std::getline and converting the tokens
+    // with std::stoi/std::stof, without the overhead of a string stream.
+    std::size_t start = 0;
+    const std::size_t line_size = line.size();
+    while (start < line_size) {
+      std::size_t end = line.find(delim, start);
+      if (end == std::string::npos) end = line_size;
       if (index == id_idx) {
-        id = std::stoi(intermediate);
+        id = parse_int_field(line, start, end);
       }
       if (index == x_idx) {
-        x = std::stof(intermediate);
+        x = parse_float_field(line, start, end);
       }
       if (index == y_idx) {
-        y = std::stof(intermediate);
+        y = parse_float_field(line, start, end);
       }
       if (index == timestamp_idx) {
-        timestamp = std::stof(intermediate);
+        timestamp = parse_float_field(line, start, end);
       }
       ++index;
+      start = end + 1;
     }
     if (prev_id == id || first_observation) {
       geom.add_point(x, y);
@@ -306,19 +350,27 @@ Trajectory CSVPointReader::read_next_trajectory() {
       prev_line = line;
     }
   }
-  if (!has_next_trajectory()) {
+  if (on_same_trajectory) {
+    // The loop ended because the file is exhausted, so every point read belongs to
+    // the trajectory prev_id. When it ended because a new id showed up, trid was set
+    // above and prev_id already holds the new id (that line is buffered in prev_line).
+    // Checking the stream here instead would relabel the trajectory with the new id
+    // whenever the buffered line is the last line of the file.
     trid = prev_id;
   }
   return Trajectory{trid, geom, timestamps};
 }
 
 bool CSVPointReader::has_next_trajectory() {
-  return ifs.peek() != EOF;
+  // A line buffered by the previous call is still a pending trajectory, even
+  // when it is the last line of the file and the stream itself is at EOF.
+  return !prev_line.empty() || ifs.peek() != EOF;
 }
 
 void CSVPointReader::reset_cursor() {
   ifs.clear();
   ifs.seekg(0, std::ios::beg);
+  prev_line.clear();
   std::string line;
   std::getline(ifs, line);
 }