From 337c3d90c55ac450bc9ec19f50cb7915e0df16ff Mon Sep 17 00:00:00 2001 From: andig Date: Mon, 6 Jul 2026 11:01:32 +0200 Subject: [PATCH 1/2] mip: run bound propagation to fixpoint via row worklist Replaces the 4-full-pass cap with a worklist seeded with all rows that re-enqueues only rows touching a tightened column. Reaches the unique fixpoint (bounds only shrink) while doing less work per call: untouched rows are never rescanned. Gate: 13/13 green; 020 6.9s -> 6.5s, 016 0.5s -> 0.2s, 018 1.2s -> 0.9s. Co-Authored-By: Claude Fable 5 --- mip/presolve.go | 195 +++++++++++++++++++++++++----------------------- 1 file changed, 103 insertions(+), 92 deletions(-) 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 } From 5240c4c6613c5d72ef441916872d22f4b6001a88 Mon Sep 17 00:00:00 2001 From: andig Date: Mon, 6 Jul 2026 11:16:34 +0200 Subject: [PATCH 2/2] mip: singleton-column presolve elimination with exact postsolve MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Costed continuous columns appearing in exactly one row are substituted out before the LP is built: a column whose cost is capped only by its row pins that row at any optimum (x = (b - rest)/a, cost folded onto the row's remaining columns), and one whose row never resists its cost sits at its bound. Folding can re-classify same-row siblings, so the pass iterates to fixpoint. Solutions, duals, reduced costs and row activities are reconstructed exactly in the caller's original space. On the evcc golden instances this finds nothing — their 616 singletons are penalty slacks (cost fights the row: a max(0, .) term no linear presolve can remove); an earlier census misread them as eliminable by ignoring ObjSense. The transform fires on minimize-form models with revenue-style singletons, e.g. via the cbc CLI. Gate: 13/13 green, trajectories bit-identical. Co-Authored-By: Claude Fable 5 --- README.md | 17 +-- mip/eliminate.go | 237 ++++++++++++++++++++++++++++++++++++++++++ mip/eliminate_test.go | 110 ++++++++++++++++++++ mip/mip.go | 14 +++ 4 files changed, 372 insertions(+), 6 deletions(-) create mode 100644 mip/eliminate.go create mode 100644 mip/eliminate_test.go 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 }