Find the way to provide more efficient methods of calculation

[o-murphy]

[dbookstaber]

This is implemented in pull request #50

[dbookstaber]

I implemented constant time-step calculation here: https://github.com/dbookstaber/py_ballistics/blob/time-step/py_ballisticcalc/trajectory_calc.py

However it looks like it will always take more iterations than the current method to achieve any particular level of precision, so I do not recommend pursuing this.

[o-murphy]

but in addition to the number of iterations, it is necessary to take into account the time for performing the calculation, it may turn out that due to a decrease in the number of calculations, the number of iterations can be increased without losing performance

we have to make a profiling and performance comparing
u can use timeit built-in module to got time some function are executing

[dbookstaber]

By "iterations" I meant the number of passes through the #region Trajectory Loop to calculate any particular trajectory. (This is without regard to the zero-finding logic.)

I think the current method of setting the time step based on the velocity is optimal. Any other approach will create excess precision at some points and lower precision at others.

[dbookstaber]

To improve efficiency, JBM recommends switching from this simple integration approach to something like RK4, which is what he uses.

[dbookstaber]

I implemented RK4 integration in this branch. It is more robust than the Euler implementation we are running, meaning that as step_size increases errors grow more slowly. However on the example I studied it's not like orders of magnitude better.

This provides a foundation for some desirable improvements, including:

1. Adaptive step_size
2. I included two different RK4 implementations, one using time and the other using range.x as the independent variable. The latter (CalcMethod.RK_dx) version could particularly benefit TrajectoryCalc.zero_angle since, as long as zero_distance % step_size == 0 then it will end exactly at the distance of interest.

[o-murphy]

Right now I don't have much time to study alternative calculation methods, so I'm completely relying on you. In any case, if we're worried about breaking something, we can always implement it as an experimental feature and enable it via the library's global settings. Like we do with gunpowder sensitivity, we can enable an alternative calculation algorithm. Use inheritance, it allows to implement new feature without massive base code changes and lost of backwards compatibility and also allows extend functionality without changing major release version

Create new module, inherit from existing backend.TrajectoryCalculator

# my_new_solver.py
from py_ballisticcalc.backend import TrajectoryCalculator

class MyNewSolver(TrajectoryCalculator):

    def trajectory(*args, **kwargs):
        # reimplement methods you need
        ...

Add global flags to trajectory_calc.py

# existing flags
_globalUsePowderSensitivity = False
def set_global_use_powder_sensitivity(value: bool) -> None:
    # pylint: disable=global-statement
    global _globalUsePowderSensitivity
    if not isinstance(value, bool):
        raise TypeError(f"set_global_use_powder_sensitivity {value=} is not a boolean")
    _globalUsePowderSensitivity = value
    
# new flag
_globalUseNewSolver = False
def set_global_use_new_solver(value: bool) -> None:
    global _globalUseNewSolver
    _globalUseNewSolver = value

Add conditional import of your new solver to interface.Calculator

from py_ballisticcalc.backend import TrajectoryCalc, get_global_use_new_solver

# Calculator class
def barrel_elevation_for_target(self, shot: Shot, target_distance: Union[float, Distance]) -> Angular:
    if get_global_use_new_solver():
        from my_new_solver import MyNewSolver
        calc = MyNewSolver()
    else:
       calc = TrajectoryCalc()

[serhiy-yevtushenko]

Hi, David (@dbookstaber). I have a question for you as expert in ballistics - could you comment on possibility of using methods like described in these and similar articles - https://iopscience.iop.org/article/10.1088/0143-0807/37/3/035001, or http://www.jaac-online.com/article/doi/10.11948/20210277 for speeding up the computations? (Search in the internet gives as well a number of other articles). Are they applicable for precise simulations used in ballistics calculator case or only for gaining insights in general behaviour of the system?

[o-murphy]

@serhiy-yevtushenko

In theory, we are considering the possibility of adding alternative calculation methods, currently the plan is to implement RK. This is something that can be implemented without a significant change in the library architecture without violating backward API compatibility.

As project owner and maintainer I think that after reworking the library architecture, there probably will be a clear opportunity to expand the functionality with methods such as the Adomian decomposition method.
However, it is worth understanding that some calculation methods will significantly reduce performance, with a slight improvement in accuracy. They may be interesting from a scientific point of view, but not so necessary for the ordinary user, whom the project is aimed at.

In addition, the library is developed with minimalism and lightness in mind, while maintaining cross-platform support and the smallest number of mandatory dependencies. Therefore, such innovations will be better implemented as optionally installed plugins.

The next stable version is planned to be v2.1.0. If you change the architecture, please increase the API version to v3

[serhiy-yevtushenko]

Hi, Dmytro. Thank you for the clarification. I believe efficiency could be considered from several viewpoints:

• precision of the computation
• speed of the computation

In my opinion, the benefit of analytical methods lays in speed of computation (compare calculations for drag-free case using formulas from Wikipedia with numerical integration like Euler or Runge-Kutta method).

Runge-Kutta method is optimizing the precision of the computations.

As I am not very experienced in ballistics, it is really interesting for me to hear David's opinion about the usage of analytical methods and their limitations. I'm able to see a lot of articles published in this area, but do not know in general their limitations and areas of applicability. Just for understanding - in some articles on analytical methods it is mentioned that computations could be easily done using a simple calculator (not a ballistic one). So, they are offering at least from a speed perspective a good potential for improvement. At the same time, I have not seen yet articles, that discuss using analytical methods for the prediction trajectories of specific projectiles with specified ballistic coefficient and G-function. Hopefully, this clarifies my question.

[dbookstaber]

@serhiy-yevtushenko, The sorts of analytic approximations discussed in those papers are at best of academic interest, and not considered useful for applied high-velocity ballistics where the drag coefficient as a function of velocity is not only variable but also must be empirically measured. Never mind accounting for variable atmospheric conditions over the trajectory, as we do here!

Applied ballistics passed the point of any speed/precision trade-off half a century ago. Runge-Kutta calculation of trajectories was "fast enough" on early 1980s-era electronics to produce solutions that exceed the precision of the inputs in fractions of a second.

(By comparison, the Euler implementation here is computationally wasteful, but on modern electronics we set step sizes that are orders of magnitude smaller than necessary for precise solutions and we don't bother fine-tuning for speed – and often don't even bother compiling the code!)

TL;DR: Basic numeric integration on these problems is so fast and precise that it's not worth making other approximations or looking for analytic shortcuts.

[serhiy-yevtushenko]

Hi, David. Thank you for the clarification. My experience of running simulation of shots with Euler using py-ballisticcalc was about 40 minutes for 900 shots - so, I would not call it blasing fast (using pure python version). I wonder what kind of performance are you getting with Runge-Kutta method

[dbookstaber]

Yeah, pure python is more of a luxury for development. When you care about speed you should:

1. Install the Cython (compiled) version. (IIRC this already increases speed over 100x.)
2. Increase _globalMaxCalcStepSizeFeet (to suit your precision requirements – 10x the default is an easy starting point if you want further speed).

We haven't completed Runge-Kutta implementation here, but it should offer a further 1-2 orders of magnitude improvement in speed over Euler.

[serhiy-yevtushenko]

Hi, David. Thank you for the tips. Could you clarify, what is the performance impact of extra-data=True? I tried now running the simulation with extensions in version 2.0.8b1, and the performance was even worse than observed before without extensions (and at the end, I got exceptions, and computation did not complete due to new exceptions)

[dbookstaber]

@serhiy-yevtushenko: We haven't ported the vertical fixes into Cython yet; I think @o-murphy will finish that shortly.

If you're not seeing orders of magnitude improvement with ext version then you must not be using the compiled code. But if you have warnings enabled and you are not seeing warnings.warn("Library running in pure python mode. " then it must at least think it is using the binaries and something else has gone wrong.

Performance impact with extra_data=True should be small but not negligible: the overhead is that instead of just recording TrajectoryData at the step frequency you request, it also does it for every horizontal 0.2ft (this is setting _globalChartResolution in next version) or (if indicated) time_step seconds of travel.

[o-murphy]

@serhiy-yevtushenko
In fact, numpy and scipy will be slower than cython in this case, one of the reasons is numpy vectors are not optimized for working with small vectors like 3d vectors and have a lot overheads on initialization. There best way could be C vectors probably

[o-murphy]

@serhiy-yevtushenko
I tried to switch from python list to c-array in for CurvePoints tooday and i got about 5-10% better performance with extra_data=True, but also i got unexpected return values instead of expected RangeError in some tests which happen not always and can't solve it yet.
so yes, the further you dig into cytonization the greater the chance of getting unexpected errors, memory leaks, and less readable and maintainable code

Now I'm on the way to partially refactor library. I wanna write a generics for calculation "Method"s, that can be used with the same interface and input data and inherit from it. the first methods will be available after refactoring is current existed (Euler) and RK4, also the cythonic code will implemented as optional installable "Method". I mean I wanna create a core where you can chose which method u want to use, or even implement your own without changing the entire lib sources

[o-murphy]

After such a refactoring, it will not matter how low-level the method is implemented, it can be completely written in another language, or even executable on a remote server, all that is needed is a wrapper that will implement the interface.

[o-murphy]

All this would allow me to focus on implementing the functionality, and let’em cook, who’s better versed math and ballistics

[serhiy-yevtushenko]

@serhiy-yevtushenko In fact, numpy and scipy will be slower than cython in this case, one of the reasons is numpy vectors are not optimized for working with small vectors like 3d vectors and have a lot overheads on initialization. There best way could be C vectors probably

I believe that usage of numpy requires reorganization from viewpoint of thinking about the code. In order to use numpy optimally, one would not try to allocate a lot of small vectors, but would rather do numpy array with dimension (N, 3), where N is the preallocated number of steps to be used. One could look at scikit-learn (from viewpoint on how they organize interfaces of machine learning methods) on good patterns to use, as scikit-learn relies both heavily on numpy and cython.

[serhiy-yevtushenko]

@serhiy-yevtushenko I tried to switch from python list to c-array in for CurvePoints tooday and i got about 5-10% better performance with extra_data=True, but also i got unexpected return values instead of expected RangeError in some tests which happen not always and can't solve it yet. so yes, the further you dig into cytonization the greater the chance of getting unexpected errors, memory leaks, and less readable and maintainable code

Now I'm on the way to partially refactor library. I wanna write a generics for calculation "Method"s, that can be used with the same interface and input data and inherit from it. the first methods will be available after refactoring is current existed (Euler) and RK4, also the cythonic code will implemented as optional installable "Method". I mean I wanna create a core where you can chose which method u want to use, or even implement your own without changing the entire lib sources

Sounds like a great plan

[serhiy-yevtushenko]

Here is the promised example for the reproduction of the cython/python differences.

This is a simplified version of script which I was running. I used it here for the computation of the maximal reachability zone of .50 BMG bullet, and reduced amount of computed angles almost 10 times (from 901 to 91)
This example requires scipy for the convex hull computation. If you don't want to add scipy, you could just remove the part of computing convex hull

import math
import time
from typing import Final, NamedTuple
from collections.abc import Callable

import numpy as np
from py_ballisticcalc import (
    DragModel,
    TableG1,
    Ammo,
    Velocity,
    Temperature,
    Weapon,
    Distance,
    Shot,
    Calculator,
    HitResult,
    Angular,
    InterfaceConfigDict,
    RangeError,
    TrajFlag,
    Weight,
)
from scipy.spatial import ConvexHull

EARTH_GRAVITY_CONSTANT: Final[float] = 9.81


def generate_reachable_region(
    calculator: Callable[[float], np.ndarray[tuple[float, float]]],
    angles: np.ndarray[float] | list[float | int],
    show_progress: bool = False,
) -> np.ndarray:
    """
    Generate all reachable points for a range of angles.
    Calculator should perform computation of trajectory till zero level, and
    return array of numpy points with trajectory points
    """
    if show_progress:
        import tqdm

        angles = tqdm.tqdm(angles, desc="Computing data for angles")
    points = []
    for angle in angles:
        trajectory = calculator(angle)
        points.extend(trajectory)
    return np.array(points)


def generate_reachable_region_points_for_angle_range(
    calculator: Callable[[float], np.ndarray[tuple[float, float]]],
    start_angle_in_degrees: float,
    end_angle_in_degrees: float,
    num_points: int,
    show_progress: bool = False,
):
    return generate_reachable_region(
        calculator,
        np.linspace(start_angle_in_degrees, end_angle_in_degrees, num_points),
        show_progress=show_progress,
    )


def extract_convex_hull_segments(
    hull: ConvexHull, points: np.ndarray
) -> list[tuple[tuple[float, float], tuple[float, float]]]:
    """
    Extract the line segments (edges) of a convex hull.

    Parameters:
        hull: ConvexHull object.
        points: Array of all points used to compute the hull.

    Returns:
        List of segments [(x1, y1, x2, y2), ...].
    """
    segments = []
    hull_points = points[hull.vertices]  # Extract vertices forming the convex hull
    for i in range(len(hull_points)):
        x1, y1 = hull_points[i]
        x2, y2 = hull_points[
            (i + 1) % len(hull_points)
        ]  # Wrap around to form the last segment
        segments.append(((x1, y1), (x2, y2)))
    return segments


def compute_convex_hull_segments(
    calculator: Callable[[float], np.ndarray[tuple[float, float]]],
    start_angle_in_degrees: float,
    end_angle_in_degrees: float,
    num_points: int,
    show_progress: bool = False,
):
    start_time = time.time()
    points = generate_reachable_region_points_for_angle_range(
        calculator,
        start_angle_in_degrees,
        end_angle_in_degrees,
        num_points,
        show_progress=show_progress,
    )
    # Compute convex hull segments
    end_time = time.time()
    print(f"There were {len(points)=:_} point computed")
    print(
        f"Points generation for convex hull {num_points=} has taken {end_time - start_time:.3f} seconds"
    )
    return extract_convex_hull_segments(ConvexHull(points), points)


def compute_weapon_reachability_convex_hull(
    calculator,
    num_points: int,
    start_angle_in_degrees: float = 0,
    end_angle_in_degrees: float = 90,
    show_progress: bool = False,
):
    segments = compute_convex_hull_segments(
        calculator,
        start_angle_in_degrees,
        end_angle_in_degrees,
        num_points,
        show_progress,
    )
    print(f"{len(segments)=:_}")
    return segments


class WeaponData(NamedTuple):
    weapon: Weapon
    ammo: Ammo

    @property
    def muzzle_velocity(self):
        return ammo.mv

    def compute_trajectory(
        self, angle_in_degrees: float, distance_in_meters: float, extra_data=True
    ):
        gun = Weapon()

        shot = Shot(
            weapon=gun,
            ammo=ammo,
            relative_angle=Angular.Degree(angle_in_degrees),
        )
        config = InterfaceConfigDict(
            cMinimumVelocity=0, cMinimumAltitude=-1, cMaximumDrop=-1
        )

        calc = Calculator(_config=config)
        try:
            shot_result = calc.fire(
                shot, Distance.Meter(distance_in_meters), extra_data=extra_data
            )
        except RangeError as e:
            if e.reason in [
                RangeError.MaximumDropReached,
                RangeError.MinimumAltitudeReached,
            ]:
                shot_result = HitResult(shot, e.incomplete_trajectory, extra=extra_data)
            else:
                raise e
        return shot_result


def generic_weapon_wrapper(
    weapon: WeaponData,
    interesting_point_finder: Callable[[HitResult], int],
    max_range_calculator: Callable[[float, float], float],
    distance_unit:Distance=Distance.Meter,
):
    def calculate_trajectory_points(
        angle_in_degrees: float,
    ) -> list[tuple[float, float]]:
        max_possible_range = max_range_calculator(
            weapon.muzzle_velocity >> Velocity.MPS, angle_in_degrees
        )
        shot = weapon.compute_trajectory(angle_in_degrees, max_possible_range)
        last_entry_index = interesting_point_finder(shot)
        if last_entry_index != -1:
            trajectory_points = shot.trajectory[:last_entry_index]
        else:
            trajectory_points = shot.trajectory
        result = []
        for p in trajectory_points:
            result.append((p.distance >> distance_unit, p.height >> distance_unit))
        return result

    return calculate_trajectory_points


def find_touch_point_index(shot: HitResult) -> int:
    for i, p in enumerate(shot.trajectory):
        if p.flag & TrajFlag.ZERO_DOWN:
            return i
    return -1


def calculate_drag_free_range(
    velocity_in_mps: float, angle_in_degrees: float, gravity: float = EARTH_GRAVITY_CONSTANT
):
    """Compute maximal range for projectile thrown with velocity in the drag-free environment"""
    angle_rad = math.radians(angle_in_degrees)
    return (velocity_in_mps ** 2 * math.sin(2 * angle_rad)) / gravity


def weapon_wrapper_for_max_distance(weapon: WeaponData, distance_unit:Distance=Distance.Meter):
    return generic_weapon_wrapper(
        weapon, find_touch_point_index, calculate_drag_free_range, distance_unit
    )


if __name__ == "__main__":
    dm = DragModel(
        0.62, TableG1, Weight.Grain(661), Distance.Inch(0.51), Distance.Inch(2.3)
    )
    ammo = Ammo(dm, Velocity.MPS(850), Temperature.Celsius(15))
    weapon = Weapon(sight_height=Distance.Centimeter(9), twist=Distance.Inch(15))

    compute_weapon_reachability_convex_hull(
        weapon_wrapper_for_max_distance(WeaponData(weapon, ammo)),
        91,
        0,
        90,
        show_progress=True,
    )

I have added option to switch usage of tqdm for showing progress on and off. I have not observed a difference in performance on my laptop.

My personal preference is to use tqdm for long-running task, as in such a way I can see, whether task is doing progress at all

Here are the run times on my laptop:
No Cython, tqdm:

UserWarning: Library running in pure python mode. For better performance install 'py_ballisticcalc.exts' binary package
  warnings.warn("Library running in pure python mode. "
Computing data for angles: 100%|██████████| 91/91 [03:28<00:00,  2.30s/it]
There were len(points)=5_882_538 point computed
Points generation for convex hull num_points=91 has taken 210.843 seconds
len(segments)=26_829

No Cython, no tqdm:

UserWarning: Library running in pure python mode. For better performance install 'py_ballisticcalc.exts' binary package
  warnings.warn("Library running in pure python mode. "
There were len(points)=5_882_538 point computed
Points generation for convex hull num_points=91 has taken 213.612 seconds
len(segments)=26_829

Cython, tqdm:

Computing data for angles: 100%|██████████| 91/91 [02:02<00:00,  1.35s/it]
There were len(points)=5_883_091 point computed
Points generation for convex hull num_points=91 has taken 124.425 seconds
len(segments)=26_837

Cython, no tqdm:

There were len(points)=5_883_091 point computed
Points generation for convex hull num_points=91 has taken 130.431 seconds
len(segments)=26_837

[serhiy-yevtushenko]

All this would allow me to focus on implementing the functionality, and let’em cook, who’s better versed math and ballistics

Short question - I would like to do an experiment with implementing a TrajectoryCalculator with ivp_solve from scipy. Which unit tests are sensible to use as a starting point and is there a recommended sequence of unit test from functionality progression viewpoint (simplest case -> more complex parameters like handling of air pressure, winds, ...)?

[dbookstaber]

Judging from my experience of using cython and writing extensions to python code on rust, the performance ranges that we observe here matches my experience of improvements that I have seen on my projects (i.e. 1,5-4 times). What I see as potential way to speed up code (which may work, or not work out) is to try using numpy vectors for data transfers and vectorized computations, and using scipy for solving ODE equations (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK23.html). Some pointers could be taken from here (https://scicomp.stackexchange.com/questions/21060/runge-kutta-simulation-for-projectile-motion-with-drag), but without own implementation of Runge-Kutta, as scipy implementations are heavily tested and optimized. Some interesting points on comparison of scipy ODE could be found in this reddit thread (https://www.reddit.com/r/CasualMath/comments/olyrxw/i_did_a_casual_survey_of_scipy_ode_solvers/?rdt=51910). I believe that this work could be done as a separate Calculator. Such an implementation may have benefits from speed perspective (although everything needs to be profiled and measured), and may have drawbacks from maintainability perspective (as it is usual with heavily optimized code).

@serhiy-yevtushenko I expected those advantages from using SciPy, and I did try to get their RK45 method working but couldn't follow far enough to get it working with reasonable performance – almost certainly due to my lack of understanding. Can you try an implementation of it? You could see a kludgy effort to hack a different calculator engine into the existing trajectory_calc.py at line 360 here. Or, today I might have a separate implementation in the current master version's rk4.py stub, but it will not be optimized or immediately offer the error measurement and control that should come for free in SciPy. I would rather spend what little free time I have on other features than reimplementing classic numerical methods....

[serhiy-yevtushenko]

Hi, @dbookstaber . I'm willing to try to get the scheme with scipy working, as I'm interested in getting better performance. What I need the help with is having the sensible set of unit tests against which I could check the implementation. I'm willing to procede as far as too see, whether it really able to bring the performance improvements

[dbookstaber]

Short question - I would like to do an experiment with implementing a TrajectoryCalculator with ivp_solve from scipy. Which unit tests are sensible to use as a starting point and is there a recommended sequence of unit test from functionality progression viewpoint (simplest case -> more complex parameters like handling of air pressure, winds, ...)?

My opinion: Start with test_trajectory.py. If you can get a basic no-wind trajectory working correctly, all of the other stuff should come for free with the existing codebase!

Note also that JBM's calculator is presumed to be correct, so you can generate arbitrary scenarios and check against his results.

[o-murphy]

@serhiy-yevtushenko
no, we don't have any comments on this at the moment. Maybe @dbookstaber has them.

Regarding any contributions to TrajectoryCalc, please don't edit TrajectoryCalc directly, don't break backward compatibility as much as possible. This creates a lot of issues that need to be resolved later, from compatibility with cython, to specific dependencies for the contribution, like numpy, scipy etc. Which definitely won't be added as a base library dependency. They can only be installed together with an extension module (like py_ballisticcalc.exts)

Use this sample for a new implementation

# my_trajectory_calc.py
# type: ignore

from typing import List

from py_ballisticcalc import TrajectoryCalc, Distance, Shot, TrajectoryData, Calculator


class MyTrajectoryCalc(TrajectoryCalc):

    ... # overrides

    def trajectory(self, shot_info: Shot, max_range: Distance, dist_step: Distance,
                   extra_data: bool = False, time_step: float = 0.0) -> List[TrajectoryData]:
        # your implementation there
        ...

class MyCalculator(Calculator):
     ... # overrides

[serhiy-yevtushenko]

My Idea is to try with derived class from my calc. The only point to the new sample - is it ok to use python 3.11 style type annotations? (I mean using list instead of List and then not importing from typing types, which already exist in python, like list, set, and other ones?

[dbookstaber]

My Idea is to try with derived class from my calc. The only point to the new sample - is it ok to use python 3.11 style type annotations? (I mean using list instead of List and then not importing from typing types, which already exist in python, like list, set, and other ones?

Dmitry is explicitly supporting 3.9-3.12, but for extensions you can use whatever version; he'll adjust as necessary when incorporated to releases.

[dbookstaber]

Going forward, Dmytro wants to separate calculation engines from the core py_ballisticcalc project. As part of this refactoring we split the RK-based calculator into https://github.com/dbookstaber/RKballistic.

[o-murphy]

Later I will add the link accordingly in README.md