Created OptimizationProblem's new LinearInequalityConstraintMatrices and LinearEqualityConstraintMatrices (#8)

* Upgraded SymbolicMath

* Bumped the version of SymbolicMath.go

* Introduced method for assembling linear inequality constraints into single matrices

* Updated main readme to take advantage of new problem structure

* wip: First attempt at creating the LinearEqualityConstraintMatrices function

* Added tests for the new Equality constraint matrices function

* Added tests for coverage of the case when there are mixtures of scalar and vector constraints

Author: kwesiRutledge

README.md

@@ -26,59 +26,40 @@ by your model.
 ## Modeling the Mathematical Program Above
 For example, to model the program above one would write the following code:
 ```
-// Constants
-modelName := "mpg-qp1"
-m := optim.NewModel(modelName)
-x := m.AddVariableVector(2)
 
-// Create Vector Constants
-c1 := optim.KVector(
-    *mat.NewVecDense(2, []float64{0.0, 1.0}),
+import (
+	...
+    "github.com/MatProGo-dev/MatProInterface.go/problem"
+	getKVector "github.com/MatProGo-dev/SymbolicMath.go/get/KVector"
+    ...
 )
 
-c2 := optim.KVector(
-    *mat.NewVecDense(2, []float64{2.0, 3.0}),
-)
+// Constants
+problemName := "mpg-qp1"
+p1 := problem.NewProblem(problemName)
+x := p1.AddVariableVector(2)
 
-// Use these to create constraints.
+// Create Vector Constants
+c1 := getKVector.From([]float64{0.0, 1.0})
+c2 := getKVector.From([]float64{2.0, 3.0})
 
-vc1, err := x.LessEq(c2)
-if err != nil {
-    t.Errorf("There was an issue creating the proper vector constraint: %v", err)
-}
+// Use these to create constraints.
+vc1 := x.LessEq(c2)
+vc2 := x.GreaterEq(c1)
 
-vc2, err := x.GreaterEq(c1)
-if err != nil {
-    t.Errorf("There was an issue creating the proper vector constraint: %v", err)
-}
+p1.Constraints = append(p1.Constraints, vc1)
+p1.Constraints = append(p1.Constraints, vc2)
 
 // Create objective
-Q1 := optim.Identity(x.Len())
+Q1 := symbolic.Identity(x.Len())
 Q1.Set(0, 1, 0.25)
 Q1.Set(1, 0, 0.25)
 Q1.Set(1, 1, 0.25)
 
-obj := optim.ScalarQuadraticExpression{
-    Q: Q1,
-    X: x,
-    L: *mat.NewVecDense(x.Len(), []float64{0, -0.97}),
-    C: 2.0,
-}
-
-// Add Constraints
-constraints := []optim.Constraint{vc1, vc2}
-for _, constr := range constraints {
-    err = m.AddConstraint(constr)
-    if err != nil {
-        t.Errorf("There was an issue adding the vector constraint to the model: %v", err)
-    }
-}
-
-// Add objective
-err = m.SetObjective(optim.Objective{obj, optim.SenseMinimize})
-if err != nil {
-    t.Errorf("There was an issue setting the objective of the Gurobi solver model: %v", err)
-}
+p1.Objective = *problem.NewObjective(
+    x.Transpose().Multiply(Q).Multiply(x),
+    problem.SenseMinimize,
+)
 
 // Solve using the solver of your choice!
 ```

go.mod

@@ -4,4 +4,4 @@ go 1.21
 
 require gonum.org/v1/gonum v0.14.0
 
-require github.com/MatProGo-dev/SymbolicMath.go v0.1.8
+require github.com/MatProGo-dev/SymbolicMath.go v0.2.1

go.sum

@@ -1,5 +1,11 @@
 github.com/MatProGo-dev/SymbolicMath.go v0.1.8 h1:lpe+6cK/2fg29WwxOykm4hKvfJeqvFUBGturC9qh5ug=
 github.com/MatProGo-dev/SymbolicMath.go v0.1.8/go.mod h1:gKbGR/6sYWi2koMUEDIPWBPi6jQPELKle0ijIM+eaHU=
+github.com/MatProGo-dev/SymbolicMath.go v0.1.9 h1:NMqvS9Bt2DWWLGxd+j3Qta4Ckq/x74gpMM7bt32om5g=
+github.com/MatProGo-dev/SymbolicMath.go v0.1.9/go.mod h1:gKbGR/6sYWi2koMUEDIPWBPi6jQPELKle0ijIM+eaHU=
+github.com/MatProGo-dev/SymbolicMath.go v0.2.0 h1:W8IREsZGeIuPAKHgJDyeCr3vLJjMr6H/O9RNFAjOVp8=
+github.com/MatProGo-dev/SymbolicMath.go v0.2.0/go.mod h1:gKbGR/6sYWi2koMUEDIPWBPi6jQPELKle0ijIM+eaHU=
+github.com/MatProGo-dev/SymbolicMath.go v0.2.1 h1:3qcLYNx3+9Ud/saS4QWxJzmDUbVvYZIWt9aX2bQ9iOk=
+github.com/MatProGo-dev/SymbolicMath.go v0.2.1/go.mod h1:gKbGR/6sYWi2koMUEDIPWBPi6jQPELKle0ijIM+eaHU=
 golang.org/x/exp v0.0.0-20230321023759-10a507213a29 h1:ooxPy7fPvB4kwsA2h+iBNHkAbp/4JxTSwCmvdjEYmug=
 golang.org/x/exp v0.0.0-20230321023759-10a507213a29/go.mod h1:CxIveKay+FTh1D0yPZemJVgC/95VzuuOLq5Qi4xnoYc=
 gonum.org/v1/gonum v0.14.0 h1:2NiG67LD1tEH0D7kM+ps2V+fXmsAnpUeec7n8tcr4S0=

problem/optimization_problem.go

@@ -5,7 +5,9 @@ import (
 
 	"github.com/MatProGo-dev/MatProInterface.go/mpiErrors"
 	"github.com/MatProGo-dev/MatProInterface.go/optim"
+	getKVector "github.com/MatProGo-dev/SymbolicMath.go/get/KVector"
 	"github.com/MatProGo-dev/SymbolicMath.go/symbolic"
+	"gonum.org/v1/gonum/mat"
 )
 
 // OptimizationProblem represents the overall constrained linear optimization model to be
@@ -340,3 +342,196 @@ func (op *OptimizationProblem) IsLinear() bool {
 	// All Checks Passed!
 	return true
 }
+
+/*
+LinearInequalityConstraintMatrices
+Description:
+
+	Returns the linear INEQUALITY constraint matrices and vectors.
+	For all linear inequality constraints, we assemble them into the form:
+		Ax <= b
+	Where A is the matrix of coefficients, x is the vector of variables, and b is the vector of constants.
+	We return A and b.
+*/
+func (op *OptimizationProblem) LinearInequalityConstraintMatrices() (symbolic.KMatrix, symbolic.KVector) {
+	// Setup
+
+	// Collect the Variables of this Problem
+	x := op.Variables
+
+	// Iterate through all constraints and collect the linear constraints
+	// into a matrix and vector.
+	scalar_constraints := make([]symbolic.ScalarConstraint, 0)
+	vector_constraints := make([]symbolic.VectorConstraint, 0)
+	for _, constraint := range op.Constraints {
+		// Skip this constraint if it is not linear
+		if !constraint.IsLinear() {
+			continue
+		}
+		// Skip this constraint if it is not an inequality
+		if constraint.ConstrSense() == symbolic.SenseEqual {
+			continue
+		}
+		switch c := constraint.(type) {
+		case symbolic.ScalarConstraint:
+			scalar_constraints = append(scalar_constraints, c)
+		case symbolic.VectorConstraint:
+			vector_constraints = append(vector_constraints, c)
+		}
+	}
+
+	// Create the matrix and vector elements from the scalar constraints
+	A_components_scalar := make([]mat.VecDense, len(scalar_constraints))
+	b_components_scalar := make([]float64, len(scalar_constraints))
+	for ii, constraint := range scalar_constraints {
+		A_components_scalar[ii], b_components_scalar[ii] = constraint.LinearInequalityConstraintRepresentation(x)
+	}
+
+	// Create the matrix and vector elements from the vector constraints
+	A_components_vector := make([]mat.Dense, len(vector_constraints))
+	b_components_vector := make([]mat.VecDense, len(vector_constraints))
+	for ii, constraint := range vector_constraints {
+		A_components_vector[ii], b_components_vector[ii] = constraint.LinearInequalityConstraintRepresentation(x)
+	}
+
+	// Assemble the matrix and vector components
+	var AOut symbolic.Expression
+	var bOut symbolic.Expression
+	scalar_constraint_matrices_exist := len(A_components_scalar) > 0
+	if scalar_constraint_matrices_exist {
+		AOut = symbolic.VecDenseToKVector(A_components_scalar[0]).Transpose()
+		for ii := 1; ii < len(A_components_scalar); ii++ {
+			AOut = symbolic.VStack(
+				AOut,
+				symbolic.VecDenseToKVector(A_components_scalar[ii]).Transpose(),
+			)
+		}
+		bOut = getKVector.From(b_components_scalar)
+	}
+
+	vector_constraint_matrices_exist := len(A_components_vector) > 0
+	if vector_constraint_matrices_exist {
+		// Create the matrix, if it doesn't already exist
+		if !scalar_constraint_matrices_exist {
+			AOut = symbolic.DenseToKMatrix(A_components_vector[0])
+			bOut = symbolic.VecDenseToKVector(b_components_vector[0])
+		} else {
+			AOut = symbolic.VStack(
+				AOut,
+				symbolic.DenseToKMatrix(A_components_vector[0]),
+			)
+			bOut = symbolic.VStack(
+				bOut,
+				symbolic.VecDenseToKVector(b_components_vector[0]),
+			)
+		}
+		for ii := 1; ii < len(A_components_vector); ii++ {
+			AOut = symbolic.VStack(
+				AOut,
+				symbolic.DenseToKMatrix(A_components_vector[ii]),
+			)
+			bOut = symbolic.VStack(
+				bOut,
+				symbolic.VecDenseToKVector(b_components_vector[ii]),
+			)
+		}
+	}
+
+	return AOut.(symbolic.KMatrix), bOut.(symbolic.KVector)
+}
+
+/*
+LinearEqualityConstraintMatrices
+Description:
+
+	Returns the linear EQUALITY constraint matrices and vectors.
+	For all linear equality constraints, we assemble them into the form:
+		Cx = d
+	Where C is the matrix of coefficients, x is the vector of variables, and d is the vector of constants.
+	We return C and d.
+*/
+func (op *OptimizationProblem) LinearEqualityConstraintMatrices() (symbolic.KMatrix, symbolic.KVector) {
+	// Setup
+
+	// Collect the Variables of this Problem
+	x := op.Variables
+
+	// Iterate through all constraints and collect the linear constraints
+	// into a matrix and vector.
+	scalar_constraints := make([]symbolic.ScalarConstraint, 0)
+	vector_constraints := make([]symbolic.VectorConstraint, 0)
+	for _, constraint := range op.Constraints {
+		// Skip this constraint if it is not linear
+		if !constraint.IsLinear() {
+			continue
+		}
+		// Skip this constraint if it is not an equality
+		if constraint.ConstrSense() != symbolic.SenseEqual {
+			continue
+		}
+		switch c := constraint.(type) {
+		case symbolic.ScalarConstraint:
+			scalar_constraints = append(scalar_constraints, c)
+		case symbolic.VectorConstraint:
+			vector_constraints = append(vector_constraints, c)
+		}
+	}
+
+	// Create the matrix and vector elements from the scalar constraints
+	C_components_scalar := make([]mat.VecDense, len(scalar_constraints))
+	d_components_scalar := make([]float64, len(scalar_constraints))
+	for ii, constraint := range scalar_constraints {
+		C_components_scalar[ii], d_components_scalar[ii] = constraint.LinearEqualityConstraintRepresentation(x)
+	}
+
+	// Create the matrix and vector elements from the vector constraints
+	C_components_vector := make([]mat.Dense, len(vector_constraints))
+	d_components_vector := make([]mat.VecDense, len(vector_constraints))
+	for ii, constraint := range vector_constraints {
+		C_components_vector[ii], d_components_vector[ii] = constraint.LinearEqualityConstraintRepresentation(x)
+	}
+
+	// Assemble the matrix and vector components
+	var COut symbolic.Expression
+	var dOut symbolic.Expression
+	scalar_constraint_matrices_exist := len(C_components_scalar) > 0
+	if scalar_constraint_matrices_exist {
+		COut = symbolic.VecDenseToKVector(C_components_scalar[0]).Transpose()
+		for ii := 1; ii < len(C_components_scalar); ii++ {
+			COut = symbolic.VStack(
+				COut,
+				symbolic.VecDenseToKVector(C_components_scalar[ii]).Transpose(),
+			)
+		}
+		dOut = getKVector.From(d_components_scalar)
+	}
+	vector_constraint_matrices_exist := len(C_components_vector) > 0
+
+	if vector_constraint_matrices_exist {
+		// Create the matrix, if it doesn't already exist
+		if !scalar_constraint_matrices_exist {
+			COut = symbolic.DenseToKMatrix(C_components_vector[0])
+			dOut = symbolic.VecDenseToKVector(d_components_vector[0])
+		} else {
+			COut = symbolic.VStack(
+				COut,
+				symbolic.DenseToKMatrix(C_components_vector[0]),
+			)
+			dOut = symbolic.VStack(
+				dOut,
+				symbolic.VecDenseToKVector(d_components_vector[0]),
+			)
+		}
+		for ii := 1; ii < len(C_components_vector); ii++ {
+			COut = symbolic.VStack(
+				COut,
+				symbolic.DenseToKMatrix(C_components_vector[ii]),
+			)
+			dOut = symbolic.VStack(
+				dOut,
+				symbolic.VecDenseToKVector(d_components_vector[ii]),
+			)
+		}
+	}
+	return COut.(symbolic.KMatrix), dOut.(symbolic.KVector)
+}