Skip to content

README.md

Latest commit

 

History

History
138 lines (102 loc) · 5.22 KB

README.md

File metadata and controls

138 lines (102 loc) · 5.22 KB

scipycpp

License: MIT C++17

Overview

scipycpp is a header-only C++ library implementing scipy's core APIs with bit-level precision alignment against Python scipy.

Built on proven C++ libraries:

  • numpycpp — numpy primitives
  • Eigen3 — linear algebra (solve, inv, det, eig, svd, cholesky)
  • pocketfft (bundled) — FFT (same library numpy/scipy use internally)
  • ckdtree (bundled) — KDTree from scipy.spatial

Quick Start

#include "scipy/core.h"

// Integration
auto [val, err] = scipy::integrate::quad([](double x){ return x*x; }, 0.0, 4.0);

// Stats
std::vector<double> pdf(3);
scipy::stats::norm_pdf(data, pdf.data(), 3);

// Linalg (Eigen3-backed)
double A[4] = {2,1,1,3}, b[2] = {5,6}, x[2];
scipy::linalg::solve(A, x, b, 2);

// FFT (pocketfft-backed — same as numpy/scipy)
scipy::fft::fft(signal, spectrum, 4);

// Spatial — KDTree with distances
scipy::spatial::KDTree<double> tree(points, n_pts, dim);
double dist; size_t idx;
tree.query(query_point, dist, idx);   // returns distance AND index

// Spatial — cross-set distance matrix (cdist)
scipy::spatial::distance::cdist(XA, mA, XB, mB, dim, dm);

// ndimage — 1D Gaussian filter
std::vector<double> out(n);
scipy::ndimage::gaussian_filter1d(src, out.data(), n, sigma);

// Signal — median filter
scipy::signal::medfilt(src, dst.data(), n, kernel_size);

// Transform — Rotation (from_matrix, from_euler, as_euler, as_matrix)
auto rot = scipy::spatial::transform::Rotation<double>::from_matrix(R);
auto euler = rot.as_euler_vec("xyz");  // returns [rx, ry, rz]

// from_euler + as_matrix — single axis (key use case: ego yaw)
double m9[9];
auto rot_z = scipy::spatial::transform::Rotation<double>::from_euler("z", yaw);
rot_z.as_matrix(m9);   // 3×3 row-major rotation matrix

// from_euler — multi-axis extrinsic / intrinsic
double angles[3] = {rx, ry, rz};
auto rot_xyz = scipy::spatial::transform::Rotation<double>::from_euler("xyz", angles);
rot_xyz.as_matrix(m9);

Dependencies

# Install dependencies
sudo dpkg -i numpycpp-dev-*.deb
sudo apt-get install libeigen3-dev

# Install scipycpp
mkdir build && cd build
cmake .. && make deb
sudo dpkg -i scipycpp-dev-*.deb
find_package(scipycpp REQUIRED)
target_link_libraries(myapp PRIVATE scipycpp::scipycpp)

Modules & ULP Alignment

All APIs are 0 ULP bit-identical to scipy for both float64/float32, including extreme values (±inf, NaN, ±0.0, subnormals, saturation inputs).

Module Backend Key APIs ULP Status
stats Cephes + numpycpp npy_exp norm.pdf, norm.cdf, norm.ppf ✅ 0 ULP
integrate sequential C++ sum trapezoid, simpson ⚠️ 0 ULP typical; ≤6 ULP uniform arrays
linalg Eigen3 partialPivLu solve ⚠️ ≤atol=1e-14
spatial pure C++ / scipy ckdtree cdist, KDTree ✅ 0 ULP
ndimage Python numpy kernel gaussian_filter1d ✅ 0 ULP
signal sort-based median medfilt ✅ 0 ULP
transform delegates to scipy Rotation.from_matrix, from_euler, as_euler, as_matrix ✅ 0 ULP

Full per-test ULP report: doc/ulp_report.csv (Auto-generated by pytest tests/test_all.py, see doc/ulp_report.md for summary)

Why some non-zero ULPs for linalg.solve?

API Max ULP Root Cause
linalg.solve ≤8.9e4 Eigen3 partialPivLu vs LAPACK gesv: different LU pivots produce different roundoff paths. Well-conditioned small matrices (2×2, identity) are bit-identical. All results within atol=1e-14 (ill-conditioned: atol=1e-10).

Integrate ULP alignment detail

API Input Max ULP Root Cause
trapezoid / simpson typical scientific data 0 ULP Sequential sum matches numpy pairwise for non-uniform data
trapezoid float64 uniform array ≤5 f64-ULP Sequential vs SIMD pairwise reorder (unavoidable without SIMD intrinsics)
trapezoid float32 uniform array ≤4 f32-ULP Same, measured in native float32 precision
simpson float64 uniform array ≤6 f64-ULP Same
simpson float32 input (any) ≤6 f32-ULP scipy computes internal sum in float32 SIMD; C++ uses sequential float32

scipy.integrate.simpson always returns float64 regardless of input dtype (Python float constants promote the result). scipy::integrate::simpson<T> matches this: it computes the intermediate sum in T (preserving scipy's precision path), then multiplies by 1.0/3.0 in double and returns double.

norm.pdf: numpycpp v1.21.2+ resolves exp via dlsym("npy_exp"), matching scipy's internal numpy math path bit-for-bit.

Testing

cd tests && make && make test
# → 299 tests (189 batch + 110 special/extreme values)
# ULP report printed to stderr + exported to doc/ulp_report.csv

See doc/ulp_report.md for the latest ULP alignment summary.

License

MIT