Janady13/missionmatch / internal/modules/staffing/optimizer.go

1package staffing 2 3import ( 4 "container/heap" 5 "context" 6 "errors" 7 "sort" 8) 9 10// The portfolio view allocates the bench across all open requirements at once to maximize total 11// assignment value. This is a MAXIMUM PROFIT UNBALANCED ASSIGNMENT, solved exactly as a minimum cost flow 12// by successive shortest paths with node potentials. The math is the classic assignment-as-min-cost-flow 13// (Tomizawa 1971; Edmonds and Karp 1972; Ahuja, Magnanti, Orlin, Network Flows, ch 9 and 12): 14// 15// - Build source -> candidate (cap 1, cost 0), candidate -> requirement (cap 1, cost = negative scaled 16// value, only for qualified pairs), requirement -> sink (cap = open slots, cost 0). 17// - Seed node potentials with one Bellman Ford pass, because maximizing makes the cost edges negative. 18// - Each iteration, find the minimum reduced cost augmenting path with Dijkstra over non negative 19// reduced costs, push one unit, and update potentials by the Dijkstra distances. Potentials keep all 20// residual reduced costs non negative, so Dijkstra stays valid for the life of the solve. 21// - STOP as soon as the cheapest augmenting path has non negative real cost: any further augmentation 22// could not increase total profit. The cost vs flow curve is convex, so this is the exact optimum. 23// 24// It is integer exact (unit capacities, integral supplies, total unimodularity, all int64 arithmetic, no 25// float in the objective) and deterministic (stable node ids from the input order, a fixed adjacency 26// order, and a priority queue that breaks ties by node id). Candidates and slots may be left unmatched, 27// which is exactly a SOURCE recommendation in the domain. The optimizer never qualifies or disqualifies a 28// candidate: the caller adds an edge only where the deterministic engine already qualified the pair. 29 30var ( 31 // ErrProblemTooLarge is returned when the candidate by requirement graph exceeds the configured size 32 // cap. An unbounded solve is a denial of service on our own API, so a too large problem is refused with 33 // a clear error rather than run. 34 ErrProblemTooLarge = errors.New("staffing: portfolio problem exceeds the configured size cap") 35 // ErrSolveBudget is returned when the solve exceeds its wall clock budget (the context deadline). The 36 // solve never hangs. 37 ErrSolveBudget = errors.New("staffing: portfolio solve exceeded the time budget") 38) 39 40// AssignEdge is a qualified candidate to requirement pairing with its scaled objective value (> 0). 41type AssignEdge struct { 42 CandidateIdx int 43 RequirementIdx int 44 Value int64 45} 46 47// PortfolioInput is the assignment problem: how many candidates, how many open slots per requirement, and 48// the qualified, valued edges. Candidate and requirement indices are 0 based and must be stable. 49type PortfolioInput struct { 50 NumCandidates int 51 Slots []int // open slots per requirement; len is the requirement count 52 Edges []AssignEdge 53} 54 55// Assignment is one candidate placed on one requirement in the optimal allocation. 56type Assignment struct { 57 CandidateIdx int 58 RequirementIdx int 59 Value int64 60} 61 62// PortfolioSolution is the optimal max profit allocation. Unassigned candidates and unfilled slots are 63// implicit (anything not in Assignments). 64type PortfolioSolution struct { 65 Assignments []Assignment 66 TotalValue int64 67} 68 69// SolvePortfolio computes the maximum profit allocation. maxEdges and maxNodes bound the problem; the 70// context deadline bounds the wall clock. An empty problem (no candidates, requirements, or edges) returns 71// an empty allocation with no error, not a failure. 72func SolvePortfolio(ctx context.Context, in PortfolioInput, maxEdges, maxNodes int) (PortfolioSolution, error) { 73 R := len(in.Slots) 74 C := in.NumCandidates 75 n := C + R + 2 76 if len(in.Edges) > maxEdges || n > maxNodes { 77 return PortfolioSolution{}, ErrProblemTooLarge 78 } 79 if C == 0 || R == 0 || len(in.Edges) == 0 { 80 return PortfolioSolution{}, nil 81 } 82 83 const source = 0 84 candNode := func(i int) int { return 1 + i } 85 reqNode := func(j int) int { return 1 + C + j } 86 sink := n - 1 87 88 g := newMCMF(n) 89 for i := 0; i < C; i++ { 90 g.addEdge(source, candNode(i), 1, 0) 91 } 92 // Add candidate to requirement edges in a fixed order so the residual graph is built identically every 93 // run (Go map iteration is randomized; the caller's slice order is the source of truth). 94 edges := append([]AssignEdge(nil), in.Edges...) 95 sort.Slice(edges, func(a, b int) bool { 96 if edges[a].CandidateIdx != edges[b].CandidateIdx { 97 return edges[a].CandidateIdx < edges[b].CandidateIdx 98 } 99 return edges[a].RequirementIdx < edges[b].RequirementIdx 100 }) 101 type ref struct { 102 node, idx, candIdx, reqIdx int 103 value int64 104 } 105 var refs []ref 106 for _, e := range edges { 107 if e.Value <= 0 || e.CandidateIdx < 0 || e.CandidateIdx >= C || e.RequirementIdx < 0 || e.RequirementIdx >= R { 108 continue // only profitable, in range, qualified edges 109 } 110 from := candNode(e.CandidateIdx) 111 idx := len(g.graph[from]) 112 g.addEdge(from, reqNode(e.RequirementIdx), 1, -e.Value) 113 refs = append(refs, ref{node: from, idx: idx, candIdx: e.CandidateIdx, reqIdx: e.RequirementIdx, value: e.Value}) 114 } 115 for j := 0; j < R; j++ { 116 if in.Slots[j] > 0 { 117 g.addEdge(reqNode(j), sink, in.Slots[j], 0) 118 } 119 } 120 121 if err := g.solveMaxProfit(ctx, source, sink); err != nil { 122 return PortfolioSolution{}, err 123 } 124 125 var sol PortfolioSolution 126 for _, r := range refs { 127 if g.graph[r.node][r.idx].cap == 0 { // the unit was pushed: the pair is assigned 128 sol.Assignments = append(sol.Assignments, Assignment{CandidateIdx: r.candIdx, RequirementIdx: r.reqIdx, Value: r.value}) 129 sol.TotalValue += r.value 130 } 131 } 132 sort.Slice(sol.Assignments, func(a, b int) bool { 133 if sol.Assignments[a].CandidateIdx != sol.Assignments[b].CandidateIdx { 134 return sol.Assignments[a].CandidateIdx < sol.Assignments[b].CandidateIdx 135 } 136 return sol.Assignments[a].RequirementIdx < sol.Assignments[b].RequirementIdx 137 }) 138 return sol, nil 139} 140 141const infCost int64 = 1 << 62 142 143type mcmfEdge struct { 144 to, rev int 145 cap int 146 cost int64 147} 148 149type mcmf struct { 150 graph [][]mcmfEdge 151} 152 153func newMCMF(n int) *mcmf { return &mcmf{graph: make([][]mcmfEdge, n)} } 154 155func (g *mcmf) addEdge(from, to, capacity int, cost int64) { 156 g.graph[from] = append(g.graph[from], mcmfEdge{to: to, rev: len(g.graph[to]), cap: capacity, cost: cost}) 157 g.graph[to] = append(g.graph[to], mcmfEdge{to: from, rev: len(g.graph[from]) - 1, cap: 0, cost: -cost}) 158} 159 160// solveMaxProfit augments unit flow along minimum cost paths while each adds profit, stopping when the 161// cheapest augmenting path no longer reduces cost. It checks the context each iteration so it never hangs. 162func (g *mcmf) solveMaxProfit(ctx context.Context, s, t int) error { 163 n := len(g.graph) 164 pot := make([]int64, n) 165 for i := range pot { 166 pot[i] = infCost 167 } 168 pot[s] = 0 169 // Bellman Ford initial potentials (cost edges are negative under maximization). O(V*E), once. 170 for iter := 0; iter < n-1; iter++ { 171 changed := false 172 for u := 0; u < n; u++ { 173 if pot[u] == infCost { 174 continue 175 } 176 for _, e := range g.graph[u] { 177 if e.cap > 0 && pot[u]+e.cost < pot[e.to] { 178 pot[e.to] = pot[u] + e.cost 179 changed = true 180 } 181 } 182 } 183 if !changed { 184 break 185 } 186 } 187 // Nodes unreachable from the source have no incoming residual edge from the reachable set, so they are 188 // never traversed; pinning their potential to 0 keeps the arithmetic clean without affecting the solve. 189 for i := range pot { 190 if pot[i] == infCost { 191 pot[i] = 0 192 } 193 } 194 195 dist := make([]int64, n) 196 parentNode := make([]int, n) 197 parentEdge := make([]int, n) 198 for { 199 select { 200 case <-ctx.Done(): 201 return ErrSolveBudget 202 default: 203 } 204 for i := range dist { 205 dist[i] = infCost 206 } 207 dist[s] = 0 208 pq := &nodePQ{{node: s, dist: 0}} 209 for pq.Len() > 0 { 210 it := heap.Pop(pq).(pqItem) 211 if it.dist > dist[it.node] { 212 continue 213 } 214 u := it.node 215 for ei := range g.graph[u] { 216 e := g.graph[u][ei] 217 if e.cap <= 0 { 218 continue 219 } 220 rc := e.cost + pot[u] - pot[e.to] // reduced cost, non negative by the potential invariant 221 nd := dist[u] + rc 222 if nd < dist[e.to] { 223 dist[e.to] = nd 224 parentNode[e.to] = u 225 parentEdge[e.to] = ei 226 heap.Push(pq, pqItem{node: e.to, dist: nd}) 227 } 228 } 229 } 230 if dist[t] == infCost { 231 break // no augmenting path remains 232 } 233 for v := 0; v < n; v++ { 234 if dist[v] < infCost { 235 pot[v] += dist[v] 236 } 237 } 238 // pot[s] stays 0, so pot[t] is the real cost of the cheapest augmenting path (the marginal cost of 239 // the next unit). Stop as soon as it is non negative: no further augmentation can add profit. 240 if pot[t]-pot[s] >= 0 { 241 break 242 } 243 // Augment one unit along the path (the bottleneck is 1, bounded by the source to candidate edge). 244 push := 1 << 30 245 for v := t; v != s; v = parentNode[v] { 246 if c := g.graph[parentNode[v]][parentEdge[v]].cap; c < push { 247 push = c 248 } 249 } 250 for v := t; v != s; v = parentNode[v] { 251 pn, pe := parentNode[v], parentEdge[v] 252 g.graph[pn][pe].cap -= push 253 g.graph[v][g.graph[pn][pe].rev].cap += push 254 } 255 } 256 return nil 257} 258 259// AssignCounterfactual is the value matrix explanation for one assignment: where else this candidate could 260// best go and who else could best serve this requirement, holding the rest of the solution fixed. Indices 261// are -1 when there is no alternative. This is the interpretable explanation, NOT a true linear program re 262// solve marginal; the two can differ, and the interpretable next best is the defensible choice for an AEDT 263// because a reviewer can verify it directly against the value matrix. 264type AssignCounterfactual struct { 265 CandidateNextReqIdx int 266 CandidateValueDelta int64 267 RequirementNextCandIdx int 268 RequirementValueDelta int64 269} 270 271// Counterfactuals computes, for each assignment, the candidate's next best requirement and the 272// requirement's next best candidate by inspecting the value matrix. Pure and deterministic (ties break to 273// the smaller index). 274func Counterfactuals(in PortfolioInput, sol PortfolioSolution) map[int]AssignCounterfactual { 275 byCand := map[int][]AssignEdge{} 276 byReq := map[int][]AssignEdge{} 277 for _, e := range in.Edges { 278 if e.Value <= 0 { 279 continue 280 } 281 byCand[e.CandidateIdx] = append(byCand[e.CandidateIdx], e) 282 byReq[e.RequirementIdx] = append(byReq[e.RequirementIdx], e) 283 } 284 out := make(map[int]AssignCounterfactual, len(sol.Assignments)) 285 for _, a := range sol.Assignments { 286 cf := AssignCounterfactual{CandidateNextReqIdx: -1, RequirementNextCandIdx: -1} 287 bestV := int64(-1) 288 for _, e := range byCand[a.CandidateIdx] { 289 if e.RequirementIdx == a.RequirementIdx { 290 continue 291 } 292 if e.Value > bestV || (e.Value == bestV && e.RequirementIdx < cf.CandidateNextReqIdx) { 293 bestV = e.Value 294 cf.CandidateNextReqIdx = e.RequirementIdx 295 } 296 } 297 if cf.CandidateNextReqIdx != -1 { 298 cf.CandidateValueDelta = a.Value - bestV 299 } 300 bestV2 := int64(-1) 301 for _, e := range byReq[a.RequirementIdx] { 302 if e.CandidateIdx == a.CandidateIdx { 303 continue 304 } 305 if e.Value > bestV2 || (e.Value == bestV2 && e.CandidateIdx < cf.RequirementNextCandIdx) { 306 bestV2 = e.Value 307 cf.RequirementNextCandIdx = e.CandidateIdx 308 } 309 } 310 if cf.RequirementNextCandIdx != -1 { 311 cf.RequirementValueDelta = a.Value - bestV2 312 } 313 out[a.CandidateIdx] = cf 314 } 315 return out 316} 317 318type pqItem struct { 319 node int 320 dist int64 321} 322 323// nodePQ is a min heap keyed by (dist, node id). The node id tiebreak makes the settle order, and thus the 324// chosen augmenting path and the final assignment, fully deterministic. 325type nodePQ []pqItem 326 327func (p nodePQ) Len() int { return len(p) } 328func (p nodePQ) Less(i, j int) bool { 329 if p[i].dist != p[j].dist { 330 return p[i].dist < p[j].dist 331 } 332 return p[i].node < p[j].node 333} 334func (p nodePQ) Swap(i, j int) { p[i], p[j] = p[j], p[i] } 335func (p *nodePQ) Push(x any) { *p = append(*p, x.(pqItem)) } 336func (p *nodePQ) Pop() any { 337 old := *p 338 n := len(old) 339 it := old[n-1] 340 *p = old[:n-1] 341 return it 342} 343