diff --git a/README.md b/README.md index 3c3ee2c..eee8251 100644 --- a/README.md +++ b/README.md @@ -67,9 +67,12 @@ It fails on any failure not listed in `testdata/pulp_known_failures.txt`. 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). +- **Presolve**: activity-based bound tightening run to fixpoint via a + row worklist, big-M coefficient tightening for binaries, + CglProbing-style binary probing (infeasibility fixing plus + integer-only merged implied bounds), and singleton-column elimination + (costed continuous singletons that pin their row or sit at a bound + are substituted out and reconstructed exactly at postsolve). - **Cuts**: Gomory mixed-integer cuts at the root, with support/dynamism hygiene, and retraction (with retries) of batches that degrade the LP numerically; rounds are budgeted in pivots (speed-invariant work @@ -134,9 +137,11 @@ It fails on any failure not listed in `testdata/pulp_known_failures.txt`. 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). +- **Partial CglPreProcess-style reductions**: singleton columns are + eliminated, but rows never are, and the evcc instances' singletons are + penalty slacks (cost fights the row — a `max(0, ·)` term no linear + presolve can remove), so node LPs stay large on them (real CBC works + on a ~4x smaller reduced model via row aggregation). - **`-mips ` warm start is parsed but not wired** to `Model.MIPStart`, so `warmStart=True` in PuLP buys nothing yet. - **No multi-threaded search.** `-threads N` is accepted, ignored. diff --git a/mip/eliminate.go b/mip/eliminate.go new file mode 100644 index 0000000..a2da8a0 --- /dev/null +++ b/mip/eliminate.go @@ -0,0 +1,237 @@ +// Singleton-column elimination (CglPreProcess-style): costed continuous +// columns appearing in one row are substituted out before the LP is built. +package mip + +import ( + "math" + + "cbcgo/problem" +) + +type elimKind byte + +const ( + elimTight elimKind = iota // cost pins the row: x = (b - rest)/a + elimFixed // row never blocks: x sits at its bound +) + +type elimRecord struct { + col, row int + kind elimKind + a float64 // column's coefficient in its row + b float64 // elimTight: row bound the column pins + val float64 // elimFixed: fixed value + obj float64 // cost at elimination time (dual postsolve shift) +} + +// reduction maps a reduced problem back to the caller's original one. +type reduction struct { + orig *problem.Problem + records []elimRecord // in elimination order; postsolve walks it backwards + colMap []int // original col index -> reduced index, -1 if eliminated +} + +// eliminateSingletons returns a reduced copy of p (p itself is untouched) +// with eligible singletons substituted out; (nil, nil) when nothing applies. +func eliminateSingletons(p *problem.Problem) (*problem.Problem, *reduction) { + if len(p.SOSs) != 0 { + return nil, nil + } + inf := problem.Inf + obj := make([]float64, len(p.Cols)) + for j := range p.Cols { + obj[j] = p.Cols[j].Obj + } + rlb := make([]float64, len(p.Rows)) + rub := make([]float64, len(p.Rows)) + touched := make([]bool, len(p.Rows)) + for ri := range p.Rows { + rlb[ri], rub[ri] = p.Rows[ri].Bounds() + } + elim := make([]bool, len(p.Cols)) + var records []elimRecord + nTight, nFixed := 0, 0 + + // folding a cost onto row siblings can re-classify them, so iterate + for changed := true; changed; { + changed = false + for j := range p.Cols { + c := &p.Cols[j] + if elim[j] || c.Integer || len(c.Idx) != 1 { + continue + } + a := c.Coef[0] + ri := c.Idx[0] + cm := obj[j] * p.ObjSense + if math.Abs(a) < 1e-9 || cm == 0 { + continue + } + d := 1.0 // objective-improving direction for x_j (minimize sense) + if cm > 0 { + d = -1 + } + rowBlocks := (a*d > 0 && rub[ri] < inf) || (a*d < 0 && rlb[ri] > -inf) + colBlocks := (d > 0 && c.UB < inf) || (d < 0 && c.LB > -inf) + switch { + case rowBlocks && !colBlocks: + // any optimum pins the row: substitute x = (b - rest)/a and + // carry the column's remaining bound onto the row + b := rub[ri] + if a*d < 0 { + b = rlb[ri] + } + lo, hi := c.LB, c.UB + if lo <= -inf { + lo = math.Inf(-1) + } + if hi >= inf { + hi = math.Inf(1) + } + nlb, nub := b-a*hi, b-a*lo + if a < 0 { + nlb, nub = b-a*lo, b-a*hi + } + if math.IsInf(nlb, -1) { + nlb = -inf + } + if math.IsInf(nub, 1) { + nub = inf + } + if nlb <= -inf && nub >= inf { + continue // fully free column: row would become vacuous + } + f := obj[j] / a + r := &p.Rows[ri] + for k, jj := range r.Idx { + if jj != j && !elim[jj] { + obj[jj] -= f * r.Coef[k] + } + } + records = append(records, elimRecord{col: j, row: ri, kind: elimTight, a: a, b: b, obj: obj[j]}) + rlb[ri], rub[ri] = nlb, nub + touched[ri], elim[j], changed = true, true, true + nTight++ + case !rowBlocks && !colBlocks: + continue // unbounded ray: leave it for the solver to report + case !rowBlocks: + // the row never resists the cost direction: x sits at its bound + v := c.UB + if d < 0 { + v = c.LB + } + if rlb[ri] > -inf { + rlb[ri] -= a * v + } + if rub[ri] < inf { + rub[ri] -= a * v + } + records = append(records, elimRecord{col: j, row: ri, kind: elimFixed, a: a, val: v, obj: obj[j]}) + touched[ri], elim[j], changed = true, true, true + nFixed++ + } + } + } + if len(records) == 0 { + return nil, nil + } + debugf("eliminate: %d singleton cols removed (%d tight, %d fixed) of %d", len(records), nTight, nFixed, len(p.Cols)) + + q := problem.New() + q.Name, q.ObjSense = p.Name, p.ObjSense + colMap := make([]int, len(p.Cols)) + for j := range p.Cols { + if elim[j] { + colMap[j] = -1 + continue + } + c := &p.Cols[j] + colMap[j] = q.AddCol(c.Name, c.LB, c.UB, obj[j], c.Integer, nil, nil) + } + for ri := range p.Rows { + r := &p.Rows[ri] + idx := make([]int, 0, len(r.Idx)) + coef := make([]float64, 0, len(r.Idx)) + for k, jj := range r.Idx { + if colMap[jj] >= 0 { + idx = append(idx, colMap[jj]) + coef = append(coef, r.Coef[k]) + } + } + nri := q.AddRow(r.Name, idx, coef, r.Sense, r.RHS) + nr := &q.Rows[nri] + nr.HasRange, nr.Range = r.HasRange, r.Range + if touched[ri] { + setRowBounds(nr, rlb[ri], rub[ri]) + } + } + return q, &reduction{orig: p, records: records, colMap: colMap} +} + +// setRowBounds rewrites a row's sense/rhs/range to represent [lb, ub]. +func setRowBounds(r *problem.Row, lb, ub float64) { + inf := problem.Inf + r.HasRange, r.Range = false, 0 + switch { + case lb == ub: + r.Sense, r.RHS = problem.EQ, lb + case lb > -inf && ub < inf: + r.Sense, r.RHS = problem.LE, ub + r.HasRange, r.Range = true, ub-lb + case ub < inf: + r.Sense, r.RHS = problem.LE, ub + default: + r.Sense, r.RHS = problem.GE, lb + } +} + +// shrinkX maps an original-space point onto the reduced column space. +func (red *reduction) shrinkX(x []float64) []float64 { + out := make([]float64, 0, len(red.colMap)) + for j, nj := range red.colMap { + if nj >= 0 { + out = append(out, x[j]) + } + } + return out +} + +// expand rewrites a reduced-space Result in the original column space, +// reconstructing eliminated columns, duals and the pinned row activities. +func (red *reduction) expand(res *Result) { + if res.X == nil { + return + } + n := len(red.orig.Cols) + x := make([]float64, n) + rc := make([]float64, n) + for j, nj := range red.colMap { + if nj >= 0 { + x[j] = res.X[nj] + if nj < len(res.ReducedCost) { + rc[j] = res.ReducedCost[nj] + } + } + } + act, price := res.RowActivity, res.RowPrice + for i := len(red.records) - 1; i >= 0; i-- { + rec := &red.records[i] + c := &red.orig.Cols[rec.col] + switch rec.kind { + case elimTight: + v := (rec.b - act[rec.row]) / rec.a + x[rec.col] = math.Min(math.Max(v, c.LB), c.UB) + act[rec.row] = rec.b + rc[rec.col] = -rec.a * price[rec.row] + price[rec.row] += rec.obj / rec.a + case elimFixed: + x[rec.col] = rec.val + act[rec.row] += rec.a * rec.val + rc[rec.col] = rec.obj - rec.a*price[rec.row] + } + } + obj := 0.0 + for j := range red.orig.Cols { + obj += red.orig.Cols[j].Obj * x[j] + } + res.X, res.ReducedCost, res.Obj = x, rc, obj +} diff --git a/mip/eliminate_test.go b/mip/eliminate_test.go new file mode 100644 index 0000000..c44b717 --- /dev/null +++ b/mip/eliminate_test.go @@ -0,0 +1,110 @@ +package mip + +import ( + "math" + "testing" + + "cbcgo/problem" +) + +// solveReducedAndExpand solves the reduced problem and maps the result back. +func solveReducedAndExpand(t *testing.T, q *problem.Problem, red *reduction) Result { + t.Helper() + res := New(q).Solve() + if res.Status != Optimal { + t.Fatalf("reduced solve status = %v", res.Status) + } + red.expand(&res) + return res +} + +// checkDualIdentity asserts rc_j == c_j - sum_i price_i * a_ij in the +// original space for every column. +func checkDualIdentity(t *testing.T, p *problem.Problem, res Result) { + t.Helper() + for j := range p.Cols { + c := &p.Cols[j] + want := c.Obj + for k, ri := range c.Idx { + want -= res.RowPrice[ri] * c.Coef[k] + } + if math.Abs(res.ReducedCost[j]-want) > 1e-7 { + t.Errorf("col %s: rc = %g, want %g", c.Name, res.ReducedCost[j], want) + } + } +} + +func TestEliminateTightSingleton(t *testing.T) { + // min -2e + g: e in [0,inf) only in e+g<=10, so any optimum pins the row + p := problem.New() + e := p.AddCol("e", 0, problem.Inf, -2, false, nil, nil) + g := p.AddCol("g", 0, 5, 1, true, nil, nil) + r1 := p.AddRow("cap", []int{e, g}, []float64{1, 1}, problem.LE, 10) + r2 := p.AddRow("dem", []int{g}, []float64{1}, problem.GE, 2) + + q, red := eliminateSingletons(p) + if red == nil || len(q.Cols) != 1 || red.colMap[e] != -1 { + t.Fatalf("expected e eliminated: %+v", red) + } + if got := q.Cols[red.colMap[g]].Obj; got != 3 { + t.Fatalf("folded g cost = %g, want 3", got) + } + res := solveReducedAndExpand(t, q, red) + if math.Abs(res.Obj-(-14)) > 1e-7 || math.Abs(res.X[e]-8) > 1e-7 || math.Abs(res.X[g]-2) > 1e-7 { + t.Fatalf("obj=%g x=%v, want obj=-14 x=[8 2]", res.Obj, res.X) + } + if math.Abs(res.RowActivity[r1]-10) > 1e-7 || math.Abs(res.RowActivity[r2]-2) > 1e-7 { + t.Fatalf("row activity = %v, want [10 2]", res.RowActivity) + } + checkDualIdentity(t, p, res) +} + +func TestEliminateFixedSingleton(t *testing.T) { + // min -x + s: s in [0,3] costs but its <= row never pushes it up: s=0 + p := problem.New() + x := p.AddCol("x", 0, 7, -1, true, nil, nil) + s := p.AddCol("s", 0, 3, 1, false, nil, nil) + r1 := p.AddRow("cap", []int{x, s}, []float64{1, 1}, problem.LE, 10) + + q, red := eliminateSingletons(p) + if red == nil || len(q.Cols) != 1 || red.colMap[s] != -1 { + t.Fatalf("expected s eliminated: %+v", red) + } + if len(red.records) != 1 || red.records[0].kind != elimFixed || red.records[0].val != 0 { + t.Fatalf("expected fixed-at-0 record: %+v", red.records) + } + res := solveReducedAndExpand(t, q, red) + if math.Abs(res.Obj-(-7)) > 1e-7 || math.Abs(res.X[x]-7) > 1e-7 || res.X[s] != 0 { + t.Fatalf("obj=%g x=%v, want obj=-7 x=[7 0]", res.Obj, res.X) + } + if math.Abs(res.RowActivity[r1]-7) > 1e-7 { + t.Fatalf("row activity = %v, want [7]", res.RowActivity) + } + checkDualIdentity(t, p, res) +} + +func TestEliminatePartialChain(t *testing.T) { + // two exports share one row; folding e2's cost flips e1 into a bounded + // penalty that must stay in the problem + p := problem.New() + e1 := p.AddCol("e1", 0, 4, -3, false, nil, nil) + e2 := p.AddCol("e2", 0, problem.Inf, -2, false, nil, nil) + g := p.AddCol("g", 0, 2, 5, true, nil, nil) + r1 := p.AddRow("cap", []int{e1, e2, g}, []float64{1, 1, 1}, problem.LE, 10) + + q, red := eliminateSingletons(p) + if red == nil || len(q.Cols) != 2 || red.colMap[e2] != -1 || red.colMap[e1] < 0 { + t.Fatalf("expected only e2 eliminated: %+v", red) + } + if got := q.Cols[red.colMap[e1]].Obj; got != -1 { + t.Fatalf("folded e1 cost = %g, want -1", got) + } + res := solveReducedAndExpand(t, q, red) + if math.Abs(res.Obj-(-24)) > 1e-7 || math.Abs(res.X[e1]-4) > 1e-7 || math.Abs(res.X[e2]-6) > 1e-7 || res.X[g] != 0 { + t.Fatalf("obj=%g x=%v, want obj=-24 x=[4 6 0]", res.Obj, res.X) + } + if math.Abs(res.RowActivity[r1]-10) > 1e-7 { + t.Fatalf("row activity = %v, want [10]", res.RowActivity) + } + checkDualIdentity(t, p, res) +} diff --git a/mip/mip.go b/mip/mip.go index 11ddea1..918e34d 100644 --- a/mip/mip.go +++ b/mip/mip.go @@ -88,6 +88,7 @@ type Model struct { Limits Limits MIPStart []float64 // optional structural start point; ints get fixed SkipProbing bool // restart passes re-derive identical probe facts + red *reduction // singleton elimination; nil when none applied live []boundOverride // bounds currently applied to LP; see solveNode rcTouched []int // columns tightened by reducedCostFix bestXSnapshot []float64 // incumbent X for the RINS neighborhood @@ -180,6 +181,10 @@ func SolveRelaxation(p *problem.Problem) Result { func (m *Model) Solve() Result { t0 := time.Now() + // restart calls pass an original-space MIP start; map it down + if m.red != nil && len(m.MIPStart) == len(m.red.orig.Cols) { + m.MIPStart = m.red.shrinkX(m.MIPStart) + } 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) @@ -200,6 +205,12 @@ func (m *Model) Solve() Result { } probe(m.P, probeDeadline) presolve(m.P) + if q, red := eliminateSingletons(m.P); red != nil { + m.P, m.red = q, red + if len(m.MIPStart) == len(red.orig.Cols) { + m.MIPStart = red.shrinkX(m.MIPStart) + } + } m.LP = simplex.Build(m.P) m.LP.Deadline = deadline } @@ -654,6 +665,9 @@ func (m *Model) Solve() Result { } res.Obj = bestInternal * m.P.ObjSense } + if m.red != nil { + m.red.expand(&res) + } return res } diff --git a/mip/presolve.go b/mip/presolve.go index a7575e1..12472b5 100644 --- a/mip/presolve.go +++ b/mip/presolve.go @@ -175,112 +175,123 @@ func setCoef(p *problem.Problem, ri, k int, v float64) { } } -// propagate tightens the working bound slices from row activity ranges, -// iterating to a fixpoint; reports false when some row proves infeasible. +// propagate tightens the working bounds via a row worklist run to fixpoint +// (order-independent: bounds only shrink); false when a row is infeasible. func propagate(p *problem.Problem, lb, ub []float64) bool { inf := problem.Inf - for range 4 { - changed := false - for ri := range p.Rows { - r := &p.Rows[ri] - rlb, rub := r.Bounds() - var minSum, maxSum float64 - var minInf, maxInf int - for k, j := range r.Idx { - a := r.Coef[k] - l, u := lb[j], ub[j] - if l <= -inf { - l = math.Inf(-1) - } - if u >= inf { - u = math.Inf(1) - } - lo, hi := a*l, a*u - if a < 0 { - lo, hi = hi, lo - } - if math.IsInf(lo, -1) { - minInf++ + nr := len(p.Rows) + inQ := make([]bool, nr) + queue := make([]int, nr) + for ri := range queue { + queue[ri] = ri + inQ[ri] = true + } + // the 1e-9 improvement floor guarantees termination; the cap only + // guards zeno chains, and a capped exit is still a valid tightening + for done := 0; len(queue) > 0 && done < 64*nr; done++ { + ri := queue[0] + queue = queue[1:] + inQ[ri] = false + r := &p.Rows[ri] + rlb, rub := r.Bounds() + var minSum, maxSum float64 + var minInf, maxInf int + for k, j := range r.Idx { + a := r.Coef[k] + l, u := lb[j], ub[j] + if l <= -inf { + l = math.Inf(-1) + } + if u >= inf { + u = math.Inf(1) + } + lo, hi := a*l, a*u + if a < 0 { + lo, hi = hi, lo + } + if math.IsInf(lo, -1) { + minInf++ + } else { + minSum += lo + } + if math.IsInf(hi, 1) { + maxInf++ + } else { + maxSum += hi + } + } + // row-level infeasibility against the activity range + scale := math.Max(1, math.Max(math.Abs(minSum), math.Abs(maxSum))) + if minInf == 0 && rub < inf && minSum > rub+1e-7*scale { + return false + } + if maxInf == 0 && rlb > -inf && maxSum < rlb-1e-7*scale { + return false + } + for k, j := range r.Idx { + a := r.Coef[k] + if a == 0 { + continue + } + l, u := lb[j], ub[j] + lf, uf := l, u + if lf <= -inf { + lf = math.Inf(-1) + } + if uf >= inf { + uf = math.Inf(1) + } + lo, hi := a*lf, a*uf + if a < 0 { + lo, hi = hi, lo + } + omin, omax := math.Inf(-1), math.Inf(1) + if minInf == 0 { + omin = minSum - lo + } else if minInf == 1 && math.IsInf(lo, -1) { + omin = minSum + } + if maxInf == 0 { + omax = maxSum - hi + } else if maxInf == 1 && math.IsInf(hi, 1) { + omax = maxSum + } + // derived bounds are rounded OUTWARD by the row's error + // scale: inward drift compounds along equality chains + out := 1e-9 * scale / math.Max(math.Abs(a), 1e-12) + nl, nu := l, u + if rub < inf && !math.IsInf(omin, -1) { + if a > 0 { + nu = math.Min(nu, (rub-omin)/a+out) } else { - minSum += lo + nl = math.Max(nl, (rub-omin)/a-out) } - if math.IsInf(hi, 1) { - maxInf++ + } + if rlb > -inf && !math.IsInf(omax, 1) { + if a > 0 { + nl = math.Max(nl, (rlb-omax)/a-out) } else { - maxSum += hi + nu = math.Min(nu, (rlb-omax)/a+out) } } - // row-level infeasibility against the activity range - scale := math.Max(1, math.Max(math.Abs(minSum), math.Abs(maxSum))) - if minInf == 0 && rub < inf && minSum > rub+1e-7*scale { - return false + if p.Cols[j].Integer { + s := 1e-7 * math.Max(1, math.Max(math.Abs(nl), math.Abs(nu))) + nl = math.Ceil(nl - s) + nu = math.Floor(nu + s) } - if maxInf == 0 && rlb > -inf && maxSum < rlb-1e-7*scale { + if nl > nu+1e-7*math.Max(1, math.Abs(nl)) { return false } - for k, j := range r.Idx { - a := r.Coef[k] - if a == 0 { - continue - } - l, u := lb[j], ub[j] - lf, uf := l, u - if lf <= -inf { - lf = math.Inf(-1) - } - if uf >= inf { - uf = math.Inf(1) - } - lo, hi := a*lf, a*uf - if a < 0 { - lo, hi = hi, lo - } - omin, omax := math.Inf(-1), math.Inf(1) - if minInf == 0 { - omin = minSum - lo - } else if minInf == 1 && math.IsInf(lo, -1) { - omin = minSum - } - if maxInf == 0 { - omax = maxSum - hi - } else if maxInf == 1 && math.IsInf(hi, 1) { - omax = maxSum - } - // derived bounds are rounded OUTWARD by the row's error - // scale: inward drift compounds along equality chains - out := 1e-9 * scale / math.Max(math.Abs(a), 1e-12) - nl, nu := l, u - if rub < inf && !math.IsInf(omin, -1) { - if a > 0 { - nu = math.Min(nu, (rub-omin)/a+out) - } else { - nl = math.Max(nl, (rub-omin)/a-out) - } - } - if rlb > -inf && !math.IsInf(omax, 1) { - if a > 0 { - nl = math.Max(nl, (rlb-omax)/a-out) - } else { - nu = math.Min(nu, (rlb-omax)/a+out) + if nl > l+1e-9 || nu < u-1e-9 { + lb[j], ub[j] = math.Max(l, nl), math.Min(u, nu) + for _, rr := range p.Cols[j].Idx { + if !inQ[rr] { + inQ[rr] = true + queue = append(queue, rr) } } - if p.Cols[j].Integer { - s := 1e-7 * math.Max(1, math.Max(math.Abs(nl), math.Abs(nu))) - nl = math.Ceil(nl - s) - nu = math.Floor(nu + s) - } - if nl > nu+1e-7*math.Max(1, math.Abs(nl)) { - return false - } - if nl > l+1e-9 || nu < u-1e-9 { - lb[j], ub[j] = math.Max(l, nl), math.Min(u, nu) - changed = true - } } } - if !changed { - return true - } } return true }