pylinkage.solver package

Submodules

pylinkage.solver.acceleration module

Per-joint numba acceleration solvers for kinematic analysis.

Each function computes the linear acceleration of a specific joint type given parent accelerations, velocities, positions, and constraints. All functions are numba-compiled for maximum performance.

Acceleration equations are derived by differentiating velocity equations, including centripetal and Coriolis terms where applicable.

pylinkage.solver.acceleration.solve_crank_acceleration(x: float, y: float, vx: float, vy: float, anchor_x: float, anchor_y: float, anchor_vx: float, anchor_vy: float, anchor_ax: float, anchor_ay: float, radius: float, omega: float, alpha: float) → tuple[float, float]

Compute crank acceleration including centripetal term.

The crank velocity is: v = anchor_v + r*ω*(-sin(θ), cos(θ)) Differentiating: a = anchor_a + r*α*(-sin(θ), cos(θ)) + r*ω²*(-cos(θ), -sin(θ))

The second term is tangential acceleration, third is centripetal.

Parameters:

x – Current X position of the crank.
y – Current Y position of the crank.
vx – Current X velocity of the crank (unused, for consistency).
vy – Current Y velocity of the crank (unused, for consistency).
anchor_x – X position of the anchor.
anchor_y – Y position of the anchor.
anchor_vx – X velocity of the anchor (unused, for consistency).
anchor_vy – Y velocity of the anchor (unused, for consistency).
anchor_ax – X acceleration of the anchor (usually 0).
anchor_ay – Y acceleration of the anchor (usually 0).
radius – Distance from anchor to crank.
omega – Angular velocity in rad/s.
alpha – Angular acceleration in rad/s².

Returns:

Acceleration (ax, ay) of the crank.

pylinkage.solver.acceleration.solve_fixed_acceleration(x: float, y: float, vx: float, vy: float, p0_x: float, p0_y: float, p0_vx: float, p0_vy: float, p0_ax: float, p0_ay: float, p1_x: float, p1_y: float, p1_vx: float, p1_vy: float, p1_ax: float, p1_ay: float, radius: float, angle: float) → tuple[float, float]

Compute fixed joint acceleration by differentiating velocity equations.

The velocity is:: vx = p0_vx - r * (dα/dt) * sin(α + γ) vy = p0_vy + r * (dα/dt) * cos(α + γ)

Differentiating gives acceleration including angular acceleration term.

Parameters:

x – Current X position of the joint.
y – Current Y position of the joint.
vx – Current X velocity of the joint.
vy – Current Y velocity of the joint.
p0_x – Position of first parent (origin).
p0_y – Position of first parent (origin).
p0_vx – Velocity of first parent.
p0_vy – Velocity of first parent.
p0_ax – Acceleration of first parent.
p0_ay – Acceleration of first parent.
p1_x – Position of second parent (angle reference).
p1_y – Position of second parent (angle reference).
p1_vx – Velocity of second parent.
p1_vy – Velocity of second parent.
p1_ax – Acceleration of second parent.
p1_ay – Acceleration of second parent.
radius – Distance from first parent.
angle – Angle offset from the parent-to-parent direction.

Returns:

Acceleration (ax, ay) of the joint, or (NaN, NaN) if singular.

pylinkage.solver.acceleration.solve_prismatic_acceleration(x: float, y: float, vx: float, vy: float, circle_x: float, circle_y: float, circle_vx: float, circle_vy: float, circle_ax: float, circle_ay: float, radius: float, line_p1_x: float, line_p1_y: float, line_p1_vx: float, line_p1_vy: float, line_p1_ax: float, line_p1_ay: float, line_p2_x: float, line_p2_y: float, line_p2_vx: float, line_p2_vy: float, line_p2_ax: float, line_p2_ay: float) → tuple[float, float]

Compute prismatic joint acceleration using implicit differentiation.

Differentiates the velocity constraint equations to get acceleration.

Parameters:

x – Current position of the joint.
y – Current position of the joint.
vx – Current velocity of the joint.
vy – Current velocity of the joint.
circle_x – Position of circle center.
circle_y – Position of circle center.
circle_vx – Velocity of circle center.
circle_vy – Velocity of circle center.
circle_ax – Acceleration of circle center.
circle_ay – Acceleration of circle center.
radius – Circle radius.
line_p1_x – Position of first line point.
line_p1_y – Position of first line point.
line_p1_vx – Velocity of first line point.
line_p1_vy – Velocity of first line point.
line_p1_ax – Acceleration of first line point.
line_p1_ay – Acceleration of first line point.
line_p2_x – Position of second line point.
line_p2_y – Position of second line point.
line_p2_vx – Velocity of second line point.
line_p2_vy – Velocity of second line point.
line_p2_ax – Acceleration of second line point.
line_p2_ay – Acceleration of second line point.

Returns:

Acceleration (ax, ay) of the joint, or (NaN, NaN) if singular.

pylinkage.solver.acceleration.solve_revolute_acceleration(x: float, y: float, vx: float, vy: float, p0_x: float, p0_y: float, p0_vx: float, p0_vy: float, p0_ax: float, p0_ay: float, p1_x: float, p1_y: float, p1_vx: float, p1_vy: float, p1_ax: float, p1_ay: float) → tuple[float, float]

Compute revolute joint acceleration using implicit differentiation.

Starting from the velocity constraint equations and differentiating again with respect to time.

From the constraint (x - p0_x)² + (y - p0_y)² = r0²: Velocity: (x - p0_x)(vx - p0_vx) + (y - p0_y)(vy - p0_vy) = 0 Acceleration (differentiating velocity constraint): (vx - p0_vx)² + (vy - p0_vy)² + (x - p0_x)(ax - p0_ax) + (y - p0_y)(ay - p0_ay) = 0

Rearranging for ax, ay gives a 2x2 linear system.

Parameters:

x – Current X position of the joint.
y – Current Y position of the joint.
vx – Current X velocity of the joint.
vy – Current Y velocity of the joint.
p0_x – Position of first parent.
p0_y – Position of first parent.
p0_vx – Velocity of first parent.
p0_vy – Velocity of first parent.
p0_ax – Acceleration of first parent.
p0_ay – Acceleration of first parent.
p1_x – Position of second parent.
p1_y – Position of second parent.
p1_vx – Velocity of second parent.
p1_vy – Velocity of second parent.
p1_ax – Acceleration of second parent.
p1_ay – Acceleration of second parent.

Returns:

Acceleration (ax, ay) of the joint, or (NaN, NaN) if singular.

pylinkage.solver.acceleration.solve_rigid_body_acceleration(x: float, y: float, p0_x: float, p0_y: float, p0_vx: float, p0_vy: float, p0_ax: float, p0_ay: float, p1_x: float, p1_y: float, p1_vx: float, p1_vy: float, p1_ax: float, p1_ay: float) → tuple[float, float]

Compute the acceleration of a third point on a rigid body.

The acceleration counterpart of solve_rigid_body_velocity(). Two points of a planar rigid body fix both its angular velocity and its angular acceleration:

a1 = a0 + α ẑ × d - ω² d ⟹ α = [dx * (p1_ay - p0_ay) - dy * (p1_ax - p0_ax)] / d²

with d = p1 - p0 (the centripetal term drops out of the cross product). The queried point then follows:

a = a0 + α ẑ × (p - p0) - ω² (p - p0)

Well-conditioned for collinear points, where intersecting the two differentiated distance constraints is singular.

Parameters:

x – Current X position of the queried point.
y – Current Y position of the queried point.
p0_x – Position of the first body point.
p0_y – Position of the first body point.
p0_vx – Velocity of the first body point.
p0_vy – Velocity of the first body point.
p0_ax – Acceleration of the first body point.
p0_ay – Acceleration of the first body point.
p1_x – Position of the second body point.
p1_y – Position of the second body point.
p1_vx – Velocity of the second body point.
p1_vy – Velocity of the second body point.
p1_ax – Acceleration of the second body point.
p1_ay – Acceleration of the second body point.

Returns:

Acceleration (ax, ay) of the queried point, or (NaN, NaN) if the two body points coincide.

pylinkage.solver.groups module

Solver functions for Assur groups.

This module provides standalone solving functions for Assur groups. These functions operate on pure geometric data (positions, distances) and delegate to the numba-optimized joint solvers.

The Assur group classes in pylinkage.assur define structure (logical properties), while this module provides the solving behavior.

pylinkage.solver.groups.solve_pp_dyad(line1_pos1: Coord, line1_pos2: Coord, line2_pos1: Coord, line2_pos2: Coord) → Coord

Solve a PP dyad position using line-line intersection.

A PP dyad consists of one internal node constrained to lie at the intersection of two lines. Each line is defined by two points.

This is used for isomers like T_R_T, T_RT_ and _TRT_.

Parameters:

line1_pos1 – Position (x, y) of the first point on line 1.
line1_pos2 – Position (x, y) of the second point on line 1.
line2_pos1 – Position (x, y) of the first point on line 2.
line2_pos2 – Position (x, y) of the second point on line 2.

Returns:

Computed (x, y) position for the internal node.

Raises:

ValueError – If the lines are parallel (no intersection).

Example

>>> pos = solve_pp_dyad(
...     line1_pos1=(0.0, 0.0),
...     line1_pos2=(4.0, 0.0),
...     line2_pos1=(2.0, -1.0),
...     line2_pos2=(2.0, 3.0),
... )
>>> # pos is (2.0, 0.0) - intersection of x-axis and line x=2

pylinkage.solver.groups.solve_rrp_dyad(anchor_pos: Coord, line_pos1: Coord, line_pos2: Coord, distance: float, hint: Coord | None = None) → Coord

Solve an RRP dyad position using circle-line intersection.

An RRP dyad consists of one internal node connected to an anchor by a link of known length, and constrained to lie on a line defined by two other nodes. The internal node position is at the intersection of the circle and line.

Parameters:

anchor_pos – Position (x, y) of the revolute anchor node.
line_pos1 – Position (x, y) of the first node defining the line.
line_pos2 – Position (x, y) of the second node defining the line.
distance – Distance from anchor to internal node.
hint – Current position for disambiguation when two solutions exist. If None, defaults to anchor_pos.

Returns:

Computed (x, y) position for the internal node.

Raises:

ValueError – If the circle and line don’t intersect (unbuildable).

Example

>>> pos = solve_rrp_dyad(
...     anchor_pos=(0.0, 0.0),
...     line_pos1=(0.0, 2.0),
...     line_pos2=(4.0, 2.0),
...     distance=2.5,
... )
>>> # pos is on the line y=2 at distance 2.5 from origin

pylinkage.solver.groups.solve_rrr_dyad(anchor0_pos: Coord, anchor1_pos: Coord, distance0: float, distance1: float, hint: Coord | None = None) → Coord

Solve an RRR dyad position using circle-circle intersection.

An RRR dyad consists of one internal node connected to two anchor nodes by links of known lengths. The internal node position is at the intersection of two circles centered at the anchors.

Parameters:

anchor0_pos – Position (x, y) of the first anchor node.
anchor1_pos – Position (x, y) of the second anchor node.
distance0 – Distance from anchor0 to internal node.
distance1 – Distance from anchor1 to internal node.
hint – Current position for disambiguation when two solutions exist. If None, defaults to anchor0_pos.

Returns:

Computed (x, y) position for the internal node.

Raises:

ValueError – If the circles don’t intersect (unbuildable configuration).

Example

>>> pos = solve_rrr_dyad(
...     anchor0_pos=(0.0, 0.0),
...     anchor1_pos=(3.0, 0.0),
...     distance0=2.0,
...     distance1=2.0,
... )
>>> # pos is approximately (1.5, 1.32...) or (1.5, -1.32...)

pylinkage.solver.groups.solve_triad(constraints: list[tuple[str, str, float]], positions: dict[str, Coord], internal_ids: tuple[str, str], hints: dict[str, Coord] | None = None) → dict[str, Coord]

Solve a triad (Class II Assur group) using Newton-Raphson.

A triad has 2 internal nodes whose positions are unknown. The constraints are distance equations between pairs of nodes (some known, some unknown). We solve the system of equations:

(x_a - x_b)² + (y_a - y_b)² = d² for each edge

where (x_a, y_a) and (x_b, y_b) are node positions and d is the edge distance. Known node positions are substituted, leaving 4 unknowns (x0, y0, x1, y1) for the 2 internal nodes.

Parameters:

constraints – List of (node_a, node_b, distance) tuples. Each defines a distance constraint between two nodes. At least 4 constraints are needed (4 unknowns).
positions – Known positions of anchor nodes.
internal_ids – The two internal node IDs (order matches unknowns).
hints – Optional initial guess positions for internal nodes. If not provided, centroid of anchors is used.

Returns:

Dict mapping internal node IDs to their computed (x, y) positions.

Raises:

ValueError – If the system cannot be solved (unbuildable).

pylinkage.solver.joints module

Per-joint numba solvers for the simulation loop.

Each function solves a specific joint type given the parent positions and constraints. All functions are numba-compiled for maximum performance.

pylinkage.solver.joints.solve_arc_crank(current_angle: float, direction: float, anchor_x: float, anchor_y: float, radius: float, angle_rate: float, arc_start: float, arc_end: float, dt: float) → tuple[float, float, float, float]

Solve arc crank position using oscillating angular motion.

The arc crank rotates around its anchor point, reversing direction when it reaches the angle limits.

Parameters:

current_angle – Current angle in radians.
direction – Current movement direction (+1.0 or -1.0).
anchor_x – X position of the anchor (rotation center).
anchor_y – Y position of the anchor (rotation center).
radius – Distance from anchor to crank.
angle_rate – Angular velocity magnitude (radians per unit time).
arc_start – Minimum angle limit in radians.
arc_end – Maximum angle limit in radians.
dt – Time step.

Returns:

Tuple of (new_x, new_y, new_angle, new_direction).

pylinkage.solver.joints.solve_crank(current_x: float, current_y: float, anchor_x: float, anchor_y: float, radius: float, angle_rate: float, dt: float) → tuple[float, float]

Solve crank position using polar rotation.

The crank rotates around its anchor point at a constant angular rate.

Parameters:

current_x – Current X position of the crank.
current_y – Current Y position of the crank.
anchor_x – X position of the anchor (rotation center).
anchor_y – Y position of the anchor (rotation center).
radius – Distance from anchor to crank.
angle_rate – Angular velocity (radians per unit time).
dt – Time step.

Returns:

New (x, y) position of the crank.

pylinkage.solver.joints.solve_fixed(p0_x: float, p0_y: float, p1_x: float, p1_y: float, radius: float, angle: float) → tuple[float, float]

Solve fixed joint using polar projection.

The fixed joint is at a fixed distance and angle from parent 0, with the angle measured relative to the line from parent 0 to parent 1.

Parameters:

p0_x – X position of first parent (origin).
p0_y – Y position of first parent (origin).
p1_x – X position of second parent (angle reference).
p1_y – Y position of second parent (angle reference).
radius – Distance from first parent.
angle – Angle offset from the parent-to-parent direction.

Returns:

New (x, y) position (always deterministic).

pylinkage.solver.joints.solve_line_line(line1_p1_x: float, line1_p1_y: float, line1_p2_x: float, line1_p2_y: float, line2_p1_x: float, line2_p1_y: float, line2_p2_x: float, line2_p2_y: float) → tuple[float, float]

Solve line-line joint using line-line intersection.

The joint is positioned at the intersection of two lines: - Line 1: passing through line1_p1 and line1_p2 - Line 2: passing through line2_p1 and line2_p2

This is deterministic when lines are not parallel (single intersection).

Parameters:

line1_p1_x – X position of first point on line 1.
line1_p1_y – Y position of first point on line 1.
line1_p2_x – X position of second point on line 1.
line1_p2_y – Y position of second point on line 1.
line2_p1_x – X position of first point on line 2.
line2_p1_y – Y position of first point on line 2.
line2_p2_x – X position of second point on line 2.
line2_p2_y – Y position of second point on line 2.

Returns:

New (x, y) position, or (NaN, NaN) if unbuildable (parallel lines).

pylinkage.solver.joints.solve_linear(current_x: float, current_y: float, circle_x: float, circle_y: float, radius: float, line_p1_x: float, line_p1_y: float, line_p2_x: float, line_p2_y: float) → tuple[float, float]

Solve linear joint using circle-line intersection.

The linear joint slides along a line while maintaining a fixed distance from a revolute anchor: - Circle: centered at revolute anchor with given radius - Line: passing through line_p1 and line_p2

When two solutions exist, picks the nearest to current position (hysteresis).

Parameters:

current_x – Current X position (for disambiguation).
current_y – Current Y position (for disambiguation).
circle_x – X position of revolute anchor (circle center).
circle_y – Y position of revolute anchor (circle center).
radius – Distance to revolute anchor (circle radius).
line_p1_x – X position of first line-defining point.
line_p1_y – Y position of first line-defining point.
line_p2_x – X position of second line-defining point.
line_p2_y – Y position of second line-defining point.

Returns:

New (x, y) position, or (NaN, NaN) if unbuildable.

pylinkage.solver.joints.solve_linear_actuator(current_extension: float, direction: float, anchor_x: float, anchor_y: float, angle: float, stroke: float, velocity: float, dt: float) → tuple[float, float, float, float]

Solve linear actuator position using oscillating linear motion.

The linear actuator moves along a line at constant velocity, reversing direction when it reaches stroke limits.

Parameters:

current_extension – Current extension from anchor (0 to stroke).
direction – Current movement direction (+1.0 or -1.0).
anchor_x – X position of the anchor (fixed end).
anchor_y – Y position of the anchor (fixed end).
angle – Direction angle in radians (from +x axis).
stroke – Maximum extension distance.
velocity – Linear velocity magnitude (units per time step).
dt – Time step.

Returns:

Tuple of (new_x, new_y, new_extension, new_direction).

pylinkage.solver.joints.solve_oscillating_cam_follower(pivot_x: float, pivot_y: float, arm_length: float, arm_angle: float) → tuple[float, float]

Solve oscillating cam follower position.

The follower is at the end of an arm rotating about a pivot point.

Parameters:

pivot_x – X position of pivot point.
pivot_y – Y position of pivot point.
arm_length – Length of the follower arm.
arm_angle – Current angle of the arm (radians from +x).

Returns:

New (x, y) position of the follower output.

pylinkage.solver.joints.solve_revolute(current_x: float, current_y: float, p0_x: float, p0_y: float, r0: float, p1_x: float, p1_y: float, r1: float) → tuple[float, float]

Solve revolute joint using circle-circle intersection.

The revolute joint is positioned at the intersection of two circles: - Circle 1: centered at parent 0 with radius r0 - Circle 2: centered at parent 1 with radius r1

When two solutions exist, picks the nearest to current position (hysteresis).

Parameters:

current_x – Current X position (for disambiguation).
current_y – Current Y position (for disambiguation).
p0_x – X position of first parent (circle 1 center).
p0_y – Y position of first parent (circle 1 center).
r0 – Distance to first parent (circle 1 radius).
p1_x – X position of second parent (circle 2 center).
p1_y – Y position of second parent (circle 2 center).
r1 – Distance to second parent (circle 2 radius).

Returns:

New (x, y) position, or (NaN, NaN) if unbuildable.

pylinkage.solver.joints.solve_translating_cam_follower(guide_x: float, guide_y: float, guide_angle: float, displacement: float) → tuple[float, float]

Solve translating cam follower position.

The follower moves along a line from the guide point at the specified angle, with displacement determining position along the line.

Parameters:

guide_x – X position of guide reference point.
guide_y – Y position of guide reference point.
guide_angle – Direction of follower motion (radians from +x).
displacement – Distance from guide along the guide axis.

Returns:

New (x, y) position of the follower.

pylinkage.solver.simulation module

Pure-numba simulation loop for linkage mechanisms.

This module provides the core simulation functions that operate on numeric arrays, avoiding Python object overhead for maximum performance.

pylinkage.solver.simulation.first_nan_step(trajectory: ndarray) → int

Find the first step with NaN positions.

Parameters:: trajectory – Array of shape (iterations, n_joints, 2).
Returns:: Index of first step with NaN, or -1 if no NaN found.

pylinkage.solver.simulation.has_nan_positions(trajectory: ndarray) → bool

Check if trajectory contains any NaN positions.

Parameters:: trajectory – Array of shape (iterations, n_joints, 2).
Returns:: True if any position is NaN (indicating unbuildable configuration).

pylinkage.solver.simulation.simulate(positions: ndarray, constraints: ndarray, joint_types: ndarray, parent_indices: ndarray, constraint_offsets: ndarray, solve_order: ndarray, iterations: int, dt: float) → ndarray

Run full simulation and return trajectory.

This is the main entry point for numba-optimized simulation. All arrays are processed in a tight loop with no Python overhead.

Parameters:

positions – Initial joint positions, shape (n_joints, 2). Will be modified to contain final positions.
constraints – Flat array of constraint values.
joint_types – Joint type code for each joint, shape (n_joints,).
parent_indices – Parent joint indices, shape (n_joints, max_parents).
constraint_offsets – Start index in constraints for each joint.
solve_order – Indices of joints to solve, in order.
iterations – Number of simulation steps to run.
dt – Time step for crank rotation.

Returns:

Trajectory array of shape (iterations, n_joints, 2) containing all joint positions at each step. If any position becomes NaN (unbuildable configuration), subsequent positions may also be NaN.

pylinkage.solver.simulation.simulate_with_kinematics(positions: ndarray, velocities: ndarray, accelerations: ndarray, constraints: ndarray, joint_types: ndarray, parent_indices: ndarray, constraint_offsets: ndarray, solve_order: ndarray, omega_values: ndarray, alpha_values: ndarray, crank_indices: ndarray, iterations: int, dt: float) → tuple[ndarray, ndarray, ndarray]

Run full simulation with velocity and acceleration computation.

This extends the basic simulate() function to also compute velocities and accelerations at each step.

Parameters:

positions – Initial joint positions, shape (n_joints, 2). Will be modified to contain final positions.
velocities – Initial velocities, shape (n_joints, 2). Will be modified to contain final velocities.
accelerations – Initial accelerations, shape (n_joints, 2). Will be modified to contain final accelerations.
constraints – Flat array of constraint values.
joint_types – Joint type code for each joint, shape (n_joints,).
parent_indices – Parent joint indices, shape (n_joints, max_parents).
constraint_offsets – Start index in constraints for each joint.
solve_order – Indices of joints to solve, in order.
omega_values – Angular velocities for crank joints (rad/s).
alpha_values – Angular accelerations for crank joints (rad/s²).
crank_indices – Indices of crank joints in the joints array.
iterations – Number of simulation steps to run.
dt – Time step for crank rotation.

Returns:

Tuple of (positions_trajectory, velocities_trajectory, accelerations_trajectory), each with shape (iterations, n_joints, 2).

pylinkage.solver.simulation.step_single(positions: ndarray, constraints: ndarray, joint_types: ndarray, parent_indices: ndarray, constraint_offsets: ndarray, solve_order: ndarray, dt: float) → None

Perform one simulation step, updating positions in-place.

Parameters:

positions – Joint positions, shape (n_joints, 2). Modified in-place.
constraints – Flat array of constraint values.
joint_types – Joint type code for each joint, shape (n_joints,).
parent_indices – Parent joint indices, shape (n_joints, max_parents).
constraint_offsets – Start index in constraints for each joint.
solve_order – Indices of joints to solve, in order.
dt – Time step for crank rotation.

pylinkage.solver.simulation.step_single_acceleration(positions: ndarray, velocities: ndarray, accelerations: ndarray, constraints: ndarray, joint_types: ndarray, parent_indices: ndarray, constraint_offsets: ndarray, solve_order: ndarray, omega_values: ndarray, alpha_values: ndarray, crank_indices: ndarray) → None

Compute accelerations for all joints after velocities are solved.

Must be called after step_single_velocity() has updated velocities.

Parameters:

positions – Joint positions, shape (n_joints, 2).
velocities – Joint velocities, shape (n_joints, 2).
accelerations – Output accelerations, shape (n_joints, 2). Modified in-place.
constraints – Flat array of constraint values.
joint_types – Joint type code for each joint, shape (n_joints,).
parent_indices – Parent joint indices, shape (n_joints, max_parents).
constraint_offsets – Start index in constraints for each joint.
solve_order – Indices of joints to solve, in order.
omega_values – Angular velocities for crank joints (rad/s).
alpha_values – Angular accelerations for crank joints (rad/s²).
crank_indices – Indices of crank joints in joints array.

pylinkage.solver.simulation.step_single_velocity(positions: ndarray, velocities: ndarray, constraints: ndarray, joint_types: ndarray, parent_indices: ndarray, constraint_offsets: ndarray, solve_order: ndarray, omega_values: ndarray, crank_indices: ndarray) → None

Compute velocities for all joints after positions are solved.

Must be called after step_single() has updated positions.

Parameters:

positions – Joint positions, shape (n_joints, 2).
velocities – Output velocities, shape (n_joints, 2). Modified in-place.
constraints – Flat array of constraint values.
joint_types – Joint type code for each joint, shape (n_joints,).
parent_indices – Parent joint indices, shape (n_joints, max_parents).
constraint_offsets – Start index in constraints for each joint.
solve_order – Indices of joints to solve, in order.
omega_values – Angular velocities for crank joints (rad/s).
crank_indices – Indices of crank joints in joints array.

pylinkage.solver.solve module

High-level solving functions for linkage mechanisms.

This module provides the primary API for solving linkage kinematics using Assur group decomposition. It bridges between the structural analysis (assur module) and the low-level numba solvers (solver.joints).

The key functions are: - solve_group(): Solve a single Assur group given anchor positions and dimensions - solve_decomposition(): Solve all groups in decomposition order

These functions are the canonical location for solving behavior. The Assur group classes themselves are pure topology (no solve methods).

pylinkage.solver.solve.solve_decomposition(result: DecompositionResult, dimensions: Dimensions, initial_positions: dict[NodeId, Coord] | None = None) → dict[NodeId, Coord]

Solve the kinematics using a decomposition result and dimensions.

This function computes the positions of all nodes by solving the Assur groups in order. It is the canonical location for decomposition-based solving.

Parameters:

result – The decomposition result with groups in solving order. Must have a valid graph reference.
dimensions – Dimensions object containing node positions and edge distances.
initial_positions – Optional override positions for nodes. If not provided, positions from dimensions are used.

Returns:

Dictionary mapping all node IDs to their computed (x, y) positions.

Raises:

UnbuildableError – If any group cannot be solved.
ValueError – If required positions are missing or graph is None.

Example

>>> from pylinkage.assur import decompose_assur_groups
>>> from pylinkage.dimensions import Dimensions
>>> result = decompose_assur_groups(graph)
>>> dims = Dimensions(node_positions={...}, edge_distances={...})
>>> positions = solve_decomposition(result, dims)
>>> print(positions)

pylinkage.solver.solve.solve_group(group: AssurGroup, positions: dict[NodeId, Coord], dimensions: Dimensions, hint_positions: dict[NodeId, Coord] | None = None) → dict[NodeId, Coord]

Solve positions for an Assur group.

This is the main dispatch function that routes to the appropriate solver based on group type (RRR, RRP, etc.).

Parameters:

group – The Assur group to solve (topology only).
positions – Known positions of anchor nodes. Must contain entries for all nodes in group.anchor_nodes.
dimensions – Dimensions object containing edge distances.
hint_positions – Optional positions to use as hints for disambiguation when there are multiple solutions.

Returns:

Dict mapping internal node IDs to their computed (x, y) positions.

Raises:

UnbuildableError – If the group cannot be solved geometrically.
NotImplementedError – If the group type is not yet supported.
ValueError – If required constraints or positions are missing.

Example

>>> from pylinkage.assur import DyadRRR, decompose_assur_groups
>>> from pylinkage.dimensions import Dimensions
>>> result = decompose_assur_groups(graph)
>>> dims = Dimensions(node_positions={...}, edge_distances={...})
>>> positions = {n: dims.get_node_position(n) for n in result.ground}
>>> for group in result.groups:
...     new_pos = solve_group(group, positions, dims)
...     positions.update(new_pos)

pylinkage.solver.types module

Numba-compatible data structures for the solver.

This module defines the numeric representation of linkages for use in the pure-numba simulation loop.

class pylinkage.solver.types.SolverData(positions: ndarray[tuple[Any, ...], dtype[float64]], constraints: ndarray[tuple[Any, ...], dtype[float64]], joint_types: ndarray[tuple[Any, ...], dtype[int32]], parent_indices: ndarray[tuple[Any, ...], dtype[int32]], constraint_offsets: ndarray[tuple[Any, ...], dtype[int32]], constraint_counts: ndarray[tuple[Any, ...], dtype[int32]], solve_order: ndarray[tuple[Any, ...], dtype[int32]], velocities: ndarray[tuple[Any, ...], dtype[float64]] | None = None, accelerations: ndarray[tuple[Any, ...], dtype[float64]] | None = None, omega_values: ndarray[tuple[Any, ...], dtype[float64]] | None = None, alpha_values: ndarray[tuple[Any, ...], dtype[float64]] | None = None, crank_indices: ndarray[tuple[Any, ...], dtype[int32]] | None = None)

Bases: object

Numba-compatible representation of a linkage.

All data is stored in contiguous numpy arrays for efficient access from numba-compiled code.

Variables:

positions (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Joint positions, shape (n_joints, 2).
constraints (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Flat array of all constraint values.
joint_types (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Joint type code for each joint, shape (n_joints,).
parent_indices (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Parent joint indices, shape (n_joints, MAX_PARENTS). Unused slots contain -1.
constraint_offsets (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Start index in constraints array for each joint.
constraint_counts (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Number of constraints for each joint.
solve_order (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Indices of joints in solving order.
velocities (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]] | None) – Joint velocities, shape (n_joints, 2). Optional.
accelerations (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]] | None) – Joint accelerations, shape (n_joints, 2). Optional.
omega_values (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]] | None) – Angular velocities for crank joints (rad/s), shape (n_cranks,).
alpha_values (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]] | None) – Angular accelerations for crank joints (rad/s²), shape (n_cranks,).
crank_indices (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]] | None) – Indices of crank joints in the joints array, shape (n_cranks,).

accelerations: ndarray[tuple[Any, ...], dtype[float64]] | None = None

alpha_values: ndarray[tuple[Any, ...], dtype[float64]] | None = None

constraint_counts: ndarray[tuple[Any, ...], dtype[int32]]

constraint_offsets: ndarray[tuple[Any, ...], dtype[int32]]

constraints: ndarray[tuple[Any, ...], dtype[float64]]

copy() → SolverData: Create a deep copy of the solver data.

crank_indices: ndarray[tuple[Any, ...], dtype[int32]] | None = None

joint_types: ndarray[tuple[Any, ...], dtype[int32]]

property n_constraints: int: Total number of constraints.

property n_joints: int: Number of joints in the linkage.

omega_values: ndarray[tuple[Any, ...], dtype[float64]] | None = None

parent_indices: ndarray[tuple[Any, ...], dtype[int32]]

positions: ndarray[tuple[Any, ...], dtype[float64]]

solve_order: ndarray[tuple[Any, ...], dtype[int32]]

velocities: ndarray[tuple[Any, ...], dtype[float64]] | None = None

pylinkage.solver.velocity module

Per-joint numba velocity solvers for kinematic analysis.

Each function computes the linear velocity of a specific joint type given the parent velocities, positions, and constraints. All functions are numba-compiled for maximum performance.

Velocity equations are derived by analytic differentiation of position equations, using implicit differentiation for constraint-based joints.

pylinkage.solver.velocity.solve_crank_velocity(x: float, y: float, anchor_x: float, anchor_y: float, anchor_vx: float, anchor_vy: float, radius: float, omega: float) → tuple[float, float]

Compute crank velocity using direct differentiation.

The crank position is: (anchor + r*cos(θ), anchor + r*sin(θ)) Differentiating: v = anchor_v + r*ω*(-sin(θ), cos(θ))

Parameters:

x – Current X position of the crank.
y – Current Y position of the crank.
anchor_x – X position of the anchor.
anchor_y – Y position of the anchor.
anchor_vx – X velocity of the anchor (usually 0).
anchor_vy – Y velocity of the anchor (usually 0).
radius – Distance from anchor to crank.
omega – Angular velocity in rad/s.

Returns:

Velocity (vx, vy) of the crank.

pylinkage.solver.velocity.solve_fixed_velocity(x: float, y: float, p0_x: float, p0_y: float, p0_vx: float, p0_vy: float, p1_x: float, p1_y: float, p1_vx: float, p1_vy: float, radius: float, angle: float) → tuple[float, float]

Compute fixed joint velocity using direct differentiation.

The fixed joint position is:: x = p0_x + r * cos(α + γ) y = p0_y + r * sin(α + γ)

Where α = atan2(p1_y - p0_y, p1_x - p0_x) is the base angle.

Differentiating:: dα/dt = [(p1_x - p0_x)(p1_vy - p0_vy) - (p1_y - p0_y)(p1_vx - p0_vx)] / d² vx = p0_vx - r * (dα/dt) * sin(α + γ) vy = p0_vy + r * (dα/dt) * cos(α + γ)

Parameters:

x – Current X position of the joint (unused, for consistency).
y – Current Y position of the joint (unused, for consistency).
p0_x – X position of first parent (origin).
p0_y – Y position of first parent (origin).
p0_vx – X velocity of first parent.
p0_vy – Y velocity of first parent.
p1_x – X position of second parent (angle reference).
p1_y – Y position of second parent (angle reference).
p1_vx – X velocity of second parent.
p1_vy – Y velocity of second parent.
radius – Distance from first parent.
angle – Angle offset from the parent-to-parent direction.

Returns:

Velocity (vx, vy) of the joint, or (NaN, NaN) if singular.

pylinkage.solver.velocity.solve_prismatic_velocity(x: float, y: float, circle_x: float, circle_y: float, circle_vx: float, circle_vy: float, radius: float, line_p1_x: float, line_p1_y: float, line_p1_vx: float, line_p1_vy: float, line_p2_x: float, line_p2_y: float, line_p2_vx: float, line_p2_vy: float) → tuple[float, float]

Compute prismatic joint velocity using implicit differentiation.

The prismatic joint satisfies: 1. Circle constraint: (x - cx)² + (y - cy)² = r² 2. Line constraint: point lies on line through p1 and p2

The line constraint in 2D:: (x - p1_x)(p2_y - p1_y) - (y - p1_y)(p2_x - p1_x) = 0

Differentiating both constraints yields a 2x2 linear system.

Parameters:

x – Current X position of the joint.
y – Current Y position of the joint.
circle_x – X position of circle center (revolute anchor).
circle_y – Y position of circle center.
circle_vx – X velocity of circle center.
circle_vy – Y velocity of circle center.
radius – Circle radius.
line_p1_x – X position of first line point.
line_p1_y – Y position of first line point.
line_p1_vx – X velocity of first line point.
line_p1_vy – Y velocity of first line point.
line_p2_x – X position of second line point.
line_p2_y – Y position of second line point.
line_p2_vx – X velocity of second line point.
line_p2_vy – Y velocity of second line point.

Returns:

Velocity (vx, vy) of the joint, or (NaN, NaN) if singular.

pylinkage.solver.velocity.solve_revolute_velocity(x: float, y: float, p0_x: float, p0_y: float, p0_vx: float, p0_vy: float, p1_x: float, p1_y: float, p1_vx: float, p1_vy: float) → tuple[float, float]

Compute revolute joint velocity using implicit differentiation.

The revolute joint satisfies two distance constraints:

(x - p0_x)² + (y - p0_y)² = r0² (x - p1_x)² + (y - p1_y)² = r1²

Differentiating with respect to time yields a 2x2 linear system:

A * [vx, vy]ᵀ = b

Where:

A = [[x - p0_x, y - p0_y],: [x - p1_x, y - p1_y]]
b = [(x - p0_x)*p0_vx + (y - p0_y)*p0_vy,: (x - p1_x)*p1_vx + (y - p1_y)*p1_vy]

Parameters:

x – Current X position of the joint.
y – Current Y position of the joint.
p0_x – X position of first parent.
p0_y – Y position of first parent.
p0_vx – X velocity of first parent.
p0_vy – Y velocity of first parent.
p1_x – X position of second parent.
p1_y – Y position of second parent.
p1_vx – X velocity of second parent.
p1_vy – Y velocity of second parent.

Returns:

Velocity (vx, vy) of the joint, or (NaN, NaN) if singular.

pylinkage.solver.velocity.solve_rigid_body_velocity(x: float, y: float, p0_x: float, p0_y: float, p0_vx: float, p0_vy: float, p1_x: float, p1_y: float, p1_vx: float, p1_vy: float) → tuple[float, float]

Compute the velocity of a third point on a rigid body.

Where solve_revolute_velocity() intersects two distance constraints, this propagates the motion of the body the point sits on. Two points of a planar rigid body fix its angular velocity:

v1 = v0 + ω ẑ × (p1 - p0) ⟹ ω = [dx * (p1_vy - p0_vy) - dy * (p1_vx - p0_vx)] / d²

with d = p1 - p0. The queried point then follows directly:

v = v0 + ω ẑ × (p - p0)

Both formulations agree wherever both apply. This one stays well-conditioned when the three points are collinear — the case that makes the constraint-intersection Jacobian singular — because it needs two distinct points on the body rather than two independent constraint directions. A ternary link carrying a coupler point on the line through its other two ports is the common example, and is a normal design rather than a degenerate one.

Parameters:

x – Current X position of the queried point.
y – Current Y position of the queried point.
p0_x – X position of the first body point.
p0_y – Y position of the first body point.
p0_vx – X velocity of the first body point.
p0_vy – Y velocity of the first body point.
p1_x – X position of the second body point.
p1_y – Y position of the second body point.
p1_vx – X velocity of the second body point.
p1_vy – Y velocity of the second body point.

Returns:

Velocity (vx, vy) of the queried point, or (NaN, NaN) if the two body points coincide.

Module contents

Pure-numba simulation solver for linkage mechanisms.

This module provides a high-performance simulation backend that uses numba JIT compilation to eliminate Python overhead in the hot loop.

The solver operates on numeric arrays rather than Python objects, achieving significant speedups for repeated simulations (e.g., during optimization).

Basic usage:

>>> from pylinkage import Linkage
>>> linkage = create_my_linkage()
>>> trajectory = linkage.step_fast(iterations=1000)
>>> # trajectory.shape == (1000, n_joints, 2)

For direct access to the solver:

>>> from pylinkage.solver import linkage_to_solver_data, simulate
>>> data = linkage_to_solver_data(linkage)
>>> trajectory = simulate(
...     data.positions, data.constraints, data.joint_types,
...     data.parent_indices, data.constraint_offsets,
...     data.solve_order, iterations=1000, dt=1.0
... )

class pylinkage.solver.SolverData(positions: ndarray[tuple[Any, ...], dtype[float64]], constraints: ndarray[tuple[Any, ...], dtype[float64]], joint_types: ndarray[tuple[Any, ...], dtype[int32]], parent_indices: ndarray[tuple[Any, ...], dtype[int32]], constraint_offsets: ndarray[tuple[Any, ...], dtype[int32]], constraint_counts: ndarray[tuple[Any, ...], dtype[int32]], solve_order: ndarray[tuple[Any, ...], dtype[int32]], velocities: ndarray[tuple[Any, ...], dtype[float64]] | None = None, accelerations: ndarray[tuple[Any, ...], dtype[float64]] | None = None, omega_values: ndarray[tuple[Any, ...], dtype[float64]] | None = None, alpha_values: ndarray[tuple[Any, ...], dtype[float64]] | None = None, crank_indices: ndarray[tuple[Any, ...], dtype[int32]] | None = None)

Bases: object

Numba-compatible representation of a linkage.

All data is stored in contiguous numpy arrays for efficient access from numba-compiled code.

Variables:

positions (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Joint positions, shape (n_joints, 2).
constraints (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Flat array of all constraint values.
joint_types (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Joint type code for each joint, shape (n_joints,).
parent_indices (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Parent joint indices, shape (n_joints, MAX_PARENTS). Unused slots contain -1.
constraint_offsets (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Start index in constraints array for each joint.
constraint_counts (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Number of constraints for each joint.
solve_order (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]]) – Indices of joints in solving order.
velocities (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]] | None) – Joint velocities, shape (n_joints, 2). Optional.
accelerations (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]] | None) – Joint accelerations, shape (n_joints, 2). Optional.
omega_values (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]] | None) – Angular velocities for crank joints (rad/s), shape (n_cranks,).
alpha_values (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]] | None) – Angular accelerations for crank joints (rad/s²), shape (n_cranks,).
crank_indices (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.int32]] | None) – Indices of crank joints in the joints array, shape (n_cranks,).

accelerations: ndarray[tuple[Any, ...], dtype[float64]] | None = None

alpha_values: ndarray[tuple[Any, ...], dtype[float64]] | None = None

constraint_counts: ndarray[tuple[Any, ...], dtype[int32]]

constraint_offsets: ndarray[tuple[Any, ...], dtype[int32]]

constraints: ndarray[tuple[Any, ...], dtype[float64]]

copy() → SolverData: Create a deep copy of the solver data.

crank_indices: ndarray[tuple[Any, ...], dtype[int32]] | None = None

joint_types: ndarray[tuple[Any, ...], dtype[int32]]

property n_constraints: int: Total number of constraints.

property n_joints: int: Number of joints in the linkage.

omega_values: ndarray[tuple[Any, ...], dtype[float64]] | None = None

parent_indices: ndarray[tuple[Any, ...], dtype[int32]]

positions: ndarray[tuple[Any, ...], dtype[float64]]

solve_order: ndarray[tuple[Any, ...], dtype[int32]]

velocities: ndarray[tuple[Any, ...], dtype[float64]] | None = None

pylinkage.solver.simulate(positions: ndarray, constraints: ndarray, joint_types: ndarray, parent_indices: ndarray, constraint_offsets: ndarray, solve_order: ndarray, iterations: int, dt: float) → ndarray

Run full simulation and return trajectory.

This is the main entry point for numba-optimized simulation. All arrays are processed in a tight loop with no Python overhead.

Parameters:

positions – Initial joint positions, shape (n_joints, 2). Will be modified to contain final positions.
constraints – Flat array of constraint values.
joint_types – Joint type code for each joint, shape (n_joints,).
parent_indices – Parent joint indices, shape (n_joints, max_parents).
constraint_offsets – Start index in constraints for each joint.
solve_order – Indices of joints to solve, in order.
iterations – Number of simulation steps to run.
dt – Time step for crank rotation.

Returns:

Trajectory array of shape (iterations, n_joints, 2) containing all joint positions at each step. If any position becomes NaN (unbuildable configuration), subsequent positions may also be NaN.

pylinkage.solver.simulate_with_kinematics(positions: ndarray, velocities: ndarray, accelerations: ndarray, constraints: ndarray, joint_types: ndarray, parent_indices: ndarray, constraint_offsets: ndarray, solve_order: ndarray, omega_values: ndarray, alpha_values: ndarray, crank_indices: ndarray, iterations: int, dt: float) → tuple[ndarray, ndarray, ndarray]

Run full simulation with velocity and acceleration computation.

This extends the basic simulate() function to also compute velocities and accelerations at each step.

Parameters:

positions – Initial joint positions, shape (n_joints, 2). Will be modified to contain final positions.
velocities – Initial velocities, shape (n_joints, 2). Will be modified to contain final velocities.
accelerations – Initial accelerations, shape (n_joints, 2). Will be modified to contain final accelerations.
constraints – Flat array of constraint values.
joint_types – Joint type code for each joint, shape (n_joints,).
parent_indices – Parent joint indices, shape (n_joints, max_parents).
constraint_offsets – Start index in constraints for each joint.
solve_order – Indices of joints to solve, in order.
omega_values – Angular velocities for crank joints (rad/s).
alpha_values – Angular accelerations for crank joints (rad/s²).
crank_indices – Indices of crank joints in the joints array.
iterations – Number of simulation steps to run.
dt – Time step for crank rotation.

Returns:

Tuple of (positions_trajectory, velocities_trajectory, accelerations_trajectory), each with shape (iterations, n_joints, 2).