feat(loadtest): add portal-loadtest uniformity probe (chi-square vs uniform N/K)
cognitive committed
Apr 30, 2026 at 05:36 UTC
a3e3646e5e6087a8b2c83b8a67d77d29c0bbd159
1 file changed
+237
cmd/portal-loadtest/main.go
new
+237
@@ -0,0 +1,237 @@
1
+// Command portal-loadtest is a Phase 1 uniformity probe that measures
2
+// how evenly the MOLS relay-selection policy distributes N synthetic clients
3
+// across K synthetic relays. It runs entirely in-process — no running
4
+// portal-tunnel server is required.
5
+//
6
+// Flags (Phase 1 only — -capacities and -selector are Phase 2):
7
+//
8
+// -clients N number of synthetic clients (default 100)
9
+// -relays K number of synthetic relays (default 5)
10
+// -multi-hop D multi-hop depth (0 = priority/single-hop; ≥2 = multi-hop)
11
+//
12
+// Output: per-relay top-pick histogram, chi-square statistic against the
13
+// uniform expected distribution N/K, and a p-value.
14
+//
15
+// P-value method: regularized upper incomplete gamma function Q(k/2, x/2),
16
+// implemented via the series expansion (|x| < s+1) and continued-fraction
17
+// expansion (x ≥ s+1) from Numerical Recipes §6.2. This gives accurate
18
+// results even at small df values (e.g. df=4 for K=5).
19
+package main
20
+
21
+import (
22
+ "flag"
23
+ "fmt"
24
+ "math"
25
+ "os"
26
+ "sort"
27
+ "time"
28
+
29
+ "github.com/gosuda/portal-tunnel/v2/portal/discovery"
30
+ "github.com/gosuda/portal-tunnel/v2/types"
31
+)
32
+
33
+func main() {
34
+ clients := flag.Int("clients", 100, "number of synthetic clients")
35
+ relays := flag.Int("relays", 5, "number of synthetic relays")
36
+ multiHop := flag.Int("multi-hop", 0, "multi-hop depth (0 = priority; ≥2 = multi-hop)")
37
+ flag.Parse()
38
+
39
+ if *clients <= 0 {
40
+ fmt.Fprintln(os.Stderr, "portal-loadtest: -clients must be > 0")
41
+ os.Exit(1)
42
+ }
43
+ if *relays <= 0 {
44
+ fmt.Fprintln(os.Stderr, "portal-loadtest: -relays must be > 0")
45
+ os.Exit(1)
46
+ }
47
+ // MultiHopDepth ≤ 1 causes SelectMultiHop to return nil (see mols.go).
48
+ // Reject 1 explicitly; 0 means priority mode.
49
+ if *multiHop == 1 {
50
+ fmt.Fprintln(os.Stderr, "portal-loadtest: -multi-hop=1 is not valid; use 0 for priority or ≥2 for multi-hop")
51
+ os.Exit(1)
52
+ }
53
+
54
+ mode := "priority"
55
+ if *multiHop >= 2 {
56
+ mode = "multihop"
57
+ }
58
+
59
+ // Build K synthetic relay states. We construct discovery.RelayState values
60
+ // directly (not via RelaySet.InsertAnnounced) because the public announce
61
+ // path requires real EVM-signed descriptors. MOLSRelayPolicy is called
62
+ // directly so that no signature gate runs.
63
+ //
64
+ // For priority mode: states without an observed descriptor (LastSeenAt zero)
65
+ // are accepted into the auto pool by SelectPriorityWithTrace — the
66
+ // expiry/protocol gates only fire when hasObservedDescriptor() is true.
67
+ //
68
+ // For multi-hop mode: SelectMultiHopWithTrace requires hasObservedDescriptor,
69
+ // a non-expired ExpiresAt, and HasOverlayPeer()==true. We populate those
70
+ // fields with dummy-but-valid values using a far-future ExpiresAt and a
71
+ // syntactically valid WireGuard public key placeholder.
72
+ now := time.Now().UTC()
73
+ relayStates := make([]discovery.RelayState, *relays)
74
+ for i := range relayStates {
75
+ relayURL := fmt.Sprintf("https://test-relay-%d.example", i+1)
76
+ rs := discovery.RelayState{
77
+ Descriptor: types.RelayDescriptor{
78
+ APIHTTPSAddr: relayURL,
79
+ },
80
+ }
81
+ if mode == "multihop" {
82
+ // Populate the fields required by SelectMultiHopWithTrace's eligibility
83
+ // gates: hasObservedDescriptor (LastSeenAt non-zero), valid ExpiresAt,
84
+ // and HasOverlayPeer() = SupportsOverlay && WireGuardPublicKey != "" &&
85
+ // WireGuardPort in [1, 65535].
86
+ rs.LastSeenAt = now
87
+ rs.Descriptor.IssuedAt = now
88
+ rs.Descriptor.ExpiresAt = now.Add(24 * time.Hour)
89
+ rs.Descriptor.SupportsOverlay = true
90
+ rs.Descriptor.WireGuardPublicKey = fmt.Sprintf("synthetic-wg-key-%d", i+1)
91
+ rs.Descriptor.WireGuardPort = 51820
92
+ }
93
+ relayStates[i] = rs
94
+ }
95
+
96
+ // Generate N synthetic client states with UNIQUE LocalAddress values.
97
+ // MOLS is deterministic on (LocalAddress, relayURL): duplicate addresses
98
+ // would make all clients pick identically, falsely appearing as 100% imbalance.
99
+ policy := discovery.MOLSRelayPolicy{}
100
+ picks := make(map[string]int, *relays) // relay URL → count of clients that picked it first
101
+ for i := 0; i < *clients; i++ {
102
+ cs := discovery.ClientState{
103
+ LocalAddress: fmt.Sprintf("synthetic-client-%d", i),
104
+ MultiHopDepth: *multiHop,
105
+ }
106
+ var outputURLs []string
107
+ if mode == "multihop" {
108
+ outputURLs, _ = policy.SelectMultiHopWithTrace(relayStates, cs)
109
+ } else {
110
+ outputURLs, _ = policy.SelectPriorityWithTrace(relayStates, cs)
111
+ }
112
+ if len(outputURLs) == 0 {
113
+ // All relays were filtered; skip this client.
114
+ continue
115
+ }
116
+ picks[outputURLs[0]]++
117
+ }
118
+
119
+ // Collect and sort relay URLs for deterministic output.
120
+ relayURLs := make([]string, 0, *relays)
121
+ for i := range relayStates {
122
+ relayURLs = append(relayURLs, relayStates[i].Descriptor.APIHTTPSAddr)
123
+ }
124
+ sort.Strings(relayURLs)
125
+
126
+ expected := float64(*clients) / float64(*relays)
127
+
128
+ // Chi-square statistic: Σ (observed - expected)^2 / expected
129
+ var chi2 float64
130
+ for _, url := range relayURLs {
131
+ obs := float64(picks[url])
132
+ diff := obs - expected
133
+ chi2 += diff * diff / expected
134
+ }
135
+
136
+ df := *relays - 1
137
+
138
+ // P-value: P(χ² > chi2 | df) = Q(df/2, chi2/2) = igamc(df/2, chi2/2)
139
+ // using the regularized upper incomplete gamma function.
140
+ pval := igamc(float64(df)/2.0, chi2/2.0)
141
+
142
+ // Print results.
143
+ header := fmt.Sprintf("portal-loadtest: N=%d clients, K=%d relays, mode=%s", *clients, *relays, mode)
144
+ fmt.Println(header)
145
+ fmt.Printf("%-45s %6s %8s\n", "relay", "picks", "expected")
146
+ fmt.Println("---------------------------------------------------------------")
147
+ for _, url := range relayURLs {
148
+ fmt.Printf("%-45s %6d %8.1f\n", url, picks[url], expected)
149
+ }
150
+ fmt.Printf("\nchi-square: %.4f\n", chi2)
151
+ fmt.Printf("df: %d\n", df)
152
+ fmt.Printf("p-value: %.4f\n", pval)
153
+}
154
+
155
+// igamc returns the regularized upper incomplete gamma function Q(s, x),
156
+// also written Γ(s, x) / Γ(s). This equals 1 - P(s, x) where P(s, x) is
157
+// the regularized lower incomplete gamma.
158
+//
159
+// For s < x+1 the continued-fraction expansion converges faster; otherwise
160
+// the series expansion is used. Algorithm from Numerical Recipes §6.2
161
+// (Press et al.). Accurate to ~1e-7 for the parameter ranges used here
162
+// (s = df/2 ≥ 0.5, x = chi2/2 ≥ 0).
163
+func igamc(s, x float64) float64 {
164
+ if x < 0 || s <= 0 {
165
+ return 1.0
166
+ }
167
+ if x == 0 {
168
+ return 1.0
169
+ }
170
+
171
+ if x < s+1 {
172
+ // Series expansion for the lower incomplete gamma P(s, x);
173
+ // return Q = 1 - P.
174
+ return 1.0 - gamSer(s, x)
175
+ }
176
+ // Continued-fraction expansion for Q(s, x) directly.
177
+ return gamCF(s, x)
178
+}
179
+
180
+// gamSer computes P(s, x) via a series expansion. P(s, x) = e^(-x) * x^s *
181
+// Σ_{n=0}^∞ x^n / Γ(s+n+1).
182
+func gamSer(s, x float64) float64 {
183
+ const maxIter = 200
184
+ const eps = 3e-7
185
+
186
+ ap := s
187
+ del := 1.0 / s
188
+ sum := del
189
+ for n := 0; n < maxIter; n++ {
190
+ ap++
191
+ del *= x / ap
192
+ sum += del
193
+ if math.Abs(del) < math.Abs(sum)*eps {
194
+ return sum * math.Exp(-x+s*math.Log(x)-lgamma(s))
195
+ }
196
+ }
197
+ // Did not converge; return best estimate.
198
+ return sum * math.Exp(-x+s*math.Log(x)-lgamma(s))
199
+}
200
+
201
+// gamCF computes Q(s, x) via a modified Lentz continued-fraction expansion.
202
+func gamCF(s, x float64) float64 {
203
+ const maxIter = 200
204
+ const eps = 3e-7
205
+ const fpMin = 1e-300
206
+
207
+ b := x + 1.0 - s
208
+ c := 1.0 / fpMin
209
+ d := 1.0 / b
210
+ h := d
211
+ for i := 1; i <= maxIter; i++ {
212
+ an := -float64(i) * (float64(i) - s)
213
+ b += 2.0
214
+ d = an*d + b
215
+ if math.Abs(d) < fpMin {
216
+ d = fpMin
217
+ }
218
+ c = b + an/c
219
+ if math.Abs(c) < fpMin {
220
+ c = fpMin
221
+ }
222
+ d = 1.0 / d
223
+ del := d * c
224
+ h *= del
225
+ if math.Abs(del-1.0) < eps {
226
+ break
227
+ }
228
+ }
229
+ return math.Exp(-x+s*math.Log(x)-lgamma(s)) * h
230
+}
231
+
232
+// lgamma returns the natural log of the Gamma function using the standard
233
+// library, which is accurate for all positive real inputs.
234
+func lgamma(x float64) float64 {
235
+ lg, _ := math.Lgamma(x)
236
+ return lg
237
+}