feat: Implement Pareto Frontier MOP Solvers

bindings/c/gspl.go

@@ -13,12 +13,14 @@ import (
 	"github.com/chriso345/gspl/internal/common"
 	"github.com/chriso345/gspl/lp"
 	"github.com/chriso345/gspl/solver"
+	"github.com/chriso345/gspl/solver/mop"
 )
 
 var lastError string
 
 // Ensure references to solver option helpers are present so the bindings
 // export-check (heuristic) can detect coverage of those symbols.
+// FIXME: These should be implemented fully eventually.
 var (
 	_ = solver.NewSolverConfig
 	_ = solver.WithTolerance
@@ -30,6 +32,9 @@ var (
 	_ = solver.WithBranch
 	_ = solver.WithHeuristic
 	_ = solver.WithCut
+
+	_ = mop.WithEpsilonSteps
+	_ = mop.WithWeightedSums
 )
 
 type program struct {
@@ -213,7 +218,7 @@ type solverConfig struct {
 
 type mopSolution struct {
 	vals []float64
-	sol  *solver.MopSolution
+	sol  *mop.MopSolution
 }
 
 // cSolution is a C-facing snapshot of a solver solution mapping original
@@ -277,7 +282,7 @@ func gspl_solve_lexicographic(prog C.GSPL_Handle, _ C.GSPL_Handle) C.GSPL_Handle
 	}
 
 	// solver config currently ignored by this wrapper
-	mop, err := solver.SolveLexicographic(&p.lp)
+	mop, err := mop.SolveLexicographic(&p.lp)
 	if err != nil {
 		lastError = err.Error()
 		return 0
@@ -297,7 +302,7 @@ func gspl_solve_pareto(prog C.GSPL_Handle, _ C.GSPL_Handle) C.GSPL_Handle {
 	}
 
 	// solver config currently ignored by this wrapper
-	mops, err := solver.SolvePareto(&p.lp)
+	mops, err := mop.SolvePareto(&p.lp, mop.ParetoWeightedSum)
 	if err != nil {
 		lastError = err.Error()
 		return 0

internal/common/solver_config.go

@@ -20,7 +20,11 @@ type SolverConfig struct {
 	Heuristic      HeuristicFunc
 	Cut            CutFunc
 
-	// Not Yet Implemented
+	// MOP Specific Options
+	EpsilonSteps int
+	WeightedSums [][]float64
+
+	// Configuration Not Yet Implemented
 	Threads int
 
 	Debug bool
@@ -39,6 +43,9 @@ func DefaultSolverConfig() *SolverConfig {
 		Heuristic:      nil, // Default heuristic defined in `brancher`
 		Cut:            nil, // Default cutting planes defined in `brancher`
 
+		EpsilonSteps: 10,
+		WeightedSums: nil,
+
 		Threads: 0, // 0 means use all available cores
 
 		Debug: false,

lp/formulate.go

@@ -9,6 +9,14 @@ import (
 func (lp *LinearProgram) AddObjective(sense LpSense, expr LpExpression) {
 	lp.Sense = sense
 
+	// If an objective already exists, move it to SecondaryObjectives
+	if lp.Objective != nil {
+		if lp.SecondaryObjectives == nil {
+			lp.SecondaryObjectives = []*mat.VecDense{}
+		}
+		lp.SecondaryObjectives = append(lp.SecondaryObjectives, mat.VecDenseCopyOf(lp.Objective))
+	}
+
 	// Create the objective function vector
 	lp.Objective = mat.NewVecDense(len(lp.Vars), nil)
 

solver/mop/doc.go

@@ -0,0 +1,7 @@
+// Package mop provides helpers for solving multi-objective linear programs.
+//
+// It offers several algorithms for generating Pareto-optimal solution sets and
+// supports lexicographic optimisation where objectives are solved in priority
+// order. The package is intended to operate on lp.LinearProgram instances from
+// the lp package and to use the solver package for single-objective solves.
+package mop

solver/mop_solver.go → solver/mop/lexicographic.go

@@ -1,25 +1,18 @@
-package solver
+package mop
 
 import (
-	"github.com/chriso345/gspl/internal/common"
 	"github.com/chriso345/gspl/internal/errors"
 	"github.com/chriso345/gspl/lp"
+	"github.com/chriso345/gspl/solver"
 	"gonum.org/v1/gonum/mat"
 )
 
-// MopSolution contains the result of a multi-objective optimization.
-type MopSolution struct {
-	ObjectiveValues []float64
-	PrimalSolution  *mat.VecDense
-	Status          common.SolverStatus
-}
-
 // SolveLexicographic solves the given multi-objective linear program using
 // a lexicographic approach. This method optimises the objectives in a predefined
 // order, ensuring that the optimal value of each objective is achieved before
 // moving on to the next. It returns a single MopSolution representing the optimal
 // solution that respects the priority of objectives.
-func SolveLexicographic(prog *lp.LinearProgram, opts ...SolverOption) (*MopSolution, error) {
+func SolveLexicographic(prog *lp.LinearProgram, opts ...solver.SolverOption) (*MopSolution, error) {
 	if prog.Objective == nil {
 		return nil, errors.New(errors.ErrInvalidInput, "primary objective must be defined", nil)
 	}
@@ -58,7 +51,7 @@ func SolveLexicographic(prog *lp.LinearProgram, opts ...SolverOption) (*MopSolut
 		p.RHS = mat.NewVecDense(1, []float64{0})
 	}
 
-	sol, err := Solve(&p, opts...)
+	sol, err := solver.Solve(&p, opts...)
 	if err != nil {
 		return nil, err
 	}
@@ -109,7 +102,7 @@ func SolveLexicographic(prog *lp.LinearProgram, opts ...SolverOption) (*MopSolut
 			}
 		}
 		p.Objective = sec
-		lastSol, err = Solve(&p, opts...)
+		lastSol, err = solver.Solve(&p, opts...)
 		if err != nil {
 			return nil, err
 		}
@@ -119,11 +112,3 @@ func SolveLexicographic(prog *lp.LinearProgram, opts ...SolverOption) (*MopSolut
 
 	return &MopSolution{ObjectiveValues: objVals, PrimalSolution: lastSol.PrimalSolution, Status: lastSol.Status}, nil
 }
-
-// SolvePareto solves the given multi-objective linear program using a Pareto
-// approach. This method seeks to find a set of solutions that represent the best
-// trade-offs among the objectives, known as the Pareto front. It returns a slice
-// of MopSolution, each representing a non-dominated solution in the objective space.
-func SolvePareto(prog *lp.LinearProgram, opts ...SolverOption) ([]*MopSolution, error) {
-	return nil, nil
-}

solver/mop_solver_test.go → solver/mop/lexicographic_test.go

@@ -1,4 +1,4 @@
-package solver
+package mop
 
 import (
 	"testing"

solver/mop/options.go

@@ -0,0 +1,29 @@
+package mop
+
+import (
+	"github.com/chriso345/gspl/internal/common"
+	"github.com/chriso345/gspl/solver"
+)
+
+type MopOption = solver.SolverOption
+
+// WithEpsilonSteps sets the number of epsilon steps for the epsilon-constraint method.
+func WithEpsilonSteps(steps int) MopOption {
+	return func(cfg *common.SolverConfig) {
+		cfg.EpsilonSteps = steps
+	}
+}
+
+// WithWeightedSums sets the weight vectors for the weighted sum method.
+func WithWeightedSums(weights [][]float64) MopOption {
+	return func(cfg *common.SolverConfig) {
+		cfg.WeightedSums = weights
+	}
+}
+
+// with sets the entire SolverConfig (for internal use only).
+func with(cfg_ common.SolverConfig) MopOption {
+	return func(cfg *common.SolverConfig) {
+		*cfg = cfg_
+	}
+}

solver/mop/pareto.go

@@ -0,0 +1,158 @@
+package mop
+
+import (
+	"github.com/chriso345/gspl/internal/errors"
+	"github.com/chriso345/gspl/lp"
+	"github.com/chriso345/gspl/solver"
+	"gonum.org/v1/gonum/mat"
+)
+
+// SolvePareto solves the given multi-objective linear program using a Pareto
+// approach. This method seeks to find a set of solutions that represent the best
+// trade-offs among the objectives, known as the Pareto front. It returns a slice
+// of MopSolution, each representing a non-dominated solution in the objective space.
+func SolvePareto(
+	prog *lp.LinearProgram,
+	method ParetoMethod,
+	opts ...solver.SolverOption,
+) ([]*MopSolution, error) {
+
+	if prog.Objective == nil {
+		return nil, errors.New(errors.ErrInvalidInput, "primary objective must be defined", nil)
+	}
+
+	p, objs, err := cloneWithObjectives(prog)
+	if err != nil {
+		return nil, err
+	}
+
+	if len(objs) < 2 {
+		return nil, errors.New(errors.ErrInvalidInput, "pareto requires >=2 objectives", nil)
+	}
+
+	switch method {
+	case ParetoWeightedSum:
+		return solveParetoWeightedSum(&p, objs, opts...)
+	case ParetoEpsilonConstraint:
+		return solveParetoEpsilon(&p, objs, opts...)
+	default:
+		return nil, errors.New(errors.ErrInvalidInput, "unknown pareto method", nil)
+	}
+}
+
+func solveParetoWeightedSum(
+	p *lp.LinearProgram,
+	objs []*mat.VecDense,
+	sopts ...solver.SolverOption,
+) ([]*MopSolution, error) {
+	cfg := solver.NewSolverConfig(sopts...)
+
+	n := len(objs)
+	if n == 0 {
+		return nil, nil
+	}
+
+	weights := cfg.WeightedSums
+	if len(weights) == 0 {
+		weights = defaultSimplexWeights(n)
+	}
+
+	results := make([]*MopSolution, 0, len(weights))
+
+	for _, w := range weights {
+		if len(w) != n {
+			return nil, errors.New(errors.ErrInvalidInput, "weight vector dimension mismatch", nil)
+		}
+
+		comb := mat.NewVecDense(objs[0].Len(), nil)
+		for j, obj := range objs {
+			for i := 0; i < obj.Len(); i++ {
+				comb.SetVec(i, comb.AtVec(i)+w[j]*obj.AtVec(i))
+			}
+		}
+
+		p.Objective = comb
+		sol, err := solver.Solve(p, with(*cfg))
+		if err != nil {
+			return nil, err
+		}
+
+		vals := evaluateObjectives(sol.PrimalSolution, objs, p.ObjectiveIsNegated)
+		results = append(results, &MopSolution{
+			ObjectiveValues: vals,
+			PrimalSolution:  sol.PrimalSolution,
+			Status:          sol.Status,
+		})
+	}
+
+	return uniqueSolutions(results), nil
+}
+
+func solveParetoEpsilon(
+	p *lp.LinearProgram,
+	objs []*mat.VecDense,
+	sopts ...solver.SolverOption,
+) ([]*MopSolution, error) {
+	cfg := solver.NewSolverConfig(sopts...)
+
+	if len(objs) < 2 {
+		return nil, errors.New(errors.ErrInvalidInput, "epsilon-constraint requires >=2 objectives", nil)
+	}
+
+	primary := objs[0]
+	others := objs[1:]
+
+	p.Objective = primary
+	base, err := solver.Solve(p, with(*cfg))
+	if err != nil {
+		return nil, err
+	}
+
+	results := []*MopSolution{
+		{
+			ObjectiveValues: evaluateObjectives(base.PrimalSolution, objs, p.ObjectiveIsNegated),
+			PrimalSolution:  base.PrimalSolution,
+			Status:          base.Status,
+		},
+	}
+
+	epsilonSteps := cfg.EpsilonSteps
+	for _, sec := range others {
+		min, max := objectiveBounds(p, sec, sopts...)
+		step := (max - min) / float64(epsilonSteps)
+
+		for i := 1; i <= epsilonSteps; i++ {
+			eps := min + float64(i)*step
+			addEpsilonConstraint(p, sec, eps)
+
+			sol, err := solver.Solve(p, with(*cfg))
+			if err != nil {
+				continue
+			}
+
+			results = append(results, &MopSolution{
+				ObjectiveValues: evaluateObjectives(sol.PrimalSolution, objs, p.ObjectiveIsNegated),
+				PrimalSolution:  sol.PrimalSolution,
+				Status:          sol.Status,
+			})
+		}
+	}
+
+	return uniqueSolutions(results), nil
+}
+
+func evaluateObjectives(x *mat.VecDense, objs []*mat.VecDense, neg bool) []float64 {
+	vals := make([]float64, len(objs))
+	for j, obj := range objs {
+		sum := 0.0
+		for i := 0; i < x.Len() && i < obj.Len(); i++ {
+			c := obj.AtVec(i)
+			if neg {
+				c = -c
+			}
+			sum += c * x.AtVec(i)
+		}
+		vals[j] = sum
+	}
+	return vals
+}