Reworked Simplify to work better on Polynomials (with ripple effects)

examples/derivative1/derivative1.go

@@ -0,0 +1,41 @@
+package main
+
+import (
+	"fmt"
+
+	"github.com/MatProGo-dev/SymbolicMath.go/symbolic"
+)
+
+/*
+derivative1.go
+Description:
+
+	This script is meant to construct a simple polynomial
+	representing an object in free fall within standard
+	earth gravity, `g`, with an initial upward velocity of `v0`.
+
+	By differentiating the position equation, you should get the
+	velocity equation.
+*/
+
+func main() {
+	// Setup
+	t := symbolic.NewVariable()
+	t.Name = "t"
+
+	g := 9.81 // m/s^2
+	v0 := 3.0 // m/2^2
+
+	// Create quadratic function
+	yPosition := symbolic.K(-0.5 * g).Multiply(
+		t.Power(2),
+	).Plus(
+		t.Multiply(v0),
+	)
+
+	// Create the derivative
+	yVelocity := yPosition.DerivativeWrt(t)
+
+	// Print the polynomial
+	fmt.Println(yVelocity.String())
+}

examples/substitute1/substitute1.go

@@ -0,0 +1,52 @@
+package main
+
+import (
+	"fmt"
+
+	getKMatrix "github.com/MatProGo-dev/SymbolicMath.go/get/KMatrix"
+	"github.com/MatProGo-dev/SymbolicMath.go/symbolic"
+)
+
+/*
+substitute1.go
+Description:
+
+	This script is meant to show how the substitution
+	method can be used to evaluate the value of your
+	expressions. In this case, we will use it to evaluate
+	a quadratic expression at the point [-3.0, 1.0].
+
+	The result should be a value of 10.0.
+*/
+
+func main() {
+	// Construct quadratic expression
+	N := 2
+	x := symbolic.NewVariableVector(N)
+	Q := getKMatrix.From(
+		symbolic.Identity(N))
+
+	// Create the quadratic polynomial
+	quadPoly := x.Transpose().Multiply(Q).Multiply(x)
+	fmt.Println("Polynomial is:", quadPoly)
+
+	// Create the "map" used for substitution
+	x0 := x[0]
+	x1 := x[1]
+
+	subMap := map[symbolic.Variable]symbolic.Expression{
+		x0: symbolic.K(-3.0),
+		x1: symbolic.K(1.0),
+	}
+
+	quadEvaluated := quadPoly.SubstituteAccordingTo(subMap)
+
+	// Print the polynomial
+	fmt.Println(
+		"Polynomial at the point x = [",
+		subMap[x0],
+		",",
+		subMap[x1],
+		"] is",
+		quadEvaluated.String())
+}

symbolic/constant_matrix.go

@@ -111,6 +111,8 @@ func (km KMatrix) Plus(e interface{}) Expression {
 	dims := km.Dims()
 	nR, nC := dims[0], dims[1]
 
+	// Algorithm
+	var out Expression
 	switch right := e.(type) {
 	case float64:
 		// Create a matrix of all elements with value right
@@ -125,14 +127,13 @@ func (km KMatrix) Plus(e interface{}) Expression {
 		var sumAsDense mat.Dense
 		sumAsDense.Add(&rightAsDense, &kmAsDense)
 
-		return DenseToKMatrix(sumAsDense)
+		out = DenseToKMatrix(sumAsDense)
 
 	case K:
-		return km.Plus(float64(right)) // Reuse float64 case
+		out = km.Plus(float64(right)) // Reuse float64 case
 
 	case Variable:
-		// Create a matrix of variables where each element has
-		// the value of the variable
+		// Create a matrix of scalar polynomials
 		var rightAsVM VariableMatrix = make([][]Variable, nR)
 		for rIndex := 0; rIndex < nR; rIndex++ {
 			rightAsVM[rIndex] = make([]Variable, nC)
@@ -141,47 +142,28 @@ func (km KMatrix) Plus(e interface{}) Expression {
 			}
 		}
 
-		return km.Plus(rightAsVM) // Reuse VariableMatrix case
+		out = km.Plus(rightAsVM) // Reuse VariableMatrix case
 
 	case mat.Dense:
-		return km.Plus(DenseToKMatrix(right)) // Reuse KMatrix case
+		out = km.Plus(DenseToKMatrix(right)) // Reuse KMatrix case
 
 	case *mat.Dense:
-		return km.Plus(*right) // Reuse mat.Dense case
-
-	case KMatrix:
-		// Create the result matrix
-		var result KMatrix = make([][]K, nR)
-		for rIndex := 0; rIndex < nR; rIndex++ {
-			result[rIndex] = make([]K, nC)
-			for cIndex := 0; cIndex < nC; cIndex++ {
-				result[rIndex][cIndex] = km[rIndex][cIndex] + right[rIndex][cIndex]
-			}
-		}
-		return result
+		out = km.Plus(*right) // Reuse mat.Dense case
 
-	case VariableMatrix:
-		// Create the result matrix
-		var result PolynomialMatrix = make([][]Polynomial, nR)
-		for rIndex := 0; rIndex < nR; rIndex++ {
-			result[rIndex] = make([]Polynomial, nC)
-			for cIndex := 0; cIndex < nC; cIndex++ {
-				result[rIndex][cIndex] = km[rIndex][cIndex].Plus(right[rIndex][cIndex]).(Polynomial)
-				// Each addition should create a polynomial
-			}
-		}
-		return result
-	case PolynomialMatrix:
-		return right.Plus(km) // Reuse PolynomialMatrix case
+	case MatrixExpression:
+		out = MatrixPlusTemplate(km, right)
+	default:
+		// If we reach this point, the input is not recognized
+		panic(
+			smErrors.UnsupportedInputError{
+				FunctionName: "KMatrix.Plus",
+				Input:        e,
+			},
+		)
 	}
 
-	// If we reach this point, the input is not recognized
-	panic(
-		smErrors.UnsupportedInputError{
-			FunctionName: "KMatrix.Plus",
-			Input:        e,
-		},
-	)
+	// Simplify and return
+	return out.AsSimplifiedExpression()
 }
 
 /*
@@ -266,42 +248,31 @@ func (km KMatrix) Multiply(e interface{}) Expression {
 		}
 	}
 
+	// Algorithm
+	var out Expression
 	switch right := e.(type) {
 	case float64:
 		// Use gonum's built-in scale function
 		kmAsDense := km.ToDense()
 		var product mat.Dense
 		product.Scale(right, &kmAsDense)
 
-		return DenseToKMatrix(product)
-
+		out = DenseToKMatrix(product)
 	case K:
-		return km.Multiply(float64(right)) // Reuse float64 case
+		out = km.Multiply(float64(right)) // Reuse float64 case
 	case Polynomial:
 		// Choose the correct output type based on the size of km
 		nR, nC := km.Dims()[0], km.Dims()[1]
-		switch {
-		case (nR == 1) && (nC == 1):
-			// If the output is a scalar, return a scalar
-			return km[0][0].Multiply(right)
-		case nC == 1:
-			// If the output is a vector, return a vector
-			var outputVec PolynomialVector = make([]Polynomial, nR)
-			for rIndex := 0; rIndex < nR; rIndex++ {
-				outputVec[rIndex] = km[rIndex][0].Multiply(right.Copy()).(Polynomial)
+		var product [][]ScalarExpression
+		for ii := 0; ii < nR; ii++ {
+			var productRow []ScalarExpression
+			for jj := 0; jj < nC; jj++ {
+				kmAsPoly := km[ii][jj].ToPolynomial()
+				productRow = append(productRow, kmAsPoly.Multiply(right).(ScalarExpression))
 			}
-			return outputVec
-		default:
-			// If the output is a matrix, return a matrix
-			var outputMat PolynomialMatrix = make([][]Polynomial, nR)
-			for rIndex := 0; rIndex < nR; rIndex++ {
-				outputMat[rIndex] = make([]Polynomial, nC)
-				for cIndex := 0; cIndex < nC; cIndex++ {
-					outputMat[rIndex][cIndex] = km[rIndex][cIndex].Multiply(right.Copy()).(Polynomial)
-				}
-			}
-			return outputMat
+			product = append(product, productRow)
 		}
+		out = ConcretizeExpression(product)
 
 	case *mat.VecDense:
 		// Use gonum's built-in multiplication function
@@ -313,77 +284,73 @@ func (km KMatrix) Multiply(e interface{}) Expression {
 		nOutput := product.Len()
 		if nOutput == 1 {
 			// If the output is a scalar, return a scalar
-			return K(product.AtVec(0))
+			out = K(product.AtVec(0))
 		} else {
 			// Otherwsie return a KVector
-			return VecDenseToKVector(product)
+			out = VecDenseToKVector(product)
 		}
 
 	case VariableVector:
 		// Choose the correct output type based on the size of km
 		nR := km.Dims()[0]
 		if nR == 1 {
 			// If the output is a scalar, return a scalar
-			var out Polynomial = K(0).ToPolynomial()
+			var prod Expression = K(0)
 			for cIndex := 0; cIndex < len(right); cIndex++ {
-				out = out.Plus(
-					right[cIndex].Multiply(km[0][cIndex]),
-				).(Polynomial)
+				prod = prod.Plus(right[cIndex].Multiply(km[0][cIndex]))
 			}
-			return out
+			out = prod
 		} else {
 			nC := km.Dims()[1]
 			// If the output is a vector, return a vector
-			var outputVec PolynomialVector = VecDenseToKVector(
-				ZerosVector(nR),
-			).ToPolynomialVector()
+			var outputVec Expression = VecDenseToKVector(ZerosVector(nR))
 			for colIndex := 0; colIndex < nC; colIndex++ {
 				kmAsDense := km.ToDense()
 				tempCol := (&kmAsDense).ColView(colIndex)
-				outputVec = outputVec.Plus(
-					right[colIndex].Multiply(tempCol),
-				).(PolynomialVector)
+				outputVec = outputVec.Plus(right[colIndex].Multiply(tempCol))
 			}
-			return outputVec
+			out = outputVec
 		}
 	case *mat.Dense:
 		// Check output dimensions
 		nOutputR := km.Dims()[0]
 		_, nOutputCols := right.Dims()
 
 		kmAsDense := km.ToDense()
-		var out mat.Dense
-		out.Mul(&kmAsDense, right)
+		var prod mat.Dense
+		prod.Mul(&kmAsDense, right)
 
 		switch {
 		case nOutputR == 1 && nOutputCols == 1:
 			// If the constant matrix is a scalar, return the scalar
-			return K(out.At(0, 0))
+			out = K(prod.At(0, 0))
 		case nOutputCols == 1:
 			// If the output is a vector, return a vector
 			var outputVec KVector = make([]K, nOutputR)
 			for rIndex := 0; rIndex < nOutputR; rIndex++ {
-				outputVec[rIndex] = K(out.At(rIndex, 0))
+				outputVec[rIndex] = K(prod.At(rIndex, 0))
 			}
-			return outputVec
+			out = outputVec
 		default:
 			// If the output is a matrix, return a matrix
-			return DenseToKMatrix(out)
+			out = DenseToKMatrix(prod)
 		}
 	case mat.Dense:
 		// Use *mat.Dense method
-		return km.Multiply(&right) // Reuse *mat.Dense case
+		out = km.Multiply(&right) // Reuse *mat.Dense case
 	case MatrixExpression:
-		return MatrixMultiplyTemplate(km, right)
+		out = MatrixMultiplyTemplate(km, right)
+	default:
+		// If we reach this point, the input is not recognized
+		panic(
+			smErrors.UnsupportedInputError{
+				FunctionName: "KMatrix.Multiply",
+				Input:        e,
+			},
+		)
 	}
 
-	// If we reach this point, the input is not recognized
-	panic(
-		smErrors.UnsupportedInputError{
-			FunctionName: "KMatrix.Multiply",
-			Input:        e,
-		},
-	)
+	return out.AsSimplifiedExpression()
 }
 
 /*