mip presolve: worklist-fixpoint propagation, singleton-column elimination

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 <file>` 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.

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
+}

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)
+}