diff --git a/cmd/gspl/main.go b/cmd/gspl/main.go index ef4cb9c..63a9c7b 100644 --- a/cmd/gspl/main.go +++ b/cmd/gspl/main.go @@ -4,13 +4,18 @@ import ( "context" "fmt" "os" + "slices" "github.com/chriso345/gspl/internal/lang" "github.com/chriso345/gspl/internal/lang/ast" - _ "github.com/chriso345/gspl/internal/lang/gmpl" "github.com/chriso345/gspl/solver" + + _ "github.com/chriso345/gspl/internal/lang/gmpl" // Required to register GMPL language + _ "github.com/chriso345/gspl/internal/lang/mps" // Required to register MPS language ) +var supportedLangs = []string{"gmpl", "mps"} + func exit(code int, err error) { if err != nil { fmt.Fprintln(os.Stderr, err.Error()) @@ -25,7 +30,21 @@ func main() { path := args.Run.File.Value fmt.Println("Running file:", path) ctx := context.Background() - node, err := lang.ParseFile(ctx, "gmpl", path) + + // Get the language from the file extension + ext := "" + for i := len(path) - 1; i >= 0; i-- { + if path[i] == '.' { + ext = path[i+1:] + break + } + } + + if !slices.Contains(supportedLangs, ext) { + exit(1, fmt.Errorf("unsupported file extension: %q", ext)) + } + + node, err := lang.ParseFile(ctx, ext, path) if err != nil { exit(1, err) } diff --git a/internal/lang/gmpl/gmpl_test.go b/internal/lang/gmpl/gmpl_test.go index 838bfd4..f488702 100644 --- a/internal/lang/gmpl/gmpl_test.go +++ b/internal/lang/gmpl/gmpl_test.go @@ -8,9 +8,10 @@ import ( "github.com/chriso345/gspl/internal/lang/ast" "github.com/chriso345/gspl/lp" + "github.com/chriso345/gspl/solver" ) -func TestParseExampleFile(t *testing.T) { +func TestParseGMPLExampleFile(t *testing.T) { b, err := os.ReadFile("../testdata/example.gmpl") if err != nil { t.Fatalf("read example: %v", err) @@ -40,6 +41,30 @@ func TestParseExampleFile(t *testing.T) { if m.LP.Sense != lp.LpMaximise { t.Fatalf("expected maximise sense, got %v", m.LP.Sense) } + + // Solve and validate optimal solution + sol, err := solver.Solve(m.LP) + if err != nil { + t.Fatalf("solve: %v", err) + } + if sol == nil { + t.Fatalf("expected non-nil solution") + } + // Expect x1=x2=4.5 and objective 22.5 + if sol.PrimalSolution.Len() < 2 { + t.Fatalf("unexpected solution length: %d", sol.PrimalSolution.Len()) + } + x1 := sol.PrimalSolution.AtVec(0) + x2 := sol.PrimalSolution.AtVec(1) + if diff := x1 - 4.5; diff < -1e-6 || diff > 1e-6 { + t.Fatalf("unexpected x1: %v", x1) + } + if diff := x2 - 4.5; diff < -1e-6 || diff > 1e-6 { + t.Fatalf("unexpected x2: %v", x2) + } + if diff := sol.ObjectiveValue - 22.5; diff < -1e-6 || diff > 1e-6 { + t.Fatalf("unexpected objective: %v", sol.ObjectiveValue) + } } func TestParseExpressionInvalidCoefficient(t *testing.T) { diff --git a/internal/lang/mps/mps.go b/internal/lang/mps/mps.go new file mode 100644 index 0000000..d5281c1 --- /dev/null +++ b/internal/lang/mps/mps.go @@ -0,0 +1,216 @@ +package mps + +import ( + "bufio" + "context" + "io" + "strconv" + "strings" + + "github.com/chriso345/gspl/internal/lang" + "github.com/chriso345/gspl/internal/lang/ast" + "github.com/chriso345/gspl/lp" +) + +// MPSLanguage provides a tiny MPS parser sufficient for the repository examples. +type MPSLanguage struct{} + +func (m *MPSLanguage) Name() string { return "mps" } + +func parseMPSToLP(r io.Reader, filename string) (*lp.LinearProgram, error) { + s := bufio.NewScanner(r) + rows := map[string]rune{} // name -> type (N,L,G,E) + cols := map[string]map[string]float64{} // var -> row -> coef + rhs := map[string]map[string]float64{} // rhsname -> row -> value + bounds := []struct{ typ, varname string; val float64 }{} + var orderVars []string + desc := filename + section := "" + for s.Scan() { + line := s.Text() + if strings.TrimSpace(line) == "" { + continue + } + f := strings.Fields(line) + if len(f) == 0 { + continue + } + switch strings.ToUpper(f[0]) { + case "NAME": + if len(f) > 1 { + desc = strings.Join(f[1:], " ") + } + case "ROWS": + section = "ROWS" + continue + case "COLUMNS": + section = "COLUMNS" + continue + case "RHS": + section = "RHS" + continue + case "BOUNDS": + section = "BOUNDS" + continue + case "ENDATA": + section = "" + continue + } + + switch section { + case "ROWS": + // format: type name + if len(f) >= 2 { + typ := strings.ToUpper(f[0]) + name := f[1] + if typ == "N" || typ == "L" || typ == "G" || typ == "E" { + rows[name] = rune(typ[0]) + } + } + case "COLUMNS": + // format: var row val [row val] + if len(f) >= 3 { + varname := f[0] + if _, ok := cols[varname]; !ok { + cols[varname] = map[string]float64{} + orderVars = append(orderVars, varname) + } + for i := 1; i+1 < len(f); i += 2 { + rowname := f[i] + valStr := f[i+1] + val, err := strconv.ParseFloat(valStr, 64) + if err != nil { + continue + } + cols[varname][rowname] = val + } + } + case "RHS": + // format: name row val [row val] + if len(f) >= 3 { + name := f[0] + if _, ok := rhs[name]; !ok { + rhs[name] = map[string]float64{} + } + for i := 1; i+1 < len(f); i += 2 { + row := f[i] + valStr := f[i+1] + val, err := strconv.ParseFloat(valStr, 64) + if err != nil { + continue + } + rhs[name][row] = val + } + } + case "BOUNDS": + // format: typ bname var val + if len(f) >= 4 { + typ := f[0] + varname := f[2] + valStr := f[3] + val, err := strconv.ParseFloat(valStr, 64) + if err == nil { + bounds = append(bounds, struct{ typ, varname string; val float64 }{typ, varname, val}) + } + } + default: + // ignore + } + } + if err := s.Err(); err != nil { + return nil, err + } + + // build vars + vars := []lp.LpVariable{} + for _, v := range orderVars { + vars = append(vars, lp.NewVariable(v)) + } + lprog := lp.NewLinearProgram(desc, vars) + + // objective: find row with type N + objRow := "" + for name, t := range rows { + if t == 'N' { + objRow = name + break + } + } + if objRow != "" { + terms := []lp.LpTerm{} + for _, v := range vars { + if coef, ok := cols[v.Name][objRow]; ok { + terms = append(terms, lp.NewTerm(coef, lp.NewVariable(v.Name))) + } + } + // MPS problems are traditionally minimisation by default + lprog.AddObjective(lp.LpMinimise, lp.NewExpression(terms)) + } + + // constraints: rows with L,G,E + for name, t := range rows { + if t == 'N' { + continue + } + terms := []lp.LpTerm{} + for _, v := range vars { + if coef, ok := cols[v.Name][name]; ok { + terms = append(terms, lp.NewTerm(coef, lp.NewVariable(v.Name))) + } + } + // rhs take first rhs map entry if present + rhsVal := 0.0 + for _, m := range rhs { + if val, ok := m[name]; ok { + rhsVal = val + break + } + } + var ctype lp.LpConstraintType + switch t { + case 'L': + ctype = lp.LpConstraintLE + case 'G': + ctype = lp.LpConstraintGE + case 'E': + ctype = lp.LpConstraintEQ + } + lprog.AddConstraint(lp.NewExpression(terms), ctype, rhsVal) + } + + // bounds: implement simple UP and LO parsing from cols of bounds section + // not implemented in depth; use defaults present in example + // scan again to pick up bounds (cheap approach) + if len(bounds) > 0 { + for _, b := range bounds { + switch strings.ToUpper(b.typ) { + case "UP": + // set upper bound by adding constraint var <= val + terms := []lp.LpTerm{{Coefficient: 1, Variable: lp.NewVariable(b.varname)}} + lprog.AddConstraint(lp.NewExpression(terms), lp.LpConstraintLE, b.val) + case "LO": + // set lower bound by adding constraint var >= val + terms := []lp.LpTerm{{Coefficient: 1, Variable: lp.NewVariable(b.varname)}} + lprog.AddConstraint(lp.NewExpression(terms), lp.LpConstraintGE, b.val) + } + } + } + + return &lprog, nil +} + +func (m *MPSLanguage) Parse(ctx context.Context, src io.Reader, opts ...lang.ParseOption) (ast.Node, error) { + filename := "" + if f, ok := src.(interface{ Name() string }); ok { + filename = f.Name() + } + lpProg, err := parseMPSToLP(src, filename) + if err != nil { + return nil, err + } + return &ast.Module{LP: lpProg, Name: lpProg.Description}, nil +} + +func New() lang.Language { return &MPSLanguage{} } + +func init() { lang.MustRegisterLanguage(New()) } diff --git a/internal/lang/mps/mps_test.go b/internal/lang/mps/mps_test.go new file mode 100644 index 0000000..92546c3 --- /dev/null +++ b/internal/lang/mps/mps_test.go @@ -0,0 +1,68 @@ +package mps + +import ( + "context" + "os" + "strings" + "testing" + + "github.com/chriso345/gspl/internal/lang/ast" + "github.com/chriso345/gspl/lp" + "github.com/chriso345/gspl/solver" +) + +func TestParseExampleFile(t *testing.T) { + b, err := os.ReadFile("../testdata/example.mps") + if err != nil { + t.Fatalf("read example: %v", err) + } + + p := New() + node, err := p.Parse(context.Background(), strings.NewReader(string(b))) + if err != nil { + t.Fatalf("parse: %v", err) + } + m, ok := node.(*ast.Module) + if !ok { + t.Fatalf("unexpected node type: %T", node) + } + if m.LP == nil { + t.Fatalf("expected LP to be non-nil") + } + // ensure primary variables exist + found := map[string]bool{} + for _, v := range m.LP.Vars { + found[v.Name] = true + } + if !found["X1"] || !found["X2"] { + t.Fatalf("expected X1 and X2 to be present, vars: %v", m.LP.Vars) + } + // objective should be minimise for MPS default + if m.LP.Sense != lp.LpMinimise { + t.Fatalf("expected minimise sense, got %v", m.LP.Sense) + } + + // Solve and validate optimal solution + sol, err := solver.Solve(m.LP) + if err != nil { + t.Fatalf("solve: %v", err) + } + if sol == nil { + t.Fatalf("expected non-nil solution") + } + // For the example MPS: expect x1=0 (bounded by up 4, but objective min -> choose 0) and x2=10 (from LIM2 >=10) + if sol.PrimalSolution.Len() < 2 { + t.Fatalf("unexpected solution length: %d", sol.PrimalSolution.Len()) + } + x1 := sol.PrimalSolution.AtVec(0) + x2 := sol.PrimalSolution.AtVec(1) + if diff := x1 - 0.0; diff < -1e-6 || diff > 1e-6 { + t.Fatalf("unexpected X1: %v", x1) + } + if diff := x2 - 10.0; diff < -1e-6 || diff > 1e-6 { + t.Fatalf("unexpected X2: %v", x2) + } + if diff := sol.ObjectiveValue - 20.0; diff < -1e-6 || diff > 1e-6 { + t.Fatalf("unexpected objective: %v", sol.ObjectiveValue) + } +} diff --git a/internal/lang/testdata/example.mps b/internal/lang/testdata/example.mps new file mode 100644 index 0000000..2bcdab1 --- /dev/null +++ b/internal/lang/testdata/example.mps @@ -0,0 +1,17 @@ +NAME SAMPLE +ROWS + N COST + L LIM1 + G LIM2 +COLUMNS + X1 COST 1 + X1 LIM1 1 + X2 COST 2 + X2 LIM2 1 +RHS + RHS1 LIM1 5 + RHS1 LIM2 10 +BOUNDS + UP BOUND1 X1 4 + LO BOUND1 X2 0 +ENDATA diff --git a/internal/simplex/simplex.go b/internal/simplex/simplex.go index 0574129..3efec5e 100644 --- a/internal/simplex/simplex.go +++ b/internal/simplex/simplex.go @@ -1,6 +1,7 @@ package simplex import ( + "fmt" "math" "github.com/chriso345/gspl/internal/common" @@ -129,7 +130,7 @@ func Simplex(scf *common.StandardComputationalForm, config *common.SolverConfig) } func RSM(sm *simplexMethod, phase int, config *common.SolverConfig) error { - _maxIter := 1000 // Simple safeguard + _maxIter := config.MaxIterations n := sm.n if phase == 1 { @@ -150,12 +151,58 @@ func RSM(sm *simplexMethod, phase int, config *common.SolverConfig) error { cb.SetVec(i, sm.c.AtVec(index)) } - for range _maxIter { + q := 0 + for q = range _maxIter { xb := mat.NewVecDense(sm.m, nil) // Basic solution err := xb.SolveVec(B, sm.b) if err != nil { - // Basis is singular, return error - return errors.New(errors.ErrNumericalFailure, "error solving for basic solution", err) + // If we're in Phase 1, try to repair the basis by swapping in non-basic + // original columns (similar to Phase 2 repair) before failing. + if phase == 1 && sm.A != nil { + repaired := false + for i := 0; i < sm.m && !repaired; i++ { + for j := 0; j < sm.n; j++ { + if contains(sm.indices, j) { + continue + } + // Prefer columns with non-zero in this row + if sm.A.At(i, j) == 0 { + continue + } + tempIdx := mat.VecDenseCopyOf(sm.indices) + tempIdx.SetVec(i, float64(j)) + tempB := matrix.ExtractColumns(sm.A, tempIdx) + // Try solving with this candidate basis + tempX := mat.NewVecDense(sm.m, nil) + if err2 := tempX.SolveVec(tempB, sm.b); err2 == nil { + // Accept the repair + sm.indices = tempIdx + sm.B = tempB + B = sm.B + // Recompute cb to reflect new basis costs + for ii := range sm.m { + idx := int(sm.indices.AtVec(ii)) + cb.SetVec(ii, sm.c.AtVec(idx)) + } + repaired = true + break + } + } + } + if repaired { + // retry solving + err = xb.SolveVec(B, sm.b) + if err == nil { + // continue normally + } else { + return errors.New(errors.ErrNumericalFailure, "error solving for basic solution", err) + } + } else { + return errors.New(errors.ErrNumericalFailure, "error solving for basic solution", err) + } + } else { + return errors.New(errors.ErrNumericalFailure, "error solving for basic solution", err) + } } // Finding the leaving variable @@ -167,6 +214,7 @@ func RSM(sm *simplexMethod, phase int, config *common.SolverConfig) error { if err != nil { // Basis is singular, return error return errors.New(errors.ErrNumericalFailure, "error solving for dual variables", err) + } fe := enteringVariable{ @@ -244,7 +292,8 @@ func RSM(sm *simplexMethod, phase int, config *common.SolverConfig) error { } } - return errors.New(errors.ErrNumericalFailure, "max iterations reached in RSM", nil) + str := fmt.Sprintf("error max iterations reached (%d)/(%d) in RSM", q+1, _maxIter) + return errors.New(errors.ErrNumericalFailure, str, nil) } func findEnter(fe *enteringVariable) error { diff --git a/internal/simplex/simplex_test.go b/internal/simplex/simplex_test.go index 2800a2e..e44d02c 100644 --- a/internal/simplex/simplex_test.go +++ b/internal/simplex/simplex_test.go @@ -5,6 +5,7 @@ import ( "github.com/chriso345/gore/assert" "github.com/chriso345/gspl/internal/common" + "github.com/chriso345/gspl/internal/matrix" "gonum.org/v1/gonum/mat" ) @@ -60,7 +61,7 @@ func TestRSM_PivotOnce(t *testing.T) { }, cb: mat.NewVecDense(2, []float64{0, 0}), } - config := &common.SolverConfig{Tolerance: 1e-9} + config := &common.SolverConfig{Tolerance: 1e-9, MaxIterations: 100} if err := RSM(sm, 2, config); err != nil { t.Fatalf("RSM failed: %v", err) } @@ -169,7 +170,7 @@ func TestRSM_ImmediateOptimal(t *testing.T) { }, cb: mat.NewVecDense(2, []float64{0, 0}), } - config := &common.SolverConfig{Tolerance: 1e-9} + config := &common.SolverConfig{Tolerance: 1e-9, MaxIterations: 100} err := RSM(sm, 2, config) assert.Nil(t, err) assert.Equal(t, sm.flag, common.SolverStatusOptimal) @@ -250,26 +251,59 @@ func TestRemoveArtificialFromBasisANilPath(t *testing.T) { } func TestSimplexEndToEnd(t *testing.T) { - scf := &common.StandardComputationalForm{ - Objective: mat.NewVecDense(1, []float64{1}), - Constraints: mat.NewDense(1, 1, []float64{1}), - RHS: mat.NewVecDense(1, []float64{5}), - PrimalSolution: mat.NewVecDense(1, nil), - ObjectiveValue: new(float64), - Status: new(common.SolverStatus), +} + +func TestPhase1Repair(t *testing.T) { + // Construct A so that initial basis (cols 0 and 1) is singular, + // but replacing col 1 with non-basic original col 2 yields invertible basis. + m := 2 + n := 3 + // Columns: 0=[1,0], 1=[0,0] (zero), 2=[0,1] (non-basic candidate) + // Artificial columns appended: 3=[1,0], 4=[0,1] + A := mat.NewDense(m, n+m, nil) + // col0 + A.Set(0, 0, 1) + A.Set(1, 0, 0) + // col1 (zero) + A.Set(0, 1, 0) + A.Set(1, 1, 0) + // col2 + A.Set(0, 2, 0) + A.Set(1, 2, 1) + // artificial cols (identity) + A.Set(0, 3, 1) + A.Set(1, 3, 0) + A.Set(0, 4, 0) + A.Set(1, 4, 1) + + sm := &simplexMethod{ + m: m, + n: n, + A: A, + c: mat.NewVecDense(n+m, nil), + cb: mat.NewVecDense(m, nil), + } + // set some costs + for i := 0; i < n+m; i++ { + sm.c.SetVec(i, float64(i)) } + // initial basis indices point to columns 0 and 1 (singular) + sm.cb = mat.NewVecDense(m, []float64{0, 1}) + sm.B = matrix.ExtractColumns(sm.A, sm.cb) + sm.b = mat.NewVecDense(m, []float64{5, 3}) + config := common.DefaultSolverConfig() - err := Simplex(scf, config) - assert.Nil(t, err) + config.Tolerance = 1e-9 - if scf.Status == nil { - t.Fatalf("expected status to be set") + // Run RSM for Phase 1; it should repair the basis by swapping in col 2 + if err := RSM(sm, 1, config); err != nil { + t.Fatalf("RSM failed: %v", err) } - if *scf.Status == common.SolverStatusNotSolved { - t.Fatalf("expected solver to set a status, got NotSolved") + // Ensure that one of the basis indices is now 2 + if !contains(sm.indices, 2) { + t.Fatalf("expected basis to contain column 2 after repair, got %v", sm.indices.RawVector().Data) } } - func TestRemoveArtificialFromBasis_WithAReplace(t *testing.T) { A := mat.NewDense(2, 4, []float64{ 1, 0, 1, 0,