diff --git a/.github/workflows/coverage1.yml b/.github/workflows/coverage1.yml new file mode 100644 index 0000000..dee4272 --- /dev/null +++ b/.github/workflows/coverage1.yml @@ -0,0 +1,40 @@ +name: Code Coverage # The name of the workflow that will appear on Github + +on: + push: + branches: [ main ] + pull_request: + # Allows you to run this workflow manually from the Actions tab + workflow_dispatch: +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ubuntu-latest, macos-latest, windows-latest] + go: [1.23] + permissions: + # Give the default GITHUB_TOKEN write permission to commit and push the + # added or changed files to the repository. + contents: write + steps: + - name: Checkout repository + uses: actions/checkout@v4 + - name: Set up Go + uses: actions/setup-go@v2 + with: + go-version: ${{ matrix.go }} + + - name: Build project + run: go install ./... + + - name: Test Coverage + run: | + go test -v -cover $(go list ./... | grep -v /examples/) -coverprofile coverage.out -coverpkg ./... + go tool cover -func coverage.out -o coverage_analysis.out + + - name: Upload coverage reports to Codecov + uses: codecov/codecov-action@v5 + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + files: ./coverage.out diff --git a/algorithms/algorithm_interface.go b/algorithms/algorithm_interface.go index 526b6e5..46625bb 100644 --- a/algorithms/algorithm_interface.go +++ b/algorithms/algorithm_interface.go @@ -6,5 +6,5 @@ import ( type AlgorithmInterface interface { // Solves the provided optimization problem. - Solve(initialState AlgorithmInternalState) (problem.Solution, error) + Solve(prob problem.OptimizationProblem) (problem.Solution, error) } diff --git a/algorithms/algorithm_type.go b/algorithms/algorithm_type.go index b8db181..30c2570 100644 --- a/algorithms/algorithm_type.go +++ b/algorithms/algorithm_type.go @@ -2,4 +2,4 @@ package algorithms type AlgorithmType int -const TypeNaive AlgorithmType = AlgorithmType(1) +const TypeNaiveTableau AlgorithmType = AlgorithmType(1) diff --git a/algorithms/internal_state.go b/algorithms/internal_state.go deleted file mode 100644 index d131895..0000000 --- a/algorithms/internal_state.go +++ /dev/null @@ -1,61 +0,0 @@ -package algorithms - -import ( - "fmt" - - "github.com/MatProGo-dev/MatProInterface.go/problem" - "github.com/MatProGo-dev/SymbolicMath.go/symbolic" - "github.com/MatProGo-dev/simplex/utils" - "gonum.org/v1/gonum/mat" -) - -type AlgorithmInternalState struct { - BasicVariables []symbolic.Variable - NonBasicVariables []symbolic.Variable - NonBasicValues *mat.VecDense - IterationCount int -} - -/* -NumberOfBasicVariables -Description: - - Returns the number of basic variables in the Internal State. -*/ -func (state *AlgorithmInternalState) NumberOfBasicVariables() int { - return len(state.BasicVariables) -} - -/* -NumberOfNonBasicVariables -Description: - - Returns the number of non-basic variables in the Internal State. -*/ -func (state *AlgorithmInternalState) NumberOfNonBasicVariables() int { - return len(state.NonBasicVariables) -} - -/* -ToTableau -Description: - - Converts the internal state to a tableau representation. -*/ -func (state *AlgorithmInternalState) ToTableau(standardForm *problem.OptimizationProblem) (utils.Tableau, error) { - // Input Processing - if standardForm == nil { - return utils.Tableau{}, fmt.Errorf( - "ToTableau: standardForm cannot be nil", - ) - } - - // Ensure that the problem is a linear program - if !standardForm.IsLinear() { - return utils.Tableau{}, fmt.Errorf( - "ToTableau: the problem is not a linear program", - ) - } - - return utils.GetInitialTableauFrom(standardForm) -} diff --git a/algorithms/naive.go b/algorithms/naive.go deleted file mode 100644 index 073cef5..0000000 --- a/algorithms/naive.go +++ /dev/null @@ -1,206 +0,0 @@ -package algorithms - -import ( - "fmt" - - "github.com/MatProGo-dev/MatProInterface.go/problem" - "github.com/MatProGo-dev/SymbolicMath.go/symbolic" - "github.com/MatProGo-dev/simplex/utils" - "gonum.org/v1/gonum/mat" -) - -type NaiveAlgorithm struct { - ProblemInStandardForm *problem.OptimizationProblem - IterationLimit int -} - -/* -ComputeFeasibleBasicSolution -Description: - - Computes a feasible solution of the BASIC variables - of the optimization problem: - maximize c^T * x - subject to A * x = b - x >= 0 - The introduction of basic and non-basic variables allows us to - rewrite the problem as: - maximize c_B^T * x_B + c_N^T * x_N - subject to A_B * x_B + A_N * x_N = b - x_B >= 0 - x_N >= 0 - Where - A_B is the matrix of coefficients of the basic variables, - A_N is the matrix of coefficients of the non-basic variables, - c_B is the vector of coefficients of the basic variables, - c_N is the vector of coefficients of the non-basic variables, - x_B is the vector of basic variables, - x_N is the vector of non-basic variables, - b is the vector of constants. - - We assume that the value of the non-basic variables is given - (i.e., they are already saved in the state). -*/ -func (algo *NaiveAlgorithm) ComputeFeasibleBasicSolution(state AlgorithmInternalState) (*mat.VecDense, error) { - // Setup - fmt.Printf("Computing feasible solution...\n") - fmt.Printf("Problem: %v\n", algo.ProblemInStandardForm) - fmt.Printf("Tableau: %v\n", algo) - nBasic := state.NumberOfBasicVariables() - - // Collect the matrices of coefficients - A, b, err := algo.ProblemInStandardForm.LinearEqualityConstraintMatrices() - if err != nil { - return nil, err - } - - // Create the matrix of coefficients of the basic variables - N, err := utils.SliceMatrixAccordingToVariableSet( - A, - algo.ProblemInStandardForm.Variables, - state.NonBasicVariables, - ) - if err != nil { - return nil, err - } - - B, err := utils.SliceMatrixAccordingToVariableSet( - A, - algo.ProblemInStandardForm.Variables, - state.BasicVariables, - ) - if err != nil { - return nil, err - } - - // Solve the system of equations - x := mat.NewVecDense(state.NumberOfBasicVariables(), nil) - - // Compute the part that comes from the rhs (b) - // xComponentFromb = B^-1 * b - var BInv *mat.Dense = mat.NewDense(nBasic, nBasic, nil) - BAsDense := B.ToDense() - fmt.Printf("A: %v\n", A) - fmt.Printf("BAsDense: %v\n", BAsDense) - err = BInv.Inverse(&BAsDense) - if err != nil { - return nil, fmt.Errorf("there was an issue inverting the matrix: %v", err) - } - fmt.Printf("BInv: %v\n", BInv) - bAsVecDense := b.ToVecDense() - fmt.Printf("bAsVecDense: %v\n", bAsVecDense) - x.MulVec(BInv, &bAsVecDense) - - fmt.Printf("x: %v\n", x) - - // Compute the part that comes from the non-basic variables - // xComponentFromXNonBasic = B^(-1) * N * x - xComponentFromXNonBasic := mat.NewVecDense(state.NumberOfBasicVariables(), nil) - BN := mat.NewDense(state.NumberOfBasicVariables(), state.NumberOfNonBasicVariables(), nil) - NAsDense := N.ToDense() - BN.Mul(BInv, &NAsDense) - xComponentFromXNonBasic.MulVec(BN, state.NonBasicValues) - x.AddVec(x, xComponentFromXNonBasic) - - return x, nil -} - -/* -ComputeObjectiveFunctionValueWithFeasibleBasicSolution -Description: - - Computes the value of the objective function - of the optimization problem: - maximize c^T * x - subject to A * x = b - x >= 0 - when the feasible solution of the BASIC variables is given as - input and the value of the non-basic variables is given - in the state of the solver. -*/ -func (algo *NaiveAlgorithm) ComputeObjectiveFunctionValueWithFeasibleBasicSolution(state AlgorithmInternalState, xBasic *mat.VecDense) (float64, error) { - // Setup - fmt.Printf("Computing objective function value...\n") - allVars := algo.ProblemInStandardForm.Variables - - // Collect the linear coefficient of the objective function - objectiveExpression := algo.ProblemInStandardForm.Objective.Expression - objectiveAsSE, tf := objectiveExpression.(symbolic.ScalarExpression) - if !tf { - return 0.0, fmt.Errorf("the objective function is not a scalar expression") - } - - // TODO(kwesi): Check that the objective function is a constant or not - - // Compute the linear coefficient of the objective function - c := objectiveAsSE.LinearCoeff(allVars) - - // Split the coefficient into the basic and non-basic variables - cB, err := utils.SliceVectorAccordingToVariableSet( - symbolic.VecDenseToKVector(c), - allVars, - state.BasicVariables, - ) - if err != nil { - return 0.0, err - } - - cN, err := utils.SliceVectorAccordingToVariableSet( - symbolic.VecDenseToKVector(c), - allVars, - state.NonBasicVariables, - ) - if err != nil { - return 0.0, err - } - - // Compute the value of the objective function - // f(x) = c_B^T * x_B + c_N^T * x_N - z := cB.Transpose().Multiply(xBasic).Plus( - cN.Transpose().Multiply(state.NonBasicValues), - ) - - zAsK, tf := z.(symbolic.K) - if !tf { - return 0.0, fmt.Errorf("the objective function is not a scalar expression") - } - - return float64(zAsK), nil -} - -func (algo *NaiveAlgorithm) Solve(initialState AlgorithmInternalState) (problem.Solution, error) { - // Setup - var stateII AlgorithmInternalState = initialState - - // Compute the feasible solution for the current choice of - // basic variables - for iter := 0; iter < algo.IterationLimit; iter++ { - // Compute the feasible Solution of the Basic variables - xBasicII, err := algo.ComputeFeasibleBasicSolution(stateII) - if err != nil { - return problem.Solution{}, err - } - - // Compute the value of the objective function - objII, err := algo.ComputeObjectiveFunctionValueWithFeasibleBasicSolution(stateII, xBasicII) - if err != nil { - return problem.Solution{}, err - } - fmt.Printf("Iteration %d: Basic Solution: %v, Objective Value: %f\n", iter, xBasicII, objII) - - // Update the state of the algorithm - stateII = stateII - - // Check if the solution is optimal - - if iter == algo.IterationLimit { - return problem.Solution{ - Objective: objII, - Status: problem.OptimizationStatus_ITERATION_LIMIT, - }, nil - } - - } - - return problem.Solution{}, nil -} diff --git a/algorithms/tableau/selection/blands_rule.go b/algorithms/tableau/selection/blands_rule.go new file mode 100644 index 0000000..8cf3821 --- /dev/null +++ b/algorithms/tableau/selection/blands_rule.go @@ -0,0 +1,108 @@ +package selection + +import ( + "fmt" + + "github.com/MatProGo-dev/simplex/utils" +) + +type BlandsRule struct{} + +/* +Description: + + SelectEnteringVariable selects an entering variable from the given tableau using the Bland's Rule. + If the result of this function is called: + `idx := SelectEnteringVariable(tableau)` + then `idx` is the index of the entering variable in `tableau.Variables`. + If no entering variable can be found (i.e., the current solution is optimal), + then -1 is returned as the variable index. +*/ +func (br BlandsRule) SelectEnteringVariable(tableau utils.Tableau) int { + // Setup + minIndex := -1 + minValue := 0.0 + + // Get the cost coefficients + costCoefficients := tableau.C() + + // Iterate through the non-basic variables and find + // the coefficient with the smallest index that is negative + for _, nonBasicVarIdx := range tableau.NonBasicVariableIndicies() { + if costCoefficients.AtVec(nonBasicVarIdx) < minValue { + minIndex = nonBasicVarIdx + minValue = costCoefficients.AtVec(nonBasicVarIdx) + } + } + + // If no entering variable is found, return -1 + return minIndex +} + +/* +SelectExitingVariable + +Description: + + SelectExitingVariable selects an exiting variable from the given tableau using the Blands Rule. + + If the result of this function is called: + `idx := SelectExitingVariable(tableau, enteringVarIdx)` + then `idx` is the index of the exiting variable in `tableau.Variables`. + If no exiting variable can be found (i.e., the problem is unbounded), + then -1 is returned as the variable index. +*/ +func (br BlandsRule) SelectExitingVariable(tableau utils.Tableau, enteringVarIdx int) int { + // Setup + minIndex := -1 + minRatio := -1.0 + + // Get the relevant matrices + A := tableau.A() + b := tableau.B() + + // Create the vector of ratios + ratios := make([]float64, tableau.NumberOfConstraints()) + for i := 0; i < tableau.NumberOfConstraints(); i++ { + if A.At(i, enteringVarIdx) > 0 { + ratios[i] = b.AtVec(i) / A.At(i, enteringVarIdx) + } else { + ratios[i] = -1.0 // Indicate that this variable cannot be used + } + } + + // Iterate through the ratios and find the smallest one + for i, ratio := range ratios { + if ratio >= 0 { // Only consider valid ratios + if minRatio < 0 || ratio < minRatio || (ratio == minRatio && tableau.BasicVariableIndicies[i] < tableau.BasicVariableIndicies[minIndex]) { + minRatio = ratio + minIndex = i + } + } + } + + // If no exiting variable is found, return -1 (this is guaranteed by the structure of the algorithm) + if minIndex == -1 { + return -1 + } + + // Return the index of the exiting variable in tableau.Variables + minIndex = tableau.BasicVariableIndicies[minIndex] + return minIndex +} + +func (br BlandsRule) SelectEnteringAndExitingVariables(tableau utils.Tableau) (int, int, error) { + // Select the entering variable + enteringVarIdx := br.SelectEnteringVariable(tableau) + if enteringVarIdx == -1 { + return -1, -1, nil // Optimal solution found, no entering variable + } + + // Select the exiting variable + exitingVarIdx := br.SelectExitingVariable(tableau, enteringVarIdx) + if exitingVarIdx == -1 { + return enteringVarIdx, -1, fmt.Errorf("BlandsRule: No exiting variable found, problem is unbounded (?)") + } + + return enteringVarIdx, exitingVarIdx, nil +} diff --git a/algorithms/tableau/state.go b/algorithms/tableau/state.go index ae3119c..84103cd 100644 --- a/algorithms/tableau/state.go +++ b/algorithms/tableau/state.go @@ -6,6 +6,7 @@ import ( "github.com/MatProGo-dev/MatProInterface.go/problem" "github.com/MatProGo-dev/SymbolicMath.go/symbolic" "github.com/MatProGo-dev/simplex/algorithms" + "github.com/MatProGo-dev/simplex/algorithms/tableau/selection" "github.com/MatProGo-dev/simplex/utils" "gonum.org/v1/gonum/mat" ) @@ -89,10 +90,16 @@ func (state *TableauAlgorithmState) CheckTerminationCondition() (bool, error) { return true, nil } -// func GetStateFromInitialProblem(prob problem.OptimizationProblem) TableauAlgorithmState { -// // Setup +func (state *TableauAlgorithmState) CurrentObjectiveValue() (float64, error) { + // Input Checking + err := state.Check() + if err != nil { + return 0.0, err + } -// } + // Setup + return state.Tableau.D(), nil +} /* GetBasicVariables @@ -142,8 +149,7 @@ Description: Returns the number of constraints as inferred by the size of the tableau. */ func (state *TableauAlgorithmState) NumberOfConstraints() int { - nRows, _ := state.Tableau.AsCompressedMatrix.Dims() - return nRows - 1 + return state.Tableau.NumberOfConstraints() } func (state *TableauAlgorithmState) XBasic() (*mat.VecDense, error) { @@ -237,6 +243,104 @@ func (state *TableauAlgorithmState) GetReducedCostVector() (*mat.VecDense, error return &finalReducedCost, nil } +func (state *TableauAlgorithmState) CalculateNextState() (TableauAlgorithmState, error) { + // Input Checking + err := state.Check() + if err != nil { + return TableauAlgorithmState{}, err + } + + // Select the pivot column and row (i.e., the entering and exiting variables in the tableau) + // Here, we use Bland's Rule to select the entering variable + // TODO(Kwesi): Make other rules available + selectionRule := selection.BlandsRule{} + enteringVarIdx, exitingVarIdx, err := selectionRule.SelectEnteringAndExitingVariables(*state.Tableau) + if err != nil { + return TableauAlgorithmState{}, fmt.Errorf("TableauAlgorithmState: Failed to select entering and exiting variables (%v)", err) + } + + // Create the new tableau + newTab, err := state.Tableau.Pivot(enteringVarIdx, exitingVarIdx) + if err != nil { + return TableauAlgorithmState{}, fmt.Errorf("TableauAlgorithmState: Failed to pivot tableau (%v)", err) + } + + // Create the new state + return TableauAlgorithmState{ + Tableau: &newTab, + IterationCount: state.IterationCount + 1, + }, nil +} + +func (state *TableauAlgorithmState) CalculateOptimalSolution() (mat.VecDense, error) { + // Input Checking + err := state.Check() + if err != nil { + return mat.VecDense{}, err + } + + // Setup + numVars := state.NumberOfVariables() + solutionVec := mat.NewVecDense(numVars, nil) + + // Create a linear system of variables consisting of: + // - The constraints + // - The non-basic variables set to zero + A, b := state.A(), state.B() + numConstraints, _ := A.Dims() + numNonBasic := state.NumberOfVariables() - state.NumberOfConstraints() + + // Augment the A and b matrices with the non-basic variable constraints + AAugmented := mat.NewDense(numConstraints+numNonBasic, numVars, nil) + AAugmented.Copy(A) + bAugmented := mat.NewVecDense(numConstraints+numNonBasic, nil) + bAugmented.CopyVec(b) + + fmt.Println("A: ", mat.Formatted(AAugmented)) + + // Add the non-basic variable constraints + nonBasicVars := state.GetNonBasicVariables() + for ii, v := range nonBasicVars { + fmt.Println("Adding non-basic variable constraint for variable: ", v) + // Find the index of the variable in the tableau + vIdxInTableau, _ := symbolic.FindInSlice(v, state.Tableau.Variables) + fmt.Println("vIdxInTableau: ", vIdxInTableau) + fmt.Println("Targeted row: ", numConstraints+ii) + AAugmented.Set(numConstraints+ii, vIdxInTableau, 1.0) + } + // b is already zero in the new rows, so we don't need to set anything in bAugmented + + // Solve the system of equations + err = solutionVec.SolveVec(AAugmented, bAugmented) + if err != nil { + return mat.VecDense{}, fmt.Errorf("TableauAlgorithmState: Failed to solve for optimal solution (%v)", err) + } + + return *solutionVec, nil +} + +func (state *TableauAlgorithmState) CreateOptimalValuesMap() (map[uint64]float64, error) { + // Input Checking + err := state.Check() + if err != nil { + return nil, err + } + + // Setup + solutionVec, err := state.CalculateOptimalSolution() + if err != nil { + return nil, err + } + + // Create the map + solutionMap := map[uint64]float64{} + for ii, v := range state.Tableau.Variables { + solutionMap[v.ID] = solutionVec.AtVec(ii) + } + + return solutionMap, nil +} + func (state *TableauAlgorithmState) ToSolution(currentStatus problem.OptimizationStatus) problem.Solution { // Construct Solution Map diff --git a/algorithms/tableau/tableau_algo.go b/algorithms/tableau/tableau_algo.go index 68e76e4..a2c25a6 100644 --- a/algorithms/tableau/tableau_algo.go +++ b/algorithms/tableau/tableau_algo.go @@ -4,24 +4,53 @@ import ( "fmt" "github.com/MatProGo-dev/MatProInterface.go/problem" + tableau_termination "github.com/MatProGo-dev/simplex/algorithms/tableau/termination" + "github.com/MatProGo-dev/simplex/utils" "gonum.org/v1/gonum/mat" ) type TableauAlgorithm struct { - ProblemInStandardForm *problem.OptimizationProblem - CurrentSolution *mat.VecDense - IterationLimit int + IterationLimit int } -func (algo *TableauAlgorithm) Solve(initialState TableauAlgorithmState) (problem.Solution, error) { +func (algo *TableauAlgorithm) CheckTerminationConditions(state TableauAlgorithmState) (tableau_termination.TerminationType, error) { + // Input Checking + err := state.Check() + if err != nil { + return tableau_termination.DidNotTerminate, err + } + + // Check If the iteration limit has been reached + if state.IterationCount >= algo.IterationLimit { + return tableau_termination.MaximumIterationsReached, nil + } + + // Check that the reduced costs are all non-negative + if state.Tableau.CanNotBeImproved() { + return tableau_termination.OptimalSolutionFound, nil + } + + return tableau_termination.DidNotTerminate, nil +} + +func (algo *TableauAlgorithm) Solve(prob problem.OptimizationProblem) (problem.Solution, error) { // Setup - var stateII TableauAlgorithmState = initialState - var status problem.OptimizationStatus = problem.OptimizationStatus_INPROGRESS + + // Create initial Tableau state from the problem + initialTableau, err := utils.GetInitialTableauFrom(&prob) + if err != nil { + return problem.Solution{}, fmt.Errorf("there was an issue creating the initial tableau: %v", err) + } + stateII := TableauAlgorithmState{ + Tableau: &initialTableau, + IterationCount: 0, + } // Loop + sol := problem.Solution{} for iter := 0; iter < algo.IterationLimit; iter++ { // Test for Termination - terminated, err := stateII.CheckTerminationCondition() + condition, err := algo.CheckTerminationConditions(stateII) if err != nil { return problem.Solution{}, fmt.Errorf( @@ -31,27 +60,44 @@ func (algo *TableauAlgorithm) Solve(initialState TableauAlgorithmState) (problem ) } - if terminated { + if condition != tableau_termination.DidNotTerminate { + sol.Status = condition.ToOptimizationStatus() + sol.Objective, err = stateII.CurrentObjectiveValue() + if err != nil { + return sol, + fmt.Errorf( + "There was an issue getting the objective value at termination: %v", + err, + ) + } + sol.Values, err = stateII.CreateOptimalValuesMap() + if err != nil { + return sol, + fmt.Errorf( + "There was an issue creating the optimal values map at termination: %v", + err, + ) + } + + // Exit the loop break } - // Compute XB, y and r - _, err = stateII.XBasic() + fmt.Println("Iteration: ", iter) + fmt.Println("Matrix: ", mat.Formatted(stateII.Tableau.AsCompressedMatrix)) + + // Update the state + stateII, err = stateII.CalculateNextState() if err != nil { return problem.Solution{}, fmt.Errorf( - "There was an issue computing the value XBasic() at iteration #%v: %v", + "There was an issue updating the state at iteration %v: %v", iter, err, ) } - // If we reach the limit, then return limit reached as status - if iter == algo.IterationLimit-1 { - status = problem.OptimizationStatus_ITERATION_LIMIT - } - } - return stateII.ToSolution(status), nil + return sol, nil } diff --git a/algorithms/tableau/termination/termination_types.go b/algorithms/tableau/termination/termination_types.go new file mode 100644 index 0000000..0a3bb80 --- /dev/null +++ b/algorithms/tableau/termination/termination_types.go @@ -0,0 +1,24 @@ +package tableau_termination + +import ( + "github.com/MatProGo-dev/MatProInterface.go/problem" +) + +type TerminationType string + +const DidNotTerminate TerminationType = "Did Not Terminate" +const MaximumIterationsReached TerminationType = "Maximum Iterations Reached" +const OptimalSolutionFound TerminationType = "Optimal Solution Found" + +func (tt TerminationType) ToOptimizationStatus() problem.OptimizationStatus { + switch tt { + case DidNotTerminate: + return problem.OptimizationStatus_INPROGRESS + case MaximumIterationsReached: + return problem.OptimizationStatus_ITERATION_LIMIT + case OptimalSolutionFound: + return problem.OptimizationStatus_OPTIMAL + default: + return problem.OptimizationStatus_INPROGRESS + } +} diff --git a/examples/readme.md b/examples/readme.md new file mode 100644 index 0000000..54a67b1 --- /dev/null +++ b/examples/readme.md @@ -0,0 +1,4 @@ +# Examples + +This directory contains a number of examples that demonstrate how to use the simplex +solver. \ No newline at end of file diff --git a/examples/sergiy_butenko1/the_simplex_method_part_ii.go b/examples/sergiy_butenko1/the_simplex_method_part_ii.go new file mode 100644 index 0000000..ba062a9 --- /dev/null +++ b/examples/sergiy_butenko1/the_simplex_method_part_ii.go @@ -0,0 +1,33 @@ +package main + +import ( + "github.com/MatProGo-dev/simplex/simplexSolver" + "github.com/MatProGo-dev/simplex/utils/examples" +) + +func main() { + // Get the example tableau + problem5 := examples.GetTestProblem5() + + // // Print the problem + // println(problem5.String()) + + // Use solver to solve the problem + solver := simplexSolver.New("Simplex Solver Example") + solver.IterationLimit = 100 + + // Solve the problem + solution, err := solver.Solve(*problem5) + if err != nil { + panic(err) + } + + // Print the solution + solutionMessage, _ := solution.Status.ToMessage() + println("Solution Status: ", solutionMessage) + println("Objective Value: ", solution.Objective) + println("Variable Values: ") + for varName, varValue := range solution.Values { + println(" ", varName, ": ", varValue) + } +} diff --git a/simplexSolver/config.go b/simplexSolver/config.go deleted file mode 100644 index 5b050be..0000000 --- a/simplexSolver/config.go +++ /dev/null @@ -1,5 +0,0 @@ -package simplexSolver - -type Configuration struct { - IterationLimit int -} diff --git a/simplexSolver/solver.go b/simplexSolver/solver.go index 923bf72..e0c9fa4 100644 --- a/simplexSolver/solver.go +++ b/simplexSolver/solver.go @@ -4,115 +4,54 @@ import ( "fmt" "github.com/MatProGo-dev/MatProInterface.go/problem" - "github.com/MatProGo-dev/SymbolicMath.go/symbolic" "github.com/MatProGo-dev/simplex/algorithms" - "github.com/MatProGo-dev/simplex/utils" - "gonum.org/v1/gonum/mat" + tableau_algorithm1 "github.com/MatProGo-dev/simplex/algorithms/tableau" ) type SimplexSolver struct { - OriginalProblem *problem.OptimizationProblem - ProblemInStandardForm *problem.OptimizationProblem - InternalState algorithms.AlgorithmInternalState - config Configuration + Name string + IterationLimit int + Algorithm algorithms.AlgorithmType } func New(name string) SimplexSolver { - // Create name for the base problem - baseProblemName := name + " Problem" return SimplexSolver{ - OriginalProblem: problem.NewProblem(baseProblemName + " (Original Problem)"), - ProblemInStandardForm: problem.NewProblem(baseProblemName + " (In Standard Form)"), + Name: name, + IterationLimit: 100, + Algorithm: algorithms.TypeNaiveTableau, } } -func For(problem *problem.OptimizationProblem, configIn Configuration) (SimplexSolver, error) { - // Create a new solver - solver := New(problem.Name + " Solver") - - // Set the original problem - original := problem - original.Name = problem.Name + " (Original Problem)" - solver.OriginalProblem = original - - // Transform the problem into the standard form where all constraints - // are equality constraints - var err error - var slackVariables []symbolic.Variable - solver.ProblemInStandardForm, slackVariables, err = solver.OriginalProblem.ToLPStandardForm1() - if err != nil { - return solver, err - } - - // Initialize Internal Solver State - solver.InternalState, err = solver.InitializeInternalState(slackVariables) - if err != nil { - return solver, err - } - - // Configure the solver - solver.config = configIn - - return solver, nil -} - -// func (solver *SimplexSolver) CurrentStateToTableau() (symbolic.KMatrix, error) { - -// } - -func (solver *SimplexSolver) InitializeInternalState(initialSlackVariables []symbolic.Variable) (algorithms.AlgorithmInternalState, error) { - // Setup - - // Initialize the internal state - out := algorithms.AlgorithmInternalState{ - BasicVariables: initialSlackVariables, - NonBasicVariables: utils.SetDifferenceOfVariables( - solver.ProblemInStandardForm.Variables, - initialSlackVariables, - ), - IterationCount: 0, - } - - nNonBasicVariables := out.NumberOfNonBasicVariables() - out.NonBasicValues = mat.NewVecDense( - nNonBasicVariables, - nil, - ) - - return out, nil -} - func (solver *SimplexSolver) CreateAlgorithm(algoType algorithms.AlgorithmType) (algorithms.AlgorithmInterface, error) { // Setup // Selection Logic switch algoType { - case algorithms.TypeNaive: - return &algorithms.NaiveAlgorithm{ - ProblemInStandardForm: solver.ProblemInStandardForm, - IterationLimit: solver.config.IterationLimit, + case algorithms.TypeNaiveTableau: + return &tableau_algorithm1.TableauAlgorithm{ + IterationLimit: solver.IterationLimit, }, nil default: - return &algorithms.NaiveAlgorithm{}, fmt.Errorf( + return &tableau_algorithm1.TableauAlgorithm{}, fmt.Errorf( "The Solve() function was given an unknown solver type: %v", algoType, ) } } -func (solver *SimplexSolver) Solve(algoType algorithms.AlgorithmType) (problem.Solution, error) { +func (solver *SimplexSolver) Solve(prob problem.OptimizationProblem) (problem.Solution, error) { // Setup // Choose Algorithm - algo, err := solver.CreateAlgorithm(algoType) + algo, err := solver.CreateAlgorithm(solver.Algorithm) if err != nil { return problem.Solution{}, fmt.Errorf( "The Solve() function was given an unknown solver type: %v", - algoType, + solver.Algorithm, ) } // Apply algorithm - return algo.Solve(solver.InternalState) + return algo.Solve(prob) } diff --git a/testing/algorithms/internal_state_test.go b/testing/algorithms/internal_state_test.go deleted file mode 100644 index 1b53fa4..0000000 --- a/testing/algorithms/internal_state_test.go +++ /dev/null @@ -1,45 +0,0 @@ -package algorithms_test - -import ( - "testing" - - "github.com/MatProGo-dev/MatProInterface.go/problem" - "github.com/MatProGo-dev/SymbolicMath.go/symbolic" - "github.com/MatProGo-dev/simplex/algorithms" - "gonum.org/v1/gonum/mat" -) - -/* -TestToTableau1 -Description: - - In this test, we verify that the ToTableau() function correctly panics when the input problem is not - a linear program. -*/ -func TestToTableau1(t *testing.T) { - // Create a non-linear optimization problem - nonLinearProblem := problem.NewProblem("TestToTableau1 Problem") - - // Create non-linear objective - x := symbolic.NewVariableVector(2) - nonLinearProblem.SetObjective( - x.Transpose().Multiply(x), - problem.SenseMinimize, - ) - - // Create an internal state object - internalState := &algorithms.AlgorithmInternalState{ - BasicVariables: []symbolic.Variable{symbolic.NewVariable(), symbolic.NewVariable()}, - NonBasicVariables: []symbolic.Variable{symbolic.NewVariable()}, - NonBasicValues: mat.NewVecDense(2, nil), - IterationCount: 0, - } - - // Expect a panic when calling ToTableau with a non-linear problem - _, err := internalState.ToTableau(nonLinearProblem) - if err == nil { - // If no error occurs, fail the test - t.Errorf("Expected an error, but got none") - } - -} diff --git a/testing/algorithms/naive_test.go b/testing/algorithms/naive_test.go deleted file mode 100644 index 60b700e..0000000 --- a/testing/algorithms/naive_test.go +++ /dev/null @@ -1,177 +0,0 @@ -package algorithms_test - -import ( - "testing" - - "github.com/MatProGo-dev/simplex/algorithms" - "github.com/MatProGo-dev/simplex/simplexSolver" - "gonum.org/v1/gonum/mat" -) - -/* -Test_NaiveAlgorithm_ComputeFeasibleSolution1 -Description: - - In this test, we verify that the ComputeFeasibleSolution() function correctly computes a feasible - solution for the basic variables. - - We use an example problem from this youtube video: - https://www.youtube.com/watch?v=QAR8zthQypc&t=483s - which is written in standard form as: - maximize 4 x1 + 3 x2 + 5 x3 - subject to - x1 + 2 x2 + 2 x3 + slack1 = 4 - 3 x1 + 4 x3 + slack2 = 12 - 2 x1 + x2 + 4 x3 + slack3 = 8 - - x1, x2, x3 >= 0 - slack1, slack2, slack3 >= 0 - - We expect the feasible solution to be: - [4, 12, 8] - for the initial tableau. -*/ -func Test_NaiveAlgorithm_ComputeFeasibleSolution1(t *testing.T) { - // Setup - problemIn := simplexSolver.GetTestProblem3() - - // Create the solver - solver, err := simplexSolver.For(problemIn, simplexSolver.Configuration{IterationLimit: 100}) - if err != nil { - t.Errorf("Expected no error, but got: %v", err) - } - - // Extract the simplex algorithm and internal state from the solver - algo, err := solver.CreateAlgorithm(algorithms.TypeNaive) - if err != nil { - t.Errorf("Expected no error when creating algorithm; received %v", err) - } - - naiveAlgo, ok := algo.(*algorithms.NaiveAlgorithm) - if !ok { - t.Errorf("Failed to convert the algorithm interface into a Naive Algorithm object (which should be possible given our inputs!)") - } - - initialState := solver.InternalState - - // Compute the feasible solution - solution, err := naiveAlgo.ComputeFeasibleBasicSolution(initialState) - if err != nil { - t.Errorf("Expected no error, but got: %v", err) - } - - // Check that the solution is correct - expectedSolution := mat.NewVecDense(3, []float64{4, 6, 8}) - if !mat.EqualApprox(solution, expectedSolution, 1e-10) { - t.Errorf("Expected feasible solution to be %v, but got %v", expectedSolution, solution) - } -} - -/* -Test_NaiveAlgorithm_SolveLoop1 -Description: - - In this test, we verify that the operations - of the solver loop are correct. - We will use a problem from this youtube video: - https://youtu.be/XMLysZSPsug?si=KMoouByHAV3TTK7h&t=377 - Our method should find that the first Basic Feasible Solution - is {3 2 4 2} and the objective value should be zero. -*/ -func Test_NaiveAlgorithm_SolveLoop1(t *testing.T) { - // Setup - problemIn := simplexSolver.GetTestProblem4() - - // Create the solver - solver, err := simplexSolver.For(problemIn, simplexSolver.Configuration{IterationLimit: 100}) - if err != nil { - t.Errorf("Expected no error, but got: %v", err) - } - - // Extract the simplex algorithm and internal state from the solver - algo, err := solver.CreateAlgorithm(algorithms.TypeNaive) - if err != nil { - t.Errorf("Expected no error when creating algorithm; received %v", err) - } - - naiveAlgo, ok := algo.(*algorithms.NaiveAlgorithm) - if !ok { - t.Errorf("Failed to convert the algorithm interface into a Naive Algorithm object (which should be possible given our inputs!)") - } - - initialState := solver.InternalState - - // Find the first basic feasible solution - basic1, err := naiveAlgo.ComputeFeasibleBasicSolution(initialState) - if err != nil { - t.Errorf("Expected no error, but got: %v", err) - } - - // Check that the solution is correct - expectedFirstBasicSolution := mat.NewVecDense(4, []float64{3, 2, 4, 2}) - if !mat.EqualApprox(basic1, expectedFirstBasicSolution, 1e-10) { - t.Errorf("Expected feasible solution to be %v, but got %v", expectedFirstBasicSolution, basic1) - } - - // Compute the value of the objective function that corresponds to the basic solution - obj1, err := naiveAlgo.ComputeObjectiveFunctionValueWithFeasibleBasicSolution(initialState, basic1) - if err != nil { - t.Errorf("Expected no error, but got: %v", err) - } - - // Check that the objective value is correct (should be zero) - expectedObj1 := 0.0 - if obj1 != expectedObj1 { - t.Errorf("Expected objective value to be %v, but got %v", expectedObj1, obj1) - } - -} - -/* -Test_NaiveAlgorithm_Solve1 -Description: - - In this test, we verify that the operations - of the NaiveAlgorithm.Solve() are correct. - We will use a problem from this youtube video: - https://youtu.be/XMLysZSPsug?si=KMoouByHAV3TTK7h&t=377 - Our method should find that the first Basic Feasible Solution - is {3 2 4 2} and the objective value should be zero. -*/ -func Test_NaiveAlgorithm_Solve1(t *testing.T) { - // Setup - problemIn := simplexSolver.GetTestProblem4() - - // Create the solver - solver, err := simplexSolver.For(problemIn, simplexSolver.Configuration{IterationLimit: 10}) - if err != nil { - t.Errorf("Expected no error, but got: %v", err) - } - - // Extract the simplex algorithm and internal state from the solver - algo, err := solver.CreateAlgorithm(algorithms.TypeNaive) - if err != nil { - t.Errorf("Expected no error when creating algorithm; received %v", err) - } - - naiveAlgo, ok := algo.(*algorithms.NaiveAlgorithm) - if !ok { - t.Errorf("Failed to convert the algorithm interface into a Naive Algorithm object (which should be possible given our inputs!)") - } - - initialState := solver.InternalState - - // Find the first basic feasible solution - _, err = naiveAlgo.Solve(initialState) - // sol1, err := naiveAlgo.Solve(initialState) - if err != nil { - t.Errorf("Unexpected error during solve loops: %v", err) - } - - // t.Errorf( - // "Objective value after %v loops: %v", - // naiveAlgo.IterationLimit, - // sol1.Objective, - // ) - -} diff --git a/testing/algorithms/stanford/state_test.go b/testing/algorithms/stanford/state_test.go index 179c41f..721995f 100644 --- a/testing/algorithms/stanford/state_test.go +++ b/testing/algorithms/stanford/state_test.go @@ -6,14 +6,14 @@ import ( "github.com/MatProGo-dev/SymbolicMath.go/symbolic" stanford_algorithm1 "github.com/MatProGo-dev/simplex/algorithms/stanford" - "github.com/MatProGo-dev/simplex/simplexSolver" "github.com/MatProGo-dev/simplex/utils" + "github.com/MatProGo-dev/simplex/utils/examples" "gonum.org/v1/gonum/mat" ) func TestStanfordAlgorithmState_NonBasicVariables1(t *testing.T) { // Setup - exampleProblem1 := simplexSolver.GetTestProblem1() + exampleProblem1 := examples.GetTestProblem1() // Create the problem in standard form problemInStandardForm, slackVariables, err := exampleProblem1.ToLPStandardForm1() @@ -77,7 +77,7 @@ Description: */ func TestStanfordAlgorithmState_BasicVariables1(t *testing.T) { // Setup - exampleProblem1 := simplexSolver.GetTestProblem1() + exampleProblem1 := examples.GetTestProblem1() // Create the problem in standard form problemInStandardForm, slackVariables, err := exampleProblem1.ToLPStandardForm1() @@ -144,7 +144,7 @@ Description: */ func TestStanfordAlgorithmState_ReducedCostVector1(t *testing.T) { // Setup - exampleProblem1 := simplexSolver.GetTestProblem2() + exampleProblem1 := examples.GetTestProblem2() // Create the problem in standard form problemInStandardForm, slackVariables, err := exampleProblem1.ToLPStandardForm1() diff --git a/testing/algorithms/tableau/selection/blands_rule_test.go b/testing/algorithms/tableau/selection/blands_rule_test.go new file mode 100644 index 0000000..038bdaf --- /dev/null +++ b/testing/algorithms/tableau/selection/blands_rule_test.go @@ -0,0 +1,107 @@ +package tableau_test + +import ( + "testing" + + "github.com/MatProGo-dev/simplex/algorithms/tableau/selection" + "github.com/MatProGo-dev/simplex/utils/examples" +) + +/* +TestBlandsRule_SelectEnteringVariable1 +Description: + + This test will verify that the Bland's Rule selection algorithm correctly selects + the entering variable for the problem shown in minute 21:42 of this youtube video: + https://www.youtube.com/watch?v=-7mCHWpQ9Fw&t=883s +*/ +func TestBlandsRule_SelectEnteringVariable1(t *testing.T) { + // Setup + + // Create the test tableau + testTableau, err := examples.GetTableauExample1() + if err != nil { + t.Errorf("Expected no error, but got: %v", err) + } + + // Create the Bland's Rule selector + selectionRule := selection.BlandsRule{} + + // Find the entering variable + enteringVarIdx := selectionRule.SelectEnteringVariable(*testTableau) + + if enteringVarIdx != 1 { + t.Errorf("Expected entering variable index to be 1, but got %d", enteringVarIdx) + } +} + +/* +TestBlandsRule_SelectExitingVariable1 +Description: + + This test will verify that the Bland's Rule selection algorithm correctly selects + the exiting variable for the problem shown in minute 21:42 of this youtube video: + https://www.youtube.com/watch?v=-7mCHWpQ9Fw&t=883s + The exiting variable selected should be the variable with index 3 (the slack variable + for the second constraint). +*/ +func TestBlandsRule_SelectExitingVariable1(t *testing.T) { + // Setup + + // Create the test tableau + testTableau, err := examples.GetTableauExample1() + if err != nil { + t.Errorf("Expected no error, but got: %v", err) + } + + // Create the Bland's Rule selector + selectionRule := selection.BlandsRule{} + + // Find the entering variable + enteringVarIdx := selectionRule.SelectEnteringVariable(*testTableau) + if enteringVarIdx != 1 { + t.Errorf("Expected entering variable index to be 1, but got %d", enteringVarIdx) + } + + // Find the exiting variable + exitingVarIdx := selectionRule.SelectExitingVariable(*testTableau, enteringVarIdx) + if exitingVarIdx != 3 { + t.Errorf("Expected exiting variable index to be 3, but got %d", exitingVarIdx) + } +} + +/* +TestBlandsRule_SelectEnteringAndExitingVariables1 +Description: + + This test will verify that the Bland's Rule selection algorithm correctly selects + the entering and exiting variables for the problem shown in minute 21:42 of this youtube video: + https://www.youtube.com/watch?v=-7mCHWpQ9Fw&t=883s + The entering variable selected should be the variable with index 1 (x2), + and the exiting variable selected should be the variable with index 3 (the slack variable + for the second constraint). +*/ +func TestBlandsRule_SelectEnteringAndExitingVariables1(t *testing.T) { + // Setup + + // Create the test tableau + testTableau, err := examples.GetTableauExample1() + if err != nil { + t.Errorf("Expected no error, but got: %v", err) + } + + // Create the Bland's Rule selector + selectionRule := selection.BlandsRule{} + + // Find the entering and exiting variables + enteringVarIdx, exitingVarIdx, err := selectionRule.SelectEnteringAndExitingVariables(*testTableau) + if err != nil { + t.Errorf("Expected no error, but got: %v", err) + } + if enteringVarIdx != 1 { + t.Errorf("Expected entering variable index to be 1, but got %d", enteringVarIdx) + } + if exitingVarIdx != 3 { + t.Errorf("Expected exiting variable index to be 3, but got %d", exitingVarIdx) + } +} diff --git a/testing/algorithms/tableau/state_test.go b/testing/algorithms/tableau/state_test.go index bdcf860..de7d526 100644 --- a/testing/algorithms/tableau/state_test.go +++ b/testing/algorithms/tableau/state_test.go @@ -1 +1,70 @@ package tableau + +import ( + "testing" + + tableau_algorithm1 "github.com/MatProGo-dev/simplex/algorithms/tableau" + "github.com/MatProGo-dev/simplex/utils" + "github.com/MatProGo-dev/simplex/utils/examples" +) + +/* +TestTableau_CalculateOptimalValue1 +Description: + + In this test, we verify that the CalculateOptimalValue() function correctly computes the optimal value + of the objective function for a known tableau. + + We use an example problem from this youtube video: + https://www.youtube.com/watch?v=-7mCHWpQ9Fw&t=883s + The optimal solution of the objective function is: + x1 = 125 + x2 = 300 + s1 = 25 + s2 = 0 + s3 = 0 + s4 = 225 +*/ +func TestTableau_CalculateOptimalSolution1(t *testing.T) { + // Setup + + // Create the test problem + testProblem := examples.GetTestProblem5() + + // Use the optimization solver to solve the problem + + // Create initial Tableau state from the problem + initialTableau, err := utils.GetInitialTableauFrom(testProblem) + if err != nil { + t.Errorf("there was an issue creating the initial tableau: %v", err) + } + state0 := tableau_algorithm1.TableauAlgorithmState{ + Tableau: &initialTableau, + IterationCount: 0, + } + + // Update state two times + state1, err := state0.CalculateNextState() + if err != nil { + t.Errorf("there was an issue calculating the next state: %v", err) + } + state2, err := state1.CalculateNextState() + if err != nil { + t.Errorf("there was an issue calculating the next state: %v", err) + } + + // Calculate the optimal solution + solVec, err := state2.CalculateOptimalSolution() + if err != nil { + t.Errorf("there was an issue calculating the optimal value: %v", err) + } + + // Check that the solution is correct + expectedSol := []float64{125.0, 300.0, 25.0, 0.0, 0.0, 225.0} + for ii, val := range expectedSol { + if solVec.AtVec(ii) != val { + t.Errorf("Expected solution value %v at index %d, but got %v", val, ii, solVec.AtVec(ii)) + } + } + +} diff --git a/testing/utils/tableau_test.go b/testing/utils/tableau_test.go index 4b1ec3c..64740c1 100644 --- a/testing/utils/tableau_test.go +++ b/testing/utils/tableau_test.go @@ -6,8 +6,9 @@ import ( "testing" "github.com/MatProGo-dev/SymbolicMath.go/symbolic" - simplexSolver "github.com/MatProGo-dev/simplex/simplexSolver" + "github.com/MatProGo-dev/simplex/algorithms/tableau/selection" "github.com/MatProGo-dev/simplex/utils" + "github.com/MatProGo-dev/simplex/utils/examples" "gonum.org/v1/gonum/mat" ) @@ -20,7 +21,7 @@ Description: */ func TestGetInitialTableau1(t *testing.T) { // Setup - problemIn := simplexSolver.GetTestProblem3() + problemIn := examples.GetTestProblem3() // Create the tableau using the initial state + problem in standard form tableau, err := utils.GetInitialTableauFrom(problemIn) @@ -79,7 +80,7 @@ Description: */ func TestComputeFeasibleSolution1(t *testing.T) { // Setup - problemIn := simplexSolver.GetTestProblem3() + problemIn := examples.GetTestProblem3() // Create the tableau tableau, err := utils.GetInitialTableauFrom(problemIn) @@ -149,3 +150,73 @@ func TestComputeFeasibleSolution1(t *testing.T) { // } // } + +/* +TestTableau_Pivot1 +Description: + + In this test, we verify that the Pivot() function correctly performs a pivot operation on the tableau. + + We use an example problem from this youtube video: + https://www.youtube.com/watch?v=-7mCHWpQ9Fw&t=883s + The initial tableau is: + [ -15 -25 0 0 0 0 0 ] + [ 1 1 1 0 0 0 450 ] + [ 0 1 0 1 0 0 300 ] + [ 4 5 0 0 1 0 2000 ] + [ 1 0 0 0 0 1 350 ] + And after pivoting using Bland's Rule, we expect the tableau to be: + [ -15 0 0 25 0 0 7500 ] + [ 1 0 1 -1 0 0 150 ] + [ 0 1 0 1 0 0 300 ] + [ 4 0 0 -5 1 0 500 ] + [ 1 0 0 0 0 1 350 ] +*/ +func TestTableau_Pivot1(t *testing.T) { + // Setup + + // Create the test tableau + testTableau, err := examples.GetTableauExample1() + if err != nil { + t.Errorf("Expected no error, but got: %v", err) + } + + // Create the Bland's Rule selector + selectionRule := selection.BlandsRule{} + + // Find the entering and exiting variables + enteringVarIdx, exitingVarIdx, err := selectionRule.SelectEnteringAndExitingVariables(*testTableau) + if err != nil { + t.Errorf("Expected no error, but got: %v", err) + } + if enteringVarIdx != 1 { + t.Errorf("Expected entering variable index to be 1, but got %d", enteringVarIdx) + } + if exitingVarIdx != 3 { + t.Errorf("Expected exiting variable index to be 3, but got %d", exitingVarIdx) + } + + // Perform the pivot operation + newTab, err := testTableau.Pivot(enteringVarIdx, exitingVarIdx) + if err != nil { + t.Errorf("Expected no error, but got: %v", err) + } + + // Compute the tableau as a dense matrix + denseTableau := newTab.AsCompressedMatrix + + // Define the expected tableau + expectedTableau := mat.NewDense(5, 7, []float64{ + -15, 0, 0, 25, 0, 0, 7500, + 1, 0, 1, -1, 0, 0, 150, + 0, 1, 0, 1, 0, 0, 300, + 4, 0, 0, -5, 1, 0, 500, + 1, 0, 0, 0, 0, 1, 350, + }) + + // Check that the tableau is correct + if !mat.EqualApprox(denseTableau, expectedTableau, 1e-10) { + t.Errorf("Expected tableau to be:\n%v\nbut got:\n%v", mat.Formatted(expectedTableau), mat.Formatted(denseTableau)) + } + +} diff --git a/simplexSolver/example_problems.go b/utils/examples/optimization_problems.go similarity index 82% rename from simplexSolver/example_problems.go rename to utils/examples/optimization_problems.go index 38b8f2c..d4e8b54 100644 --- a/simplexSolver/example_problems.go +++ b/utils/examples/optimization_problems.go @@ -1,4 +1,4 @@ -package simplexSolver +package examples import ( "github.com/MatProGo-dev/MatProInterface.go/problem" @@ -186,3 +186,39 @@ func GetTestProblem4() *problem.OptimizationProblem { return out } + +/* +GetTestProblem5 +Description: + + Returns the LP from the Youtube video: + https://www.youtube.com/watch?v=-7mCHWpQ9Fw&t=883s +*/ +func GetTestProblem5() *problem.OptimizationProblem { + // Setup + out := problem.NewProblem("TestProblem5") + + // Create variables + x := out.AddVariableVector(2) + + // Create Basic Objective + c := getKVector.From([]float64{15.0, 25.0}) + out.SetObjective( + c.Transpose().Multiply(x), + problem.SenseMaximize, + ) + + // Create Constraints + A := getKMatrix.From([][]float64{ + {1.0, 1.0}, + {0.0, 1.0}, + {4.0, 5.0}, + {1.0, 0.0}, + }) + b := getKVector.From([]float64{450.0, 300.0, 2000.0, 350.0}) + out.Constraints = append(out.Constraints, A.Multiply(x).LessEq(b)) + out.Constraints = append(out.Constraints, x.AtVec(0).GreaterEq(0.0)) + out.Constraints = append(out.Constraints, x.AtVec(1).GreaterEq(0.0)) + + return out +} diff --git a/utils/examples/tableau_examples.go b/utils/examples/tableau_examples.go new file mode 100644 index 0000000..e06f83e --- /dev/null +++ b/utils/examples/tableau_examples.go @@ -0,0 +1,41 @@ +package examples + +import ( + getKMatrix "github.com/MatProGo-dev/SymbolicMath.go/get/KMatrix" + "github.com/MatProGo-dev/SymbolicMath.go/symbolic" + "github.com/MatProGo-dev/simplex/utils" +) + +/* +GetTableauExample1 +Description: + + Returns a simple tableau from the Youtube video: + https://www.youtube.com/watch?v=-7mCHWpQ9Fw&t=883s +*/ +func GetTableauExample1() (*utils.Tableau, error) { + // Setup + + // Create variables + x := symbolic.NewVariableVector(2) + s := symbolic.NewVariableVector(4) + + // Create the condensed tableau matrix directly + condensedKMatrix := getKMatrix.From([][]float64{ + {-15, -25, 0.0, 0.0, 0.0, 0.0, 0.0}, + {1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 450.0}, + {0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 300.0}, + {4.0, 5.0, 0.0, 0.0, 1.0, 0.0, 2000.0}, + {1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 350.0}, + }) + condensed := condensedKMatrix.ToDense() + + // Create the tableau + out := utils.Tableau{ + Variables: append(x, s...), + BasicVariableIndicies: []int{2, 3, 4, 5}, + AsCompressedMatrix: &condensed, + } + + return &out, out.Check() +} diff --git a/utils/tableau.go b/utils/tableau.go index 7a948ce..6562ae6 100644 --- a/utils/tableau.go +++ b/utils/tableau.go @@ -2,6 +2,7 @@ package utils import ( "fmt" + "math" "github.com/MatProGo-dev/MatProInterface.go/problem" getKMatrix "github.com/MatProGo-dev/SymbolicMath.go/get/KMatrix" @@ -11,10 +12,23 @@ import ( ) // The Tableau Representation of a linear program in standard form +// Specifically, the tableau represents the problem: +// +// minimize c^T * x + d +// subject to A * x = b +// x >= 0 +// +// It represents the problem by storing the following information: +// - The list of all variables in the problem (including slack variables) +// - The indicies of the basic variables in the list of all variables +// - A Matrix representing the tableau as follows +// +// | c^T | d | +// | A | b | type Tableau struct { Variables []symbolic.Variable BasicVariableIndicies []int // The basic variables in order of their connection to the constraint rows - AsCompressedMatrix *mat.Dense // The compressed matrix contains all + AsCompressedMatrix *mat.Dense // The compressed matrix contains all of the information } /* @@ -158,7 +172,11 @@ func (tableau *Tableau) C() *mat.VecDense { topRow := compressedMatrixCopy.RowView(0) topRowAsVecDense, _ := topRow.(*mat.VecDense) - return mat.NewVecDense(nTableauCols-1, topRowAsVecDense.RawVector().Data) + // Return the C vector (excluding the last entry which is the constant term) + out := mat.NewVecDense(nTableauCols-1, nil) + out.CopyVec(topRowAsVecDense.SliceVec(0, nTableauCols-1)) + + return out } /* @@ -254,6 +272,18 @@ func (tableau *Tableau) CNonBasic() (*mat.VecDense, error) { return &cNonBasicAsVecDense, nil } +func (tableau *Tableau) D() float64 { + // Check that tableau is valid + err := tableau.Check() + if err != nil { + panic(err) + } + + // Return the value of d + _, nTableauCols := tableau.AsCompressedMatrix.Dims() + return tableau.AsCompressedMatrix.At(0, nTableauCols-1) +} + func (tableau *Tableau) NonBasicVariableIndicies() []int { // Input Processing err := tableau.Check() @@ -290,6 +320,11 @@ func (tableau *Tableau) NonBasicVariables() []symbolic.Variable { return out } +func (tableau *Tableau) NumberOfConstraints() int { + nRows, _ := tableau.AsCompressedMatrix.Dims() + return nRows - 1 +} + func (tableau *Tableau) BasicVariables() []symbolic.Variable { // Input Processing err := tableau.Check() @@ -330,7 +365,7 @@ func GetInitialTableauFrom(problemIn *problem.OptimizationProblem) (Tableau, err // Transform the problem into the standard form where all constraints // are equality constraints - problemInStandardForm, slackVariables, err := problemIn.ToLPStandardForm1() + problemInStandardForm, slackVariables, err := problemIn.ToLPStandardForm2() if err != nil { return Tableau{}, err } @@ -343,18 +378,24 @@ func GetInitialTableauFrom(problemIn *problem.OptimizationProblem) (Tableau, err } // Create the matrix + // [ A | b ] A, b, err := problemInStandardForm.LinearEqualityConstraintMatrices() if err != nil { return Tableau{}, err } Ab := symbolic.HStack(A, b) + // Create the cost vector + // [ -c^T | -d ] objectiveExpression := problemInStandardForm.Objective.Expression.(symbolic.ScalarExpression) c := objectiveExpression.LinearCoeff(problemInStandardForm.Variables) + d := objectiveExpression.Constant() + + c.ScaleVec(-1.0, &c) // We negate c because we are converting from a maximization to a minimization problem var cExtended *mat.VecDense = mat.NewVecDense(c.Len()+1, nil) cExtended.CopyVec(&c) - cExtended.SetVec(c.Len(), 0.0) + cExtended.SetVec(c.Len(), -d) tableauMatCondensed := symbolic.VStack( symbolic.VecDenseToKVector(*cExtended).Transpose(), @@ -393,6 +434,28 @@ func (tableau *Tableau) AllObjectiveRowEntriesAreLessThanOrEqualToZero() bool { return true } +/* +CanNotBeImproved +Description: + + Returns true if the tableau's objective row indicates that + the objective function can not be improved by changing + any of the non-basic variables. +*/ +func (tableau *Tableau) CanNotBeImproved() bool { + // Get the coefficients of the non-basic variables + c := tableau.C() + + // Check if all coefficients are less than or equal to zero + for ii := 0; ii < c.Len(); ii++ { + if c.AtVec(ii) < 0 { + return false + } + } + + return true +} + /* ComputeFeasibleSolution Description: @@ -499,44 +562,92 @@ func (tableau *Tableau) NumberOfNonBasicVariables() int { } /* -SelectPivotColumn +Pivot Description: - Selects the pivot column based on the tableau. - This is the column from the non-basic variables that leads to the - largest increase in the objective function value. - It returns the index of the pivot column in the tableau. + Performs a pivot operation on the tableau. + This operation produces a new tableau where the entering variable + becomes a basic variable and the exiting variable becomes a non-basic variable. */ -func (tableau *Tableau) SelectPivotColumn() (int, error) { +func (tableau *Tableau) Pivot(enteringVarIdx int, exitingVarIdx int) (Tableau, error) { // Input Processing err := tableau.Check() if err != nil { - return -1, fmt.Errorf("SelectPivotColumn: %v", err) + return Tableau{}, fmt.Errorf("Pivot: %v", err) + } + + // Check that the entering variable is not already a basic variable + if foundIdx, _ := symbolic.FindInSlice(enteringVarIdx, tableau.BasicVariableIndicies); foundIdx != -1 { + return Tableau{}, fmt.Errorf("Pivot: the entering variable is already a basic variable") + } + + // Check that the exiting variable is a basic variable + if foundIdx, _ := symbolic.FindInSlice(exitingVarIdx, tableau.BasicVariableIndicies); foundIdx == -1 { + return Tableau{}, fmt.Errorf("Pivot: the exiting variable is not a basic variable") } // Setup - nNonBasic := tableau.NumberOfNonBasicVariables() - if nNonBasic == 0 { - return -1, fmt.Errorf("SelectPivotColumn: there are no non-basic variables") + nRows, nCols := tableau.AsCompressedMatrix.Dims() + newTableauMat := mat.NewDense(nRows, nCols, nil) + newTableauMat.Copy(tableau.AsCompressedMatrix) + + exitingConstraintIdx, err := symbolic.FindInSlice(exitingVarIdx, tableau.BasicVariableIndicies) + if err != nil { + return Tableau{}, fmt.Errorf("Pivot: could not find exiting variable in basic variable indicies (%v)", err) } - // Get the coefficients of the non-basic variables - objectiveExpression := getKVector.From(tableau.C()).Transpose().Multiply( - symbolic.VariableVector(tableau.Variables), - ) - objectiveSE, _ := objectiveExpression.(symbolic.ScalarExpression) + // Perform the pivot operation + // - "Normalize" the pivot element (i.e., the element at A[exitingVarIdx, enteringVarIdx]) + normalizingFactor := 1.0 / tableau.AsCompressedMatrix.At(exitingConstraintIdx+1, enteringVarIdx) + newPivotRow := mat.NewVecDense(nCols, nil) + for jj := 0; jj < nCols; jj++ { + newPivotRow.SetVec(jj, tableau.AsCompressedMatrix.At(exitingConstraintIdx+1, jj)*normalizingFactor) + } + newTableauMat.SetRow(exitingConstraintIdx+1, newPivotRow.RawVector().Data) + + // - Zero out the other entries in the entering variable column + for ii := 0; ii < nRows; ii++ { + // Skip the pivot row (i.e., the row of the entering variable) + if ii == exitingConstraintIdx+1 { + continue + } + // Skip any row that is already zero + if math.Abs(newTableauMat.At(ii, enteringVarIdx)) < 1e-14 { + continue + } + + // Otherwise, determine the factor needed to zero out this entry + factorToZeroOut := -1.0 * newTableauMat.At(ii, enteringVarIdx) + rowToAdd := mat.NewVecDense(nCols, nil) + rowToAdd.ScaleVec(factorToZeroOut, newPivotRow) + rowToAdd.AddVec(rowToAdd, newTableauMat.RowView(ii)) - c := objectiveSE.LinearCoeff(tableau.NonBasicVariables()) + // Update the row in the new tableau + newTableauMat.SetRow(ii, rowToAdd.RawVector().Data) + } - // Find the index of the pivot column - pivotColIndex := -1 - maxValue := 0.0 - for ii := 0; ii < nNonBasic; ii++ { - if c.AtVec(ii) > maxValue { - maxValue = c.AtVec(ii) - pivotColIndex = ii + // Update the list of basic variable indicies + newBasicVariableIndicies := make([]int, len(tableau.BasicVariableIndicies)) + copy(newBasicVariableIndicies, tableau.BasicVariableIndicies) + for ii, bvIdx := range tableau.BasicVariableIndicies { + if bvIdx == exitingVarIdx { + newBasicVariableIndicies[ii] = enteringVarIdx + break } } - return pivotColIndex, nil + // Create the new tableau + newTableau := Tableau{ + Variables: tableau.Variables, + BasicVariableIndicies: newBasicVariableIndicies, + AsCompressedMatrix: newTableauMat, + } + + // Check the new tableau for validity + err = newTableau.Check() + if err != nil { + return Tableau{}, fmt.Errorf("Pivot: the resulting tableau is invalid (%v)", err) + } + + return newTableau, nil }