Python ENnUI Demo¶

Demonstrate calling compiled ennui with nanobind bindings from python as pyennui. Compare results against pure-python implementation both for accuracy and computational cost.

In [1]:

# Setup
import sys
import os

# C++ compiled with nanobind bindings
# sys.path.append(os.path.join("..", "..", "build", "bin", "Release"))
import pyennui

# Pure-python implementation
sys.path.append(os.path.join("..", "..", "python", "reference"))
import pure_ennui

import numpy as np
from timeit import default_timer as timer

error_abs = lambda x, y: np.linalg.norm(x - y, np.inf)
error_rel = lambda x, y: np.linalg.norm(
    (x - y) / (np.maximum(np.abs(x), np.abs(y))), np.inf
)
print_error_abs = lambda str, x, y: print(
    "Rel. err., %s : %.2e" % (str, error_abs(x, y))
)
print_error_rel = lambda str, x, y: print(
    "Abs. err., %s : %.2e" % (str, error_rel(x, y))
)
print_time = lambda str, start, end: print("%s : %e" % (str, end - start))

tol = 1e-14

# Setup import sys import os # C++ compiled with nanobind bindings # sys.path.append(os.path.join("..", "..", "build", "bin", "Release")) import pyennui # Pure-python implementation sys.path.append(os.path.join("..", "..", "python", "reference")) import pure_ennui import numpy as np from timeit import default_timer as timer error_abs = lambda x, y: np.linalg.norm(x - y, np.inf) error_rel = lambda x, y: np.linalg.norm( (x - y) / (np.maximum(np.abs(x), np.abs(y))), np.inf ) print_error_abs = lambda str, x, y: print( "Rel. err., %s : %.2e" % (str, error_abs(x, y)) ) print_error_rel = lambda str, x, y: print( "Abs. err., %s : %.2e" % (str, error_rel(x, y)) ) print_time = lambda str, start, end: print("%s : %e" % (str, end - start)) tol = 1e-14

Help¶

Automated help displays available methods and arguments

In [2]:

help(pyennui)

help(pyennui)

Help on module pyennui:

NAME
    pyennui

SUBMODULES
    geodetic
    math
    mechanization

DATA
    normalize_SO3_return_arg = <nanobind.nb_func object>
        normalize_SO3_return_arg(arg0: numpy.ndarray[dtype=float64, writable=False, shape=(3, 3), order='*'], arg1: numpy.ndarray[dtype=float64, shape=(3, 3), order='*'], /) -> None

FILE
    c:\users\mrwalke\appdata\local\anaconda3\envs\ennui_dev\lib\site-packages\pyennui.pyd

Validate gravitation model¶

gravitation_ecef returns gravitation vector as a function of ECEF coordinates. Contrast against reference and pure-python implementations.

In [3]:

WhiteHouse_ECEF = np.array(
    (1.1150423452941689e06, -4.8438122981491517e06, 3.9835202164462707e06)
)
WhiteHouse_gamma = np.array(
    (-1.7170260919766687e00, 7.4588665943134185e00, -6.1541304311837033e00)
)
p = WhiteHouse_ECEF

N = 1000

# Pure-python
start = timer()
for _ in range(N):
    gamma_pure = pure_ennui.geodetic.gravitation_ecef(p)
end = timer()
print_time("Timing, pure-python", start, end)

# C++
start = timer()
for _ in range(N):
    gamma = pyennui.geodetic.gravitation_ecef(p)
end = timer()
print_time("Timing, C++        ", start, end)

print_error_abs("C++ vs. pure-py", gamma, gamma_pure)
print_error_abs("C++ vs. ref.   ", gamma, WhiteHouse_gamma)
assert error_rel(gamma, WhiteHouse_gamma) < tol, "C++ does not match reference."

WhiteHouse_ECEF = np.array( (1.1150423452941689e06, -4.8438122981491517e06, 3.9835202164462707e06) ) WhiteHouse_gamma = np.array( (-1.7170260919766687e00, 7.4588665943134185e00, -6.1541304311837033e00) ) p = WhiteHouse_ECEF N = 1000 # Pure-python start = timer() for _ in range(N): gamma_pure = pure_ennui.geodetic.gravitation_ecef(p) end = timer() print_time("Timing, pure-python", start, end) # C++ start = timer() for _ in range(N): gamma = pyennui.geodetic.gravitation_ecef(p) end = timer() print_time("Timing, C++ ", start, end) print_error_abs("C++ vs. pure-py", gamma, gamma_pure) print_error_abs("C++ vs. ref. ", gamma, WhiteHouse_gamma) assert error_rel(gamma, WhiteHouse_gamma) < tol, "C++ does not match reference."

Timing, pure-python : 4.455900e-03
Timing, C++         : 1.283100e-03
Rel. err., C++ vs. pure-py : 0.00e+00
Rel. err., C++ vs. ref.    : 1.78e-15

Matrix handling¶

Demonstrate passing matrix through use of normalize_SO3 which orthonormalizes a rotation matrix using singular-value decomposition.

In [4]:

m = np.eye(3)
m[0, 1] = 1e-3
print_error_abs("pre-orthonormalization ", m @ m.transpose(), np.eye(3))
m = pyennui.math.normalize_SO3(m)
print_error_abs("post-orthonormalization", m @ m.transpose(), np.eye(3))
assert error_abs(m @ m.transpose(), np.eye(3)) < tol, "C++ does not match reference."

m = np.eye(3) m[0, 1] = 1e-3 print_error_abs("pre-orthonormalization ", m @ m.transpose(), np.eye(3)) m = pyennui.math.normalize_SO3(m) print_error_abs("post-orthonormalization", m @ m.transpose(), np.eye(3)) assert error_abs(m @ m.transpose(), np.eye(3)) < tol, "C++ does not match reference."

Rel. err., pre-orthonormalization  : 1.00e-03
Rel. err., post-orthonormalization : 1.11e-16

Comparable implementation returning as an argument.

In [5]:

# Using a return argument
m = np.eye(3)
m[0, 1] = 1e-3
rslt = np.ndarray((3, 3))

print_error_abs("pre-orthonorm ", m @ m.transpose(), np.eye(3))
pyennui.normalize_SO3_return_arg(m, rslt)
m = rslt
print_error_abs("post-orthonorm", m @ m.transpose(), np.eye(3))

# Using a return argument m = np.eye(3) m[0, 1] = 1e-3 rslt = np.ndarray((3, 3)) print_error_abs("pre-orthonorm ", m @ m.transpose(), np.eye(3)) pyennui.normalize_SO3_return_arg(m, rslt) m = rslt print_error_abs("post-orthonorm", m @ m.transpose(), np.eye(3))

Rel. err., pre-orthonorm  : 1.00e-03
Rel. err., post-orthonorm : 1.11e-16

State propagation¶

Define initial state, inertial measurements, elapsed time, and final state as computed using a reference algorithm.

In [6]:

WhiteHouse_mean_prop = {
    "prior": {
        "position": np.array(
            [1.115042345294169e06, -4.843812298149152e06, 3.983520216446271e06]
        ),
        "velocity": np.array(
            [1.700252783993267e00, 5.799253612971604e00, -7.143100966293305e00]
        ),
        "attitude": np.matrix(
            [
                [
                    -1.689390579588263e-01,
                    -4.453768089569571e-01,
                    -8.792605374627603e-01,
                ],
                [8.000355898459490e-01, -5.830041132072308e-01, 1.415954058693114e-01],
                [-5.756758199506289e-01, -6.795187282384265e-01, 4.548094637289362e-01],
            ]
        ),
    },
    "specific_force": np.array(
        [-2.351012554013630e00, 1.857770394070348e00, 1.940270045514844e01]
    ),
    "angular_rate": np.array(
        [6.941949163919531e00, -3.383664771751461e00, 7.607424972048551e00]
    ),
    "dt": 1.000000000000000e-03,
    "posterior": {
        "position": np.array(
            [1.115042346984841e06, -4.843812292346284e06, 3.983520209304588e06]
        ),
        "velocity": np.array(
            [1.681091028114939e00, 5.806482921506987e00, -7.140263743663021e00]
        ),
        "attitude": np.matrix(
            [
                [
                    -1.753142992423087e-01,
                    -4.501584286525450e-01,
                    -8.755696920258544e-01,
                ],
                [7.960624873695160e-01, -5.880875442140073e-01, 1.429599823146225e-01],
                [-5.792662709706455e-01, -6.719452577802820e-01, 4.614543941305069e-01],
            ]
        ),
    },
}

# Pre-allocate output
posterior = {
    "position": np.empty((3)),
    "velocity": np.empty((3)),
    "attitude": np.empty((3, 3)),
}

# Pre-compute gravitation vector
gamma = pyennui.geodetic.gravitation_ecef(WhiteHouse_mean_prop["prior"]["position"])

WhiteHouse_mean_prop = { "prior": { "position": np.array( [1.115042345294169e06, -4.843812298149152e06, 3.983520216446271e06] ), "velocity": np.array( [1.700252783993267e00, 5.799253612971604e00, -7.143100966293305e00] ), "attitude": np.matrix( [ [ -1.689390579588263e-01, -4.453768089569571e-01, -8.792605374627603e-01, ], [8.000355898459490e-01, -5.830041132072308e-01, 1.415954058693114e-01], [-5.756758199506289e-01, -6.795187282384265e-01, 4.548094637289362e-01], ] ), }, "specific_force": np.array( [-2.351012554013630e00, 1.857770394070348e00, 1.940270045514844e01] ), "angular_rate": np.array( [6.941949163919531e00, -3.383664771751461e00, 7.607424972048551e00] ), "dt": 1.000000000000000e-03, "posterior": { "position": np.array( [1.115042346984841e06, -4.843812292346284e06, 3.983520209304588e06] ), "velocity": np.array( [1.681091028114939e00, 5.806482921506987e00, -7.140263743663021e00] ), "attitude": np.matrix( [ [ -1.753142992423087e-01, -4.501584286525450e-01, -8.755696920258544e-01, ], [7.960624873695160e-01, -5.880875442140073e-01, 1.429599823146225e-01], [-5.792662709706455e-01, -6.719452577802820e-01, 4.614543941305069e-01], ] ), }, } # Pre-allocate output posterior = { "position": np.empty((3)), "velocity": np.empty((3)), "attitude": np.empty((3, 3)), } # Pre-compute gravitation vector gamma = pyennui.geodetic.gravitation_ecef(WhiteHouse_mean_prop["prior"]["position"])

In [7]:

# Call C++ state-propagation

N = 1000

# Pure-python
start = timer()
for _ in range(N):
    pyennui.mechanization.ecef.fwd_pva_S03(
        WhiteHouse_mean_prop["prior"]["position"],
        WhiteHouse_mean_prop["prior"]["velocity"],
        WhiteHouse_mean_prop["prior"]["attitude"],
        gamma,
        WhiteHouse_mean_prop["specific_force"],
        WhiteHouse_mean_prop["angular_rate"],
        WhiteHouse_mean_prop["dt"],
        posterior["position"],
        posterior["velocity"],
        posterior["attitude"],
    )

end = timer()
print_time("Timing, C++", start, end)

# Print results and error if tol exceeded
for key in posterior:
    rslt = posterior[key]
    expected = WhiteHouse_mean_prop["posterior"][key]
    print_error_rel("%s" % key[0], rslt, expected)
    assert error_rel(rslt, expected) < tol, "C++ does not match reference."
print("Success!")

# Pure-python
start = timer()
for _ in range(N):
    [posterior["position"], posterior["velocity"], posterior["attitude"]] = (
        pure_ennui.mechanization.ecef.fwd_pva_S03(
            WhiteHouse_mean_prop["prior"]["position"],
            WhiteHouse_mean_prop["prior"]["velocity"],
            WhiteHouse_mean_prop["prior"]["attitude"],
            gamma,
            WhiteHouse_mean_prop["specific_force"],
            WhiteHouse_mean_prop["angular_rate"],
            WhiteHouse_mean_prop["dt"],
        )
    )

end = timer()
print_time("Timing, pure-python", start, end)

# Print results and error if tol exceeded
for key in posterior:
    rslt = posterior[key]
    expected = WhiteHouse_mean_prop["posterior"][key]
    print_error_rel("%s" % key[0], rslt, expected)
    assert error_rel(rslt, expected) < tol, "C++ does not match reference."
print("Success!")

# Call C++ state-propagation N = 1000 # Pure-python start = timer() for _ in range(N): pyennui.mechanization.ecef.fwd_pva_S03( WhiteHouse_mean_prop["prior"]["position"], WhiteHouse_mean_prop["prior"]["velocity"], WhiteHouse_mean_prop["prior"]["attitude"], gamma, WhiteHouse_mean_prop["specific_force"], WhiteHouse_mean_prop["angular_rate"], WhiteHouse_mean_prop["dt"], posterior["position"], posterior["velocity"], posterior["attitude"], ) end = timer() print_time("Timing, C++", start, end) # Print results and error if tol exceeded for key in posterior: rslt = posterior[key] expected = WhiteHouse_mean_prop["posterior"][key] print_error_rel("%s" % key[0], rslt, expected) assert error_rel(rslt, expected) < tol, "C++ does not match reference." print("Success!") # Pure-python start = timer() for _ in range(N): [posterior["position"], posterior["velocity"], posterior["attitude"]] = ( pure_ennui.mechanization.ecef.fwd_pva_S03( WhiteHouse_mean_prop["prior"]["position"], WhiteHouse_mean_prop["prior"]["velocity"], WhiteHouse_mean_prop["prior"]["attitude"], gamma, WhiteHouse_mean_prop["specific_force"], WhiteHouse_mean_prop["angular_rate"], WhiteHouse_mean_prop["dt"], ) ) end = timer() print_time("Timing, pure-python", start, end) # Print results and error if tol exceeded for key in posterior: rslt = posterior[key] expected = WhiteHouse_mean_prop["posterior"][key] print_error_rel("%s" % key[0], rslt, expected) assert error_rel(rslt, expected) < tol, "C++ does not match reference." print("Success!")

Timing, C++ : 3.068800e-03
Abs. err., p : 2.34e-16
Abs. err., v : 3.96e-16
Abs. err., a : 5.26e-16
Success!
Timing, pure-python : 1.227745e-01
Abs. err., p : 2.34e-16
Abs. err., v : 3.96e-16
Abs. err., a : 5.22e-16
Success!