From 243963063cacd83c00dc79404b6bf80c957b9de6 Mon Sep 17 00:00:00 2001 From: awe Date: Mon, 21 Sep 2026 13:39:04 +0300 Subject: [PATCH] make gpr update --- .gitignore | 2 + Makefile | 28 +- README.md | 41 +- .../gpr_cpp/CMakeLists.txt | 78 + data_acq_and_processing/gpr_cpp/README.md | 87 + .../gpr_cpp/build_cmake.sh | 87 + .../gpr_cpp/examples/standalone_bp.cpp | 145 + .../gpr_cpp/include/gpr/ascan.hpp | 58 + .../gpr_cpp/include/gpr/config.hpp | 102 + .../gpr_cpp/include/gpr/device.hpp | 52 + .../gpr_cpp/include/gpr/geometry.hpp | 28 + .../gpr_cpp/include/gpr/nvrtc_kernel.hpp | 25 + .../gpr_cpp/include/gpr/sar3d.hpp | 44 + .../gpr_cpp/include/gpr/types.hpp | 25 + data_acq_and_processing/gpr_cpp/src/ascan.cpp | 198 + .../gpr_cpp/src/device.cpp | 34 + .../gpr_cpp/src/geometry.cpp | 68 + .../gpr_cpp/src/kernel_src.hpp | 72 + .../gpr_cpp/src/nvrtc_kernel.cpp | 48 + data_acq_and_processing/gpr_cpp/src/sar3d.cpp | 111 + .../third_party/pocketfft/pocketfft_hdronly.h | 3853 +++++++++++++++++ requirements.txt | 5 + 22 files changed, 5161 insertions(+), 30 deletions(-) create mode 100644 data_acq_and_processing/gpr_cpp/CMakeLists.txt create mode 100644 data_acq_and_processing/gpr_cpp/README.md create mode 100644 data_acq_and_processing/gpr_cpp/build_cmake.sh create mode 100644 data_acq_and_processing/gpr_cpp/examples/standalone_bp.cpp create mode 100644 data_acq_and_processing/gpr_cpp/include/gpr/ascan.hpp create mode 100644 data_acq_and_processing/gpr_cpp/include/gpr/config.hpp create mode 100644 data_acq_and_processing/gpr_cpp/include/gpr/device.hpp create mode 100644 data_acq_and_processing/gpr_cpp/include/gpr/geometry.hpp create mode 100644 data_acq_and_processing/gpr_cpp/include/gpr/nvrtc_kernel.hpp create mode 100644 data_acq_and_processing/gpr_cpp/include/gpr/sar3d.hpp create mode 100644 data_acq_and_processing/gpr_cpp/include/gpr/types.hpp create mode 100644 data_acq_and_processing/gpr_cpp/src/ascan.cpp create mode 100644 data_acq_and_processing/gpr_cpp/src/device.cpp create mode 100644 data_acq_and_processing/gpr_cpp/src/geometry.cpp create mode 100644 data_acq_and_processing/gpr_cpp/src/kernel_src.hpp create mode 100644 data_acq_and_processing/gpr_cpp/src/nvrtc_kernel.cpp create mode 100644 data_acq_and_processing/gpr_cpp/src/sar3d.cpp create mode 100644 data_acq_and_processing/gpr_cpp/third_party/pocketfft/pocketfft_hdronly.h diff --git a/.gitignore b/.gitignore index deb8359..7bdfab8 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,8 @@ __pycache__/ *.py[codz] *$py.class build* +# build* иначе прячет скрипт сборки вендоренного gpr_cpp (его build/ остаётся игнорируемым) +!data_acq_and_processing/gpr_cpp/build_cmake.sh .codex .DS_Store # C extensions diff --git a/Makefile b/Makefile index a627502..84e4c0d 100644 --- a/Makefile +++ b/Makefile @@ -8,15 +8,17 @@ VISA_LDFLAGS := $(LDFLAGS) -lvisa BUILD_DIR := build BIN_DIR := $(BUILD_DIR)/bin -# --- gpr_cpp: external C++/CUDA 3D SAR library (VNA_AlGo). Build it first with -# `bash $(GPR_CPP_DIR)/build_cmake.sh`; it is NOT built by this Makefile. CUDA -# headers and libnvrtc come from the pip nvidia-* packages inside the VNA_AlGo -# venv (same detection as build_cmake.sh); when those are absent (e.g. Jetson, -# where CUDA comes with JetPack) the system toolkit at /usr/local/cuda is used. -# libcuda is always the system driver. -GPR_CPP_DIR ?= $(HOME)/Documents/VNA_AlGo/gpr_cpp +# --- gpr_cpp: vendored C++/CUDA 3D SAR library (originally from VNA_AlGo). +# Built automatically by the $(GPR_CPP_LIB) rule below via its own CMake +# (build_cmake.sh). CUDA headers and libnvrtc come from the pip nvidia-* +# packages inside this repo's venv (same detection as build_cmake.sh); when +# those are absent (e.g. Jetson, where CUDA comes with JetPack) the system +# toolkit at /usr/local/cuda is used. libcuda is always the system driver. +GPR_CPP_DIR ?= data_acq_and_processing/gpr_cpp GPR_CPP_LIB ?= $(GPR_CPP_DIR)/build/libgpr_core_cpp.a -GPR_VENV ?= $(GPR_CPP_DIR)/../.venv +# Absolute so the derived -rpath (libnvrtc dir) survives running the binary +# from outside the repo root. +GPR_VENV ?= $(CURDIR)/.venv GPR_CUDA_INC_DIR ?= $(dir $(firstword $(wildcard $(GPR_VENV)/lib/python*/site-packages/nvidia/*/include/nvrtc.h))) # cuda.h can live in a DIFFERENT pip package than nvrtc.h (Jetson: cuda_runtime # vs cuda_nvrtc); the unified cu13 package on dev machines has both in one dir. @@ -181,6 +183,14 @@ $(BIN_DIR)/data_processor: $(DATA_PROCESSOR_OBJS) $(GPR_CPP_LIB) @mkdir -p $(BIN_DIR) $(CXX) $(DATA_PROCESSOR_OBJS) -o $@ $(DATA_PROCESSOR_LDFLAGS) +# Vendored gpr_cpp static lib: delegate to its CMake build (also produces the +# standalone_bp example used by sar3d_parity cross-checks). +GPR_CPP_SRC_FILES := $(wildcard $(GPR_CPP_DIR)/src/*) $(wildcard $(GPR_CPP_DIR)/include/gpr/*) \ + $(wildcard $(GPR_CPP_DIR)/examples/*) $(GPR_CPP_DIR)/CMakeLists.txt $(GPR_CPP_DIR)/build_cmake.sh + +$(GPR_CPP_LIB): $(GPR_CPP_SRC_FILES) + GPR_BUILD_PYTHON=OFF VENV=$(abspath $(GPR_VENV)) bash $(GPR_CPP_DIR)/build_cmake.sh + # The collector compiles into its own build/kamil/ tree with collector-only # includes (no VISA), so the shared driver/config sources it reuses never collide # with the orchestrator's objects of the same name. Basenames are unique, so a @@ -202,4 +212,4 @@ $(BUILD_DIR)/%.o: %.cpp -include $(DEPFILES) clean: - rm -rf $(BUILD_DIR) + rm -rf $(BUILD_DIR) $(GPR_CPP_DIR)/build diff --git a/README.md b/README.md index 5e3b5cb..aef7fb1 100644 --- a/README.md +++ b/README.md @@ -60,7 +60,7 @@ controls plus a pixel-identical live plot through an embedded web UI. | Path | Contents | | --- | --- | -| `data_acq_and_processing/` | C++ pipeline: `sweep_orchestrator/`, `preprocessing/`, `processing/` (including `locator/`), shared `common_cpp/` (config, IPC rings), and vendored `third_party/` (Eigen, nlohmann/json). | +| `data_acq_and_processing/` | C++ pipeline: `sweep_orchestrator/`, `preprocessing/`, `processing/` (including `locator/`), shared `common_cpp/` (config, IPC rings), the vendored `gpr_cpp/` C++/CUDA 3D SAR library, and vendored `third_party/` (Eigen, nlohmann/json). | | `python_app/` | Python application: `gui/`, `orchestration/` (process supervisor, shared-memory readers, config), `hardware_full/` (device drivers and services), `webui/`, `storage/`, `models/`, `workflows/`, `scripts/` (producers and tools), and `tests/`. | | `run_config_examples/` | Ready-made configuration profiles, one per hardware setup. | | `deploy/` | systemd daemon installer for Raspberry Pi appliances. | @@ -75,12 +75,12 @@ controls plus a pixel-identical live plot through an embedded web UI. - A C++20 compiler and `make` — Eigen and nlohmann/json are vendored, so no extra C++ dependencies are needed. - **Python 3.10 or newer**, in a virtual environment. - For the `sar3d` processor (3D SAR): an NVIDIA GPU with the driver installed - (`libcuda`) and the external `gpr_cpp` library from the VNA_AlGo project. It - must be built **before** `make` (see below); its location is configurable via - the `GPR_CPP_DIR` make variable (default `$HOME/Documents/VNA_AlGo/gpr_cpp`). - CUDA headers/libnvrtc are taken from the pip `nvidia-*` packages in the - VNA_AlGo venv when present, otherwise from the system toolkit at - `/usr/local/cuda` (Jetson/JetPack, apt `cuda-toolkit`). + (`libcuda`) and `cmake` (the vendored `gpr_cpp` library at + `data_acq_and_processing/gpr_cpp/` builds via its own CMake, invoked + automatically by `make`). CUDA headers/libnvrtc are taken from the pip + `nvidia-*` packages in this repo's venv (installed by `requirements.txt` on + x86_64), otherwise from the system toolkit at `/usr/local/cuda` + (Jetson/JetPack, apt `cuda-toolkit`). ## Build and run @@ -90,19 +90,22 @@ From the repository root: python3 -m venv .venv .venv/bin/python -m pip install --upgrade pip .venv/bin/python -m pip install -r requirements.txt -# 3D SAR library (external, CMake + NVRTC; CUDA headers come from the pip -# nvidia-* packages in the VNA_AlGo venv, or /usr/local/cuda as a fallback): -bash ~/Documents/VNA_AlGo/gpr_cpp/build_cmake.sh make # build the C++ binaries into build/bin/ ``` -To verify the sar3d processor against the VNA_AlGo reference on the recorded +`make` also builds the vendored 3D SAR library +(`data_acq_and_processing/gpr_cpp/`, CMake + NVRTC — needs `cmake`; CUDA +headers come from the pip `nvidia-*` packages in `.venv`, or `/usr/local/cuda` +as a fallback). + +To verify the sar3d processor against the gpr_cpp reference on the recorded dataset, build and run the on-demand parity tool: ```bash make sar3d_parity ./build/bin/sar3d_parity -# compare the printed peak/argmax with gpr_cpp's ./build/standalone_bp on the same arguments +# compare the printed peak/argmax with +# ./data_acq_and_processing/gpr_cpp/build/standalone_bp on the same arguments ``` Run the desktop GUI: @@ -127,18 +130,14 @@ that pin the stock `./start.sh` provisioning works; the extra steps are only for the CUDA-backed `sar3d` processor: ```bash -# 1. gpr_cpp against the JetPack system CUDA (no pip nvidia-* packages needed): -git clone ~/Documents/VNA_AlGo -GPR_BUILD_PYTHON=OFF bash ~/Documents/VNA_AlGo/gpr_cpp/build_cmake.sh -# build_cmake.sh falls back to /usr/local/cuda automatically; override with -# CUDA_INC=... CUDA_LIB=... if the toolkit lives elsewhere. - -# 2. copy the replay dataset (for dataset-replay profiles) and point +# 1. copy the replay dataset (for dataset-replay profiles) and point # radar.replay.dataset_dir in the profile at the local copy — the path in # run_config_examples/*.json is machine-specific. -# 3. normal launch; make picks up /usr/local/cuda when the VNA_AlGo venv has -# no pip nvidia-* packages: +# 2. normal launch; the vendored gpr_cpp builds against the JetPack system CUDA +# (requirements.txt installs no pip nvidia-* packages on aarch64, and +# build_cmake.sh falls back to /usr/local/cuda automatically; override with +# CUDA_INC=... CUDA_LIB=... if the toolkit lives elsewhere): ./start.sh --profile run_config_examples/run_config_replay_vna.example.json ``` diff --git a/data_acq_and_processing/gpr_cpp/CMakeLists.txt b/data_acq_and_processing/gpr_cpp/CMakeLists.txt new file mode 100644 index 0000000..59afccb --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/CMakeLists.txt @@ -0,0 +1,78 @@ +# 3.22 = apt cmake на Ubuntu 22.04 (Jetson/JetPack); фич новее здесь нет. +cmake_minimum_required(VERSION 3.22) +project(gpr_cpp LANGUAGES CXX) # CXX-only: CUDA компилируется в рантайме через NVRTC + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif() +set(CMAKE_CXX_FLAGS_RELEASE "-O3") + +# bindings/pymodule.cpp не вендорен в radar_system — включение потребует вернуть его из VNA_AlGo +option(GPR_BUILD_PYTHON "Build pybind11 module" OFF) +option(GPR_BUILD_TESTS "Build C++ parity tests" OFF) +option(GPR_BUILD_EXAMPLES "Build standalone C++ example" ON) + +# Python окружение (venv) — источник pybind11 и pip CUDA-библиотек (nvrtc, cudart). +if(NOT DEFINED PYTHON_EXECUTABLE) + set(PYTHON_EXECUTABLE python3) +endif() + +# CUDA include/lib из pip-пакетов nvidia-*; задаются -DGPR_CUDA_INC / -DGPR_CUDA_LIB +# (детект делает build_cmake.sh). Заголовки: nvrtc.h, cuda.h (Driver API). +# GPR_CUDA_INC может быть СПИСКОМ каталогов (через ';'): в раздельных pip-пакетах +# (Jetson: cuda_nvrtc + cuda_runtime) nvrtc.h и cuda.h лежат в разных include/. +if(NOT DEFINED GPR_CUDA_INC OR NOT DEFINED GPR_CUDA_LIB) + message(FATAL_ERROR "Задайте -DGPR_CUDA_INC= и -DGPR_CUDA_LIB=") +endif() + +# libcuda (драйвер) — системный +find_library(CUDA_DRIVER_LIB cuda PATHS /usr/lib/x86_64-linux-gnu /usr/lib/aarch64-linux-gnu + /usr/lib REQUIRED) +# libnvrtc из pip +find_library(NVRTC_LIB NAMES nvrtc libnvrtc.so.13 libnvrtc.so.12 nvrtc.so + PATHS ${GPR_CUDA_LIB} REQUIRED NO_DEFAULT_PATH) + +# PocketFFT (header-only, vendored) +add_library(pocketfft INTERFACE) +target_include_directories(pocketfft INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/third_party/pocketfft) + +# --- ядро библиотеки (usable без Python) --- +add_library(gpr_core_cpp STATIC + src/geometry.cpp + src/ascan.cpp + src/device.cpp + src/nvrtc_kernel.cpp + src/sar3d.cpp) +target_include_directories(gpr_core_cpp PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR}/include + ${GPR_CUDA_INC}) +target_link_libraries(gpr_core_cpp PUBLIC pocketfft ${NVRTC_LIB} ${CUDA_DRIVER_LIB}) +# rpath на каталог pip libnvrtc — чтобы бинарник/модуль грузились без LD_LIBRARY_PATH +set_target_properties(gpr_core_cpp PROPERTIES POSITION_INDEPENDENT_CODE ON + BUILD_RPATH "${GPR_CUDA_LIB}" INSTALL_RPATH "${GPR_CUDA_LIB}") + +if(GPR_BUILD_PYTHON) + execute_process(COMMAND ${PYTHON_EXECUTABLE} -m pybind11 --cmakedir + OUTPUT_VARIABLE PB11_DIR OUTPUT_STRIP_TRAILING_WHITESPACE RESULT_VARIABLE PB11_RC) + if(PB11_RC EQUAL 0) + find_package(pybind11 CONFIG REQUIRED HINTS ${PB11_DIR}) + else() + find_package(pybind11 CONFIG REQUIRED) + endif() + pybind11_add_module(_gpr_cpp bindings/pymodule.cpp) + target_link_libraries(_gpr_cpp PRIVATE gpr_core_cpp) + set_target_properties(_gpr_cpp PROPERTIES BUILD_RPATH "${GPR_CUDA_LIB}" INSTALL_RPATH "${GPR_CUDA_LIB}") +endif() + +if(GPR_BUILD_TESTS) + add_executable(gpr_test_parity tests/cpp/test_parity.cpp) + target_link_libraries(gpr_test_parity PRIVATE gpr_core_cpp) +endif() + +if(GPR_BUILD_EXAMPLES) + add_executable(standalone_bp examples/standalone_bp.cpp) + target_link_libraries(standalone_bp PRIVATE gpr_core_cpp) + set_target_properties(standalone_bp PROPERTIES BUILD_RPATH "${GPR_CUDA_LIB}" INSTALL_RPATH "${GPR_CUDA_LIB}") +endif() diff --git a/data_acq_and_processing/gpr_cpp/README.md b/data_acq_and_processing/gpr_cpp/README.md new file mode 100644 index 0000000..681b23f --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/README.md @@ -0,0 +1,87 @@ +# gpr_cpp — C++/CUDA счёт-ядро MIMO GPR 3D SAR + +> Вендорено из репозитория VNA_AlGo (`gpr_cpp/`, коммит `544dc59`); python-биндинги +> (`bindings/`, pybind11) при переносе не взяты — сборка здесь только C++ +> (`GPR_BUILD_PYTHON=OFF`). Собирается из корневого Makefile radar_system +> (цель `$(GPR_CPP_LIB)` вызывает `build_cmake.sh`). + +Порт вычислительного ядра из Python/CuPy (`gpr_core`) в самостоятельную +C++/CUDA библиотеку для интеграции в C++ проект. +Алгоритм — см. `ALGORITHM.md` в репозитории VNA_AlGo. + +## Что внутри (граница C++ ↔ Python) + +**В C++ (`libgpr_core_cpp`):** геометрия, препроцессинг A-сканов (bg-subtract, +фазовая коррекция движения, k0-вставка + Ханн + IFFT через PocketFFT, нормировка), +3D SAR апертурный объём (одно CUDA-ядро на всю апертуру). Возвращает **сырые** +`acc`(complex64)/`cnt`(float32). + +**Остаётся в Python:** загрузка `.npy`, финализация объёма (abs+gate+gaussian+ +нормировка), CLEAN-детекция, трек, оценка скорости, plotly. + +## Ключевое решение: CUDA без nvcc (NVRTC + Driver API) + +На целевых машинах (x86+T600, Jetson) нет полного CUDA toolkit — только runtime + +**NVRTC** (как у cupy). Поэтому ядро **не компилируется nvcc'ом**: его исходник +встроен строкой ([src/kernel_src.hpp](src/kernel_src.hpp)), компилируется в рантайме +через NVRTC в **cubin под точную арх устройства** (запрос compute capability → +`sm_XY`, без JIT PTX), грузится и запускается через **CUDA Driver API**. Хост-часть +собирается обычным g++. Плюс: без toolkit и без привязки к SM — один код для +x86 (sm_75) и Jetson (sm_87). + +## Сборка + +```bash +bash build_cmake.sh # авто-детект pip CUDA-путей; собирает всё +# артефакты в build/: libgpr_core_cpp.a, _gpr_cpp*.so, standalone_bp +``` +Требуется: g++ ≥13, cmake ≥3.24, pybind11 (`pip install pybind11`), pip-пакеты +`nvidia-*` (nvrtc+cudart, идут с cupy), системный `libcuda`. rpath на каталог +libnvrtc прописан — бинарники/модуль грузятся без `LD_LIBRARY_PATH`. + +## Публичный C++ API ([include/gpr/](include/gpr/)) + +```cpp +#include "gpr/config.hpp" // GeometryConfig/BPConfig/MotionConfig/SARConfig (дефолты = config.py) +#include "gpr/ascan.hpp" // preprocess_snapshot(...) -> AScanSet +#include "gpr/sar3d.hpp" // SarEngine::aperture_volume(frames, grid) -> VolumeRaw{acc,cnt} +#include "gpr/geometry.hpp" // frame_offsets, frame_weights, build_grid_scene, point_to_scene + +gpr::SarEngine engine(geom, bp, mot, sar); // компилит ядро под устройство +auto vol = engine.aperture_volume(frames, grid); // один запуск ядра на апертуру +// vol.acc / vol.cnt — сырые; финализацию (gaussian/normalize) делает вызывающий +``` +Пример без Python: [examples/standalone_bp.cpp](examples/standalone_bp.cpp) — грузит +S21 из `.npy`, строит объём, печатает максимум. +```bash +./build/standalone_bp 0038_id39_ns11361438137486 +``` + +## Python-модуль + +```python +import _gpr_cpp # из gpr_cpp/build +eng = _gpr_cpp.SarEngine(_gpr_cpp.GeometryConfig(), _gpr_cpp.BPConfig(), + _gpr_cpp.MotionConfig(), _gpr_cpp.SARConfig()) +a = _gpr_cpp.preprocess_snapshot(s21, freq, background, mot, geom, bp, True) # -> AScanSet +acc, cnt = eng.aperture_volume([(a, dx, dy, dz, w), ...], gx, gy, gz) # numpy +``` + +## Паритет (проверено на x86 T600) + +| Этап | Порог | Факт | +|---|---|---| +| config дефолты vs dataclass | точно | OK | +| A-сканы raw H vs numpy (float64) | <1e-9 | 1.7e-12 | +| A-сканы норм. H (FFT-ULP через медиану) | <1e-4 | 2.5e-6 | +| **3D объём C++ ядро vs cupy** (float32) | <5e-3 | **3e-7, argmax идентичен** | + +Тесты: [../tests/test_parity_cpp.py](../tests/test_parity_cpp.py) (C++ vs cupy). + +## Статус + +- [x] x86 (T600, sm_75, CUDA 13): собрано, паритет, standalone работает. +- [ ] Jetson (Orin NX, sm_87, CUDA 12.6): порт тривиален (NVRTC сам берёт sm_87); + собрать `build_cmake.sh` на Jetson с его venv. +- Опц.: ветка `impl='cpp'` в `gpr_core/sar3d.py` (адаптер `gpr_cpp_backend.py`), + кэш device-A-сканов (сейчас каждый вызов грузит стек на GPU). diff --git a/data_acq_and_processing/gpr_cpp/build_cmake.sh b/data_acq_and_processing/gpr_cpp/build_cmake.sh new file mode 100644 index 0000000..1af3377 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/build_cmake.sh @@ -0,0 +1,87 @@ +#!/usr/bin/env bash +# CMake-сборка libgpr_core_cpp + standalone (без nvcc; NVRTC/driver). +# +# Источник CUDA (nvrtc.h / libnvrtc), по убыванию приоритета: +# 1. env CUDA_INC / CUDA_LIB — заданы явно; +# 2. pip-пакеты nvidia-* в venv ($VENV, по умолчанию ../../.venv — venv +# radar_system) — дев-машина; +# 3. системный toolkit /usr/local/cuda — Jetson/JetPack, apt cuda-toolkit. +# +# GPR_BUILD_PYTHON=ON вернул бы pybind-модуль, но bindings/ не вендорен в +# radar_system (см. README) — по умолчанию OFF. +set -eu +HERE="$(cd "$(dirname "$0")" && pwd)" +VENV="${VENV:-$HERE/../../.venv}" +PY="$VENV/bin/python" +GPR_BUILD_PYTHON="${GPR_BUILD_PYTHON:-OFF}" + +CUDA_INC="${CUDA_INC:-}" +CUDA_LIB="${CUDA_LIB:-}" + +if [ -z "$CUDA_INC" ] || [ -z "$CUDA_LIB" ]; then + if [ -x "$PY" ]; then + # Печатается только когда найдены nvrtc.h, cuda.h И библиотека — частичная + # детекция не может подсунуть мусорный путь. nvrtc.h и cuda.h могут лежать + # в РАЗНЫХ каталогах (раздельные pip-пакеты cuda_nvrtc/cuda_runtime на + # Jetson) — тогда CUDA_INC становится ';'-списком для CMake. + PIP_CUDA=$("$PY" - <<'P' 2>/dev/null || true +import glob, os, site +sp = site.getsitepackages()[0] +nvrtc_hdr = glob.glob(os.path.join(sp, 'nvidia', '**', 'include', 'nvrtc.h'), recursive=True) +cuda_hdr = glob.glob(os.path.join(sp, 'nvidia', '**', 'include', 'cuda.h'), recursive=True) +lib = sorted(glob.glob(os.path.join(sp, 'nvidia', '**', 'lib', 'libnvrtc.so*'), recursive=True)) +if nvrtc_hdr and cuda_hdr and lib: + inc_dirs = [] + for hdr in (nvrtc_hdr[0], cuda_hdr[0]): + d = os.path.dirname(hdr) + if d not in inc_dirs: + inc_dirs.append(d) + print(';'.join(inc_dirs)) + print(os.path.dirname(lib[0])) +P +) + if [ -n "$PIP_CUDA" ]; then + CUDA_INC=$(printf '%s\n' "$PIP_CUDA" | sed -n 1p) + CUDA_LIB=$(printf '%s\n' "$PIP_CUDA" | sed -n 2p) + fi + fi +fi + +if [ -z "$CUDA_INC" ] || [ -z "$CUDA_LIB" ]; then + if [ -e /usr/local/cuda/include/nvrtc.h ]; then + CUDA_INC=/usr/local/cuda/include + CUDA_LIB=/usr/local/cuda/lib64 + echo "[build_cmake] pip nvidia-* не найдены; использую системный CUDA: $CUDA_INC" + else + echo "[build_cmake] CUDA не найден: нет ни pip nvidia-* в $VENV, ни /usr/local/cuda." >&2 + echo "[build_cmake] Задайте CUDA_INC= и CUDA_LIB=." >&2 + exit 1 + fi +fi + +PYTHON_ARG=() +if [ "$GPR_BUILD_PYTHON" = "ON" ] && [ -x "$PY" ]; then + PYTHON_ARG=(-DPYTHON_EXECUTABLE="$PY") +fi + +# cmake: сначала из venv (pip install cmake — свежий без sudo), потом системный. +CMAKE="$VENV/bin/cmake" +if [ ! -x "$CMAKE" ]; then + CMAKE=$(command -v cmake || true) +fi +if [ -z "$CMAKE" ]; then + echo "[build_cmake] cmake не найден. Поставьте его в venv: \"$VENV/bin/pip\" install cmake" >&2 + echo "[build_cmake] (или системно: sudo apt install cmake)" >&2 + exit 1 +fi + +"$CMAKE" -S "$HERE" -B "$HERE/build" \ + ${PYTHON_ARG[@]+"${PYTHON_ARG[@]}"} \ + -DGPR_BUILD_PYTHON="$GPR_BUILD_PYTHON" \ + -DGPR_CUDA_INC="$CUDA_INC" \ + -DGPR_CUDA_LIB="$CUDA_LIB" \ + -DCMAKE_BUILD_TYPE=Release "$@" 2>&1 | tail -8 +"$CMAKE" --build "$HERE/build" -j4 2>&1 | tail -15 +echo "--- артефакты ---" +ls -la "$HERE/build/libgpr_core_cpp.a" "$HERE/build/standalone_bp" "$HERE"/build/_gpr_cpp*.so 2>/dev/null | awk '{print $5,$9}' +echo "CUDA_LIB=$CUDA_LIB (для LD_LIBRARY_PATH при запуске)" diff --git a/data_acq_and_processing/gpr_cpp/examples/standalone_bp.cpp b/data_acq_and_processing/gpr_cpp/examples/standalone_bp.cpp new file mode 100644 index 0000000..c263c4c --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/examples/standalone_bp.cpp @@ -0,0 +1,145 @@ +// Автономный пример: C++ без Python. Грузит S21 снапшота из .npy, +// строит A-сканы и 3D SAR-объём апертуры, печатает максимум |acc/cnt|. +// +// Использование: +// standalone_bp +// где preprocessed_dir содержит папки снапшотов с i{rx}_o{tx}_s21.npy и _freq.npy. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "gpr/ascan.hpp" +#include "gpr/config.hpp" +#include "gpr/geometry.hpp" +#include "gpr/sar3d.hpp" + +using namespace gpr; + +// --- минимальный загрузчик numpy .npy (v1.0, 1-D complex128 или float64) --- +struct NpyArray { std::vector f64; std::vector c128v; bool is_complex = false; }; + +static NpyArray load_npy(const std::string& path) { + std::ifstream in(path, std::ios::binary); + if (!in) throw std::runtime_error("не открыть " + path); + char magic[6]; in.read(magic, 6); + uint8_t ver[2]; in.read((char*)ver, 2); + uint16_t hlen; in.read((char*)&hlen, 2); + std::string hdr(hlen, '\0'); in.read(&hdr[0], hlen); + auto sp = hdr.find("'shape':"); + auto lp = hdr.find('(', sp), rp = hdr.find(')', sp); + long n = std::stol(hdr.substr(lp + 1, rp - lp - 1)); + NpyArray a; + auto has = [&](const char* s) { return hdr.find(s) != std::string::npos; }; + if (has("c8")) { // complex64 → c128 + a.is_complex = true; a.c128v.resize(n); + std::vector> tmp(n); + in.read((char*)tmp.data(), n * sizeof(std::complex)); + for (long i = 0; i < n; ++i) a.c128v[i] = c128(tmp[i].real(), tmp[i].imag()); + } else if (has("c16")) { // complex128 + a.is_complex = true; a.c128v.resize(n); + in.read((char*)a.c128v.data(), n * sizeof(c128)); + } else if (has("f4")) { // float32 → double + a.f64.resize(n); + std::vector tmp(n); in.read((char*)tmp.data(), n * sizeof(float)); + for (long i = 0; i < n; ++i) a.f64[i] = tmp[i]; + } else { // float64 + a.f64.resize(n); in.read((char*)a.f64.data(), n * sizeof(double)); + } + return a; +} + +int main(int argc, char** argv) { + if (argc < 3) { printf("usage: %s \n", argv[0]); return 1; } + std::string root = argv[1], center = argv[2]; + + GeometryConfig geom; geom.finalize(); + BPConfig bp; MotionConfig mot; SARConfig sar; + auto pairs = geom.pair_keys(); + auto tx = geom.output_idx, rx = geom.input_idx; + std::sort(tx.begin(), tx.end()); std::sort(rx.begin(), rx.end()); + + // список снапшотов, выбор апертуры вокруг center + std::vector snaps; + if (DIR* d = opendir(root.c_str())) { + for (dirent* e; (e = readdir(d));) { + std::string n = e->d_name; + if (!n.empty() && isdigit(n[0])) snaps.push_back(n); + } + closedir(d); + } + std::sort(snaps.begin(), snaps.end()); + int ci = -1; for (int i = 0; i < (int)snaps.size(); ++i) if (snaps[i] == center) ci = i; + if (ci < 0) { printf("центр %s не найден\n", center.c_str()); return 1; } + int lo = std::max(0, ci - sar.frame_radius), hi = std::min((int)snaps.size(), ci + sar.frame_radius + 1); + std::vector aperture(snaps.begin() + lo, snaps.begin() + hi); + printf("апертура: %d кадров вокруг %s\n", (int)aperture.size(), center.c_str()); + + auto load_snapshot = [&](const std::string& snap, + std::map>& s21, + std::map>& freq) { + for (int p = 0; p < (int)pairs.size(); ++p) { + int i_tx = pairs[p].first, i_rx = pairs[p].second; + int phys_tx = tx[i_tx], phys_rx = rx[i_rx]; + char fn[512]; + snprintf(fn, sizeof(fn), "%s/%s/i%d_o%d_s21.npy", root.c_str(), snap.c_str(), phys_rx, phys_tx); + s21[{i_tx, i_rx}] = load_npy(fn).c128v; + snprintf(fn, sizeof(fn), "%s/%s/i%d_o%d_freq.npy", root.c_str(), snap.c_str(), phys_rx, phys_tx); + freq[{i_tx, i_rx}] = load_npy(fn).f64; + } + }; + + auto ts = compute_pair_timestamps(mot, geom); + SarEngine engine(geom, bp, mot, sar); + + // держим AScanSet живыми + std::vector> keep; + std::vector frames; + std::vector frame_times; + for (auto& snap : aperture) { + // время из суффикса _nsNNNN + auto pos = snap.rfind("_ns"); + double t = (pos != std::string::npos) ? std::stod(snap.substr(pos + 3)) * 1e-9 : 0.0; + frame_times.push_back(t); + } + double t_ref = frame_times[ci - lo]; + auto offs = frame_offsets(frame_times, t_ref, mot); + auto wts = frame_weights(offs, sar.edge_frame_weight); + + for (int f = 0; f < (int)aperture.size(); ++f) { + std::map> s21; + std::map> freq; + load_snapshot(aperture[f], s21, freq); + auto a = std::make_shared( + preprocess_snapshot(s21, freq, nullptr, mot, geom, bp, ts, true)); + keep.push_back(a); + frames.push_back(ApertureFrame{a.get(), offs[f], wts[f]}); + } + + // центр сетки: под центральным снапшотом на глубине 5.6 м (сцена) + auto center_scene = point_to_scene(0.0, 0.0, 5.6, mot, sar); + Grid grid = build_grid_scene(center_scene, sar); + printf("сетка %dx%dx%d, счёт объёма...\n", grid.NX(), grid.NY(), grid.NZ()); + + VolumeRaw vol = engine.aperture_volume(frames, grid); + + // финализация: |acc/cnt|, максимум + double vmax = 0; size_t amax = 0; + for (size_t i = 0; i < vol.acc.size(); ++i) { + double v = std::abs(vol.acc[i]) / (vol.cnt[i] + 1e-12); + if (v > vmax) { vmax = v; amax = i; } + } + int NX = vol.NX, NZ = vol.NZ; + int ix = amax % NX, iz = (amax / NX) % NZ, iy = amax / ((size_t)NX * NZ); + printf("МАКСИМУМ объёма: |acc/cnt|=%.4e в вокселе (ix=%d,iy=%d,iz=%d)\n", vmax, ix, iy, iz); + printf(" X=%.1f см, Y=%.1f см, Z=%.1f см\n", + grid.x[ix] * 100, grid.y[iy] * 100, grid.z[iz] * 100); + printf("OK: автономный C++ путь отработал без Python.\n"); + return 0; +} diff --git a/data_acq_and_processing/gpr_cpp/include/gpr/ascan.hpp b/data_acq_and_processing/gpr_cpp/include/gpr/ascan.hpp new file mode 100644 index 0000000..e15d5f5 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/include/gpr/ascan.hpp @@ -0,0 +1,58 @@ +// Препроцессинг S21 → комплексные A-сканы (host, float64). Порт preprocess.py. +#pragma once +#include +#include +#include +#include + +#include "gpr/config.hpp" +#include "gpr/types.hpp" + +namespace gpr { + +using PairKey = std::pair; + +// Тайминги пары (нужны slot/n_slots/t_center для фазовой коррекции). +struct PairTiming { + int tx_interleave_slot = 0; + int tx_interleave_slots = 1; + double t_center_s = 0.0; +}; +using PairTimestamps = std::map; + +PairTimestamps compute_pair_timestamps(const MotionConfig&, const GeometryConfig&); + +// A-сканы одного снапшота, общая ось времени у всех пар (assert). +struct AScanSet { + int n_pairs = 0; + int n_fft = 0; + double t0 = 0.0, dt = 0.0; // общая ось: t_k = t0 + k*dt + std::vector pairs; // порядок как в geom.pair_keys() + std::vector H; // [n_pairs * n_fft] float64 + const c128* pair_ptr(int p) const { return H.data() + (size_t)p * n_fft; } +}; + +// Порт compute_pair_ascans: bg-subtract → phase correction → k0-IFFT → норма. +// s21/freq/background — по ключам (i_tx,i_rx); background может быть nullptr. +AScanSet preprocess_snapshot( + const std::map>& s21, + const std::map>& freq, + const std::map>* background, + const MotionConfig&, const GeometryConfig&, const BPConfig&, + const PairTimestamps&, bool normalize = true); + +// Отдельные шаги (для тестов/переиспользования): +std::vector apply_intra_sweep_phase_correction( + const std::vector& s21, const std::vector& freq_full, + const PairTiming&, const MotionConfig&, double wave_speed, + double f_start, double f_stop, int n_slots_default); + +// Возвращает t (ось), h (комплексный A-скан). n_fft через out-параметр. +void compute_ascan_bp(const std::vector& s21, const std::vector& freq, + double f_start, double f_stop, bool window, int oversample, + std::vector& t_out, std::vector& h_out); + +double pair_norm_scale(const std::vector& h, const std::vector& t, + double wave_speed, const BPConfig&); + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/include/gpr/config.hpp b/data_acq_and_processing/gpr_cpp/include/gpr/config.hpp new file mode 100644 index 0000000..ca16e6d --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/include/gpr/config.hpp @@ -0,0 +1,102 @@ +// Конфигурация счёт-ядра. Дефолты 1:1 повторяют gpr_core/config.py. +#pragma once +#include +#include +#include +#include +#include +#include + +namespace gpr { + +struct GeometryConfig { + std::vector input_idx{0, 1, 2, 3}; + std::vector output_idx{0, 1}; + // (x,y,z) метры; порядок по физическим индексам как в config.py + std::vector> tx_positions{{0.785, 0.0, 0.0}, {-0.785, 0.0, 0.0}}; + std::vector> rx_positions{ + {-0.42, 0.0, 0.0}, {0.48, 0.0, 0.0}, {0.18, 0.0, 0.0}, {-0.12, 0.0, 0.0}}; + double eps_r = 1.0; + double bp_plane_y = 0.0; + + // --- производные (заполняются finalize) --- + std::vector x_tx, y_tx, z_tx; // отсортированы по phys-индексу + std::vector x_rx, y_rx, z_rx; + int n_tx = 0, n_rx = 0; + double v = 0.0; + + void finalize() { + auto sort_xyz = [](std::vector idx, const std::vector>& pos, + std::vector& X, std::vector& Y, std::vector& Z) { + std::vector order(idx.size()); + for (size_t i = 0; i < idx.size(); ++i) order[i] = (int)i; + std::sort(order.begin(), order.end(), [&](int a, int b) { return idx[a] < idx[b]; }); + X.clear(); Y.clear(); Z.clear(); + for (int o : order) { X.push_back(pos[o][0]); Y.push_back(pos[o][1]); Z.push_back(pos[o][2]); } + }; + sort_xyz(output_idx, tx_positions, x_tx, y_tx, z_tx); + sort_xyz(input_idx, rx_positions, x_rx, y_rx, z_rx); + n_tx = (int)output_idx.size(); + n_rx = (int)input_idx.size(); + v = 3e8 / std::sqrt(eps_r); + } + + std::vector> pair_keys() const { + std::vector> p; + for (int i = 0; i < n_tx; ++i) + for (int j = 0; j < n_rx; ++j) p.emplace_back(i, j); + return p; + } +}; + +struct BPConfig { + double f_start = 20 * 1e8; // 2 ГГц + double f_stop = 6.0e9; + double min_depth = 3.2; + double max_depth = 10.0; + double grid_x_pad = 2.0; + double z_min = 0.1; + int nx = 300; + int nz = 300; + int oversample = 8; + bool window = true; + bool pair_normalize = true; + double pair_norm_percentile = 50.0; + double pair_norm_depth_min = 3.2; + double pair_norm_depth_max = 10.0; + double pair_norm_eps = 1e-15; + double comp_range_power = 0.25; + double comp_weight_max = 5.0; + double comp_ref_depth = 5.0; + double smooth_sigma = 1.5; + bool bg_subtract = true; + double cf_eps = 1e-12; +}; + +struct MotionConfig { + double speed_m_s = 0.854; + double look_angle_deg = 6.0; + double pitch_deg = 6.0; + double tx_sweep_time_s = 0.079; + double tx_switch_time_s = 1e-5; + // frequency_timing_mode всегда interleaved_tx_by_frequency + double direction_sign = +1.0; + bool apply_freq_phase_correction = true; + double motion_phase_sign = -1.0; + + double frame_time_s(int n_tx) const { + return n_tx * tx_sweep_time_s + std::max(0, n_tx - 1) * tx_switch_time_s; + } +}; + +struct SARConfig { + int frame_radius = 6; + double edge_frame_weight = 0.75; + std::array half_width_m{1.8, 1.0, 1.8}; + std::array grid{150, 120, 160}; // NX, NY, NZ + bool coord_mode_ground = true; // 'ground' vs 'radar_local' + bool use_compensation = true; + bool normalize_each_frame = true; +}; + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/include/gpr/device.hpp b/data_acq_and_processing/gpr_cpp/include/gpr/device.hpp new file mode 100644 index 0000000..f3ce243 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/include/gpr/device.hpp @@ -0,0 +1,52 @@ +// CUDA Driver API RAII: контекст, поток, буфер устройства. +#pragma once +#include + +#include +#include +#include + +namespace gpr { + +void cu_check(CUresult r, const char* what); + +// Инициализация драйвера + primary context на device 0. Возвращает (major,minor) CC. +struct CudaContext { + CUdevice dev = 0; + CUcontext ctx = nullptr; + int cc_major = 0, cc_minor = 0; + CudaContext(); + ~CudaContext(); + std::string arch_flag() const; // "--gpu-architecture=sm_XY" +}; + +template +class DeviceBuffer { + CUdeviceptr ptr_ = 0; + size_t n_ = 0; + + public: + DeviceBuffer() = default; + explicit DeviceBuffer(size_t n) { alloc(n); } + ~DeviceBuffer() { free(); } + DeviceBuffer(const DeviceBuffer&) = delete; + DeviceBuffer& operator=(const DeviceBuffer&) = delete; + DeviceBuffer(DeviceBuffer&& o) noexcept : ptr_(o.ptr_), n_(o.n_) { o.ptr_ = 0; o.n_ = 0; } + + void alloc(size_t n) { + free(); + if (n) cu_check(cuMemAlloc(&ptr_, n * sizeof(T)), "cuMemAlloc"); + n_ = n; + } + void free() { + if (ptr_) cuMemFree(ptr_); + ptr_ = 0; n_ = 0; + } + void upload(const T* host, size_t n) { cu_check(cuMemcpyHtoD(ptr_, host, n * sizeof(T)), "HtoD"); } + void download(T* host, size_t n) const { cu_check(cuMemcpyDtoH(host, ptr_, n * sizeof(T)), "DtoH"); } + CUdeviceptr ptr() const { return ptr_; } + CUdeviceptr* addr() { return &ptr_; } + size_t size() const { return n_; } +}; + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/include/gpr/geometry.hpp b/data_acq_and_processing/gpr_cpp/include/gpr/geometry.hpp new file mode 100644 index 0000000..53c4a3f --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/include/gpr/geometry.hpp @@ -0,0 +1,28 @@ +// Геометрия: опорное расхождение пары, смещения/веса кадров, построение сетки. +#pragma once +#include + +#include "gpr/config.hpp" +#include "gpr/types.hpp" + +namespace gpr { + +// 1/geo_ref пары (порт pair_geo_ref, скаляр float64). +double pair_geo_ref(const GeometryConfig&, const BPConfig&, int i_tx, int i_rx, double yy); + +// Смещения кадров относительно опорного (порт sar_frame_offsets). +// frame_times[k] — время кадра (сек, от ns-таймстемпа); t_ref — время опорного. +std::vector frame_offsets(const std::vector& frame_times, double t_ref, + const MotionConfig&); + +// Косинусные веса апертуры (порт sar_frame_weights); dt берётся относительно центра. +std::vector frame_weights(const std::vector& offsets, double edge_weight); + +// Сетка вокруг точки сцены (порт sar_build_grid_scene). +Grid build_grid_scene(const std::array& center_scene, const SARConfig&); + +// point_to_scene: оси радара → сцена (поворот на pitch, если ground). +std::array point_to_scene(double x, double y, double z, const MotionConfig&, + const SARConfig&); + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/include/gpr/nvrtc_kernel.hpp b/data_acq_and_processing/gpr_cpp/include/gpr/nvrtc_kernel.hpp new file mode 100644 index 0000000..de194af --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/include/gpr/nvrtc_kernel.hpp @@ -0,0 +1,25 @@ +// Компиляция строки CUDA-ядра через NVRTC → cubin под арх устройства → CUfunction. +#pragma once +#include + +#include + +namespace gpr { + +class NvrtcKernel { + CUmodule module_ = nullptr; + CUfunction func_ = nullptr; + + public: + NvrtcKernel() = default; + ~NvrtcKernel(); + NvrtcKernel(const NvrtcKernel&) = delete; + NvrtcKernel& operator=(const NvrtcKernel&) = delete; + + // Компилирует src в cubin под arch_flag ("--gpu-architecture=sm_XY"), + // загружает модуль и берёт функцию func_name. + void compile(const char* src, const std::string& arch_flag, const char* func_name); + CUfunction func() const { return func_; } +}; + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/include/gpr/sar3d.hpp b/data_acq_and_processing/gpr_cpp/include/gpr/sar3d.hpp new file mode 100644 index 0000000..31850f3 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/include/gpr/sar3d.hpp @@ -0,0 +1,44 @@ +// 3D SAR: сборка взвешенного объёма апертуры за один запуск ядра. +#pragma once +#include +#include + +#include "gpr/ascan.hpp" +#include "gpr/config.hpp" +#include "gpr/device.hpp" +#include "gpr/geometry.hpp" +#include "gpr/nvrtc_kernel.hpp" +#include "gpr/types.hpp" + +namespace gpr { + +// Один кадр апертуры: A-сканы + смещение + вес. +struct ApertureFrame { + const AScanSet* ascans; + FrameOffset offset; + double weight; +}; + +// Сырой результат объёма (Python делает finish: abs+gate+gaussian+normalize). +struct VolumeRaw { + std::vector acc; // [NY*NZ*NX] + std::vector cnt; // [NY*NZ*NX] + int NX = 0, NY = 0, NZ = 0; +}; + +class SarEngine { + public: + SarEngine(const GeometryConfig&, const BPConfig&, const MotionConfig&, const SARConfig&); + // Весь объём апертуры за один запуск ядра; возвращает сырые acc/cnt (host). + VolumeRaw aperture_volume(const std::vector& frames, const Grid&); + + private: + GeometryConfig geom_; + BPConfig bp_; + MotionConfig mot_; + SARConfig sar_; + std::unique_ptr ctx_; + std::unique_ptr kernel_; +}; + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/include/gpr/types.hpp b/data_acq_and_processing/gpr_cpp/include/gpr/types.hpp new file mode 100644 index 0000000..ab02ce3 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/include/gpr/types.hpp @@ -0,0 +1,25 @@ +// Базовые типы счёт-ядра. +#pragma once +#include +#include +#include + +namespace gpr { + +using c64 = std::complex; // complex64, layout == cuComplex (2×float) +using c128 = std::complex; // float64 для препроцессинга/IFFT + +// Ground-scene сетка (порт dict из geometry.py: x,y,z оси в метрах). +struct Grid { + std::vector x, y, z; + int NX() const { return (int)x.size(); } + int NY() const { return (int)y.size(); } + int NZ() const { return (int)z.size(); } +}; + +// Смещение антенн кадра апертуры (порт sar_frame_offsets). +struct FrameOffset { + double dx_m = 0.0, dy_m = 0.0, dz_m = 0.0, dt_s = 0.0, horizontal_m = 0.0; +}; + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/src/ascan.cpp b/data_acq_and_processing/gpr_cpp/src/ascan.cpp new file mode 100644 index 0000000..3ef5512 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/src/ascan.cpp @@ -0,0 +1,198 @@ +#include "gpr/ascan.hpp" + +#include +#include +#include +#include + +#include "pocketfft_hdronly.h" + +namespace gpr { + +PairTimestamps compute_pair_timestamps(const MotionConfig& mot, const GeometryConfig& geom) { + // interleaved_tx_by_frequency: слот Tx = порядок в отсортированном output_idx. + std::vector tx = geom.output_idx, rx = geom.input_idx; + std::sort(tx.begin(), tx.end()); + std::sort(rx.begin(), rx.end()); + int n_slots = (int)tx.size(); + double frame_center = 0.5 * n_slots * mot.tx_sweep_time_s; + PairTimestamps ts; + for (int slot = 0; slot < (int)tx.size(); ++slot) { + int i_tx = slot; // индекс в отсортированном == логический + for (int j = 0; j < (int)rx.size(); ++j) { + ts[{i_tx, j}] = PairTiming{slot, n_slots, frame_center}; + } + } + return ts; +} + +static double median_of_diff(const std::vector& f) { + std::vector d(f.size() - 1); + for (size_t i = 0; i + 1 < f.size(); ++i) d[i] = f[i + 1] - f[i]; + std::sort(d.begin(), d.end()); + size_t n = d.size(); + return (n % 2) ? d[n / 2] : 0.5 * (d[n / 2 - 1] + d[n / 2]); +} + +std::vector apply_intra_sweep_phase_correction( + const std::vector& s21, const std::vector& freq_full, + const PairTiming& info, const MotionConfig& cfg, double wave_speed, + double f_start, double f_stop, int n_slots_default) { + std::vector out = s21; + int n_full = (int)freq_full.size(); + // индексы в полосе + std::vector idx; + for (int i = 0; i < n_full; ++i) + if (freq_full[i] >= f_start && freq_full[i] <= f_stop) idx.push_back(i); + if (idx.empty()) throw std::runtime_error("phase corr: пустая полоса"); + if (n_full < 2) return out; + + int n_slots = info.tx_interleave_slots ? info.tx_interleave_slots : n_slots_default; + int slot = info.tx_interleave_slot; + double dt_base = cfg.tx_sweep_time_s / (n_full - 1); + // t_all по всему свипу; t_center — середина ПОЛНОГО свипа (как в preprocess.py) + double t_first = (double)(n_slots * 0 + slot) * dt_base; + double t_last = (double)(n_slots * (n_full - 1) + slot) * dt_base; + double t_center = 0.5 * (t_first + t_last); + double theta = cfg.look_angle_deg * M_PI / 180.0; + if (!cfg.apply_freq_phase_correction) return out; + for (int i : idx) { + double t_abs = (double)(n_slots * i + slot) * dt_base; + double dt_intra = t_abs - t_center; + double delta_range = cfg.direction_sign * cfg.speed_m_s * std::cos(theta) * dt_intra; + double dtau = 2.0 * delta_range / wave_speed; + double phi = 2.0 * M_PI * freq_full[i] * dtau; + double ang = cfg.motion_phase_sign * phi; + out[i] *= c128(std::cos(ang), std::sin(ang)); + } + return out; +} + +static int next_pow2(int x) { + int p = 1; + while (p < x) p <<= 1; + return p; +} + +void compute_ascan_bp(const std::vector& s21, const std::vector& freq, + double f_start, double f_stop, bool window, int oversample, + std::vector& t_out, std::vector& h_out) { + std::vector freq_cut; + std::vector s_cut; + for (size_t i = 0; i < freq.size(); ++i) + if (freq[i] >= f_start && freq[i] <= f_stop) { + freq_cut.push_back(freq[i]); + s_cut.push_back(s21[i]); + } + if (freq_cut.size() < 2) throw std::runtime_error("ascan: <2 точек в полосе"); + + double df = median_of_diff(freq_cut); + int n = (int)freq_cut.size(); + int k0 = (int)std::nearbyint(freq_cut[0] / df); // round-half-to-even как python round + int min_len = 2 * (k0 + n - 1); + int n_fft_base = next_pow2(min_len); + int n_fft = n_fft_base * oversample; + double dt = 1.0 / ((double)n_fft * df); + + // окно Ханна (numpy np.hanning: 0.5 - 0.5*cos(2*pi*i/(N-1))) + std::vector H(n_fft, c128(0.0, 0.0)); + for (int i = 0; i < n; ++i) { + c128 v = s_cut[i]; + if (window) { + double w = 0.5 - 0.5 * std::cos(2.0 * M_PI * i / (n - 1)); + v *= w; + } + H[k0 + i] = v; + } + + // IFFT через pocketfft (backward + 1/n_fft == numpy ifft) + h_out.assign(n_fft, c128(0.0, 0.0)); + pocketfft::shape_t shape{(size_t)n_fft}; + pocketfft::stride_t stride{(ptrdiff_t)sizeof(c128)}; + pocketfft::c2c(shape, stride, stride, {0}, pocketfft::BACKWARD, + H.data(), h_out.data(), (double)(1.0 / n_fft)); + + t_out.resize(n_fft); + for (int i = 0; i < n_fft; ++i) t_out[i] = i * dt; +} + +double pair_norm_scale(const std::vector& h, const std::vector& t, + double wave_speed, const BPConfig& bp) { + std::vector vals; + vals.reserve(h.size()); + for (size_t i = 0; i < h.size(); ++i) { + double z = t[i] * wave_speed / 2.0; + if (z >= bp.pair_norm_depth_min && z <= bp.pair_norm_depth_max) + vals.push_back(std::abs(h[i])); + } + if (vals.empty()) + for (auto& c : h) vals.push_back(std::abs(c)); + std::sort(vals.begin(), vals.end()); + // numpy.percentile linear interpolation + double q = bp.pair_norm_percentile / 100.0; + double pos = q * (vals.size() - 1); + int lo = (int)std::floor(pos); + int hi = std::min(lo + 1, (int)vals.size() - 1); + double scale = vals[lo] + (pos - lo) * (vals[hi] - vals[lo]); + if (!std::isfinite(scale) || scale <= bp.pair_norm_eps) scale = 1.0; + return scale; +} + +AScanSet preprocess_snapshot( + const std::map>& s21, + const std::map>& freq, + const std::map>* background, + const MotionConfig& mot, const GeometryConfig& geom, const BPConfig& bp, + const PairTimestamps& ts, bool normalize) { + AScanSet out; + out.pairs = geom.pair_keys(); + out.n_pairs = (int)out.pairs.size(); + + std::vector> Hs; + std::vector> Ts; + for (auto& key : out.pairs) { + auto its = s21.find(key); + if (its == s21.end()) throw std::runtime_error("preprocess: нет пары в s21"); + std::vector s = its->second; // копия + if (background) { + auto ib = background->find(key); + if (ib != background->end()) + for (size_t i = 0; i < s.size(); ++i) s[i] -= ib->second[i]; + } + const std::vector& fr = freq.at(key); + auto itt = ts.find(key); + if (itt != ts.end()) + s = apply_intra_sweep_phase_correction(s, fr, itt->second, mot, geom.v, + bp.f_start, bp.f_stop, geom.n_tx); + std::vector t; + std::vector h; + compute_ascan_bp(s, fr, bp.f_start, bp.f_stop, bp.window, bp.oversample, t, h); + Hs.push_back(std::move(h)); + Ts.push_back(std::move(t)); + } + + // общая ось времени — проверка + int n_fft = (int)Ts[0].size(); + double t0 = Ts[0][0]; + double dt = (Ts[0].back() - t0) / (n_fft - 1); + for (auto& t : Ts) + if ((int)t.size() != n_fft || std::abs(t[0] - t0) > 1e-18) + throw std::runtime_error("A-сканы: разная ось времени у пар"); + + if (normalize) { + for (int p = 0; p < out.n_pairs; ++p) { + double scale = pair_norm_scale(Hs[p], Ts[p], geom.v, bp); + for (auto& c : Hs[p]) c /= (scale + bp.pair_norm_eps); + } + } + + out.n_fft = n_fft; + out.t0 = t0; + out.dt = dt; + out.H.resize((size_t)out.n_pairs * n_fft); + for (int p = 0; p < out.n_pairs; ++p) + std::copy(Hs[p].begin(), Hs[p].end(), out.H.begin() + (size_t)p * n_fft); + return out; +} + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/src/device.cpp b/data_acq_and_processing/gpr_cpp/src/device.cpp new file mode 100644 index 0000000..c60432c --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/src/device.cpp @@ -0,0 +1,34 @@ +#include "gpr/device.hpp" + +#include + +namespace gpr { + +void cu_check(CUresult r, const char* what) { + if (r != CUDA_SUCCESS) { + const char* s = nullptr; + cuGetErrorString(r, &s); + throw std::runtime_error(std::string("CUDA error in ") + what + ": " + (s ? s : "?")); + } +} + +CudaContext::CudaContext() { + cu_check(cuInit(0), "cuInit"); + cu_check(cuDeviceGet(&dev, 0), "cuDeviceGet"); + cu_check(cuDeviceGetAttribute(&cc_major, CU_DEVICE_ATTRIBUTE_COMPUTE_CAPABILITY_MAJOR, dev), + "cc_major"); + cu_check(cuDeviceGetAttribute(&cc_minor, CU_DEVICE_ATTRIBUTE_COMPUTE_CAPABILITY_MINOR, dev), + "cc_minor"); + cu_check(cuDevicePrimaryCtxRetain(&ctx, dev), "primaryCtxRetain"); + cu_check(cuCtxSetCurrent(ctx), "ctxSetCurrent"); +} + +CudaContext::~CudaContext() { + if (ctx) cuDevicePrimaryCtxRelease(dev); +} + +std::string CudaContext::arch_flag() const { + return "--gpu-architecture=sm_" + std::to_string(cc_major) + std::to_string(cc_minor); +} + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/src/geometry.cpp b/data_acq_and_processing/gpr_cpp/src/geometry.cpp new file mode 100644 index 0000000..7576228 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/src/geometry.cpp @@ -0,0 +1,68 @@ +#include "gpr/geometry.hpp" + +#include + +namespace gpr { + +double pair_geo_ref(const GeometryConfig& g, const BPConfig& bp, int i_tx, int i_rx, double yy) { + double xc = (g.x_tx[i_tx] + g.x_rx[i_rx]) / 2.0; + double Rtx = std::sqrt(std::pow(xc - g.x_tx[i_tx], 2) + std::pow(yy - g.y_tx[i_tx], 2) + + std::pow(bp.comp_ref_depth - g.z_tx[i_tx], 2)); + double Rrx = std::sqrt(std::pow(xc - g.x_rx[i_rx], 2) + std::pow(yy - g.y_rx[i_rx], 2) + + std::pow(bp.comp_ref_depth - g.z_rx[i_rx], 2)); + return 1.0 / (Rtx * Rrx + 1e-12) + 1e-30; +} + +std::vector frame_offsets(const std::vector& frame_times, double t_ref, + const MotionConfig& cfg) { + std::vector out(frame_times.size()); + double pitch = cfg.pitch_deg * M_PI / 180.0; + for (size_t k = 0; k < frame_times.size(); ++k) { + double t = frame_times[k] - t_ref; + double h = cfg.direction_sign * cfg.speed_m_s * t; + out[k].dx_m = 0.0; + out[k].dy_m = h * std::sin(pitch); + out[k].dz_m = h * std::cos(pitch); + out[k].dt_s = t; + out[k].horizontal_m = h; + } + return out; +} + +std::vector frame_weights(const std::vector& offs, double edge_weight) { + size_t n = offs.size(); + std::vector w(n, 1.0); + double edge = std::min(std::max(edge_weight, 0.0), 1.0); + if (n == 0 || edge >= 1.0) return w; + double dmax = 0.0; + for (auto& o : offs) dmax = std::max(dmax, std::abs(o.dt_s)); + if (dmax <= 1e-12) return w; + for (size_t i = 0; i < n; ++i) { + double r = std::min(std::abs(offs[i].dt_s) / dmax, 1.0); + w[i] = edge + (1.0 - edge) * 0.5 * (1.0 + std::cos(M_PI * r)); + } + return w; +} + +std::array point_to_scene(double x, double y, double z, const MotionConfig& mot, + const SARConfig& sar) { + if (!sar.coord_mode_ground) return {x, y, z}; + double p = mot.pitch_deg * M_PI / 180.0; + double ys = y * std::cos(p) - z * std::sin(p); + double zs = y * std::sin(p) + z * std::cos(p); + return {x, ys, zs}; +} + +Grid build_grid_scene(const std::array& c, const SARConfig& sar) { + Grid g; + auto lin = [](double a, double b, int n, std::vector& out) { + out.resize(n); + for (int i = 0; i < n; ++i) out[i] = a + (b - a) * i / (double)(n - 1); + }; + lin(c[0] - sar.half_width_m[0], c[0] + sar.half_width_m[0], sar.grid[0], g.x); + lin(c[1] - sar.half_width_m[1], c[1] + sar.half_width_m[1], sar.grid[1], g.y); + lin(c[2] - sar.half_width_m[2], c[2] + sar.half_width_m[2], sar.grid[2], g.z); + return g; +} + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/src/kernel_src.hpp b/data_acq_and_processing/gpr_cpp/src/kernel_src.hpp new file mode 100644 index 0000000..eb6621a --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/src/kernel_src.hpp @@ -0,0 +1,72 @@ +// Исходник CUDA-ядра для NVRTC (порт bp_aperture_contribution из gpr_core/kernels.py). +// Self-contained complex (без cupy-заголовков). +#pragma once + +namespace gpr { + +inline const char* APERTURE_KERNEL_SRC = R"CUDA( +struct cf { float x, y; }; +__device__ __forceinline__ cf operator+(cf a, cf b){ return {a.x+b.x, a.y+b.y}; } +__device__ __forceinline__ cf operator-(cf a, cf b){ return {a.x-b.x, a.y-b.y}; } +__device__ __forceinline__ cf operator*(float s, cf a){ return {s*a.x, s*a.y}; } + +// Весь взвешенный объём апертуры за один запуск: один поток на воксель, +// внутренний цикл по n_contrib = n_frames*n_pairs. acc/cnt пишутся один раз. +extern "C" __global__ void bp_aperture_contribution( + const cf* __restrict__ H, + const float t0, const float inv_dt, const int n_fft, const float inv_v, + const float* __restrict__ xt, const float* __restrict__ yt, const float* __restrict__ zt, + const float* __restrict__ xr, const float* __restrict__ yr, const float* __restrict__ zr, + const float* __restrict__ inv_geo_ref, const float* __restrict__ wgt, + const float* __restrict__ gx, const float* __restrict__ gy, const float* __restrict__ gz, + const int NX, const int NY, const int NZ, const int n_contrib, + const float cos_p, const float sin_p, + const float min_depth, const float max_depth, + const float comp_pow, const float w_max, const int use_comp, + cf* __restrict__ acc, float* __restrict__ cnt) +{ + const long total = (long)NX * NY * NZ; + for (long gid = blockIdx.x * blockDim.x + threadIdx.x; + gid < total; gid += (long)gridDim.x * blockDim.x) { + int ix = gid % NX; + int iz = (gid / NX) % NZ; + int iy = gid / ((long)NX * NZ); + float x = gx[ix]; + float ys = gy[iy], zs = gz[iz]; + float y_rad = ys * cos_p + zs * sin_p; + float z_rad = -ys * sin_p + zs * cos_p; + bool depth_ok = (z_rad >= min_depth) && (z_rad <= max_depth); + + cf a = {0.f, 0.f}; + float c = 0.f; + for (int i = 0; i < n_contrib; ++i) { + float w = wgt[i]; + float dxt = x - xt[i], dyt = y_rad - yt[i], dzt = z_rad - zt[i]; + float dxr = x - xr[i], dyr = y_rad - yr[i], dzr = z_rad - zr[i]; + float Rt = sqrtf(dxt*dxt + dyt*dyt + dzt*dzt); + float Rr = sqrtf(dxr*dxr + dyr*dyr + dzr*dzr); + float pos = ((Rt + Rr) * inv_v - t0) * inv_dt; + bool inside = (pos >= 0.f) && (pos <= (float)(n_fft - 1)); + if (!inside) continue; + int k = (int)floorf(pos); + if (k < 0) k = 0; + if (k > n_fft - 2) k = n_fft - 2; + float f = pos - (float)k; + const cf* h = H + (long)i * n_fft; + cf hv = h[k] + f * (h[k+1] - h[k]); + if (use_comp) { + float ratio = (1.f/(Rt*Rr + 1e-12f) + 1e-30f) * inv_geo_ref[i]; + float denom = (comp_pow == 0.25f) ? sqrtf(sqrtf(ratio)) : powf(ratio, comp_pow); + float ww = 1.f/(denom + 1e-12f); + ww = fminf(fmaxf(ww, 0.f), w_max); + hv = ww * hv; + } + if (depth_ok) { a = a + (w * hv); c += w; } + } + acc[gid] = a; + cnt[gid] = c; + } +} +)CUDA"; + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/src/nvrtc_kernel.cpp b/data_acq_and_processing/gpr_cpp/src/nvrtc_kernel.cpp new file mode 100644 index 0000000..0fce0d6 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/src/nvrtc_kernel.cpp @@ -0,0 +1,48 @@ +#include "gpr/nvrtc_kernel.hpp" + +#include + +#include +#include +#include + +#include "gpr/device.hpp" + +namespace gpr { + +static void nvrtc_check(nvrtcResult r, const char* what) { + if (r != NVRTC_SUCCESS) + throw std::runtime_error(std::string("NVRTC error in ") + what + ": " + + nvrtcGetErrorString(r)); +} + +NvrtcKernel::~NvrtcKernel() { + if (module_) cuModuleUnload(module_); +} + +void NvrtcKernel::compile(const char* src, const std::string& arch_flag, const char* func_name) { + nvrtcProgram prog; + nvrtc_check(nvrtcCreateProgram(&prog, src, "gpr_kernel.cu", 0, nullptr, nullptr), + "createProgram"); + std::vector opts{arch_flag.c_str()}; + nvrtcResult cr = nvrtcCompileProgram(prog, (int)opts.size(), opts.data()); + if (cr != NVRTC_SUCCESS) { + size_t logsz = 0; + nvrtcGetProgramLogSize(prog, &logsz); + std::string log(logsz, '\0'); + if (logsz) nvrtcGetProgramLog(prog, &log[0]); + nvrtcDestroyProgram(&prog); + throw std::runtime_error("NVRTC compile failed:\n" + log); + } + // cubin под точную арх (без JIT PTX — драйвер грузит SASS напрямую) + size_t cubsz = 0; + nvrtc_check(nvrtcGetCUBINSize(prog, &cubsz), "getCUBINSize"); + std::string cubin(cubsz, '\0'); + nvrtc_check(nvrtcGetCUBIN(prog, &cubin[0]), "getCUBIN"); + nvrtcDestroyProgram(&prog); + + cu_check(cuModuleLoadDataEx(&module_, cubin.data(), 0, nullptr, nullptr), "moduleLoad"); + cu_check(cuModuleGetFunction(&func_, module_, func_name), "getFunction"); +} + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/src/sar3d.cpp b/data_acq_and_processing/gpr_cpp/src/sar3d.cpp new file mode 100644 index 0000000..e32e265 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/src/sar3d.cpp @@ -0,0 +1,111 @@ +#include "gpr/sar3d.hpp" + +#include +#include +#include + +#include "kernel_src.hpp" + +namespace gpr { + +SarEngine::SarEngine(const GeometryConfig& g, const BPConfig& bp, const MotionConfig& mot, + const SARConfig& sar) + : geom_(g), bp_(bp), mot_(mot), sar_(sar) { + geom_.finalize(); + ctx_ = std::make_unique(); + kernel_ = std::make_unique(); + kernel_->compile(APERTURE_KERNEL_SRC, ctx_->arch_flag(), "bp_aperture_contribution"); +} + +VolumeRaw SarEngine::aperture_volume(const std::vector& frames, const Grid& grid) { + if (frames.empty()) throw std::runtime_error("aperture_volume: нет кадров"); + auto pairs = geom_.pair_keys(); + int n_pairs = (int)pairs.size(); + int n_frames = (int)frames.size(); + int n_contrib = n_frames * n_pairs; + + const AScanSet& a0 = *frames[0].ascans; + int n_fft = a0.n_fft; + double t0 = a0.t0, dt = a0.dt; + + int NX = grid.NX(), NY = grid.NY(), NZ = grid.NZ(); + size_t nvox = (size_t)NX * NY * NZ; + + // Плоские массивы по (кадр, пара) + std::vector H_stack((size_t)n_contrib * n_fft); + std::vector xt(n_contrib), yt(n_contrib), zt(n_contrib); + std::vector xr(n_contrib), yr(n_contrib), zr(n_contrib); + std::vector igr(n_contrib), wgt(n_contrib); + int idx = 0; + for (int f = 0; f < n_frames; ++f) { + const AScanSet& a = *frames[f].ascans; + if (a.n_fft != n_fft || std::abs(a.t0 - t0) > 1e-18) + throw std::runtime_error("кадры апертуры: разная ось времени"); + double dx = frames[f].offset.dx_m, dy = frames[f].offset.dy_m, dz = frames[f].offset.dz_m; + float w = (float)frames[f].weight; + for (int p = 0; p < n_pairs; ++p) { + int i = pairs[p].first, j = pairs[p].second; + // A-скан пары p → в стек (каст c128→c64) + const c128* src = a.pair_ptr(p); + c64* dst = H_stack.data() + (size_t)idx * n_fft; + for (int k = 0; k < n_fft; ++k) dst[k] = c64((float)src[k].real(), (float)src[k].imag()); + xt[idx] = (float)(geom_.x_tx[i] + dx); yt[idx] = (float)(geom_.y_tx[i] + dy); + zt[idx] = (float)(geom_.z_tx[i] + dz); + xr[idx] = (float)(geom_.x_rx[j] + dx); yr[idx] = (float)(geom_.y_rx[j] + dy); + zr[idx] = (float)(geom_.z_rx[j] + dz); + igr[idx] = (float)(1.0 / pair_geo_ref(geom_, bp_, i, j, geom_.bp_plane_y)); + wgt[idx] = w; + ++idx; + } + } + + std::vector gx(grid.x.begin(), grid.x.end()); + std::vector gy(grid.y.begin(), grid.y.end()); + std::vector gz(grid.z.begin(), grid.z.end()); + + float cos_p = 1.f, sin_p = 0.f; + if (sar_.coord_mode_ground) { + double p = mot_.pitch_deg * M_PI / 180.0; + cos_p = (float)std::cos(p); sin_p = (float)std::sin(p); + } + + // Заливка на устройство + DeviceBuffer dH(H_stack.size()); dH.upload(H_stack.data(), H_stack.size()); + DeviceBuffer dxt(n_contrib), dyt(n_contrib), dzt(n_contrib); + DeviceBuffer dxr(n_contrib), dyr(n_contrib), dzr(n_contrib); + DeviceBuffer digr(n_contrib), dwgt(n_contrib); + dxt.upload(xt.data(), n_contrib); dyt.upload(yt.data(), n_contrib); dzt.upload(zt.data(), n_contrib); + dxr.upload(xr.data(), n_contrib); dyr.upload(yr.data(), n_contrib); dzr.upload(zr.data(), n_contrib); + digr.upload(igr.data(), n_contrib); dwgt.upload(wgt.data(), n_contrib); + DeviceBuffer dgx(NX), dgy(NY), dgz(NZ); + dgx.upload(gx.data(), NX); dgy.upload(gy.data(), NY); dgz.upload(gz.data(), NZ); + DeviceBuffer dacc(nvox); + DeviceBuffer dcnt(nvox); + + float f_t0 = (float)t0, f_inv_dt = (float)(1.0 / dt), f_inv_v = (float)(1.0 / geom_.v); + float min_d = (float)bp_.min_depth, max_d = (float)bp_.max_depth; + float comp_pow = (float)bp_.comp_range_power, w_max = (float)bp_.comp_weight_max; + int use_comp = sar_.use_compensation ? 1 : 0; + + void* args[] = { + dH.addr(), &f_t0, &f_inv_dt, &n_fft, &f_inv_v, + dxt.addr(), dyt.addr(), dzt.addr(), dxr.addr(), dyr.addr(), dzr.addr(), + digr.addr(), dwgt.addr(), dgx.addr(), dgy.addr(), dgz.addr(), + &NX, &NY, &NZ, &n_contrib, &cos_p, &sin_p, &min_d, &max_d, + &comp_pow, &w_max, &use_comp, dacc.addr(), dcnt.addr()}; + + int block = 256; + int gridDim = (int)std::min(65535, ((long)nvox + block - 1) / block); + cu_check(cuLaunchKernel(kernel_->func(), gridDim, 1, 1, block, 1, 1, 0, 0, args, nullptr), + "launch aperture"); + cu_check(cuCtxSynchronize(), "sync"); + + VolumeRaw out; + out.NX = NX; out.NY = NY; out.NZ = NZ; + out.acc.resize(nvox); out.cnt.resize(nvox); + dacc.download(out.acc.data(), nvox); + dcnt.download(out.cnt.data(), nvox); + return out; +} + +} // namespace gpr diff --git a/data_acq_and_processing/gpr_cpp/third_party/pocketfft/pocketfft_hdronly.h b/data_acq_and_processing/gpr_cpp/third_party/pocketfft/pocketfft_hdronly.h new file mode 100644 index 0000000..de6c8d6 --- /dev/null +++ b/data_acq_and_processing/gpr_cpp/third_party/pocketfft/pocketfft_hdronly.h @@ -0,0 +1,3853 @@ +/* +This file is part of pocketfft. + +Copyright (C) 2010-2024 Max-Planck-Society +Copyright (C) 2019-2020 Peter Bell + +For the odd-sized DCT-IV transforms: + Copyright (C) 2003, 2007-14 Matteo Frigo + Copyright (C) 2003, 2007-14 Massachusetts Institute of Technology + +For the prev_good_size search: + Copyright (C) 2024 Tan Ping Liang, Peter Bell + +For the safeguards against integer overflow in good_size search: + Copyright (C) 2024 Cris Luengo + +Authors: Martin Reinecke, Peter Bell + +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. +* Redistributions in binary form must reproduce the above copyright notice, this + list of conditions and the following disclaimer in the documentation and/or + other materials provided with the distribution. +* Neither the name of the copyright holder nor the names of its contributors may + be used to endorse or promote products derived from this software without + specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#ifndef POCKETFFT_HDRONLY_H +#define POCKETFFT_HDRONLY_H + +#ifndef __cplusplus +#error This file is C++ and requires a C++ compiler. +#endif + +#if !(__cplusplus >= 201103L || (defined(_MSVC_LANG) && _MSVC_LANG >= 201103L)) +#error This file requires at least C++11 support. +#endif + +#ifndef POCKETFFT_NAMESPACE +#define POCKETFFT_NAMESPACE pocketfft +#endif + +#ifndef POCKETFFT_CACHE_SIZE +#define POCKETFFT_CACHE_SIZE 0 +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#if POCKETFFT_CACHE_SIZE!=0 +#include +#include +#endif + +#ifndef POCKETFFT_NO_MULTITHREADING +#include +#include +#include +#include +#include +#include +#include + +#ifdef POCKETFFT_PTHREADS +# include +#endif +#endif + +#if defined(__GNUC__) +#define POCKETFFT_NOINLINE __attribute__((noinline)) +#define POCKETFFT_RESTRICT __restrict__ +#elif defined(_MSC_VER) +#define POCKETFFT_NOINLINE __declspec(noinline) +#define POCKETFFT_RESTRICT __restrict +#else +#define POCKETFFT_NOINLINE +#define POCKETFFT_RESTRICT +#endif + +namespace POCKETFFT_NAMESPACE { + +namespace detail { +using std::size_t; +using std::ptrdiff_t; + +// Always use std:: for functions +template T cos(T) = delete; +template T sin(T) = delete; +template T sqrt(T) = delete; + +using shape_t = std::vector; +using stride_t = std::vector; + +constexpr bool FORWARD = true, + BACKWARD = false; + +// only enable vector support for gcc>=5.0 and clang>=5.0 +#ifndef POCKETFFT_NO_VECTORS +#define POCKETFFT_NO_VECTORS +#if defined(__INTEL_COMPILER) +// do nothing. This is necessary because this compiler also sets __GNUC__. +#elif defined(__clang__) +// AppleClang has their own version numbering +#ifdef __apple_build_version__ +# if (__clang_major__ > 9) || (__clang_major__ == 9 && __clang_minor__ >= 1) +# undef POCKETFFT_NO_VECTORS +# endif +#elif __clang_major__ >= 5 +# undef POCKETFFT_NO_VECTORS +#endif +#elif defined(__GNUC__) +#if __GNUC__>=5 +#undef POCKETFFT_NO_VECTORS +#endif +#endif +#endif + +template struct VLEN { static constexpr size_t val=1; }; + +#ifndef POCKETFFT_NO_VECTORS +#if (defined(__AVX512F__)) +template<> struct VLEN { static constexpr size_t val=16; }; +template<> struct VLEN { static constexpr size_t val=8; }; +#elif (defined(__AVX__)) +template<> struct VLEN { static constexpr size_t val=8; }; +template<> struct VLEN { static constexpr size_t val=4; }; +#elif (defined(__SSE2__)) +template<> struct VLEN { static constexpr size_t val=4; }; +template<> struct VLEN { static constexpr size_t val=2; }; +#elif (defined(__VSX__)) +template<> struct VLEN { static constexpr size_t val=4; }; +template<> struct VLEN { static constexpr size_t val=2; }; +#elif (defined(__ARM_NEON__) || defined(__ARM_NEON)) +template<> struct VLEN { static constexpr size_t val=4; }; +template<> struct VLEN { static constexpr size_t val=2; }; +#elif (defined(__riscv_vector) || defined(__riscv_v_intrinsic)) +template<> struct VLEN { static constexpr size_t val=8; }; +template<> struct VLEN { static constexpr size_t val=4; }; +#else +#define POCKETFFT_NO_VECTORS +#endif +#endif + +// std::aligned_alloc is a bit cursed ... it doesn't exist on MacOS < 10.15 +// and in musl, and other OSes seem to have even more peculiarities. +// Let's unconditionally work around it for now. +#if defined(POCKETFFT_USE_POSIX_MEMALIGN) && (defined(__APPLE__) || defined(__unix__)) +// Use posix_memalign on POSIX systems when explicitly enabled. +// The portable aligned_alloc below stores metadata at ptr[-1], which conflicts +// with ASAN's heap redzone and causes intermittent bus errors. +inline void *aligned_alloc(size_t align, size_t size) + { + align = std::max(align, sizeof(void*)); // posix_memalign requires align >= sizeof(void*) + void *ptr = nullptr; + if (posix_memalign(&ptr, align, size) != 0) + throw std::bad_alloc(); + return ptr; + } +inline void aligned_dealloc(void *ptr) + { free(ptr); } +#else +# if 0 +//#if (__cplusplus >= 201703L) && (!defined(__MINGW32__)) && (!defined(_MSC_VER)) && (__MAC_OS_X_VERSION_MIN_REQUIRED >= MAC_OS_X_VERSION_10_15) +inline void *aligned_alloc(size_t align, size_t size) + { + // aligned_alloc() requires that the requested size is a multiple of "align" + void *ptr = ::aligned_alloc(align,(size+align-1)&(~(align-1))); + if (!ptr) throw std::bad_alloc(); + return ptr; + } +inline void aligned_dealloc(void *ptr) + { free(ptr); } +#else // portable emulation +inline void *aligned_alloc(size_t align, size_t size) + { + align = std::max(align, alignof(max_align_t)); + void *ptr = malloc(size+align); + if (!ptr) throw std::bad_alloc(); + void *res = reinterpret_cast + ((reinterpret_cast(ptr) & ~(uintptr_t(align-1))) + uintptr_t(align)); + (reinterpret_cast(res))[-1] = ptr; + return res; + } +inline void aligned_dealloc(void *ptr) + { if (ptr) free((reinterpret_cast(ptr))[-1]); } +# endif +#endif + +template class arr + { + private: + T *p; + size_t sz; + +#if defined(POCKETFFT_NO_VECTORS) + static T *ralloc(size_t num) + { + if (num==0) return nullptr; + void *res = malloc(num*sizeof(T)); + if (!res) throw std::bad_alloc(); + return reinterpret_cast(res); + } + static void dealloc(T *ptr) + { free(ptr); } +#else + static T *ralloc(size_t num) + { + if (num==0) return nullptr; + void *ptr = aligned_alloc(64, num*sizeof(T)); + return static_cast(ptr); + } + static void dealloc(T *ptr) + { aligned_dealloc(ptr); } +#endif + + public: + arr() : p(nullptr), sz(0) {} + explicit arr(size_t n) : p(ralloc(n)), sz(n) {} + arr(arr &&other) noexcept + : p(other.p), sz(other.sz) + { other.p=nullptr; other.sz=0; } + ~arr() { dealloc(p); } + + void resize(size_t n) + { + if (n==sz) return; + dealloc(p); + p = ralloc(n); + sz = n; + } + + T &operator[](size_t idx) { return p[idx]; } + const T &operator[](size_t idx) const { return p[idx]; } + + T *data() { return p; } + const T *data() const { return p; } + + size_t size() const { return sz; } + }; + +template struct cmplx { + T r, i; + cmplx() = default; + cmplx(T r_, T i_) : r(r_), i(i_) {} + void Set(T r_, T i_) { r=r_; i=i_; } + void Set(T r_) { r=r_; i=T(0); } + cmplx &operator+= (const cmplx &other) + { r+=other.r; i+=other.i; return *this; } + templatecmplx &operator*= (T2 other) + { r*=other; i*=other; return *this; } + templatecmplx &operator*= (const cmplx &other) + { + T tmp = r*other.r - i*other.i; + i = r*other.i + i*other.r; + r = tmp; + return *this; + } + templatecmplx &operator+= (const cmplx &other) + { r+=other.r; i+=other.i; return *this; } + templatecmplx &operator-= (const cmplx &other) + { r-=other.r; i-=other.i; return *this; } + template auto operator* (const T2 &other) const + -> cmplx + { return {r*other, i*other}; } + template auto operator+ (const cmplx &other) const + -> cmplx + { return {r+other.r, i+other.i}; } + template auto operator- (const cmplx &other) const + -> cmplx + { return {r-other.r, i-other.i}; } + template auto operator* (const cmplx &other) const + -> cmplx + { return {r*other.r-i*other.i, r*other.i + i*other.r}; } + template auto special_mul (const cmplx &other) const + -> cmplx + { + using Tres = cmplx; + return fwd ? Tres(r*other.r+i*other.i, i*other.r-r*other.i) + : Tres(r*other.r-i*other.i, r*other.i+i*other.r); + } +}; +template inline void PM(T &a, T &b, T c, T d) + { a=c+d; b=c-d; } +template inline void PMINPLACE(T &a, T &b) + { T t = a; a+=b; b=t-b; } +template inline void MPINPLACE(T &a, T &b) + { T t = a; a-=b; b=t+b; } +template cmplx conj(const cmplx &a) + { return {a.r, -a.i}; } +template void special_mul (const cmplx &v1, const cmplx &v2, cmplx &res) + { + res = fwd ? cmplx(v1.r*v2.r+v1.i*v2.i, v1.i*v2.r-v1.r*v2.i) + : cmplx(v1.r*v2.r-v1.i*v2.i, v1.r*v2.i+v1.i*v2.r); + } + +template void ROT90(cmplx &a) + { auto tmp_=a.r; a.r=-a.i; a.i=tmp_; } +template void ROTX90(cmplx &a) + { auto tmp_= fwd ? -a.r : a.r; a.r = fwd ? a.i : -a.i; a.i=tmp_; } + +// +// twiddle factor section +// +template class sincos_2pibyn + { + private: + using Thigh = typename std::conditional<(sizeof(T)>sizeof(double)), T, double>::type; + size_t N, mask, shift; + arr> v1, v2; + + static cmplx calc(size_t x, size_t n, Thigh ang) + { + x<<=3; + if (x<4*n) // first half + { + if (x<2*n) // first quadrant + { + if (x(std::cos(Thigh(x)*ang), std::sin(Thigh(x)*ang)); + return cmplx(std::sin(Thigh(2*n-x)*ang), std::cos(Thigh(2*n-x)*ang)); + } + else // second quadrant + { + x-=2*n; + if (x(-std::sin(Thigh(x)*ang), std::cos(Thigh(x)*ang)); + return cmplx(-std::cos(Thigh(2*n-x)*ang), std::sin(Thigh(2*n-x)*ang)); + } + } + else + { + x=8*n-x; + if (x<2*n) // third quadrant + { + if (x(std::cos(Thigh(x)*ang), -std::sin(Thigh(x)*ang)); + return cmplx(std::sin(Thigh(2*n-x)*ang), -std::cos(Thigh(2*n-x)*ang)); + } + else // fourth quadrant + { + x-=2*n; + if (x(-std::sin(Thigh(x)*ang), -std::cos(Thigh(x)*ang)); + return cmplx(-std::cos(Thigh(2*n-x)*ang), -std::sin(Thigh(2*n-x)*ang)); + } + } + } + + public: + POCKETFFT_NOINLINE explicit sincos_2pibyn(size_t n) + : N(n) + { + constexpr auto pi = 3.141592653589793238462643383279502884197L; + Thigh ang = Thigh(0.25L*pi/n); + size_t nval = (n+2)/2; + shift = 1; + while((size_t(1)< operator[](size_t idx) const + { + if (2*idx<=N) + { + auto x1=v1[idx&mask], x2=v2[idx>>shift]; + return cmplx(T(x1.r*x2.r-x1.i*x2.i), T(x1.r*x2.i+x1.i*x2.r)); + } + idx = N-idx; + auto x1=v1[idx&mask], x2=v2[idx>>shift]; + return cmplx(T(x1.r*x2.r-x1.i*x2.i), -T(x1.r*x2.i+x1.i*x2.r)); + } + }; + +struct util // hack to avoid duplicate symbols + { + static POCKETFFT_NOINLINE size_t largest_prime_factor (size_t n) + { + size_t res=1; + while ((n&1)==0) + { res=2; n>>=1; } + for (size_t x=3; x*x<=n; x+=2) + while ((n%x)==0) + { res=x; n/=x; } + if (n>1) res=n; + return res; + } + + static POCKETFFT_NOINLINE double cost_guess (size_t n) + { + constexpr double lfp=1.1; // penalty for non-hardcoded larger factors + size_t ni=n; + double result=0.; + while ((n&1)==0) + { result+=2; n>>=1; } + for (size_t x=3; x*x<=n; x+=2) + while ((n%x)==0) + { + result+= (x<=5) ? double(x) : lfp*double(x); // penalize larger prime factors + n/=x; + } + if (n>1) result+=(n<=5) ? double(n) : lfp*double(n); + return result*double(ni); + } + + /* inner workings of good_size_cmplx() */ + template + static POCKETFFT_NOINLINE UIntT good_size_cmplx_typed(UIntT n) + { + static_assert(std::numeric_limits::is_integer && (!std::numeric_limits::is_signed), + "type must be unsigned integer"); + if (n<=12) return n; + if (n>std::numeric_limits::max()/11/2) + { + // The algorithm below doesn't work for this value, the multiplication can overflow. + if (sizeof(UIntT)(n); + if (res<=std::numeric_limits::max()) + return static_cast(res); + } + // Otherwise, this size is ridiculously large, people shouldn't be computing FFTs this large. + throw std::runtime_error("FFT size is too large."); + } + + UIntT bestfac=2*n; + for (UIntT f11=1; f11n) + { + if (x>=1; + } + else + return n; + } + } + return bestfac; + } + /* returns the smallest composite of 2, 3, 5, 7 and 11 which is >= n */ + static POCKETFFT_NOINLINE size_t good_size_cmplx(size_t n) + { + return good_size_cmplx_typed(n); + } + /* returns the smallest composite of 2, 3, 5, 7 and 11 which is >= n + and a multiple of required_factor. */ + static POCKETFFT_NOINLINE size_t good_size_cmplx(size_t n, + size_t required_factor) + { + if (required_factor<1) + throw std::runtime_error("required factor must not be 0"); + return good_size_cmplx((n+required_factor-1)/required_factor) * required_factor; + } + + /* inner workings of good_size_real() */ + template + static POCKETFFT_NOINLINE UIntT good_size_real_typed(UIntT n) + { + static_assert(std::numeric_limits::is_integer && (!std::numeric_limits::is_signed), + "type must be unsigned integer"); + if (n<=6) return n; + if (n>std::numeric_limits::max()/5/2) + { + // The algorithm below doesn't work for this value, the multiplication can overflow. + if (sizeof(UIntT)(n); + if (res<=std::numeric_limits::max()) + return static_cast(res); + } + // Otherwise, this size is ridiculously large, people shouldn't be computing FFTs this large. + throw std::runtime_error("FFT size is too large."); + } + + UIntT bestfac=2*n; + for (UIntT f5=1; f5n) + { + if (x>=1; + } + else + return n; + } + } + return bestfac; + } + /* returns the smallest composite of 2, 3, 5 which is >= n */ + static POCKETFFT_NOINLINE size_t good_size_real(size_t n) + { + return good_size_real_typed(n); + } + /* returns the smallest composite of 2, 3, 5 which is >= n + and a multiple of required_factor. */ + static POCKETFFT_NOINLINE size_t good_size_real(size_t n, + size_t required_factor) + { + if (required_factor<1) + throw std::runtime_error("required factor must not be 0"); + return good_size_real((n+required_factor-1)/required_factor) * required_factor; + } + + /* inner workings of prev_good_size_cmplx() */ + template + static POCKETFFT_NOINLINE UIntT prev_good_size_cmplx_typed(UIntT n) + { + static_assert(std::numeric_limits::is_integer && (!std::numeric_limits::is_signed), + "type must be unsigned integer"); + if (n<=12) return n; + if (n>std::numeric_limits::max()/11) + { + // The algorithm below doesn't work for this value, the multiplication can overflow. + if (sizeof(UIntT)(n); + if (res<=std::numeric_limits::max()) + return static_cast(res); + } + // Otherwise, this size is ridiculously large, people shouldn't be computing FFTs this large. + throw std::runtime_error("FFT size is too large."); + } + + UIntT bestfound = 1; + for (UIntT f11 = 1;f11 <= n; f11 *= 11) + for (UIntT f117 = f11; f117 <= n; f117 *= 7) + for (UIntT f1175 = f117; f1175 <= n; f1175 *= 5) + { + UIntT x = f1175; + while (x*2 <= n) x *= 2; + if (x > bestfound) bestfound = x; + while (true) + { + if (x * 3 <= n) x *= 3; + else if (x % 2 == 0) x /= 2; + else break; + + if (x > bestfound) bestfound = x; + } + } + return bestfound; + } + /* returns the largest composite of 2, 3, 5, 7 and 11 which is <= n */ + static POCKETFFT_NOINLINE size_t prev_good_size_cmplx(size_t n) + { + return prev_good_size_cmplx_typed(n); + } + + /* inner workings of prev_good_size_real() */ + template + static POCKETFFT_NOINLINE UIntT prev_good_size_real_typed(UIntT n) + { + static_assert(std::numeric_limits::is_integer && (!std::numeric_limits::is_signed), + "type must be unsigned integer"); + if (n<=6) return n; + if (n>std::numeric_limits::max()/5) + { + // The algorithm below doesn't work for this value, the multiplication can overflow. + if (sizeof(UIntT)(n); + if (res<=std::numeric_limits::max()) + return static_cast(res); + } + // Otherwise, this size is ridiculously large, people shouldn't be computing FFTs this large. + throw std::runtime_error("FFT size is too large."); + } + + UIntT bestfound = 1; + for (UIntT f5 = 1; f5 <= n; f5 *= 5) + { + UIntT x = f5; + while (x*2 <= n) x *= 2; + if (x > bestfound) bestfound = x; + while (true) + { + if (x * 3 <= n) x *= 3; + else if (x % 2 == 0) x /= 2; + else break; + + if (x > bestfound) bestfound = x; + } + } + return bestfound; + } + /* returns the largest composite of 2, 3, 5 which is <= n */ + static POCKETFFT_NOINLINE size_t prev_good_size_real(size_t n) + { + return prev_good_size_real_typed(n); + } + + static size_t prod(const shape_t &shape) + { + size_t res=1; + for (auto sz: shape) + res*=sz; + return res; + } + + static POCKETFFT_NOINLINE void sanity_check(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, bool inplace) + { + auto ndim = shape.size(); + if (ndim<1) throw std::runtime_error("ndim must be >= 1"); + if ((stride_in.size()!=ndim) || (stride_out.size()!=ndim)) + throw std::runtime_error("stride dimension mismatch"); + if (inplace && (stride_in!=stride_out)) + throw std::runtime_error("stride mismatch"); + } + + static POCKETFFT_NOINLINE void sanity_check(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, bool inplace, + const shape_t &axes) + { + sanity_check(shape, stride_in, stride_out, inplace); + auto ndim = shape.size(); + shape_t tmp(ndim,0); + for (auto ax : axes) + { + if (ax>=ndim) throw std::invalid_argument("bad axis number"); + if (++tmp[ax]>1) throw std::invalid_argument("axis specified repeatedly"); + } + } + + static POCKETFFT_NOINLINE void sanity_check(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, bool inplace, + size_t axis) + { + sanity_check(shape, stride_in, stride_out, inplace); + if (axis>=shape.size()) throw std::invalid_argument("bad axis number"); + } + +#ifdef POCKETFFT_NO_MULTITHREADING + static size_t thread_count (size_t /*nthreads*/, const shape_t &/*shape*/, + size_t /*axis*/, size_t /*vlen*/) + { return 1; } +#else + static size_t thread_count (size_t nthreads, const shape_t &shape, + size_t axis, size_t vlen) + { + if (nthreads==1) return 1; + size_t size = prod(shape); + size_t parallel = size / (shape[axis] * vlen); + if (shape[axis] < 1000) + parallel /= 4; + size_t max_threads = nthreads == 0 ? + std::thread::hardware_concurrency() : nthreads; + return std::max(size_t(1), std::min(parallel, max_threads)); + } +#endif + }; + +namespace threading { + +#ifdef POCKETFFT_NO_MULTITHREADING + +constexpr inline size_t thread_id() { return 0; } +constexpr inline size_t num_threads() { return 1; } + +template +void thread_map(size_t /* nthreads */, Func f) + { f(); } + +#else + +inline size_t &thread_id() + { + static thread_local size_t thread_id_=0; + return thread_id_; + } +inline size_t &num_threads() + { + static thread_local size_t num_threads_=1; + return num_threads_; + } +static const size_t max_threads = std::max(1u, std::thread::hardware_concurrency()); + +class latch + { + std::atomic num_left_; + std::mutex mut_; + std::condition_variable completed_; + using lock_t = std::unique_lock; + + public: + explicit latch(size_t n): num_left_(n) {} + + void count_down() + { + lock_t lock(mut_); + if (--num_left_) + return; + completed_.notify_all(); + } + + void wait() + { + lock_t lock(mut_); + completed_.wait(lock, [this]{ return is_ready(); }); + } + bool is_ready() { return num_left_ == 0; } + }; + +template class concurrent_queue + { + std::queue q_; + std::mutex mut_; + std::atomic size_; + using lock_t = std::lock_guard; + + public: + + void push(T val) + { + lock_t lock(mut_); + ++size_; + q_.push(std::move(val)); + } + + bool try_pop(T &val) + { + if (size_ == 0) return false; + lock_t lock(mut_); + // Queue might have been emptied while we acquired the lock + if (q_.empty()) return false; + + val = std::move(q_.front()); + --size_; + q_.pop(); + return true; + } + + bool empty() const { return size_==0; } + }; + +// C++ allocator with support for over-aligned types +template struct aligned_allocator + { + using value_type = T; + template + explicit aligned_allocator(const aligned_allocator&) {} + aligned_allocator() = default; + + T *allocate(size_t n) + { + void* mem = aligned_alloc(alignof(T), n*sizeof(T)); + return static_cast(mem); + } + + void deallocate(T *p, size_t /*n*/) + { aligned_dealloc(p); } + }; + +class thread_pool + { + // A reasonable guess, probably close enough for most hardware + static constexpr size_t cache_line_size = 64; + struct alignas(cache_line_size) worker + { + std::thread thread; + std::condition_variable work_ready; + std::mutex mut; + std::atomic_flag busy_flag = ATOMIC_FLAG_INIT; + std::function work; + + void worker_main( + std::atomic &shutdown_flag, + std::atomic &unscheduled_tasks, + concurrent_queue> &overflow_work) + { + using lock_t_inner = std::unique_lock; + bool expect_work = true; + while (!shutdown_flag || expect_work) + { + std::function local_work; + if (expect_work || unscheduled_tasks == 0) + { + lock_t_inner lock(mut); + // Wait until there is work to be executed + work_ready.wait(lock, [&]{ return (work || shutdown_flag); }); + local_work.swap(work); + expect_work = false; + } + + bool marked_busy = false; + if (local_work) + { + marked_busy = true; + local_work(); + } + + if (!overflow_work.empty()) + { + if (!marked_busy && busy_flag.test_and_set()) + { + expect_work = true; + continue; + } + marked_busy = true; + + while (overflow_work.try_pop(local_work)) + { + --unscheduled_tasks; + local_work(); + } + } + + if (marked_busy) busy_flag.clear(); + } + } + }; + + concurrent_queue> overflow_work_; + std::mutex mut_; + std::vector> workers_; + std::atomic shutdown_; + std::atomic unscheduled_tasks_; + using lock_t = std::lock_guard; + + void create_threads() + { + lock_t lock(mut_); + size_t nthreads=workers_.size(); + for (size_t i=0; ibusy_flag.clear(); + worker_i->work = nullptr; + worker_i->thread = std::thread([worker_i, this] + { + worker_i->worker_main(shutdown_, unscheduled_tasks_, overflow_work_); + }); + } + catch (...) + { + shutdown_locked(); + throw; + } + } + } + + void shutdown_locked() + { + shutdown_ = true; + for (auto &worker_i : workers_) + worker_i.work_ready.notify_all(); + + for (auto &worker_i : workers_) + if (worker_i.thread.joinable()) + worker_i.thread.join(); + } + + public: + explicit thread_pool(size_t nthreads): + workers_(nthreads) + { create_threads(); } + + thread_pool(): thread_pool(max_threads) {} + + ~thread_pool() { shutdown(); } + + void submit(std::function work) + { + lock_t lock(mut_); + if (shutdown_) + throw std::runtime_error("Work item submitted after shutdown"); + + ++unscheduled_tasks_; + + // First check for any idle workers and wake those + for (auto &worker_i : workers_) + if (!worker_i.busy_flag.test_and_set()) + { + --unscheduled_tasks_; + { + lock_t lock_inner(worker_i.mut); + worker_i.work = std::move(work); + } + worker_i.work_ready.notify_one(); + return; + } + + // If no workers were idle, push onto the overflow queue for later + overflow_work_.push(std::move(work)); + } + + void shutdown() + { + lock_t lock(mut_); + shutdown_locked(); + } + + void restart() + { + shutdown_ = false; + create_threads(); + } + }; + +inline thread_pool & get_pool() + { + static thread_pool pool; +#ifdef POCKETFFT_PTHREADS + static std::once_flag f; + std::call_once(f, + []{ + pthread_atfork( + +[]{ get_pool().shutdown(); }, // prepare + +[]{ get_pool().restart(); }, // parent + +[]{ get_pool().restart(); } // child + ); + }); +#endif + + return pool; + } + +/** Map a function f over nthreads */ +template +void thread_map(size_t nthreads, Func f) + { + if (nthreads == 0) + nthreads = max_threads; + + if (nthreads == 1) + { f(); return; } + + auto & pool = get_pool(); + latch counter(nthreads); + std::exception_ptr ex; + std::mutex ex_mut; + for (size_t i=0; i lock(ex_mut); + ex = std::current_exception(); + } + counter.count_down(); + }); + } + counter.wait(); + if (ex) + std::rethrow_exception(ex); + } + +#endif + +} + +// +// complex FFTPACK transforms +// + +template class cfftp + { + private: + struct fctdata + { + size_t fct; + cmplx *tw, *tws; + }; + + size_t length; + arr> mem; + std::vector fact; + + void add_factor(size_t factor) + { fact.push_back({factor, nullptr, nullptr}); } + +template void pass2 (size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const cmplx * POCKETFFT_RESTRICT wa) const + { + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+2*c)]; }; + auto WA = [wa, ido](size_t x, size_t i) + { return wa[i-1+x*(ido-1)]; }; + + if (ido==1) + for (size_t k=0; k(CC(i,0,k)-CC(i,1,k),WA(0,i),CH(i,k,1)); + } + } + } + +#define POCKETFFT_PREP3(idx) \ + T t0 = CC(idx,0,k), t1, t2; \ + PM (t1,t2,CC(idx,1,k),CC(idx,2,k)); \ + CH(idx,k,0)=t0+t1; +#define POCKETFFT_PARTSTEP3a(u1,u2,twr,twi) \ + { \ + T ca=t0+t1*twr; \ + T cb{-t2.i*twi, t2.r*twi}; \ + PM(CH(0,k,u1),CH(0,k,u2),ca,cb) ;\ + } +#define POCKETFFT_PARTSTEP3b(u1,u2,twr,twi) \ + { \ + T ca=t0+t1*twr; \ + T cb{-t2.i*twi, t2.r*twi}; \ + special_mul(ca+cb,WA(u1-1,i),CH(i,k,u1)); \ + special_mul(ca-cb,WA(u2-1,i),CH(i,k,u2)); \ + } +template void pass3 (size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const cmplx * POCKETFFT_RESTRICT wa) const + { + constexpr T0 tw1r=-0.5, + tw1i= (fwd ? -1: 1) * T0(0.8660254037844386467637231707529362L); + + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+3*c)]; }; + auto WA = [wa, ido](size_t x, size_t i) + { return wa[i-1+x*(ido-1)]; }; + + if (ido==1) + for (size_t k=0; k void pass4 (size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const cmplx * POCKETFFT_RESTRICT wa) const + { + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+4*c)]; }; + auto WA = [wa, ido](size_t x, size_t i) + { return wa[i-1+x*(ido-1)]; }; + + if (ido==1) + for (size_t k=0; k(t4); + PM(CH(0,k,0),CH(0,k,2),t2,t3); + PM(CH(0,k,1),CH(0,k,3),t1,t4); + } + else + for (size_t k=0; k(t4); + PM(CH(0,k,0),CH(0,k,2),t2,t3); + PM(CH(0,k,1),CH(0,k,3),t1,t4); + } + for (size_t i=1; i(t4); + CH(i,k,0) = t2+t3; + special_mul(t1+t4,WA(0,i),CH(i,k,1)); + special_mul(t2-t3,WA(1,i),CH(i,k,2)); + special_mul(t1-t4,WA(2,i),CH(i,k,3)); + } + } + } + +#define POCKETFFT_PREP5(idx) \ + T t0 = CC(idx,0,k), t1, t2, t3, t4; \ + PM (t1,t4,CC(idx,1,k),CC(idx,4,k)); \ + PM (t2,t3,CC(idx,2,k),CC(idx,3,k)); \ + CH(idx,k,0).r=t0.r+t1.r+t2.r; \ + CH(idx,k,0).i=t0.i+t1.i+t2.i; + +#define POCKETFFT_PARTSTEP5a(u1,u2,twar,twbr,twai,twbi) \ + { \ + T ca,cb; \ + ca.r=t0.r+twar*t1.r+twbr*t2.r; \ + ca.i=t0.i+twar*t1.i+twbr*t2.i; \ + cb.i=twai*t4.r twbi*t3.r; \ + cb.r=-(twai*t4.i twbi*t3.i); \ + PM(CH(0,k,u1),CH(0,k,u2),ca,cb); \ + } + +#define POCKETFFT_PARTSTEP5b(u1,u2,twar,twbr,twai,twbi) \ + { \ + T ca,cb; \ + ca.r=t0.r+twar*t1.r+twbr*t2.r; \ + ca.i=t0.i+twar*t1.i+twbr*t2.i; \ + cb.i=twai*t4.r twbi*t3.r; \ + cb.r=-(twai*t4.i twbi*t3.i); \ + special_mul(ca+cb,WA(u1-1,i),CH(i,k,u1)); \ + special_mul(ca-cb,WA(u2-1,i),CH(i,k,u2)); \ + } +template void pass5 (size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const cmplx * POCKETFFT_RESTRICT wa) const + { + constexpr T0 tw1r= T0(0.3090169943749474241022934171828191L), + tw1i= (fwd ? -1: 1) * T0(0.9510565162951535721164393333793821L), + tw2r= T0(-0.8090169943749474241022934171828191L), + tw2i= (fwd ? -1: 1) * T0(0.5877852522924731291687059546390728L); + + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+5*c)]; }; + auto WA = [wa, ido](size_t x, size_t i) + { return wa[i-1+x*(ido-1)]; }; + + if (ido==1) + for (size_t k=0; k(da,WA(u1-1,i),CH(i,k,u1)); \ + special_mul(db,WA(u2-1,i),CH(i,k,u2)); \ + } + +template void pass7(size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const cmplx * POCKETFFT_RESTRICT wa) const + { + constexpr T0 tw1r= T0(0.6234898018587335305250048840042398L), + tw1i= (fwd ? -1 : 1) * T0(0.7818314824680298087084445266740578L), + tw2r= T0(-0.2225209339563144042889025644967948L), + tw2i= (fwd ? -1 : 1) * T0(0.9749279121818236070181316829939312L), + tw3r= T0(-0.9009688679024191262361023195074451L), + tw3i= (fwd ? -1 : 1) * T0(0.433883739117558120475768332848359L); + + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+7*c)]; }; + auto WA = [wa, ido](size_t x, size_t i) + { return wa[i-1+x*(ido-1)]; }; + + if (ido==1) + for (size_t k=0; k void ROTX45(T &a) const + { + constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L); + if (fwd) + { auto tmp_=a.r; a.r=hsqt2*(a.r+a.i); a.i=hsqt2*(a.i-tmp_); } + else + { auto tmp_=a.r; a.r=hsqt2*(a.r-a.i); a.i=hsqt2*(a.i+tmp_); } + } +template void ROTX135(T &a) const + { + constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L); + if (fwd) + { auto tmp_=a.r; a.r=hsqt2*(a.i-a.r); a.i=hsqt2*(-tmp_-a.i); } + else + { auto tmp_=a.r; a.r=hsqt2*(-a.r-a.i); a.i=hsqt2*(tmp_-a.i); } + } + +template void pass8 (size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const cmplx * POCKETFFT_RESTRICT wa) const + { + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+8*c)]; }; + auto WA = [wa, ido](size_t x, size_t i) + { return wa[i-1+x*(ido-1)]; }; + + if (ido==1) + for (size_t k=0; k(a3); + + ROTX90(a7); + PMINPLACE(a5,a7); + ROTX45(a5); + ROTX135(a7); + + PM(a0,a4,CC(0,0,k),CC(0,4,k)); + PM(a2,a6,CC(0,2,k),CC(0,6,k)); + PM(CH(0,k,0),CH(0,k,4),a0+a2,a1); + PM(CH(0,k,2),CH(0,k,6),a0-a2,a3); + ROTX90(a6); + PM(CH(0,k,1),CH(0,k,5),a4+a6,a5); + PM(CH(0,k,3),CH(0,k,7),a4-a6,a7); + } + else + for (size_t k=0; k(a3); + + ROTX90(a7); + PMINPLACE(a5,a7); + ROTX45(a5); + ROTX135(a7); + + PM(a0,a4,CC(0,0,k),CC(0,4,k)); + PM(a2,a6,CC(0,2,k),CC(0,6,k)); + PM(CH(0,k,0),CH(0,k,4),a0+a2,a1); + PM(CH(0,k,2),CH(0,k,6),a0-a2,a3); + ROTX90(a6); + PM(CH(0,k,1),CH(0,k,5),a4+a6,a5); + PM(CH(0,k,3),CH(0,k,7),a4-a6,a7); + } + for (size_t i=1; i(a7); + PMINPLACE(a1,a3); + ROTX90(a3); + PMINPLACE(a5,a7); + ROTX45(a5); + ROTX135(a7); + PM(a0,a4,CC(i,0,k),CC(i,4,k)); + PM(a2,a6,CC(i,2,k),CC(i,6,k)); + PMINPLACE(a0,a2); + CH(i,k,0) = a0+a1; + special_mul(a0-a1,WA(3,i),CH(i,k,4)); + special_mul(a2+a3,WA(1,i),CH(i,k,2)); + special_mul(a2-a3,WA(5,i),CH(i,k,6)); + ROTX90(a6); + PMINPLACE(a4,a6); + special_mul(a4+a5,WA(0,i),CH(i,k,1)); + special_mul(a4-a5,WA(4,i),CH(i,k,5)); + special_mul(a6+a7,WA(2,i),CH(i,k,3)); + special_mul(a6-a7,WA(6,i),CH(i,k,7)); + } + } + } + + +#define POCKETFFT_PREP11(idx) \ + T t1 = CC(idx,0,k), t2, t3, t4, t5, t6, t7, t8, t9, t10, t11; \ + PM (t2,t11,CC(idx,1,k),CC(idx,10,k)); \ + PM (t3,t10,CC(idx,2,k),CC(idx, 9,k)); \ + PM (t4,t9 ,CC(idx,3,k),CC(idx, 8,k)); \ + PM (t5,t8 ,CC(idx,4,k),CC(idx, 7,k)); \ + PM (t6,t7 ,CC(idx,5,k),CC(idx, 6,k)); \ + CH(idx,k,0).r=t1.r+t2.r+t3.r+t4.r+t5.r+t6.r; \ + CH(idx,k,0).i=t1.i+t2.i+t3.i+t4.i+t5.i+t6.i; + +#define POCKETFFT_PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,out1,out2) \ + { \ + T ca = t1 + t2*x1 + t3*x2 + t4*x3 + t5*x4 +t6*x5, \ + cb; \ + cb.i=y1*t11.r y2*t10.r y3*t9.r y4*t8.r y5*t7.r; \ + cb.r=-(y1*t11.i y2*t10.i y3*t9.i y4*t8.i y5*t7.i ); \ + PM(out1,out2,ca,cb); \ + } +#define POCKETFFT_PARTSTEP11a(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5) \ + POCKETFFT_PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,CH(0,k,u1),CH(0,k,u2)) +#define POCKETFFT_PARTSTEP11(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5) \ + { \ + T da,db; \ + POCKETFFT_PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,da,db) \ + special_mul(da,WA(u1-1,i),CH(i,k,u1)); \ + special_mul(db,WA(u2-1,i),CH(i,k,u2)); \ + } + +template void pass11 (size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const cmplx * POCKETFFT_RESTRICT wa) const + { + constexpr T0 tw1r= T0(0.8412535328311811688618116489193677L), + tw1i= (fwd ? -1 : 1) * T0(0.5406408174555975821076359543186917L), + tw2r= T0(0.4154150130018864255292741492296232L), + tw2i= (fwd ? -1 : 1) * T0(0.9096319953545183714117153830790285L), + tw3r= T0(-0.1423148382732851404437926686163697L), + tw3i= (fwd ? -1 : 1) * T0(0.9898214418809327323760920377767188L), + tw4r= T0(-0.6548607339452850640569250724662936L), + tw4i= (fwd ? -1 : 1) * T0(0.7557495743542582837740358439723444L), + tw5r= T0(-0.9594929736144973898903680570663277L), + tw5i= (fwd ? -1 : 1) * T0(0.2817325568414296977114179153466169L); + + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+11*c)]; }; + auto WA = [wa, ido](size_t x, size_t i) + { return wa[i-1+x*(ido-1)]; }; + + if (ido==1) + for (size_t k=0; k void passg (size_t ido, size_t ip, + size_t l1, T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const cmplx * POCKETFFT_RESTRICT wa, + const cmplx * POCKETFFT_RESTRICT csarr) const + { + const size_t cdim=ip; + size_t ipph = (ip+1)/2; + size_t idl1 = ido*l1; + + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + auto CC = [cc,ido,cdim](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+cdim*c)]; }; + auto CX = [cc, ido, l1](size_t a, size_t b, size_t c) -> T& + { return cc[a+ido*(b+l1*c)]; }; + auto CX2 = [cc, idl1](size_t a, size_t b) -> T& + { return cc[a+idl1*b]; }; + auto CH2 = [ch, idl1](size_t a, size_t b) -> const T& + { return ch[a+idl1*b]; }; + + arr> wal(ip); + wal[0] = cmplx(1., 0.); + for (size_t i=1; i(csarr[i].r,fwd ? -csarr[i].i : csarr[i].i); + + for (size_t k=0; kip) iwal-=ip; + cmplx xwal=wal[iwal]; + iwal+=l; if (iwal>ip) iwal-=ip; + cmplx xwal2=wal[iwal]; + for (size_t ik=0; ikip) iwal-=ip; + cmplx xwal=wal[iwal]; + for (size_t ik=0; ik(x1,wa[idij],CX(i,k,j)); + idij=(jc-1)*(ido-1)+i-1; + special_mul(x2,wa[idij],CX(i,k,jc)); + } + } + } + } + +template void pass_all(T c[], T0 fct) const + { + if (length==1) { c[0]*=fct; return; } + size_t l1=1; + arr ch(length); + T *p1=c, *p2=ch.data(); + + for(size_t k1=0; k1 (ido, l1, p1, p2, fact[k1].tw); + else if(ip==8) + pass8(ido, l1, p1, p2, fact[k1].tw); + else if(ip==2) + pass2(ido, l1, p1, p2, fact[k1].tw); + else if(ip==3) + pass3 (ido, l1, p1, p2, fact[k1].tw); + else if(ip==5) + pass5 (ido, l1, p1, p2, fact[k1].tw); + else if(ip==7) + pass7 (ido, l1, p1, p2, fact[k1].tw); + else if(ip==11) + pass11 (ido, l1, p1, p2, fact[k1].tw); + else + { + passg(ido, ip, l1, p1, p2, fact[k1].tw, fact[k1].tws); + std::swap(p1,p2); + } + std::swap(p1,p2); + l1=l2; + } + if (p1!=c) + { + if (fct!=1.) + for (size_t i=0; i void exec(T c[], T0 fct, bool fwd) const + { fwd ? pass_all(c, fct) : pass_all(c, fct); } + + private: + POCKETFFT_NOINLINE void factorize() + { + size_t len=length; + while ((len&7)==0) + { add_factor(8); len>>=3; } + while ((len&3)==0) + { add_factor(4); len>>=2; } + if ((len&1)==0) + { + len>>=1; + // factor 2 should be at the front of the factor list + add_factor(2); + std::swap(fact[0].fct, fact.back().fct); + } + for (size_t divisor=3; divisor*divisor<=len; divisor+=2) + while ((len%divisor)==0) + { + add_factor(divisor); + len/=divisor; + } + if (len>1) add_factor(len); + } + + size_t twsize() const + { + size_t twsize=0, l1=1; + for (size_t k=0; k11) + twsize+=ip; + l1*=ip; + } + return twsize; + } + + void comp_twiddle() + { + sincos_2pibyn twiddle(length); + size_t l1=1; + size_t memofs=0; + for (size_t k=0; k11) + { + fact[k].tws=mem.data()+memofs; + memofs+=ip; + for (size_t j=0; j class rfftp + { + private: + struct fctdata + { + size_t fct; + T0 *tw, *tws; + }; + + size_t length; + arr mem; + std::vector fact; + + void add_factor(size_t factor) + { fact.push_back({factor, nullptr, nullptr}); } + +/* (a+ib) = conj(c+id) * (e+if) */ +template inline void MULPM + (T1 &a, T1 &b, T2 c, T2 d, T3 e, T3 f) const + { a=c*e+d*f; b=c*f-d*e; } + +template void radf2 (size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa) const + { + auto WA = [wa,ido](size_t x, size_t i) { return wa[i+x*(ido-1)]; }; + auto CC = [cc,ido,l1](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+l1*c)]; }; + auto CH = [ch,ido](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+2*c)]; }; + + for (size_t k=0; k void radf3(size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa) const + { + constexpr T0 taur=-0.5, taui=T0(0.8660254037844386467637231707529362L); + + auto WA = [wa,ido](size_t x, size_t i) { return wa[i+x*(ido-1)]; }; + auto CC = [cc,ido,l1](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+l1*c)]; }; + auto CH = [ch,ido](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+3*c)]; }; + + for (size_t k=0; k void radf4(size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa) const + { + constexpr T0 hsqt2=T0(0.707106781186547524400844362104849L); + + auto WA = [wa,ido](size_t x, size_t i) { return wa[i+x*(ido-1)]; }; + auto CC = [cc,ido,l1](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+l1*c)]; }; + auto CH = [ch,ido](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+4*c)]; }; + + for (size_t k=0; k void radf5(size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa) const + { + constexpr T0 tr11= T0(0.3090169943749474241022934171828191L), + ti11= T0(0.9510565162951535721164393333793821L), + tr12= T0(-0.8090169943749474241022934171828191L), + ti12= T0(0.5877852522924731291687059546390728L); + + auto WA = [wa,ido](size_t x, size_t i) { return wa[i+x*(ido-1)]; }; + auto CC = [cc,ido,l1](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+l1*c)]; }; + auto CH = [ch,ido](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+5*c)]; }; + + for (size_t k=0; k void radfg(size_t ido, size_t ip, size_t l1, + T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr) const + { + const size_t cdim=ip; + size_t ipph=(ip+1)/2; + size_t idl1 = ido*l1; + + auto CC = [cc,ido,cdim](size_t a, size_t b, size_t c) -> T& + { return cc[a+ido*(b+cdim*c)]; }; + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> const T& + { return ch[a+ido*(b+l1*c)]; }; + auto C1 = [cc,ido,l1] (size_t a, size_t b, size_t c) -> T& + { return cc[a+ido*(b+l1*c)]; }; + auto C2 = [cc,idl1] (size_t a, size_t b) -> T& + { return cc[a+idl1*b]; }; + auto CH2 = [ch,idl1] (size_t a, size_t b) -> T& + { return ch[a+idl1*b]; }; + + if (ido>1) + { + for (size_t j=1, jc=ip-1; j=ip) iang-=ip; + T0 ar1=csarr[2*iang], ai1=csarr[2*iang+1]; + iang+=l; if (iang>=ip) iang-=ip; + T0 ar2=csarr[2*iang], ai2=csarr[2*iang+1]; + iang+=l; if (iang>=ip) iang-=ip; + T0 ar3=csarr[2*iang], ai3=csarr[2*iang+1]; + iang+=l; if (iang>=ip) iang-=ip; + T0 ar4=csarr[2*iang], ai4=csarr[2*iang+1]; + for (size_t ik=0; ik=ip) iang-=ip; + T0 ar1=csarr[2*iang], ai1=csarr[2*iang+1]; + iang+=l; if (iang>=ip) iang-=ip; + T0 ar2=csarr[2*iang], ai2=csarr[2*iang+1]; + for (size_t ik=0; ik=ip) iang-=ip; + T0 ar=csarr[2*iang], ai=csarr[2*iang+1]; + for (size_t ik=0; ik void radb2(size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa) const + { + auto WA = [wa,ido](size_t x, size_t i) { return wa[i+x*(ido-1)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+2*c)]; }; + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + + for (size_t k=0; k void radb3(size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa) const + { + constexpr T0 taur=-0.5, taui=T0(0.8660254037844386467637231707529362L); + + auto WA = [wa,ido](size_t x, size_t i) { return wa[i+x*(ido-1)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+3*c)]; }; + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + + for (size_t k=0; k void radb4(size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa) const + { + constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L); + + auto WA = [wa,ido](size_t x, size_t i) { return wa[i+x*(ido-1)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+4*c)]; }; + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + + for (size_t k=0; k void radb5(size_t ido, size_t l1, + const T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa) const + { + constexpr T0 tr11= T0(0.3090169943749474241022934171828191L), + ti11= T0(0.9510565162951535721164393333793821L), + tr12= T0(-0.8090169943749474241022934171828191L), + ti12= T0(0.5877852522924731291687059546390728L); + + auto WA = [wa,ido](size_t x, size_t i) { return wa[i+x*(ido-1)]; }; + auto CC = [cc,ido](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+5*c)]; }; + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + + for (size_t k=0; k void radbg(size_t ido, size_t ip, size_t l1, + T * POCKETFFT_RESTRICT cc, T * POCKETFFT_RESTRICT ch, + const T0 * POCKETFFT_RESTRICT wa, const T0 * POCKETFFT_RESTRICT csarr) const + { + const size_t cdim=ip; + size_t ipph=(ip+1)/ 2; + size_t idl1 = ido*l1; + + auto CC = [cc,ido,cdim](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+cdim*c)]; }; + auto CH = [ch,ido,l1](size_t a, size_t b, size_t c) -> T& + { return ch[a+ido*(b+l1*c)]; }; + auto C1 = [cc,ido,l1](size_t a, size_t b, size_t c) -> const T& + { return cc[a+ido*(b+l1*c)]; }; + auto C2 = [cc,idl1](size_t a, size_t b) -> T& + { return cc[a+idl1*b]; }; + auto CH2 = [ch,idl1](size_t a, size_t b) -> T& + { return ch[a+idl1*b]; }; + + for (size_t k=0; kip) iang-=ip; + T0 ar1=csarr[2*iang], ai1=csarr[2*iang+1]; + iang+=l; if(iang>ip) iang-=ip; + T0 ar2=csarr[2*iang], ai2=csarr[2*iang+1]; + iang+=l; if(iang>ip) iang-=ip; + T0 ar3=csarr[2*iang], ai3=csarr[2*iang+1]; + iang+=l; if(iang>ip) iang-=ip; + T0 ar4=csarr[2*iang], ai4=csarr[2*iang+1]; + for (size_t ik=0; ikip) iang-=ip; + T0 ar1=csarr[2*iang], ai1=csarr[2*iang+1]; + iang+=l; if(iang>ip) iang-=ip; + T0 ar2=csarr[2*iang], ai2=csarr[2*iang+1]; + for (size_t ik=0; ikip) iang-=ip; + T0 war=csarr[2*iang], wai=csarr[2*iang+1]; + for (size_t ik=0; ik void copy_and_norm(T *c, T *p1, T0 fct) const + { + if (p1!=c) + { + if (fct!=1.) + for (size_t i=0; i void exec(T c[], T0 fct, bool r2hc) const + { + if (length==1) { c[0]*=fct; return; } + size_t nf=fact.size(); + arr ch(length); + T *p1=c, *p2=ch.data(); + + if (r2hc) + for(size_t k1=0, l1=length; k1>=2; } + if ((len%2)==0) + { + len>>=1; + // factor 2 should be at the front of the factor list + add_factor(2); + std::swap(fact[0].fct, fact.back().fct); + } + for (size_t divisor=3; divisor*divisor<=len; divisor+=2) + while ((len%divisor)==0) + { + add_factor(divisor); + len/=divisor; + } + if (len>1) add_factor(len); + } + + size_t twsize() const + { + size_t twsz=0, l1=1; + for (size_t k=0; k5) twsz+=2*ip; + l1*=ip; + } + return twsz; + } + + void comp_twiddle() + { + sincos_2pibyn twid(length); + size_t l1=1; + T0 *ptr=mem.data(); + for (size_t k=0; k5) // special factors required by *g functions + { + fact[k].tws=ptr; ptr+=2*ip; + fact[k].tws[0] = 1.; + fact[k].tws[1] = 0.; + for (size_t i=2, ic=2*ip-2; i<=ic; i+=2, ic-=2) + { + fact[k].tws[i ] = twid[i/2*(length/ip)].r; + fact[k].tws[i+1] = twid[i/2*(length/ip)].i; + fact[k].tws[ic] = twid[i/2*(length/ip)].r; + fact[k].tws[ic+1] = -twid[i/2*(length/ip)].i; + } + } + l1*=ip; + } + } + + public: + POCKETFFT_NOINLINE explicit rfftp(size_t length_) + : length(length_) + { + if (length==0) throw std::runtime_error("zero-length FFT requested"); + if (length==1) return; + factorize(); + mem.resize(twsize()); + comp_twiddle(); + } +}; + +// +// complex Bluestein transforms +// + +template class fftblue + { + private: + size_t n, n2; + cfftp plan; + arr> mem; + cmplx *bk, *bkf; + + template void fft(cmplx c[], T0 fct) const + { + arr> akf(n2); + + /* initialize a_k and FFT it */ + for (size_t m=0; m(c[m],bk[m],akf[m]); + auto zero = akf[0]*T0(0); + for (size_t m=n; m(bkf[0]); + for (size_t m=1; m<(n2+1)/2; ++m) + { + akf[m] = akf[m].template special_mul(bkf[m]); + akf[n2-m] = akf[n2-m].template special_mul(bkf[m]); + } + if ((n2&1)==0) + akf[n2/2] = akf[n2/2].template special_mul(bkf[n2/2]); + + /* inverse FFT */ + plan.exec (akf.data(),1.,false); + + /* multiply by b_k */ + for (size_t m=0; m(bk[m])*fct; + } + + public: + POCKETFFT_NOINLINE explicit fftblue(size_t length) + : n(length), n2(util::good_size_cmplx(n*2-1)), plan(n2), mem(n+n2/2+1), + bk(mem.data()), bkf(mem.data()+n) + { + /* initialize b_k */ + sincos_2pibyn tmp(2*n); + bk[0].Set(1, 0); + + size_t coeff=0; + for (size_t m=1; m=2*n) coeff-=2*n; + bk[m] = tmp[coeff]; + } + + /* initialize the zero-padded, Fourier transformed b_k. Add normalisation. */ + arr> tbkf(n2); + T0 xn2 = T0(1)/T0(n2); + tbkf[0] = bk[0]*xn2; + for (size_t m=1; m void exec(cmplx c[], T0 fct, bool fwd) const + { fwd ? fft(c,fct) : fft(c,fct); } + + template void exec_r(T c[], T0 fct, bool fwd) + { + arr> tmp(n); + if (fwd) + { + auto zero = T0(0)*c[0]; + for (size_t m=0; m(tmp.data(),fct); + c[0] = tmp[0].r; + std::copy_n (&tmp[1].r, n-1, &c[1]); + } + else + { + tmp[0].Set(c[0],c[0]*0); + std::copy_n (c+1, n-1, &tmp[1].r); + if ((n&1)==0) tmp[n/2].i=T0(0)*c[0]; + for (size_t m=1; 2*m(tmp.data(),fct); + for (size_t m=0; m class pocketfft_c + { + private: + std::unique_ptr> packplan; + std::unique_ptr> blueplan; + size_t len; + + public: + POCKETFFT_NOINLINE explicit pocketfft_c(size_t length) + : len(length) + { + if (length==0) throw std::runtime_error("zero-length FFT requested"); + size_t tmp = (length<50) ? 0 : util::largest_prime_factor(length); + if (tmp*tmp <= length) + { + packplan=std::unique_ptr>(new cfftp(length)); + return; + } + double comp1 = util::cost_guess(length); + double comp2 = 2*util::cost_guess(util::good_size_cmplx(2*length-1)); + comp2*=1.5; /* fudge factor that appears to give good overall performance */ + if (comp2>(new fftblue(length)); + else + packplan=std::unique_ptr>(new cfftp(length)); + } + + template POCKETFFT_NOINLINE void exec(cmplx c[], T0 fct, bool fwd) const + { packplan ? packplan->exec(c,fct,fwd) : blueplan->exec(c,fct,fwd); } + + size_t length() const { return len; } + }; + +// +// flexible (FFTPACK/Bluestein) real-valued 1D transform +// + +template class pocketfft_r + { + private: + std::unique_ptr> packplan; + std::unique_ptr> blueplan; + size_t len; + + public: + POCKETFFT_NOINLINE explicit pocketfft_r(size_t length) + : len(length) + { + if (length==0) throw std::runtime_error("zero-length FFT requested"); + size_t tmp = (length<50) ? 0 : util::largest_prime_factor(length); + if (tmp*tmp <= length) + { + packplan=std::unique_ptr>(new rfftp(length)); + return; + } + double comp1 = 0.5*util::cost_guess(length); + double comp2 = 2*util::cost_guess(util::good_size_cmplx(2*length-1)); + comp2*=1.5; /* fudge factor that appears to give good overall performance */ + if (comp2>(new fftblue(length)); + else + packplan=std::unique_ptr>(new rfftp(length)); + } + + template POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool fwd) const + { packplan ? packplan->exec(c,fct,fwd) : blueplan->exec_r(c,fct,fwd); } + + size_t length() const { return len; } + }; + + +// +// sine/cosine transforms +// + +template class T_dct1 + { + private: + pocketfft_r fftplan; + + public: + POCKETFFT_NOINLINE explicit T_dct1(size_t length) + : fftplan(2*(length-1)) {} + + template POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho, + int /*type*/, bool /*cosine*/) const + { + constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L); + size_t N=fftplan.length(), n=N/2+1; + if (ortho) + { c[0]*=sqrt2; c[n-1]*=sqrt2; } + arr tmp(N); + tmp[0] = c[0]; + for (size_t i=1; i class T_dst1 + { + private: + pocketfft_r fftplan; + + public: + POCKETFFT_NOINLINE explicit T_dst1(size_t length) + : fftplan(2*(length+1)) {} + + template POCKETFFT_NOINLINE void exec(T c[], T0 fct, + bool /*ortho*/, int /*type*/, bool /*cosine*/) const + { + size_t N=fftplan.length(), n=N/2-1; + arr tmp(N); + tmp[0] = tmp[n+1] = c[0]*0; + for (size_t i=0; i class T_dcst23 + { + private: + pocketfft_r fftplan; + std::vector twiddle; + + public: + POCKETFFT_NOINLINE explicit T_dcst23(size_t length) + : fftplan(length), twiddle(length) + { + sincos_2pibyn tw(4*length); + for (size_t i=0; i POCKETFFT_NOINLINE void exec(T c[], T0 fct, bool ortho, + int type, bool cosine) const + { + constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L); + size_t N=length(); + size_t NS2 = (N+1)/2; + if (type==2) + { + if (!cosine) + for (size_t k=1; k class T_dcst4 + { + private: + size_t N; + std::unique_ptr> fft; + std::unique_ptr> rfft; + arr> C2; + + public: + POCKETFFT_NOINLINE explicit T_dcst4(size_t length) + : N(length), + fft((N&1) ? nullptr : new pocketfft_c(N/2)), + rfft((N&1)? new pocketfft_r(N) : nullptr), + C2((N&1) ? 0 : N/2) + { + if ((N&1)==0) + { + sincos_2pibyn tw(16*N); + for (size_t i=0; i POCKETFFT_NOINLINE void exec(T c[], T0 fct, + bool /*ortho*/, int /*type*/, bool cosine) const + { + size_t n2 = N/2; + if (!cosine) + for (size_t k=0, kc=N-1; k y(N); + { + size_t i=0, m=n2; + for (; mexec(y.data(), fct, true); + { + auto SGN = [](size_t i) + { + constexpr T0 sqrt2=T0(1.414213562373095048801688724209698L); + return (i&2) ? -sqrt2 : sqrt2; + }; + c[n2] = y[0]*SGN(n2+1); + size_t i=0, i1=1, k=1; + for (; k> y(n2); + for(size_t i=0; iexec(y.data(), fct, true); + for(size_t i=0, ic=n2-1; i std::shared_ptr get_plan(size_t length) + { +#if POCKETFFT_CACHE_SIZE==0 + return std::make_shared(length); +#else + constexpr size_t nmax=POCKETFFT_CACHE_SIZE; + static std::array, nmax> cache; + static std::array last_access{{0}}; + static size_t access_counter = 0; + static std::mutex mut; + + auto find_in_cache = [&]() -> std::shared_ptr + { + for (size_t i=0; ilength()==length)) + { + // no need to update if this is already the most recent entry + if (last_access[i]!=access_counter) + { + last_access[i] = ++access_counter; + // Guard against overflow + if (access_counter == 0) + last_access.fill(0); + } + return cache[i]; + } + + return nullptr; + }; + + { + std::lock_guard lock(mut); + auto p = find_in_cache(); + if (p) return p; + } + auto plan = std::make_shared(length); + { + std::lock_guard lock(mut); + auto p = find_in_cache(); + if (p) return p; + + size_t lru = 0; + for (size_t i=1; i class cndarr: public arr_info + { + protected: + const char *d; + + public: + cndarr(const void *data_, const shape_t &shape_, const stride_t &stride_) + : arr_info(shape_, stride_), + d(reinterpret_cast(data_)) {} + const T &operator[](ptrdiff_t ofs) const + { return *reinterpret_cast(d+ofs); } + }; + +template class ndarr: public cndarr + { + public: + ndarr(void *data_, const shape_t &shape_, const stride_t &stride_) + : cndarr::cndarr(const_cast(data_), shape_, stride_) + {} + T &operator[](ptrdiff_t ofs) + { return *reinterpret_cast(const_cast(cndarr::d+ofs)); } + }; + +template class multi_iter + { + private: + shape_t pos; + const arr_info &iarr, &oarr; + ptrdiff_t p_ii, p_i[N], str_i, p_oi, p_o[N], str_o; + size_t idim, rem; + + void advance_i() + { + for (int i_=int(pos.size())-1; i_>=0; --i_) + { + auto i = size_t(i_); + if (i==idim) continue; + p_ii += iarr.stride(i); + p_oi += oarr.stride(i); + if (++pos[i] < iarr.shape(i)) + return; + pos[i] = 0; + p_ii -= ptrdiff_t(iarr.shape(i))*iarr.stride(i); + p_oi -= ptrdiff_t(oarr.shape(i))*oarr.stride(i); + } + } + + public: + multi_iter(const arr_info &iarr_, const arr_info &oarr_, size_t idim_) + : pos(iarr_.ndim(), 0), iarr(iarr_), oarr(oarr_), p_ii(0), + str_i(iarr.stride(idim_)), p_oi(0), str_o(oarr.stride(idim_)), + idim(idim_), rem(iarr.size()/iarr.shape(idim)) + { + auto nshares = threading::num_threads(); + if (nshares==1) return; + if (nshares==0) throw std::runtime_error("can't run with zero threads"); + auto myshare = threading::thread_id(); + if (myshare>=nshares) throw std::runtime_error("impossible share requested"); + size_t nbase = rem/nshares; + size_t additional = rem%nshares; + size_t lo = myshare*nbase + ((myshare=0; --i_) + { + auto i = size_t(i_); + p += arr.stride(i); + if (++pos[i] < arr.shape(i)) + return; + pos[i] = 0; + p -= ptrdiff_t(arr.shape(i))*arr.stride(i); + } + } + ptrdiff_t ofs() const { return p; } + size_t remaining() const { return rem; } + }; + +class rev_iter + { + private: + shape_t pos; + const arr_info &arr; + std::vector rev_axis; + std::vector rev_jump; + size_t last_axis, last_size; + shape_t shp; + ptrdiff_t p, rp; + size_t rem; + + public: + rev_iter(const arr_info &arr_, const shape_t &axes) + : pos(arr_.ndim(), 0), arr(arr_), rev_axis(arr_.ndim(), 0), + rev_jump(arr_.ndim(), 1), p(0), rp(0) + { + for (auto ax: axes) + rev_axis[ax]=1; + last_axis = axes.back(); + last_size = arr.shape(last_axis)/2 + 1; + shp = arr.shape(); + shp[last_axis] = last_size; + rem=1; + for (auto i: shp) + rem *= i; + } + void advance() + { + --rem; + for (int i_=int(pos.size())-1; i_>=0; --i_) + { + auto i = size_t(i_); + p += arr.stride(i); + if (!rev_axis[i]) + rp += arr.stride(i); + else + { + rp -= arr.stride(i); + if (rev_jump[i]) + { + rp += ptrdiff_t(arr.shape(i))*arr.stride(i); + rev_jump[i] = 0; + } + } + if (++pos[i] < shp[i]) + return; + pos[i] = 0; + p -= ptrdiff_t(shp[i])*arr.stride(i); + if (rev_axis[i]) + { + rp -= ptrdiff_t(arr.shape(i)-shp[i])*arr.stride(i); + rev_jump[i] = 1; + } + else + rp -= ptrdiff_t(shp[i])*arr.stride(i); + } + } + ptrdiff_t ofs() const { return p; } + ptrdiff_t rev_ofs() const { return rp; } + size_t remaining() const { return rem; } + }; + +template struct VTYPE {}; +template using vtype_t = typename VTYPE::type; + +#ifndef POCKETFFT_NO_VECTORS +template<> struct VTYPE + { + using type = float __attribute__ ((vector_size (VLEN::val*sizeof(float)))); + }; +template<> struct VTYPE + { + using type = double __attribute__ ((vector_size (VLEN::val*sizeof(double)))); + }; +template<> struct VTYPE + { + using type = long double __attribute__ ((vector_size (VLEN::val*sizeof(long double)))); + }; +#endif + +template arr alloc_tmp(const shape_t &shape, + size_t axsize, size_t elemsize) + { + auto othersize = util::prod(shape)/axsize; + auto tmpsize = axsize*((othersize>=VLEN::val) ? VLEN::val : 1); + return arr(tmpsize*elemsize); + } +template arr alloc_tmp(const shape_t &shape, + const shape_t &axes, size_t elemsize) + { + size_t fullsize=util::prod(shape); + size_t tmpsize=0; + for (size_t i=0; i=VLEN::val) ? VLEN::val : 1); + if (sz>tmpsize) tmpsize=sz; + } + return arr(tmpsize*elemsize); + } + +template void copy_input(const multi_iter &it, + const cndarr> &src, cmplx> *POCKETFFT_RESTRICT dst) + { + for (size_t i=0; i void copy_input(const multi_iter &it, + const cndarr &src, vtype_t *POCKETFFT_RESTRICT dst) + { + for (size_t i=0; i void copy_input(const multi_iter &it, + const cndarr &src, T *POCKETFFT_RESTRICT dst) + { + if (dst == &src[it.iofs(0)]) return; // in-place + for (size_t i=0; i void copy_output(const multi_iter &it, + const cmplx> *POCKETFFT_RESTRICT src, ndarr> &dst) + { + for (size_t i=0; i void copy_output(const multi_iter &it, + const vtype_t *POCKETFFT_RESTRICT src, ndarr &dst) + { + for (size_t i=0; i void copy_output(const multi_iter &it, + const T *POCKETFFT_RESTRICT src, ndarr &dst) + { + if (src == &dst[it.oofs(0)]) return; // in-place + for (size_t i=0; i struct add_vec { using type = vtype_t; }; +template struct add_vec> + { using type = cmplx>; }; +template using add_vec_t = typename add_vec::type; + +template +POCKETFFT_NOINLINE void general_nd(const cndarr &in, ndarr &out, + const shape_t &axes, T0 fct, size_t nthreads, const Exec & exec, + const bool allow_inplace=true) + { + std::shared_ptr plan; + + for (size_t iax=0; iaxlength())) + plan = get_plan(len); + + threading::thread_map( + util::thread_count(nthreads, in.shape(), axes[iax], VLEN::val), + [&] { + constexpr auto vlen = VLEN::val; + auto storage = alloc_tmp(in.shape(), len, sizeof(T)); + const auto &tin(iax==0? in : out); + multi_iter it(tin, out, axes[iax]); +#ifndef POCKETFFT_NO_VECTORS + if (vlen>1) + while (it.remaining()>=vlen) + { + it.advance(vlen); + auto tdatav = reinterpret_cast *>(storage.data()); + exec(it, tin, out, tdatav, *plan, fct); + } +#endif + while (it.remaining()>0) + { + it.advance(1); + auto buf = allow_inplace && it.stride_out() == sizeof(T) ? + &out[it.oofs(0)] : reinterpret_cast(storage.data()); + exec(it, tin, out, buf, *plan, fct); + } + }); // end of parallel region + fct = T0(1); // factor has been applied, use 1 for remaining axes + } + } + +struct ExecC2C + { + bool forward; + + template void operator () ( + const multi_iter &it, const cndarr> &in, + ndarr> &out, T * buf, const pocketfft_c &plan, T0 fct) const + { + copy_input(it, in, buf); + plan.exec(buf, fct, forward); + copy_output(it, buf, out); + } + }; + +template void copy_hartley(const multi_iter &it, + const vtype_t *POCKETFFT_RESTRICT src, ndarr &dst) + { + for (size_t j=0; j void copy_hartley(const multi_iter &it, + const T *POCKETFFT_RESTRICT src, ndarr &dst) + { + dst[it.oofs(0)] = src[0]; + size_t i=1, i1=1, i2=it.length_out()-1; + for (i=1; i void copy_FHT(const multi_iter &it, + const vtype_t *POCKETFFT_RESTRICT src, ndarr &dst) + { + for (size_t j=0; j void copy_FHT(const multi_iter &it, + const T *POCKETFFT_RESTRICT src, ndarr &dst) + { + dst[it.oofs(0)] = src[0]; + size_t i=1, i1=1, i2=it.length_out()-1; + for (i=1; i void operator () ( + const multi_iter &it, const cndarr &in, ndarr &out, + T * buf, const pocketfft_r &plan, T0 fct) const + { + copy_input(it, in, buf); + plan.exec(buf, fct, true); + copy_hartley(it, buf, out); + } + }; +struct ExecFHT + { + template void operator () ( + const multi_iter &it, const cndarr &in, ndarr &out, + T * buf, const pocketfft_r &plan, T0 fct) const + { + copy_input(it, in, buf); + plan.exec(buf, fct, true); + copy_FHT(it, buf, out); + } + }; + +struct ExecDcst + { + bool ortho; + int type; + bool cosine; + + template + void operator () (const multi_iter &it, const cndarr &in, + ndarr &out, T * buf, const Tplan &plan, T0 fct) const + { + copy_input(it, in, buf); + plan.exec(buf, fct, ortho, type, cosine); + copy_output(it, buf, out); + } + }; + +template POCKETFFT_NOINLINE void general_r2c( + const cndarr &in, ndarr> &out, size_t axis, bool forward, T fct, + size_t nthreads) + { + auto plan = get_plan>(in.shape(axis)); + size_t len=in.shape(axis); + threading::thread_map( + util::thread_count(nthreads, in.shape(), axis, VLEN::val), + [&] { + constexpr auto vlen = VLEN::val; + auto storage = alloc_tmp(in.shape(), len, sizeof(T)); + multi_iter it(in, out, axis); +#ifndef POCKETFFT_NO_VECTORS + if (vlen>1) + while (it.remaining()>=vlen) + { + it.advance(vlen); + auto tdatav = reinterpret_cast *>(storage.data()); + copy_input(it, in, tdatav); + plan->exec(tdatav, fct, true); + for (size_t j=0; j0) + { + it.advance(1); + auto tdata = reinterpret_cast(storage.data()); + copy_input(it, in, tdata); + plan->exec(tdata, fct, true); + out[it.oofs(0)].Set(tdata[0]); + size_t i=1, ii=1; + if (forward) + for (; i POCKETFFT_NOINLINE void general_c2r( + const cndarr> &in, ndarr &out, size_t axis, bool forward, T fct, + size_t nthreads) + { + auto plan = get_plan>(out.shape(axis)); + size_t len=out.shape(axis); + threading::thread_map( + util::thread_count(nthreads, in.shape(), axis, VLEN::val), + [&] { + constexpr auto vlen = VLEN::val; + auto storage = alloc_tmp(out.shape(), len, sizeof(T)); + multi_iter it(in, out, axis); +#ifndef POCKETFFT_NO_VECTORS + if (vlen>1) + while (it.remaining()>=vlen) + { + it.advance(vlen); + auto tdatav = reinterpret_cast *>(storage.data()); + for (size_t j=0; jexec(tdatav, fct, false); + copy_output(it, tdatav, out); + } +#endif + while (it.remaining()>0) + { + it.advance(1); + auto tdata = reinterpret_cast(storage.data()); + tdata[0]=in[it.iofs(0)].r; + { + size_t i=1, ii=1; + if (forward) + for (; iexec(tdata, fct, false); + copy_output(it, tdata, out); + } + }); // end of parallel region + } + +struct ExecR2R + { + bool r2h, forward; + + template void operator () ( + const multi_iter &it, const cndarr &in, ndarr &out, T * buf, + const pocketfft_r &plan, T0 fct) const + { + copy_input(it, in, buf); + if ((!r2h) && forward) + for (size_t i=2; i void c2c(const shape_t &shape, const stride_t &stride_in, + const stride_t &stride_out, const shape_t &axes, bool forward, + const std::complex *data_in, std::complex *data_out, T fct, + size_t nthreads=1) + { + if (util::prod(shape)==0) return; + util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes); + cndarr> ain(data_in, shape, stride_in); + ndarr> aout(data_out, shape, stride_out); + general_nd>(ain, aout, axes, fct, nthreads, ExecC2C{forward}); + } + +template void dct(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes, + int type, const T *data_in, T *data_out, T fct, bool ortho, size_t nthreads=1) + { + if ((type<1) || (type>4)) throw std::invalid_argument("invalid DCT type"); + if (util::prod(shape)==0) return; + util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes); + cndarr ain(data_in, shape, stride_in); + ndarr aout(data_out, shape, stride_out); + const ExecDcst exec{ortho, type, true}; + if (type==1) + general_nd>(ain, aout, axes, fct, nthreads, exec); + else if (type==4) + general_nd>(ain, aout, axes, fct, nthreads, exec); + else + general_nd>(ain, aout, axes, fct, nthreads, exec); + } + +template void dst(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes, + int type, const T *data_in, T *data_out, T fct, bool ortho, size_t nthreads=1) + { + if ((type<1) || (type>4)) throw std::invalid_argument("invalid DST type"); + if (util::prod(shape)==0) return; + util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes); + cndarr ain(data_in, shape, stride_in); + ndarr aout(data_out, shape, stride_out); + const ExecDcst exec{ortho, type, false}; + if (type==1) + general_nd>(ain, aout, axes, fct, nthreads, exec); + else if (type==4) + general_nd>(ain, aout, axes, fct, nthreads, exec); + else + general_nd>(ain, aout, axes, fct, nthreads, exec); + } + +template void r2c(const shape_t &shape_in, + const stride_t &stride_in, const stride_t &stride_out, size_t axis, + bool forward, const T *data_in, std::complex *data_out, T fct, + size_t nthreads=1) + { + if (util::prod(shape_in)==0) return; + util::sanity_check(shape_in, stride_in, stride_out, false, axis); + cndarr ain(data_in, shape_in, stride_in); + shape_t shape_out(shape_in); + shape_out[axis] = shape_in[axis]/2 + 1; + ndarr> aout(data_out, shape_out, stride_out); + general_r2c(ain, aout, axis, forward, fct, nthreads); + } + +template void r2c(const shape_t &shape_in, + const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes, + bool forward, const T *data_in, std::complex *data_out, T fct, + size_t nthreads=1) + { + if (util::prod(shape_in)==0) return; + util::sanity_check(shape_in, stride_in, stride_out, false, axes); + r2c(shape_in, stride_in, stride_out, axes.back(), forward, data_in, data_out, + fct, nthreads); + if (axes.size()==1) return; + + shape_t shape_out(shape_in); + shape_out[axes.back()] = shape_in[axes.back()]/2 + 1; + auto newaxes = shape_t{axes.begin(), --axes.end()}; + c2c(shape_out, stride_out, stride_out, newaxes, forward, data_out, data_out, + T(1), nthreads); + } + +template void c2r(const shape_t &shape_out, + const stride_t &stride_in, const stride_t &stride_out, size_t axis, + bool forward, const std::complex *data_in, T *data_out, T fct, + size_t nthreads=1) + { + if (util::prod(shape_out)==0) return; + util::sanity_check(shape_out, stride_in, stride_out, false, axis); + shape_t shape_in(shape_out); + shape_in[axis] = shape_out[axis]/2 + 1; + cndarr> ain(data_in, shape_in, stride_in); + ndarr aout(data_out, shape_out, stride_out); + general_c2r(ain, aout, axis, forward, fct, nthreads); + } + +template void c2r(const shape_t &shape_out, + const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes, + bool forward, const std::complex *data_in, T *data_out, T fct, + size_t nthreads=1) + { + if (util::prod(shape_out)==0) return; + if (axes.size()==1) + return c2r(shape_out, stride_in, stride_out, axes[0], forward, + data_in, data_out, fct, nthreads); + util::sanity_check(shape_out, stride_in, stride_out, false, axes); + auto shape_in = shape_out; + shape_in[axes.back()] = shape_out[axes.back()]/2 + 1; + auto nval = util::prod(shape_in); + stride_t stride_inter(shape_in.size()); + stride_inter.back() = sizeof(cmplx); + for (int i=int(shape_in.size())-2; i>=0; --i) + stride_inter[size_t(i)] = + stride_inter[size_t(i+1)]*ptrdiff_t(shape_in[size_t(i+1)]); + arr> tmp(nval); + auto newaxes = shape_t{axes.begin(), --axes.end()}; + c2c(shape_in, stride_in, stride_inter, newaxes, forward, data_in, tmp.data(), + T(1), nthreads); + c2r(shape_out, stride_inter, stride_out, axes.back(), forward, + tmp.data(), data_out, fct, nthreads); + } + +template void r2r_fftpack(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes, + bool real2hermitian, bool forward, const T *data_in, T *data_out, T fct, + size_t nthreads=1) + { + if (util::prod(shape)==0) return; + util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes); + cndarr ain(data_in, shape, stride_in); + ndarr aout(data_out, shape, stride_out); + general_nd>(ain, aout, axes, fct, nthreads, + ExecR2R{real2hermitian, forward}); + } + +template void r2r_separable_hartley(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes, + const T *data_in, T *data_out, T fct, size_t nthreads=1) + { + if (util::prod(shape)==0) return; + util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes); + cndarr ain(data_in, shape, stride_in); + ndarr aout(data_out, shape, stride_out); + general_nd>(ain, aout, axes, fct, nthreads, ExecHartley{}, + false); + } + +template void r2r_genuine_hartley(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes, + const T *data_in, T *data_out, T fct, size_t nthreads=1) + { + if (util::prod(shape)==0) return; + if (axes.size()==1) + return r2r_separable_hartley(shape, stride_in, stride_out, axes, data_in, + data_out, fct, nthreads); + util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes); + shape_t tshp(shape); + tshp[axes.back()] = tshp[axes.back()]/2+1; + arr> tdata(util::prod(tshp)); + stride_t tstride(shape.size()); + tstride.back()=sizeof(std::complex); + for (size_t i=tstride.size()-1; i>0; --i) + tstride[i-1]=tstride[i]*ptrdiff_t(tshp[i]); + r2c(shape, stride_in, tstride, axes, true, data_in, tdata.data(), fct, nthreads); + cndarr> atmp(tdata.data(), tshp, tstride); + ndarr aout(data_out, shape, stride_out); + simple_iter iin(atmp); + rev_iter iout(aout, axes); + while(iin.remaining()>0) + { + auto v = atmp[iin.ofs()]; + aout[iout.ofs()] = v.r+v.i; + aout[iout.rev_ofs()] = v.r-v.i; + iin.advance(); iout.advance(); + } + } + +template void r2r_separable_fht(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes, + const T *data_in, T *data_out, T fct, size_t nthreads=1) + { + if (util::prod(shape)==0) return; + util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes); + cndarr ain(data_in, shape, stride_in); + ndarr aout(data_out, shape, stride_out); + general_nd>(ain, aout, axes, fct, nthreads, ExecFHT{}, + false); + } + +template void r2r_genuine_fht(const shape_t &shape, + const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes, + const T *data_in, T *data_out, T fct, size_t nthreads=1) + { + if (util::prod(shape)==0) return; + if (axes.size()==1) + return r2r_separable_fht(shape, stride_in, stride_out, axes, data_in, + data_out, fct, nthreads); + util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes); + shape_t tshp(shape); + tshp[axes.back()] = tshp[axes.back()]/2+1; + arr> tdata(util::prod(tshp)); + stride_t tstride(shape.size()); + tstride.back()=sizeof(std::complex); + for (size_t i=tstride.size()-1; i>0; --i) + tstride[i-1]=tstride[i]*ptrdiff_t(tshp[i]); + r2c(shape, stride_in, tstride, axes, true, data_in, tdata.data(), fct, nthreads); + cndarr> atmp(tdata.data(), tshp, tstride); + ndarr aout(data_out, shape, stride_out); + simple_iter iin(atmp); + rev_iter iout(aout, axes); + while(iin.remaining()>0) + { + auto v = atmp[iin.ofs()]; + aout[iout.ofs()] = v.r-v.i; + aout[iout.rev_ofs()] = v.r+v.i; + iin.advance(); iout.advance(); + } + } + +} // namespace detail + +using detail::FORWARD; +using detail::BACKWARD; +using detail::shape_t; +using detail::stride_t; +using detail::c2c; +using detail::c2r; +using detail::r2c; +using detail::r2r_fftpack; +using detail::r2r_separable_hartley; +using detail::r2r_genuine_hartley; +using detail::r2r_separable_fht; +using detail::r2r_genuine_fht; +using detail::dct; +using detail::dst; + +} // namespace pocketfft + +#undef POCKETFFT_NOINLINE +#undef POCKETFFT_RESTRICT + +#endif // POCKETFFT_HDRONLY_H diff --git a/requirements.txt b/requirements.txt index e525707..6522b48 100644 --- a/requirements.txt +++ b/requirements.txt @@ -12,3 +12,8 @@ pyqtgraph>=0.13.7 rpi-hardware-pwm>=0.2.2,<1 fastapi>=0.110 # embedded web UI (browser view + controls), served in-process uvicorn[standard]>=0.29 # ASGI server for the web UI, run on a background thread +# CUDA headers + libnvrtc for building the vendored gpr_cpp library (3D SAR). +# x86_64 dev machines only: on Jetson (aarch64) CUDA comes with JetPack at +# /usr/local/cuda and build_cmake.sh falls back to it automatically. +nvidia-cuda-nvrtc==13.3.33; platform_machine == "x86_64" # nvrtc.h + libnvrtc +nvidia-cuda-runtime==13.3.29; platform_machine == "x86_64" # cuda.h (Driver API)