| 1 | // Command portal-loadtest is a Phase 1 uniformity probe that measures |
| 2 | // how evenly MOLS relay selection 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. Selection functions are 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 SelectPriority; the expiry/protocol |
| 66 | // gates only fire when hasObservedDescriptor() is true. |
| 67 | // |
| 68 | // For multi-hop mode: SelectMultiHop requires hasObservedDescriptor, a |
| 69 | // 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 SelectMultiHop'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 | picks := make(map[string]int, *relays) // relay URL → count of clients that picked it first |
| 100 | for i := 0; i < *clients; i++ { |
| 101 | cs := discovery.RouteState{ |
| 102 | LocalAddress: fmt.Sprintf("synthetic-client-%d", i), |
| 103 | MultiHopDepth: *multiHop, |
| 104 | } |
| 105 | var outputURLs []string |
| 106 | if mode == "multihop" { |
| 107 | outputURLs = discovery.SelectMultiHop(relayStates, cs) |
| 108 | } else { |
| 109 | outputURLs = discovery.SelectPriority(relayStates, cs) |
| 110 | } |
| 111 | if len(outputURLs) == 0 { |
| 112 | // All relays were filtered; skip this client. |
| 113 | continue |
| 114 | } |
| 115 | picks[outputURLs[0]]++ |
| 116 | } |
| 117 | |
| 118 | // Collect and sort relay URLs for deterministic output. |
| 119 | relayURLs := make([]string, 0, *relays) |
| 120 | for i := range relayStates { |
| 121 | relayURLs = append(relayURLs, relayStates[i].Descriptor.APIHTTPSAddr) |
| 122 | } |
| 123 | sort.Strings(relayURLs) |
| 124 | |
| 125 | expected := float64(*clients) / float64(*relays) |
| 126 | |
| 127 | // Chi-square statistic: Σ (observed - expected)^2 / expected |
| 128 | var chi2 float64 |
| 129 | for _, url := range relayURLs { |
| 130 | obs := float64(picks[url]) |
| 131 | diff := obs - expected |
| 132 | chi2 += diff * diff / expected |
| 133 | } |
| 134 | |
| 135 | df := *relays - 1 |
| 136 | |
| 137 | // P-value: P(χ² > chi2 | df) = Q(df/2, chi2/2) = igamc(df/2, chi2/2) |
| 138 | // using the regularized upper incomplete gamma function. |
| 139 | pval := igamc(float64(df)/2.0, chi2/2.0) |
| 140 | |
| 141 | // Print results. |
| 142 | header := fmt.Sprintf("portal-loadtest: N=%d clients, K=%d relays, mode=%s", *clients, *relays, mode) |
| 143 | fmt.Println(header) |
| 144 | fmt.Printf("%-45s %6s %8s\n", "relay", "picks", "expected") |
| 145 | fmt.Println("---------------------------------------------------------------") |
| 146 | for _, url := range relayURLs { |
| 147 | fmt.Printf("%-45s %6d %8.1f\n", url, picks[url], expected) |
| 148 | } |
| 149 | fmt.Printf("\nchi-square: %.4f\n", chi2) |
| 150 | fmt.Printf("df: %d\n", df) |
| 151 | fmt.Printf("p-value: %.4f\n", pval) |
| 152 | } |
| 153 | |
| 154 | // igamc returns the regularized upper incomplete gamma function Q(s, x), |
| 155 | // also written Γ(s, x) / Γ(s). This equals 1 - P(s, x) where P(s, x) is |
| 156 | // the regularized lower incomplete gamma. |
| 157 | // |
| 158 | // For s < x+1 the continued-fraction expansion converges faster; otherwise |
| 159 | // the series expansion is used. Algorithm from Numerical Recipes §6.2 |
| 160 | // (Press et al.). Accurate to ~1e-7 for the parameter ranges used here |
| 161 | // (s = df/2 ≥ 0.5, x = chi2/2 ≥ 0). |
| 162 | func igamc(s, x float64) float64 { |
| 163 | if x < 0 || s <= 0 { |
| 164 | return 1.0 |
| 165 | } |
| 166 | if x == 0 { |
| 167 | return 1.0 |
| 168 | } |
| 169 | |
| 170 | if x < s+1 { |
| 171 | // Series expansion for the lower incomplete gamma P(s, x); |
| 172 | // return Q = 1 - P. |
| 173 | return 1.0 - gamSer(s, x) |
| 174 | } |
| 175 | // Continued-fraction expansion for Q(s, x) directly. |
| 176 | return gamCF(s, x) |
| 177 | } |
| 178 | |
| 179 | // gamSer computes P(s, x) via a series expansion. P(s, x) = e^(-x) * x^s * |
| 180 | // Σ_{n=0}^∞ x^n / Γ(s+n+1). |
| 181 | func gamSer(s, x float64) float64 { |
| 182 | const maxIter = 200 |
| 183 | const eps = 3e-7 |
| 184 | |
| 185 | ap := s |
| 186 | del := 1.0 / s |
| 187 | sum := del |
| 188 | for n := 0; n < maxIter; n++ { |
| 189 | ap++ |
| 190 | del *= x / ap |
| 191 | sum += del |
| 192 | if math.Abs(del) < math.Abs(sum)*eps { |
| 193 | return sum * math.Exp(-x+s*math.Log(x)-lgamma(s)) |
| 194 | } |
| 195 | } |
| 196 | // Did not converge; return best estimate. |
| 197 | return sum * math.Exp(-x+s*math.Log(x)-lgamma(s)) |
| 198 | } |
| 199 | |
| 200 | // gamCF computes Q(s, x) via a modified Lentz continued-fraction expansion. |
| 201 | func gamCF(s, x float64) float64 { |
| 202 | const maxIter = 200 |
| 203 | const eps = 3e-7 |
| 204 | const fpMin = 1e-300 |
| 205 | |
| 206 | b := x + 1.0 - s |
| 207 | c := 1.0 / fpMin |
| 208 | d := 1.0 / b |
| 209 | h := d |
| 210 | for i := 1; i <= maxIter; i++ { |
| 211 | an := -float64(i) * (float64(i) - s) |
| 212 | b += 2.0 |
| 213 | d = an*d + b |
| 214 | if math.Abs(d) < fpMin { |
| 215 | d = fpMin |
| 216 | } |
| 217 | c = b + an/c |
| 218 | if math.Abs(c) < fpMin { |
| 219 | c = fpMin |
| 220 | } |
| 221 | d = 1.0 / d |
| 222 | del := d * c |
| 223 | h *= del |
| 224 | if math.Abs(del-1.0) < eps { |
| 225 | break |
| 226 | } |
| 227 | } |
| 228 | return math.Exp(-x+s*math.Log(x)-lgamma(s)) * h |
| 229 | } |
| 230 | |
| 231 | // lgamma returns the natural log of the Gamma function using the standard |
| 232 | // library, which is accurate for all positive real inputs. |
| 233 | func lgamma(x float64) float64 { |
| 234 | lg, _ := math.Lgamma(x) |
| 235 | return lg |
| 236 | } |