Copapy

Copapy is a Python framework for deterministic, low-latency real-time computation with automatic differentiation support, targeting hardware applications - for example in the fields of robotics, aerospace, SDR, embedded systems and control systems in general.

GPU frameworks like PyTorch, JAX and TensorFlow jump-started the development in the field of AI. With the right balance of flexibility and performance, they allow for fast iteration of new ideas while still being performant enough to test or even use them in production.

This is exactly what Copapy aims for - but in the field of embedded realtime computation. While making use of the ergonomics of Python, the tooling, and the general Python ecosystem, Copapy runs seamlessly optimized machine code. Despite being highly portable, the copy-and-patch compiler allows for effortless and fast deployment without any dependencies beyond Python. It's designed to feel like writing Python scripts with a shallow learning curve, but under the hood it produces high-performance, statically typed and memory-safe code with a minimized set of possible runtime errors¹. To maximize productivity, the framework provides detailed type hints to catch most errors even before compilation.

Embedded systems come with a variety of CPU architectures. The copy-and-patch compiler already supports the most common ones², and porting it to new architectures is straightforward if a C compiler for the target architecture is available³. The generated code depends only on the CPU architecture. The generated binaries neither perform system calls nor rely on external libraries like libc. This makes Copapy both highly deterministic and easy to deploy on different realtime operating systems (RTOS) or bare metal.

The main features can be summarized as:

• Fast to write & easy to read
• Memory and type safety with a minimal set of runtime errors
• Deterministic execution
• Automatic differentiation for efficient realtime optimization (reverse-mode)
• Optimized machine code for x86_64, 32 Bit ARM (Cortex-A and Cortex-M) and AArch64
• Highly portable to new architectures
• Small Python package with minimal dependencies and no cross-compile toolchain required

Execution of the compiled code is managed by a runner application. The runner is implemented in C and handles I/O and communication with the Copapy framework. The overall design emphasizes minimal complexity of the runner to simplify portability, since this part must be adapted for the individual hardware/application. Because patching of memory addresses is done by the runner, the different architecture-specific relocation types are unified to an architecture-independent format by Copapy before sending the patch instructions to the runner. This keeps the runner implementation as minimal as possible.

The design targets either an architecture with a realtime-patched Linux kernel - where the runner uses the same CPU and memory as Linux but executes in a realtime thread - or a setup where even higher determinism is required. In such cases, the runner can be executed on a separate crossover MCU running on bare metal or a RTOS.

The Copapy framework also includes a runner as a Python module built from the same C code. This allows frictionless testing of code and might be valuable for using Copapy in conventional application development.

Current state

While hardware I/O is obviously a core aspect of the project it is not yet available. However the computation engine is already fully functional - for all above mentioned target architectures - and available for testing and experimentation simply by installing the package. The project focuses now on integration into the first demonstration hardware platform.

Furthermore in development are currently:

• Array stencils for handling large arrays and generating SIMD-optimized code - e.g., for machine vision and neural network applications
• Constant regrouping for further symbolic optimization of the computation graph

Despite missing SIMD-optimization, benchmark performance shows promising numbers. The following chart plots the results in comparison to NumPy 2.3.5:

For the benchmark (tests/benchmark.py) timings for 30,000 iterations of calculating the term sum((v1 + i) @ v2 for i in range(10)) were measured on a Ryzen 5 3400G. The vectors v1 and v2 both have lengths of v_size, which was varied from 10 to 500 according to the chart. For the NumPy case the i in range(10) loop was vectorized like this: np.sum((v1 + i) @ v2) with i being an NDArray of shape [10, 1]. The number of calculated scalar operations is the same for both implementations. Copapy benefits from lower overhead by calling a single function from Python per iteration, whereas the NumPy variant requires three. Interestingly, the chart shows no indication that for increasing v_size the calling overhead for NumPy will be compensated by faster SIMD instructions. Note that in this benchmark the Copapy case does not move any data between Python and the compiled code.

Furthermore, for many applications Copapy performance will benefit by reducing the actual number of operations significantly compared to a NumPy implementation, by precomputing constant values known at compile time and benefiting from sparsity. Multiplying by zero (e.g., in a diagonal matrix) eliminates a whole branch in the computation graph. Operations with no effect, like multiplications by 1 or additions with zero, get eliminated at compile time.

For testing and using Copapy to speed up computations in conventional Python programs there is also the @cp.jit decorator available, to compile functions on first use and cache the compiled version for later calls:

import copapy as cp

@cp.jit
def calculation(x: float, y: float) -> float:
    return sum(x ** 2 + y ** 2 + i for i in range(10))

# Compile and run:
result1 = calculation(2.5, 1.2)

# Run cached compiled version:
result2 = calculation(3.1, 4.7)

Note that cp.jit is not currently highly optimized for data transfer between Python and the compiled code.

Install

To install Copapy, you can use pip. Precompiled wheels are available for Linux (x86_64, AArch64, ARMv7, i686), Windows (x86_64) and macOS (x86_64, AArch64):

pip install copapy

Examples

Basic example

A very simple example program using Copapy can look like this:

import copapy as cp

# Define variables
a = cp.value(0.25)
b = cp.value(0.87)

# Define computations
c = a + b * 2.0
d = c ** 2 + cp.sin(a)
e = cp.sqrt(b)

# Create a target (default is local), compile and run
tg = cp.Target()
tg.compile(c, d, e)
tg.run()

# Read the results
print("Result c:", tg.read_value(c))
print("Result d:", tg.read_value(d))
print("Result e:", tg.read_value(e))

IMU Sensor fusion

This example demonstrates how concisely the Madgwick AHRS algorithm for IMU sensor fusion can be expressed in Copapy. While the original C implementation contains more than 130 lines of code (excluding blank lines and comments), the algorithm can be expressed in about 10 lines with Copapy. One reason is Copapy's automatic differentiation, which computes the Madgwick correction gradient with cp.grad():

import copapy as cp
from copapy import vector, quaternion

def madgwick_imu_update(q: quaternion, gyro: vector[float],
                        accel: vector[float]) -> quaternion:
    BETA: float = 0.1
    SAMPLE_INTERVAL: float = 0.01  # seconds
    objective = q.rotate_vector(vector([0.0, 0.0, 1.0])) - accel.normalize()
    cost = 0.5 * objective.dot(objective)
    gradient = cp.grad(cost, q).normalize()
    gyro_quat = quaternion(0.0, *gyro)
    q_dot_gyro = 0.5 * (q @ gyro_quat)
    q_dot = q_dot_gyro - BETA * gradient
    return (q + q_dot * SAMPLE_INTERVAL).normalize()

# Usage
q = quaternion(cp.value(0.7071), cp.value(0.7071), cp.value(0.0), cp.value(0.0))
gyro = vector([0.01, 0.02, 0.015])
accel = vector([0.0, 0.0, 1.0])

new_q = madgwick_imu_update(q, gyro, accel)

tg = cp.Target()
tg.compile(new_q)
tg.run()

print(f"Updated orientation: {tg.read_value(new_q)}")

See tests/test_quaternion.py for the documented version.

Inverse kinematics

Another example using autograd in Copapy, here implementing gradient descent to solve an inverse kinematics problem for a two-joint 2D arm:

import copapy as cp

# Arm lengths
l1, l2 = 1.8, 2.0

# Target position
target = cp.vector([0.7, 0.7])

# Learning rate for iterative adjustment
alpha = 0.1

def forward_kinematics(theta1, theta2):
    """Return positions of joint and end-effector."""
    joint = cp.vector([l1 * cp.cos(theta1), l1 * cp.sin(theta1)])
    end_effector = joint + cp.vector([l2 * cp.cos(theta1 + theta2),
                                     l2 * cp.sin(theta1 + theta2)])
    return joint, end_effector

# Start values
theta = cp.vector([cp.value(0.0), cp.value(0.0)])

# Iterative inverse kinematics
for _ in range(48):
    joint, effector = forward_kinematics(theta[0], theta[1])
    error = ((target - effector) ** 2).sum()

    theta -= alpha * cp.grad(error, theta)

tg = cp.Target()
tg.compile(error, theta, joint)
tg.run()

print(f"Joint angles: {tg.read_value(theta)}")
print(f"Joint position: {tg.read_value(joint)}")
print(f"End-effector position: {tg.read_value(effector)}")
print(f"quadratic error = {tg.read_value(error)}")

Joint angles: [-0.7221821546554565, 2.6245293617248535]
Joint position: [1.3509329557418823, -1.189529299736023]
End-effector position: [0.6995794177055359, 0.7014330625534058]
quadratic error = 2.2305819129542215e-06

How it works

The compilation step starts with tracing the Python code to generate an acyclic directed graph (DAG) of variables and operations. The code can contain functions, closures, branching, and so on, but conditional branching is only allowed when the condition is known at tracing time (a cp.iif function exists to work around this). In the next step, this DAG is optimized and linearized into a sequence of operations. Each operation is mapped to a precompiled stencil or a combination of several stencils. A stencil is a piece of machine code with placeholders for memory addresses pointing to other code or data. The compiler generates patch instructions that fill these placeholders with the correct memory addresses.

After compilation, the binary code built from the stencils, the constant data, and the patch instructions is handed to the runner for execution. The runner allocates memory for code and data, copies both into place, applies the patch instructions, and finally executes the code.

The C code for a very simple stencil can look like this:

add_float_float(float arg1, float arg2) {
    result_float_float(arg1 + arg2, arg2);
}

The call to the dummy function result_float_float ensures that the compiler keeps the result and the second operand in registers for later use. The dummy function acts as a placeholder for the next stencil. Copapy uses two virtual registers, which map on most relevant architectures to actual hardware registers. Data that cannot be kept in a register is stored in statically allocated heap memory. Stack memory may be used inside some stencils, but its usage is essentially fixed and independent of the Copapy program, so total memory requirements are known at compile time.

The machine code for the function above, compiled for x86_64, looks like this:

0000000000000000 <add_float_float>:
   0:    f3 0f 58 c1              addss  %xmm1,%xmm0
   4:    e9 00 00 00 00           jmp    9 <.LC1+0x1>
            5: R_X86_64_PLT32    result_float_float-0x4

Based on the relocation entry for the jmp to the symbol result_float_float, the jmp instruction is stripped when it is the last instruction in a stencil. Thus, a Copapy addition operation results in a single instruction. For stencils containing multiple branch exits, only the final jmp is removed; the others are patched to jump to the next stencil.

For more complex operations - where inlining is less useful - stencils call a non-stencil function, such as in this example:

0000000000000000 <sin_float>:
   0:    48 83 ec 08              sub    $0x8,%rsp
   4:    e8 00 00 00 00           call   9 <sin_float+0x9>
            5: R_X86_64_PLT32    sinf-0x4
   9:    48 83 c4 08              add    $0x8,%rsp
   d:    e9 00 00 00 00           jmp    12 <.LC0+0x2>
            e: R_X86_64_PLT32    result_float-0x4

Unlike stencils, non-stencil functions like sinf are not stripped and do not need to be tail-call-optimizable. These functions can be provided as C code and compiled together with the stencils or can be object files like in the case of sinf compiled from C and assembly code and merged into the stencil object files. Math functions like sinf are currently provided by the MUSL C library, with architecture-specific optimizations.

Non-stencil functions and constants used by them are stored together with the stencils in ELF object files for each supported CPU architecture. The required non-stencil functions and constants are bundled during compilation. The compiler includes only the data and code required for a specific Copapy program.

The Copapy compilation process is independent of the actual instruction set. It relies purely on relocation entries and symbol metadata from the ELF file generated by the C compiler.

Developer Guide

Feedback and contributions are welcome - please open an issue or submit a pull request on GitHub.

To get started with development, first clone the repository:

git clone https://github.com/Nonannet/copapy.git
cd copapy

You may set up a virtual environment:

python -m venv .venv
source .venv/bin/activate  # On Windows: `.venv\Scripts\activate`

Build and install the package and dev dependencies:

pip install -e .[dev]

If the build fails because no suitable C compiler is installed, you can either install one or use the binary package from PyPI:

pip install copapy[dev]

When running pytest, it will use the binary components from PyPI, but all Python code is executed from the local repository.

To run all tests, you need the stencil object files and the compiled runner. You can download them from GitHub or build them yourself with gcc.

Download the latest binaries from GitHub:

python tools/get_binaries.py

Build the binaries from source on Linux:

bash tools/build.sh

Run the tests:

pytest

License

This project is licensed under the MIT license - see the LICENSE file for details.

1. Errors like divide-by-zero are currently still possible. The feasibility of tracking value ranges in the type system is under investigation to enable compile-time checks. ↩

2. Supported architectures: x86_64, x86 (32 Bit), AArch64, ARMv6/7 (non-Thumb) and ARMv7 Thumb for Cortex-A and Cortex-M. ↩

3. The compiler must support tail-call optimization (TCO). Currently, GCC is supported. Porting to a new architecture requires implementing a subset of relocation types used by that architecture. ↩