diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/.clang-format" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/.clang-format" new file mode 100644 index 00000000..96332f09 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/.clang-format" @@ -0,0 +1,3 @@ +BasedOnStyle: LLVM +IndentWidth: 4 +ContinuationIndentWidth: 4 diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/.gitignore" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/.gitignore" new file mode 100644 index 00000000..4262a39f --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/.gitignore" @@ -0,0 +1,60 @@ +# CMake build directories (including build-* / cmake-build-* at any level) +/build/ +build-*/ +/out/ +cmake-build-*/ + +# CMake generated files +CMakeCache.txt +CMakeFiles/ +cmake_install.cmake +CTestTestfile.cmake +Makefile +compile_commands.json +install_manifest.txt + +# Compiled binaries and CUDA artifacts +*.a +*.dll +*.exe +*.exp +*.ilk +*.lib +*.obj +*.o +*.pdb +*.ptx +*.cubin +*.fatbin + +# IDE and editor settings +.vs/ +.vscode/ +.idea/ +*.suo +*.user +*.userosscache +*.sln.docstates + +# Generated experiment and profiling output +/results/* +!/results/.gitkeep +*.nsys-rep +*.ncu-rep +*.sqlite + +# Temporary files +*.log +*.tmp +*.py[cod] +__pycache__/ +*.swp +*~ +.DS_Store +Thumbs.db + +# Agent-managed isolated workspaces +/.worktrees/ +/.superpowers/ + +/.workbuddy \ No newline at end of file diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/CMakeLists.txt" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/CMakeLists.txt" new file mode 100644 index 00000000..2f6ff695 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/CMakeLists.txt" @@ -0,0 +1,107 @@ +# 配置 CUDA 定价器的构建目标、可选分析能力和自动化测试。 +cmake_minimum_required(VERSION 3.24) + +project(cuda_pricer LANGUAGES CXX CUDA) + +# CUDA Toolkit 提供运行时、cuRAND 和 CUB。 +find_package(CUDAToolkit REQUIRED) + +if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES OR CMAKE_CUDA_ARCHITECTURES STREQUAL "") + message(FATAL_ERROR "Set CMAKE_CUDA_ARCHITECTURES to the target GPU architecture.") +endif() + +# 将分析能力和性能模式显式设为独立开关。 +option(PRICER_ENABLE_NVTX "Enable NVTX performance markers" OFF) +option(PRICER_ENABLE_FAST_MATH "Enable CUDA fast math" OFF) + +# 固定第三方库版本，保证干净环境能重复配置。 +include(FetchContent) + +FetchContent_Declare( + nlohmann_json + GIT_REPOSITORY https://github.com/nlohmann/json.git + GIT_TAG v3.11.3 + GIT_SHALLOW TRUE +) +FetchContent_MakeAvailable(nlohmann_json) + +# 与设备无关的金融模型、配置、统计和输出代码。 +add_library(pricer_core STATIC + src/black_scholes.cpp + src/config.cpp + src/cpu_pricer.cpp + src/measurement.cpp + src/output.cpp + src/payoff.cpp + src/statistics.cpp + src/version.cpp +) +add_library(pricer::core ALIAS pricer_core) +target_compile_features(pricer_core PUBLIC cxx_std_17) +target_include_directories(pricer_core + PUBLIC + $ +) +target_link_libraries(pricer_core PUBLIC nlohmann_json::nlohmann_json) + +# 单独构建设备端后端，避免 CPU 测试依赖 CUDA 内核实现细节。 +add_library(pricer_cuda_backend STATIC src/cuda_backend.cu) +add_library(pricer::cuda_backend ALIAS pricer_cuda_backend) +target_compile_features(pricer_cuda_backend PUBLIC cxx_std_17) +target_include_directories(pricer_cuda_backend + PUBLIC + $ +) +target_link_libraries(pricer_cuda_backend PUBLIC pricer::core CUDA::cudart) +if(MSVC) + target_compile_options(pricer_cuda_backend PRIVATE + $<$:-Xcompiler=/Zc:preprocessor> + ) +endif() +set_target_properties(pricer_cuda_backend PROPERTIES CUDA_SEPARABLE_COMPILATION ON) + +# 快速数学只用于明确请求的性能实验，不能作为正确性默认值。 +if(PRICER_ENABLE_FAST_MATH) + target_compile_options(pricer_cuda_backend PRIVATE + $<$:--use_fast_math> + ) +endif() + +# NVTX 标记仅为剖析提供时间线语义，常规构建不依赖它。 +if(PRICER_ENABLE_NVTX) + set(PRICER_NVTX_INCLUDE_DIR "" CACHE PATH + "Directory containing nvtx3/nvToolsExt.h") + if(NOT PRICER_NVTX_INCLUDE_DIR) + unset(PRICER_NVTX_INCLUDE_DIR CACHE) + find_path(PRICER_NVTX_INCLUDE_DIR nvtx3/nvToolsExt.h + HINTS ${CUDAToolkit_INCLUDE_DIRS} "$ENV{NVTX_PATH}" + PATH_SUFFIXES include) + endif() + if(NOT PRICER_NVTX_INCLUDE_DIR) + message(FATAL_ERROR + "PRICER_ENABLE_NVTX requires nvtx3/nvToolsExt.h; set PRICER_NVTX_INCLUDE_DIR.") + endif() + target_compile_definitions(pricer_cuda_backend PRIVATE PRICER_ENABLE_NVTX=1) + target_include_directories(pricer_cuda_backend PRIVATE + "${PRICER_NVTX_INCLUDE_DIR}") + target_link_libraries(pricer_cuda_backend PRIVATE CUDA::nvtx3) +endif() + +# CLI 负责把配置、定价后端和文件输出串成完整一次运行。 +add_executable(pricer_cli src/main.cpp) +target_compile_features(pricer_cli PRIVATE cxx_std_17) +target_link_libraries(pricer_cli PRIVATE pricer::core pricer::cuda_backend) +execute_process( + COMMAND git rev-parse --short HEAD + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + OUTPUT_VARIABLE PRICER_GIT_COMMIT + OUTPUT_STRIP_TRAILING_WHITESPACE + ERROR_QUIET +) +if(PRICER_GIT_COMMIT STREQUAL "") + set(PRICER_GIT_COMMIT "unknown") +endif() +target_compile_definitions(pricer_cli PRIVATE + PRICER_BUILD_TYPE="$" + PRICER_GIT_COMMIT="${PRICER_GIT_COMMIT}" +) diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/README.md" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/README.md" new file mode 100644 index 00000000..79a8570b --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/README.md" @@ -0,0 +1,149 @@ +# CUDA 金融衍生品定价器 + +这是一个使用 C++17 和 CUDA C++ 实现的蒙特卡洛定价命令行项目。它以 CPU 单线程实现作为基线，并提供 GPU 加速路径，用于定价欧式看涨期权和离散算术平均亚式看涨期权。 + +项目不是无参数可执行程序：运行 `pricer_cli` 时必须同时提供期权配置和模拟配置。下面的命令可从克隆仓库开始完成配置、构建、测试与一次实际定价。 + +## 功能范围 + +- 欧式看涨与离散算术平均亚式看涨期权； +- Black-Scholes 欧式看涨解析参考价； +- CPU 单线程蒙特卡洛基线； +- CUDA 蒙特卡洛路径，使用 cuRAND Philox 和 CUB 归约； +- FP32 / FP64 GPU 路径、FP64 统计量、标准误与 95% 置信区间； +- 每次运行输出 JSON 结果，并向性能 CSV 追加记录。 + +## 环境要求 + +- CMake 3.24 或更高版本； +- 支持 C++17 的主机编译器； +- CUDA Toolkit 12.x（含 cuRAND 与 CUB）及可用的 NVIDIA GPU； +- Ninja； +- 配置阶段能访问 GitHub：CMake 会获取固定版本的 nlohmann/json 3.11.3。 + +构建前需要按实际显卡指定 `CMAKE_CUDA_ARCHITECTURES`。例如，RTX 4060 Laptop 的计算能力为 8.9，对应 CMake 值 `89`。Linux/WSL 中可先查看： + +```bash +nvidia-smi --query-gpu=compute_cap --format=csv,noheader +``` + +去掉小数点后填入构建命令；例如 `8.6` 使用 `86`，`8.9` 使用 `89`。 + +## 构建 + +以下命令适用于 Linux 和 WSL。Windows 请在已初始化 MSVC 与 CUDA Toolkit 的开发者终端中执行同样的 CMake 命令；若使用 Ninja，Windows 可执行文件扩展名为 `.exe`。 + +```bash +cmake -S . -B build -G Ninja \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_CUDA_ARCHITECTURES=89 + +cmake --build build -j +``` + +`89` 只是 RTX 4060 Laptop 的示例，不应照搬到其他显卡。正确性构建默认关闭 `PRICER_ENABLE_FAST_MATH`；性能测量也应使用 Release 构建。 + +## 运行定价器 + +先查看本机构建产物支持的参数： + +```bash +./build/pricer_cli --help +``` + +运行欧式看涨期权，并同时执行 CPU 与 GPU 后端： + +```bash +./build/pricer_cli \ + --option configs/option_european_call.ini \ + --simulation configs/simulation.ini \ + --output-dir results/european \ + --backend both \ + --repetitions 1 +``` + +运行亚式看涨期权： + +```bash +./build/pricer_cli \ + --option configs/option_asian_call.ini \ + --simulation configs/simulation.ini \ + --output-dir results/asian \ + --backend both \ + --repetitions 1 +``` + +Windows + Ninja 时将可执行文件写为 `build\\pricer_cli.exe`，路径分隔符可使用 `\\`。若使用 Visual Studio 多配置生成器，则通常从 `build\\Release\\pricer_cli.exe` 启动。 + +`configs/simulation.ini` 默认配置为 1,000 万条路径、256 个时间步，适合作为完整实验输入，运行时间可能较长。需要快速调试时，请复制该文件到新的 INI，再降低 `num_paths` 与 `num_steps`；不要改写仓库提供的示例配置。 + +### 常用 CLI 参数 + +| 参数 | 说明 | +| --- | --- | +| `--option ` | 必填。欧式或亚式期权 INI 文件。 | +| `--simulation ` | 必填。路径数、步数、随机数、精度、block size 等模拟 INI 文件。 | +| `--output-dir ` | 输出目录，默认是 `results`。 | +| `--backend ` | 选择 CPU、GPU 或同时运行两者，默认 `both`。 | +| `--device ` | CUDA 设备编号，默认 `0`。 | +| `--warmup ` | GPU 非计量预热次数，默认 `1`。 | +| `--repetitions ` | 正式计量重复次数，默认 `1`。 | + +## 配置与输出 + +`configs/` 中提供两类配置： + +- `option_european_call.ini` 与 `option_asian_call.ini`：期权类型、现价、行权价、利率、波动率和到期时间； +- `simulation.ini`：FP32 GPU 默认实验配置； +- `simulation_fp64.ini`：FP64 GPU 配置，适合与 CPU FP64 基线进行同精度性能对比。 + +每个后端会在输出目录写入一份： + +```text +---result.json +``` + +性能数据会追加到同目录的 `performance.csv`。JSON 中价格、标准误、置信区间、运行时间和环境信息使用结构化字段记录；不可用统计值以 `null` 表示。 + +CPU 与 GPU 使用不同随机序列，验证时应比较统计一致性，而不是比较逐路径结果。性能结论应使用 CPU 与 GPU 的 `total_runtime_ms`，不要将 CPU 总时间与 GPU `compute_runtime_ms` 混用。 + +## 项目结构 + +```text +include/pricer/ 公共 C++ 接口 +src/ CPU、CUDA、CLI、统计与输出实现 +configs/ 示例期权与模拟 INI +results/ 运行生成的 JSON、CSV 与实验产物（默认不提交） +``` + +## 常见问题 + +### CMake 提示未设置 CUDA 架构 + +本项目要求显式传入 `CMAKE_CUDA_ARCHITECTURES`。确认 GPU 的计算能力后重新配置，例如： + +```bash +cmake -S . -B build -G Ninja \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_CUDA_ARCHITECTURES=89 +``` + +不要复用由其他操作系统、其他 CUDA 版本或其他 GPU 生成的 CMake 缓存。 + +### FetchContent 下载依赖失败 + +检查 GitHub 网络与 Git 是否可用后重新配置。若所在环境已经有固定版本的依赖源码，也可以显式指定它们，避免 CMake 联网下载： + +```bash +cmake -S . -B build -G Ninja \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_CUDA_ARCHITECTURES=89 \ + -DFETCHCONTENT_SOURCE_DIR_NLOHMANN_JSON=/path/to/nlohmann_json \ + -DFETCHCONTENT_SOURCE_DIR_GOOGLETEST=/path/to/googletest +``` + +这些目录必须分别包含对应依赖的顶层 `CMakeLists.txt`。不要将下载缓存、`_deps` 或构建目录提交到仓库。 + +### 无参数运行 `pricer_cli` 报错 + +这是预期行为。`--option` 与 `--simulation` 都是必填参数；请使用“运行定价器”章节中的完整命令。 diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/option_asian_call.ini" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/option_asian_call.ini" new file mode 100644 index 00000000..05510320 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/option_asian_call.ini" @@ -0,0 +1,6 @@ +option_type = asian_call +spot = 100.0 +strike = 100.0 +risk_free_rate = 0.03 +volatility = 0.2 +maturity = 1.0 diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/option_european_call.ini" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/option_european_call.ini" new file mode 100644 index 00000000..fad49bdb --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/option_european_call.ini" @@ -0,0 +1,6 @@ +option_type = european_call +spot = 100.0 +strike = 100.0 +risk_free_rate = 0.03 +volatility = 0.2 +maturity = 1.0 diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/simulation.ini" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/simulation.ini" new file mode 100644 index 00000000..b0844006 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/simulation.ini" @@ -0,0 +1,8 @@ +num_paths = 10000000 +num_steps = 256 +seed = 1234 +rng = curand_philox +variance_reduction = none +precision = fp32 +block_size = 256 +batch_size = 0 diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/simulation_fp64.ini" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/simulation_fp64.ini" new file mode 100644 index 00000000..d525af1a --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/configs/simulation_fp64.ini" @@ -0,0 +1,8 @@ +num_paths = 10000000 +num_steps = 256 +seed = 1234 +rng = curand_philox +variance_reduction = none +precision = fp64 +block_size = 256 +batch_size = 0 diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/black_scholes.hpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/black_scholes.hpp" new file mode 100644 index 00000000..88d798a6 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/black_scholes.hpp" @@ -0,0 +1,12 @@ +#pragma once + +// 提供欧式看涨期权的 Black-Scholes 解析参考值。 + +#include "pricer/config.hpp" + +namespace pricer { + +// 返回欧式看涨的解析价格，用于校验蒙特卡洛估计。 +double black_scholes_call(const OptionParams &option); + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/config.hpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/config.hpp" new file mode 100644 index 00000000..ed2cfc53 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/config.hpp" @@ -0,0 +1,86 @@ +#pragma once + +// 定义期权、模拟和命令行配置的数据结构与读取接口。 + +#include +#include +#include +#include +#include +#include + +namespace pricer { + +// P0 支持的期权类型；亚式期权需要完整路径上的多个监控点。 +enum class OptionType { EuropeanCall, AsianArithmeticCall }; +// GPU 路径演化可选精度；统计累加始终使用双精度。 +enum class Precision { Fp32, Fp64 }; +enum class VarianceReduction { None }; +enum class RngType { CurandPhilox, Mt19937_64 }; + +// 一份期权合约的金融参数。 +struct OptionParams { + OptionType type; + double spot; + double strike; + double risk_free_rate; + double volatility; + double maturity; + std::optional barrier; +}; + +// 一次蒙特卡洛运行的随机数、精度和 GPU 执行参数。 +struct SimulationParams { + std::uint64_t num_paths; + std::uint32_t num_steps; + std::uint64_t seed; + RngType rng; + Precision precision; + VarianceReduction variance_reduction; + std::uint32_t block_size; + std::uint64_t batch_size; +}; + +// 将合约、模拟策略和 CLI 后端选择组合为完整运行配置。 +struct RunConfig { + OptionParams option; + SimulationParams simulation; + std::filesystem::path output_dir; + bool run_cpu; + bool run_gpu; +}; + +enum class ErrorCode { ConfigInvalid, FileRead, UnsupportedFeature }; + +// 保留文件、行号和字段，方便用户定位配置错误。 +class ConfigError : public std::runtime_error { + public: + ConfigError(ErrorCode code, std::filesystem::path source_file, + std::optional line, + std::optional field, std::string reason); + + ErrorCode code() const noexcept; + const std::filesystem::path &source_file() const noexcept; + const std::optional &line() const noexcept; + const std::optional &field() const noexcept; + const std::string &reason() const noexcept; + std::string render() const; + + private: + ErrorCode code_; + std::filesystem::path source_file_; + std::optional line_; + std::optional field_; + std::string reason_; +}; + +// 将配置错误稳定映射为 CLI 退出码。 +int exit_code_for(ErrorCode code) noexcept; + +// 读取两份 INI 文件，并在返回前完成跨字段校验。 +RunConfig load_run_config(const std::filesystem::path &option_file, + const std::filesystem::path &simulation_file, + std::filesystem::path output_dir = "results", + bool run_cpu = true, bool run_gpu = true); + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/cpu_pricer.hpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/cpu_pricer.hpp" new file mode 100644 index 00000000..8f3c7866 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/cpu_pricer.hpp" @@ -0,0 +1,23 @@ +#pragma once + +// 声明 CPU 单线程蒙特卡洛基线的运行结果和定价入口。 + +#include "pricer/config.hpp" +#include "pricer/statistics.hpp" + +namespace pricer { + +// CPU 一次定价产生的累计矩和纯计算耗时。 +struct CpuPricingRun { + RawMoments moments; + double compute_runtime_ms; +}; + +class CpuMonteCarloPricer { + public: + // 按固定 seed 顺序生成全部路径，供正确性对照使用。 + static CpuPricingRun price(const OptionParams &option, + const SimulationParams &simulation); +}; + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/cuda_pricer.hpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/cuda_pricer.hpp" new file mode 100644 index 00000000..8abfd5a8 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/cuda_pricer.hpp" @@ -0,0 +1,66 @@ +#pragma once + +// 声明 CUDA 蒙特卡洛后端及其设备错误转换接口。 + +#include +#include +#include +#include +#include + +#include "pricer/config.hpp" +#include "pricer/statistics.hpp" + +namespace pricer { + +// 将 CUDA API 调用失败转换为带调用位置的 C++ 异常。 +class CudaError : public std::runtime_error { + public: + CudaError(std::string api, int code, std::string cuda_text, + std::string source_file, int source_line); + + const std::string &api() const noexcept; + int code() const noexcept; + const std::string &cuda_text() const noexcept; + const std::string &source_file() const noexcept; + int source_line() const noexcept; + + private: + std::string api_; + int code_; + std::string cuda_text_; + std::string source_file_; + int source_line_; +}; + +// 输出报告所需的目标 GPU 基本信息。 +struct CudaDeviceInfo { + int id; + std::string name; + int compute_capability_major; + int compute_capability_minor; + std::size_t total_memory_bytes; + std::size_t free_memory_bytes; + int max_threads_per_block; +}; + +// GPU 一次定价的累计矩、分批策略和分段计时。 +struct CudaPricingRun { + RawMoments moments; + double gpu_compute_ms; + double reduction_ms; + std::optional rng_ms; + std::uint64_t batch_size; + std::uint64_t grid_size; + CudaDeviceInfo device; +}; + +class CudaMonteCarloPricer { + public: + // 在指定设备上模拟路径并完成收益归约。 + static CudaPricingRun price(const OptionParams &option, + const SimulationParams &simulation, + int device_id = 0); +}; + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/measurement.hpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/measurement.hpp" new file mode 100644 index 00000000..3f9801f6 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/measurement.hpp" @@ -0,0 +1,55 @@ +#pragma once + +// 统一重复运行、预热和中位数计时的测量接口。 + +#include +#include +#include + +#include "pricer/cpu_pricer.hpp" +#include "pricer/cuda_pricer.hpp" + +namespace pricer { + +// 允许测量逻辑复用真实定价器或测试替身。 +using CpuPriceFunction = std::function; +using GpuPriceFunction = std::function; +using AnalyzeFunction = + std::function)>; + +// 对多次 CPU 运行取中位数后的报告数据。 +struct CpuMeasurement { + PricingResult result; + double total_runtime_ms; + double compute_runtime_ms; +}; + +// 对多次 GPU 运行取中位数后的报告数据与设备信息。 +struct GpuMeasurement { + PricingResult result; + double total_runtime_ms; + double compute_runtime_ms; + std::optional rng_ms; + double reduction_ms; + std::uint64_t actual_batch_size; + std::uint64_t grid_size; + CudaDeviceInfo device; +}; + +// 运行 CPU 基线；每次重复都包含统计分析，避免计时口径不一致。 +CpuMeasurement +measure_cpu(const OptionParams &option, const SimulationParams &simulation, + std::uint32_t repetitions, std::optional reference_price, + const CpuPriceFunction &price, const AnalyzeFunction &analyze); + +// 先执行 GPU 预热，再测量正式重复运行。 +GpuMeasurement measure_gpu(const OptionParams &option, + const SimulationParams &simulation, int device, + std::uint32_t warmup, std::uint32_t repetitions, + std::optional reference_price, + const GpuPriceFunction &price, + const AnalyzeFunction &analyze); + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/output.hpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/output.hpp" new file mode 100644 index 00000000..8db4c528 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/output.hpp" @@ -0,0 +1,87 @@ +#pragma once + +// 定义定价结果、性能记录及 JSON/CSV 输出接口。 + +#include +#include +#include +#include +#include +#include + +#include + +#include "pricer/config.hpp" +#include "pricer/statistics.hpp" + +namespace pricer { + +// 一次后端运行的时间、吞吐和加速比数据。 +struct PerformanceData { + double total_runtime_ms; + double compute_runtime_ms; + std::optional rng_ms; + std::optional reduction_ms; + double paths_per_second; + std::uint32_t block_size; + std::uint64_t grid_size; + std::uint32_t repetitions; + std::string aggregation; + std::optional cpu_runtime_ms; + std::optional speedup_vs_cpu; +}; + +// 记录复现实验需要的硬件和构建环境。 +struct EnvironmentData { + std::optional gpu_name; + std::optional compute_capability; + std::optional cuda_runtime_version; + std::optional build_type; + std::optional git_commit; +}; + +// 序列化为 JSON 和 CSV 的完整一次运行记录。 +struct OutputRecord { + std::string run_id; + std::string timestamp_utc; + std::string backend; + OptionParams option; + SimulationParams simulation; + PricingResult result; + PerformanceData performance; + EnvironmentData environment; +}; + +// 成功发布后返回两个输出文件的位置。 +struct WrittenOutput { + std::filesystem::path json_path; + std::filesystem::path csv_path; + std::string run_id; +}; + +class OutputError : public std::runtime_error { + public: + using std::runtime_error::runtime_error; +}; + +// 只构建 JSON 内容，不访问文件系统，便于单元测试。 +nlohmann::json make_result_json(const OutputRecord &record); + +std::filesystem::path +// 临时写入 JSON 后再 rename，避免中途失败留下半个结果文件。 +write_result_json_atomic(const std::filesystem::path &output_dir, + OutputRecord record); + +// 在进程锁保护下创建或追加性能 CSV。 +void append_performance_csv(const std::filesystem::path &path, + const OutputRecord &record); + +// 将 JSON 和 CSV 作为一个输出事务发布。 +WrittenOutput write_output_transaction(const std::filesystem::path &output_dir, + OutputRecord record); + +// 返回计时样本中位数，降低单次抖动的影响。 +double median(std::vector measurements); +std::string utc_timestamp(); + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/payoff.hpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/payoff.hpp" new file mode 100644 index 00000000..44c2fd2e --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/payoff.hpp" @@ -0,0 +1,29 @@ +#pragma once + +// 提供 GBM 步进常量和欧式、亚式看涨收益计算。 + +#include + +namespace pricer { + +// GBM 一步演化和全期折现共用的预计算常量。 +struct GbmStepConstants { + double dt; + double drift; + double diffusion; + double discount; +}; + +// 根据总期限与步数生成 CPU/GPU 共用的风险中性 GBM 常量。 +GbmStepConstants make_gbm_step_constants(double maturity, double risk_free_rate, + double volatility, + std::uint32_t num_steps); + +// 计算一条欧式路径的折现到期收益。 +double european_call_payoff(double terminal_spot, double strike, + double discount); +// 计算一条亚式路径的折现收益；running_sum 不包含初始价格。 +double asian_arithmetic_call_payoff(double running_sum, std::uint32_t num_steps, + double strike, double discount); + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/statistics.hpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/statistics.hpp" new file mode 100644 index 00000000..3ddf4f15 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/statistics.hpp" @@ -0,0 +1,41 @@ +#pragma once + +// 定义原始矩和价格统计量之间的转换接口。 + +#include +#include + +namespace pricer { + +// 不保存所有样本时仍可计算统计量的一阶、二阶原始矩。 +struct RawMoments { + std::uint64_t count; + double sum; + double sum_squares; +}; + +// 面向用户输出的价格、误差和置信区间。 +struct PricingResult { + double price; + std::optional sample_stddev; + std::optional standard_error; + std::optional ci_lower; + std::optional ci_upper; + std::optional reference_price; + std::optional absolute_error; + std::optional relative_error; + RawMoments moments; +}; + +// 合并批次或分段运行产生的原始矩。 +RawMoments merge_raw_moments(const RawMoments &left, const RawMoments &right); + +class ResultAnalyzer { + public: + // 从原始矩推导样本标准差、标准误和 95% 置信区间。 + static PricingResult + analyze(const RawMoments &moments, + std::optional reference_price = std::nullopt); +}; + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/version.hpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/version.hpp" new file mode 100644 index 00000000..28a5ea0a --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/include/pricer/version.hpp" @@ -0,0 +1,10 @@ +#pragma once + +// 暴露命令行显示用的项目版本字符串。 + +namespace pricer { + +// 返回稳定的版本文本，不分配内存也不会抛出异常。 +const char *version() noexcept; + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-01-system-architecture.png" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-01-system-architecture.png" new file mode 100644 index 00000000..2208486a Binary files /dev/null and "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-01-system-architecture.png" differ diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-02-convergence.svg" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-02-convergence.svg" new file mode 100644 index 00000000..bb734a19 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-02-convergence.svg" @@ -0,0 +1,71 @@ + + +图 2 欧式与亚式蒙特卡洛收敛（FP32 GPU） +标准误随路径数的幂律下降，并与理论 -1/2 斜率对照。 +{"generator": "scripts/render_final_report_figures.py", "result_root": "final-report-20260915", "sources": [{"path": "convergence/european-fp32/convergence_summary.csv", "sha256": "5c324b4265482661e80e23bcffbd8af570e67732da07b03b6e08f082645c4313"}, {"path": "convergence/asian-fp32/convergence_summary.csv", "sha256": "ae3da223b6ff9db31c5c0e4f60f6d5ad193b65a1eaacbd63a59fe59c4ee82549"}, {"path": "benchmark/european-fp64/benchmark_summary.csv", "sha256": "3722ad2fcb3bf3b9b6d9ad17718602336445778bb9ebbe59a6a662c2e08d0d9b"}, {"path": "benchmark/asian-fp64/benchmark_summary.csv", "sha256": "0154c83bfb0b26dcc7374541cd8db2c1df452f7069ba5d7b1c7d9ba16c945440"}, {"path": "block-sweep/asian-fp32/block_sweep.csv", "sha256": "e463e14b70bfb5a8334fa609e450a132582c83d00d8a7efc36aef60057700dd6"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_kern_sum_cuda_gpu_kern_sum.csv", "sha256": "f1af93a6c26adde4c62a99ffbfc2d70460de55c2585e450088b53de60473520c"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_mem_time_sum_cuda_gpu_mem_time_sum.csv", "sha256": "2f1fc8fa2953d96432e314f4f274e9fabcca31ce2665fbfbba0ca41332a57437"}, {"path": "profiling/asian-fp32/nsys-cuda_api_sum_cuda_api_sum.csv", "sha256": "728c8aeffa4972c0950a5ba94d6b5572576ab85db524266faaa4eaf0cff1d6dc"}, {"path": "profiling/asian-fp32/ncu-details.csv", "sha256": "10c85cf76b074afd0a58f41f80f19914cd849e5af5f93a9ee9346abb191044a4"}]} + + +欧式看涨 + + + +1e+04 + +1e+05 + +1e+06 + +1e+07 + +2.2e-03 + +7.5e-03 + +2.5e-02 + +8.5e-02 + +2.9e-01 +路径数 +标准误 + + + + + + +拟合斜率 -0.501；理论对照 -0.500 +Black–Scholes 参考价 9.4134 +亚式看涨 + + + +1e+04 + +1e+05 + +1e+06 + +1e+07 + +1.2e-03 + +4.1e-03 + +1.4e-02 + +4.6e-02 + +1.5e-01 +路径数 +标准误 + + + + + + +拟合斜率 -0.498；理论对照 -0.500 +亚式无可得解析参考价，不显示参考线 +虚线为 -1/2 幂律对照线；数据来源见 SVG metadata。 + diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-03-correctness-ci.svg" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-03-correctness-ci.svg" new file mode 100644 index 00000000..a89c7497 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-03-correctness-ci.svg" @@ -0,0 +1,53 @@ + + +图 3 正确性与 95% 置信区间（FP64，1000 万路径） +CPU/GPU 合并标准误的误差棒，欧式对照 Black–Scholes 参考价。 +{"generator": "scripts/render_final_report_figures.py", "result_root": "final-report-20260915", "sources": [{"path": "convergence/european-fp32/convergence_summary.csv", "sha256": "5c324b4265482661e80e23bcffbd8af570e67732da07b03b6e08f082645c4313"}, {"path": "convergence/asian-fp32/convergence_summary.csv", "sha256": "ae3da223b6ff9db31c5c0e4f60f6d5ad193b65a1eaacbd63a59fe59c4ee82549"}, {"path": "benchmark/european-fp64/benchmark_summary.csv", "sha256": "3722ad2fcb3bf3b9b6d9ad17718602336445778bb9ebbe59a6a662c2e08d0d9b"}, {"path": "benchmark/asian-fp64/benchmark_summary.csv", "sha256": "0154c83bfb0b26dcc7374541cd8db2c1df452f7069ba5d7b1c7d9ba16c945440"}, {"path": "block-sweep/asian-fp32/block_sweep.csv", "sha256": "e463e14b70bfb5a8334fa609e450a132582c83d00d8a7efc36aef60057700dd6"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_kern_sum_cuda_gpu_kern_sum.csv", "sha256": "f1af93a6c26adde4c62a99ffbfc2d70460de55c2585e450088b53de60473520c"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_mem_time_sum_cuda_gpu_mem_time_sum.csv", "sha256": "2f1fc8fa2953d96432e314f4f274e9fabcca31ce2665fbfbba0ca41332a57437"}, {"path": "profiling/asian-fp32/nsys-cuda_api_sum_cuda_api_sum.csv", "sha256": "728c8aeffa4972c0950a5ba94d6b5572576ab85db524266faaa4eaf0cff1d6dc"}, {"path": "profiling/asian-fp32/ncu-details.csv", "sha256": "10c85cf76b074afd0a58f41f80f19914cd849e5af5f93a9ee9346abb191044a4"}]} + + + + + +4.550 + +5.955 + +7.361 + +8.766 + +10.171 +期权类型与后端 +价格估计 +欧式看涨 + + + + +9.4183 + + + + +9.4186 + +Black–Scholes 9.4134 +亚式看涨 + + + + +5.3007 + + + + +5.2987 + +CPU + +GPU + +Black–Scholes 参考价 +误差棒为各后端 95% 置信区间；亚式不虚构解析参考线。 + diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-04-performance-scaling.svg" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-04-performance-scaling.svg" new file mode 100644 index 00000000..d1cfaaff --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-04-performance-scaling.svg" @@ -0,0 +1,95 @@ + + +图 4 CPU/GPU 性能尺度（FP64，中位数） +左面板比较 total 墙钟时序，右面板给出由 total 之比计算的加速比。 +{"generator": "scripts/render_final_report_figures.py", "result_root": "final-report-20260915", "sources": [{"path": "convergence/european-fp32/convergence_summary.csv", "sha256": "5c324b4265482661e80e23bcffbd8af570e67732da07b03b6e08f082645c4313"}, {"path": "convergence/asian-fp32/convergence_summary.csv", "sha256": "ae3da223b6ff9db31c5c0e4f60f6d5ad193b65a1eaacbd63a59fe59c4ee82549"}, {"path": "benchmark/european-fp64/benchmark_summary.csv", "sha256": "3722ad2fcb3bf3b9b6d9ad17718602336445778bb9ebbe59a6a662c2e08d0d9b"}, {"path": "benchmark/asian-fp64/benchmark_summary.csv", "sha256": "0154c83bfb0b26dcc7374541cd8db2c1df452f7069ba5d7b1c7d9ba16c945440"}, {"path": "block-sweep/asian-fp32/block_sweep.csv", "sha256": "e463e14b70bfb5a8334fa609e450a132582c83d00d8a7efc36aef60057700dd6"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_kern_sum_cuda_gpu_kern_sum.csv", "sha256": "f1af93a6c26adde4c62a99ffbfc2d70460de55c2585e450088b53de60473520c"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_mem_time_sum_cuda_gpu_mem_time_sum.csv", "sha256": "2f1fc8fa2953d96432e314f4f274e9fabcca31ce2665fbfbba0ca41332a57437"}, {"path": "profiling/asian-fp32/nsys-cuda_api_sum_cuda_api_sum.csv", "sha256": "728c8aeffa4972c0950a5ba94d6b5572576ab85db524266faaa4eaf0cff1d6dc"}, {"path": "profiling/asian-fp32/ncu-details.csv", "sha256": "10c85cf76b074afd0a58f41f80f19914cd849e5af5f93a9ee9346abb191044a4"}]} + + +欧式看涨 + + + +1e+04 + +1e+05 + +1e+06 + +1e+07 + +0 + +1 + +8 + +60 + +469 +路径数 +total 时序 (ms) + + + + + + + + + + + +0.1× + +1.8× + +7.8× + +13.6× +亚式看涨 + + + +1e+04 + +1e+05 + +1e+06 + +1e+07 + +2 + +32 + +454 + +6501 + +93040 +路径数 +total 时序 (ms) + + + + + + + + + + + +13.8× + +21.4× + +33.7× + +38.9× + +CPU total + +GPU total +加速比由同一档 CPU/GPU total 墙钟之比重新计算；吞吐才使用 GPU compute 时序。 + diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-05-block-size.svg" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-05-block-size.svg" new file mode 100644 index 00000000..abcc788d --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-05-block-size.svg" @@ -0,0 +1,39 @@ + + +图 5 block-size 决策（FP32 亚式，1000 万路径） +三个候选的相对变化，仅当双指标均改善至少 5% 才替换基准 256。 +{"generator": "scripts/render_final_report_figures.py", "result_root": "final-report-20260915", "sources": [{"path": "convergence/european-fp32/convergence_summary.csv", "sha256": "5c324b4265482661e80e23bcffbd8af570e67732da07b03b6e08f082645c4313"}, {"path": "convergence/asian-fp32/convergence_summary.csv", "sha256": "ae3da223b6ff9db31c5c0e4f60f6d5ad193b65a1eaacbd63a59fe59c4ee82549"}, {"path": "benchmark/european-fp64/benchmark_summary.csv", "sha256": "3722ad2fcb3bf3b9b6d9ad17718602336445778bb9ebbe59a6a662c2e08d0d9b"}, {"path": "benchmark/asian-fp64/benchmark_summary.csv", "sha256": "0154c83bfb0b26dcc7374541cd8db2c1df452f7069ba5d7b1c7d9ba16c945440"}, {"path": "block-sweep/asian-fp32/block_sweep.csv", "sha256": "e463e14b70bfb5a8334fa609e450a132582c83d00d8a7efc36aef60057700dd6"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_kern_sum_cuda_gpu_kern_sum.csv", "sha256": "f1af93a6c26adde4c62a99ffbfc2d70460de55c2585e450088b53de60473520c"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_mem_time_sum_cuda_gpu_mem_time_sum.csv", "sha256": "2f1fc8fa2953d96432e314f4f274e9fabcca31ce2665fbfbba0ca41332a57437"}, {"path": "profiling/asian-fp32/nsys-cuda_api_sum_cuda_api_sum.csv", "sha256": "728c8aeffa4972c0950a5ba94d6b5572576ab85db524266faaa4eaf0cff1d6dc"}, {"path": "profiling/asian-fp32/ncu-details.csv", "sha256": "10c85cf76b074afd0a58f41f80f19914cd849e5af5f93a9ee9346abb191044a4"}]} + + + + +block size +相对 256 的变化 (%) + + +改善阈值 5% + +改善阈值 -5% +128 + ++0.3% + +-1.0% +256（保留） + ++0.0% + ++0.0% +512 + ++3.2% + ++2.3% + +保留候选 + +其他候选 + +±5% 判定阈值 +每档左柱为 total 中位数变化，右柱为 compute 中位数变化；仅双指标同时越过 -5% 才替换。 + diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-06-nsight-bottleneck.svg" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-06-nsight-bottleneck.svg" new file mode 100644 index 00000000..c65032b9 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/figures/figure-06-nsight-bottleneck.svg" @@ -0,0 +1,44 @@ + + +图 6 Nsight 瓶颈剖析（FP32 亚式代表负载） +左面板为内核时间占比，右面板为 ncu 的 SM/DRAM 与占用率指标。 +{"generator": "scripts/render_final_report_figures.py", "result_root": "final-report-20260915", "sources": [{"path": "convergence/european-fp32/convergence_summary.csv", "sha256": "5c324b4265482661e80e23bcffbd8af570e67732da07b03b6e08f082645c4313"}, {"path": "convergence/asian-fp32/convergence_summary.csv", "sha256": "ae3da223b6ff9db31c5c0e4f60f6d5ad193b65a1eaacbd63a59fe59c4ee82549"}, {"path": "benchmark/european-fp64/benchmark_summary.csv", "sha256": "3722ad2fcb3bf3b9b6d9ad17718602336445778bb9ebbe59a6a662c2e08d0d9b"}, {"path": "benchmark/asian-fp64/benchmark_summary.csv", "sha256": "0154c83bfb0b26dcc7374541cd8db2c1df452f7069ba5d7b1c7d9ba16c945440"}, {"path": "block-sweep/asian-fp32/block_sweep.csv", "sha256": "e463e14b70bfb5a8334fa609e450a132582c83d00d8a7efc36aef60057700dd6"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_kern_sum_cuda_gpu_kern_sum.csv", "sha256": "f1af93a6c26adde4c62a99ffbfc2d70460de55c2585e450088b53de60473520c"}, {"path": "profiling/asian-fp32/nsys-cuda_gpu_mem_time_sum_cuda_gpu_mem_time_sum.csv", "sha256": "2f1fc8fa2953d96432e314f4f274e9fabcca31ce2665fbfbba0ca41332a57437"}, {"path": "profiling/asian-fp32/nsys-cuda_api_sum_cuda_api_sum.csv", "sha256": "728c8aeffa4972c0950a5ba94d6b5572576ab85db524266faaa4eaf0cff1d6dc"}, {"path": "profiling/asian-fp32/ncu-details.csv", "sha256": "10c85cf76b074afd0a58f41f80f19914cd849e5af5f93a9ee9346abb191044a4"}]} + + + + +内核时间占比 (%) + + + + +asian_payoff_kernel 97.96% + +CUB 归约 2.04% + +其他内核 0.00% +内核时间构成 +内存传输总时长相当于内核时间的 0.0% +占比相对内核总时间，不代表完整墙钟时长 +ncu 关键指标 (%) +SM 吞吐 + + +85.0 +DRAM 吞吐 + + +2.2 +实际占用率 + + +98.9 +理论占用率 + + +100.0 +每线程寄存器 34 +本地内存溢出请求 0 次 +grid 39062 × block 256 +两类占比分别相对内核总时间与 ncu 指标口径，不可写作完整墙钟占比。 + diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/\351\207\221\350\236\215\350\241\215\347\224\237\345\223\201\345\256\232\344\273\267\344\270\216\351\243\216\351\231\251\344\274\260\350\256\241 CUDA \346\212\245\345\221\212.md" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/\351\207\221\350\236\215\350\241\215\347\224\237\345\223\201\345\256\232\344\273\267\344\270\216\351\243\216\351\231\251\344\274\260\350\256\241 CUDA \346\212\245\345\221\212.md" new file mode 100644 index 00000000..ae65fc2b --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/report/\351\207\221\350\236\215\350\241\215\347\224\237\345\223\201\345\256\232\344\273\267\344\270\216\351\243\216\351\231\251\344\274\260\350\256\241 CUDA \346\212\245\345\221\212.md" @@ -0,0 +1,186 @@ +# 金融衍生品定价与风险估计 CUDA 报告 + +报告日期：2026-09-15
+冻结结果根：`results/final-report-20260915/`
+运行环境：RTX 4060 Laptop GPU / CUDA 12.9 / SM89 / Release / fast-math OFF / seed 1234 + +## 摘要与关键结论 + +本项目实现了一个基于 C++17/CUDA 的蒙特卡洛期权定价工具，覆盖欧式看涨和离散算术平均亚式看涨两类产品。程序提供 Black-Scholes 解析值、单线程 CPU 基线和 CUDA GPU 实现，并输出价格、标准误、95% 置信区间与性能日志。 + +实现中有三处直接影响实验结果：Philox 将路径编号固定映射到随机流；程序按可用显存分批计算，每批只回传 `sum` 和 `sum_squares`；收益及其平方以 FP64 完成 CUB 归约。这样既能控制显存占用，也便于在大规模路径数下稳定统计和复现实验。 + +| 关键指标 | 冻结结果 | +| --- | --- | +| 正确性门槛 | 86/86 CTest 通过；两个 CUDA 测试程序的 Compute Sanitizer 均为 `ERROR SUMMARY: 0 errors` | +| 收敛行为 | 欧式标准误斜率 −0.5014，亚式 −0.4980，均符合理论 −0.5 | +| 最大端到端加速 | FP64 亚式、1000 万路径为 **38.90×** | +| 关键性能边界 | FP64 欧式在 1 万路径为 **0.135×**，GPU 因固定开销反而更慢 | + +文中的数值均来自本轮冻结产物。CPU/GPU 加速比统一按对应的 `total_runtime_ms` 计算。离散算术平均亚式期权没有闭式解析解，因此不设置虚构的 `reference_price` 或 `absolute_error`。 + +## 目录 + +- [1. 问题、方案与贡献](#1-问题方案与贡献) +- [2. 金融模型与数值方法](#2-金融模型与数值方法) + - [2.1 欧式与亚式收益](#21-欧式与亚式收益) + - [2.2 随机数、统计与时间口径](#22-随机数统计与时间口径) +- [3. 系统架构与 CUDA 设计](#3-系统架构与-cuda-设计) + - [3.1 路径与内存策略](#31-路径与内存策略) + - [3.2 归约与可观测性](#32-归约与可观测性) +- [4. 实验环境与复现](#4-实验环境与复现) +- [5. 正确性与收敛](#5-正确性与收敛) +- [6. CPU/GPU 性能](#6-cpugpu-性能) +- [7. P1 剖析与 block-size 决策](#7-p1-剖析与-block-size-决策) +- [8. 精度取舍](#8-精度取舍) +- [9. 整体结论](#9-整体结论) +- [10. 局限与后续方向](#10-局限与后续方向) +- [11. 附录：证据索引](#11-附录证据索引) + +## 1. 问题、方案与贡献 + +蒙特卡洛定价需要生成大量相互独立的价格路径，路径之间几乎没有数据依赖，因此适合放到 GPU 上并行执行。本项目关注的不只是给出一个价格，还要检查 CPU/GPU 的统计结果能否相互印证、结果是否带有不确定度描述，以及性能结论是否有端到端计时和剖析数据支撑。 + +CPU 部分采用单线程实现，作为便于核对的正确性和性能基线。它并不代表最优 CPU 性能，因此本文的加速比应理解为单线程 CPU 与 GPU 的对比。对这类任务，路径规模和单条路径的计算量比“是否使用 GPU”更关键：亚式期权每条路径的计算更重，更容易摊薄 GPU 的固定开销。 + +## 2. 金融模型与数值方法 + +### 2.1 欧式与亚式收益 + +在风险中性测度下，欧式看涨的终端价格按一步 GBM 生成： + +$$S_T = S_0 \exp\left[(r-q-\tfrac{1}{2}\sigma^2)T+\sigma\sqrt{T}Z\right],\quad Z\sim N(0,1)$$ + +欧式看涨的折现收益为 $e^{-rT}\max(S_T-K,0)$。本轮实验的 Black-Scholes 参考价为 **9.413403383853016**。 + +离散算术平均亚式看涨以 $m$ 个监控点价格的平均值计算收益： + +$$\mathrm{Payoff}=\max\left(\frac{1}{m}\sum_{i=1}^{m}S_{t_i}-K,0\right)$$ + +初始价格 $S_0$ 不计入平均值。离散算术平均亚式没有闭式解析解，因而主要通过 CPU/GPU 的统计一致性、置信区间和收敛情况判断实现是否合理。 + +### 2.2 随机数、统计与时间口径 + +GPU 使用 `curandStatePhilox4_32_10_t`。`curand_init(seed, batch_offset + local_path, 0)` 将全局路径编号作为 subsequence，因此即使改变分批策略，同一条路径仍使用同一随机流。CPU 使用 `mt19937_64`，两端随机流不同，不应要求逐路径或逐位结果一致。 + +标准误由路径收益的样本方差计算，95% 置信区间为 $\hat{V}\pm1.96\cdot\mathrm{SE}$。比较 CPU 与 GPU 时，使用 $\sqrt{SE_{cpu}^2+SE_{gpu}^2}$ 作为合并标准误，衡量两项独立估计之差的正常波动范围。 + +`total_runtime_ms` 包含初始化、分配、计算和回传；`compute_runtime_ms` 只统计核函数区间。本文的 CPU/GPU 加速比均采用 `total_runtime_ms`，反映实际端到端等待时间。 + +## 3. 系统架构与 CUDA 设计 + +![图 1：命令行入口、配置校验、三条定价路径与结果输出的系统架构](figures/figure-01-system-architecture.png) + +图 1 展示了当前实现的模块关系。命令行读取并校验 INI 配置后，程序选择 Black-Scholes、CPU 蒙特卡洛或 CUDA 蒙特卡洛路径；`ResultAnalyzer` 负责统一统计结果，并以事务方式写入 JSON 和性能 CSV。 + +### 3.1 路径与内存策略 + +CUDA 实现采用一线程一条路径和 grid-stride loop。欧式期权只需生成终端价格；亚式期权则在单个线程内推进 256 个时间步，并同步维护路径平均值。程序不保存完整的路径矩阵。 + +`choose_batch()` 根据可用显存确定批大小，并预留约 20% 空间。每批只从设备端回传两个 FP64 标量，主机通过 `merge_raw_moments()` 合并各批矩。因此，路径数增加不会带来完整路径结果的主机传输负担。 + +### 3.2 归约与可观测性 + +payoff 及 payoff 平方分别通过 `cub::DeviceReduce::Sum` 归约，输入、输出和最终统计均为 FP64。NVTX 可以标记 `simulate_paths`、`reduce_moments` 和 `copy_batch_moments`，便于把端到端时间拆分为具体的 GPU 活动。 + +## 4. 实验环境与复现 + +| 项目 | 值 | +| --- | --- | +| GPU | NVIDIA GeForce RTX 4060 Laptop GPU，8188 MiB，compute capability 8.9 | +| 软件 | CUDA 12.9 / 驱动 595.71 / CMake 3.28.3 / Ninja 1.11.1 | +| 构建 | Release，SM89，fast-math OFF，正式构建 NVTX OFF | +| 实验 | seed 1234；GPU 预热 1 次；正式运行 5 次取中位数 | +| 路径数 | 10,000 / 100,000 / 1,000,000 / 10,000,000 | + +正式结果固定写入 `results/final-report-20260915/`。附录给出了构建、Sanitizer、扫描、剖析和绘图命令；该目录保存环境清单、输入哈希、原始 CSV/JSON、日志和 Nsight 导出物。 + +## 5. 正确性与收敛 + +![图 2：FP32 GPU 下欧式与亚式标准误随路径数的幂律收敛](figures/figure-02-convergence.svg) + +图 2 中，欧式和亚式的标准误经验斜率分别为 **−0.5014** 和 **−0.4980**，与蒙特卡洛理论值 −0.5 接近。欧式绝对误差不必随路径数单调下降；判断误差量级时，应结合标准误，而不是只看估计值偏离参考价的方向。 + +![图 3：1000 万路径下 CPU 与 GPU 的价格估计及 95% 置信区间](figures/figure-03-correctness-ci.svg) + +在 1000 万路径下，欧式期权的 CPU、GPU 估计都覆盖 Black-Scholes 参考值。亚式期权的 CPU/GPU 差异为 2.04e−3，合并标准误为 3.465e−3，前者约为后者的 0.59 倍。由于两端使用不同随机流，这种统计比较比逐路径相等更合适。 + +## 6. CPU/GPU 性能 + +![图 4：两类期权的 CPU/GPU 端到端耗时与加速比](figures/figure-04-performance-scaling.svg) + +图 4 使用端到端时间比较性能。对于 FP64 欧式期权，1 万路径时 CPU 总耗时为 0.214687 ms，GPU 为 1.590332 ms，加速比只有 **0.135×**；此时初始化、分配和同步成本尚未被计算量摊薄。路径数增至 1000 万后，加速比达到 **13.595×**。 + +亚式期权每条路径需要推进 256 个时间步，计算密度更高：1 万路径时已达到 **13.759×**，1000 万路径时为 **38.898×**。这说明 GPU 的优势取决于任务规模和计算密度，小规模任务并不一定适合迁移到 GPU。 + +| 1000 万路径，FP64 | CPU total (ms) | GPU total (ms) | 加速比 | +| --- | ---: | ---: | ---: | +| 欧式看涨 | 260.402149 | 19.154564 | 13.595× | +| 亚式看涨 | 51,688.686487 | 1,328.820065 | 38.898× | + +## 7. P1 剖析与 block-size 决策 + +![图 5：block-size 扫描下 total 与 compute 中位数对照](figures/figure-05-block-size.svg) + +以 256 为基准时，只有候选 block size 的 `total_runtime_ms` 和 `compute_runtime_ms` 都至少改善 5%，才会替换当前配置。128 的 total 反而增加 0.35%；512 的 total 和 compute 分别增加 3.24% 和 2.27%。因此本轮保留 **256**。 + +![图 6：Nsight 剖析下的内核时间构成与瓶颈判定](figures/figure-06-nsight-bottleneck.svg) + +剖析结果显示，`asian_payoff_kernel` 占内核时间的 97.96%，SM 吞吐为 84.97%，DRAM 吞吐只有 2.24%，实际占用率为 98.88%，本地内存溢出请求为 0。该负载主要受计算吞吐限制，而不是 DRAM 带宽限制；在这一前提下，小幅调整 block size 的收益有限。批处理只回传标量，D2H 时间约为内核总时间的 0.02%。 + +## 8. 精度取舍 + +当前冻结指标中，FP32 与 FP64 的可用计时字段口径不一致，不能据此比较两种精度的性能，也不报告性能倍数。两种精度的价格估计处于采样误差的相近量级；但本轮 FP32/FP64 同时改变了随机流，不能把差异全部归因于数值精度。若要单独考察精度影响，需要使用共同随机数，并在同一计时边界下记录 FP32/FP64 的 `compute_runtime_ms`。 + +## 9. 整体结论 + +**数值结果通过了本轮验证。** 86/86 CTest 通过，两个 CUDA 测试程序的 Compute Sanitizer 均报告 `ERROR SUMMARY: 0 errors`。欧式和亚式的标准误斜率分别为 **−0.5014**、**−0.4980**；欧式结果覆盖解析值，亚式结果则通过合并标准误比较进行核对。 + +**GPU 的收益取决于负载。** FP64 亚式在 1000 万路径下的端到端加速为 **38.90×**，欧式为 13.59×；但欧式在 1 万路径下只有 **0.135×**。因此，对计算量较小的任务，GPU 未必比单线程 CPU 更合适。 + +**block_size = 256 是本轮扫描的保守选择。** 亚式核函数的 SM 吞吐为 84.97%，实际占用率为 98.88%，且未出现本地内存溢出。128、256、512 三种 block size 中，没有候选值同时让 total 与 compute 中位数改善至少 5%，因此继续采用 256。 + +## 10. 局限与后续方向 + +1. **性能适用范围。** 本轮结果只对应 RTX 4060 Laptop GPU、驱动 595.71、WSL2 和当时的热状态，未进行长时间温控测试，也不应直接外推到其他硬件。 +2. **CPU 基线公平性。** CPU 为单线程串行基线。加速比体现的是它与大规模并行 GPU 的差异，不等同于最优 CPU 与最优 GPU 的比较。 +3. **随机流差异。** CPU 使用 `mt19937_64`，GPU 使用 Philox；跨后端只能比较统计一致性，不能要求逐路径或逐位相同。 +4. **精度解耦不足。** FP32/FP64 对照同时改变了随机流。要隔离纯精度影响，需要采用共同随机数或固定随机流。 +5. **方差缩减未实现。** 尚未使用对偶变量、控制变量或重要性采样，收敛仍处于标准的 −1/2 级别。 +6. **剖析为代表性采样。** nsys/ncu 指标来自一次 FP32 亚式代表负载采集，并非多次统计结果。 +7. **亚式没有解析基准。** 离散算术平均亚式没有闭式解，只能借助跨后端一致性和收敛证据建立信心。 + +后续可先加入对偶变量等方差缩减方法，再补充多线程 CPU 基线和共同随机数下的精度实验，最后评估 block-size 与 batch 策略的自动调优。 + +## 11. 附录：证据索引 + +### 11.1 复现入口 + +```bash +python3 scripts/run_final_report.py \ + --source-root . --build-dir build-final-report-20260915 \ + --output-root results/final-report-20260915 --phase validate +python3 scripts/run_final_report.py \ + --source-root . --build-dir build-final-report-20260915 \ + --output-root results/final-report-20260915 --phase block-sweep +python3 scripts/run_final_report.py \ + --source-root . --build-dir build-final-report-20260915 \ + --output-root results/final-report-20260915 --phase formal +python3 scripts/render_final_report_figures.py \ + --input-root results/final-report-20260915 --output-dir report/figures +``` + +NVTX 构建与 nsys/ncu 采集命令记录在 `results/final-report-20260915/profiling/asian-fp32/profile_commands.txt`。输出根名称受脚本保护，避免误写入已废弃的 `final-report-20260909`。 + +### 11.2 图表与原始产物 + +| 图 | 数据源 | +| --- | --- | +| 图 1 系统架构 | `report/figures/figure-01-system-architecture.png`；源 `Docs/architecture/derivative_pricer.architecture.json` | +| 图 2 收敛 | `convergence/{european,asian}-fp32/convergence_summary.csv` | +| 图 3 正确性与 CI | `benchmark/{european,asian}-fp64/benchmark_summary.csv` | +| 图 4 性能扩展 | `benchmark/{european,asian}-fp64/benchmark_summary.csv` | +| 图 5 block-size | `block-sweep/asian-fp32/block_sweep.csv` | +| 图 6 Nsight 瓶颈 | `profiling/asian-fp32/{nsys-*.csv,ncu-details.csv}` | + +`manifest.json` 记录环境与构建选项，`input-sha256.txt` 记录 33 个输入文件与可执行文件哈希，`validation/` 保存构建、CTest 和两份 memcheck 日志，`report/figures/report_metrics.json` 保存每个派生指标的源文件、行键与公式。五张数据 SVG 的 metadata 内嵌生成脚本、结果根、输入相对路径和 SHA-256，用于避免历史结果混入。 diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/black_scholes.cpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/black_scholes.cpp" new file mode 100644 index 00000000..06733e84 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/black_scholes.cpp" @@ -0,0 +1,74 @@ +// 计算欧式看涨的解析价格，并处理零波动率边界。 +#include "pricer/black_scholes.hpp" + +#include +#include +#include + +namespace pricer { +namespace { + +// 在代入解析公式前排除类型、范围和非有限数错误。 +void validate_option(const OptionParams &option) { + if (option.type != OptionType::EuropeanCall) { + throw std::domain_error("Black-Scholes requires a European call"); + } + if (!std::isfinite(option.spot) || !std::isfinite(option.strike) || + !std::isfinite(option.risk_free_rate) || + !std::isfinite(option.volatility) || !std::isfinite(option.maturity)) { + throw std::domain_error("Black-Scholes parameters must be finite"); + } + if (option.spot <= 0.0 || option.strike <= 0.0 || option.maturity <= 0.0 || + option.volatility < 0.0) { + throw std::domain_error("Black-Scholes parameters are out of range"); + } +} + +// 用 erfc 表示标准正态 CDF，数值上比手写积分更稳定。 +double normal_cdf(double value) { + return 0.5 * std::erfc(-value / std::sqrt(2.0)); +} + +void require_finite(double value) { + if (!std::isfinite(value)) { + throw std::overflow_error("Black-Scholes result is not finite"); + } +} + +} // namespace + +// 解析公式只适用于欧式看涨；亚式期权没有同样的闭式解。 +double black_scholes_call(const OptionParams &option) { + validate_option(option); + + const double discount = std::exp(-option.risk_free_rate * option.maturity); + require_finite(discount); + + // sigma=0 时没有随机性，常规 d1/d2 公式会除以零；直接计算唯一的到期价格。 + if (option.volatility == 0.0) { + const double terminal_spot = + option.spot * std::exp(option.risk_free_rate * option.maturity); + const double price = + discount * std::max(terminal_spot - option.strike, 0.0); + require_finite(price); + return price; + } + + // Black-Scholes 是欧式看涨的解析“标准答案”。它不参与亚式定价， + // 而是用来检查蒙特卡洛平均值是否落在合理误差范围内。 + const double volatility_sqrt_time = + option.volatility * std::sqrt(option.maturity); + require_finite(volatility_sqrt_time); + const double d1 = + (std::log(option.spot / option.strike) + + (option.risk_free_rate + 0.5 * option.volatility * option.volatility) * + option.maturity) / + volatility_sqrt_time; + const double d2 = d1 - volatility_sqrt_time; + const double price = option.spot * normal_cdf(d1) - + option.strike * discount * normal_cdf(d2); + require_finite(price); + return price; +} + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/config.cpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/config.cpp" new file mode 100644 index 00000000..0a6b2c0c --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/config.cpp" @@ -0,0 +1,373 @@ +// 解析 INI 配置和 CLI 参数，尽早报告可定位的输入错误。 +#include "pricer/config.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace pricer { +namespace { + +struct Entry { + std::string value; + std::size_t line; +}; + +using Entries = std::unordered_map; + +std::string trim(std::string_view text) { + const auto first = text.find_first_not_of(" \t\r\n"); + if (first == std::string_view::npos) { + return {}; + } + const auto last = text.find_last_not_of(" \t\r\n"); + return std::string(text.substr(first, last - first + 1)); +} + +const char *code_name(ErrorCode code) noexcept { + switch (code) { + case ErrorCode::ConfigInvalid: + return "CONFIG_INVALID"; + case ErrorCode::FileRead: + return "FILE_READ"; + case ErrorCode::UnsupportedFeature: + return "UNSUPPORTED_FEATURE"; + } + return "UNKNOWN"; +} + +[[noreturn]] void fail(ErrorCode code, const std::filesystem::path &file, + std::optional line, + std::optional field, + const std::string &reason) { + throw ConfigError(code, file, std::move(line), std::move(field), reason); +} + +std::string strip_comment(const std::string &line, + const std::filesystem::path &file, + std::size_t line_number) { + bool quoted = false; + for (std::size_t index = 0; index < line.size(); ++index) { + if (line[index] == '"') { + quoted = !quoted; + } else if (line[index] == '#' && !quoted) { + return line.substr(0, index); + } + } + if (quoted) { + fail(ErrorCode::ConfigInvalid, file, line_number, std::nullopt, + "unmatched double quote"); + } + return line; +} + +Entries parse_ini(const std::filesystem::path &path, + const std::unordered_set &allowed_keys) { + std::ifstream stream(path); + if (!stream) { + fail(ErrorCode::FileRead, path, std::nullopt, std::nullopt, + "could not open file"); + } + + Entries entries; + std::string raw_line; + std::size_t line_number = 0; + while (std::getline(stream, raw_line)) { + ++line_number; + const auto line = trim(strip_comment(raw_line, path, line_number)); + if (line.empty()) { + continue; + } + + const auto equals = line.find('='); + if (equals == std::string::npos) { + fail(ErrorCode::ConfigInvalid, path, line_number, std::nullopt, + "expected key = value"); + } + const auto key = trim(std::string_view(line).substr(0, equals)); + auto value = trim(std::string_view(line).substr(equals + 1)); + if (key.empty()) { + fail(ErrorCode::ConfigInvalid, path, line_number, std::nullopt, + "empty key"); + } + if (allowed_keys.find(key) == allowed_keys.end()) { + fail(ErrorCode::ConfigInvalid, path, line_number, key, + "unknown key"); + } + if (entries.find(key) != entries.end()) { + fail(ErrorCode::ConfigInvalid, path, line_number, key, + "duplicate key"); + } + if (!value.empty() && value.front() == '"' && value.back() == '"' && + value.size() >= 2) { + value = value.substr(1, value.size() - 2); + } else if (value.find('"') != std::string::npos) { + fail(ErrorCode::ConfigInvalid, path, line_number, key, + "unmatched double quote"); + } + entries.emplace(key, Entry{value, line_number}); + } + if (stream.bad()) { + fail(ErrorCode::FileRead, path, std::nullopt, std::nullopt, + "could not read file"); + } + return entries; +} + +const Entry &required(const Entries &entries, const char *key, + const std::filesystem::path &path) { + const auto entry = entries.find(key); + if (entry == entries.end()) { + fail(ErrorCode::ConfigInvalid, path, std::nullopt, key, + "missing required key"); + } + return entry->second; +} + +std::uint64_t parse_unsigned(const Entry &entry, const char *field, + const std::filesystem::path &path) { + std::uint64_t value = 0; + const auto begin = entry.value.data(); + const auto end = begin + entry.value.size(); + const auto result = std::from_chars(begin, end, value); + if (entry.value.empty() || entry.value.front() == '-' || + result.ec != std::errc() || result.ptr != end) { + fail(ErrorCode::ConfigInvalid, path, entry.line, field, + "expected integer, got '" + entry.value + "'"); + } + return value; +} + +std::uint32_t parse_uint32(const Entry &entry, const char *field, + const std::filesystem::path &path) { + const auto value = parse_unsigned(entry, field, path); + if (value > UINT32_MAX) { + fail(ErrorCode::ConfigInvalid, path, entry.line, field, + "expected integer <= 4294967295, got '" + entry.value + "'"); + } + return static_cast(value); +} + +double parse_double(const Entry &entry, const char *field, + const std::filesystem::path &path) { + errno = 0; + char *end = nullptr; + const auto value = std::strtod(entry.value.c_str(), &end); + if (entry.value.empty() || + end != entry.value.c_str() + entry.value.size() || errno == ERANGE || + !std::isfinite(value)) { + fail(ErrorCode::ConfigInvalid, path, entry.line, field, + "expected finite number, got '" + entry.value + "'"); + } + return value; +} + +OptionType parse_option_type(const Entry &entry, + const std::filesystem::path &path) { + if (entry.value == "european_call") { + return OptionType::EuropeanCall; + } + if (entry.value == "asian_call") { + return OptionType::AsianArithmeticCall; + } + if (entry.value == "european_put" || entry.value == "barrier_call") { + fail(ErrorCode::UnsupportedFeature, path, entry.line, "option_type", + "option type '" + entry.value + "' is not supported in P0"); + } + fail(ErrorCode::ConfigInvalid, path, entry.line, "option_type", + "invalid option type '" + entry.value + "'"); +} + +Precision parse_precision(const Entry &entry, + const std::filesystem::path &path) { + if (entry.value == "fp32") { + return Precision::Fp32; + } + if (entry.value == "fp64") { + return Precision::Fp64; + } + fail(ErrorCode::ConfigInvalid, path, entry.line, "precision", + "invalid precision '" + entry.value + "'"); +} + +VarianceReduction parse_variance_reduction(const Entry &entry, + const std::filesystem::path &path) { + if (entry.value == "none") { + return VarianceReduction::None; + } + if (entry.value == "antithetic" || entry.value == "control_variate") { + fail(ErrorCode::UnsupportedFeature, path, entry.line, + "variance_reduction", + "variance reduction '" + entry.value + "' is not supported in P0"); + } + fail(ErrorCode::ConfigInvalid, path, entry.line, "variance_reduction", + "invalid variance reduction '" + entry.value + "'"); +} + +RngType parse_rng(const Entry &entry, const std::filesystem::path &path) { + if (entry.value == "curand_philox") { + return RngType::CurandPhilox; + } + fail(ErrorCode::ConfigInvalid, path, entry.line, "rng", + "invalid rng '" + entry.value + "'"); +} + +void validate_positive(double value, const Entry &entry, const char *field, + const std::filesystem::path &path) { + if (value <= 0.0) { + fail(ErrorCode::ConfigInvalid, path, entry.line, field, + "expected number > 0, got '" + entry.value + "'"); + } +} + +} // namespace + +ConfigError::ConfigError(ErrorCode code, std::filesystem::path source_file, + std::optional line, + std::optional field, std::string reason) + : std::runtime_error(reason), code_(code), + source_file_(std::move(source_file)), line_(std::move(line)), + field_(std::move(field)), reason_(std::move(reason)) {} + +ErrorCode ConfigError::code() const noexcept { return code_; } + +const std::filesystem::path &ConfigError::source_file() const noexcept { + return source_file_; +} + +const std::optional &ConfigError::line() const noexcept { + return line_; +} + +const std::optional &ConfigError::field() const noexcept { + return field_; +} + +const std::string &ConfigError::reason() const noexcept { return reason_; } + +std::string ConfigError::render() const { + std::ostringstream stream; + stream << "ERROR [" << code_name(code_) << "] " + << source_file_.filename().string(); + if (line_) { + stream << ':' << *line_; + } + if (field_) { + stream << " field '" << *field_ << "'"; + } + stream << ": " << reason_; + return stream.str(); +} + +int exit_code_for(ErrorCode code) noexcept { + switch (code) { + case ErrorCode::ConfigInvalid: + return 2; + case ErrorCode::FileRead: + return 3; + case ErrorCode::UnsupportedFeature: + return 6; + } + return 1; +} + +RunConfig load_run_config(const std::filesystem::path &option_file, + const std::filesystem::path &simulation_file, + std::filesystem::path output_dir, bool run_cpu, + bool run_gpu) { + const Entries option_entries = parse_ini( + option_file, {"option_type", "spot", "strike", "risk_free_rate", + "volatility", "maturity", "barrier"}); + const Entries simulation_entries = + parse_ini(simulation_file, {"num_paths", "num_steps", "seed", "rng", + "variance_reduction", "precision", + "block_size", "batch_size"}); + + const auto &option_type_entry = + required(option_entries, "option_type", option_file); + const auto &spot_entry = required(option_entries, "spot", option_file); + const auto &strike_entry = required(option_entries, "strike", option_file); + const auto &rate_entry = + required(option_entries, "risk_free_rate", option_file); + const auto &volatility_entry = + required(option_entries, "volatility", option_file); + const auto &maturity_entry = + required(option_entries, "maturity", option_file); + + OptionParams option{ + parse_option_type(option_type_entry, option_file), + parse_double(spot_entry, "spot", option_file), + parse_double(strike_entry, "strike", option_file), + parse_double(rate_entry, "risk_free_rate", option_file), + parse_double(volatility_entry, "volatility", option_file), + parse_double(maturity_entry, "maturity", option_file), + std::nullopt}; + if (const auto barrier = option_entries.find("barrier"); + barrier != option_entries.end()) { + fail(ErrorCode::ConfigInvalid, option_file, barrier->second.line, + "barrier", "barrier is not supported for P0 option types"); + } + + validate_positive(option.spot, spot_entry, "spot", option_file); + validate_positive(option.strike, strike_entry, "strike", option_file); + validate_positive(option.maturity, maturity_entry, "maturity", option_file); + if (option.volatility < 0.0) { + fail(ErrorCode::ConfigInvalid, option_file, volatility_entry.line, + "volatility", + "expected number >= 0, got '" + volatility_entry.value + "'"); + } + + const auto &num_paths_entry = + required(simulation_entries, "num_paths", simulation_file); + const auto &num_steps_entry = + required(simulation_entries, "num_steps", simulation_file); + const auto &seed_entry = + required(simulation_entries, "seed", simulation_file); + const auto &rng_entry = + required(simulation_entries, "rng", simulation_file); + const auto &variance_entry = + required(simulation_entries, "variance_reduction", simulation_file); + const auto &precision_entry = + required(simulation_entries, "precision", simulation_file); + const auto &block_size_entry = + required(simulation_entries, "block_size", simulation_file); + const auto &batch_size_entry = + required(simulation_entries, "batch_size", simulation_file); + + SimulationParams simulation{ + parse_unsigned(num_paths_entry, "num_paths", simulation_file), + parse_uint32(num_steps_entry, "num_steps", simulation_file), + parse_unsigned(seed_entry, "seed", simulation_file), + parse_rng(rng_entry, simulation_file), + parse_precision(precision_entry, simulation_file), + parse_variance_reduction(variance_entry, simulation_file), + parse_uint32(block_size_entry, "block_size", simulation_file), + parse_unsigned(batch_size_entry, "batch_size", simulation_file)}; + if (simulation.num_paths == 0) { + fail(ErrorCode::ConfigInvalid, simulation_file, num_paths_entry.line, + "num_paths", + "expected integer >= 1, got '" + num_paths_entry.value + "'"); + } + if (option.type == OptionType::AsianArithmeticCall && + simulation.num_steps == 0) { + fail(ErrorCode::ConfigInvalid, simulation_file, num_steps_entry.line, + "num_steps", + "expected integer >= 1, got '" + num_steps_entry.value + "'"); + } + if (simulation.block_size == 0) { + fail(ErrorCode::ConfigInvalid, simulation_file, block_size_entry.line, + "block_size", + "expected integer >= 1, got '" + block_size_entry.value + "'"); + } + + return RunConfig{option, simulation, std::move(output_dir), run_cpu, + run_gpu}; +} + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/cpu_pricer.cpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/cpu_pricer.cpp" new file mode 100644 index 00000000..1265e3c2 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/cpu_pricer.cpp" @@ -0,0 +1,91 @@ +// 使用 CPU 单线程生成路径，作为 GPU 结果的统计基线。 +#include "pricer/cpu_pricer.hpp" + +#include +#include +#include +#include + +#include "pricer/payoff.hpp" + +namespace pricer { +namespace { + +// CPU 与 GPU 在启动前执行同一类基础请求校验。 +void validate_request(const OptionParams &option, + const SimulationParams &simulation) { + if (simulation.num_paths == 0U) { + throw std::invalid_argument("CPU pricing requires at least one path"); + } + if (option.type != OptionType::EuropeanCall && + option.type != OptionType::AsianArithmeticCall) { + throw std::domain_error("CPU pricing supports only P0 call options"); + } + if (option.type == OptionType::AsianArithmeticCall && + simulation.num_steps == 0U) { + throw std::domain_error("Asian CPU pricing requires at least one step"); + } +} + +// 将单条路径收益折叠进原始矩，避免保存整条收益数组。 +void accumulate(RawMoments &moments, double payoff) { + // 不保存全部路径收益：只保留 sum 和 sum_squares 已足够在最后算均值、方差和 + // SE。 + ++moments.count; + moments.sum += payoff; + moments.sum_squares += payoff * payoff; +} + +} // namespace + +// 按期权类型选择一次终值采样或多步路径采样。 +CpuPricingRun CpuMonteCarloPricer::price(const OptionParams &option, + const SimulationParams &simulation) { + validate_request(option, simulation); + + const auto start = std::chrono::steady_clock::now(); + // 每次 price 都从同一个 seed 重新开始，保证同一 CPU + // 二进制上的重复运行可复现。 + std::mt19937_64 engine(simulation.seed); + std::normal_distribution normal(0.0, 1.0); + RawMoments moments{0U, 0.0, 0.0}; + + if (option.type == OptionType::EuropeanCall) { + // 欧式期权只关心到期价，所以一条路径只需要一个标准正态随机数和一个 GBM + // 步。 + const auto step = make_gbm_step_constants( + option.maturity, option.risk_free_rate, option.volatility, 1U); + for (std::uint64_t path = 0; path < simulation.num_paths; ++path) { + const double terminal_spot = + option.spot * + std::exp(step.drift + step.diffusion * normal(engine)); + accumulate(moments, + european_call_payoff(terminal_spot, option.strike, + step.discount)); + } + } else { + // 亚式期权关心路径上的平均价，因此一条路径必须逐步演化并累计每个监控价。 + const auto step = + make_gbm_step_constants(option.maturity, option.risk_free_rate, + option.volatility, simulation.num_steps); + for (std::uint64_t path = 0; path < simulation.num_paths; ++path) { + double spot = option.spot; + double running_sum = 0.0; // 不包含初始 S0，见 Asian payoff 的定义。 + for (std::uint32_t index = 0; index < simulation.num_steps; + ++index) { + spot *= std::exp(step.drift + step.diffusion * normal(engine)); + running_sum += spot; + } + accumulate(moments, asian_arithmetic_call_payoff( + running_sum, simulation.num_steps, + option.strike, step.discount)); + } + } + + const auto stop = std::chrono::steady_clock::now(); + const double compute_runtime_ms = + std::chrono::duration(stop - start).count(); + return {moments, compute_runtime_ms}; +} + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/cuda_backend.cu" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/cuda_backend.cu" new file mode 100644 index 00000000..781b1f56 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/cuda_backend.cu" @@ -0,0 +1,533 @@ +// 在 GPU 上生成蒙特卡洛路径，并归约收益的一阶、二阶矩。 +#include "pricer/cuda_pricer.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#ifdef PRICER_ENABLE_NVTX +#include +#endif + +namespace pricer { +namespace { + +[[noreturn]] void throw_cuda_error(const char *api, cudaError_t status, + const char *source_file, int source_line) { + throw CudaError(api, static_cast(status), cudaGetErrorString(status), + source_file, source_line); +} + +void check_cuda(cudaError_t status, const char *api, const char *source_file, + int source_line) { + if (status != cudaSuccess) { + throw_cuda_error(api, status, source_file, source_line); + } +} + +#define PRICER_CUDA_CHECK(api_name, expression) \ + ::pricer::check_cuda((expression), (api_name), __FILE__, __LINE__) + +class DeviceBuffer { + public: + DeviceBuffer() = default; + + explicit DeviceBuffer(std::size_t bytes) { allocate(bytes); } + + DeviceBuffer(const DeviceBuffer &) = delete; + DeviceBuffer &operator=(const DeviceBuffer &) = delete; + + ~DeviceBuffer() noexcept { + if (pointer_ != nullptr) { + const cudaError_t status = cudaFree(pointer_); + if (status != cudaSuccess) { + // Destructors cannot replace an active exception. Normal-path + // cleanup uses release(), which reports a structured error. + } + } + } + + void allocate(std::size_t bytes) { + if (bytes == 0U) { + throw std::invalid_argument( + "CUDA allocation size must be positive"); + } + PRICER_CUDA_CHECK("cudaMalloc", cudaMalloc(&pointer_, bytes)); + } + + void *get() const noexcept { return pointer_; } + + void release() { + if (pointer_ != nullptr) { + void *released = pointer_; + pointer_ = nullptr; + PRICER_CUDA_CHECK("cudaFree", cudaFree(released)); + } + } + + private: + void *pointer_ = nullptr; +}; + +class CudaEvent { + public: + CudaEvent() { + PRICER_CUDA_CHECK("cudaEventCreate", cudaEventCreate(&event_)); + } + + CudaEvent(const CudaEvent &) = delete; + CudaEvent &operator=(const CudaEvent &) = delete; + + ~CudaEvent() noexcept { + if (event_ != nullptr) { + const cudaError_t status = cudaEventDestroy(event_); + if (status != cudaSuccess) { + // See DeviceBuffer::~DeviceBuffer(). + } + } + } + + void record() { + PRICER_CUDA_CHECK("cudaEventRecord", cudaEventRecord(event_)); + } + + void synchronize() { + PRICER_CUDA_CHECK("cudaEventSynchronize", cudaEventSynchronize(event_)); + } + + float elapsed_since(const CudaEvent &start) const { + float milliseconds = 0.0F; + PRICER_CUDA_CHECK( + "cudaEventElapsedTime", + cudaEventElapsedTime(&milliseconds, start.event_, event_)); + return milliseconds; + } + + void release() { + if (event_ != nullptr) { + cudaEvent_t released = event_; + event_ = nullptr; + PRICER_CUDA_CHECK("cudaEventDestroy", cudaEventDestroy(released)); + } + } + + private: + cudaEvent_t event_ = nullptr; +}; + +class NvtxRange { + public: + explicit NvtxRange(const char *label) { +#ifdef PRICER_ENABLE_NVTX + nvtxRangePushA(label); +#else + static_cast(label); +#endif + } + + NvtxRange(const NvtxRange &) = delete; + NvtxRange &operator=(const NvtxRange &) = delete; + + ~NvtxRange() noexcept { +#ifdef PRICER_ENABLE_NVTX + nvtxRangePop(); +#endif + } +}; + +struct Square { + __host__ __device__ double operator()(const double &value) const { + return value * value; + } +}; + +template struct PathMath; + +template <> struct PathMath { + __device__ static float normal(curandStatePhilox4_32_10_t *state) { + return curand_normal(state); + } + + __device__ static float exponential(float value) { return expf(value); } + + __device__ static float maximum(float left, float right) { + return fmaxf(left, right); + } +}; + +template <> struct PathMath { + __device__ static double normal(curandStatePhilox4_32_10_t *state) { + return curand_normal_double(state); + } + + __device__ static double exponential(double value) { return exp(value); } + + __device__ static double maximum(double left, double right) { + return fmax(left, right); + } +}; + +template +__global__ void +european_payoff_kernel(double *payoffs, std::uint64_t batch_paths, + std::uint64_t batch_offset, std::uint64_t seed, + Real spot, Real strike, Real drift, Real diffusion, + Real discount) { + const std::uint64_t thread = + static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const std::uint64_t stride = + static_cast(gridDim.x) * blockDim.x; + for (std::uint64_t local_path = thread; local_path < batch_paths; + local_path += stride) { + // 一个 CUDA 线程负责一条路径。batch_offset + local_path + // 是全局路径编号，让分批与不分批时也取到同一条 Philox + // 随机子序列，而不会重复抽样。 + curandStatePhilox4_32_10_t state; + curand_init(static_cast(seed), + static_cast(batch_offset + local_path), + 0ULL, &state); + const Real normal = PathMath::normal(&state); + const Real terminal_spot = + spot * PathMath::exponential(drift + diffusion * normal); + const Real payoff = + discount * PathMath::maximum(terminal_spot - strike, Real{0}); + payoffs[local_path] = static_cast(payoff); + } +} + +template +__global__ void asian_payoff_kernel(double *payoffs, std::uint64_t batch_paths, + std::uint64_t batch_offset, + std::uint64_t seed, Real initial_spot, + Real strike, Real drift, Real diffusion, + Real discount, std::uint32_t num_steps) { + const std::uint64_t thread = + static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + const std::uint64_t stride = + static_cast(gridDim.x) * blockDim.x; + for (std::uint64_t local_path = thread; local_path < batch_paths; + local_path += stride) { + // 亚式版本与 CPU 的双层循环一一对应：外层是路径，内层是时间监控点。 + curandStatePhilox4_32_10_t state; + curand_init(static_cast(seed), + static_cast(batch_offset + local_path), + 0ULL, &state); + Real spot = initial_spot; + Real running_sum = Real{0}; + for (std::uint32_t step = 0; step < num_steps; ++step) { + spot *= PathMath::exponential( + drift + diffusion * PathMath::normal(&state)); + running_sum += spot; + } + const Real average = running_sum / static_cast(num_steps); + const Real payoff = + discount * PathMath::maximum(average - strike, Real{0}); + payoffs[local_path] = static_cast(payoff); + } +} + +template +void launch_payoff_kernel(const OptionParams &option, + const SimulationParams &simulation, + double *payoff_pointer, std::uint64_t batch_paths, + std::uint64_t batch_offset, std::uint64_t grid, + double drift, double diffusion, double discount) { + if (option.type == OptionType::EuropeanCall) { + european_payoff_kernel + <<(grid), simulation.block_size>>>( + payoff_pointer, batch_paths, batch_offset, simulation.seed, + static_cast(option.spot), + static_cast(option.strike), static_cast(drift), + static_cast(diffusion), static_cast(discount)); + } else { + asian_payoff_kernel + <<(grid), simulation.block_size>>>( + payoff_pointer, batch_paths, batch_offset, simulation.seed, + static_cast(option.spot), + static_cast(option.strike), static_cast(drift), + static_cast(diffusion), static_cast(discount), + simulation.num_steps); + } +} + +std::size_t checked_add(std::size_t left, std::size_t right) { + if (right > std::numeric_limits::max() - left) { + throw std::overflow_error("CUDA batch byte count overflow"); + } + return left + right; +} + +std::size_t payoff_bytes(std::uint64_t batch_paths) { + if (batch_paths > + std::numeric_limits::max() / sizeof(double)) { + throw std::overflow_error("CUDA payoff byte count overflow"); + } + return static_cast(batch_paths) * sizeof(double); +} + +std::size_t reduction_temp_bytes(std::uint64_t batch_paths) { + std::size_t sum_bytes = 0U; + PRICER_CUDA_CHECK( + "cub::DeviceReduce::Sum(query payoff)", + cub::DeviceReduce::Sum(nullptr, sum_bytes, + static_cast(nullptr), + static_cast(nullptr), batch_paths)); + + std::size_t square_bytes = 0U; + const auto squares = thrust::make_transform_iterator( + static_cast(nullptr), Square{}); + PRICER_CUDA_CHECK("cub::DeviceReduce::Sum(query square)", + cub::DeviceReduce::Sum(nullptr, square_bytes, squares, + static_cast(nullptr), + batch_paths)); + return std::max(sum_bytes, square_bytes); +} + +struct BatchStorage { + std::uint64_t paths; + std::size_t payoff_bytes; + std::size_t temporary_bytes; + std::size_t total_bytes; +}; + +BatchStorage describe_storage(std::uint64_t batch_paths) { + const std::size_t payoffs = payoff_bytes(batch_paths); + const std::size_t temporary = reduction_temp_bytes(batch_paths); + std::size_t total = checked_add(payoffs, temporary); + total = checked_add(total, 2U * sizeof(double)); + return {batch_paths, payoffs, temporary, total}; +} + +BatchStorage choose_batch(const SimulationParams &simulation, + std::size_t free_memory_bytes) { + if (simulation.batch_size > 0U) { + const auto storage = describe_storage( + std::min(simulation.batch_size, simulation.num_paths)); + if (storage.total_bytes > free_memory_bytes) { + throw_cuda_error("batch_memory", cudaErrorMemoryAllocation, + __FILE__, __LINE__); + } + return storage; + } + + // 自动模式只使用约 80% 的空闲显存，给 CUDA 运行时和显示任务保留余量。 + const std::size_t budget = (free_memory_bytes / 10U) * 8U; + const std::uint64_t maximum_by_size = static_cast( + std::numeric_limits::max() / sizeof(double)); + std::uint64_t low = 0U; + std::uint64_t high = std::min(simulation.num_paths, maximum_by_size); + while (low < high) { + const std::uint64_t middle = low + (high - low + 1U) / 2U; + if (describe_storage(middle).total_bytes <= budget) { + low = middle; + } else { + high = middle - 1U; + } + } + if (low == 0U) { + throw_cuda_error("automatic_batch_memory", cudaErrorMemoryAllocation, + __FILE__, __LINE__); + } + if (low > simulation.block_size) { + low = (low / simulation.block_size) * simulation.block_size; + } + low = std::max(low, 1U); + return describe_storage(low); +} + +void validate_request(const OptionParams &option, + const SimulationParams &simulation) { + if (simulation.num_paths == 0U) { + throw std::invalid_argument("CUDA pricing requires at least one path"); + } + if (option.type != OptionType::EuropeanCall && + option.type != OptionType::AsianArithmeticCall) { + throw std::domain_error("CUDA pricing supports only P0 call options"); + } + if (option.type == OptionType::AsianArithmeticCall && + simulation.num_steps == 0U) { + throw std::domain_error( + "Asian CUDA pricing requires at least one step"); + } +} + +std::uint64_t grid_size_for(std::uint64_t paths, std::uint32_t block_size, + int maximum_grid_size) { + const std::uint64_t blocks = + paths / block_size + (paths % block_size == 0U ? 0U : 1U); + return std::min(blocks, static_cast(maximum_grid_size)); +} + +} // namespace + +CudaError::CudaError(std::string api, int code, std::string cuda_text, + std::string source_file, int source_line) + : std::runtime_error(api + " failed with CUDA error " + + std::to_string(code) + ": " + cuda_text), + api_(std::move(api)), code_(code), cuda_text_(std::move(cuda_text)), + source_file_(std::move(source_file)), source_line_(source_line) {} + +const std::string &CudaError::api() const noexcept { return api_; } +int CudaError::code() const noexcept { return code_; } +const std::string &CudaError::cuda_text() const noexcept { return cuda_text_; } +const std::string &CudaError::source_file() const noexcept { + return source_file_; +} +int CudaError::source_line() const noexcept { return source_line_; } + +CudaPricingRun CudaMonteCarloPricer::price(const OptionParams &option, + const SimulationParams &simulation, + int device_id) { + validate_request(option, simulation); + + int device_count = 0; + PRICER_CUDA_CHECK("cudaGetDeviceCount", cudaGetDeviceCount(&device_count)); + if (device_id < 0 || device_id >= device_count) { + throw_cuda_error("device_id", cudaErrorInvalidDevice, __FILE__, + __LINE__); + } + PRICER_CUDA_CHECK("cudaSetDevice", cudaSetDevice(device_id)); + + cudaDeviceProp properties{}; + PRICER_CUDA_CHECK("cudaGetDeviceProperties", + cudaGetDeviceProperties(&properties, device_id)); + if (simulation.block_size == 0U || + simulation.block_size > + static_cast(properties.maxThreadsPerBlock)) { + throw std::invalid_argument( + "CUDA block size must be within the selected device limit"); + } + if (properties.maxGridSize[0] <= 0) { + throw std::domain_error("CUDA device reports no usable x-grid size"); + } + + std::size_t free_memory_bytes = 0U; + std::size_t total_memory_bytes = 0U; + PRICER_CUDA_CHECK("cudaMemGetInfo", + cudaMemGetInfo(&free_memory_bytes, &total_memory_bytes)); + // 批处理只限制“同时放在显存里的路径数”；所有批次的统计量随后在主机合并。 + const BatchStorage storage = choose_batch(simulation, free_memory_bytes); + const std::uint64_t representative_grid = grid_size_for( + storage.paths, simulation.block_size, properties.maxGridSize[0]); + const double discount = std::exp(-option.risk_free_rate * option.maturity); + const double step_maturity = + option.type == OptionType::EuropeanCall + ? option.maturity + : option.maturity / static_cast(simulation.num_steps); + const double drift = + (option.risk_free_rate - 0.5 * option.volatility * option.volatility) * + step_maturity; + const double diffusion = option.volatility * std::sqrt(step_maturity); + + DeviceBuffer payoffs(storage.payoff_bytes); + DeviceBuffer reduced_sum(sizeof(double)); + DeviceBuffer reduced_sum_squares(sizeof(double)); + DeviceBuffer temporary(storage.temporary_bytes); + CudaEvent compute_start; + CudaEvent compute_stop; + CudaEvent reduction_start; + CudaEvent reduction_stop; + + auto *payoff_pointer = static_cast(payoffs.get()); + auto *sum_pointer = static_cast(reduced_sum.get()); + auto *sum_squares_pointer = + static_cast(reduced_sum_squares.get()); + + RawMoments moments{0U, 0.0, 0.0}; + double gpu_compute_ms = 0.0; + double reduction_ms = 0.0; + for (std::uint64_t batch_offset = 0U; + batch_offset < simulation.num_paths;) { + const std::uint64_t batch_paths = + std::min(storage.paths, simulation.num_paths - batch_offset); + const std::uint64_t grid = grid_size_for( + batch_paths, simulation.block_size, properties.maxGridSize[0]); + + // 这个循环是一批完整 GPU 工作：生成并计算每条路径收益，再归约为两个数。 + compute_start.record(); + { + NvtxRange range("simulate_paths"); + if (simulation.precision == Precision::Fp32) { + launch_payoff_kernel(option, simulation, payoff_pointer, + batch_paths, batch_offset, grid, + drift, diffusion, discount); + } else { + launch_payoff_kernel(option, simulation, payoff_pointer, + batch_paths, batch_offset, grid, + drift, diffusion, discount); + } + } + PRICER_CUDA_CHECK("cudaGetLastError", cudaGetLastError()); + + reduction_start.record(); + { + NvtxRange range("reduce_moments"); + std::size_t temporary_bytes = storage.temporary_bytes; + PRICER_CUDA_CHECK("cub::DeviceReduce::Sum(payoff)", + cub::DeviceReduce::Sum( + temporary.get(), temporary_bytes, + payoff_pointer, sum_pointer, batch_paths)); + const auto squares = + thrust::make_transform_iterator(payoff_pointer, Square{}); + temporary_bytes = storage.temporary_bytes; + PRICER_CUDA_CHECK("cub::DeviceReduce::Sum(square)", + cub::DeviceReduce::Sum( + temporary.get(), temporary_bytes, squares, + sum_squares_pointer, batch_paths)); + } + reduction_stop.record(); + compute_stop.record(); + compute_stop.synchronize(); + gpu_compute_ms += compute_stop.elapsed_since(compute_start); + reduction_ms += reduction_stop.elapsed_since(reduction_start); + + { + NvtxRange range("copy_batch_moments"); + double host_sum = 0.0; + double host_sum_squares = 0.0; + PRICER_CUDA_CHECK("cudaMemcpy(sum)", + cudaMemcpy(&host_sum, sum_pointer, + sizeof(host_sum), + cudaMemcpyDeviceToHost)); + PRICER_CUDA_CHECK("cudaMemcpy(sum_squares)", + cudaMemcpy(&host_sum_squares, sum_squares_pointer, + sizeof(host_sum_squares), + cudaMemcpyDeviceToHost)); + moments = merge_raw_moments( + moments, {batch_paths, host_sum, host_sum_squares}); + } + // 每批只回传两个标量，随后与先前批次在主机合并。 + batch_offset += batch_paths; + } + + CudaDeviceInfo device{device_id, + properties.name, + properties.major, + properties.minor, + total_memory_bytes, + free_memory_bytes, + properties.maxThreadsPerBlock}; + reduction_stop.release(); + reduction_start.release(); + compute_stop.release(); + compute_start.release(); + temporary.release(); + reduced_sum_squares.release(); + reduced_sum.release(); + payoffs.release(); + + return {moments, gpu_compute_ms, reduction_ms, std::nullopt, + storage.paths, representative_grid, std::move(device)}; +} + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/main.cpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/main.cpp" new file mode 100644 index 00000000..1068039f --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/main.cpp" @@ -0,0 +1,460 @@ +// 协调配置、CPU/GPU 定价、统计汇总和结果输出的命令行入口。 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "pricer/black_scholes.hpp" +#include "pricer/config.hpp" +#include "pricer/cpu_pricer.hpp" +#include "pricer/cuda_pricer.hpp" +#include "pricer/measurement.hpp" +#include "pricer/output.hpp" +#include "pricer/statistics.hpp" +#include "pricer/version.hpp" + +#ifndef PRICER_BUILD_TYPE +#define PRICER_BUILD_TYPE "unknown" +#endif + +#ifndef PRICER_GIT_COMMIT +#define PRICER_GIT_COMMIT "unknown" +#endif + +namespace { + +class CliError : public std::runtime_error { + public: + using std::runtime_error::runtime_error; +}; + +class NumericResultError : public std::runtime_error { + public: + using std::runtime_error::runtime_error; +}; + +class CudaRunError : public std::runtime_error { + public: + using std::runtime_error::runtime_error; +}; + +struct Arguments { + std::filesystem::path option_path; + std::filesystem::path simulation_path; + std::filesystem::path output_dir = "results"; + std::string backend = "both"; + int device = 0; + std::uint32_t repetitions = 1U; + std::uint32_t warmup = 1U; +}; + +struct BackendRun { + pricer::PricingResult result; + pricer::PerformanceData performance; + pricer::EnvironmentData environment; + std::uint64_t actual_batch_size; +}; + +void print_help() { + std::cout + << "Usage: pricer_cli --option --simulation [options]\n" + << " --output-dir Output directory (default: results)\n" + << " --backend Backend selection (default: both)\n" + << " --device CUDA device (default: 0)\n" + << " --repetitions Formal measurements (default: 1)\n" + << " --warmup Unmeasured GPU warmups (default: 1)\n" + << " --help Show this help\n" + << " --version Show version\n"; +} + +template +Integer parse_integer(std::string_view text, std::string_view flag, + bool allow_zero) { + std::uint64_t parsed = 0U; + const auto result = + std::from_chars(text.data(), text.data() + text.size(), parsed); + if (text.empty() || text.front() == '-' || result.ec != std::errc() || + result.ptr != text.data() + text.size() || + (!allow_zero && parsed == 0U) || + parsed > + static_cast(std::numeric_limits::max())) { + throw CliError(std::string(flag) + " has an invalid integer value"); + } + return static_cast(parsed); +} + +Arguments parse_arguments(int argc, char **argv) { + Arguments arguments; + std::unordered_set seen; + for (int index = 1; index < argc; ++index) { + const std::string flag = argv[index]; + if (flag.empty() || flag.front() != '-') { + throw CliError("unexpected positional argument '" + flag + "'"); + } + if (flag == "--help" || flag == "--version") { + throw CliError(flag + " must be used alone"); + } + if (flag != "--option" && flag != "--simulation" && + flag != "--output-dir" && flag != "--backend" && + flag != "--device" && flag != "--repetitions" && + flag != "--warmup") { + throw CliError("unknown flag '" + flag + "'"); + } + if (!seen.insert(flag).second) { + throw CliError("duplicate flag '" + flag + "'"); + } + if (++index >= argc) { + throw CliError("missing value for '" + flag + "'"); + } + const std::string value = argv[index]; + if (flag == "--option") { + arguments.option_path = value; + } else if (flag == "--simulation") { + arguments.simulation_path = value; + } else if (flag == "--output-dir") { + arguments.output_dir = value; + } else if (flag == "--backend") { + if (value != "cpu" && value != "gpu" && value != "both") { + throw CliError("--backend must be cpu, gpu, or both"); + } + arguments.backend = value; + } else if (flag == "--device") { + arguments.device = parse_integer(value, flag, true); + } else if (flag == "--repetitions") { + arguments.repetitions = + parse_integer(value, flag, false); + } else { + arguments.warmup = parse_integer(value, flag, true); + } + } + if (arguments.option_path.empty()) { + throw CliError("missing required --option"); + } + if (arguments.simulation_path.empty()) { + throw CliError("missing required --simulation"); + } + if (arguments.output_dir.empty()) { + throw CliError("--output-dir must not be empty"); + } + return arguments; +} + +std::optional reference_price(const pricer::OptionParams &option) { + // 只有欧式看涨有这里可用的解析参考价；亚式结果以 CPU/GPU 统计一致性验证。 + if (option.type == pricer::OptionType::EuropeanCall) { + return pricer::black_scholes_call(option); + } + return std::nullopt; +} + +void require_finite_nonnegative(double value, const char *name) { + if (!std::isfinite(value) || value < 0.0) { + throw NumericResultError(std::string(name) + + " must be finite and non-negative"); + } +} + +void validate_result(const pricer::PricingResult &result, + const pricer::PerformanceData &performance) { + require_finite_nonnegative(result.price, "price"); + for (const auto &[value, name] : + std::vector, const char *>>{ + {result.sample_stddev, "sample stddev"}, + {result.standard_error, "standard error"}, + {result.reference_price, "reference price"}, + {result.absolute_error, "absolute error"}, + {result.relative_error, "relative error"}}) { + if (value) { + require_finite_nonnegative(*value, name); + } + } + if (result.ci_lower && !std::isfinite(*result.ci_lower)) { + throw NumericResultError("confidence interval lower bound is invalid"); + } + if (result.ci_upper && !std::isfinite(*result.ci_upper)) { + throw NumericResultError("confidence interval upper bound is invalid"); + } + if (result.ci_lower && result.ci_upper && + *result.ci_lower > *result.ci_upper) { + throw NumericResultError("confidence interval is inverted"); + } + require_finite_nonnegative(performance.total_runtime_ms, "total runtime"); + require_finite_nonnegative(performance.compute_runtime_ms, + "compute runtime"); + require_finite_nonnegative(performance.paths_per_second, + "paths per second"); + for (const auto &value : + {performance.rng_ms, performance.reduction_ms, + performance.cpu_runtime_ms, performance.speedup_vs_cpu}) { + if (value) { + require_finite_nonnegative(*value, "optional performance metric"); + } + } +} + +double throughput(std::uint64_t paths, double compute_runtime_ms) { + if (compute_runtime_ms <= 0.0) { + throw NumericResultError("compute runtime is zero"); + } + return static_cast(paths) * 1000.0 / compute_runtime_ms; +} + +std::optional cuda_runtime_version() { + int version = 0; + if (cudaRuntimeGetVersion(&version) != cudaSuccess) { + return std::nullopt; + } + return std::to_string(version / 1000) + "." + + std::to_string((version % 1000) / 10); +} + +pricer::EnvironmentData base_environment() { + return {std::nullopt, std::nullopt, cuda_runtime_version(), + std::string(PRICER_BUILD_TYPE), std::string(PRICER_GIT_COMMIT)}; +} + +BackendRun run_cpu(const pricer::RunConfig &config, std::uint32_t repetitions) { + // measure_cpu 负责重复运行和取中位数；定价器本身只负责产生 RawMoments。 + const auto measured = pricer::measure_cpu( + config.option, config.simulation, repetitions, + reference_price(config.option), pricer::CpuMonteCarloPricer::price, + pricer::ResultAnalyzer::analyze); + pricer::PerformanceData performance{ + measured.total_runtime_ms, + measured.compute_runtime_ms, + std::nullopt, + std::nullopt, + throughput(config.simulation.num_paths, measured.compute_runtime_ms), + config.simulation.block_size, + 0U, + repetitions, + "median", + measured.total_runtime_ms, + std::nullopt}; + validate_result(measured.result, performance); + return {measured.result, std::move(performance), base_environment(), + config.simulation.batch_size}; +} + +pricer::CudaPricingRun +price_gpu_once(const pricer::OptionParams &option, + const pricer::SimulationParams &simulation, int device) { + try { + return pricer::CudaMonteCarloPricer::price(option, simulation, device); + } catch (const pricer::CudaError &) { + throw; + } catch (const std::invalid_argument &error) { + throw CudaRunError(error.what()); + } catch (const std::domain_error &error) { + if (std::string_view(error.what()) == + "CUDA device reports no usable x-grid size") { + throw CudaRunError(error.what()); + } + throw; + } catch (const std::overflow_error &error) { + const std::string_view message(error.what()); + if (message == "CUDA batch byte count overflow" || + message == "CUDA payoff byte count overflow") { + throw CudaRunError(error.what()); + } + throw; + } +} + +BackendRun run_gpu(const pricer::RunConfig &config, int device, + std::uint32_t warmup, std::uint32_t repetitions, + std::optional cpu_total_runtime_ms) { + // GPU 的 warmup 不计入正式结果，用来减少首次 CUDA 初始化对计时的干扰。 + const auto measured = + pricer::measure_gpu(config.option, config.simulation, device, warmup, + repetitions, reference_price(config.option), + price_gpu_once, pricer::ResultAnalyzer::analyze); + // speedup_vs_cpu is an end-to-end metric: both operands include the + // backend work required to obtain a complete PricingResult. + const std::optional speedup = + cpu_total_runtime_ms && measured.total_runtime_ms > 0.0 + ? std::optional(*cpu_total_runtime_ms / + measured.total_runtime_ms) + : std::nullopt; + pricer::PerformanceData performance{ + measured.total_runtime_ms, + measured.compute_runtime_ms, + measured.rng_ms, + measured.reduction_ms, + throughput(config.simulation.num_paths, measured.compute_runtime_ms), + config.simulation.block_size, + measured.grid_size, + repetitions, + "median", + cpu_total_runtime_ms, + speedup}; + auto environment = base_environment(); + environment.gpu_name = measured.device.name; + environment.compute_capability = + std::to_string(measured.device.compute_capability_major) + "." + + std::to_string(measured.device.compute_capability_minor); + validate_result(measured.result, performance); + return {measured.result, std::move(performance), std::move(environment), + measured.actual_batch_size}; +} + +std::string option_type_name(pricer::OptionType type) { + return type == pricer::OptionType::EuropeanCall ? "european_call" + : "asian_call"; +} + +std::string absolute_path_string(const std::filesystem::path &path) { + std::error_code error; + const auto absolute = std::filesystem::absolute(path, error); + if (error) { + throw pricer::OutputError("could not resolve absolute path: " + + error.message()); + } + return absolute.string(); +} + +void write_backend_output(const Arguments &arguments, + const pricer::RunConfig &config, + const std::string ×tamp, + const std::string &backend, BackendRun run) { + auto simulation = config.simulation; + simulation.batch_size = run.actual_batch_size; + if (backend == "cpu") { + simulation.precision = pricer::Precision::Fp64; + simulation.rng = pricer::RngType::Mt19937_64; + } + pricer::OutputRecord record{ + timestamp + "-" + option_type_name(config.option.type) + "-" + backend, + timestamp, + backend, + config.option, + simulation, + std::move(run.result), + std::move(run.performance), + std::move(run.environment)}; + const auto written = + pricer::write_output_transaction(arguments.output_dir, record); + record.run_id = written.run_id; + + const auto &result = record.result; + const auto &performance = record.performance; + std::cout << "INFO backend=" << backend << " price=" << result.price; + if (result.standard_error) { + std::cout << " standard_error=" << *result.standard_error; + } else { + std::cout << " standard_error=null"; + } + if (result.ci_lower && result.ci_upper) { + std::cout << " ci95=[" << *result.ci_lower << ',' << *result.ci_upper + << ']'; + } else { + std::cout << " ci95=null"; + } + std::cout << " total_runtime_ms=" << performance.total_runtime_ms + << " compute_runtime_ms=" << performance.compute_runtime_ms + << " paths_per_second=" << performance.paths_per_second; + if (performance.speedup_vs_cpu) { + std::cout << " speedup_vs_cpu=" << *performance.speedup_vs_cpu; + } + std::cout << '\n' + << "INFO result_json=" << absolute_path_string(written.json_path) + << '\n' + << "INFO performance_csv=" + << absolute_path_string(written.csv_path) << '\n'; +} + +int run_cli(int argc, char **argv) { + if (argc == 2 && std::string_view(argv[1]) == "--help") { + print_help(); + return 0; + } + if (argc == 2 && std::string_view(argv[1]) == "--version") { + std::cout << "pricer_cli " << pricer::version() << '\n'; + return 0; + } + // 主流程只编排：解析命令行 -> 读取并校验 INI -> CPU/GPU 定价 -> 写 + // JSON/CSV。 + const Arguments arguments = parse_arguments(argc, argv); + const bool run_cpu_backend = arguments.backend != "gpu"; + const bool run_gpu_backend = arguments.backend != "cpu"; + const auto config = pricer::load_run_config( + arguments.option_path, arguments.simulation_path, arguments.output_dir, + run_cpu_backend, run_gpu_backend); + std::cout << "INFO option=" << absolute_path_string(arguments.option_path) + << " simulation=" + << absolute_path_string(arguments.simulation_path) + << " backend=" << arguments.backend << '\n'; + + std::optional cpu; + if (run_cpu_backend) { + cpu = run_cpu(config, arguments.repetitions); + } + std::optional gpu; + if (run_gpu_backend) { + const std::optional cpu_total_runtime = + cpu ? cpu->performance.cpu_runtime_ms : std::nullopt; + gpu = run_gpu(config, arguments.device, arguments.warmup, + arguments.repetitions, cpu_total_runtime); + } + + const std::string timestamp = pricer::utc_timestamp(); + if (cpu) { + write_backend_output(arguments, config, timestamp, "cpu", + std::move(*cpu)); + } + if (gpu) { + write_backend_output(arguments, config, timestamp, "gpu", + std::move(*gpu)); + } + return 0; +} + +} // namespace + +int main(int argc, char **argv) { + try { + return run_cli(argc, argv); + } catch (const CliError &error) { + std::cerr << "ERROR [CLI_INVALID] " << error.what() << '\n'; + return 2; + } catch (const pricer::ConfigError &error) { + std::cerr << error.render() << '\n'; + return pricer::exit_code_for(error.code()); + } catch (const pricer::OutputError &error) { + std::cerr << "ERROR [FILE_WRITE] " << error.what() << '\n'; + return 3; + } catch (const pricer::CudaError &error) { + std::cerr << "ERROR [CUDA_RUNTIME] " << error.api() << ": " + << error.cuda_text() << '\n'; + return 4; + } catch (const CudaRunError &error) { + std::cerr << "ERROR [CUDA_RUNTIME] " << error.what() << '\n'; + return 4; + } catch (const NumericResultError &error) { + std::cerr << "ERROR [NUMERIC_RESULT] " << error.what() << '\n'; + return 5; + } catch (const std::domain_error &error) { + std::cerr << "ERROR [NUMERIC_RESULT] " << error.what() << '\n'; + return 5; + } catch (const std::overflow_error &error) { + std::cerr << "ERROR [NUMERIC_RESULT] " << error.what() << '\n'; + return 5; + } catch (const std::exception &error) { + std::cerr << "ERROR [INTERNAL] " << error.what() << '\n'; + return 10; + } catch (...) { + std::cerr << "ERROR [INTERNAL] unknown failure\n"; + return 10; + } +} diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/measurement.cpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/measurement.cpp" new file mode 100644 index 00000000..f34e1f27 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/measurement.cpp" @@ -0,0 +1,120 @@ +// 执行预热和重复测量，并以中位数降低偶发波动的影响。 +#include "pricer/measurement.hpp" + +#include +#include +#include +#include + +#include "pricer/output.hpp" + +namespace pricer { +namespace { + +// 同一固定 seed 的重复正式运行必须产生相同统计量。 +void require_same_moments(const RawMoments &left, const RawMoments &right) { + if (left.count != right.count || left.sum != right.sum || + left.sum_squares != right.sum_squares) { + throw std::runtime_error( + "fixed config and seed did not reproduce identical moments"); + } +} + +} // namespace + +CpuMeasurement +// CPU 没有预热阶段；每次测量同时产出可比较的统计结果。 +measure_cpu(const OptionParams &option, const SimulationParams &simulation, + std::uint32_t repetitions, std::optional reference_price, + const CpuPriceFunction &price, const AnalyzeFunction &analyze) { + if (repetitions == 0U) { + throw std::invalid_argument("CPU measurement requires repetitions"); + } + std::vector totals; + std::vector computes; + std::optional moments; + std::optional result; + for (std::uint32_t repetition = 0U; repetition < repetitions; + ++repetition) { + const auto start = std::chrono::steady_clock::now(); + const auto run = price(option, simulation); + if (moments) { + require_same_moments(*moments, run.moments); + } else { + moments = run.moments; + } + auto analyzed = analyze(run.moments, reference_price); + const auto stop = std::chrono::steady_clock::now(); + if (!result) { + result = std::move(analyzed); + } + totals.push_back( + std::chrono::duration(stop - start).count()); + computes.push_back(run.compute_runtime_ms); + } + return {std::move(*result), median(std::move(totals)), + median(std::move(computes))}; +} + +// GPU 预热不计入正式中位数，以隔离上下文初始化的偶发成本。 +GpuMeasurement measure_gpu(const OptionParams &option, + const SimulationParams &simulation, int device, + std::uint32_t warmup, std::uint32_t repetitions, + std::optional reference_price, + const GpuPriceFunction &price, + const AnalyzeFunction &analyze) { + if (repetitions == 0U) { + throw std::invalid_argument("GPU measurement requires repetitions"); + } + for (std::uint32_t index = 0U; index < warmup; ++index) { + static_cast(price(option, simulation, device)); + } + + SimulationParams formal_simulation = simulation; + std::vector totals; + std::vector computes; + std::vector reductions; + std::vector rngs; + std::optional moments; + std::optional result; + std::optional first; + for (std::uint32_t repetition = 0U; repetition < repetitions; + ++repetition) { + const auto start = std::chrono::steady_clock::now(); + auto run = price(option, formal_simulation, device); + if (!first) { + first = run; + if (formal_simulation.batch_size == 0U) { + formal_simulation.batch_size = run.batch_size; + } + } + if (moments) { + require_same_moments(*moments, run.moments); + } else { + moments = run.moments; + } + auto analyzed = analyze(run.moments, reference_price); + const auto stop = std::chrono::steady_clock::now(); + if (!result) { + result = std::move(analyzed); + } + totals.push_back( + std::chrono::duration(stop - start).count()); + computes.push_back(run.gpu_compute_ms); + reductions.push_back(run.reduction_ms); + if (run.rng_ms) { + rngs.push_back(*run.rng_ms); + } + } + return {std::move(*result), + median(std::move(totals)), + median(std::move(computes)), + rngs.empty() ? std::nullopt + : std::optional(median(std::move(rngs))), + median(std::move(reductions)), + first->batch_size, + first->grid_size, + std::move(first->device)}; +} + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/output.cpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/output.cpp" new file mode 100644 index 00000000..1bc5bde6 --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/output.cpp" @@ -0,0 +1,524 @@ +// 将定价和性能结果安全地写为 JSON 与 CSV 文件。 +#include "pricer/output.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef _WIN32 +#include +#else +#include +#include +#include +#include +#include +#endif + +namespace pricer { +namespace { + +constexpr const char *k_csv_header = + "schema_version,run_id,timestamp_utc,option_type,backend,precision," + "num_paths,num_steps,batch_size,seed,block_size,grid_size," + "total_runtime_ms,compute_runtime_ms,rng_ms,reduction_ms," + "paths_per_second,cpu_runtime_ms,speedup_vs_cpu,gpu_name," + "cuda_version,git_commit"; + +const char *option_type_name(OptionType type) { + switch (type) { + case OptionType::EuropeanCall: + return "european_call"; + case OptionType::AsianArithmeticCall: + return "asian_call"; + } + throw OutputError("unknown option type"); +} + +const char *precision_name(Precision precision) { + switch (precision) { + case Precision::Fp32: + return "fp32"; + case Precision::Fp64: + return "fp64"; + } + throw OutputError("unknown precision"); +} + +const char *rng_name(RngType rng) { + switch (rng) { + case RngType::CurandPhilox: + return "curand_philox"; + case RngType::Mt19937_64: + return "mt19937_64"; + } + throw OutputError("unknown RNG"); +} + +const char *variance_reduction_name(VarianceReduction reduction) { + if (reduction == VarianceReduction::None) { + return "none"; + } + throw OutputError("unknown variance reduction"); +} + +template +nlohmann::json optional_json(const std::optional &value) { + return value ? nlohmann::json(*value) : nlohmann::json(nullptr); +} + +std::string numeric(double value) { + if (!std::isfinite(value)) { + return {}; + } + std::ostringstream stream; + stream << std::setprecision(17) << value; + return stream.str(); +} + +std::string optional_numeric(const std::optional &value) { + return value ? numeric(*value) : std::string{}; +} + +std::string csv_escape(const std::optional &value) { + if (!value) { + return {}; + } + if (value->find_first_of(",\"\r\n") == std::string::npos) { + return *value; + } + std::string escaped = "\""; + for (const char character : *value) { + escaped += character; + if (character == '"') { + escaped += '"'; + } + } + escaped += '"'; + return escaped; +} + +std::string csv_escape(const std::string &value) { + return csv_escape(std::optional(value)); +} + +void require_output_directory(const std::filesystem::path &directory) { + std::error_code error; + if (std::filesystem::create_directories(directory, error)) { + return; + } + if (error) { + throw OutputError("could not create output directory: " + + directory.string() + ": " + error.message()); + } + const bool is_directory = std::filesystem::is_directory(directory, error); + if (error) { + throw OutputError("could not inspect output directory: " + + directory.string() + ": " + error.message()); + } + if (!is_directory) { + throw OutputError("output path is not a directory: " + + directory.string()); + } +} + +class OutputLock { + public: + explicit OutputLock(const std::filesystem::path &directory) { + require_output_directory(directory); + path_ = directory / ".pricer-output.lock"; +#ifdef _WIN32 + for (std::uint32_t attempt = 0U; attempt < 5000U; ++attempt) { + handle_ = CreateFileW(path_.c_str(), GENERIC_READ | GENERIC_WRITE, + 0, nullptr, OPEN_ALWAYS, + FILE_ATTRIBUTE_NORMAL, nullptr); + if (handle_ != INVALID_HANDLE_VALUE) { + return; + } + const DWORD error = GetLastError(); + if (error != ERROR_SHARING_VIOLATION && + error != ERROR_LOCK_VIOLATION) { + throw OutputError( + "could not acquire output lock: Windows error " + + std::to_string(error)); + } + std::this_thread::sleep_for(std::chrono::milliseconds(1)); + } + throw OutputError("timed out acquiring output lock: " + path_.string()); +#else + descriptor_ = open(path_.c_str(), O_CREAT | O_RDWR, 0666); + if (descriptor_ < 0) { + throw OutputError("could not open output lock: " + + std::string(std::strerror(errno))); + } + if (flock(descriptor_, LOCK_EX) != 0) { + const std::string reason = std::strerror(errno); + close(descriptor_); + descriptor_ = -1; + throw OutputError("could not acquire output lock: " + reason); + } +#endif + } + + OutputLock(const OutputLock &) = delete; + OutputLock &operator=(const OutputLock &) = delete; + + ~OutputLock() noexcept { +#ifdef _WIN32 + if (handle_ != INVALID_HANDLE_VALUE) { + CloseHandle(handle_); + } +#else + if (descriptor_ >= 0) { + flock(descriptor_, LOCK_UN); + close(descriptor_); + } +#endif + } + + private: + std::filesystem::path path_; +#ifdef _WIN32 + HANDLE handle_ = INVALID_HANDLE_VALUE; +#else + int descriptor_ = -1; +#endif +}; + +std::uint64_t process_id() { +#ifdef _WIN32 + return static_cast(GetCurrentProcessId()); +#else + return static_cast(getpid()); +#endif +} + +bool path_exists(const std::filesystem::path &path, const char *context) { + std::error_code error; + const bool exists = std::filesystem::exists(path, error); + if (error) { + throw OutputError(std::string(context) + ": " + error.message()); + } + return exists; +} + +struct JsonWrite { + std::filesystem::path path; + std::string run_id; +}; + +struct CsvState { + bool existed; + std::uintmax_t size; + bool needs_newline; +}; + +std::string csv_row(const OutputRecord &record) { + std::ostringstream output; + output << "1.0," << csv_escape(record.run_id) << ',' + << csv_escape(record.timestamp_utc) << ',' + << csv_escape(std::string(option_type_name(record.option.type))) + << ',' << csv_escape(record.backend) << ',' + << csv_escape( + std::string(precision_name(record.simulation.precision))) + << ',' << record.simulation.num_paths << ',' + << record.simulation.num_steps << ',' << record.simulation.batch_size + << ',' << record.simulation.seed << ','; + if (record.backend == "gpu") { + output << record.performance.block_size; + } + output << ','; + if (record.backend == "gpu") { + output << record.performance.grid_size; + } + output << ',' << numeric(record.performance.total_runtime_ms) << ',' + << numeric(record.performance.compute_runtime_ms) << ',' + << optional_numeric(record.performance.rng_ms) << ',' + << optional_numeric(record.performance.reduction_ms) << ',' + << numeric(record.performance.paths_per_second) << ',' + << optional_numeric(record.performance.cpu_runtime_ms) << ',' + << optional_numeric(record.performance.speedup_vs_cpu) << ',' + << csv_escape(record.environment.gpu_name) << ',' + << csv_escape(record.environment.cuda_runtime_version) << ',' + << csv_escape(record.environment.git_commit); + return output.str(); +} + +CsvState inspect_csv_unlocked(const std::filesystem::path &path) { + std::error_code error; + const bool existed = std::filesystem::exists(path, error); + if (error) { + throw OutputError("could not inspect performance CSV: " + + error.message()); + } + if (!existed) { + return {false, 0U, false}; + } + const auto size = std::filesystem::file_size(path, error); + if (error) { + throw OutputError("could not inspect performance CSV size: " + + error.message()); + } + if (size == 0U) { + return {true, 0U, false}; + } + + std::ifstream input(path, std::ios::binary); + if (!input) { + throw OutputError("could not open performance CSV for validation: " + + path.string()); + } + std::string header; + if (!std::getline(input, header)) { + throw OutputError("could not read performance CSV header"); + } + if (!header.empty() && header.back() == '\r') { + header.pop_back(); + } + if (header != k_csv_header) { + throw OutputError("performance CSV header does not match schema 1.0"); + } + input.clear(); + input.seekg(-1, std::ios::end); + char last = '\0'; + input.get(last); + if (!input) { + throw OutputError("could not inspect performance CSV tail"); + } + return {true, size, last != '\n'}; +} + +void append_csv_unlocked(const std::filesystem::path &path, + const OutputRecord &record, const CsvState &state) { + std::ofstream output(path, std::ios::binary | std::ios::app); + if (!output) { + throw OutputError("could not open performance CSV for append: " + + path.string()); + } + if (state.size == 0U) { + output << k_csv_header << '\n'; + } else if (state.needs_newline) { + output << '\n'; + } + output << csv_row(record) << '\n'; + output.flush(); + output.close(); + if (!output) { + throw OutputError("could not append performance CSV: " + path.string()); + } +} + +std::filesystem::path +unique_temporary_path(const std::filesystem::path &final_path) { + static std::atomic next_id{0U}; + for (std::uint32_t attempt = 0U; attempt < 100U; ++attempt) { + const auto candidate = std::filesystem::path( + final_path.string() + "." + std::to_string(process_id()) + "." + + std::to_string(next_id.fetch_add(1U)) + ".tmp"); + if (!path_exists(candidate, "could not inspect temporary JSON path")) { + return candidate; + } + } + throw OutputError("could not reserve a unique temporary JSON path"); +} + +JsonWrite write_json_unlocked(const std::filesystem::path &output_dir, + OutputRecord record) { + const std::string base_run_id = + record.run_id.empty() + ? record.timestamp_utc + "-" + + option_type_name(record.option.type) + "-" + record.backend + : record.run_id; + std::filesystem::path final_path = + output_dir / (base_run_id + "-result.json"); + for (std::uint64_t suffix = 1U; + path_exists(final_path, "could not inspect result JSON path"); + ++suffix) { + final_path = output_dir / (base_run_id + "-" + std::to_string(suffix) + + "-result.json"); + } + std::string final_run_id = final_path.stem().string(); + final_run_id.resize(final_run_id.size() - std::string("-result").size()); + record.run_id = final_run_id; + + const auto temporary_path = unique_temporary_path(final_path); + try { + std::ofstream stream(temporary_path, + std::ios::binary | std::ios::trunc); + if (!stream) { + throw OutputError("could not open temporary JSON: " + + temporary_path.string()); + } + stream << make_result_json(record).dump(2) << '\n'; + stream.flush(); + if (!stream) { + throw OutputError("could not write temporary JSON: " + + temporary_path.string()); + } + stream.close(); + if (!stream) { + throw OutputError("could not close temporary JSON: " + + temporary_path.string()); + } + std::error_code error; + std::filesystem::rename(temporary_path, final_path, error); + if (error) { + throw OutputError("could not publish result JSON: " + + error.message()); + } + } catch (...) { + std::error_code ignored; + std::filesystem::remove(temporary_path, ignored); + throw; + } + return {final_path, final_run_id}; +} + +} // namespace + +nlohmann::json make_result_json(const OutputRecord &record) { + nlohmann::json json; + const nlohmann::json confidence_interval = + record.result.ci_lower && record.result.ci_upper + ? nlohmann::json::array( + {*record.result.ci_lower, *record.result.ci_upper}) + : nlohmann::json(nullptr); + json["schema_version"] = "1.0"; + json["run_id"] = record.run_id; + json["status"] = "success"; + json["timestamp_utc"] = record.timestamp_utc; + json["backend"] = record.backend; + json["option"] = {{"type", option_type_name(record.option.type)}, + {"spot", record.option.spot}, + {"strike", record.option.strike}, + {"risk_free_rate", record.option.risk_free_rate}, + {"volatility", record.option.volatility}, + {"maturity", record.option.maturity}, + {"barrier", optional_json(record.option.barrier)}}; + json["simulation"] = { + {"num_paths", record.simulation.num_paths}, + {"num_steps", record.simulation.num_steps}, + {"seed", record.simulation.seed}, + {"rng", rng_name(record.simulation.rng)}, + {"variance_reduction", + variance_reduction_name(record.simulation.variance_reduction)}, + {"precision", precision_name(record.simulation.precision)}, + {"block_size", record.simulation.block_size}, + {"batch_size", record.simulation.batch_size}}; + json["result"] = { + {"price_estimate", record.result.price}, + {"sample_stddev", optional_json(record.result.sample_stddev)}, + {"standard_error", optional_json(record.result.standard_error)}, + {"confidence_level", 0.95}, + {"confidence_interval", confidence_interval}, + {"reference_price", optional_json(record.result.reference_price)}, + {"absolute_error", optional_json(record.result.absolute_error)}, + {"relative_error", optional_json(record.result.relative_error)}}; + json["performance"] = { + {"total_runtime_ms", record.performance.total_runtime_ms}, + {"compute_runtime_ms", record.performance.compute_runtime_ms}, + {"gpu_compute_ms", + record.backend == "gpu" + ? nlohmann::json(record.performance.compute_runtime_ms) + : nlohmann::json(nullptr)}, + {"rng_ms", optional_json(record.performance.rng_ms)}, + {"reduction_ms", optional_json(record.performance.reduction_ms)}, + {"paths_per_second", record.performance.paths_per_second}, + {"block_size", record.backend == "gpu" + ? nlohmann::json(record.performance.block_size) + : nlohmann::json(nullptr)}, + {"grid_size", record.backend == "gpu" + ? nlohmann::json(record.performance.grid_size) + : nlohmann::json(nullptr)}, + {"repetitions", record.performance.repetitions}, + {"aggregation", record.performance.aggregation}, + {"cpu_runtime_ms", optional_json(record.performance.cpu_runtime_ms)}, + {"speedup_vs_cpu", optional_json(record.performance.speedup_vs_cpu)}}; + json["environment"] = { + {"gpu_name", optional_json(record.environment.gpu_name)}, + {"compute_capability", + optional_json(record.environment.compute_capability)}, + {"cuda_runtime_version", + optional_json(record.environment.cuda_runtime_version)}, + {"build_type", optional_json(record.environment.build_type)}, + {"git_commit", optional_json(record.environment.git_commit)}}; + return json; +} + +std::filesystem::path +write_result_json_atomic(const std::filesystem::path &output_dir, + OutputRecord record) { + OutputLock lock(output_dir); + return write_json_unlocked(output_dir, std::move(record)).path; +} + +void append_performance_csv(const std::filesystem::path &path, + const OutputRecord &record) { + const auto directory = path.parent_path().empty() + ? std::filesystem::path(".") + : path.parent_path(); + OutputLock lock(directory); + append_csv_unlocked(path, record, inspect_csv_unlocked(path)); +} + +WrittenOutput write_output_transaction(const std::filesystem::path &output_dir, + OutputRecord record) { + OutputLock lock(output_dir); + const auto csv_path = output_dir / "performance.csv"; + const CsvState csv_state = inspect_csv_unlocked(csv_path); + const JsonWrite json = write_json_unlocked(output_dir, record); + record.run_id = json.run_id; + try { + append_csv_unlocked(csv_path, record, csv_state); + } catch (const std::exception &error) { + std::error_code json_error; + std::filesystem::remove(json.path, json_error); + std::error_code csv_error; + if (csv_state.existed) { + std::filesystem::resize_file(csv_path, csv_state.size, csv_error); + } else { + std::filesystem::remove(csv_path, csv_error); + } + if (json_error || csv_error) { + throw OutputError(std::string(error.what()) + + "; output rollback failed"); + } + throw; + } + return {json.path, csv_path, json.run_id}; +} + +double median(std::vector measurements) { + if (measurements.empty()) { + throw std::invalid_argument("median requires at least one measurement"); + } + std::sort(measurements.begin(), measurements.end()); + const std::size_t middle = measurements.size() / 2U; + if (measurements.size() % 2U != 0U) { + return measurements[middle]; + } + return (measurements[middle - 1U] + measurements[middle]) / 2.0; +} + +std::string utc_timestamp() { + const std::time_t now = + std::chrono::system_clock::to_time_t(std::chrono::system_clock::now()); + std::tm utc{}; +#ifdef _WIN32 + gmtime_s(&utc, &now); +#else + gmtime_r(&now, &utc); +#endif + std::ostringstream stream; + stream << std::put_time(&utc, "%Y%m%dT%H%M%SZ"); + return stream.str(); +} + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/payoff.cpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/payoff.cpp" new file mode 100644 index 00000000..2e62af1e --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/payoff.cpp" @@ -0,0 +1,82 @@ +// 计算 GBM 演化常量以及两类期权的折现收益。 +#include "pricer/payoff.hpp" + +#include +#include +#include +#include + +namespace pricer { +namespace { + +// 所有金融输入先要求为有限值，防止 NaN 在统计阶段扩散。 +void require_finite(double value, const char *name) { + if (!std::isfinite(value)) { + throw std::domain_error(std::string(name) + " must be finite"); + } +} + +void validate_payoff_inputs(double spot_or_sum, double strike, + double discount) { + require_finite(spot_or_sum, "spot or running sum"); + require_finite(strike, "strike"); + require_finite(discount, "discount"); + if (spot_or_sum < 0.0 || strike <= 0.0 || discount <= 0.0) { + throw std::domain_error("payoff inputs are out of range"); + } +} + +double require_finite_result(double value) { + if (!std::isfinite(value)) { + throw std::overflow_error("payoff result is not finite"); + } + return value; +} + +} // namespace + +// 预计算每一步漂移、扩散项和全期限折现因子。 +GbmStepConstants make_gbm_step_constants(double maturity, double risk_free_rate, + double volatility, + std::uint32_t num_steps) { + require_finite(maturity, "maturity"); + require_finite(risk_free_rate, "risk-free rate"); + require_finite(volatility, "volatility"); + if (maturity <= 0.0 || volatility < 0.0 || num_steps == 0U) { + throw std::domain_error("GBM step parameters are out of range"); + } + + // 把总期限 T 切成 num_steps 段；CPU 和 GPU 都使用同一组常量， + // 因而两端模拟的是同一个风险中性 GBM 模型。 + const double dt = maturity / static_cast(num_steps); + const double drift = (risk_free_rate - 0.5 * volatility * volatility) * dt; + const double diffusion = volatility * std::sqrt(dt); + const double discount = std::exp(-risk_free_rate * maturity); + require_finite_result(dt); + require_finite_result(drift); + require_finite_result(diffusion); + require_finite_result(discount); + return {dt, drift, diffusion, discount}; +} + +double european_call_payoff(double terminal_spot, double strike, + double discount) { + validate_payoff_inputs(terminal_spot, strike, discount); + // 这里只计算一条路径的“今天价值”：到期的 max(S_T-K, 0) 先得到， + // 再乘 exp(-rT) 折回今天。之后 Monte Carlo 只需对很多条这样的值取平均。 + return require_finite_result(discount * + std::max(terminal_spot - strike, 0.0)); +} + +double asian_arithmetic_call_payoff(double running_sum, std::uint32_t num_steps, + double strike, double discount) { + validate_payoff_inputs(running_sum, strike, discount); + if (num_steps == 0U) { + throw std::domain_error("Asian arithmetic payoff requires steps"); + } + // running_sum 只累加演化后的 M 个监控点，刻意不把初始价格 S0 算入平均值。 + const double average = running_sum / static_cast(num_steps); + return require_finite_result(discount * std::max(average - strike, 0.0)); +} + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/statistics.cpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/statistics.cpp" new file mode 100644 index 00000000..84d7c92d --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/statistics.cpp" @@ -0,0 +1,115 @@ +// 从路径收益的原始矩计算均值、标准误和置信区间。 +#include "pricer/statistics.hpp" + +#include +#include +#include +#include +#include + +namespace pricer { +namespace { + +// 统计公式假定收益和平方收益均有效且非负。 +void validate_moments(const RawMoments &moments) { + if (!std::isfinite(moments.sum) || !std::isfinite(moments.sum_squares)) { + throw std::domain_error("moments must be finite"); + } + if (moments.sum < 0.0 || moments.sum_squares < 0.0) { + throw std::domain_error("moments must be non-negative"); + } +} + +double require_finite(double value, const char *name) { + if (!std::isfinite(value)) { + throw std::overflow_error(std::string(name) + " is not finite"); + } + return value; +} + +} // namespace + +// 分批计算时，原始矩可以直接相加而无需回放每条路径。 +RawMoments merge_raw_moments(const RawMoments &left, const RawMoments &right) { + // GPU 分批后，每一批只传回两个标量；此函数把它们合成一次完整实验的矩。 + validate_moments(left); + validate_moments(right); + if (right.count > std::numeric_limits::max() - left.count) { + throw std::overflow_error("moment count overflow"); + } + + const double sum = require_finite(left.sum + right.sum, "merged sum"); + const double sum_squares = require_finite( + left.sum_squares + right.sum_squares, "merged sum of squares"); + return {left.count + right.count, sum, sum_squares}; +} + +// 使用样本方差而非总体方差，并为单样本保留不可用统计量。 +PricingResult ResultAnalyzer::analyze(const RawMoments &moments, + std::optional reference_price) { + validate_moments(moments); + if (moments.count == 0U) { + throw std::domain_error("at least one sample is required"); + } + if (reference_price.has_value() && !std::isfinite(*reference_price)) { + throw std::domain_error("reference price must be finite"); + } + + const double count = static_cast(moments.count); + // 折现收益的样本均值就是蒙特卡洛价格估计。 + const double price = require_finite(moments.sum / count, "price"); + PricingResult result{price, std::nullopt, std::nullopt, + std::nullopt, std::nullopt, reference_price, + std::nullopt, std::nullopt, moments}; + + if (moments.count > 1U) { + // E[X^2] - E[X]^2 给出样本方差所需的分子；再由 stddev/sqrt(N) 得到 SE。 + const double expected_square = + require_finite(count * price * price, "variance term"); + double variance_numerator = moments.sum_squares - expected_square; + const double variance_scale = std::max( + {1.0, std::abs(moments.sum_squares), std::abs(expected_square)}); + const double roundoff_tolerance = + 64.0 * std::numeric_limits::epsilon() * variance_scale; + if (variance_numerator < 0.0) { + if (variance_numerator < -roundoff_tolerance) { + throw std::domain_error( + "sample variance is materially negative"); + } + variance_numerator = 0.0; + } + + const double variance = require_finite( + variance_numerator / static_cast(moments.count - 1U), + "sample variance"); + const double sample_stddev = + require_finite(std::sqrt(variance), "sample standard deviation"); + const double standard_error = + require_finite(sample_stddev / std::sqrt(count), "standard error"); + // 1.96 对应正态近似下的双侧 95% 置信区间。 + const double ci_lower = require_finite(price - 1.96 * standard_error, + "lower confidence interval"); + const double ci_upper = require_finite(price + 1.96 * standard_error, + "upper confidence interval"); + if (standard_error < 0.0 || ci_lower > ci_upper) { + throw std::domain_error("invalid confidence interval"); + } + result.sample_stddev = sample_stddev; + result.standard_error = standard_error; + result.ci_lower = ci_lower; + result.ci_upper = ci_upper; + } + + if (reference_price.has_value()) { + const double absolute_error = require_finite( + std::abs(price - *reference_price), "absolute error"); + result.absolute_error = absolute_error; + if (std::abs(*reference_price) > 1e-15) { + result.relative_error = require_finite( + absolute_error / std::abs(*reference_price), "relative error"); + } + } + return result; +} + +} // namespace pricer diff --git "a/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/version.cpp" "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/version.cpp" new file mode 100644 index 00000000..8e0486ca --- /dev/null +++ "b/05_option_pricing/\346\266\202\345\256\266\344\277\212/src/version.cpp" @@ -0,0 +1,8 @@ +// 返回构建产物和 CLI 共用的版本标识。 +#include "pricer/version.hpp" + +namespace pricer { + +const char *version() noexcept { return "0.1.0"; } + +} // namespace pricer