From c0997609a594f7b0c8d269fa8819ca27b5d44db8 Mon Sep 17 00:00:00 2001 From: jasmaa Date: Wed, 5 Sep 2018 22:14:06 -0400 Subject: [PATCH 1/5] initial implementation --- .../code/python/gaussian_elimination.py | 75 +++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 contents/gaussian_elimination/code/python/gaussian_elimination.py diff --git a/contents/gaussian_elimination/code/python/gaussian_elimination.py b/contents/gaussian_elimination/code/python/gaussian_elimination.py new file mode 100644 index 000000000..d7c644397 --- /dev/null +++ b/contents/gaussian_elimination/code/python/gaussian_elimination.py @@ -0,0 +1,75 @@ +import numpy as np + +def gaussian_elimination(A): + # Step 1: Go by number of pivots with end being the min number of rows or cols + for pivot_i in range(min(A.shape[0], A.shape[1])): + + # Step 2: Swap row with highest element in col + max_i = np.argmax(abs(A[pivot_i:,pivot_i])) + pivot_i + + temp = A[pivot_i,:].copy() + A[pivot_i,:] = A[max_i,:] + A[max_i,:] = temp + + if A[pivot_i, pivot_i] == 0: + # Skip eval on singular matrix + continue + + for row in range(pivot_i+1, A.shape[0]): + # Step 3: Get fraction + frac = -A[row,pivot_i] / A[pivot_i,pivot_i] + # Step 4: Add rows + A[row,:] = A[row,:] + frac*A[pivot_i,:] + +""" +Assumes A is already ref +""" +def gauss_jordon_elimination(A): + # Find and set each pivot to one via row scaling + col = 0 + for row in range(A.shape[0]): + while A[row,col] == 0: + col += 1 + if col >= A.shape[1]: + break + if col >= A.shape[1]: + continue + A[row,:] = A[row,:] / A[row,col] + + # kill rows above + for r in range(row): + A[r,:] = A[r,:] - A[r,col] * A[row,:] + + +""" +Assumes A has a unique solution and A in ref +""" +def back_sub(A): + sol = np.zeros(A.shape[0]).T + for pivot_i in range(A.shape[0]-1, -1, -1): + s = 0 + for col in range(pivot_i+1, A.shape[1]-1): + s += A[pivot_i,col] * sol[col] + sol[pivot_i] = (A[pivot_i,A.shape[1]-1] - s) / A[pivot_i,pivot_i] + return sol + + +# MAIN +A = np.matrix('0. 2 1 -8; 1 -2 -3 0; -1 1 2 3') +#A = np.matrix('1. -2 -6 12; 2 4 12 -17; 1 -4 -12 22') +print(A) +print() + +gaussian_elimination(A) +print("ref") +print(A) +print() + +print("back sub") +print(back_sub(A)) +print() + +gauss_jordon_elimination(A) +print("rref") +print(A) +print() From a5d48102e6c21c07a1da05379be9e2048904d91a Mon Sep 17 00:00:00 2001 From: jasmaa Date: Thu, 6 Sep 2018 20:44:19 -0400 Subject: [PATCH 2/5] Debug ref and rref --- .../code/python/gaussian_elimination.py | 84 ++++++++++++------- 1 file changed, 52 insertions(+), 32 deletions(-) diff --git a/contents/gaussian_elimination/code/python/gaussian_elimination.py b/contents/gaussian_elimination/code/python/gaussian_elimination.py index d7c644397..bc20a6517 100644 --- a/contents/gaussian_elimination/code/python/gaussian_elimination.py +++ b/contents/gaussian_elimination/code/python/gaussian_elimination.py @@ -1,42 +1,52 @@ import numpy as np def gaussian_elimination(A): - # Step 1: Go by number of pivots with end being the min number of rows or cols - for pivot_i in range(min(A.shape[0], A.shape[1])): + + pivot_row = 0 + + # Step 1: Go by column + for pivot_col in range(min(A.shape[0], A.shape[1])): # Step 2: Swap row with highest element in col - max_i = np.argmax(abs(A[pivot_i:,pivot_i])) + pivot_i + max_i = np.argmax(abs(A[pivot_row:,pivot_col])) + pivot_row - temp = A[pivot_i,:].copy() - A[pivot_i,:] = A[max_i,:] + temp = A[pivot_row,:].copy() + A[pivot_row,:] = A[max_i,:] A[max_i,:] = temp - if A[pivot_i, pivot_i] == 0: - # Skip eval on singular matrix + # Skip on singular matrix, not actually a pivot + if A[pivot_row, pivot_col] == 0: continue - for row in range(pivot_i+1, A.shape[0]): + # Steps 3 & 4: Zero out elements below pivot + for r in range(pivot_row+1, A.shape[0]): # Step 3: Get fraction - frac = -A[row,pivot_i] / A[pivot_i,pivot_i] + frac = -A[r,pivot_col] / A[pivot_row,pivot_col] # Step 4: Add rows - A[row,:] = A[row,:] + frac*A[pivot_i,:] + A[r,:] = A[r,:] + frac*A[pivot_row,:] + + pivot_row += 1 + """ Assumes A is already ref """ def gauss_jordon_elimination(A): - # Find and set each pivot to one via row scaling + col = 0 + + # Scan for pivots for row in range(A.shape[0]): - while A[row,col] == 0: + while col < A.shape[1] and A[row,col] == 0: col += 1 - if col >= A.shape[1]: - break + if col >= A.shape[1]: continue + + # Set each pivot to one via row scaling A[row,:] = A[row,:] / A[row,col] - # kill rows above + # Zero out elements above pivot for r in range(row): A[r,:] = A[r,:] - A[r,col] * A[row,:] @@ -44,32 +54,42 @@ def gauss_jordon_elimination(A): """ Assumes A has a unique solution and A in ref """ -def back_sub(A): +def back_substitution(A): + sol = np.zeros(A.shape[0]).T + + # Go by pivots along diagonal for pivot_i in range(A.shape[0]-1, -1, -1): s = 0 for col in range(pivot_i+1, A.shape[1]-1): s += A[pivot_i,col] * sol[col] sol[pivot_i] = (A[pivot_i,A.shape[1]-1] - s) / A[pivot_i,pivot_i] + return sol -# MAIN -A = np.matrix('0. 2 1 -8; 1 -2 -3 0; -1 1 2 3') -#A = np.matrix('1. -2 -6 12; 2 4 12 -17; 1 -4 -12 22') -print(A) -print() +def main(): + #A = np.matrix('0. 2 1 -8; 1 -2 -3 0; -1 1 2 3') + #A = np.matrix('1. -2 -6 12 1; 2 4 12 -17 1; 1 -4 -12 22 1') + #A = np.matrix('1 1 1 1 1; 0 1 0 0 0; 0 0 0 0 1; 0 0 0 1 0') + A = np.matrix('2. 3 4 6; 1 2 3 4; 3 -4 0 10') + + print("Original") + print(A, "\n") + + gaussian_elimination(A) + print("Gaussian elimination, REF") + print(A, "\n") + + print("Back subsitution") + print(back_substitution(A)) + print() + + gauss_jordon_elimination(A) + print("Gauss-Jordan, RREF") + print(A, "\n") -gaussian_elimination(A) -print("ref") -print(A) -print() -print("back sub") -print(back_sub(A)) -print() +if __name__ == "__main__": + main() -gauss_jordon_elimination(A) -print("rref") -print(A) -print() From 8aa5fa98d639d4255e19719dd1ab55e77ae67fe3 Mon Sep 17 00:00:00 2001 From: jasmaa Date: Thu, 6 Sep 2018 20:50:53 -0400 Subject: [PATCH 3/5] Add to gitbook --- .../code/python/gaussian_elimination.py | 3 --- contents/gaussian_elimination/gaussian_elimination.md | 8 ++++++++ 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/contents/gaussian_elimination/code/python/gaussian_elimination.py b/contents/gaussian_elimination/code/python/gaussian_elimination.py index bc20a6517..f03eaef46 100644 --- a/contents/gaussian_elimination/code/python/gaussian_elimination.py +++ b/contents/gaussian_elimination/code/python/gaussian_elimination.py @@ -69,9 +69,6 @@ def back_substitution(A): def main(): - #A = np.matrix('0. 2 1 -8; 1 -2 -3 0; -1 1 2 3') - #A = np.matrix('1. -2 -6 12 1; 2 4 12 -17 1; 1 -4 -12 22 1') - #A = np.matrix('1 1 1 1 1; 0 1 0 0 0; 0 0 0 0 1; 0 0 0 1 0') A = np.matrix('2. 3 4 6; 1 2 3 4; 3 -4 0 10') print("Original") diff --git a/contents/gaussian_elimination/gaussian_elimination.md b/contents/gaussian_elimination/gaussian_elimination.md index d48b46df4..aeed20f81 100644 --- a/contents/gaussian_elimination/gaussian_elimination.md +++ b/contents/gaussian_elimination/gaussian_elimination.md @@ -360,6 +360,8 @@ In code, this looks like: [import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs) {% sample lang="hs" %} [import:10-36, lang:"haskell"](code/haskell/gaussianElimination.hs) +{% sample lang="py" %} +[import:3-28, lang:"python"](code/python/gaussian_elimination.py) {% endmethod %} Now, to be clear: this algorithm creates an upper-triangular matrix. @@ -393,6 +395,8 @@ This code does not exist yet in rust, so here's Julia code (sorry for the inconv [import:67-93, lang:"julia"](code/julia/gaussian_elimination.jl) {% sample lang="hs" %} [import:38-46, lang:"haskell"](code/haskell/gaussianElimination.hs) +{% sample lang="py" %} +[import:31-51, lang:"python"](code/python/gaussian_elimination.py) {% endmethod %} ## Back-substitution @@ -423,6 +427,8 @@ In code, this involves keeping a rolling sum of all the values we substitute in [import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs) {% sample lang="hs" %} [import:48-53, lang:"haskell"](code/haskell/gaussianElimination.hs) +{% sample lang="py" %} +[import:54-68, lang:"python"](code/python/gaussian_elimination.py) {% endmethod %} ## Conclusions @@ -445,6 +451,8 @@ As for what's next... Well, we are in for a treat! The above algorithm clearly h [import, lang:"rust"](code/rust/gaussian_elimination.rs) {% sample lang="hs" %} [import, lang:"haskell"](code/haskell/gaussianElimination.hs) +{% sample lang="py" %} +[import, lang:"python"](code/python/gaussian_elimination.py) {% endmethod %} From 70a487b89f593c753cfdf01774338b571807aa02 Mon Sep 17 00:00:00 2001 From: jasmaa Date: Fri, 14 Sep 2018 12:16:45 -0400 Subject: [PATCH 4/5] Make requested changes --- .../code/python/gaussian_elimination.py | 53 +++++++++---------- .../gaussian_elimination.md | 4 +- 2 files changed, 27 insertions(+), 30 deletions(-) diff --git a/contents/gaussian_elimination/code/python/gaussian_elimination.py b/contents/gaussian_elimination/code/python/gaussian_elimination.py index f03eaef46..5fe36548e 100644 --- a/contents/gaussian_elimination/code/python/gaussian_elimination.py +++ b/contents/gaussian_elimination/code/python/gaussian_elimination.py @@ -8,82 +8,79 @@ def gaussian_elimination(A): for pivot_col in range(min(A.shape[0], A.shape[1])): # Step 2: Swap row with highest element in col - max_i = np.argmax(abs(A[pivot_row:,pivot_col])) + pivot_row + max_i = np.argmax(abs(A[pivot_row:, pivot_col])) + pivot_row - temp = A[pivot_row,:].copy() - A[pivot_row,:] = A[max_i,:] - A[max_i,:] = temp + temp = A[pivot_row, :].copy() + A[pivot_row, :] = A[max_i, :] + A[max_i, :] = temp - # Skip on singular matrix, not actually a pivot + # Skip on singular matrix, not actually a pivot if A[pivot_row, pivot_col] == 0: continue # Steps 3 & 4: Zero out elements below pivot - for r in range(pivot_row+1, A.shape[0]): + for r in range(pivot_row + 1, A.shape[0]): # Step 3: Get fraction - frac = -A[r,pivot_col] / A[pivot_row,pivot_col] + frac = -A[r, pivot_col] / A[pivot_row, pivot_col] # Step 4: Add rows - A[r,:] = A[r,:] + frac*A[pivot_row,:] + A[r, :] += frac * A[pivot_row, :] pivot_row += 1 -""" -Assumes A is already ref -""" -def gauss_jordon_elimination(A): +# Assumes A is already ref +def gauss_jordan_elimination(A): col = 0 # Scan for pivots for row in range(A.shape[0]): - while col < A.shape[1] and A[row,col] == 0: + while col < A.shape[1] and A[row, col] == 0: col += 1 if col >= A.shape[1]: continue # Set each pivot to one via row scaling - A[row,:] = A[row,:] / A[row,col] + A[row, :] /= A[row, col] # Zero out elements above pivot for r in range(row): - A[r,:] = A[r,:] - A[r,col] * A[row,:] + A[r, :] -= A[r, col] * A[row, :] -""" -Assumes A has a unique solution and A in ref -""" +# Assumes A has a unique solution and A in ref def back_substitution(A): sol = np.zeros(A.shape[0]).T # Go by pivots along diagonal - for pivot_i in range(A.shape[0]-1, -1, -1): + for pivot_i in range(A.shape[0] - 1, -1, -1): s = 0 - for col in range(pivot_i+1, A.shape[1]-1): - s += A[pivot_i,col] * sol[col] - sol[pivot_i] = (A[pivot_i,A.shape[1]-1] - s) / A[pivot_i,pivot_i] + for col in range(pivot_i + 1, A.shape[1] - 1): + s += A[pivot_i, col] * sol[col] + sol[pivot_i] = (A[pivot_i, A.shape[1] - 1] - s) / A[pivot_i, pivot_i] return sol def main(): - A = np.matrix('2. 3 4 6; 1 2 3 4; 3 -4 0 10') + A = np.array([[2, 3, 4, 6], + [1, 2, 3, 4,], + [3, -4, 0, 10]], dtype=float) print("Original") print(A, "\n") gaussian_elimination(A) - print("Gaussian elimination, REF") + print("Gaussian elimination") print(A, "\n") print("Back subsitution") - print(back_substitution(A)) - print() + print(back_substitution(A), "\n") - gauss_jordon_elimination(A) - print("Gauss-Jordan, RREF") + gauss_jordan_elimination(A) + print("Gauss-Jordan") print(A, "\n") diff --git a/contents/gaussian_elimination/gaussian_elimination.md b/contents/gaussian_elimination/gaussian_elimination.md index aeed20f81..7976a0c7e 100644 --- a/contents/gaussian_elimination/gaussian_elimination.md +++ b/contents/gaussian_elimination/gaussian_elimination.md @@ -396,7 +396,7 @@ This code does not exist yet in rust, so here's Julia code (sorry for the inconv {% sample lang="hs" %} [import:38-46, lang:"haskell"](code/haskell/gaussianElimination.hs) {% sample lang="py" %} -[import:31-51, lang:"python"](code/python/gaussian_elimination.py) +[import:31-49, lang:"python"](code/python/gaussian_elimination.py) {% endmethod %} ## Back-substitution @@ -428,7 +428,7 @@ In code, this involves keeping a rolling sum of all the values we substitute in {% sample lang="hs" %} [import:48-53, lang:"haskell"](code/haskell/gaussianElimination.hs) {% sample lang="py" %} -[import:54-68, lang:"python"](code/python/gaussian_elimination.py) +[import:52-64, lang:"python"](code/python/gaussian_elimination.py) {% endmethod %} ## Conclusions From 290fa377047f8df610f243b323f8e9269932cc37 Mon Sep 17 00:00:00 2001 From: jasmaa Date: Sat, 22 Sep 2018 21:56:37 -0400 Subject: [PATCH 5/5] Change ref to row echelon form --- .../gaussian_elimination/code/python/gaussian_elimination.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/contents/gaussian_elimination/code/python/gaussian_elimination.py b/contents/gaussian_elimination/code/python/gaussian_elimination.py index 5fe36548e..3f11cd761 100644 --- a/contents/gaussian_elimination/code/python/gaussian_elimination.py +++ b/contents/gaussian_elimination/code/python/gaussian_elimination.py @@ -28,7 +28,7 @@ def gaussian_elimination(A): pivot_row += 1 -# Assumes A is already ref +# Assumes A is already row echelon form def gauss_jordan_elimination(A): col = 0 @@ -49,7 +49,7 @@ def gauss_jordan_elimination(A): A[r, :] -= A[r, col] * A[row, :] -# Assumes A has a unique solution and A in ref +# Assumes A has a unique solution and A in row echelon form def back_substitution(A): sol = np.zeros(A.shape[0]).T