From 4dff70b18b84f25949360e0d0c5119783f7b714f Mon Sep 17 00:00:00 2001 From: PJ Date: Sun, 16 Aug 2026 17:45:26 +0530 Subject: [PATCH] feat(analyze): add the seed-paired signed-rank comparison and record the holm family --paired contrasts two arms running the same seeds seed by seed with the wilcoxon signed-rank test rather than treating them as two independent samples, reporting the per-seed differences, the sign, a12 within pairs and the seeds usable in one arm only. --question names the family holm corrected within, and the family size is recorded next to the p-values rather than left to the reader to reconstruct. --- cmd/internal-tools/analyze/analysis.go | 66 +++- cmd/internal-tools/analyze/main.go | 18 +- cmd/internal-tools/analyze/paired.go | 276 +++++++++++++++++ cmd/internal-tools/analyze/paired_test.go | 356 ++++++++++++++++++++++ cmd/internal-tools/analyze/report.go | 28 ++ 5 files changed, 735 insertions(+), 9 deletions(-) create mode 100644 cmd/internal-tools/analyze/paired.go create mode 100644 cmd/internal-tools/analyze/paired_test.go diff --git a/cmd/internal-tools/analyze/analysis.go b/cmd/internal-tools/analyze/analysis.go index f9f26c5..dadda02 100644 --- a/cmd/internal-tools/analyze/analysis.go +++ b/cmd/internal-tools/analyze/analysis.go @@ -1,6 +1,7 @@ package main import ( + "fmt" "math" "slices" "time" @@ -51,17 +52,57 @@ type pairwiseResult struct { } type analysis struct { - GeneratedAt time.Time `json:"generated_at"` - Outcome string `json:"outcome"` - Arms []armSummary `json:"arms"` - LogRank *logRankResult `json:"log_rank"` - Pairwise []pairwiseResult `json:"pairwise"` - Notes []string `json:"notes,omitempty"` + GeneratedAt time.Time `json:"generated_at"` + Outcome string `json:"outcome"` + // Question names the family Holm corrects within. The correction is applied + // across the comparisons of one research question and never across the + // paper, so the family a p-value was adjusted in has to be recorded next to + // it rather than left to the reader to reconstruct. + Question string `json:"question,omitempty"` + HolmFamilySize int `json:"holm_family_size"` + Arms []armSummary `json:"arms"` + LogRank *logRankResult `json:"log_rank"` + Pairwise []pairwiseResult `json:"pairwise"` + Paired *pairedComparison `json:"paired,omitempty"` + Notes []string `json:"notes,omitempty"` } const outcomeDescription = "steps to first violation, right-censored at the step budget" func analyse(arms []arm, now time.Time) analysis { + result, testable := baseAnalysis(arms, now) + if len(testable) >= 2 { + result.Pairwise = comparePairs(testable) + result.HolmFamilySize = countCorrected(result.Pairwise) + } + return result +} + +// analysePaired is the seed-matched design of the actuation ablation: two arms +// running the same seeds, contrasted seed by seed rather than as two +// independent samples. +func analysePaired(arms []arm, now time.Time) (analysis, error) { + result, testable := baseAnalysis(arms, now) + if len(testable) != 2 { + return analysis{}, fmt.Errorf("a paired comparison needs exactly two arms with usable runs, found %d", len(testable)) + } + comparison, err := pairArms(testable[0], testable[1]) + if err != nil { + return analysis{}, err + } + if comparison.Pairs == 0 { + return analysis{}, fmt.Errorf("arms %q and %q share no seed with a usable run in both", + testable[0].Name, testable[1].Name) + } + if !math.IsNaN(comparison.PValue) { + comparison.HolmPValue = holm([]float64{comparison.PValue})[0] + result.HolmFamilySize = 1 + } + result.Paired = &comparison + return result, nil +} + +func baseAnalysis(arms []arm, now time.Time) (analysis, []arm) { result := analysis{GeneratedAt: now, Outcome: outcomeDescription} for _, current := range arms { result.Arms = append(result.Arms, summarize(current)) @@ -86,9 +127,18 @@ func analyse(arms []arm, now time.Time) analysis { } test := logRank(names, groups) result.LogRank = &test - result.Pairwise = comparePairs(testable) } - return result + return result, testable +} + +func countCorrected(pairs []pairwiseResult) int { + corrected := 0 + for _, pair := range pairs { + if !math.IsNaN(pair.PValue) { + corrected++ + } + } + return corrected } func comparePairs(arms []arm) []pairwiseResult { diff --git a/cmd/internal-tools/analyze/main.go b/cmd/internal-tools/analyze/main.go index 96d980a..47e8d2d 100644 --- a/cmd/internal-tools/analyze/main.go +++ b/cmd/internal-tools/analyze/main.go @@ -25,6 +25,9 @@ Usage: Each directory is one produced by the campaign tool and must hold campaign.json and runs.jsonl. Directories sharing an arm label are pooled and must agree on the step budget. + +One invocation is one research question: Holm corrects across the comparisons it +produces and across nothing else. ` type stringList []string @@ -48,8 +51,12 @@ func run(arguments []string, stdout, stderr io.Writer) error { } var directories stringList var jsonPath string + var question string + var paired bool flagSet.Var(&directories, "campaign", "campaign directory to read; repeat for more, or pass them as arguments") flagSet.StringVar(&jsonPath, "json", "", "write the machine-readable summary here, or - for stdout") + flagSet.StringVar(&question, "question", "", "the research question these campaigns answer; Holm corrects within one invocation, and this records which family that was") + flagSet.BoolVar(&paired, "paired", false, "the two arms ran the same seeds: contrast them seed by seed with the Wilcoxon signed-rank test instead of the rank-sum test") if err := flagSet.Parse(arguments); err != nil { return err } @@ -73,7 +80,16 @@ func run(arguments []string, stdout, stderr io.Writer) error { if err != nil { return err } - result := analyse(arms, time.Now().UTC()) + var result analysis + if paired { + result, err = analysePaired(arms, time.Now().UTC()) + if err != nil { + return err + } + } else { + result = analyse(arms, time.Now().UTC()) + } + result.Question = question writeReport(result, stdout) if jsonPath == "" { diff --git a/cmd/internal-tools/analyze/paired.go b/cmd/internal-tools/analyze/paired.go new file mode 100644 index 0000000..792ed20 --- /dev/null +++ b/cmd/internal-tools/analyze/paired.go @@ -0,0 +1,276 @@ +package main + +import ( + "fmt" + "math" + "slices" +) + +type signedRankResult struct { + Pairs int `json:"pairs"` + NonZero int `json:"non_zero_pairs"` + Statistic float64 `json:"signed_rank_v"` + PValue float64 `json:"p_value"` + Exact bool `json:"exact"` +} + +// exactSignedRankLimit matches R's wilcox.test: the exact null distribution is +// used only below this many non-zero differences, and only when nothing is tied. +const exactSignedRankLimit = 50 + +// signedRank is the two-sided Wilcoxon signed-rank test over paired +// differences. Zero differences are dropped before ranking and the statistic is +// the sum of the ranks carried by the positive differences, which is the +// quantity R's wilcox.test calls V. Wilcoxon (1945), "Individual Comparisons by +// Ranking Methods", Biometrics Bulletin 1(6), 80-83. +func signedRank(differences []float64) signedRankResult { + result := signedRankResult{ + Pairs: len(differences), + Statistic: math.NaN(), + PValue: math.NaN(), + } + var magnitudes []float64 + var positive []bool + for _, difference := range differences { + if difference == 0 { + continue + } + magnitudes = append(magnitudes, math.Abs(difference)) + positive = append(positive, difference > 0) + } + result.NonZero = len(magnitudes) + if result.NonZero == 0 { + return result + } + ranks, tieGroups := midRanks(magnitudes) + statistic := 0.0 + for index, rank := range ranks { + if positive[index] { + statistic += rank + } + } + result.Statistic = statistic + + droppedZeros := len(differences) != result.NonZero + if len(tieGroups) == 0 && !droppedZeros && result.NonZero < exactSignedRankLimit { + result.Exact = true + result.PValue = exactSignedRankTwoSided(statistic, result.NonZero) + return result + } + result.PValue = normalSignedRankTwoSided(statistic, result.NonZero, tieGroups) + return result +} + +// normalSignedRankTwoSided follows the large-sample branch of R's wilcox.test: +// +// mean = n(n+1)/4 +// variance = n(n+1)(2n+1)/24 - sum(t^3 - t)/48 +// +// where t runs over the sizes of the groups tied on the absolute difference. +// The 0.5 shift toward the null mean is the continuity correction. +func normalSignedRankTwoSided(statistic float64, count int, tieGroups []int) float64 { + variance := signedRankVariance(count, tieGroups) + if variance <= 0 { + return 1 + } + size := float64(count) + centered := statistic - size*(size+1)/4 + correction := 0.0 + switch { + case centered > 0: + correction = 0.5 + case centered < 0: + correction = -0.5 + } + z := (centered - correction) / math.Sqrt(variance) + tail := math.Min(standardNormalUpperTail(z), standardNormalUpperTail(-z)) + return math.Min(2*tail, 1) +} + +func signedRankVariance(count int, tieGroups []int) float64 { + size := float64(count) + tieAdjustment := 0.0 + for _, group := range tieGroups { + tied := float64(group) + tieAdjustment += tied*tied*tied - tied + } + return size*(size+1)*(2*size+1)/24 - tieAdjustment/48 +} + +// exactSignedRankTwoSided doubles the smaller exact tail, as R's wilcox.test +// does. +func exactSignedRankTwoSided(statistic float64, count int) float64 { + if statistic > float64(count)*float64(count+1)/4 { + return math.Min(2*exactSignedRankUpperTail(statistic, count), 1) + } + return math.Min(2*exactSignedRankLowerTail(statistic, count), 1) +} + +// exactSignedRankUpperTail is P(V >= statistic) under the null with no ties. +func exactSignedRankUpperTail(statistic float64, count int) float64 { + counts := exactSignedRankCounts(count) + total, tail := 0.0, 0.0 + for value, weight := range counts { + total += weight + if float64(value) >= statistic { + tail += weight + } + } + return tail / total +} + +func exactSignedRankLowerTail(statistic float64, count int) float64 { + counts := exactSignedRankCounts(count) + total, tail := 0.0, 0.0 + for value, weight := range counts { + total += weight + if float64(value) <= statistic { + tail += weight + } + } + return tail / total +} + +// exactSignedRankCounts returns the number of sign assignments producing each +// value of V from 0 to n(n+1)/2. V is the sum of the ranks held by the positive +// differences, so the count is a subset-sum tally over the ranks 1 to n. +func exactSignedRankCounts(count int) []float64 { + high := count * (count + 1) / 2 + table := make([]float64, high+1) + table[0] = 1 + for rank := 1; rank <= count; rank++ { + for sum := high; sum >= rank; sum-- { + if table[sum-rank] != 0 { + table[sum] += table[sum-rank] + } + } + } + return table +} + +// pairedComparison is the seed-matched contrast the actuation ablation reports. +// The difference is the first arm's steps to first violation less the second's, +// so a positive median means the second arm reached its first violation sooner, +// and Sign carries that direction as a number the decision rule can read. +type pairedComparison struct { + First string `json:"first"` + Second string `json:"second"` + Pairs int `json:"pairs"` + UnpairedSeeds []int64 `json:"unpaired_seeds,omitempty"` + MedianDifference float64 `json:"median_step_difference"` + Sign int `json:"sign"` + FirstSooner int `json:"first_sooner"` + SecondSooner int `json:"second_sooner"` + Tied int `json:"tied"` + // A12 is the within-pair form of the Vargha-Delaney effect size, the share + // of matched seeds on which the first arm took more steps, counting a tie as + // half. A matched design has no reason to compare the two arms as pooled + // bags of runs when each seed has a partner. + A12 float64 `json:"a12_within_pairs"` + Statistic float64 `json:"signed_rank_v"` + PValue float64 `json:"p_value"` + HolmPValue float64 `json:"holm_p_value"` + Exact bool `json:"exact"` +} + +// pairArms matches the two arms by seed and contrasts them pair by pair. +// Censored runs enter at the step budget, the same convention the unpaired +// comparison uses. A seed usable in one arm and not the other is named rather +// than dropped silently, because that is a host that lost a run and it is what +// the campaign manifest exists to make visible. +func pairArms(first, second arm) (pairedComparison, error) { + firstBySeed, err := usableBySeed(first) + if err != nil { + return pairedComparison{}, err + } + secondBySeed, err := usableBySeed(second) + if err != nil { + return pairedComparison{}, err + } + + comparison := pairedComparison{ + First: first.Name, + Second: second.Name, + A12: math.NaN(), + Statistic: math.NaN(), + PValue: math.NaN(), + HolmPValue: math.NaN(), + } + var differences []float64 + for _, seed := range sortedSeeds(firstBySeed, secondBySeed) { + left, inFirst := firstBySeed[seed] + right, inSecond := secondBySeed[seed] + if !inFirst || !inSecond { + comparison.UnpairedSeeds = append(comparison.UnpairedSeeds, seed) + continue + } + difference := observationOf(left, first.Budget).Steps - observationOf(right, second.Budget).Steps + differences = append(differences, difference) + switch { + case difference < 0: + comparison.FirstSooner++ + case difference > 0: + comparison.SecondSooner++ + default: + comparison.Tied++ + } + } + comparison.Pairs = len(differences) + if comparison.Pairs == 0 { + return comparison, nil + } + + comparison.MedianDifference = medianOf(differences) + switch { + case comparison.MedianDifference > 0: + comparison.Sign = 1 + case comparison.MedianDifference < 0: + comparison.Sign = -1 + } + comparison.A12 = (float64(comparison.SecondSooner) + 0.5*float64(comparison.Tied)) / float64(comparison.Pairs) + test := signedRank(differences) + comparison.Statistic = test.Statistic + comparison.PValue = test.PValue + comparison.Exact = test.Exact + return comparison, nil +} + +func usableBySeed(current arm) (map[int64]classifiedRun, error) { + bySeed := map[int64]classifiedRun{} + for _, item := range current.Runs { + if item.ExcludedBecause != "" { + continue + } + if _, repeated := bySeed[item.Seed]; repeated { + return nil, fmt.Errorf("arm %q has more than one usable run for seed %d: a seed-matched "+ + "comparison cannot choose between them", current.Name, item.Seed) + } + bySeed[item.Seed] = item + } + return bySeed, nil +} + +func sortedSeeds(sets ...map[int64]classifiedRun) []int64 { + var seeds []int64 + seen := map[int64]bool{} + for _, set := range sets { + for seed := range set { + if seen[seed] { + continue + } + seen[seed] = true + seeds = append(seeds, seed) + } + } + slices.Sort(seeds) + return seeds +} + +func medianOf(values []float64) float64 { + sorted := slices.Sorted(slices.Values(values)) + middle := len(sorted) / 2 + if len(sorted)%2 == 1 { + return sorted[middle] + } + return (sorted[middle-1] + sorted[middle]) / 2 +} diff --git a/cmd/internal-tools/analyze/paired_test.go b/cmd/internal-tools/analyze/paired_test.go new file mode 100644 index 0000000..56d8ac6 --- /dev/null +++ b/cmd/internal-tools/analyze/paired_test.go @@ -0,0 +1,356 @@ +package main + +import ( + "math" + "testing" +) + +// Hollander and Wolfe (1973), 29f: Hamilton depression scale factor +// measurements on nine patients, first at admission and again after tranquilizer +// treatment. R's wilcox.test help page uses exactly these vectors as its paired +// example and reports +// +// wilcox.test(x, y, paired = TRUE, alternative = "greater") +// ## V = 40, p-value = 0.01953 +var ( + depressionAtAdmission = []float64{1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30} + depressionAfterOneWeek = []float64{0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29} +) + +func differencesOf(first, second []float64) []float64 { + differences := make([]float64, len(first)) + for index := range first { + differences[index] = first[index] - second[index] + } + return differences +} + +func TestSignedRank_MatchesPublishedDepressionResult(t *testing.T) { + result := signedRank(differencesOf(depressionAtAdmission, depressionAfterOneWeek)) + + if result.Statistic != 40 { + t.Errorf("statistic %v, want 40", result.Statistic) + } + if !result.Exact { + t.Error("expected the exact null distribution for nine untied differences") + } + upper := exactSignedRankUpperTail(40, 9) + if math.Abs(upper-0.01953) > 5e-6 { + t.Errorf("one-sided p-value %.6f, want 0.01953", upper) + } + if math.Abs(result.PValue-2*0.01953125) > 1e-9 { + t.Errorf("two-sided p-value %.6f, want %.6f", result.PValue, 2*0.01953125) + } +} + +// Reversing the pairs mirrors the statistic about n(n+1)/2 and leaves the +// two-sided p-value alone, which R reports as V = 5 on the same data. +func TestSignedRank_ReversedPairsMirrorTheStatistic(t *testing.T) { + forward := signedRank(differencesOf(depressionAtAdmission, depressionAfterOneWeek)) + reversed := signedRank(differencesOf(depressionAfterOneWeek, depressionAtAdmission)) + if reversed.Statistic != 5 { + t.Errorf("reversed statistic %v, want 5", reversed.Statistic) + } + if math.Abs(reversed.PValue-forward.PValue) > 1e-12 { + t.Errorf("reversed p-value %v, want %v", reversed.PValue, forward.PValue) + } +} + +// The exact null distribution must be a proper distribution: 2^n sign +// assignments in total, symmetric about n(n+1)/4. +func TestExactSignedRankCounts_FormAProperSymmetricDistribution(t *testing.T) { + counts := exactSignedRankCounts(8) + total := 0.0 + for _, count := range counts { + total += count + } + if total != 256 { + t.Errorf("counts sum to %v, want 2^8 = 256", total) + } + for index := range counts { + if counts[index] != counts[len(counts)-1-index] { + t.Errorf("count at %d is %v but %v at the mirrored point", index, counts[index], counts[len(counts)-1-index]) + } + } +} + +// The tie-corrected variance is checked against the exact permutation variance +// of the statistic, computed here by enumerating every sign assignment over the +// observed midranks. That is an independent calculation rather than a second +// call into the implementation under test. +func TestSignedRankVariance_MatchesExactPermutationVariance(t *testing.T) { + cases := [][]float64{ + {1, -2, 3, -4, 5, 6, -7, 8}, + {12, -12, 12, 12, -5, 5, 30, -30}, + {-400, 400, 400, -400, 400, 400, 400, 400}, + {3, 3, 3, -3, -3, 7, 7, 9, 9}, + } + for _, differences := range cases { + magnitudes := make([]float64, len(differences)) + for index, difference := range differences { + magnitudes[index] = math.Abs(difference) + } + ranks, tieGroups := midRanks(magnitudes) + mean, variance := permutationMomentsOfSignedRank(ranks) + size := float64(len(ranks)) + if expected := size * (size + 1) / 4; math.Abs(mean-expected) > 1e-9 { + t.Errorf("permutation mean %v for %v, want %v", mean, differences, expected) + } + if got := signedRankVariance(len(ranks), tieGroups); math.Abs(got-variance) > 1e-9 { + t.Errorf("variance %v for %v, want the permutation variance %v", got, differences, variance) + } + } +} + +// With ties present the normal approximation is the only branch available, so +// it is checked against the exact permutation p-value of the same statistic on +// the same data. +func TestSignedRank_TiedDifferencesTrackTheExactPermutationPValue(t *testing.T) { + cases := [][]float64{ + {40, 40, 40, -12, 33, 40, 40, -3, 40, 21, 40, 40}, + {-5, -5, -5, -5, 9, 9, 2, 2, -1, -1, 40, 40}, + {40, 40, -12, 33, 40, -3, 21, 40, 15, -9}, + } + for _, differences := range cases { + result := signedRank(differences) + if result.Exact { + t.Errorf("%v used the exact null distribution despite ties", differences) + } + exact := permutationSignedRankTwoSided(differences) + if math.Abs(result.PValue-exact) > 0.03 { + t.Errorf("normal approximation p %.4f for %v, want near the permutation p %.4f", + result.PValue, differences, exact) + } + } +} + +// Every difference the same size is the degenerate end of the tie correction, +// and the normal approximation is genuinely poor there: the exact randomization +// p-value on this sample is 0.3438 against the approximation's 0.2273. The tool +// keeps R's formula rather than the randomization p-value so that a reviewer +// running wilcox.test on the same differences reads the same number, and the +// expected value here is that published formula worked through by hand: +// +// n = 10, one tied group of 10, V = 7 * 5.5 = 38.5 +// mean = 10 * 11 / 4 = 27.5 +// variance = 10 * 11 * 21 / 24 - (10^3 - 10) / 48 = 96.25 - 20.625 = 75.625 +// z = (38.5 - 27.5 - 0.5) / sqrt(75.625) +func TestSignedRank_EveryDifferenceTheSameSizeFollowsTheDocumentedFormula(t *testing.T) { + differences := []float64{7, 7, 7, 7, 7, 7, 7, -7, -7, -7} + result := signedRank(differences) + + if result.Statistic != 38.5 { + t.Errorf("statistic %v, want 38.5", result.Statistic) + } + if got := signedRankVariance(10, []int{10}); math.Abs(got-75.625) > 1e-12 { + t.Errorf("variance %v, want 75.625", got) + } + expected := 2 * standardNormalUpperTail(10.5/math.Sqrt(75.625)) + if math.Abs(result.PValue-expected) > 1e-12 { + t.Errorf("p-value %v, want %v", result.PValue, expected) + } + if randomization := permutationSignedRankTwoSided(differences); math.Abs(randomization-0.3438) > 5e-4 { + t.Errorf("randomization p-value %.4f, want 0.3438", randomization) + } +} + +// R drops zero differences before ranking and tests what is left, so a pair +// where both arms took the same number of steps carries no direction and must +// not be ranked as though it did. +func TestSignedRank_ZeroDifferencesAreDropped(t *testing.T) { + result := signedRank([]float64{0, 0, 3, -1, 2}) + if result.Pairs != 5 || result.NonZero != 3 { + t.Errorf("pairs %d non-zero %d, want 5 and 3", result.Pairs, result.NonZero) + } + // Ranks over |{3, 1, 2}| are 3, 1, 2, and the positive differences hold 3 + // and 2. + if result.Statistic != 5 { + t.Errorf("statistic %v, want 5", result.Statistic) + } + if result.Exact { + t.Error("used the exact null distribution despite dropped zeros") + } +} + +func TestSignedRank_EveryDifferenceZero(t *testing.T) { + result := signedRank([]float64{0, 0, 0}) + if result.NonZero != 0 { + t.Errorf("non-zero pairs %d, want 0", result.NonZero) + } + if !math.IsNaN(result.PValue) || !math.IsNaN(result.Statistic) { + t.Errorf("result %+v, want everything undefined", result) + } +} + +// permutationMomentsOfSignedRank enumerates every sign assignment and returns +// the mean and variance of the statistic over them. +func permutationMomentsOfSignedRank(ranks []float64) (float64, float64) { + values := signedRankPermutationValues(ranks) + mean := 0.0 + for _, value := range values { + mean += value + } + mean /= float64(len(values)) + variance := 0.0 + for _, value := range values { + variance += (value - mean) * (value - mean) + } + return mean, variance / float64(len(values)) +} + +// permutationSignedRankTwoSided is the exact randomization p-value: the +// proportion of sign assignments whose statistic is at least as far from the +// null mean as the observed one. +func permutationSignedRankTwoSided(differences []float64) float64 { + magnitudes := make([]float64, 0, len(differences)) + observed := 0.0 + for _, difference := range differences { + if difference == 0 { + continue + } + magnitudes = append(magnitudes, math.Abs(difference)) + } + ranks, _ := midRanks(magnitudes) + position := 0 + for _, difference := range differences { + if difference == 0 { + continue + } + if difference > 0 { + observed += ranks[position] + } + position++ + } + size := float64(len(ranks)) + mean := size * (size + 1) / 4 + values := signedRankPermutationValues(ranks) + extreme := 0 + for _, value := range values { + if math.Abs(value-mean) >= math.Abs(observed-mean)-1e-9 { + extreme++ + } + } + return float64(extreme) / float64(len(values)) +} + +func signedRankPermutationValues(ranks []float64) []float64 { + values := make([]float64, 0, 1< 1e-12 { + t.Errorf("a12 within pairs %v, want %v", comparison.A12, 2.5/3) + } +} + +// A run excluded as missing data cannot be paired against anything, and the +// seed it came from has to be named rather than silently shrinking the sample. +func TestPairArms_NamesSeedsUsableInOneArmOnly(t *testing.T) { + pre := arm{Name: "pre", Budget: 40, Runs: []classifiedRun{ + violatingRun(1, 30, 30, "p"), + {Seed: 2, ExcludedBecause: reasonTimedOut}, + violatingRun(3, 20, 20, "p"), + }} + post := arm{Name: "post", Budget: 40, Runs: []classifiedRun{ + violatingRun(1, 10, 10, "p"), + violatingRun(2, 11, 11, "p"), + }} + + comparison, err := pairArms(pre, post) + if err != nil { + t.Fatal(err) + } + if comparison.Pairs != 1 { + t.Fatalf("%d pairs, want 1", comparison.Pairs) + } + if len(comparison.UnpairedSeeds) != 2 || comparison.UnpairedSeeds[0] != 2 || comparison.UnpairedSeeds[1] != 3 { + t.Errorf("unpaired seeds %v, want [2 3]", comparison.UnpairedSeeds) + } +} + +func TestPairArms_RefusesTwoUsableRunsForOneSeed(t *testing.T) { + pooled := arm{Name: "pre", Budget: 40, Runs: []classifiedRun{ + violatingRun(1, 30, 30, "p"), + violatingRun(1, 12, 12, "p"), + }} + post := arm{Name: "post", Budget: 40, Runs: []classifiedRun{violatingRun(1, 10, 10, "p")}} + + if _, err := pairArms(pooled, post); err == nil { + t.Fatal("paired two arms where one seed ran twice") + } +} + +// The paired comparison is the ablation's decision rule, so the direction it +// reports has to survive the arms being passed the other way round. +func TestPairArms_DirectionReversesWithTheArms(t *testing.T) { + pre := arm{Name: "pre", Budget: 400, Runs: []classifiedRun{ + cleanRun(1, 400), cleanRun(2, 400), violatingRun(3, 380, 380, "p"), + cleanRun(4, 400), violatingRun(5, 350, 350, "p"), + }} + post := arm{Name: "post", Budget: 400, Runs: []classifiedRun{ + violatingRun(1, 40, 40, "p"), violatingRun(2, 90, 90, "p"), violatingRun(3, 60, 60, "p"), + violatingRun(4, 120, 120, "p"), violatingRun(5, 30, 30, "p"), + }} + + forward, err := pairArms(pre, post) + if err != nil { + t.Fatal(err) + } + reversed, err := pairArms(post, pre) + if err != nil { + t.Fatal(err) + } + if forward.Sign != 1 || reversed.Sign != -1 { + t.Errorf("signs %+d and %+d, want +1 then -1", forward.Sign, reversed.Sign) + } + if forward.MedianDifference != -reversed.MedianDifference { + t.Errorf("median differences %v and %v, want opposites", forward.MedianDifference, reversed.MedianDifference) + } + if math.Abs(forward.PValue-reversed.PValue) > 1e-12 { + t.Errorf("p-values %v and %v, want the same two-sided value", forward.PValue, reversed.PValue) + } + if math.Abs(forward.A12+reversed.A12-1) > 1e-12 { + t.Errorf("a12 %v and %v, want them to sum to 1", forward.A12, reversed.A12) + } +} diff --git a/cmd/internal-tools/analyze/report.go b/cmd/internal-tools/analyze/report.go index 63bcff3..f9f8d75 100644 --- a/cmd/internal-tools/analyze/report.go +++ b/cmd/internal-tools/analyze/report.go @@ -94,6 +94,10 @@ func writeReport(result analysis, out io.Writer) { }) } + if result.Paired != nil { + writePaired(out, *result.Paired) + } + if len(result.Pairwise) > 0 { fmt.Fprintln(out, "\npairwise wilcoxon rank-sum, censored runs held at the budget") fmt.Fprintln(out, "a12 above 0.5 means the first arm takes more steps to its first violation") @@ -112,11 +116,35 @@ func writeReport(result analysis, out io.Writer) { }) } + if result.HolmFamilySize > 0 { + family := "this invocation" + if result.Question != "" { + family = result.Question + } + fmt.Fprintf(out, "\nholm correction applied within %s, over %d comparison(s)\n", family, result.HolmFamilySize) + } + for _, note := range result.Notes { fmt.Fprintf(out, "\nnote: %s\n", note) } } +func writePaired(out io.Writer, comparison pairedComparison) { + fmt.Fprintf(out, "\npaired per-seed difference, %s minus %s, censored runs held at the budget\n", + comparison.First, comparison.Second) + fmt.Fprintf(out, "%d seed pair(s): %s sooner in %d, %s sooner in %d, tied in %d\n", + comparison.Pairs, comparison.First, comparison.FirstSooner, + comparison.Second, comparison.SecondSooner, comparison.Tied) + fmt.Fprintf(out, "median difference %+.1f steps, sign %+d, a12 within pairs %.3f\n", + comparison.MedianDifference, comparison.Sign, comparison.A12) + fmt.Fprintf(out, "wilcoxon signed-rank v %.1f, p %s, holm p %s\n", + comparison.Statistic, formatPValue(comparison.PValue), formatPValue(comparison.HolmPValue)) + if len(comparison.UnpairedSeeds) > 0 { + fmt.Fprintf(out, "%d seed(s) usable in one arm only and left out of the pairing: %v\n", + len(comparison.UnpairedSeeds), comparison.UnpairedSeeds) + } +} + func sortedKeys(counts map[string]int) []string { return slices.Sorted(maps.Keys(counts)) }