diff --git a/README.md b/README.md index fa0a94b..3c3ee2c 100644 --- a/README.md +++ b/README.md @@ -62,18 +62,28 @@ It fails on any failure not listed in `testdata/pulp_known_failures.txt`. target instances it trades the tuned root trajectory for no bound gain. - **Basis factorization**: singleton triangularization with a sparse-LU - kernel and product-form (eta) updates — O(nnz) pivots, periodic - refactorization. No dense inverse. + kernel (in-place row elimination, counted arenas — allocation-free on + the hot path) and product-form (eta) updates — O(nnz) pivots, periodic + refactorization. No dense inverse. Per-factor data is int32-compacted + and arena-consolidated, and all solve scratch is shared per LP, so a + refactorization costs a handful of allocations instead of a dozen. - **Presolve**: iterated activity-based bound tightening, big-M coefficient tightening for binaries, and CglProbing-style binary probing (infeasibility fixing plus integer-only merged implied bounds). - **Cuts**: Gomory mixed-integer cuts at the root, with support/dynamism - hygiene, bound-driven round control, and retraction (with retries) of - batches that degrade the LP numerically; probing implication cuts + hygiene, and retraction (with retries) of batches that degrade the LP + numerically; rounds are budgeted in pivots (speed-invariant work + units), so the cut set is a deterministic function of the problem and + survives engine speed changes — wall clock only remains as a safety + cap; probing implication cuts (CglProbing as a cut generator) on large instances, with slackened implied bounds so propagation drift can never cut off the optimum; TwoMIR-lite cuts (sparse pairwise tableau-row aggregations through the - MIR derivation) on large instances; slack cuts are dropped before the + MIR derivation) on large instances; single-row MIR cuts + (CglMixedIntegerRounding2-style, VUB/bound substitution with a divisor + search) seed the first two rounds on large instances — later MIR + rounds are measured-negative (row bloat, face drift); slack cuts are + dropped before the tree so only root-active rows ride into node re-solves. - **Branch and bound**: best-first with depth-first plunging, warm-started child bases, node-level bound propagation on branching, monotone bound @@ -87,12 +97,14 @@ It fails on any failure not listed in `testdata/pulp_known_failures.txt`. the incumbent fixes the variable at the node without spending a branch; pseudocost selection deeper (reliability-branching shape). - **Heuristics**: caller-provided MIP start (`mip.Model.MIPStart`, - completed before the cut loop so reduced-cost fixing bites), 1-opt - incumbent polish (CbcHeuristicLocal-style binary flips via warm dual - re-solves), face walk (least-degradation dive along the LP-optimal + completed via a warm child solve before the cut loop so reduced-cost + fixing bites; the trivial start is deliberately not polished — measured + as pure pivot burn), 1-opt incumbent polish on real tree incumbents + (CbcHeuristicLocal-style binary flips via warm dual re-solves), face walk (least-degradation dive along the LP-optimal face — proves optimality outright on degenerate alternate-optima instances), RENS, feasibility pump, batch rounding dive, RINS-lite - with exponential failure backoff; heuristic bursts are time-boxed + warm-started from the node's own basis, with exponential failure + backoff; heuristic bursts are time-boxed (root MaxTime/3, deeper MaxTime/8) and skipped once the incumbent sits near the node bound, so the tree keeps its budget. - **Anti-degeneracy**: EXPAND (Gill et al.) on the primal ratio test — @@ -106,9 +118,9 @@ It fails on any failure not listed in `testdata/pulp_known_failures.txt`. ## Missing vs. real CBC -- **Cut families beyond GMI, probing and pairwise TwoMIR**: no knapsack - cover, clique, flow-cover, or lift-and-project cuts; no cuts below - the root. A sound multi-row c-MIR generator (equality-chain +- **Cut families beyond GMI, probing, single-row MIR and pairwise + TwoMIR**: no knapsack cover, clique, flow-cover, or lift-and-project + cuts; no cuts below the root. A sound multi-row c-MIR generator (equality-chain aggregation with exact variable-bound substitution, property-tested) exists in `mip/cmir.go` but stays unwired: measured on the target instances it separates nothing that GMI + probing have not already @@ -117,10 +129,11 @@ It fails on any failure not listed in `testdata/pulp_known_failures.txt`. violated row (plus the revisit tabu) instead of dual steepest edge, and runs unperturbed. The Clp-style engine (DSE, Harris ratio test, cost perturbation — `simplex/dual2.go`) is property-tested but gated - off as measured-negative on the target instances. FTRAN/BTRAN are - always dense-vector solves; Clp's hypersparse triangular solves are - not implemented (measured ~20% result density here, so the payoff is - bounded). + off as measured-negative on the target instances. BTRAN is hypersparse + (activation-graph guarded, bitwise-identical to the dense path) for + mostly-zero right-hand sides and FTRAN skips zero pivots naturally, + but there is no Clp-style hypersparse bookkeeping across the eta file + (measured ~20% result density bounds the further payoff). - **No CglPreProcess-style reductions**: presolve tightens bounds and coefficients but never eliminates rows/columns, so node LPs stay large (real CBC works on a ~4x smaller reduced model for these instances). diff --git a/mip/cmir.go b/mip/cmir.go index 580d6a6..1386dcc 100644 --- a/mip/cmir.go +++ b/mip/cmir.go @@ -135,6 +135,38 @@ func linkVar(p *problem.Problem, r1, r2 int) int { // divisor/complement search, and adds violated cuts. Returns rows added. var cmirDbg struct{ windows, bails, derives, noBin, f0rej, violrej int } +// rowMIRCuts derives single-row MIR cuts (CglMixedIntegerRounding2-style) +// from the first nRows original rows: no aggregation, just VUB/bound +// substitution and the divisor search on each row in both directions. +func (mo *Model) rowMIRCuts(x []float64, nRows int) int { + p := mo.P + vubs := mo.detectVUBs() + cmirDbg = struct{ windows, bails, derives, noBin, f0rej, violrej int }{} + added := 0 + agg := map[int]float64{} + for ri := 0; ri < nRows && added < maxCutsPer; ri++ { + r := &p.Rows[ri] + clear(agg) + for pos, j := range r.Idx { + agg[j] += r.Coef[pos] + } + rl, ru := r.Bounds() + if ru < problem.Inf { + if mo.cmirDerive(agg, ru, 1, vubs, x) { + added++ + } + } + if rl > -problem.Inf { + if mo.cmirDerive(agg, rl, -1, vubs, x) { + added++ + } + } + } + debugf("rowmir: vubs=%d derives=%d noBin=%d f0rej=%d violrej=%d added=%d", + len(vubs), cmirDbg.derives, cmirDbg.noBin, cmirDbg.f0rej, cmirDbg.violrej, added) + return added +} + func (mo *Model) cmirCuts(x []float64) int { vubs := mo.detectVUBs() if len(vubs) == 0 { diff --git a/mip/mip.go b/mip/mip.go index d201742..11ddea1 100644 --- a/mip/mip.go +++ b/mip/mip.go @@ -95,6 +95,9 @@ type Model struct { // scratch for node-level bound propagation (propagatedChild) propLB, propUB, propL0, propU0 []float64 + // debug-only pivot attribution per heuristic (SOLVER_DEBUG) + dbgFW, dbgFP, dbgRINS, dbgDive, dbgNodeCold int64 + // per-column pseudocosts: observed bound gain per unit fraction psUp, psDn []float64 psUpN, psDnN []int @@ -176,6 +179,11 @@ func SolveRelaxation(p *problem.Problem) Result { } func (m *Model) Solve() Result { + t0 := time.Now() + mark := func(phase string) { + st := m.LP.Stats + debugf("phase: %s at %v (solves %d, pivots %d)", phase, time.Since(t0).Round(time.Millisecond), st.Solves, st.Phase1+st.Phase2+st.Dual) + } m.live = nil deadline := time.Time{} if m.Limits.MaxTime > 0 { @@ -201,15 +209,21 @@ func (m *Model) Solve() Result { var startObj float64 var startX, startAct, startRC, startPrice []float64 haveStart := false + // the mipstart block's root solve doubles as cut round 0's when + // reducedCostFix left the LP untouched (preSolvedPivots keeps the cut + // pivot ledger identical to solving it inside the loop) + var preSolved *simplex.State + var preSolvedPivots int64 if len(m.MIPStart) == len(m.P.Cols) && len(m.P.SOSs) == 0 { + preBase := m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual if status, st, _ := m.LP.ColdSolve(); status == simplex.Optimal { + preSolvedPivots = m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual - preBase rx, _, rrc, _ := m.LP.Solution(st) robj := m.LP.InternalObjective(st) - if obj, x, act, rc, price, ok := m.completePoint(&node{brCol: -1}, m.MIPStart); ok { - if pObj, px, pAct, pRC, pPrice, better := m.polishIncumbent(&node{brCol: -1}, obj, x, deadline); better { - debugf("polish: mipstart %g -> %g", obj, pObj) - obj, x, act, rc, price = pObj, px, pAct, pRC, pPrice - } + // polishing the trivial start measured as pure waste: it burned + // 45k pivots on 020 while reducedCostFix fixed 0 columns at any + // cutoff far above the optimum; tree heuristics polish for real + if obj, x, act, rc, price, ok := m.completePoint(&node{brCol: -1}, m.MIPStart, st); ok { startObj, startX, startAct, startRC, startPrice = obj, x, act, rc, price haveStart = true debugf("mipstart: pre-cut incumbent obj=%g", obj) @@ -219,12 +233,19 @@ func (m *Model) Solve() Result { m.LP.SetBound(ov.idx, m.P.Cols[ov.idx].LB, m.P.Cols[ov.idx].UB) } m.live = nil + if len(m.rcTouched) == 0 { + preSolved = st // bounds fully restored: st still solves this LP + } } } + mark("mipstart done") // GMI cut rounds tighten the root while its bound still moves; capped // at a fifth of the budget so the tree always gets its time origRows := len(m.P.Rows) + // solved state of the CURRENT m.LP, nil whenever the LP was rebuilt + // after the solve — lets the slack-drop pass skip a duplicate ColdSolve + var rootSt *simplex.State // restart passes inherit the first pass's cuts: go straight to the tree if len(m.P.SOSs) == 0 && !m.SkipProbing { cutDeadline := deadline @@ -236,17 +257,48 @@ func (m *Model) Solve() Result { flat := 0 // bound improvement can pause a round and resume lastBatch := -1 // first row index of the previous round's cuts poisoned := 0 + // pivots are speed-invariant work units: budgeting rounds by pivots + // keeps the cut set deterministic across engine speed changes; the + // wall-clock box stays as a safety cap only + pivots := func() int64 { return m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual } + cutBase := pivots() + if preSolved != nil { + cutBase -= preSolvedPivots // the reused solve's work stays on the ledger + } + cutPivotBudget := int64(10 * m.LP.NumRows()) for round := range maxCutRound { if !cutDeadline.IsZero() && time.Now().After(cutDeadline) { break } - status, st, _ := m.LP.ColdSolve() + if pivots()-cutBase > cutPivotBudget { + // the final batch was never validated by a re-solve: drop + // it, matching the poison path's retraction semantics + if lastBatch >= 0 { + truncateRows(m.P, lastBatch) + stats := m.LP.Stats + m.LP = simplex.Build(m.P) + m.LP.Stats = stats + rootSt = nil + } + debugf("cuts: pivot budget spent after %d rounds", round) + break + } + var status simplex.Status + var st *simplex.State + if round == 0 && preSolved != nil { + status, st = simplex.Optimal, preSolved + } else { + status, st, _ = m.LP.ColdSolve() + } + rootSt = nil if status != simplex.Optimal { // the last batch poisoned the LP (degenerate grind): // retract it and continue with the last good relaxation if lastBatch >= 0 { truncateRows(m.P, lastBatch) + stats := m.LP.Stats m.LP = simplex.Build(m.P) + m.LP.Stats = stats // counters survive the retraction rebuild m.LP.Deadline = cutDeadline debugf("cuts: retracted poison batch, back to %d rows", lastBatch) lastBatch = -1 @@ -257,6 +309,7 @@ func (m *Model) Solve() Result { m.LP.Deadline = deadline break } + rootSt = st obj := m.LP.InternalObjective(st) if round > 0 && obj-prevObj < math.Max(1e-7, 1e-9*math.Abs(obj)) { if flat++; flat >= 2 { @@ -275,12 +328,18 @@ func (m *Model) Solve() Result { gmi := m.gomoryCuts(st) // probing cuts pay off on big fixed-charge instances; on small // ones they poison node re-solves that branching closes anyway - prb := 0 + prb, mir := 0, 0 if m.LP.NumRows() > 1500 { prb = m.probingCuts(x, cutDeadline) + // single-row MIR only seeds rounds 0-1: repeated rounds + // bloat the rows and blur the face the walk needs + if round <= 1 { + mir = m.rowMIRCuts(x, origRows) + } } - added := gmi + prb - debugf("cuts: round %d added %d rows (gmi %d, probing %d, bound %g)", round, added, gmi, prb, obj) + added := gmi + prb + mir + debugf("cuts: round %d added %d rows (gmi %d, probing %d, mir %d, bound %g, pivots %d)", + round, added, gmi, prb, mir, obj, m.LP.Stats.Phase1+m.LP.Stats.Phase2+m.LP.Stats.Dual) if added == 0 { break } @@ -288,14 +347,20 @@ func (m *Model) Solve() Result { m.LP = simplex.Build(m.P) m.LP.Stats = stats // carry pivot counters across rebuilds m.LP.Deadline = cutDeadline + rootSt = nil } m.LP.Deadline = deadline } // keep only cuts tight at the root: slack rows bloat every node re-solve if len(m.P.Rows) > origRows { - if status, st, _ := m.LP.ColdSolve(); status == simplex.Optimal { - _, rowAct, _, _ := m.LP.Solution(st) + if rootSt == nil { + if status, st, _ := m.LP.ColdSolve(); status == simplex.Optimal { + rootSt = st + } + } + if rootSt != nil { + _, rowAct, _, _ := m.LP.Solution(rootSt) if dropped := dropSlackCuts(m.P, origRows, rowAct); dropped > 0 { stats := m.LP.Stats m.LP = simplex.Build(m.P) @@ -306,6 +371,7 @@ func (m *Model) Solve() Result { } } + mark("cuts+drop done") pq := &nodeHeap{{bound: math.Inf(-1), brCol: -1}} heap.Init(pq) @@ -467,7 +533,7 @@ func (m *Model) Solve() Result { // no incumbent yet: caller's MIP start first, then heuristics // (children above already captured their bounds) if !hasIncumbent && nodeCount == 1 && len(m.MIPStart) == len(m.P.Cols) { - if hObj, hx, hAct, hRC, hPrice, ok := m.completePoint(nd, m.MIPStart); ok { + if hObj, hx, hAct, hRC, hPrice, ok := m.completePoint(nd, m.MIPStart, endState); ok { debugf("mipstart: accepted obj=%g", hObj) newIncumbent(hObj, hx, hAct, hRC, hPrice) } else { @@ -493,16 +559,22 @@ func (m *Model) Solve() Result { if hasIncumbent { cutoff = bestInternal } + fwBase := m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual hObj, hx, hAct, hRC, hPrice, ok := m.faceWalk(nd, x, endState, obj, cutoff) + m.dbgFW += m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual - fwBase if ok { debugf("facewalk: integral vertex on face obj=%g", hObj) } else { hObj, hx, hAct, hRC, hPrice, ok = m.rensImprove(nd, x) } if !ok { + fpBase := m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual hObj, hx, hAct, hRC, hPrice, ok = m.feasibilityPump(nd, x, endState) + m.dbgFP += m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual - fpBase if !ok { + dvBase := m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual hObj, hx, hAct, hRC, hPrice, ok = m.diveForIncumbent(nd, x, endState) + m.dbgDive += m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual - dvBase } } m.LP.Deadline = savedDL @@ -515,7 +587,8 @@ func (m *Model) Solve() Result { // LP-complete the rest; accept only strict improvements. // Failures back the frequency off (CBC adaptive frequency) if hasIncumbent && nodeCount%(64< %g", bestInternal, hObj) if pObj, px, pAct, pRC, pPrice, better := m.polishIncumbent(nd, hObj, hx, deadline); better { @@ -526,6 +599,7 @@ func (m *Model) Solve() Result { } else { rinsFails++ } + m.dbgRINS += m.LP.Stats.Phase1 + m.LP.Stats.Phase2 + m.LP.Stats.Dual - rnBase } nd, backtrack = near, far } @@ -549,6 +623,7 @@ func (m *Model) Solve() Result { proven := hasIncumbent && !m.improves(remaining, bestInternal) debugf("mip exit: nodes=%d best=%g remaining=%g proven=%v stopped=%v proofLost=%v heap=%d", nodeCount, bestInternal, remaining, proven, stopped, proofLost, pq.Len()) + debugf("pivot split: facewalk %d, fpump %d, dive %d, rins %d, coldfallbacks %d", m.dbgFW, m.dbgFP, m.dbgDive, m.dbgRINS, m.dbgNodeCold) res := Result{HasIncumbent: hasIncumbent, NodeCount: nodeCount} switch { @@ -683,9 +758,11 @@ func (m *Model) solveNode(nd *node) (status simplex.Status, x, rowAct, rc, price // a degenerate warm start can stall at the iteration cap; a fresh // basis usually solves the same node quickly if status == simplex.IterLimit && (m.LP.Deadline.IsZero() || time.Now().Before(m.LP.Deadline)) { + m.dbgNodeCold++ status, endState, _ = m.LP.ColdSolve() } } else { + m.dbgNodeCold++ status, endState, _ = m.LP.ColdSolve() } if status == simplex.Optimal { @@ -738,7 +815,8 @@ func (m *Model) feasibilityPump(nd *node, x []float64, endState *simplex.State) debugf("fp: integral at iter %d", iter) // polish: real objective, all integers fixed at their values m.LP.SwapCost(orig) - return m.completePoint(nd, curX) + // fp's basis optimized the L1 cost, not the real one: cold here + return m.completePoint(nd, curX, nil) } if prev != nil && slicesEqual(prev, rounded) { debugf("fp: cycle at iter %d", iter) @@ -771,7 +849,7 @@ func slicesEqual(a, b []float64) bool { // rinsImprove fixes integers where the incumbent and the node LP agree and // dives the few disagreeing ones (CBC's RINS neighborhood, LP-approximated). -func (m *Model) rinsImprove(nd *node, x []float64) (obj float64, dx, rowAct, rc, price []float64, ok bool) { +func (m *Model) rinsImprove(nd *node, x []float64, st *simplex.State) (obj float64, dx, rowAct, rc, price []float64, ok bool) { best := m.bestXSnapshot if best == nil { return 0, nil, nil, nil, nil, false @@ -790,7 +868,9 @@ func (m *Model) rinsImprove(nd *node, x []float64) (obj float64, dx, rowAct, rc, v := math.Max(lb, math.Min(ub, bv)) ovs = append(ovs, boundOverride{j, v, v}) } - child := &node{overrides: ovs, bound: nd.bound, depth: nd.depth + 1} + // warm-start from the node's own solved basis like every other child + // solve: the fixings are ordinary bound overrides for the dual repair + child := &node{overrides: ovs, parentState: st, bound: nd.bound, depth: nd.depth + 1} status, cx, _, _, _, _, cState := m.solveNode(child) if status != simplex.Optimal { return 0, nil, nil, nil, nil, false @@ -957,7 +1037,7 @@ func (m *Model) polishIncumbent(nd *node, obj float64, x []float64, deadline tim // completePoint fixes every integer column at its rounded value from point // (a structural vector) and LP-solves the continuous completion. -func (m *Model) completePoint(nd *node, point []float64) (obj float64, x, rowAct, rc, price []float64, ok bool) { +func (m *Model) completePoint(nd *node, point []float64, st *simplex.State) (obj float64, x, rowAct, rc, price []float64, ok bool) { ovs := make([]boundOverride, len(nd.overrides), len(nd.overrides)+len(m.P.Cols)) copy(ovs, nd.overrides) for j, c := range m.P.Cols { @@ -968,7 +1048,7 @@ func (m *Model) completePoint(nd *node, point []float64) (obj float64, x, rowAct v := math.Max(lb, math.Min(ub, math.Round(point[j]))) ovs = append(ovs, boundOverride{j, v, v}) } - child := &node{overrides: ovs, bound: nd.bound, depth: nd.depth + 1} + child := &node{overrides: ovs, parentState: st, bound: nd.bound, depth: nd.depth + 1} status, cx, cAct, cRC, cPrice, cObj, _ := m.solveNode(child) if status != simplex.Optimal { return 0, nil, nil, nil, nil, false @@ -1053,6 +1133,10 @@ func (m *Model) reducedCostFix(rootX, rootRC []float64, rootObj, best float64) { if slack < 0 { return } + nFixed := len(m.rcTouched) + defer func() { + debugf("rcfix: cutoff %g fixed %d cols (total %d)", best, len(m.rcTouched)-nFixed, len(m.rcTouched)) + }() live := make(map[int]bool, len(m.live)) for _, ov := range m.live { live[ov.idx] = true diff --git a/mip/presolve.go b/mip/presolve.go index cb9b1f5..a7575e1 100644 --- a/mip/presolve.go +++ b/mip/presolve.go @@ -289,6 +289,11 @@ func propagate(p *problem.Problem, lb, ub []float64) bool { // infeasible side fixes the binary; otherwise merged implied bounds apply. // Effort is time-boxed: partial results are valid tightenings. func probe(p *problem.Problem, deadline time.Time) { + nFix, nTight, nProbed := 0, 0, 0 + t0 := time.Now() + defer func() { + debugf("probe: %d binaries probed, %d fixed, %d bounds tightened in %v", nProbed, nFix, nTight, time.Since(t0)) + }() n := len(p.Cols) base := make([]float64, 2*n) lb0, ub0 := make([]float64, n), make([]float64, n) @@ -303,6 +308,7 @@ func probe(p *problem.Problem, deadline time.Time) { if !deadline.IsZero() && time.Now().After(deadline) { return } + nProbed++ copy(lb0, base[:n]) copy(ub0, base[n:]) copy(lb1, base[:n]) @@ -317,9 +323,11 @@ func probe(p *problem.Problem, deadline time.Time) { case !feas0: p.Cols[j].LB = 1 base[j] = 1 + nFix++ case !feas1: p.Cols[j].UB = 0 base[n+j] = 0 + nFix++ default: // merged implied bounds only for integer columns: continuous // merges compound propagation drift into a collapse cascade @@ -330,10 +338,12 @@ func probe(p *problem.Problem, deadline time.Time) { if l := math.Min(lb0[i], lb1[i]); l > base[i]+1e-9 { p.Cols[i].LB = l base[i] = l + nTight++ } if u := math.Max(ub0[i], ub1[i]); u < base[n+i]-1e-9 { p.Cols[i].UB = u base[n+i] = u + nTight++ } } } diff --git a/simplex/factor.go b/simplex/factor.go index 5f48ec7..1310be7 100644 --- a/simplex/factor.go +++ b/simplex/factor.go @@ -13,7 +13,7 @@ const ( // triPivot resolves one variable by substitution: basis position pos pivots // on row row with diagonal a. type triPivot struct { - pos, row int + pos, row int32 a float64 } @@ -39,27 +39,44 @@ type factor struct { kRows []int // kernel rows (original indices) kPos []int // kernel basis positions klu *sparseLU // sparse LU of the kernel - rowKIdx []int // row -> kernel row index, -1 otherwise + rowKIdx []int32 // row -> kernel row index, -1 otherwise - // solve scratch, reused across calls (solver is single-threaded); - // per-call make() showed up as ~8% madvise in profiles - xScratch, kScratch []float64 + // btran activation graph in solve order (bwd, kernel, reversed fwd): + // a pivot runs only when its w entry or an earlier written row feeds it + posSeq []int32 // basis position -> pivot seq + rdrHead []int32 // row -> [rdrHead[r], rdrHead[r+1]) into rdrList + rdrList []int32 // pivot seqs whose column reads that row + + ws *factorWS // shared per-LP solve scratch (solver is single-threaded) } -// factorWS holds factorize's reusable working arrays; nil means allocate -// fresh (the solver is single-threaded, so one per LP suffices). +// factorWS holds factorize's reusable working arrays plus the solve-time +// scratch shared by every factor of one LP; nil means allocate fresh (the +// solver is single-threaded, so one per LP suffices). type factorWS struct { rowCount, colCount []int rowActive, colActive []bool rowCols [][]int32 queue []int + + // solve scratch: xScratch/kScratch are per-call temporaries; ySolve + // stays all-zero between btranSparse calls (written-rows restore); + // actGen+gen are the epoch-stamped activation flags + xScratch, kScratch, ySolve []float64 + actGen []int32 + written []int32 + gen int32 } func newFactorWS(m int) *factorWS { + fbuf := make([]float64, 3*m) + ibuf := make([]int32, 2*m+1) return &factorWS{ rowCount: make([]int, m), colCount: make([]int, m), rowActive: make([]bool, m), colActive: make([]bool, m), rowCols: make([][]int32, m), queue: make([]int, 0, m), + xScratch: fbuf[:m], kScratch: fbuf[m : 2*m], ySolve: fbuf[2*m:], + actGen: ibuf[: m+1 : m+1], written: ibuf[m+1:][:0], } } @@ -133,7 +150,7 @@ func factorize(m int, cols []int32, tblRow [][]int32, tblVal [][]float64, ws *fa if math.Abs(a) < pivotTol { continue // leave to the kernel } - f.fwd = append(f.fwd, triPivot{pos, r, a}) + f.fwd = append(f.fwd, triPivot{int32(pos), int32(r), a}) rowActive[r], colActive[pos] = false, false colRow, colVal := tblRow[cols[pos]], tblVal[cols[pos]] for k, rr := range colRow { @@ -183,7 +200,7 @@ func factorize(m int, cols []int32, tblRow [][]int32, tblVal [][]float64, ws *fa if row < 0 { continue } - f.bwd = append(f.bwd, triPivot{pos, row, a}) + f.bwd = append(f.bwd, triPivot{int32(pos), int32(row), a}) rowActive[row], colActive[pos] = false, false for _, p := range rowCols[row] { if colActive[p] { @@ -196,13 +213,13 @@ func factorize(m int, cols []int32, tblRow [][]int32, tblVal [][]float64, ws *fa } // kernel: whatever remains, inverted densely - f.rowKIdx = make([]int, m) + f.rowKIdx = make([]int32, m) for i := range f.rowKIdx { f.rowKIdx[i] = -1 } for r := range m { if rowActive[r] { - f.rowKIdx[r] = len(f.kRows) + f.rowKIdx[r] = int32(len(f.kRows)) f.kRows = append(f.kRows, r) } } @@ -245,15 +262,86 @@ func factorize(m int, cols []int32, tblRow [][]int32, tblVal [][]float64, ws *fa return nil } } - buf := make([]float64, m+k) - f.xScratch, f.kScratch = buf[:m], buf[m:] + f.ws = ws + f.buildBtranGraph() return f } +// buildBtranGraph assigns each pivot its btran solve-order seq and builds +// the row -> reader-pivots CSR used to activate only reachable pivots. +func (f *factor) buildBtranGraph() { + m, nb, nf := f.m, len(f.bwd), len(f.fwd) + kernelSeq := int32(nb) + // one arena for the graph's fixed-size int32 arrays + arena := make([]int32, m+(m+1)) + f.posSeq = arena[:m:m] + head := arena[m:] + + for i, tp := range f.bwd { + f.posSeq[tp.pos] = int32(i) + } + for _, pos := range f.kPos { + f.posSeq[pos] = kernelSeq + } + for i, tp := range f.fwd { + f.posSeq[tp.pos] = int32(nb + 1 + (nf - 1 - i)) + } + + countTri := func(tri []triPivot) { + for _, tp := range tri { + for _, r := range f.tblRow[f.cols[tp.pos]] { + if r != tp.row { + head[r+1]++ + } + } + } + } + countTri(f.fwd) + countTri(f.bwd) + for _, pos := range f.kPos { + for _, r := range f.tblRow[f.cols[pos]] { + if f.rowKIdx[r] < 0 { + head[r+1]++ + } + } + } + for r := range m { + head[r+1] += head[r] + } + // cursorless CSR fill: advance head[r] itself, then shift it back + list := make([]int32, head[m]) + fillTri := func(tri []triPivot, seqOf func(i int) int32) { + for i, tp := range tri { + s := seqOf(i) + for _, r := range f.tblRow[f.cols[tp.pos]] { + if r != tp.row { + list[head[r]] = s + head[r]++ + } + } + } + } + fillTri(f.bwd, func(i int) int32 { return int32(i) }) + fillTri(f.fwd, func(i int) int32 { return int32(nb + 1 + (nf - 1 - i)) }) + for _, pos := range f.kPos { + for _, r := range f.tblRow[f.cols[pos]] { + if f.rowKIdx[r] < 0 { + list[head[r]] = kernelSeq + head[r]++ + } + } + } + for r := m; r > 0; r-- { + head[r] = head[r-1] + } + head[0] = 0 + f.rdrHead, f.rdrList = head, list +} + // ftranFactor solves B*x = v in place: v becomes x indexed by basis // position (x[pos] = multiplier of basic column pos). func (f *factor) ftran(v []float64) { - x := f.xScratch // by basis position + x := f.ws.xScratch // by basis position clear(x) // forward: row singletons for _, tp := range f.fwd { @@ -270,7 +358,7 @@ func (f *factor) ftran(v []float64) { // kernel k := len(f.kRows) if k > 0 { - kv := f.kScratch + kv := f.ws.kScratch[:k] for ki, r := range f.kRows { kv[ki] = v[r] } @@ -305,14 +393,26 @@ func (f *factor) ftran(v []float64) { // btran solves B^T*y = w in place: w is indexed by basis position on entry, // y by row on exit. Ops are the adjoints of ftran's, in reverse order. func (f *factor) btran(w []float64) { - y := f.xScratch // by row + nnz := 0 + for _, v := range w[:f.m] { + if v != 0 { + nnz++ + } + } + // hypersparse pay-off needs a mostly-zero rhs; activation bookkeeping + // costs about one extra column pass per active pivot + if nnz*10 <= f.m { + f.btranSparse(w) + return + } + y := f.ws.xScratch // by row clear(y) // adjoint of backward pass, in forward order for _, tp := range f.bwd { s := w[tp.pos] cr, cv := f.tblRow[f.cols[tp.pos]], f.tblVal[f.cols[tp.pos]] for k, r := range cr { - if int(r) != tp.row { + if r != tp.row { s -= cv[k] * y[r] } } @@ -321,7 +421,7 @@ func (f *factor) btran(w []float64) { // adjoint of kernel: solve K^T y_K = (w_K - cols^T y_known) k := len(f.kRows) if k > 0 { - kw := f.kScratch + kw := f.ws.kScratch[:k] for ki, pos := range f.kPos { s := w[pos] cr, cv := f.tblRow[f.cols[pos]], f.tblVal[f.cols[pos]] @@ -343,7 +443,7 @@ func (f *factor) btran(w []float64) { s := w[tp.pos] cr, cv := f.tblRow[f.cols[tp.pos]], f.tblVal[f.cols[tp.pos]] for k, r := range cr { - if int(r) != tp.row { + if r != tp.row { s -= cv[k] * y[r] } } @@ -352,6 +452,89 @@ func (f *factor) btran(w []float64) { copy(w, y) } +// btranSparse is btran gathering only activation-reachable pivots; skipped +// pivots would compute exactly 0, so results match the dense path bitwise. +func (f *factor) btranSparse(w []float64) { + y := f.ws.ySolve // by row, all-zero on entry (restored at exit) + f.ws.gen++ + gen := f.ws.gen + act := f.ws.actGen + written := f.ws.written[:0] + for pos := range f.m { + if w[pos] != 0 { + act[f.posSeq[pos]] = gen + } + } + mark := func(row int32) { + for _, t := range f.rdrList[f.rdrHead[row]:f.rdrHead[row+1]] { + act[t] = gen + } + } + for i, tp := range f.bwd { + if act[i] != gen { + continue + } + s := w[tp.pos] + cr, cv := f.tblRow[f.cols[tp.pos]], f.tblVal[f.cols[tp.pos]] + for k, r := range cr { + if r != tp.row { + s -= cv[k] * y[r] + } + } + if yr := s / tp.a; yr != 0 { + y[tp.row] = yr + written = append(written, tp.row) + mark(tp.row) + } + } + nb := len(f.bwd) + if k := len(f.kRows); k > 0 && act[nb] == gen { + kw := f.ws.kScratch[:k] + for ki, pos := range f.kPos { + s := w[pos] + cr, cv := f.tblRow[f.cols[pos]], f.tblVal[f.cols[pos]] + for kk, r := range cr { + if f.rowKIdx[r] < 0 { + s -= cv[kk] * y[r] + } + } + kw[ki] = s + } + f.klu.solveT(kw) + for ki, r := range f.kRows { + if yr := kw[ki]; yr != 0 { + y[r] = yr + written = append(written, int32(r)) + mark(int32(r)) + } + } + } + nf := len(f.fwd) + for si := range nf { + if act[nb+1+si] != gen { + continue + } + tp := f.fwd[nf-1-si] + s := w[tp.pos] + cr, cv := f.tblRow[f.cols[tp.pos]], f.tblVal[f.cols[tp.pos]] + for k, r := range cr { + if r != tp.row { + s -= cv[k] * y[r] + } + } + if yr := s / tp.a; yr != 0 { + y[tp.row] = yr + written = append(written, tp.row) + mark(tp.row) + } + } + copy(w, y) + for _, r := range written { + y[r] = 0 + } + f.ws.written = written +} + // applyEtas maps x = Binv_factor*v to Binv_current*v (FTRAN direction). func applyEtas(etas []*eta, x []float64) { for _, e := range etas { diff --git a/simplex/factor_test.go b/simplex/factor_test.go index cdafa75..00a3371 100644 --- a/simplex/factor_test.go +++ b/simplex/factor_test.go @@ -173,3 +173,72 @@ func TestFactorEtaUpdate(t *testing.T) { } } } + +// TestBtranSparseMatchesDense checks the activation-guarded btran against +// the dense path bitwise on sparse right-hand sides (unit vectors included). +func TestBtranSparseMatchesDense(t *testing.T) { + rng := rand.New(rand.NewSource(11)) + for trial := range 40 { + m := 150 + rng.Intn(150) + colRow, colVal := randomBasis(rng, m) + f := factorize(m, identCols(m), colRow, colVal, nil) + if f == nil { + continue + } + for sub := range 20 { + w := make([]float64, m) + nnz := 1 + rng.Intn(m/12) // always below the m/10 cutover + for range nnz { + w[rng.Intn(m)] = rng.NormFloat64() + } + dense := append([]float64(nil), w...) + f.btranSparse(w) + got := append([]float64(nil), w...) + // dense reference: run the tail of btran directly + y := make([]float64, m) + for _, tp := range f.bwd { + s := dense[tp.pos] + cr, cv := f.tblRow[f.cols[tp.pos]], f.tblVal[f.cols[tp.pos]] + for k, r := range cr { + if r != tp.row { + s -= cv[k] * y[r] + } + } + y[tp.row] = s / tp.a + } + if k := len(f.kRows); k > 0 { + kw := make([]float64, k) + for ki, pos := range f.kPos { + s := dense[pos] + cr, cv := f.tblRow[f.cols[pos]], f.tblVal[f.cols[pos]] + for kk, r := range cr { + if f.rowKIdx[r] < 0 { + s -= cv[kk] * y[r] + } + } + kw[ki] = s + } + f.klu.solveT(kw) + for ki, r := range f.kRows { + y[r] = kw[ki] + } + } + for i := len(f.fwd) - 1; i >= 0; i-- { + tp := f.fwd[i] + s := dense[tp.pos] + cr, cv := f.tblRow[f.cols[tp.pos]], f.tblVal[f.cols[tp.pos]] + for k, r := range cr { + if r != tp.row { + s -= cv[k] * y[r] + } + } + y[tp.row] = s / tp.a + } + for i := range m { + if got[i] != y[i] && !(got[i] == 0 && y[i] == 0) { + t.Fatalf("trial %d/%d m=%d row %d: sparse %v dense %v", trial, sub, m, i, got[i], y[i]) + } + } + } + } +} diff --git a/simplex/simplex.go b/simplex/simplex.go index e54ce6d..a94840a 100644 --- a/simplex/simplex.go +++ b/simplex/simplex.go @@ -89,6 +89,9 @@ type LP struct { // run/recomputeBasics per-call vectors, reused (single-threaded) runCost, runY, runA, residual []float64 + // pivot's eta staging, copied into exact-size arrays per eta + etaIdxWS []int32 + etaValWS []float64 // Stats accumulates pivot counts across solves (diagnostics only). Stats struct { @@ -1016,17 +1019,13 @@ func (lp *LP) pivot(st *State, q int, dir float64, a []float64, t float64, leave if t < 0 { t = 0 } - // update all basic values - for i := range lp.m { - if a[i] == 0 { - continue - } - st.value[st.basicOf[i]] -= a[i] * dir * t - } - oldVal := st.value[q] - st.value[q] = oldVal + dir*t - if isFlip { + for i := range lp.m { + if a[i] == 0 { + continue + } + st.value[st.basicOf[i]] -= a[i] * dir * t + } if dir > 0 { st.status[q] = atUpper st.value[q] = lp.ub[q] @@ -1037,6 +1036,25 @@ func (lp *LP) pivot(st *State, q int, dir float64, a []float64, t float64, leave return } + // one pass over alpha updates basic values and stages the eta support; + // exact-size copies follow (etas outlive the pivot, scratch does not) + if cap(lp.etaIdxWS) < lp.m { + lp.etaIdxWS = make([]int32, 0, lp.m) + lp.etaValWS = make([]float64, 0, lp.m) + } + wsIdx, wsVal := lp.etaIdxWS[:0], lp.etaValWS[:0] + for i, v := range a[:lp.m] { + if v == 0 { + continue + } + st.value[st.basicOf[i]] -= v * dir * t + if math.Abs(v) > etaDropTol { + wsIdx = append(wsIdx, int32(i)) + wsVal = append(wsVal, v) + } + } + st.value[q] += dir * t + leaving := st.basicOf[leaveRow] // snap the leaving variable exactly onto the bound it reached to avoid // floating-point drift. @@ -1049,22 +1067,10 @@ func (lp *LP) pivot(st *State, q int, dir float64, a []float64, t float64, leave st.value[leaving] = ub } - // product-form update: alpha is exactly the eta for this basis change; - // counting first sizes the arrays exactly (etas outlive the pivot) - nnz := 0 - for _, v := range a { - if math.Abs(v) > etaDropTol { - nnz++ - } - } - idx := make([]int32, 0, nnz) - val := make([]float64, 0, nnz) - for i, v := range a { - if math.Abs(v) > etaDropTol { - idx = append(idx, int32(i)) - val = append(val, v) - } - } + idx := make([]int32, len(wsIdx)) + val := make([]float64, len(wsVal)) + copy(idx, wsIdx) + copy(val, wsVal) st.etas = append(st.etas, &eta{r: leaveRow, idx: idx, val: val, ar: a[leaveRow]}) st.basicOf[leaveRow] = q st.status[q] = basic diff --git a/simplex/sparselu.go b/simplex/sparselu.go index b2a4802..d62eab5 100644 --- a/simplex/sparselu.go +++ b/simplex/sparselu.go @@ -24,6 +24,16 @@ type sparseLU struct { func luFactorize(k int, colIdx [][]int32, colVal [][]float64) *sparseLU { rowIdx := make([][]int32, k) rowVal := make([][]float64, k) + cnt := make([]int, k) + for c := range k { + for _, r := range colIdx[c] { + cnt[r]++ + } + } + for r := range k { + rowIdx[r] = make([]int32, 0, cnt[r]) + rowVal[r] = make([]float64, 0, cnt[r]) + } for c := range k { for t, r := range colIdx[c] { rowIdx[r] = append(rowIdx[r], int32(c)) @@ -48,6 +58,7 @@ func luFactorize(k int, colIdx [][]int32, colVal [][]float64) *sparseLU { } acc := make([]float64, k) + mark := make([]bool, k) // pivot-row cols seen in the current row's scan touched := make([]int32, 0, k) for step := range k { @@ -120,40 +131,36 @@ func luFactorize(k int, colIdx [][]int32, colVal [][]float64) *sparseLU { lu.lIdx[step] = append(lu.lIdx[step], int32(rr)) lu.lVal[step] = append(lu.lVal[step], mult) - nIdx := make([]int32, 0, len(rowIdx[rr])+len(touched)) - nVal := make([]float64, 0, len(rowIdx[rr])+len(touched)) - for t, c := range rowIdx[rr] { + // rewrite the row in place: the write index never passes the + // read index, and fill-in only appends after the scan + row, val := rowIdx[rr], rowVal[rr] + w := 0 + for t, c := range row { if c == pc || colUsed[c] { continue } - v := rowVal[rr][t] + v := val[t] if a := acc[c]; a != 0 { v -= mult * a + mark[c] = true } if math.Abs(v) > etaDropTol { - nIdx = append(nIdx, c) - nVal = append(nVal, v) + row[w], val[w] = c, v + w++ } } - // fill-in: pivot-row cols absent from row rr + row, val = row[:w], val[:w] + // fill-in: pivot-row cols absent from the original row rr for _, c := range touched { - if a := acc[c]; a != 0 { - found := false - for _, cc := range rowIdx[rr] { - if cc == c { - found = true - break - } - } - if !found { - if v := -mult * a; math.Abs(v) > etaDropTol { - nIdx = append(nIdx, c) - nVal = append(nVal, v) - } + if a := acc[c]; a != 0 && !mark[c] { + if v := -mult * a; math.Abs(v) > etaDropTol { + row = append(row, c) + val = append(val, v) } } + mark[c] = false } - rowIdx[rr], rowVal[rr] = nIdx, nVal + rowIdx[rr], rowVal[rr] = row, val } for _, c := range touched { acc[c] = 0