diff --git a/symbolic/constant_vector.go b/symbolic/constant_vector.go index ef9d1da..bcbb093 100644 --- a/symbolic/constant_vector.go +++ b/symbolic/constant_vector.go @@ -145,7 +145,8 @@ func (kv KVector) Plus(rightIn interface{}) Expression { // Constants kvLen := kv.Len() - // Management + // Algorithm + var out Expression switch right := rightIn.(type) { case float64: // Check to see if the output is a vector or a scalar @@ -159,42 +160,26 @@ func (kv KVector) Plus(rightIn interface{}) Expression { eAsVec.ScaleVec(right, &tempOnes) // Add the values - return kv.Plus(VecDenseToKVector(eAsVec)) - case K, Variable, Monomial, Polynomial: - // Create a new polynomial vector - var out []ScalarExpression - for _, element := range kv { - out = append(out, element.Plus(right).(ScalarExpression)) - } - return ConcretizeVectorExpression(out) + out = kv.Plus(VecDenseToKVector(eAsVec)) + case Expression: + out = VectorPlusTemplate(kv, right) case *mat.VecDense: - return kv.Plus(VecDenseToKVector(*right)) // Convert to KVector + out = kv.Plus(VecDenseToKVector(*right)) // Convert to KVector case mat.VecDense: - return kv.Plus(VecDenseToKVector(right)) // Convert to KVector - - case KVector: - // Compute Addition - var result mat.VecDense - kvAsVec := kv.ToVecDense() - eAsVec := right.ToVecDense() - result.AddVec(&kvAsVec, &eAsVec) - - return VecDenseToKVector(result) - - case VariableVector: - return right.Plus(kv) - - case MonomialVector: - return right.Plus(kv) - - case PolynomialVector: - return right.Plus(kv) + out = kv.Plus(VecDenseToKVector(right)) // Convert to KVector default: - errString := fmt.Sprintf("Unrecognized expression type %T for addition of KVector kv.Plus(%v)!", right, right) - panic(fmt.Errorf(errString)) + panic( + smErrors.UnsupportedInputError{ + FunctionName: "KVector.Plus", + Input: rightIn, + }, + ) } + + // Simplify and return + return out.AsSimplifiedExpression() } /* diff --git a/symbolic/monomial_vector.go b/symbolic/monomial_vector.go index 1563f2a..65b44a1 100644 --- a/symbolic/monomial_vector.go +++ b/symbolic/monomial_vector.go @@ -150,92 +150,10 @@ func (mv MonomialVector) Plus(term1 interface{}) Expression { switch right := term1.(type) { case float64: out = mv.Plus(K(right)) - case K: - // Convert the scalar to a scalar vector - tempVD := OnesVector(mv.Len()) - tempVD.ScaleVec(float64(right), &tempVD) - - out = mv.Plus(VecDenseToKVector(tempVD)) - case Monomial: - // Check to see if all elements of the monomial vector, - // are all monomials like the input monomial. - monomialVectorMatches := true - for _, monomial := range mv { - if !monomial.MatchesFormOf(right) { - monomialVectorMatches = false - } - } - - if monomialVectorMatches { - // If all elements of the monomial vector are monomials like the input monomial, - // then simply add the coefficients. and return a monomial vector. - var mvOut MonomialVector - for _, monomial := range mv { - mvOut = append(mvOut, monomial.Plus(right).(Monomial)) - } - out = mvOut - } else { - // Otherwise, create a polynomial vector - var pv PolynomialVector - for _, monomial := range mv { - pv = append(pv, monomial.Plus(right).(Polynomial)) - } - out = pv.Simplify() - } - case KVector: - if mv.IsConstant() { - // If monomial vector is really a constant vector, - // then don't convert down but simply update the coefficients. - var kvOut KVector = VecDenseToKVector(mv.Constant()) - out = kvOut.Plus(right) - } else { - // Create a polynomial vector - var ve []ScalarExpression - for ii, monomial := range mv { - ve = append(ve, monomial.Plus(right[ii]).(ScalarExpression)) - } - out = ConcretizeVectorExpression(ve) - } - case MonomialVector: - // Check to see if all elements of the monomial vector, - // are all monomials like the input monomial. - monomialVectorMatches := true - for ii, monomial := range mv { - if !monomial.MatchesFormOf(right[ii]) { - monomialVectorMatches = false - } - } - - if monomialVectorMatches { - // If all elements of the monomial vector are monomials like the input monomial, - // then simply add the coefficients. and return a monomial vector. - var mvOut MonomialVector - for ii, monomial := range mv { - mvOut = append(mvOut, monomial.Plus(right[ii]).(Monomial)) - } - out = mvOut - } else { - // Otherwise, create a polynomial vector - var pv PolynomialVector - for ii, monomial := range mv { - sumII := monomial.Plus(right[ii]) - switch sumII.(type) { - case Monomial: - pv = append(pv, sumII.(Monomial).ToPolynomial()) - case Polynomial: - pv = append(pv, sumII.(Polynomial)) - default: - panic( - fmt.Errorf( - "Unexpected type of sumII in the Plus() method: %T (%v)", - sumII, sumII, - ), - ) - - } - } - out = pv.Simplify() - } + case int: + out = mv.Plus(K(float64(right))) + case Expression: + out = VectorPlusTemplate(mv, right) default: // If the right hand side is an unsupported type, then panic panic( diff --git a/symbolic/polynomial_vector.go b/symbolic/polynomial_vector.go index df2cb85..45dc581 100644 --- a/symbolic/polynomial_vector.go +++ b/symbolic/polynomial_vector.go @@ -196,46 +196,12 @@ func (pv PolynomialVector) Plus(e interface{}) Expression { switch right := e.(type) { case float64: out = pv.Plus(K(right)) - case K: - pvCopy := pv - - // Algorithm - var sum []ScalarExpression - for _, polynomial := range pvCopy { - tempSum := polynomial.Plus(right) - sum = append(sum, tempSum.(ScalarExpression)) - } - out = ConcretizeExpression(sum) - case Variable: - pvCopy := pv - for ii, polynomial := range pv { - sum := polynomial.Plus(right) - pvCopy[ii] = sum.(Polynomial) - } - out = pvCopy - - case Polynomial: - pvCopy := pv - - // Algorithm - for ii, polynomial := range pv { - sum := polynomial.Plus(right) - pvCopy[ii] = sum.(Polynomial) - } - out = pvCopy - case VectorExpression: - pvCopy := pv - - // Cast right - rightAsVector, _ := ToVectorExpression(right) - - // Algorithm - var sum []ScalarExpression - for ii, polynomial := range pvCopy { - tempSum := polynomial.Plus(rightAsVector.AtVec(ii)) - sum = append(sum, tempSum.(ScalarExpression)) - } - out = ConcretizeExpression(sum) + // The generic Expression case is sufficient because VectorPlusTemplate + // is designed to handle all supported vector types (e.g., PolynomialVector, + // VectorLinearExpr, etc.) and correctly implements addition for them. + // This ensures that all previously supported types are handled safely. + case Expression: + out = VectorPlusTemplate(pv, right) default: // Default response is a panic panic( diff --git a/testing/symbolic/constant_matrix_test.go b/testing/symbolic/constant_matrix_test.go index 0c2c559..c922b55 100644 --- a/testing/symbolic/constant_matrix_test.go +++ b/testing/symbolic/constant_matrix_test.go @@ -14,8 +14,10 @@ import ( "testing" getKMatrix "github.com/MatProGo-dev/SymbolicMath.go/get/KMatrix" + getKVector "github.com/MatProGo-dev/SymbolicMath.go/get/KVector" "github.com/MatProGo-dev/SymbolicMath.go/smErrors" "github.com/MatProGo-dev/SymbolicMath.go/symbolic" + "gonum.org/v1/gonum/mat" ) /* @@ -745,6 +747,41 @@ func TestKMatrix_Multiply10(t *testing.T) { } } +/* +TestKMatrix_Multiply11 +Description: + + Tests that the Multiply() method properly multiplies a large KMatrix + that was causing issues when trying to debug the tricky gonum LP: + https://github.com/gonum/gonum/issues/1914 + The KMatrix should have 57 variables. +*/ +func TestKMatrix_Multiply11(t *testing.T) { + // Setup + varCount := 57 + + // Create the variables + x := symbolic.NewVariableVector(varCount) + + // - Linear Equality Constraints + b := getKVector.From( + []float64{38, 5, 2, 33, 28, 14, 2, 48, 8, 133, 117, 34, 48, 20, 16, 50, 30, 75, 10, 40, 6, 70, 10, 5, 167, 13, 2, 118, 12, 98, 67, 157, 55, 2, 4, 4, 25, 4}, + ) + + Abuf := []float64{ + 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.95, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.95} + A := mat.NewDense(38, varCount, Abuf) + + // Create constraint with NO panics + defer func() { + if r := recover(); r != nil { + t.Errorf("Did not expect panic; received %v", r) + } + }() + + getKMatrix.From(A).Multiply(x).Eq(b) +} + /* TestKMatrix_Transpose1 Description: