1 package profile
2 3 import (
4 "sort"
5 6 "git.mleku.dev/mleku/dendrite/pkg/lattice"
7 "git.mleku.dev/mleku/dendrite/pkg/ratio"
8 )
9 10 // Stats holds derived statistical measures computed from a Profile.
11 // All values use exact rational arithmetic for determinism.
12 type Stats struct {
13 // PathEntropy is the Shannon entropy of the path frequency distribution.
14 // Higher means more diverse bonding across the lattice.
15 PathEntropy ratio.Ratio `json:"path_entropy"`
16 17 // SurprisalVariance is the variance of -log2(p) across bonded paths.
18 // Human text has higher variance (bursty surprisal); AI text is smoother.
19 SurprisalVariance ratio.Ratio `json:"surprisal_variance"`
20 21 // BurstinessGini is the Gini coefficient of bond counts across nodes.
22 // Measures inequality in how bonds distribute across the lattice.
23 // Human text is burstier (higher Gini); AI text distributes more evenly.
24 BurstinessGini ratio.Ratio `json:"burstiness_gini"`
25 26 // VertexCoverage is the fraction of lattice nodes that received at least
27 // one bond during the growth session.
28 VertexCoverage ratio.Ratio `json:"vertex_coverage"`
29 30 // AvgWalkDistance is the mean number of walk steps before bonding.
31 AvgWalkDistance ratio.Ratio `json:"avg_walk_distance"`
32 33 // BondRate is bonds / tokens ingested.
34 BondRate ratio.Ratio `json:"bond_rate"`
35 36 // NewVertexRate is new vertices / tokens ingested.
37 NewVertexRate ratio.Ratio `json:"new_vertex_rate"`
38 39 // TransitionEntropy is the Shannon entropy of the bigram transition
40 // frequency distribution. Measures sequential diversity.
41 TransitionEntropy ratio.Ratio `json:"transition_entropy"`
42 }
43 44 // Compute derives statistics from a raw profile.
45 // latticeSize is the total number of nodes in the lattice at time of profiling.
46 func Compute(p *Profile, latticeSize int) Stats {
47 var s Stats
48 49 if p.TokensIngested == 0 {
50 return s
51 }
52 53 // Bond rate.
54 s.BondRate = ratio.New(p.BondEvents, p.TokensIngested)
55 56 // New vertex rate.
57 s.NewVertexRate = ratio.New(p.NewVertices, p.TokensIngested)
58 59 // Vertex coverage.
60 if latticeSize > 0 {
61 s.VertexCoverage = ratio.New(int64(len(p.PathFreq)), int64(latticeSize))
62 }
63 64 // Average walk distance.
65 if p.BondEvents > 0 {
66 totalSteps := int64(0)
67 for steps, count := range p.WalkDistHist {
68 totalSteps += int64(steps) * count
69 }
70 s.AvgWalkDistance = ratio.New(totalSteps, p.BondEvents)
71 }
72 73 // Path entropy and surprisal variance.
74 s.PathEntropy, s.SurprisalVariance = computeEntropyAndVariance(p.PathFreq, p.BondEvents)
75 76 // Burstiness (Gini coefficient of path frequencies).
77 s.BurstinessGini = computeGini(p.PathFreq)
78 79 // Transition entropy.
80 totalTransitions := int64(0)
81 transitionCounts := make(map[[2]string]int64, len(p.TransitionFreq))
82 for k, v := range p.TransitionFreq {
83 transitionCounts[k] = v
84 totalTransitions += v
85 }
86 if totalTransitions > 0 {
87 // Convert to a flat frequency map for entropy computation.
88 flatFreq := make(map[int]int64, len(transitionCounts))
89 idx := 0
90 for _, v := range transitionCounts {
91 flatFreq[idx] = v
92 idx++
93 }
94 s.TransitionEntropy, _ = computeEntropyAndVarianceGeneric(flatFreq, totalTransitions)
95 }
96 97 return s
98 }
99 100 // computeEntropyAndVariance computes Shannon entropy and surprisal variance
101 // from a frequency distribution. Uses rational arithmetic throughout.
102 //
103 // Entropy = -sum(p_i * log2(p_i)) approximated as:
104 // H = log2(N) - (1/N) * sum(f_i * log2(f_i))
105 // where f_i are frequencies and N = sum(f_i).
106 //
107 // Since exact log2 of rationals produces irrationals, we approximate using
108 // integer log2 (floor). This introduces bounded error but preserves the
109 // comparison ordering: distributions with higher true entropy will have
110 // higher approximated entropy.
111 func computeEntropyAndVariance(freqs map[lattice.NodeID]int64, total int64) (entropy, variance ratio.Ratio) {
112 if total == 0 || len(freqs) == 0 {
113 return ratio.Zero, ratio.Zero
114 }
115 116 // Compute entropy approximation: log2(total) - (1/total) * sum(f * log2(f))
117 totalR := ratio.FromInt(total)
118 log2Total := log2Scaled(total)
119 120 sumFLogF := ratio.Zero
121 for _, f := range freqs {
122 if f > 0 {
123 sumFLogF = sumFLogF.Add(ratio.FromInt(f).Mul(log2Scaled(f)))
124 }
125 }
126 127 entropy = log2Total.Sub(sumFLogF.Div(totalR))
128 if entropy.IsNegative() {
129 entropy = ratio.Zero
130 }
131 132 // Surprisal variance: Var(-log2(p_i)) where p_i = f_i/total.
133 // Mean surprisal is entropy H.
134 // Var = E[S^2] - H^2 where S_i = -log2(p_i) = log2(total) - log2(f_i).
135 sumS2 := ratio.Zero
136 for _, f := range freqs {
137 if f > 0 {
138 s := log2Total.Sub(log2Scaled(f))
139 s2 := s.Mul(s)
140 p := ratio.New(f, total)
141 sumS2 = sumS2.Add(p.Mul(s2))
142 }
143 }
144 variance = sumS2.Sub(entropy.Mul(entropy))
145 if variance.IsNegative() {
146 variance = ratio.Zero
147 }
148 149 return entropy, variance
150 }
151 152 // computeEntropyAndVarianceGeneric is the same as computeEntropyAndVariance
153 // but works with int-keyed maps (for transition frequencies).
154 func computeEntropyAndVarianceGeneric(freqs map[int]int64, total int64) (entropy, variance ratio.Ratio) {
155 if total == 0 || len(freqs) == 0 {
156 return ratio.Zero, ratio.Zero
157 }
158 159 totalR := ratio.FromInt(total)
160 log2Total := log2Scaled(total)
161 162 sumFLogF := ratio.Zero
163 for _, f := range freqs {
164 if f > 0 {
165 sumFLogF = sumFLogF.Add(ratio.FromInt(f).Mul(log2Scaled(f)))
166 }
167 }
168 169 entropy = log2Total.Sub(sumFLogF.Div(totalR))
170 if entropy.IsNegative() {
171 entropy = ratio.Zero
172 }
173 174 sumS2 := ratio.Zero
175 for _, f := range freqs {
176 if f > 0 {
177 s := log2Total.Sub(log2Scaled(f))
178 s2 := s.Mul(s)
179 p := ratio.New(f, total)
180 sumS2 = sumS2.Add(p.Mul(s2))
181 }
182 }
183 variance = sumS2.Sub(entropy.Mul(entropy))
184 if variance.IsNegative() {
185 variance = ratio.Zero
186 }
187 188 return entropy, variance
189 }
190 191 // computeGini computes the Gini coefficient of a frequency distribution.
192 // Gini = (2 * sum_i(i * x_i)) / (n * sum(x_i)) - (n+1)/n
193 // where x_i are the sorted values and i is 1-based rank.
194 func computeGini(freqs map[lattice.NodeID]int64) ratio.Ratio {
195 if len(freqs) == 0 {
196 return ratio.Zero
197 }
198 199 // Extract and sort values.
200 vals := make([]int64, 0, len(freqs))
201 for _, v := range freqs {
202 vals = append(vals, v)
203 }
204 sort.Slice(vals, func(i, j int) bool { return vals[i] < vals[j] })
205 206 n := int64(len(vals))
207 total := int64(0)
208 weightedSum := int64(0)
209 for i, v := range vals {
210 total += v
211 weightedSum += int64(i+1) * v
212 }
213 214 if total == 0 || n == 0 {
215 return ratio.Zero
216 }
217 218 // Gini = (2 * weightedSum) / (n * total) - (n + 1) / n
219 term1 := ratio.New(2*weightedSum, n*total)
220 term2 := ratio.New(n+1, n)
221 g := term1.Sub(term2)
222 if g.IsNegative() {
223 return ratio.Zero
224 }
225 return g
226 }
227 228 // log2Scale is the fixed-point scaling factor for log2 computations.
229 // Using 1024 (2^10) gives 10 bits of fractional precision while
230 // keeping all arithmetic in exact rationals.
231 const log2Scale = 1024
232 233 // log2Scaled returns an approximation of log2(n) as a ratio with
234 // denominator log2Scale. For n > 0, computes floor(log2(n)) as the
235 // integer part, then refines the fractional part using bisection
236 // on the residual: after extracting the integer bits k such that
237 // 2^k <= n < 2^(k+1), the fractional part is log2(n / 2^k) which
238 // lies in [0, 1). We approximate this by testing whether
239 // n^2 >= 2^(2k+1) repeatedly at increasing resolution.
240 //
241 // Returns ratio.Zero for n <= 0.
242 func log2Scaled(n int64) ratio.Ratio {
243 if n <= 0 {
244 return ratio.Zero
245 }
246 if n == 1 {
247 return ratio.Zero
248 }
249 250 // Integer part: floor(log2(n)).
251 k := int64(0)
252 v := n
253 for v > 1 {
254 v >>= 1
255 k++
256 }
257 258 // Fractional part via repeated squaring / bisection.
259 // We compute frac * log2Scale as an integer.
260 // Start with remainder r = n, base = 2^k.
261 // At each step, square r, check if r^2 >= base*2, if so
262 // the next bit of the fractional log is 1.
263 frac := int64(0)
264 // r represents n / 2^k as a fraction in [1, 2).
265 // We track r_num / r_den where initially r = n / 2^k.
266 rNum := n
267 rDen := int64(1) << k
268 if rDen <= 0 {
269 // Overflow protection for very large k.
270 return ratio.New(k*log2Scale, log2Scale)
271 }
272 273 for bit := int64(log2Scale / 2); bit > 0; bit >>= 1 {
274 // Square: r = r^2 / 2 (normalize to keep in [1, 2) range).
275 // Actually: r' = r^2. If r' >= 2, then this fractional bit is 1
276 // and r' = r' / 2.
277 rNum = rNum * rNum
278 rDen = rDen * rDen
279 280 // Check overflow — if numbers get too big, stop refining.
281 if rNum < 0 || rDen < 0 {
282 break
283 }
284 285 // If r^2 >= 2 (i.e., rNum/rDen >= 2, i.e., rNum >= 2*rDen):
286 if rNum >= 2*rDen {
287 frac += bit
288 // r = r / 2: keep rNum, double rDen.
289 rDen *= 2
290 }
291 }
292 293 return ratio.New(k*log2Scale+frac, log2Scale)
294 }
295