Mug and Spoon PCBO

This note describes the PCBO that runs in code/mug_spoon/pcbo3d.cpp (binary pcbo3d_v37 and later), in the configuration that performed best on 22 Oct 2026. Every definition below is the one the code computes. Parameter names in typewriter are the environment variables of pcbo3d; their values in the best configuration are collected in Section 10.

The problem is stated generically: two rigid bodies, finitely many contact constraints, a static-equilibrium residual. The mug and spoon is the instance.

1. Spaces

1.1 Configuration space

A configuration is a pair of rigid poses, the mug $M$ and the spoon $S$:

with $p$ positions and $q$ unit quaternions ($\pm q$ identified). Gravity and the table are invariant under the planar Euclidean group

acting on both bodies at once. The physically distinct configurations form the quotient $X/G$, of dimension $12 - 3 = 9$.

Gauge-free distance. For $a, b \in X$ let $g_{ab} \in G$ align $b$ to $a$: it rotates $b$ about $e_z$ by the yaw of $R_{q_M^b} R_{q_M^a}^{\top}$ and translates it so that the mugs' horizontal positions coincide (gauge_align). With $L = 0.1\,$m per radian (L_TWIST),

This is gdist. All distances between landed configurations ("images") below are $d$.

1.2 Constraints, strata and the feasible set

The model has finitely many geometric features $j = 1,\dots,F$ ($F = 413$): gaps $g_j(x)$ between pairs of geoms (mug parts, spoon parts, the table). They are signed distances, exact for the curved pairs. A feature is active when $g_j(x) \le b$ ($b = 2\cdot10^{-5}$, band). A configuration is valid when no gap is below $-\epsilon_{\text{pen}}$ ($\epsilon_{\text{pen}} = 2\cdot 10^{-7}$) and at least one active contact joins mug and spoon.

For a set $H$ of held features, the contact stratum is

Equilibrium residual. At a valid $x$ with active contacts $c = 1,\dots,n_c$, each with normal $n_c$ and friction coefficient $\mu = 10$, the circular friction cone is replaced by the inscribed polygonal cone with $K = 16$ generators $\{f_{c,k}\}$. The residual is the non-negative least-squares value

where $W_{c,k} \in \mathbb{R}^{12}$ is the wrench of generator $f_{c,k}$ on both bodies (force and torque about each centre of mass) and $b$ is the wrench that balances gravity, in units $F^\ast = 0.3 \cdot 9.81\,$N and $L^\ast = 0.1\,$m. It is computed by Lawson–Hanson NNLS with a KKT check; $V_{\text{in}}$ is the residual norm. Because the cone is inscribed, $V_{\text{in}} = 0$ certifies an equilibrium. A circumscribed cone gives $V_{\text{out}} \le V$, which certifies infeasibility when positive.

• Exact: $V_{\text{in}}(x) \le 10^{-10}$.
• Feasible set: $\mathcal{F} = \{x \text{ valid} : V_{\text{in}}(x) \le 10^{-10}\}$.
• Flat (robust) set: $\mathcal{F}_{\text{flat}} = \{x \in \mathcal{F} : \Phi(x) \ge \tfrac12\}$. Here $\Phi(x)$ is the share of $K_f = 8$ random unit tangent directions $u$ of the active stratum (gauge excluded) along which the step $x \mapsto \text{proj}_{\Sigma_H}(x \oplus \rho u)$, with $\rho = 10^{-5}$, stays exact (flat_share).

The objects of interest are the connected components of $\mathcal{F}_{\text{flat}}/G$ ("scenes"). The landed flat states have intrinsic dimension about 3-4: the Levina–Bickel MLE gives 3.6-4.1 on all flat landings and 2.8-3.1 within the largest component. This matches strata of dimension $9 - r$ for held ranks $r = 5,\dots,8$.

1.3 The search space and the landing map

PCBO does not move in $X$. It moves in the input space of the landing law:

a 10-dimensional manifold with boundary. A point $z = (q_M, h, q_S, p_S)$ fixes the mug's orientation, its lowest point at height $h$ above the table, the spoon's orientation and origin. The mug's horizontal position is the origin, which fixes the translational gauge. The reference law $\mu_0$ on $Z$ is uniform: Haar on both $SO(3)$ factors, uniform on $[0, \bar h]$ and on $B$.

Chart. At $z$ the normalised coordinates of a nearby $z'$ are

(in_diff; $\log$ is the rotation vector). A step $s \in \mathbb{R}^{10}$ is applied by $z \oplus s$: left multiplication by $\exp(\pi s_{\text{rot}})$ on each orientation, and reflection at the boundary for $h$ and $p_S$ (in_step). Norms $|\cdot|$ of steps are taken in these coordinates.

Landing map $\Lambda_r : Z \to X \cup \{\varnothing\}$, for a target rank $r \in \{5,6,7,8\}$:

1. Form the initial pose $x_0(z)$.
2. Order the features by increasing gap at $x_0$, and hold, in that order, those whose gradient raises the rank of the held set together with the gauge directions. A row raises it when its component orthogonal to the span exceeds $0.05$. Continue until rank $r$; every penetrating feature is also held.
3. Project $x_0$ onto $\Sigma_H$ by damped Gauss–Newton (minimum-norm steps).
4. If new penetrations appear, add them to $H$ and project again, up to 4 rounds.
5. If no round succeeds, the landing fails: $\Lambda_r(z) = \varnothing$, at no oracle cost.

$\Lambda_r$ is piecewise smooth and many-to-one: whole basins of $Z$ collapse onto one image. This fact governs most of the design below. Each particle $i$ keeps a fixed rank $r_i = 5 + (i \bmod 4)$; write $\Lambda_i = \Lambda_{r_i}$ and $x_i = \Lambda_i(z_i)$ (the input pose if the landing failed).

2. Objective and oracle cost

Every flat state, found long ago or just now, is a global minimiser ($E = 0$). That is why plain PCBO concentrates on what it has found (Section 9).

Cost. One oracle call per evaluated landing. For an exact landing, add the flatness test, run sequentially (FLAT_EARLY=1): the same directions are drawn in the same order, and the test stops at 4 exact steps (flat) or 5 failed ones (not flat). Steps whose projection fails are not evaluated. The verdict is identical to the full 8-step test (checked: 1,329 of 1,329 exact landings agree). An exact landing costs about $4.8$ calls instead of $9$. Failed landings and all geometry cost nothing. Budgets are counted in oracle calls.

3. The metric

The metric needs only a metric space $(X, d_X)$, a group $G$ of isometries that leaves the constraints and the objective unchanged, and the objective $E$. It satisfies the metric axioms, and it keeps three jobs apart: geometry (the metric), which points are solutions (the objective, used only in the consensus weights, Section 4), and where the next evaluation should go (diversity, Section 6).

3.1 Quotient placement metric

the root-mean-square displacement of each body's matter ($\rho_B$ its mass distribution, $c_B$ its centre of mass, $\Sigma_B$ its second moment; every body weighted equally). The table's symmetry group $G = SE(2)$ (rotation about the vertical, horizontal translation) is divided out exactly:

a planar Procrustes problem with a closed-form minimiser. $d_G$ is a metric on scenes: symmetric, zero only for symmetric copies of the same scene, and the triangle inequality holds exactly because the minimum over $G$ is global.

3.2 Contact term

with $H$ the set of active contact constraints (gap within the contact band) and $\lambda$ the bodies' mean RMS radius (4.66 cm): gaining or losing a contact counts as much as moving a body by its own size. As the $\ell^2$ sum of $d_G$ and the square root of a Hamming distance, $d$ is a metric, unchanged by the table's symmetry. On the solution set, the corresponding path metric is an extended metric whose infinite distances are exactly the disconnected components.

3.3 Diversity and the archive

Let $A_t$ be every configuration evaluated so far. The novelty of a proposal $x$ is $\nu_t(x) = \min_{y \in A_t} d(x, y)$, and the archive's spacing $\varepsilon_t$ is the median distance of its points to their nearest neighbours. A proposal is evaluated with probability $\min\{1, (\nu_t(x)/\varepsilon_t)^m\}$, $m = 2$ (Section 6): an $\varepsilon$-cover of the solution set built online. The crowding of an evaluated state $j$ is $G_j = -\log(nn_j/\varepsilon_t)$, $nn_j$ its distance to its nearest other archive point (used by the consensus, Section 4). In the dynamics of Section 5 the conformal factor is 1, and the noise floor is $\varepsilon_t$ (converted to inputs by the length scale $s$), so it shrinks as the cover gets finer.

4. Consensus points

Feasibility first, crowding second (CB2O, Fornasier et al.). Each particle's consensus averages only feasible particles, those at exact equilibria ($E \le 10^{-10}$), with weights $\exp(-\alpha\,[D_{ij}^2 + G_j])$, $G_j$ the crowding of Section 3.3, so the consensus leans towards sparsely explored solutions; when no feasible particle is within reach, the weights of Section 4.1 apply. $\alpha$ is chosen per particle so that each consensus averages $8N/128$ particles (62 for $N = 1000$), which keeps the kernel fixed as $N$ grows.

4.1 Polarised weights (Fornasier–Sun, arXiv 2403.06662, eq. (16))

Build the $k$-nearest-neighbour graph ($k = 8$) on the images $\{x_j\}$. Edge lengths are conformal,

and $D_{ij}$ is the shortest-path (graph-geodesic) distance, following Bernstein–de Silva–Langford–Tenenbaum. Particle $i$ weights particle $j$ by

with $\alpha = 10$, $\vartheta = 0.05$ (VTH), $\varkappa = 0.01$ (VKAP), normalised over $j$. As in PCBO, $\arg\min_j A(j,i)$ is close to the minimiser nearest $i$.

4.2 Continuous consensus point (CONS_CONT=1)

The weighted mean is taken in the chart at $z_i$ (one step of a weighted Karcher mean):

It is a point of $Z$, not one of the particles. A medoid consensus can only move to existing particles, so repulsion between medoids cannot reach empty space. The continuous point can.

4.3 Balanced repulsion of the consensus points (Massimo's (b), in the form of Fornasier–Hütter)

The consensus points $C = (c_1,\dots,c_N)$ minimise

Up to the normalisation $N$ against $N-1$ and an additive constant, the bracket is $\tfrac{N}{2}$ times the energy distance between the empirical measure of the consensus points and a target $\omega$:

the MMD with the conditionally negative definite kernel $-|\cdot|^p$, $1 \le p < 2$. By Fornasier–Hütter (arXiv 1310.1120), the balanced form, attraction to $\omega$ and repulsion among the points with the same kernel, has $\omega$ as its unique minimiser in the limit. Repulsion alone has no target. The target is Massimo's suggestion (d), "uniform on the hull of the found solutions":

where $\mathcal{Z}^{\text{flat}}_t$ are the inputs that produced flat landings so far (a pool of 256 points per iteration). $\mathcal{B}$ is decreased by 20 Jacobi gradient steps of size $\eta = 0.05$, each coordinate clipped:

Write $\ell_i = |c_i \ominus z_i|$ and $u_i = (c_i \ominus z_i)/\ell_i$.

5. The dynamics

5.1 The SDE (chart coordinates, frozen metric)

with $\lambda = 1$, $\sigma = 0.2$ and the noise floor $\tilde\kappa_t = \varepsilon_t / 0.2\,\text{m}$, the archive spacing of Section 3.3 in input units: it shrinks as the cover of the solution set gets finer (Fornasier's adaptive floor).

• Drift: towards the particle's own repelled consensus point, as in PCBO.
• Anisotropic noise (ANISO=1; Carrillo et al. 2019, Fornasier et al. 2021): each coordinate's amplitude is its distance to the consensus point, plus a floor. A particle that shares a find's mug pose varies mainly the spoon.
• The floor is Massimo's suggestion (a), $\sigma(|X - X_\alpha| + \kappa)$, with $\kappa$ set from the join radius: $\kappa_{\text{img}} = 0.125\,$m, where the planner joins connected pairs about half the time.

5.2 The time step (Euler–Maruyama proposal)

With $\Delta t = 0.15$ and $\xi \sim \mathcal{N}(0, I_{10})$,

6. The screen before the oracle call

A proposal is first moved along the particle's contact stratum, which needs no oracle call, and is then accepted for evaluation with probability $\min\{1, (\nu_t(x')/\varepsilon_t)^m\}$, $m = 2$: the distance $\nu_t$ of the proposed state $x'$ to everything evaluated so far, against the archive's spacing $\varepsilon_t$ (Section 3.3). If it is rejected, the noise is redrawn with the same drift, at most 50 times; after 50 rejections the last proposal is used. Only accepted moves cost oracle calls, so the screen decides where the calls go: towards states unlike any evaluated so far.

7. One iteration

for each particle i (in parallel):
    land z_i at rank r_i; if valid: one oracle call -> E_i;
        if exact: sequential flatness test (<= 8 calls) -> flat or not
add the evaluated states to the archive A_t (nearest-neighbour distances, spacing eps_t), the flat states and the flat inputs Z^flat
crowding G_i <- -log(nn_i / eps_t)
graph: kNN on images, conformal edges, all-pairs shortest paths D
weights omega_ij ~ exp(-alpha (E_j (E_i + theta) / varkappa + D_ij^2))
m_i <- z_i (+) sum_j omega_ij (z_j (-) z_i)
C <- 20 gradient steps on B(C)  (pull to m_i; energy distance to omega, the hull of Z^flat)
for each particle i:
    repeat up to 50 times:
        s <- drift + anisotropic noise (floor kappa), clipped
        rescale s so that the image moves ~ min(T_max, s ||s_raw||)
        land z_i (+) s geometrically; accept with probability min(1, (nu_t / eps_t)^m)
    z_i <- z_i (+) s

8. Parameters of the configuration

9. How the components are counted

Components are counted by checked joins: straight segments of exact equilibria, and a planner that slides along the contact strata and releases one contact at a time. Two of PCBO's flat states are in one component when a chain of such joins connects them. The joins may pass through 10,386 flat states of random landings and 8,771 stable poses of Majid's, which serve only as stepping stones and are not counted (joins within 12 cm, up to 4 attempts per pair of groups, failed joins retried in reverse, passes until nothing more joins). Two states in different components are separated at the planner's resolution; separation is not proved.

10. Open points

• Separation. Two components are distinct at the planner's resolution; no certificate proves them separate.
• The contact term at large $N$. Its value at 1,000 particles is being measured by the same run with and without it.
• Model. With $\mu = 10$ the friction cones have a half-angle of about 84°, and the mug's wall is 24 flat boxes, so a few equilibria (a spoon held by its shaft tip in a crease between two wall facets) exist only in the model.
• Mean field. With the screen and the slides along the strata, the scheme is a Markov chain whose proposal is the Euler step of Section 5; its mean-field limit is not established.

Code: code/mug_spoon/pcbo3d.cpp, landing3d.hpp, strata3d.hpp, join3d.hpp, oracle3d.hpp, geom3d.hpp. Results: code/results/mug_spoon/ (pcbo_mechanics/, flat_early/, synthesis/). Log: Docs/majid/mug_spoon_results.md.