phase_a_layout.py

#!/usr/bin/env python3
"""
phase_a_layout.py -- lay out a 16-patch Phase A test plot so that yard gradients cancel.
Run: python3 phase_a_layout.py

Four on/off factors give 2^4 = 16 recipes, one patch each, in a 4x4 grid (rows run upslope to
downslope, columns west to east). Recipe numbers 0-15 are placed by a 4x4 magic square: a recipe's four
bits say which factors are on. In every 4x4 magic square each bit is on in exactly two cells of every
row and every column (Fitting 1931, Jahresbericht DMV 40, p. 190), so slope, sun and drainage
differences between rows or between columns cancel out of every single-factor comparison.

The square used is affine: each bit is an XOR of row and column bits. Then every other effect
(interactions too) is either balanced in all rows and columns or tied entirely to the rows or to the
columns. That is Yates's (1937) confounding of a factorial in a row-column layout, the "quasi-Latin
square", here at its smallest size. The magic-square property itself is more than the design needs:
row and column balance is what cancels gradients.

The factors are the founder's choice: edit FACTORS. The defaults come from the 2026-10-22 handoff.
The layout suits the plot tests in MISSION §8: the till-in-place patches (Route A) and the
floating-shoe test bed. It does not replace the lab cube tests.
"""
import itertools

FACTORS = [("Roman hot-mix lime", "none", "hot-mixed quicklime"),     # (name, level off, level on)
           ("rice husk ash", "none", "husk ash"),
           ("grid face", "top", "bottom"),
           ("cure", "method A", "method B")]
KEEP = (0, 1)        # the interaction that must stay clear of rows and columns: lime x husk ash
                     # (lime + pozzolan is the binder; each alone is not)

# An affine square: value bit k = (r_k . row bits) XOR (c_k . column bits) XOR constant. Bit k is
# balanced on rows iff c_k != 0, on columns iff r_k != 0, on both diagonals iff r_k != c_k.
FORMS = [(r, c) for r in range(1, 4) for c in range(1, 4) if r != c]
def dot(a, b):
    return bin(a & b).count("1") & 1
def square(forms, const):
    return [[sum((dot(r, i) ^ dot(c, j) ^ ((const >> k) & 1)) << k for k, (r, c) in enumerate(forms))
             for j in range(4)] for i in range(4)]
def effect_form(forms, S):                   # an effect (set of factors) is the XOR of their forms
    r = c = 0
    for k in S:
        r ^= forms[k][0]; c ^= forms[k][1]
    return r, c
def status(forms, S):
    r, c = effect_form(forms, S)
    return "balanced" if r and c else "tied to rows" if r else "tied to columns"

best = None
for forms in itertools.permutations(FORMS, 4):
    span = {(0, 0)}
    for f in forms:
        span |= {(a ^ f[0], b ^ f[1]) for a, b in span}
    if len(span) < 16 or status(forms, KEEP) != "balanced":
        continue                               # need four independent forms (16 distinct recipes)
    score = sum(status(forms, S) == "balanced" for S in itertools.combinations(range(4), 2))
    if best is None or score > best[0]:
        best = (score, forms)
forms = best[1]
M = square(forms, 0)
lines = M + [list(c) for c in zip(*M)] + [[M[i][i] for i in range(4)], [M[i][3 - i] for i in range(4)]]
assert sorted(v for r in M for v in r) == list(range(16)) and all(sum(l) == 30 for l in lines)
assert all(sum((v >> k) & 1 for v in l) == 2 for l in lines for k in range(4))

print("PHASE A TEST PLOT: 16 patches, one per recipe (a 4x4 magic square)")
print("rows run upslope -> downslope, columns west -> east; each cell is a recipe number\n")
for r in M:
    print("   " + "  ".join(f"{v:2d}" for v in r))
print("\nrecipe = 4 bits:  " + ";  ".join(f"bit {k} = {n} ({off} / {on})" for k, (n, off, on) in enumerate(FACTORS)))
print("\npatches:")
for i, r in enumerate(M):
    for j, v in enumerate(r):
        print(f"   row {i+1} col {j+1}  #{v:2d}  " +
              ", ".join(FACTORS[k][2] if (v >> k) & 1 else FACTORS[k][1] for k in range(4)))
print("\nevery factor is 'on' in exactly 2 patches of every row, column and diagonal: True")

print("\nall 15 effects (main effects, then interactions):")
lost = 0
for n in range(1, 5):
    for S in itertools.combinations(range(4), n):
        st = status(forms, S)
        lost += st != "balanced"
        name = " x ".join(FACTORS[k][0] for k in S)
        note = {"balanced": "estimable, free of row and column effects",
                "tied to rows": "confounded with the slope (rows)",
                "tied to columns": "confounded with east-west (columns)"}[st]
        print(f"   {name}: {note}")
free = 15 - lost - 4 - sum(status(forms, S) == "balanced" for S in itertools.combinations(range(4), 2))
print(f"\n{lost} of 15 effects are tied to rows or columns (the price of removing both gradients).")
print(f"Analysis: estimate each effect as mean(on) - mean(off). With one replicate only {free} degree(s)")
print("of freedom are left for error, and only if 3- and 4-factor interactions are negligible. Judge")
print("the effects on a half-normal plot (Daniel 1959), or lay a second plot with another square to")
print("get a real error estimate. If only the slope matters, block on rows alone: 3 effects are lost, not 6.")