Nonlinear Operators

🧮 Elementwise Math

sim_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

sim_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.sim_add_elementwise_math_operator(ctx, continuous, nil, quantized, {
    mode = "floor"
})

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

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

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

-- Step function (threshold at 0.5)
ooc.sim_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.sim_add_elementwise_math_operator(ctx, a, b, mask, {
    mode = "gt",
    rhs_source = "field"
})

-- Equality test with tolerance
ooc.sim_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.sim_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.sim_add_elementwise_math_operator(ctx, density, nil, mask, {
    mode = "step",
    threshold = 0.1,
    true_value = 1.0,
    false_value = 0.0
})

---

🔢 Complex Math

sim_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

sim_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.sim_add_complex_math_operator(ctx, z, nil, out, {
    mode = "exp"
})

-- Complex multiplication with constant (0.5 + 0.5i)
ooc.sim_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.sim_add_complex_math_operator(ctx, a, b, product, {
    mode = "mul",
    rhs_source = "field"
})

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

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

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

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

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

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

-- Complex power: z^(0.5) (square root alternative)
ooc.sim_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.sim_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

sim_add_chaos_map_operator(ctx, input, out, opts)

Iterate discrete chaotic maps on a complex state field. The real and imaginary components encode the two-dimensional phase space. Useful for generating deterministic chaos, studying nonlinear dynamics, and creating complex textures.

Method Signature

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

Returns: Operator handle (userdata)

Mathematical Formulation

Standard (Chirikov) Map:

The standard map is a canonical example of Hamiltonian chaos, encoding angle $\theta$ in the real part and momentum $p$ in the imaginary part:

$$\begin{aligned} p_{n+1} &= p_n + K \sin(\alpha \theta_n) \\ \theta_{n+1} &= \theta_n + p_{n+1} \pmod{2\pi} \end{aligned}$$

where $K$ is the chaos parameter and $\alpha$ is angle_scale.

Kick modes:

• kick_drift: Update $p$ first, then $\theta$ (default)
• drift_kick: Update $\theta$ first, then $p$
• symmetric: Half-kick, drift, half-kick (symplectic)

Ikeda Map:

Models light in a nonlinear optical resonator:

$$z_{n+1} = A + u \cdot z_n \cdot e^{i(|z_n|^2 b + a)}$$

where:

• $A$ is the offset (ikeda_offset_re + i * ikeda_offset_im)
• $u$ is the contraction factor
• $a$ is the phase bias
• $b$ is the nonlinearity strength

Exponential Map:

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

where $s$ is the complex scale and $c$ is the complex constant.

Blending:

$$z_{\text{out}} = (1 - \beta) z_{\text{old}} + \beta \cdot z_{\text{new}}$$

where $\beta$ is the blend parameter.

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

Standard Map Parameters:

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

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$

Boundary & Initial Conditions

• Operates elementwise; no spatial boundary handling
• Initial conditions determine the trajectory; small changes can lead to vastly different outcomes (sensitive dependence)
• Standard map: Typically initialize with $\theta \in [0, 2\pi)$ and $p$ near zero
• Ikeda map: Typically initialize near the fixed point or attractor
• Exponential map: Initialize within bounded region to avoid overflow

Stability & Convergence

• Standard map:

  • $K < K_c \approx 0.9716$ produces mostly regular (KAM tori) behavior
  • $K > K_c$ transitions to widespread chaos
  • $K \approx 1$ shows mixed regular/chaotic phase space
  • Large $K$ values produce fully chaotic behavior

• Ikeda map:

  • $u < 1$ ensures contraction (bounded attractor)
  • Multiple attractors possible; initial conditions determine which is reached
  • Larger $b$ increases complexity of the attractor

• Exponential map:

  • Can escape to infinity; consider clamping or careful parameter choice
  • Julia sets exist for certain $c$ values (e.g., $c \approx -0.65$)

Performance Notes

• iterations_per_step > 1 multiplies computation per timestep
• Standard map uses one sin evaluation per iteration
• Ikeda map requires magnitude calculation and complex exponential
• Exponential map requires complex exponential (expensive)
• blend < 1 adds interpolation overhead but smooths transitions

Examples

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

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

-- Standard map with regular (integrable) motion
ooc.sim_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.sim_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.sim_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
})

-- Exponential map (Julia set exploration)
ooc.sim_add_chaos_map_operator(ctx, state, out, {
    map_type = "exponential",
    exp_c_re = -0.8,
    exp_c_im = 0.156
})

-- Exponential map with scaling
ooc.sim_add_chaos_map_operator(ctx, state, out, {
    map_type = "exponential",
    exp_scale_re = 0.5,
    exp_scale_im = 0.1,
    exp_c_re = -0.65,
    exp_c_im = 0.0
})

-- Multiple iterations per step with blending (smooth animation)
ooc.sim_add_chaos_map_operator(ctx, state, out, {
    map_type = "standard",
    k = 0.95,
    iterations_per_step = 4,
    blend = 0.75
})

-- Fast iteration for texture generation
ooc.sim_add_chaos_map_operator(ctx, state, out, {
    map_type = "ikeda",
    ikeda_u = 0.9,
    ikeda_b = 6.0,
    iterations_per_step = 10,
    blend = 1.0
})

---

🔁 Hysteretic

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

Method Signature

sim_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 z + (1-\alpha) \cdot x$$

where $z$ is the internal hysteretic variable and $\alpha$, $A$, $\beta$, $\gamma$, $n$ are shape parameters.

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
• 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
• rate_limit > 0 adds conditional logic per element

Examples

-- Schmitt trigger with ±0.3 thresholds
ooc.sim_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.sim_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.sim_add_hysteretic_operator(ctx, position, force, {
    mode = "play",
    play_radius = 0.1,
    input_gain = 100.0  -- spring constant
})

-- Bouc-Wen for structural hysteresis
ooc.sim_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.sim_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.sim_add_hysteretic_operator(ctx, ac_signal, envelope, {
    mode = "play",
    input_mode = "abs",
    play_radius = 0.05,
    smooth = 0.8
})