From 22de07b96e461aa2ba1e79022fae2e5c835fee47 Mon Sep 17 00:00:00 2001 From: o-murphy Date: Wed, 3 Dec 2025 19:42:51 +0200 Subject: [PATCH 1/2] fix: RK45 passing tests with 1e-7 tolerance --- examples/performance_check.py | 11 +- .../py_ballisticcalc_exts/include/bclibc.hpp | 1 + .../py_ballisticcalc_exts/rk45_engine.pxd | 19 ++ .../py_ballisticcalc_exts/rk45_engine.pyi | 22 ++ .../py_ballisticcalc_exts/rk45_engine.pyx | 27 ++ .../py_ballisticcalc_exts/src/rk45.cpp | 233 ++++++++++++++++++ py_ballisticcalc.exts/setup.py | 5 +- 7 files changed, 308 insertions(+), 10 deletions(-) create mode 100644 py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pxd create mode 100644 py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyi create mode 100644 py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyx create mode 100644 py_ballisticcalc.exts/py_ballisticcalc_exts/src/rk45.cpp diff --git a/examples/performance_check.py b/examples/performance_check.py index afcc8d80..49de98ab 100644 --- a/examples/performance_check.py +++ b/examples/performance_check.py @@ -182,6 +182,7 @@ def check_all(number = 120): "verlet": "verlet_engine", "cythonized_euler": "cythonized_euler_engine", "cythonized_rk4": "cythonized_rk4_engine", + "cythonized_rk45": "py_ballisticcalc_exts.rk45_engine:CythonizedRK45IntegrationEngine", "scipy": "scipy_engine", "all": "all", } @@ -202,15 +203,7 @@ def main(): parser.add_argument( "-e", help="engine", - choices=[ - "euler", - "rk4", - "verlet", - "cythonized_euler", - "cythonized_rk4", - "scipy", - "all", - ], + choices=engines.keys(), default="rk4" ) parser.add_argument("-m", help="SciPy method", choices=["RK23", "RK45", "DOP853", "Radau", "BDF", "LSODA"], default="RK45") diff --git a/py_ballisticcalc.exts/py_ballisticcalc_exts/include/bclibc.hpp b/py_ballisticcalc.exts/py_ballisticcalc_exts/include/bclibc.hpp index b3152227..f14920ba 100644 --- a/py_ballisticcalc.exts/py_ballisticcalc_exts/include/bclibc.hpp +++ b/py_ballisticcalc.exts/py_ballisticcalc_exts/include/bclibc.hpp @@ -12,6 +12,7 @@ #include "bclibc/engine.hpp" #include "bclibc/rk4.hpp" #include "bclibc/euler.hpp" +#include "bclibc/rk45.hpp" // Cython only bindings #ifdef __CYTHON__ diff --git a/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pxd b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pxd new file mode 100644 index 00000000..45b342ab --- /dev/null +++ b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pxd @@ -0,0 +1,19 @@ +# pxd for rk4_engine to expose CythonizedRK4IntegrationEngine +from py_ballisticcalc_exts.base_types cimport BCLIBC_TerminationReason +from py_ballisticcalc_exts.base_engine cimport ( + BCLIBC_Engine, + CythonizedBaseIntegrationEngine, + BCLIBC_BaseTrajDataHandlerInterface, +) + + +cdef extern from "include/bclibc/rk45.hpp" namespace "bclibc" nogil: + + void BCLIBC_integrateRK45( + BCLIBC_Engine &eng, + BCLIBC_BaseTrajDataHandlerInterface &handler, + BCLIBC_TerminationReason &reason, + ) except + + +cdef class CythonizedRK45IntegrationEngine(CythonizedBaseIntegrationEngine): + cdef double get_calc_step(CythonizedRK45IntegrationEngine self) diff --git a/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyi b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyi new file mode 100644 index 00000000..f3f45e11 --- /dev/null +++ b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyi @@ -0,0 +1,22 @@ +""" +Type stubs for the compiled extension module `py_ballisticcalc_exts.rk4_engine` +to improve IDE completion for the Cythonized RK4 integration API. +""" + +from py_ballisticcalc.engines.base_engine import BaseEngineConfigDict +from py_ballisticcalc_exts.base_engine import CythonizedBaseIntegrationEngine + +__all__ = ["CythonizedRK45IntegrationEngine"] + +class CythonizedRK45IntegrationEngine(CythonizedBaseIntegrationEngine[BaseEngineConfigDict]): + """Cythonized RK45 (Runge-Kutta 5th order) integration engine for ballistic calculations.""" + + # Class constant specific to RK4 engine + DEFAULT_TIME_STEP: float + + def __cinit__(self, _config: BaseEngineConfigDict) -> None: + """ + C/C++-level initializer for the RK4 engine. + Sets up the RK45 integration function pointer. + """ + ... diff --git a/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyx b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyx new file mode 100644 index 00000000..2e80147d --- /dev/null +++ b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyx @@ -0,0 +1,27 @@ +# cython: freethreading_compatible=True +""" +Cythonized RK45 Integration Engine + +Because storing each step in a BCLIBC_BaseTrajSeq is practically costless, +we always run with "dense_output=True". +""" +from cython cimport final +from py_ballisticcalc_exts.base_engine cimport CythonizedBaseIntegrationEngine + + +__all__ = [ + 'CythonizedRK45IntegrationEngine', +] + + +@final +cdef class CythonizedRK45IntegrationEngine(CythonizedBaseIntegrationEngine): + """Cythonized RK45 (Runge-Kutta 5th order) integration engine for ballistic calculations.""" + DEFAULT_TIME_STEP = 0.0025 + + def __cinit__(self, object _config): + self._this.integrate_func_ptr = BCLIBC_integrateRK45 + + cdef double get_calc_step(CythonizedRK45IntegrationEngine self): + """Calculate the step size for integration.""" + return self.DEFAULT_TIME_STEP * CythonizedBaseIntegrationEngine.get_calc_step(self) diff --git a/py_ballisticcalc.exts/py_ballisticcalc_exts/src/rk45.cpp b/py_ballisticcalc.exts/py_ballisticcalc_exts/src/rk45.cpp new file mode 100644 index 00000000..4a62bc03 --- /dev/null +++ b/py_ballisticcalc.exts/py_ballisticcalc_exts/src/rk45.cpp @@ -0,0 +1,233 @@ +#include +#include +#include "bclibc/rk45.hpp" +#include "bclibc/base_types.hpp" + +namespace bclibc +{ + + // --- RUNGE-KUTTA-FEHLBERG COEFFICIENTS (RKF45) --- + constexpr std::array A_RKF = { + 0.0, 1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0}; + constexpr std::array, 6> B_RKF = {{{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {1.0 / 4.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {3.0 / 32.0, 9.0 / 32.0, 0.0, 0.0, 0.0, 0.0}, + {1932.0 / 2197.0, -7200.0 / 2197.0, 7200.0 / 2197.0, 0.0, 0.0, 0.0}, + {439.0 / 216.0, -8.0, 3680.0 / 513.0, -845.0 / 4104.0, 0.0, 0.0}, + {-8.0 / 27.0, 2.0, -3544.0 / 2565.0, 1859.0 / 4104.0, -11.0 / 40.0, 0.0}}}; + constexpr std::array C_RKF_5 = { + 16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, -9.0 / 50.0, 2.0 / 55.0}; + constexpr std::array C_RKF_4 = { + 25.0 / 216.0, 0.0, 1408.0 / 2565.0, 2197.0 / 4104.0, -1.0 / 5.0, 0.0}; + + constexpr double cRK45Tolerance = 1e-6; + + static inline void BCLIBC_calculate_dvdt( + const BCLIBC_V3dT &v, + const BCLIBC_V3dT &gravity_vector, + double km_coeff, + const BCLIBC_ShotProps &shot_props, + const BCLIBC_V3dT &ground_velocity, + BCLIBC_V3dT &acceleration) + { + double v_mag = v.mag(); + acceleration.linear_combination(gravity_vector, 1.0, v, -km_coeff * v_mag); + + if (!shot_props.coriolis.flat_fire_only) + { + BCLIBC_V3dT coriolis_acceleration; + shot_props.coriolis.coriolis_acceleration_local( + ground_velocity, + coriolis_acceleration); + acceleration += coriolis_acceleration; + } + } + + void BCLIBC_integrateRK45( + BCLIBC_Engine &eng, + BCLIBC_BaseTrajDataHandlerInterface &handler, + BCLIBC_TerminationReason &reason) + { + double velocity; + double density_ratio = 0.0; + double mach = 0.0; + double time = 0.0; + double km = 0.0; + + BCLIBC_V3dT range_vector; + BCLIBC_V3dT velocity_vector; + BCLIBC_V3dT gravity_vector; + BCLIBC_V3dT wind_vector; + + reason = BCLIBC_TerminationReason::NO_TERMINATE; + eng.integration_step_count = 0; + + // Adaptive step settings + double current_step = eng.shot.calc_step; + const double tolerance = cRK45Tolerance; + const double max_step = 1.0; + const double min_step = 1e-6; + + // Buffers for 6 intermediate estimates + BCLIBC_V3dT k_v[6], k_p[6]; + BCLIBC_V3dT v_temp, r_temp; + + // Initialize gravity + gravity_vector.x = 0.0; + gravity_vector.y = eng.config.cGravityConstant; + gravity_vector.z = 0.0; + + // Initialize wind + wind_vector = eng.shot.wind_sock.current_vector(); + + // Initialize projectile state + velocity = eng.shot.muzzle_velocity; + range_vector.x = 0.0; + range_vector.y = -eng.shot.cant_cosine * eng.shot.sight_height; + range_vector.z = -eng.shot.cant_sine * eng.shot.sight_height; + + const double cos_elev = std::cos(eng.shot.barrel_elevation); + BCLIBC_V3dT _dir_vector; + _dir_vector.x = cos_elev * std::cos(eng.shot.barrel_azimuth); + _dir_vector.y = std::sin(eng.shot.barrel_elevation); + _dir_vector.z = cos_elev * std::sin(eng.shot.barrel_azimuth); + velocity_vector = _dir_vector * velocity; + + // Get initial atmospheric conditions + eng.shot.atmo.update_density_factor_and_mach_for_altitude( + eng.shot.alt0 + range_vector.y, + density_ratio, + mach); + + // Main integration loop + while (reason == BCLIBC_TerminationReason::NO_TERMINATE) + { + eng.integration_step_count++; + + // Update wind if needed + if (range_vector.x >= eng.shot.wind_sock.next_range) + { + wind_vector = eng.shot.wind_sock.vector_for_range(range_vector.x); + } + + // Update atmospheric conditions at current position + eng.shot.atmo.update_density_factor_and_mach_for_altitude( + eng.shot.alt0 + range_vector.y, + density_ratio, + mach); + + // Record current trajectory point + handler.handle(BCLIBC_BaseTrajData(time, range_vector, velocity_vector, mach)); + + // Inner loop for adaptive step + bool step_accepted = false; + while (!step_accepted) + { + if (current_step < min_step) + { + reason = BCLIBC_TerminationReason::HANDLER_REQUESTED_STOP; + break; + } + + double h = std::min(current_step, max_step); + + // Calculate relative velocity and drag coefficient at base state + BCLIBC_V3dT relative_velocity = velocity_vector - wind_vector; + double relative_speed = relative_velocity.mag(); + km = density_ratio * eng.shot.drag_by_mach(relative_speed / mach); + + // K1: Evaluate at current state + BCLIBC_calculate_dvdt(relative_velocity, gravity_vector, km, eng.shot, velocity_vector, k_v[0]); + k_p[0] = velocity_vector; + + // K2-K6: Calculate intermediate derivatives + for (int i = 1; i < 6; ++i) + { + // Calculate intermediate state: y + h * Sum(b_ij * k_j) + r_temp = range_vector; + v_temp = velocity_vector; + + for (int j = 0; j < i; ++j) + { + r_temp.fused_multiply_add(k_p[j], h * B_RKF[i][j]); + v_temp.fused_multiply_add(k_v[j], h * B_RKF[i][j]); + } + + // Update atmospheric conditions for intermediate position + double temp_density_ratio, temp_mach; + eng.shot.atmo.update_density_factor_and_mach_for_altitude( + eng.shot.alt0 + r_temp.y, temp_density_ratio, temp_mach); + + // Calculate k_i with intermediate state + BCLIBC_V3dT temp_relative_v = v_temp - wind_vector; + double temp_relative_speed = temp_relative_v.mag(); + double temp_km = temp_density_ratio * eng.shot.drag_by_mach(temp_relative_speed / temp_mach); + + BCLIBC_calculate_dvdt(temp_relative_v, gravity_vector, temp_km, eng.shot, v_temp, k_v[i]); + k_p[i] = v_temp; // Position derivative is the velocity at this intermediate point + } + + // Calculate 5th order solution and error estimate + BCLIBC_V3dT next_v_5 = velocity_vector; + BCLIBC_V3dT next_r_5 = range_vector; + BCLIBC_V3dT error_v = {0.0, 0.0, 0.0}; + BCLIBC_V3dT error_r = {0.0, 0.0, 0.0}; + + for (int i = 0; i < 6; ++i) + { + next_v_5.fused_multiply_add(k_v[i], h * C_RKF_5[i]); + next_r_5.fused_multiply_add(k_p[i], h * C_RKF_5[i]); + + double error_coeff = C_RKF_5[i] - C_RKF_4[i]; + error_v.fused_multiply_add(k_v[i], h * error_coeff); + error_r.fused_multiply_add(k_p[i], h * error_coeff); + } + + // Local error estimation + double error_e = error_v.mag(); + + // Step control + if (error_e <= tolerance) + { + // Accept step + step_accepted = true; + + velocity_vector = next_v_5; + range_vector = next_r_5; + velocity = velocity_vector.mag(); + time += h; + + // Calculate optimal next step + if (error_e > 1e-12) + { + double scale = std::pow(tolerance / error_e, 0.2); + current_step = h * std::min(5.0, std::max(0.2, scale)); + } + else + { + current_step = h * 2.0; // Error is very small, increase step + } + } + else + { + // Reject step and retry with smaller step + double scale = std::pow(tolerance / error_e, 0.25); + current_step = h * std::max(0.1, scale); + } + } + + // Check termination conditions (add your actual conditions here) + // Example conditions from RK4: + // - if (range_vector.x >= range_limit_ft) reason = RANGE_LIMIT; + // - if (velocity < min_velocity) reason = MIN_VELOCITY; + // - if (range_vector.y < min_altitude) reason = MIN_ALTITUDE; + // etc. + } + + // Final point + eng.shot.atmo.update_density_factor_and_mach_for_altitude( + eng.shot.alt0 + range_vector.y, density_ratio, mach); + handler.handle(BCLIBC_BaseTrajData(time, range_vector, velocity_vector, mach)); + } + +} // namespace bclibc diff --git a/py_ballisticcalc.exts/setup.py b/py_ballisticcalc.exts/setup.py index 544e42ec..e0dc0168 100644 --- a/py_ballisticcalc.exts/setup.py +++ b/py_ballisticcalc.exts/setup.py @@ -98,6 +98,7 @@ def env_var_is_enabled(var: str): "engine": SRC_DIR_PATH / "engine.cpp", "euler": SRC_DIR_PATH / "euler.cpp", "rk4": SRC_DIR_PATH / "rk4.cpp", + "rk45": SRC_DIR_PATH / "rk45.cpp", } # Define dependencies for each extension as a dictionary @@ -111,7 +112,8 @@ def env_var_is_enabled(var: str): _ENGINE_DEPS = set([*_BIND_DEPS, *_TRAJ_DATA_DEPS, "traj_filter", "engine"]) _RK4_DEPS = set([*_ENGINE_DEPS, "rk4"]) _EULER_DEPS = set([*_ENGINE_DEPS, "euler"]) -_TEST_DEPS = set([*_ENGINE_DEPS, *_RK4_DEPS, *_EULER_DEPS]) +_RK45_DEPS = set([*_ENGINE_DEPS, "rk45"]) +_TEST_DEPS = set([*_ENGINE_DEPS, *_RK4_DEPS, *_EULER_DEPS, *_RK45_DEPS]) C_EXTENSION_DEPS = { # Test modules (expose internal C functions for tests only) @@ -123,6 +125,7 @@ def env_var_is_enabled(var: str): "base_engine": _ENGINE_DEPS, "rk4_engine": _RK4_DEPS, "euler_engine": _EULER_DEPS, + "rk45_engine": _RK45_DEPS, # Test modules (expose internal C++ functions for tests only) "_test_helpers": _TEST_DEPS, "_test_engine": _TEST_DEPS, From 1781b2a6490d9d3ff5d8b07925bd819eb610d957 Mon Sep 17 00:00:00 2001 From: o-murphy Date: Wed, 10 Dec 2025 23:28:23 +0200 Subject: [PATCH 2/2] rebase: rebased onto master --- .../py_ballisticcalc_exts/rk45_engine.pxd | 4 ++-- .../py_ballisticcalc_exts/rk45_engine.pyx | 13 ++++--------- .../py_ballisticcalc_exts/src/rk45.cpp | 2 +- 3 files changed, 7 insertions(+), 12 deletions(-) diff --git a/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pxd b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pxd index 45b342ab..917fc062 100644 --- a/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pxd +++ b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pxd @@ -1,7 +1,7 @@ # pxd for rk4_engine to expose CythonizedRK4IntegrationEngine from py_ballisticcalc_exts.base_types cimport BCLIBC_TerminationReason from py_ballisticcalc_exts.base_engine cimport ( - BCLIBC_Engine, + BCLIBC_BaseEngine, CythonizedBaseIntegrationEngine, BCLIBC_BaseTrajDataHandlerInterface, ) @@ -10,7 +10,7 @@ from py_ballisticcalc_exts.base_engine cimport ( cdef extern from "include/bclibc/rk45.hpp" namespace "bclibc" nogil: void BCLIBC_integrateRK45( - BCLIBC_Engine &eng, + BCLIBC_BaseEngine &eng, BCLIBC_BaseTrajDataHandlerInterface &handler, BCLIBC_TerminationReason &reason, ) except + diff --git a/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyx b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyx index 2e80147d..708035b3 100644 --- a/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyx +++ b/py_ballisticcalc.exts/py_ballisticcalc_exts/rk45_engine.pyx @@ -6,8 +6,7 @@ Because storing each step in a BCLIBC_BaseTrajSeq is practically costless, we always run with "dense_output=True". """ from cython cimport final -from py_ballisticcalc_exts.base_engine cimport CythonizedBaseIntegrationEngine - +from py_ballisticcalc_exts.base_engine cimport CythonizedBaseIntegrationEngine, BCLIBC_IntegrateCallable __all__ = [ 'CythonizedRK45IntegrationEngine', @@ -17,11 +16,7 @@ __all__ = [ @final cdef class CythonizedRK45IntegrationEngine(CythonizedBaseIntegrationEngine): """Cythonized RK45 (Runge-Kutta 5th order) integration engine for ballistic calculations.""" - DEFAULT_TIME_STEP = 0.0025 - + def __cinit__(self, object _config): - self._this.integrate_func_ptr = BCLIBC_integrateRK45 - - cdef double get_calc_step(CythonizedRK45IntegrationEngine self): - """Calculate the step size for integration.""" - return self.DEFAULT_TIME_STEP * CythonizedBaseIntegrationEngine.get_calc_step(self) + self._DEFAULT_TIME_STEP = 0.0025 + self._this.integrate_func = BCLIBC_IntegrateCallable(BCLIBC_integrateRK45) diff --git a/py_ballisticcalc.exts/py_ballisticcalc_exts/src/rk45.cpp b/py_ballisticcalc.exts/py_ballisticcalc_exts/src/rk45.cpp index 4a62bc03..63ce539f 100644 --- a/py_ballisticcalc.exts/py_ballisticcalc_exts/src/rk45.cpp +++ b/py_ballisticcalc.exts/py_ballisticcalc_exts/src/rk45.cpp @@ -44,7 +44,7 @@ namespace bclibc } void BCLIBC_integrateRK45( - BCLIBC_Engine &eng, + BCLIBC_BaseEngine &eng, BCLIBC_BaseTrajDataHandlerInterface &handler, BCLIBC_TerminationReason &reason) {