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, ARM32 (Cortex-A and Cortex-M), ARM64, RISC-V (32 and 64 Bit) and TriCore (AURIX)
• 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

The computation engine is fully functional - for all above mentioned target architectures - and available for testing and experimentation simply by installing the package.

Hardware I/O and real-time execution are available with the runner for the Zephyr RTOS. It executes periodic tasks with deterministic timing (logical execution time): a task reads its inputs at its release and publishes its outputs at the end of its period. The host is connected over a serial port (e.g. USB) by cp.SerialTarget to deploy the tasks, start and stop them, and to read and write values while they are running. Device I/O channels - GPIO, PWM, quadrature encoder and ADC - are bound to the values of the tasks. If the tasks need more CPU time than available, the tasks with the lowest priority are degraded to keep the connection to the host responsive. Supported boards are currently the Teensy 4.1 (Cortex-M7).

Values of periodic tasks are signals in concept: one value per run of the task. cp.delay, cp.integrate, cp.rate and cp.feedback operate on these signals without side effects. Their state is kept in the memory of the task, cp.dt is the period of the task. Orientations as quaternions are integrated from angular velocities by cp.integrate_rotation(gyro), cp.angular_rate(q) is the angular velocity of an orientation:

setpoint = cp.value(1.0)
position = cp.value(0.0)  # e.g. bound to an encoder with tg.bind
error = setpoint - position
control = 2.0 * error + 5.0 * cp.integrate(error) + 0.1 * cp.rate(error)
smoothed = cp.feedback(lambda prev: prev + 0.2 * (control - prev), init=0.0)

tg = cp.LocalTarget(tick=1e-3)
tg.deploy(cp.Task(smoothed, period=1e-3))
tg.step(0.1)  # 100 runs in simulated time

Array stencils for handling large arrays - e.g., for machine vision and neural network applications - are available and use auto-vectorizing of the C compiler. Automatic differentiation covers element-wise array operations, reductions and matrix products, but not yet convolutions, pooling, sorting, transposing and slicing of arrays.

Furthermore in development is currently constant regrouping for further symbolic optimization of the computation graph.

Benchmark performance shows promising numbers. The following chart compares three ways of running the same computation with Copapy and NumPy 2.0.1:

For the benchmark (tests/benchmark.py) the median time per iteration of calculating the term sum((v1 + i) @ v2 for i in range(10)) was measured on a Ryzen 5 3400G (15 runs of 30,000 iterations each). 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 all implementations. The Copapy variants are:

• Copapy scalar: cp.vector(..., packed=False), one scalar stencil per element. The code grows with v_size but they profit from sparcity (no benefit in this benchmark).
• Copapy packed: cp.vector(..., packed=True), array stencils that are auto-vectorized by the C compiler. They compile much faster for large arrays. Without the packed argument, vectors with more than 64 elements are packed automatically.
• Copapy jit: the term as a function decorated with @cp.jit, called with the two NumPy arrays. Unlike the other variants, this includes copying the input arrays and reads back the result.

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 numpy as np
import copapy as cp

@cp.jit
def rms_error(measured, reference):
    error = measured - reference
    return cp.sqrt((error * error).sum() / error.size)

@cp.jit
def normalize(v):
    return cp.map_to_numpy(v / cp.sqrt((v * v).sum()))

measured = np.array([1.0, 2.1, 2.9, 4.2], dtype=np.float32)
reference = np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float32)

# Compile and run:
rms = rms_error(measured, reference)  # 0.1225

# Run cached compiled version:
rms = rms_error(measured * 2, reference)

# Array result as NumPy array:
unit = normalize(reference)  # array([0.183, 0.365, 0.548, 0.730], dtype=float32)

Arguments can be numbers, sequences of numbers, or arrays with the buffer protocol, such as NumPy arrays. A function is compiled again for each new combination of argument types and array shapes. NumPy arrays of byte, 32 bit int or 32 bit float are copied to the compiled code as they are, other numeric types are converted to 32 bit. Array results are returned as nested lists, or as NumPy arrays if marked with cp.map_to_numpy. Such a NumPy array is a view of the memory of the compiled code and changes with the next call, unless copy=True is passed.

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: vector((-0.7229549288749695, 2.6246931552886963))
Joint position: vector((1.3498780727386475, -1.1907261610031128))
End-effector position: vector((0.6998459100723267, 0.7006908655166626))
quadratic error = 5.010516588299652e-07

Real-time tasks on a device

Periodic tasks run in real time on a device with the Zephyr runner. This example deploys a position controller to a Teensy 4.1. Every millisecond it reads an encoder and sets the duty cycle of a PWM output, while the host changes the setpoint:

import time
import copapy as cp

# Connect to a Teensy 4.1 running the Zephyr runner (port found by USB id)
tg = cp.SerialTarget(device=cp.devices.teensy41)

# Position controller: encoder in, PWM duty cycle out
setpoint = cp.value(0.0)  # written by the host while running
position = cp.value(0)    # read from the encoder at each release
tg.bind(position, tg.io.encoder('enc0'))

error = setpoint - position * 0.01
duty = cp.clamp(0.5 + error * 2.0, 0.0, 1.0)
tg.bind(duty, tg.io.pwm('pwm2'))

controller = cp.Task(duty, period=1e-3)  # released every millisecond
tg.deploy(controller)
tg.run()  # real-time execution on the device

tg.write_value(setpoint, 0.1)
time.sleep(1.0)
print(f"Duty cycle: {tg.read_value(duty)}")
print(f"Jobs run: {tg.stats(controller)['runs']}")
tg.stop()

Duty cycle: 0.699999988079071
Jobs run: 1002

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. For the Cortex-M targets with a single precision FPU most math functions are based on CMSIS-DSP instead, which avoids the software emulated double precision arithmetic of the MUSL functions.

Array operations use the same mechanism. Instead of passing operands in registers, an array stencil accesses its operands and its result in heap memory through placeholder symbols, which are patched to the actual addresses. The loop over the elements is in a non-stencil function shared by all stencils of that operation:

add_floatarr_floatarr(void) {
    aux_arr_add_float_float_vv(REF(ref_arg0), REF(ref_arg1), REF(ref_out), *(int *)REF(ref_arg2));
    result_float_ref();
}

Here ref_arg0 and ref_arg1 are the two arrays, ref_arg2 is the number of elements and ref_out is the result array. The name of the dummy function result_float_ref() tells the compiler the element type of the result.

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

The runner for the Zephyr RTOS (targets/zephyr) is built in a Zephyr workspace (tested with Zephyr 4.4 and the Zephyr SDK, ZEPHYR_BASE set):

bash tools/build_zephyr.sh mps2/an500  # Cortex-M7 in qemu-system-arm
bash tools/build_zephyr.sh teensy41    # Teensy 4.1, runner on the USB port

The tests tests/test_runner_zephyr.py and tests/test_serial_target.py run it in qemu (path of qemu-system-arm in the environment variable QEMU_SYSTEM_ARM) or on a board (serial port in COPARUN_SERIAL_PORT).

License

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

Footnotes

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), ARMv7 Thumb for Cortex-A and Cortex-M (M4 and M7), ARMv8.1-M Thumb with Helium (Cortex-M55 and M85), RISC-V (RV32/RV64 with hardware floating point) and TriCore 1.6.1 (AURIX TC2xx). ↩

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. ↩