Janady13/missionmatch / internal/platform/learn/logistic.go

1// Package learn is the pure Go numerical core of the MAG learning loop: it fits interpretable win 2// probability coefficients from labeled outcomes and calibrates them, with NO database and NO Python. 3// float64 lives ONLY inside these offline solvers (an offline numerical optimization, never a request path 4// score); every coefficient the caller stores and every number the platform surfaces is re quantized to 5// shopspring/decimal by the caller (cmd/train-model), so the served path stays decimal exact. 6package learn 7 8import ( 9 "errors" 10 "math" 11) 12 13// ErrNoRows is returned when a fit is asked over an empty dataset. 14var ErrNoRows = errors.New("learn: no rows") 15 16// LogisticModel is a fitted ridge penalized logistic regression on STANDARDIZED features. Serving computes 17// sigmoid(Intercept + sum_j Coeffs[j] * (x_j - Means[j]) / Scales[j]). Storing the means and scales lets the 18// serve path standardize a new feature vector exactly as training did. 19type LogisticModel struct { 20 Intercept float64 21 Coeffs []float64 // one per feature, on the standardized scale 22 Means []float64 // feature means used to standardize 23 Scales []float64 // feature standard deviations used to standardize 24 Iterations int // Newton iterations to convergence 25} 26 27// FitLogistic fits a ridge penalized logistic regression by iteratively reweighted least squares (Newton 28// Raphson on the penalized log likelihood). x is n rows by p features, y is n labels in {0,1}, l2 is the 29// ridge strength (a zero mean Gaussian prior on the non intercept coefficients, so a thin dataset shrinks 30// toward zero rather than diverging; the intercept is never penalized). Features are standardized internally. 31// Deterministic: the same input always yields the same coefficients. 32func FitLogistic(x [][]float64, y []int, l2 float64) (LogisticModel, error) { 33 n := len(x) 34 if n == 0 { 35 return LogisticModel{}, ErrNoRows 36 } 37 if len(y) != n { 38 return LogisticModel{}, errors.New("learn: x and y length mismatch") 39 } 40 p := len(x[0]) 41 if p == 0 { 42 return LogisticModel{}, errors.New("learn: no features") 43 } 44 45 means := make([]float64, p) 46 scales := make([]float64, p) 47 for j := 0; j < p; j++ { 48 var s float64 49 for i := 0; i < n; i++ { 50 s += x[i][j] 51 } 52 means[j] = s / float64(n) 53 } 54 for j := 0; j < p; j++ { 55 var v float64 56 for i := 0; i < n; i++ { 57 d := x[i][j] - means[j] 58 v += d * d 59 } 60 sd := math.Sqrt(v / float64(n)) 61 if sd < 1e-12 { 62 sd = 1 // a constant feature gets scale 1, so its standardized value is exactly 0 and it contributes nothing 63 } 64 scales[j] = sd 65 } 66 z := make([][]float64, n) 67 for i := 0; i < n; i++ { 68 z[i] = make([]float64, p) 69 for j := 0; j < p; j++ { 70 z[i][j] = (x[i][j] - means[j]) / scales[j] 71 } 72 } 73 74 // A tiny ridge floor keeps the Hessian positive definite even when a feature is constant (a zero 75 // standardized column) or the data is perfectly separable (coefficients would otherwise diverge). It is 76 // small enough not to bias a normal fit and is standard numerical hygiene. 77 effL2 := l2 78 if effL2 < 1e-8 { 79 effL2 = 1e-8 80 } 81 82 // beta index 0 is the intercept, index j+1 is the coefficient of standardized feature j. 83 dim := p + 1 84 beta := make([]float64, dim) 85 const maxIter = 100 86 const tol = 1e-9 87 iters := 0 88 for it := 0; it < maxIter; it++ { 89 iters = it + 1 90 g := make([]float64, dim) // gradient 91 h := make([][]float64, dim) // Hessian 92 for a := range h { 93 h[a] = make([]float64, dim) 94 } 95 for i := 0; i < n; i++ { 96 eta := beta[0] 97 for j := 0; j < p; j++ { 98 eta += beta[j+1] * z[i][j] 99 } 100 pi := sigmoid(eta) 101 w := pi * (1 - pi) 102 if w < 1e-10 { 103 w = 1e-10 // keep the weight matrix positive definite when a point is near certain 104 } 105 r := float64(y[i]) - pi 106 // the design row is xi = [1, z_i0, z_i1, ...]; g += xi*r and h += w * xi xi^T 107 g[0] += r 108 h[0][0] += w 109 for a := 0; a < p; a++ { 110 g[a+1] += z[i][a] * r 111 h[0][a+1] += w * z[i][a] 112 h[a+1][0] += w * z[i][a] 113 for b := 0; b < p; b++ { 114 h[a+1][b+1] += w * z[i][a] * z[i][b] 115 } 116 } 117 } 118 // ridge on the non intercept coefficients only 119 for a := 1; a < dim; a++ { 120 g[a] -= effL2 * beta[a] 121 h[a][a] += effL2 122 } 123 delta, err := solveLinear(h, g) 124 if err != nil { 125 return LogisticModel{}, err 126 } 127 var maxStep float64 128 for a := 0; a < dim; a++ { 129 beta[a] += delta[a] 130 if s := math.Abs(delta[a]); s > maxStep { 131 maxStep = s 132 } 133 } 134 if maxStep < tol { 135 break 136 } 137 } 138 return LogisticModel{Intercept: beta[0], Coeffs: beta[1:], Means: means, Scales: scales, Iterations: iters}, nil 139} 140 141// Predict returns the model probability in [0,1] for a RAW (unstandardized) feature vector. 142func (m LogisticModel) Predict(features []float64) float64 { 143 eta := m.Intercept 144 for j := range m.Coeffs { 145 x := 0.0 146 if j < len(features) && j < len(m.Means) { 147 x = (features[j] - m.Means[j]) / m.Scales[j] 148 } 149 eta += m.Coeffs[j] * x 150 } 151 return sigmoid(eta) 152} 153 154// sigmoid is the numerically stable logistic function. 155func sigmoid(x float64) float64 { 156 if x >= 0 { 157 return 1 / (1 + math.Exp(-x)) 158 } 159 e := math.Exp(x) 160 return e / (1 + e) 161} 162 163// solveLinear solves h x = g for a small dense system by Gaussian elimination with partial pivoting. h is 164// copied, not mutated. It errors on a singular system. 165func solveLinear(h [][]float64, g []float64) ([]float64, error) { 166 n := len(g) 167 m := make([][]float64, n) 168 for i := 0; i < n; i++ { 169 m[i] = make([]float64, n+1) 170 copy(m[i], h[i]) 171 m[i][n] = g[i] 172 } 173 for col := 0; col < n; col++ { 174 piv := col 175 for r := col + 1; r < n; r++ { 176 if math.Abs(m[r][col]) > math.Abs(m[piv][col]) { 177 piv = r 178 } 179 } 180 if math.Abs(m[piv][col]) < 1e-14 { 181 return nil, errors.New("learn: singular system") 182 } 183 m[col], m[piv] = m[piv], m[col] 184 for r := 0; r < n; r++ { 185 if r == col { 186 continue 187 } 188 f := m[r][col] / m[col][col] 189 for c := col; c <= n; c++ { 190 m[r][c] -= f * m[col][c] 191 } 192 } 193 } 194 x := make([]float64, n) 195 for i := 0; i < n; i++ { 196 x[i] = m[i][n] / m[i][i] 197 } 198 return x, nil 199} 200 201// FitPlatt fits Platt scaling: a one dimensional logistic map from raw model scores to calibrated 202// probabilities, returning a and b such that calibrated = sigmoid(a*score + b). It reuses FitLogistic with 203// the score as the single feature and recovers a and b on the raw score scale. Platt is the cold start 204// calibrator: parametric, robust when the label set is thin, later replaced by isotonic past ~1000 labels. 205func FitPlatt(scores []float64, y []int) (a, b float64, err error) { 206 x := make([][]float64, len(scores)) 207 for i := range scores { 208 x[i] = []float64{scores[i]} 209 } 210 m, e := FitLogistic(x, y, 0) 211 if e != nil { 212 return 0, 0, e 213 } 214 // standardized eta = Intercept + Coeff*(s-mean)/scale = (Coeff/scale)*s + (Intercept - Coeff*mean/scale) 215 a = m.Coeffs[0] / m.Scales[0] 216 b = m.Intercept - m.Coeffs[0]*m.Means[0]/m.Scales[0] 217 return a, b, nil 218} 219