Operator syntax

The operator API provides a more compositional interface to gratopy’s projection operators. It is currently experimental and primarily focused on parallel-beam Radon transforms.

The legacy gratopy.ProjectionSettings API remains the main documented interface for the full feature set of gratopy. The operator API complements it with a syntax that is often more convenient when working with operator algebra, adjoint operators, and experimental kernels.

Warning

The operator API is experimental. Backward-incompatible changes may be introduced without a full deprecation cycle while the interface and internal abstractions are still settling.

Extent placeholders such as gratopy.utilities.ExtentPlaceholder are supported experimentally for Radon operators when exactly one of the image or detector extents is a placeholder. Passing placeholders for both extents at once remains unsupported and raises NotImplementedError.

Current scope

The current operator API supports in particular:

• the gratopy.operator.projection.Radon operator,

• adjoints via T,

• operator composition and arithmetic,

• custom OpenCL kernels via gratopy.operator.opencl.OpenCLKernelSpec.

At the moment, the operator API should be understood as Radon-only for execution. The exported Fanbeam class is currently a non-executable placeholder for the future operator port; use the legacy API for fanbeam projections.

Quick example

A Radon transform and its adjoint can be used as follows:

import numpy as np
import pyopencl as cl
import pyopencl.array as clarray
import gratopy

ctx = cl.create_some_context(interactive=False)
queue = cl.CommandQueue(ctx)

Nx = 128
img = np.zeros((Nx, Nx), dtype=np.float32)

R = gratopy.operator.Radon(image_domain=Nx, angles=180)

sino = R.apply_to(img, queue=queue)
backprojection = R.T.apply_to(sino)

The same operations can be written with operator syntax when the arrays already reside on the OpenCL device:

device_img = clarray.to_device(queue, img)
sino = R * device_img
backprojection = R.T * sino

Queue selection is deterministic and does not depend on earlier applications. An explicit queue= takes precedence; otherwise the queue is inferred from a device argument or device output. Consequently, every application to a NumPy array must either receive queue= explicitly or receive a caller-provided pyopencl.array.Array output. The multiplication shorthand has no place to pass a queue and is therefore intended for device arrays.

Output reuse and events

Applications allocate a result when output is omitted. Allocation-sensitive iterative code can provide a compatible device output instead:

output = clarray.empty(queue, R.output_shape, dtype=np.float32)
R.apply_to(device_img, output=output)

Compositions forward output to their final operation, while sums write their first summand directly into it before accumulating the remaining terms. Composite intermediates may still be allocated internally. Reusing outputs is recommended in long iterative loops; retaining every newly returned output necessarily retains the corresponding device memory.

OpenCL execution is asynchronous. Passing return_event=True returns (result, events) with the result’s current event list in addition to recording those events on the result array.

Detailed geometry example

The operator API also supports explicit geometry helper objects from gratopy.utilities. This is often the clearest way to specify image extent, detector geometry, shifts, and angular sampling explicitly.

import numpy as np
import pyopencl as cl
import gratopy
from gratopy.utilities import Angles, Detectors, ImageDomain

ctx = cl.create_some_context(interactive=False)
queue = cl.CommandQueue(ctx)

img = np.zeros((192, 128), dtype=np.float32)

image_domain = ImageDomain(
    size=(192, 128),
    extent=3.0,
    center=(0.1, -0.2),
)

angles = Angles.uniform_interval(
    start=0.0,
    end=np.pi / 2,
    number=120,
)

detectors = Detectors(
    number=220,
    extent=3.0,
    center=0.15,
    reversed=False,
)

R = gratopy.operator.Radon(
    image_domain=image_domain,
    angles=angles,
    detectors=detectors,
)

sino = R.apply_to(img, queue=queue)
backproj = R.T.apply_to(sino)

The same setup can also be written directly inline when constructing the operator:

R = gratopy.operator.Radon(
    image_domain=ImageDomain(size=(192, 128), extent=3.0, center=(0.1, -0.2)),
    angles=Angles.uniform_interval(0.0, np.pi / 2, 120),
    detectors=Detectors(number=220, extent=3.0, center=0.15),
)

This explicit style is particularly useful when experimenting with geometry in Python code, because image domain, angles, and detector settings become immutable first-class values that can be safely reused by multiple operators. To change a detector or image setting, construct a new value, for example with dataclasses.replace(). Angles makes private copies of its input arrays and exposes them read-only so subsequent changes to caller-owned arrays cannot invalidate an operator’s cached geometry.

Extent placeholders

For Radon operators, one physical extent can be inferred from the other by using gratopy.utilities.ExtentPlaceholder. This is useful when one wants either the smallest detector covering a fixed image domain, or the largest image domain covered by a fixed detector.

For example, the detector extent can be inferred from a fixed image extent:

from gratopy.utilities import Detectors, ExtentPlaceholder, ImageDomain

R = gratopy.operator.Radon(
    image_domain=ImageDomain(size=128, extent=2.0),
    angles=180,
    detectors=Detectors(number=200, extent=ExtentPlaceholder.FULL),
)

Conversely, the image extent can be inferred from a fixed detector extent:

R = gratopy.operator.Radon(
    image_domain=ImageDomain(size=128, extent=ExtentPlaceholder.FULL),
    angles=180,
    detectors=Detectors(number=200, extent=2.0),
)

Only one side may use an extent placeholder at a time. Passing placeholders for both the image and detector extents is unsupported and raises NotImplementedError. If the requested placeholder semantics are geometrically impossible for the supplied centers and fixed extent, construction raises ValueError.

Adjoint convention

The Radon adjoint uses the same weighted discretization as the legacy API. Angular quadrature weights from gratopy.utilities.Angles are included in the backprojection kernel. Thus R.T denotes the adjoint with respect to gratopy’s physical image and sinogram pairings; it is not generally the plain Euclidean transpose of the unweighted forward-projection matrix.

Operator algebra

Operators inherit from gratopy.operator.base.Operator, which builds non-mutating expression nodes for arithmetic, scaling, adjoints, and composition. Expression nodes retain references to their original operands; they do not copy concrete operators or their runtime caches. For example, one can form a Gram operator

G = R.T * R

and apply it to an image:

gram_img = G.apply_to(img, queue=queue)

This is one of the main motivations for the operator interface: projection operators can be combined with a syntax that mirrors the underlying linear algebra. The multiplication syntax is intentionally overloaded:

• A * B composes two operators,

• alpha * A and A * alpha scale an operator,

• A * x applies an operator to a non-operator argument.

Addition and subtraction construct sum expressions. Algebra creates dedicated expression nodes and never mutates or copies concrete leaves.

Norm estimation

Every operator provides gratopy.operator.base.Operator.norm_estimate(). The default "poweriteration" algorithm applies power iteration to A.T * A and supports OpenCL operators via an explicit queue:

estimate = R.norm_estimate(queue=queue, number_iterations=30)

The alternative "naive" algorithm combines leaf values structurally using the triangle inequality for sums and submultiplicativity of operator norms for compositions. These inequalities preserve certified upper bounds, but unknown leaf norms are currently obtained from finite power iterations and are not themselves certified upper bounds. The resulting combined value is therefore a heuristic unless certified bounds are available for every leaf.

Class structure

The operator implementation is intentionally layered in a small number of classes.

As an experimental interface, the operator API may still change in backward-incompatible ways without a full deprecation cycle while the design is settling.

gratopy.operator.base.Operator

Provides the common operator interface and constructs dedicated adjoint, scale, sum, and composition expression nodes. Concrete operands are shared across the expression tree rather than copied.

gratopy.operator.opencl._OpenCLOperator

Internal helper base for OpenCL-backed operators. It implements shared execution plumbing such as queue inference, array coercion, output allocation, program cacheing, and kernel lookup.

gratopy.operator.projection.Radon

Concrete Radon transform operator. It owns its operator state, geometry-specific preparation, and binds these to the OpenCL kernels.

This separation keeps the generic algebra in gratopy.operator.base backend-agnostic while concentrating OpenCL-specific behavior in a gratopy- specific internal layer.

Compiled-program lifecycle

Compiled OpenCL programs are shared across concrete operators using the same context, kernel source, build options, and template-expansion mode. The global registry retains these program bundles weakly; each operator that has actually used a bundle retains a strong runtime lease. Consequently, equivalent live operators share compilation work, while a bundle is released automatically after its last operator lease disappears. Adjoint and composite expressions retain their concrete leaves and therefore participate in the same lifecycle.

Kernel instances are local to each calling thread because OpenCL kernel arguments are mutable. Threads share the compiled program but not argument state. This protects kernel argument setup, but it does not yet constitute a guarantee that complete operators can be applied concurrently: other lazy runtime caches still require a dedicated thread-safety pass.

The registry can be invalidated explicitly when required:

from gratopy.operator import invalidate_kernel_cache

invalidate_kernel_cache()

Live operators acquire newly compiled bundles on their next application. Invocations already in progress may finish with their existing program. Invalidation removes bundles from future lookup; an existing operator may keep its old lease alive until its next application or until the operator is released.

Custom kernels

One goal of the new operator API is to make kernel experimentation easier. Kernel sources can be specified via gratopy.operator.opencl.OpenCLKernelSpec.

A custom kernel spec can be passed directly to the operator:

from gratopy.operator import OpenCLKernelSpec

spec = OpenCLKernelSpec.from_path("scratch/my_radon.cl", base_name="radon")
R = gratopy.operator.Radon(image_domain=128, angles=180, kernel_spec=spec)

This allows experimenting with alternative kernels while keeping the Python-side operator interface unchanged.

Subclassing _OpenCLOperator

For experimental custom operators, the internal class gratopy.operator.opencl._OpenCLOperator provides a default apply_to() implementation. Although _OpenCLOperator is internal, its documented hook methods form the intended customization surface for OpenCL-backed operators.

The default execution pipeline performs, in order:

1. queue inference,

2. coercion of array-like inputs to pyopencl.array.Array,

3. direction-aware input validation,

4. output allocation (if needed),

5. output validation,

6. forward or adjoint kernel lookup,

7. kernel invocation.

Scalar multiplication is represented by a dedicated expression node instead of mutating or copying the concrete OpenCL operator.

Subclasses can adapt this behavior mostly via hooks instead of overriding the entire method.

Important hooks are:

• gratopy.operator.opencl._OpenCLOperator._default_kernel_spec() for the default kernel bundle,

• gratopy.operator.opencl._OpenCLOperator._expected_output_shape() for shape inference,

• gratopy.operator.opencl._OpenCLOperator._validate_argument() and gratopy.operator.opencl._OpenCLOperator._validate_output() for validation,

• gratopy.operator.opencl._OpenCLOperator._get_kernel() for choosing the compiled kernel,

• gratopy.operator.opencl._OpenCLOperator._kernel_arguments() for supplying additional kernel arguments beyond output and input buffers,

• gratopy.operator.opencl._OpenCLOperator._global_shape() for customizing the OpenCL launch shape.

In simple cases, a custom operator only needs to provide a kernel spec and static input/output shapes. The shared OpenCL implementation dispatches the forward and adjoint kernels through apply_to() and apply_adjoint_to(); T is an adjoint expression wrapper around the same concrete operator and therefore shares its runtime caches.

Limitations and status

The operator API is still evolving. In particular:

• only gratopy.operator.projection.Radon currently has an executable operator implementation; Fanbeam remains a placeholder,

• extent placeholders are currently supported experimentally for Radon operators when exactly one of the image or detector extents is a placeholder,

• higher-level solver interfaces are still centered around the legacy API.

For the full and mature feature set of gratopy, the legacy API documented in Getting started and Reference manual remains the main reference.