Nonlinear Operators

Elementwise Math 🧮

add_elementwise_math_operator(ctx, lhs, rhs, out, opts)

Perform elementwise discrete math operations including floor, fractional part, modulo, step functions, comparisons, and conditional selection. Essential for creating masks, quantization, and control flow in field processing.

Method Signature

add_elementwise_math_operator(ctx, lhs, [rhs], output, [options]) -> operator

Returns: Operator handle (userdata)

Note: The rhs field is optional and can be nil when using rhs_source = "constant".

Mathematical Formulation

Rounding Operations:

• floor: $\text{out} = \lfloor x \rfloor$
• fract: $\text{out} = x - \lfloor x \rfloor \in [0, 1)$

Modulo:

$$\text{out} = x - y \cdot \lfloor x / y \rfloor$$

with guard against division by values smaller than epsilon.

Step Function:

$$\text{out} = \begin{cases} \text{true\_value} & \text{if } x \geq \text{threshold} \\ \text{false\_value} & \text{otherwise} \end{cases}$$

Comparisons:

• eq: $\text{out} = (|x - y| \leq \epsilon) \;?\; \text{true\_value} : \text{false\_value}$
• lt: $\text{out} = (x < y) \;?\; \text{true\_value} : \text{false\_value}$
• gt: $\text{out} = (x > y) \;?\; \text{true\_value} : \text{false\_value}$

Conditional Select:

$$\text{out} = \begin{cases} y & \text{if } x \geq \text{threshold} \\ \text{false\_value} & \text{otherwise} \end{cases}$$

Parameters

Core Parameters:

Parameter	Type	Default	Description
mode	enum	"floor"	Operation to perform (see below)
rhs_source	enum	"constant"	Source for rhs: field or constant
rhs_constant	double	1.0	Constant rhs value when rhs_source = "constant"
accumulate	boolean	false	Add to output instead of overwriting
scale_by_dt	boolean	true	Scale accumulated writes by dt

Mode options: floor, fract, mod, step, eq, lt, gt, select

Comparison & Selection Parameters:

Parameter	Type	Default	Range	Description
threshold	double	0.0	unbounded	Threshold for step/select comparisons
epsilon	double	1e-6	≥0	Tolerance for eq comparisons and mod guards
true_value	double	1.0	unbounded	Value emitted when condition is true
false_value	double	0.0	unbounded	Value emitted when condition is false

Boundary & Initial Conditions

• Operates elementwise; no boundary handling required
• For binary operations with rhs_source = "field", both fields must have matching dimensions
• Complex fields: operations apply to real part only (use Complex Math for full complex operations)

Stability & Convergence

• All operations are unconditionally stable
• floor/fract: Discontinuous at integers; may cause aliasing in smooth fields
• mod: Protected against division by zero via epsilon guard
• step: Discontinuous at threshold; consider smoothstep for continuous alternatives
• Comparisons produce discrete 0/1-like outputs; consider smoothing for gradual transitions

Performance Notes

• Lightweight pointwise operations; O(N) complexity
• Floor/fract use fast hardware intrinsics when available
• Comparisons with rhs_source = "constant" avoid second field read
• Select mode is branch-free on most architectures

Examples

-- Floor quantization
ooc.add_elementwise_math_operator(ctx, continuous, nil, quantized, {
    mode = "floor"
})

-- Extract fractional part (sawtooth from ramp)
ooc.add_elementwise_math_operator(ctx, phase, nil, sawtooth, {
    mode = "fract"
})

-- Modulo 3 wrapping
ooc.add_elementwise_math_operator(ctx, lhs, nil, out, {
    mode = "mod",
    rhs_constant = 3.0,
    rhs_source = "constant"
})

-- Modulo with another field
ooc.add_elementwise_math_operator(ctx, lhs, divisor, out, {
    mode = "mod",
    rhs_source = "field"
})

-- Step function (threshold at 0.5)
ooc.add_elementwise_math_operator(ctx, signal, nil, binary, {
    mode = "step",
    threshold = 0.5,
    true_value = 1.0,
    false_value = 0.0
})

-- Field-vs-field comparison (greater than)
ooc.add_elementwise_math_operator(ctx, a, b, mask, {
    mode = "gt",
    rhs_source = "field"
})

-- Equality test with tolerance
ooc.add_elementwise_math_operator(ctx, measured, expected, match, {
    mode = "eq",
    rhs_source = "field",
    epsilon = 0.01
})

-- Conditional select (choose values where mask >= 0.5)
ooc.add_elementwise_math_operator(ctx, mask, values, out, {
    mode = "select",
    threshold = 0.5,
    rhs_source = "field",
    false_value = 0.0
})

-- Create binary mask from continuous field
ooc.add_elementwise_math_operator(ctx, density, nil, mask, {
    mode = "step",
    threshold = 0.1,
    true_value = 1.0,
    false_value = 0.0
})

---

Complex Math 🔢

add_complex_math_operator(ctx, lhs, rhs, out, opts)

Perform elementwise complex arithmetic, transcendental functions, and component extraction. Supports full complex number operations with configurable scaling, bias, and phase wrapping.

Method Signature

add_complex_math_operator(ctx, lhs, [rhs], output, [options]) -> operator

Returns: Operator handle (userdata)

Note: The rhs field is optional for unary operations or when using rhs_source = "constant".

Mathematical Formulation

Binary Arithmetic:

For complex numbers $z_1 = a + bi$ and $z_2 = c + di$:

• add: $z_1 + z_2 = (a+c) + (b+d)i$
• sub: $z_1 - z_2 = (a-c) + (b-d)i$
• mul: $z_1 \cdot z_2 = (ac - bd) + (ad + bc)i$
• div: $z_1 / z_2 = \frac{ac + bd}{c^2 + d^2} + \frac{bc - ad}{c^2 + d^2}i$
• pow: $z_1^{z_2} = e^{z_2 \ln(z_1)}$

Unary Transcendental:

• exp: $e^z = e^a(\cos b + i \sin b)$
• log: $\ln z = \ln|z| + i \arg(z)$
• sqrt: $\sqrt{z} = \sqrt{|z|} e^{i \arg(z)/2}$

Trigonometric:

• sin: $\sin z = \sin a \cosh b + i \cos a \sinh b$
• cos: $\cos z = \cos a \cosh b - i \sin a \sinh b$
• tan: $\tan z = \sin z / \cos z$
• sinh/cosh/tanh: Hyperbolic analogs

Component Extraction:

• conj: $\bar{z} = a - bi$
• abs: $|z| = \sqrt{a^2 + b^2}$
• arg: $\arg(z) = \text{atan2}(b, a)$
• real: $\text{Re}(z) = a$
• imag: $\text{Im}(z) = b$
• neg: $-z = -a - bi$

Processing Chain:

$$\text{out} = \text{op}(\text{lhs\_scale} \cdot z_1, \text{rhs\_scale} \cdot z_2) + \text{bias}$$

Parameters

Core Parameters:

Parameter	Type	Default	Description
mode	enum	"add"	Operation to perform (see below)
rhs_source	enum	"constant"	Source for rhs: field or constant
accumulate	boolean	false	Add to output instead of overwriting
scale_by_dt	boolean	true	Scale accumulated writes by dt

Mode options:

• Binary: add, sub, mul, div, pow
• Unary transcendental: exp, log, sqrt
• Trigonometric: sin, cos, tan, sinh, cosh, tanh
• Component: conj, abs, arg, real, imag, neg

Constant RHS:

Parameter	Type	Default	Description
rhs_constant_re	double	0.0	Real part of constant rhs
rhs_constant_im	double	0.0	Imaginary part of constant rhs

Scaling & Bias:

Parameter	Type	Default	Description
lhs_scale_re	double	1.0	Real scale applied to lhs
lhs_scale_im	double	0.0	Imaginary scale applied to lhs
rhs_scale_re	double	1.0	Real scale applied to rhs
rhs_scale_im	double	0.0	Imaginary scale applied to rhs
bias_re	double	0.0	Real bias added after operation
bias_im	double	0.0	Imaginary bias added after operation

Numerical Guards:

Parameter	Type	Default	Range	Description
epsilon	double	1e-9	≥0	Tolerance for log/division guards

Output Configuration:

Parameter	Type	Default	Description
output_component	enum	"real"	Component written to real outputs: real, imag, magnitude, phase
phase_wrap	enum	"none"	Phase wrapping: none, signed, unsigned, unit

Phase wrap options:

• none: No wrapping (raw atan2 output)
• signed: Wrap to $[-\pi, \pi]$
• unsigned: Wrap to $[0, 2\pi]$
• unit: Normalize to $[0, 1]$

Boundary & Initial Conditions

• Operates elementwise; no boundary handling required
• For binary operations with rhs_source = "field", both fields must have matching dimensions
• Division and log operations are protected against zero by epsilon

Stability & Convergence

• All operations are numerically stable within floating-point precision
• Division: Protected against division by zero; denominator clamped to epsilon
• Log: Protected against log(0); magnitude clamped to epsilon
• Pow: May produce large values for complex exponents; consider clamping outputs
• Exp: Can overflow for large real parts; typical float range is $e^{88}$

Performance Notes

• Unary operations require single field read
• Binary operations with constant rhs avoid second field read
• Trigonometric and transcendental functions use hardware intrinsics when available
• Component extraction (abs, arg, real, imag) is highly optimized
• Consider using phase_wrap = "unit" for normalized phase comparisons

Examples

-- Complex exponential
ooc.add_complex_math_operator(ctx, z, nil, out, {
    mode = "exp"
})

-- Complex multiplication with constant (0.5 + 0.5i)
ooc.add_complex_math_operator(ctx, z, nil, out, {
    mode = "mul",
    rhs_source = "constant",
    rhs_constant_re = 0.5,
    rhs_constant_im = 0.5
})

-- Field-vs-field multiplication
ooc.add_complex_math_operator(ctx, a, b, product, {
    mode = "mul",
    rhs_source = "field"
})

-- Extract magnitude
ooc.add_complex_math_operator(ctx, z, nil, mag, {
    mode = "abs"
})

-- Extract phase normalized to [0, 1]
ooc.add_complex_math_operator(ctx, z, nil, phase, {
    mode = "arg",
    phase_wrap = "unit"
})

-- Extract phase in [-π, π]
ooc.add_complex_math_operator(ctx, z, nil, phase, {
    mode = "arg",
    phase_wrap = "signed"
})

-- Complex conjugate
ooc.add_complex_math_operator(ctx, z, nil, z_conj, {
    mode = "conj"
})

-- Field-vs-field division with epsilon guard
ooc.add_complex_math_operator(ctx, numerator, denominator, quotient, {
    mode = "div",
    rhs_source = "field",
    epsilon = 1e-6
})

-- Complex logarithm
ooc.add_complex_math_operator(ctx, z, nil, log_z, {
    mode = "log",
    epsilon = 1e-12
})

-- Complex power: z^(0.5) (square root alternative)
ooc.add_complex_math_operator(ctx, z, nil, sqrt_z, {
    mode = "pow",
    rhs_source = "constant",
    rhs_constant_re = 0.5,
    rhs_constant_im = 0.0
})

-- Scale and bias: (2+i)*z + (0.1 + 0.2i)
ooc.add_complex_math_operator(ctx, z, nil, out, {
    mode = "mul",
    rhs_source = "constant",
    rhs_constant_re = 1.0,
    rhs_constant_im = 0.0,
    lhs_scale_re = 2.0,
    lhs_scale_im = 1.0,
    bias_re = 0.1,
    bias_im = 0.2
})

---

Chaos Map 🌀

add_chaos_map_operator(ctx, input, out, opts)

Iterate discrete chaotic maps on a real or complex state field. Complex fields use the real and imaginary components as the two phase-space coordinates; real fields seed the imaginary coordinate at zero and write back the updated real component.

Method Signature

add_chaos_map_operator(ctx, input, output, [options]) -> operator

Returns: Operator handle (userdata)

Mathematical Formulation

At each application, the operator advances the selected map family

$$z_{n+1} = F(z_n; \theta)$$

and optionally blends the result with the previous state:

$$z_{\text{out}} = (1 - \beta) z_n + \beta z_{n+1}$$

where $\beta$ is blend. iterations_per_step applies the selected map multiple times per simulation step before the final blend.

Standard (Chirikov) Map:

$$\begin{aligned} p_{n+1} &= p_n + K \sin(s x_n) \\ x_{n+1} &= x_n + p_{n+1} \end{aligned}$$

where $K$ is the chaos parameter and $s$ is angle_scale. kick_mode chooses whether the map is applied as kick-drift, drift-kick, or a symmetric half-kick split.

Ikeda Map:

$$t = a - \frac{b}{1 + x_n^2 + y_n^2}, \qquad z_{n+1} = o + u z_n e^{i t}$$

where $o = o_x + i o_y$ is the complex offset, $u$ is the contraction factor, and $a$, $b$ control the nonlinear phase shift.

Exponential Map:

$$z_{n+1} = e^{s z_n} + c$$

with complex scale $s$ and offset $c$.

Quadratic Map:

$$z_{n+1} = a z_n^2 + b z_n + c$$

Henon Map:

$$\begin{aligned} x_{n+1} &= o_x + g_x x_n + g_y y_n - a x_n^2 \\ y_{n+1} &= o_y + b x_n \end{aligned}$$

Lozi Map:

$$\begin{aligned} x_{n+1} &= o_x + g_x x_n + g_y y_n - a |x_n| \\ y_{n+1} &= o_y + b x_n \end{aligned}$$

When lozi_abs_epsilon > 0, the absolute value is softened to $\sqrt{x_n^2 + \varepsilon^2}$.

Tinkerbell Map:

$$\begin{aligned} x_{n+1} &= o_x + g_{x^2}x_n^2 + g_{y^2}y_n^2 + a x_n + b y_n \\ y_{n+1} &= o_y + g_{xy}x_n y_n + c x_n + d y_n \end{aligned}$$

When present, k_field, u_field, a_field, b_field, c_field, and d_field override the corresponding constants per sample. wrap_* constrains the real and imaginary coordinates after each update, and escape_* defines how runaway states are handled once |z| exceeds escape_radius.

Parameters

Core Parameters:

Parameter	Type	Default	Range	Description
map_type	enum	"standard"	see below	Chaotic map family
iterations_per_step	integer	1	≥1	Map iterations per simulation step
blend	double	1.0	[0, 1]	Blend factor (1 = full update)

Map type options: standard, ikeda, exponential, quadratic, henon, lozi, tinkerbell

Standard Map Parameters:

Parameter	Type	Default	Range	Description
kick_mode	enum	"kick_drift"	see below	Kick/drift ordering
k	double	1.0	unbounded	Chaos parameter $K$
angle_scale	double	1.0	unbounded	Scale factor $s$ for the sine argument

Kick mode options: kick_drift, drift_kick, symmetric

Ikeda Map Parameters:

Parameter	Type	Default	Range	Description
ikeda_u	double	0.9	[0, 1]	Contraction factor
ikeda_a	double	0.4	unbounded	Phase bias
ikeda_b	double	6.0	unbounded	Nonlinearity strength
ikeda_offset_re	double	1.0	unbounded	Real part of offset $A$
ikeda_offset_im	double	0.0	unbounded	Imaginary part of offset $A$

Exponential Map Parameters:

Parameter	Type	Default	Range	Description
exp_scale_re	double	1.0	unbounded	Real part of scale $s$
exp_scale_im	double	0.0	unbounded	Imaginary part of scale $s$
exp_c_re	double	0.0	unbounded	Real part of constant $c$
exp_c_im	double	0.0	unbounded	Imaginary part of constant $c$

Quadratic Map Parameters:

Parameter	Type	Default	Range	Description
quad_a_re	double	1.0	unbounded	Real part of quadratic coefficient $a$
quad_a_im	double	0.0	unbounded	Imaginary part of quadratic coefficient $a$
quad_b_re	double	0.0	unbounded	Real part of linear coefficient $b$
quad_b_im	double	0.0	unbounded	Imaginary part of linear coefficient $b$
quad_c_re	double	0.0	unbounded	Real part of constant $c$
quad_c_im	double	0.0	unbounded	Imaginary part of constant $c$

Henon Map Parameters:

Parameter	Type	Default	Range	Description
henon_a	double	1.4	unbounded	Quadratic coefficient $a$
henon_b	double	0.3	unbounded	Coupling coefficient $b$
henon_x_gain	double	0.0	unbounded	Linear X gain $g_x$
henon_y_gain	double	1.0	unbounded	Cross-coupling gain $g_y$
henon_offset_re	double	1.0	unbounded	X offset $o_x$
henon_offset_im	double	0.0	unbounded	Y offset $o_y$

Lozi Map Parameters:

Parameter	Type	Default	Range	Description
lozi_a	double	1.7	unbounded	Absolute-value coefficient $a$
lozi_b	double	0.5	unbounded	Coupling coefficient $b$
lozi_x_gain	double	0.0	unbounded	Linear X gain $g_x$
lozi_y_gain	double	1.0	unbounded	Cross-coupling gain $g_y$
lozi_offset_re	double	1.0	unbounded	X offset $o_x$
lozi_offset_im	double	0.0	unbounded	Y offset $o_y$
lozi_abs_epsilon	double	0.0	≥0	Smooth-absolute epsilon (0 uses exact abs)

Tinkerbell Map Parameters:

Parameter	Type	Default	Range	Description
tinkerbell_a	double	0.9	unbounded	Linear X coefficient $a$
tinkerbell_b	double	-0.6013	unbounded	Linear Y coefficient $b$ in the X update
tinkerbell_c	double	2.0	unbounded	Linear X coefficient $c$ in the Y update
tinkerbell_d	double	0.5	unbounded	Linear Y coefficient $d$ in the Y update
tinkerbell_x2_gain	double	1.0	unbounded	Quadratic X gain $g_{x^2}$
tinkerbell_y2_gain	double	-1.0	unbounded	Quadratic Y gain $g_{y^2}$
tinkerbell_xy_gain	double	2.0	unbounded	Cross term gain $g_{xy}$
tinkerbell_offset_re	double	0.0	unbounded	X offset $o_x$
tinkerbell_offset_im	double	0.0	unbounded	Y offset $o_y$

Per-Element Override Fields:

Parameter	Type	Default	Range	Description
k_field	field handle	nil	—	Optional per-sample override for k on the standard map
u_field	field handle	nil	—	Optional per-sample override for ikeda_u
a_field	field handle	nil	—	Optional per-sample override for the map-specific a coefficient
b_field	field handle	nil	—	Optional per-sample override for the map-specific b coefficient
c_field	field handle	nil	—	Optional per-sample override for the map-specific c coefficient
d_field	field handle	nil	—	Optional per-sample override for the map-specific d coefficient

Domain & Escape Controls:

Parameter	Type	Default	Range	Description
wrap_mode_re	enum	"none"	see below	Value-space wrapping mode for the real component
wrap_mode_im	enum	"none"	see below	Value-space wrapping mode for the imaginary component
wrap_min_re	double	0.0	unbounded	Lower real-component wrap bound
wrap_max_re	double	2π	unbounded	Upper real-component wrap bound
wrap_min_im	double	0.0	unbounded	Lower imaginary-component wrap bound
wrap_max_im	double	2π	unbounded	Upper imaginary-component wrap bound
escape_mode	enum	"none"	see below	Action to take when `
escape_radius	double	0.0	≥0	Escape threshold (<= 0 disables escape handling)
escape_reset_re	double	0.0	unbounded	Real component used by escape_mode = "reset"
escape_reset_im	double	0.0	unbounded	Imaginary component used by escape_mode = "reset"

Wrap mode options: none, periodic, clamp, mirror

Escape mode options: none, clamp, reset, nan

Boundary & Initial Conditions

• Operates elementwise; no spatial boundary handling is involved
• Initial conditions determine the trajectory; nearby states can diverge rapidly
• Real input fields evolve on the real axis (Im(z) = 0) and write back only the updated real component
• wrap_* acts on state coordinates, not on spatial sampling
• escape_radius <= 0 disables escape handling entirely

Stability & Convergence

• Standard map: Larger k and higher iterations_per_step generally increase chaotic mixing
• Ikeda map: |ikeda_u| < 1 promotes bounded attractors; larger ikeda_b increases nonlinear folding
• Quadratic, exponential, and tinkerbell maps: Can diverge quickly; wrap_* and escape_* are useful safety valves
• Lozi: lozi_abs_epsilon > 0 smooths the non-differentiable corner at the origin
• blend < 1 damps abrupt state changes but also changes the raw map dynamics

Performance Notes

• iterations_per_step > 1 scales cost linearly with the number of iterations
• Standard, Henon, and Lozi are comparatively cheap polynomial/elementary maps
• Ikeda and exponential families require complex exponentials and are the most expensive per iteration
• Per-element override fields add extra field reads
• wrap_* and escape_* introduce lightweight branching but can prevent runaway values from poisoning later stages

Examples

-- Standard map with K = 1 (near chaotic threshold)
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "standard",
    k = 1.0,
    kick_mode = "symmetric"
})

-- Standard map with strong chaos
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "standard",
    k = 5.0,
    kick_mode = "kick_drift"
})

-- Standard map with regular (integrable) motion
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "standard",
    k = 0.5,  -- below critical value
    angle_scale = 1.0
})

-- Ikeda map (classic parameters)
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "ikeda",
    ikeda_u = 0.9,
    ikeda_a = 0.4,
    ikeda_b = 6.0
})

-- Ikeda map with offset
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "ikeda",
    ikeda_u = 0.85,
    ikeda_a = 0.4,
    ikeda_b = 6.0,
    ikeda_offset_re = 1.0,
    ikeda_offset_im = 0.5
})

-- Quadratic complex map
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "quadratic",
    quad_a_re = 1.0,
    quad_c_re = -0.8,
    quad_c_im = 0.156
})

-- Exponential map (escape-clamped)
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "exponential",
    exp_c_re = -0.8,
    exp_c_im = 0.156,
    escape_mode = "clamp",
    escape_radius = 8.0
})

-- Henon map with explicit linear gains
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "henon",
    henon_a = 1.4,
    henon_b = 0.3,
    henon_y_gain = 1.0
})

-- Lozi map with smoothed absolute value near the origin
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "lozi",
    lozi_a = 1.7,
    lozi_b = 0.5,
    lozi_abs_epsilon = 0.01
})

-- Tinkerbell attractor with softened updates
ooc.add_chaos_map_operator(ctx, state, out, {
    map_type = "tinkerbell",
    tinkerbell_a = 0.9,
    tinkerbell_b = -0.6013,
    tinkerbell_c = 2.0,
    tinkerbell_d = 0.5,
    blend = 0.8
})

-- Spatially varying standard-map chaos parameter
ooc.add_chaos_map_operator(ctx, phase_space, next_phase_space, {
    map_type = "standard",
    k_field = local_k,
    angle_scale = 1.0,
    iterations_per_step = 2
})

---

Hysteretic 🔁

add_hysteretic_operator(ctx, input, output, opts)

Stateful hysteresis operator supporting Schmitt trigger, play model, and Bouc-Wen dynamics. Models memory-dependent nonlinear transformations common in mechanical systems, electronic circuits, and material science. Complex fields are supported and are processed component-wise with independent state on the real and imaginary channels.

Method Signature

add_hysteretic_operator(ctx, input, output, [options]) -> operator

Returns: Operator handle (userdata)

Mathematical Formulation

Schmitt Trigger Mode:

Binary output with hysteresis band:

$$y_n = \begin{cases} y_{\text{high}} & \text{if } x_n > T_{\text{high}} \\ y_{\text{low}} & \text{if } x_n < T_{\text{low}} \\ y_{n-1} & \text{otherwise} \end{cases}$$

where $T_{\text{low}}$ and $T_{\text{high}}$ define the hysteresis thresholds.

Play Model:

Deadband hysteresis where output only changes when input exceeds accumulated play:

$$y_n = \begin{cases} x_n - r & \text{if } x_n > y_{n-1} + r \\ x_n + r & \text{if } x_n < y_{n-1} - r \\ y_{n-1} & \text{otherwise} \end{cases}$$

where $r$ is the play_radius.

Bouc-Wen Model:

Differential hysteresis model for structural mechanics:

$$\dot{z} = A\dot{x} - \beta|\dot{x}||z|^{n-1}z - \gamma\dot{x}|z|^n$$ $$y = \alpha \cdot x + (1-\alpha) \cdot z$$

where $z$ is the internal hysteretic variable and $\alpha$, $A$, $\beta$, $\gamma$, $n$ are shape parameters. The input is first preprocessed by input_gain, input_bias, and input_mode; complex fields are processed component-wise.

For a complex sample $u = u_r + i u_i$, the operator evaluates the same hysteresis law independently on $u_r$ and $u_i$, then recombines the result:

$$H(u) = H(u_r) + i H(u_i)$$

Parameters

Core Parameters:

Parameter	Type	Default	Description
mode	enum	"schmitt"	Hysteresis model: schmitt, play, bouc_wen
threshold_mode	enum	"bounds"	Threshold specification: bounds or center_width
input_mode	enum	"direct"	Input preprocessing: direct, abs, squared
input_gain	double	1.0	Gain applied to input before processing
input_bias	double	0.0	Bias added to input
output_gain	double	1.0	Gain applied to output
output_bias	double	0.0	Bias added to output
accumulate	boolean	false	Add to output instead of overwriting
scale_by_dt	boolean	true	Scale accumulated writes by dt

Threshold Parameters (bounds mode):

Parameter	Type	Default	Description
threshold_low	double	-0.5	Lower threshold for hysteresis band
threshold_high	double	0.5	Upper threshold for hysteresis band

Threshold Parameters (center_width mode):

Parameter	Type	Default	Description
threshold_center	double	0.0	Center of hysteresis band
threshold_width	double	1.0	Width of hysteresis band (≥0)

Schmitt Trigger Parameters:

Parameter	Type	Default	Description
output_low	double	-1.0	Output when below lower threshold
output_high	double	1.0	Output when above upper threshold

Play Model Parameters:

Parameter	Type	Default	Description
play_radius	double	0.0	Deadband radius for play hysteresis

Bouc-Wen Parameters:

Parameter	Type	Default	Description
bw_alpha	double	0.1	Ratio of post-yield to pre-yield stiffness
bw_A	double	1.0	Controls tangent stiffness
bw_beta	double	0.5	Shape parameter (softening/hardening)
bw_gamma	double	0.5	Shape parameter (pinching)
bw_n	double	2.0	Exponent controlling sharpness
bw_z_clamp	double	0.0	Clamp on internal variable z (0 = no clamp)
bw_xdot_clamp	double	0.0	Clamp on velocity estimate (0 = no clamp)

Smoothing & Rate Limiting:

Parameter	Type	Default	Description
smooth	double	0.0	Exponential smoothing factor [0, 1]
rate_limit	double	0.0	Maximum
state_min	double	-1e6	Minimum internal state value
state_max	double	1e6	Maximum internal state value

Initialization:

Parameter	Type	Default	Description
initialize_from_input	boolean	true	Initialize state from first input sample
initial_output	double	0.0	Initial output if not from input
initial_input	double	0.0	Initial input reference
initial_z	double	0.0	Initial Bouc-Wen internal variable

Boundary & Initial Conditions

• Operates elementwise; no spatial boundary handling
• Internal state is maintained per-element across timesteps
• Complex fields keep separate hysteresis state for the real and imaginary channels
• Use initialize_from_input = true for automatic state initialization
• For deterministic behavior, set explicit initial values

Stability & Convergence

• Schmitt mode: Unconditionally stable; discrete state transitions
• Play mode: Unconditionally stable; bounded by input range
• Bouc-Wen mode: Stability depends on parameter choices:
  • $\beta + \gamma > 0$ ensures bounded dissipation
  • Large $n$ values can cause stiff dynamics; reduce timestep accordingly
  • Use bw_z_clamp to prevent runaway internal states

Performance Notes

• Maintains per-element internal state; memory usage scales with field size
• Schmitt and play modes are computationally lightweight
• Bouc-Wen requires velocity estimation (finite difference) and is more expensive
• Complex fields approximately double the scalar state/work because the channels are evolved independently
• rate_limit > 0 adds conditional logic per element

Examples

-- Schmitt trigger with ±0.3 thresholds
ooc.add_hysteretic_operator(ctx, input, output, {
    mode = "schmitt",
    threshold_low = -0.3,
    threshold_high = 0.3,
    output_low = 0.0,
    output_high = 1.0
})

-- Center-width threshold specification
ooc.add_hysteretic_operator(ctx, input, output, {
    mode = "schmitt",
    threshold_mode = "center_width",
    threshold_center = 0.5,
    threshold_width = 0.2  -- thresholds at 0.4 and 0.6
})

-- Play model with deadband
ooc.add_hysteretic_operator(ctx, position, force, {
    mode = "play",
    play_radius = 0.1,
    input_gain = 100.0  -- spring constant
})

-- Bouc-Wen for structural hysteresis
ooc.add_hysteretic_operator(ctx, displacement, restoring_force, {
    mode = "bouc_wen",
    bw_alpha = 0.05,
    bw_A = 1.0,
    bw_beta = 0.5,
    bw_gamma = 0.5,
    bw_n = 2.0
})

-- Smooth Schmitt with rate limiting (debouncing)
ooc.add_hysteretic_operator(ctx, noisy_signal, clean_signal, {
    mode = "schmitt",
    threshold_low = -0.1,
    threshold_high = 0.1,
    smooth = 0.9,
    rate_limit = 10.0  -- max 10 units/second change rate
})

-- Magnitude-based hysteresis (rectified input)
ooc.add_hysteretic_operator(ctx, ac_signal, envelope, {
    mode = "play",
    input_mode = "abs",
    play_radius = 0.05,
    smooth = 0.8
})

-- Complex I/Q hysteresis (applied independently to real and imaginary parts)
ooc.add_hysteretic_operator(ctx, iq_signal, iq_shaped, {
    mode = "schmitt",
    threshold_low = -0.2,
    threshold_high = 0.2,
    output_low = -1.0,
    output_high = 1.0
})

---

Neural Infer 🧠

add_neural_infer_operator(ctx, input, out, opts)

Run an inference-only learned operator on a field while carrying explicit metadata about determinism, device placement, precision, and accepted tensor shape.

Experimental: Model adapters, execution contracts, and fallback behavior may evolve as the neural runtime is refined.

Method Signature

add_neural_infer_operator(ctx, input, output, [options]) -> operator

Returns: Operator handle (userdata)

Mathematical Formulation

Conceptually, the operator applies a learned map:

$$\hat{y} = F_{\theta}(\tilde{x})$$

with optional affine input/output transforms

$$\tilde{x} = \begin{cases} s_{\text{in}} x + b_{\text{in}} & \text{when input normalization is enabled} \\ x & \text{otherwise} \end{cases}$$

and writes

$$y = \begin{cases} \Delta t \cdot \hat{y} & \text{when dt scaling is enabled} \\ \hat{y} & \text{otherwise} \end{cases}$$

Here, normalize_input controls the affine preprocessing and scale_by_dt controls the optional $\Delta t$ writeback scaling.

If no backend model or custom inference callback is available, the current runtime falls back to an affine passthrough using the configured normalization/denormalization settings.

Parameters

Inference Contract:

Parameter	Type	Default	Description
model_id	string	"model"	Model identifier resolved by the backend adapter
determinism_policy	enum	"inherit"	inherit, strict, best_effort, off
device_requirement	enum	"any"	any, cpu_only, accelerator_preferred, accelerator_required
precision_mode	enum	"default"	default, fp32, fp64, mixed, fp16, bf16

Shape Constraints:

Parameter	Type	Default	Description
min_rank	integer	1	Minimum accepted input rank
max_rank	integer	0	Maximum accepted input rank (0 disables the upper bound)
channel_axis	integer	255	Channel axis (255 = auto)
channels_last	boolean	true	Auto-axis policy when channel_axis = 255
allow_complex_input	boolean	true	Allow complex-valued input tensors/fields

Normalization & Writeback:

Parameter	Type	Default	Description
normalize_input	boolean	false	Apply affine preprocessing before inference
input_scale	double	1.0	Input normalization scale
input_bias	double	0.0	Input normalization bias
output_scale	double	1.0	Output denormalization scale
output_bias	double	0.0	Output denormalization bias
accumulate	boolean	false	Add prediction into the output field
scale_by_dt	boolean	true	Scale writes by substep dt

Boundary & Initial Conditions

• Input and output fields must have identical shape and storage type
• Rank and channel layout must satisfy the configured shape constraints
• Complex inputs are permitted only when allow_complex_input = true

Stability & Convergence

• Numerical stability depends on the learned model rather than a fixed analytic formula
• strict determinism may reduce backend choices but improves reproducibility
• scale_by_dt = true is appropriate when the model predicts a rate-like increment rather than a full state

Performance Notes

• Allocates a temporary prediction buffer each step
• Backend selection and precision can materially affect throughput
• The operator is nonlocal from the scheduler’s perspective even though the I/O shape matches pointwise field updates

Examples

-- Deterministic surrogate inference
ooc.add_neural_infer_operator(ctx, state, prediction, {
    model_id = "burgers_surrogate_v1",
    determinism_policy = "strict",
    precision_mode = "fp32"
})

-- Explicit tensor layout contract
ooc.add_neural_infer_operator(ctx, input, output, {
    model_id = "vision_block",
    min_rank = 3,
    max_rank = 3,
    channel_axis = 0,
    channels_last = false
})

-- Accumulate model output as an RHS term
ooc.add_neural_infer_operator(ctx, state, rhs, {
    model_id = "residual_term",
    accumulate = true,
    scale_by_dt = true
})

---

Neural Hybrid 🔀

add_neural_hybrid_operator(ctx, input, out, opts)

Combine an analytic baseline with a learned residual prediction. This is useful when you want the model to correct, not replace, a known signal path.

Experimental: Model adapters, execution contracts, and blend semantics may evolve as the neural runtime is refined.

Method Signature

add_neural_hybrid_operator(ctx, input, output, [options]) -> operator

Returns: Operator handle (userdata)

Mathematical Formulation

The operator forms a hybrid write from the input field $x$ and a model prediction $\hat{r}$:

$$y = g_a x + g_r \hat{r}$$

where:

• $g_a$ is analytic_gain
• $g_r$ is residual_gain
• $\hat{r} = F_{\theta}(\tilde{x})$

As with add_neural_infer_operator, the runtime can normalize inputs before inference and apply output scaling/bias afterward. If no model backend is available, the residual path currently falls back to an affine-transformed copy of the input.

Parameters

Hybrid Weights:

Parameter	Type	Default	Description
analytic_gain	double	1.0	Scale applied to the direct analytic/input contribution
residual_gain	double	1.0	Scale applied to the learned residual contribution

Inference Contract:

Parameter	Type	Default	Description
model_id	string	"model"	Residual model identifier
determinism_policy	enum	"inherit"	inherit, strict, best_effort, off
device_requirement	enum	"any"	any, cpu_only, accelerator_preferred, accelerator_required
precision_mode	enum	"default"	default, fp32, fp64, mixed, fp16, bf16

Shape Constraints & Normalization:

Parameter	Type	Default	Description
min_rank	integer	1	Minimum accepted input rank
max_rank	integer	0	Maximum accepted input rank (0 disables the upper bound)
channel_axis	integer	255	Channel axis (255 = auto)
channels_last	boolean	true	Auto-axis policy when channel_axis = 255
allow_complex_input	boolean	true	Allow complex-valued input tensors/fields
normalize_input	boolean	false	Apply affine preprocessing before residual inference
input_scale	double	1.0	Input normalization scale
input_bias	double	0.0	Input normalization bias
output_scale	double	1.0	Output denormalization scale on the residual path
output_bias	double	0.0	Output denormalization bias on the residual path
accumulate	boolean	false	Add the hybrid output into the destination field
scale_by_dt	boolean	true	Scale writes by substep dt

Boundary & Initial Conditions

• Input and output fields must match in shape and scalar/complex storage
• The same shape-contract rules as Neural Infer apply
• Complex outputs combine the analytic and residual paths component-wise

Stability & Convergence

• analytic_gain and residual_gain directly control how strongly the learned correction can perturb the baseline
• Start with small residual_gain values when introducing a new surrogate into an existing simulation
• scale_by_dt = true is appropriate when the hybrid output represents an increment rather than a full-state overwrite

Performance Notes

• Same backend and buffer costs as Neural Infer, plus the final blend with the input field
• Useful for residual-learning workflows where the analytic path remains interpretable
• Determinism, device, and precision choices are exposed explicitly so the runtime contract stays inspectable

Examples

-- Residual correction on top of the current state
ooc.add_neural_hybrid_operator(ctx, state, next_state, {
    model_id = "navier_stokes_residual_v2",
    analytic_gain = 1.0,
    residual_gain = 0.2
})

-- CPU-only deterministic hybrid for reproducible tests
ooc.add_neural_hybrid_operator(ctx, u, out, {
    model_id = "test_residual",
    determinism_policy = "strict",
    device_requirement = "cpu_only",
    residual_gain = 0.1
})

-- Accumulate learned correction into an RHS buffer
ooc.add_neural_hybrid_operator(ctx, state, rhs, {
    model_id = "forcing_residual",
    analytic_gain = 0.0,
    residual_gain = 1.0,
    accumulate = true,
    scale_by_dt = true
})