Measurement robustness

cmd/rtcompare-example/main.go

@@ -1,96 +1,151 @@
+// Command rtcompare-example compares two median implementations and reports
+// whether the difference between them is one this machine can actually resolve.
+//
+// The comparison itself is the least interesting part. What the example is
+// really about is everything around it: sizing the batches so that a difference
+// below the clock's resolution is measurable at all, finding out what the
+// harness invents on its own before trusting it, and noticing that the two
+// candidates are not equally well behaved. rtcompare.Compare does all of that in
+// one call and reports what it decided, which is what the first half of this
+// program shows. The second half takes the same measurements apart by hand, for
+// when the summary is not enough.
 package main
 
 import (
 	"fmt"
-	"runtime"
+	"os"
 
 	"github.com/TomTonic/rtcompare"
 )
 
+// sink absorbs the results of measured work. Assigning to a package-level
+// variable keeps a candidate honest; it is cheap insurance rather than a
+// necessity, since the Go compiler does not delete a loop merely because
+// nothing reads what it computes.
+var sink float64
+
+const arraySize = 50
+
 func main() {
-	const (
-		N              = 50     // size of input array
-		repeats        = 101    // number of timing samples
-		innerLoops     = 2000   // number of median calls per timing sample
-		precisionLevel = 10_000 // bootstrap repetitions
-	)
-
-	rng := rtcompare.NewDPRNG()
-
-	// Initialitze two working arrays
-	workArrayMedian := make([]float64, N)
-	safeState := rng.State
-	fillArray(workArrayMedian, rng) // rng is passed by value here so we should not need to safeguard its state
-	if safeState != rng.State {
-		panic("rng state was modified unexpectedly")
+	// Both candidates refresh their input inside the loop, because QuickMedian
+	// mutates what it is given and the two must do the same work per operation
+	// to be comparable. That refresh is measured along with the candidate and
+	// dilutes the difference between them; see the note on attenuation in the
+	// Collect documentation. It cannot be hoisted into Setup, because it has to
+	// happen per operation rather than per batch.
+	quick := medianCandidate("QuickMedian", rtcompare.QuickMedian)
+	sorting := medianCandidate("Median", rtcompare.Median)
+
+	// Everything is left at its default: the batches are sized so the clock
+	// contributes at most a tenth of a percent, both candidates are validated
+	// against themselves, the measurement order is interleaved, and the
+	// resampling scheme is chosen from the dependence actually observed.
+	fmt.Println("Comparing QuickMedian against Median. This validates the")
+	fmt.Println("harness against each candidate first, so it takes a few seconds.")
+
+	report, err := rtcompare.Compare(quick, sorting, rtcompare.CompareOptions{
+		Collect:    rtcompare.CollectOptions{GCBetween: true},
+		Thresholds: []float64{0.05, 0.10, 0.20, 0.30},
+	})
+	if err != nil {
+		fail("comparing: %v", err)
 	}
-	workArrayQuick := make([]float64, N)
-	fillArray(workArrayQuick, rng)
-
-	// Warm-up both methods
-	_ = rtcompare.Median(workArrayMedian)
-	_ = rtcompare.QuickMedian(workArrayQuick)
-
-	// Collect timing samples
-	var timesMedian []float64
-	var timesQuick []float64
-
-	for range repeats {
-		// Set rng to a new state for each timing sample
-		rng = rtcompare.NewDPRNG()
-
-		// make sure to avoid GC noise
-		runtime.GC()
-
-		// Measure Median
-		t1 := rtcompare.SampleTime()
-		// we neet to measure multiple iterations of the function to make sure the time measurement
-		// is not polluted by the timer's resolution too much (quantization noise)
-		for range innerLoops {
-			// Refresh the data in the working array - this function has constant runtime.
-			// Even though Median does not mutate its input we need to do this for the results to be comparable.
-			fillArray(workArrayMedian, rng)
-			_ = rtcompare.Median(workArrayMedian)
-		}
-		t2 := rtcompare.SampleTime()
-		durMedian := float64(rtcompare.DiffTimeStamps(t1, t2)) / float64(innerLoops)
-		timesMedian = append(timesMedian, durMedian)
-
-		// the Median function allocates memory, so we trigger a GC cycle again to reduce noise
-		runtime.GC()
-
-		// Measure QuickMedian
-		t3 := rtcompare.SampleTime()
-		for range innerLoops {
-			// Refresh the data in the working array - this function has constant runtime.
-			// This is necessary as QuickMedian mutates its input. On the other hand, it does not allocate extra memory.
-			fillArray(workArrayQuick, rng)
-			_ = rtcompare.QuickMedian(workArrayQuick)
+
+	fmt.Printf("\n%s\n", report)
+
+	if !report.Resolved {
+		fmt.Println("\nNothing has been resolved; the run says only that any difference is small.")
+		return
+	}
+
+	// The report carries the evidence as well as the verdict, so the pieces are
+	// there when the summary is not enough.
+	fmt.Println("\n--- the evidence behind that verdict ---")
+
+	// How large is the difference, and how precisely is that known? This is the
+	// question to ask when no threshold was given to you.
+	e := report.Estimate
+	fmt.Printf("\nEstimated difference: %s\n", e)
+	fmt.Printf("  point estimate %+.2f%%, bootstrap median %+.2f%% (a large gap would mean the\n"+
+		"  statistic behaves awkwardly on this data and the interval deserves suspicion)\n",
+		e.Delta*100, e.BootstrapMedian*100)
+	fmt.Printf("  the interval %s zero, so a difference %s established\n",
+		yesNo(e.Excludes(0), "excludes", "includes"),
+		yesNo(e.Excludes(0), "has been", "has not been"))
+
+	// What the harness does to identical code, which is what the result above
+	// has to be read against.
+	fmt.Printf("\nNoise floor across both candidates: %.2f%%\n", report.NoiseFloor*100)
+	for _, v := range []struct {
+		name string
+		val  rtcompare.HarnessValidation
+	}{{"QuickMedian", report.ValidationA}, {"Median", report.ValidationB}} {
+		fmt.Printf("\n%s:\n%s\n", v.name, v.val)
+	}
+
+	// A trend across a run is invisible to the bootstrap, which treats the
+	// samples as an unordered bag, so Compare tests for it separately.
+	for _, d := range []struct {
+		name string
+		rep  rtcompare.DriftReport
+	}{{"QuickMedian", report.DriftA}, {"Median", report.DriftB}} {
+		if d.rep.N > 0 && d.rep.Drifted(0.05) {
+			fmt.Printf("\nnote: %s %s\n", d.name, d.rep)
 		}
-		t4 := rtcompare.SampleTime()
-		durQuick := float64(rtcompare.DiffTimeStamps(t3, t4)) / float64(innerLoops)
-		timesQuick = append(timesQuick, durQuick)
 	}
 
-	// Compare the timing distributions using bootstrap
-	speedups := []float64{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0} // relative speedups to test
-	results, err := rtcompare.CompareSamples(timesQuick, timesMedian, speedups, precisionLevel)
-	if err != nil {
-		panic(err)
+	if report.BlockLength > 1 {
+		fmt.Printf("\nThe samples were correlated with their neighbours (%+.3f), so they were\n"+
+			"resampled in blocks of %d rather than one at a time. Median allocates a copy on\n"+
+			"every call, which tends to put it there.\n", report.Autocorrelation, report.BlockLength)
+	}
+
+	fmt.Println("\nConfidence that QuickMedian beats Median by at least:")
+	for _, t := range []float64{0.05, 0.10, 0.20, 0.30} {
+		fmt.Printf("  %6.2f%%   %7.2f%%\n", t*100, report.Confidence[t]*100)
 	}
 
-	// Report results
-	fmt.Println("⏱️ Runtime comparison: QuickMedian vs. Median for arrays of size", N)
-	for _, r := range results {
-		fmt.Printf("Speedup ≥ %.2f%% → Confidence: %.3f%%\n", r.RelativeSpeedupSampleAvsSampleB*100.0, r.Confidence*100.0)
+	// sink is never meant to be read for its value, only written to during the
+	// batches above so the compiler cannot prove the loop's result is unused
+	// and optimize it away. That makes it invisible to static analysis in a
+	// package main, which can see the whole program and correctly notice that
+	// nothing downstream ever looks at it — unlike the same pattern in this
+	// module's own test files, where a library package's unused-variable check
+	// is deliberately more conservative. This line is the fix: a read that is
+	// itself immediately discarded, costing nothing at runtime, but enough to
+	// tell the linter what the accumulation already told the compiler.
+	_ = sink
+}
+
+// medianCandidate wraps one median implementation as a candidate.
+func medianCandidate(name string, median func([]float64) float64) rtcompare.Candidate {
+	return rtcompare.Candidate{
+		Name: name,
+		Batch: func(n uint64) {
+			rng := rtcompare.NewDPRNG(0x5EED)
+			work := make([]float64, arraySize)
+			var acc float64
+			for range n {
+				// Constant cost in the array's length, so it adds the same
+				// amount to both candidates.
+				for i := range work {
+					work[i] = rng.Float64()
+				}
+				acc += median(work)
+			}
+			sink += acc
+		},
 	}
 }
 
-// fillArray fills the given array with random float64 values using the provided DPRNG.
-// The function modifies the contents of the array in place.
-// This function has constant runtime for an array of a fixed size as rtcompare.DPRNG generates values in constant time.
-func fillArray(array []float64, rng rtcompare.DPRNG) {
-	for i := range array {
-		array[i] = rng.Float64()
+func yesNo(cond bool, yes, no string) string {
+	if cond {
+		return yes
 	}
+	return no
+}
+
+func fail(format string, args ...any) {
+	fmt.Fprintf(os.Stderr, "rtcompare-example: "+format+"\n", args...)
+	os.Exit(1)
 }