-
Notifications
You must be signed in to change notification settings - Fork 0
/
Q3.py
375 lines (321 loc) · 11.9 KB
/
Q3.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
"""
111901030
Mayank Singla
Coding Assignment 3 - Q3
"""
# %%
from random import uniform
import numpy as np
import matplotlib.pyplot as plt
def handleError(method):
"""
Decorator Factory function.
Returns a decorator that normally calls the method of a class by forwarding all its arguments to the method.
It surrounds the method calling in try-except block to handle errors gracefully.
"""
def decorator(ref, *args, **kwargs):
"""
Decorator function that surrounds the method of a class in try-except block and call the methods and handles error gracefully.
"""
try:
# Return the same value as that of the method if any
return method(ref, *args, **kwargs)
except Exception as err:
print(type(err))
print(err)
return decorator
class RowVectorFloat:
"""
Represents a row vector of float values.
"""
@handleError
def _validateListValues(self, lst):
"""
Validates the values in the list.
Returns True if correct.
"""
if not isinstance(lst, list):
raise Exception("Invalid input - Expected list")
for i in lst:
if not isinstance(i, float) and not isinstance(i, int):
raise Exception(
f"Invalid type of value received {type(i)}.\nExpected float or int."
)
return True
@handleError
def _validateIndex(self, index):
"""
Validates the input index.
Returns True if correct.
"""
if not isinstance(index, int):
raise Exception(
f"Invalid type of index received {type(index)}.\nExpected int."
)
n = len(self.vec)
if index >= n or index < (-n):
raise Exception(f"Index out of range.")
return True
@handleError
def __init__(self, lst):
"""
Initializes the row vector with the values in the list.
"""
if not self._validateListValues(lst):
return
# Creating a new copy of the list
self.vec = list(lst)
@handleError
def __str__(self):
"""
Returns the string representation of the row vector.
"""
return " ".join(
f"{i:.2f}" if isinstance(i, float) else str(i) for i in self.vec
)
@handleError
def __len__(self):
"""
Returns the length of the row vector.
"""
return len(self.vec)
@handleError
def __getitem__(self, index):
"""
Returns the value at the given index.
"""
if not self._validateIndex(index):
return
return self.vec[index]
@handleError
def __setitem__(self, index, value):
"""
Sets the value at the given index
"""
if not self._validateIndex(index):
return
elif not self._validateListValues([value]):
return
self.vec[index] = value
@handleError
def __add__(self, rv):
"""
Adds two row vectors.
Operator looks for __add__ in left operand.
"""
if not isinstance(rv, RowVectorFloat):
raise Exception("Invalid input - Expected RowVectorFloat")
elif len(self) != len(rv):
raise Exception("Invalid input - Expected same length vectors")
ans = [self.vec[i] + rv.vec[i] for i in range(len(self))]
return RowVectorFloat(ans)
@handleError
def __radd__(self, rv):
"""
Adds two row vectors.
Operator looks for __add__ in right operand.
"""
return self.__add__(rv)
@handleError
def __mul__(self, scalar):
"""
Multiplies a row vector with a scalar.
Operator looks for __mul__ in left operand.
"""
if not isinstance(scalar, (int, float)):
raise Exception("Invalid input - Expected scalar")
ans = [
self.vec[i] * scalar if self.vec[i] != 0 else 0.00 for i in range(len(self))
]
return RowVectorFloat(ans)
@handleError
def __rmul__(self, scalar):
"""
Multiplies a row vector with a scalar.
Operator looks for __mul__ in right operand.
"""
return self.__mul__(scalar)
class SquareMatrixFloat:
"""
Represents a square matrix
"""
@handleError
def __init__(self, n):
"""
Initializes the square matrix as a list of RowVectorFloat objects.
"""
if not isinstance(n, int):
raise Exception("Invalid input - Expected int")
self.n = n
self.mat = [RowVectorFloat([0] * n) for _ in range(n)]
@handleError
def __str__(self):
"""
Returns the string representation of the square matrix.
"""
ans = "The matrix is:\n"
for i in range(self.n):
ans += str(self.mat[i]) + "\n"
return ans
@handleError
def sampleSymmetric(self):
"""
Samples a random symmetric matrix.
"""
for i in range(self.n):
for j in range(i):
self.mat[i][j] = uniform(0, 1)
self.mat[j][i] = self.mat[i][j]
self.mat[i][i] = uniform(0, self.n)
@handleError
def toRowEchelonForm(self):
"""
Converts the matrix to its row echelon form.
"""
pivotRow, pivotCol = 0, 0
while pivotRow < self.n and pivotCol < self.n:
# Find the row with first non-zero entry in the pivot column
nonZeroRow = pivotRow
while nonZeroRow < self.n and self.mat[nonZeroRow][pivotCol] == 0:
nonZeroRow += 1
# If no non-zero entry found, move to next pivot column
if nonZeroRow == self.n:
pivotCol += 1
continue
# Swap the pivot row with the nonZeroRow
if nonZeroRow != pivotRow:
self.mat[pivotRow], self.mat[nonZeroRow] = (
self.mat[nonZeroRow],
self.mat[pivotRow],
)
# Multiply each element in the pivot row by the inverse of the pivot, so the pivot equals 1.
self.mat[pivotRow] = self.mat[pivotRow] * (1 / self.mat[pivotRow][pivotCol])
# Make the pivot entry 1 (to avoid python precision issues)
self.mat[pivotRow][pivotCol] = 1.00
# Add multiples of the pivot row to each of the lower rows, so every element in the pivot column of the lower rows equals 0.
for i in range(pivotRow + 1, self.n):
if self.mat[i][pivotCol] != 0:
self.mat[i] = (
self.mat[i] + (-self.mat[i][pivotCol]) * self.mat[pivotRow]
)
self.mat[i][pivotCol] = 0.00 # To avoid python precision issues
# Move to the next column and next row
pivotCol += 1
pivotRow += 1
return self
@handleError
def isDRDominant(self):
"""
Checks if the matrix is Diagonally row dominant.
Here, I am checking strictly for DRDominant though it was not mentioned in the question, but it is required for checking the correctness of Jacobi method.
"""
for i in range(self.n):
sumRemRow = -self.mat[i][i]
for j in range(len(self.mat[i])):
sumRemRow += self.mat[i][j]
if self.mat[i][i] <= sumRemRow:
return False
return True
@handleError
def _validateListValues(self, lst):
"""
Validates the values in the list.
Returns True if correct.
"""
if not isinstance(lst, list):
raise Exception("Invalid input - Expected list")
for i in lst:
if not isinstance(i, float) and not isinstance(i, int):
raise Exception(
f"Invalid type of value received {type(i)}.\nExpected float or int."
)
return True
@handleError
def _iterativeMethod(self, b, m, isJacobi=True):
"""
Takes a list (denoting vector `b`) and number of iterations `m` as its arguments, and performs `m` iterations of the Jacobi / Gauss-Siedel iterative procedure.
Return the final iteration value, and value of the term `||Ax⁽ᵏ⁾ - b||₂` from all the m iterations
"""
if not self._validateListValues(b):
return
if not isinstance(m, int) or m < 1:
raise Exception(f"Invalid m received {m}.\nExpected a positive int.")
if isJacobi and not self.isDRDominant():
raise Exception("Not solving because convergence is not guranteed.")
# Converting the class matrix to numpy array
A = []
for i in range(self.n):
el = []
for j in range(self.n):
el.append(self.mat[i][j])
A.append(el)
A = np.array(A)
# Initialize the iteration vector
prevX = [0] * self.n
x = [0] * self.n
ansError = []
# Doing m iterations
for _ in range(m):
# Computing each xᵢ
for i in range(self.n):
sumRemRow = 0
for j in range(self.n):
if isJacobi:
if j != i:
sumRemRow += self.mat[i][j] * prevX[j]
else:
if j < i:
sumRemRow += self.mat[i][j] * x[j]
elif j > i:
sumRemRow += self.mat[i][j] * prevX[j]
x[i] = (b[i] - sumRemRow) / self.mat[i][i]
# ||Ax⁽ᵏ⁾ - b||₂
ansError.append(np.linalg.norm(A @ np.array(x) - np.array(b)))
prevX = x[:]
ansX = prevX[:] # Final iteration value
return (ansError, ansX)
@handleError
def jSolve(self, b, m):
"""
Takes a list (denoting vector `b`) and number of iterations `m` as its arguments, and performs `m` iterations of the Jacobi iterative procedure.
Return the final iteration value, and value of the term `||Ax⁽ᵏ⁾ - b||₂` from all the m iterations
"""
return self._iterativeMethod(b, m, isJacobi=True)
@handleError
def gsSolve(self, b, m):
"""
Takes a list (denoting vector `b`) and number of iterations `m` as its arguments, and performs `m` iterations of the Gauss-Siedel iterative procedure.
Return the final iteration value, and value of the term `||Ax⁽ᵏ⁾ - b||₂` from all the m iterations
"""
return self._iterativeMethod(b, m, isJacobi=False)
@handleError
def visualizeRateOfConvergence(self):
"""
Visualize the rate of convergence of Jacobi and Gauss-Siedel methods of a linear system with a diagonally dominant square symmetric matrix.
"""
# Constructing a diagonally dominant square symmetric matrix
self.sampleSymmetric()
while not self.isDRDominant():
self.sampleSymmetric()
numIterations = 50 # Number of iterations used
b = list(range(self.n)) # Vector b
# Solving the system using Jacobi method
jacobiError, jacobiX = self.jSolve(b, numIterations)
# Solving the system using Gauss-Siedel method
gsError, gsX = self.gsSolve(b, numIterations)
xpoints = list(range(1, numIterations + 1))
# Plotting the rate of convergence
plt.title(
f"Rate of convergence of Jacobi and Gauss-Siedel methods\n({numIterations} iterations, {self.n}x{self.n} matrix)"
)
plt.xlabel("Number of iterations")
plt.ylabel("Error : ||Ax⁽ᵏ⁾ - b||₂")
plt.plot(xpoints, jacobiError, c="b", label="Jacobi Method")
plt.plot(xpoints, gsError, c="r", label="Gauss-Siedel Method")
plt.legend()
plt.grid()
plt.show()
if __name__ == "__main__":
# Sample Test Case
s = SquareMatrixFloat(10)
s.visualizeRateOfConvergence()