Tarea 1

Nombre: Miguel Rodrigo Castañeda Santiago

175840

Clase: Array

Descripción: Clase para algebra lineal



In [1]:

    
class Array:
    def __init__(self, list_of_rows):
        "Constructor y validador"
        # obtener dimensiones
        self.data = list_of_rows
        nrow = len(list_of_rows)
        #  ___caso vector: redimensionar correctamente
        if not isinstance(list_of_rows[0], list):
            nrow = 1
            self.data = [[x] for x in list_of_rows]
        # ahora las columnas deben estar bien aunque sea un vector
        ncol = len(self.data[0])
        self.shape = (nrow, ncol)
        # validar tamano correcto de filas
        if any([len(r) != ncol for r in self.data]):
            raise Exception("Las filas deben ser del mismo tamano")
        for r in range(nrow):
            for c in range(ncol):
                if not isinstance(self.data[r][c],(int, float, complex)):
                    raise Exception("Los valores deben ser numéricos")

    def __add__(self, other):
        "Suma de matrices"
        if isinstance(other, Array):
            if self.shape != other.shape:
                raise Exception("Las dimensiones son distintas!")
            rows, cols = self.shape
            newArray = Array([[0. for c in range(cols)] for r in range(rows)])
            for r in range(rows):
                for c in range(cols):
                    newArray.data[r][c] = self.data[r][c] + other.data[r][c]
            return newArray
        elif isinstance(other, (int, float, complex)): # en caso de que el lado derecho sea solo un numero
            rows, cols = self.shape
            newArray = Array([[0. for c in range(cols)] for r in range(rows)])
            for r in range(rows):
                for c in range(cols):
                    newArray.data[r][c] = self.data[r][c] + other
            return newArray
        else:
            return NotImplemented # es un tipo de error particular usado en estos metodos

    def __mul__(self, other):
        if isinstance(other, (int,float, complex)): # Multiplicacion escalar
            rows, cols = self.shape
            newArray = Array.zeros(self.shape)
            for r in range(rows):
                for c in range(cols):
                    newArray.data[r][c] = self.data[r][c] * other
            return newArray
        if isinstance(other, Array):
            rows, cols = self.shape
            rowsOther, colsOther = other.shape
            if cols != rowsOther:
                raise Exception ("Las matrices no se pueden multiplicar")
            newArray = Array.zeros((rows,colsOther))
            for i in range(rows):
                for j in range(colsOther):
                    for k in range(cols):
                        newArray.data[i][j] += self.data[i][k] * other.data[k][j]
            return newArray
        else:
            return NotImplemented

    def __rmul__(self, other):
        return self.__mul__(other)

    def __sub__(self, other):
        "Resta de matrices"
        if isinstance(other, Array):
            if self.shape != other.shape:
                raise Exception("Las dimensiones son distintas!")
            rows, cols = self.shape
            newArray = Array([[0. for c in range(cols)] for r in range(rows)])
            for r in range(rows):
                for c in range(cols):
                    newArray.data[r][c] = self.data[r][c] - other.data[r][c]
            return newArray
        elif isinstance(other, (int, float, complex)):
            rows, cols = self.shape
            newArray = Array([[0. for c in range(cols)] for r in range(rows)])
            for r in range(rows):
                for c in range(cols):
                    newArray.data[r][c] = self.data[r][c] - other
            return newArray
        else:
            return NotImplemented

    def __radd__(self, other):
        return self.__add__(other)

    def __rsub__(self, other):
        return self.__sub__(other)

    def __repr__(self):
        return self.__class__.__name__ + "(" + self.__str__() + ")"

    def __str__(self):
        rows, cols = self.shape
        resultado = ""
        if rows > 1:
            resultado += "["
        for r in range(rows):
            resultado += "["
            for c in range(cols):
                resultado += str(self.data[r][c]) + " "
            if r + 1 == rows:
                resultado += "]"
            else:
                resultado += "]\n"
        if rows > 1:
            resultado += "]"
        return resultado

    def __getitem__(self, idx):
        return self.data[idx[0]][idx[1]]

    def __setitem__(self, idx, new_value):
        self.data[idx[0]][idx[1]] = new_value

    def zeros(shape):
        if not isinstance(shape,tuple):
            raise "Error en el tipo de dato"
        if len(shape) != 2:
            raise "Debe tener la forma (nrow, ncol) "
        nrow, ncol  = shape[0], shape[1]
        newArray = Array([[0. for c in range(ncol)] for r in range(nrow)])
        return newArray

    def eye(n):
        nrow, ncol = n,n
        shape = (n,n)
        newArray = Array.zeros(shape)
        for r in range(nrow):
            for c in range(ncol):
                if r == c:
                    newArray[r,c] = 1.
        return newArray

    def transpose(self):
        nrow , ncol = self.shape
        newArray = Array.zeros((ncol,nrow))
        for r in range(nrow):
            for c in range(ncol):
                newArray[c,r] = self.data[r][c]
        return newArray

    def forward_subs(L,b):
        rows,cols = L.shape
        if rows != cols:
            raise Exception("La matriz L debe ser cuadrada ...")

        rowb, colb = b.shape
        if colb != 1:
            raise Exception("No implementada la funcion de resolver más de un vector a la vez ...")

        n = rows
        for i in range(n - 2):
            for j in range(1,n):
                if L[i,j] != 0:
                    raise Exception("La matriz no es de la forma triangular inferior ...")

        x = Array.zeros((rowb,colb))

        for m in range(rowb):
            valor = 0.
            for i in range(0,m):
                valor += L[m,i] * x[i,0]
            x[m,0] = (b[m,0] - valor)/L[m,m]

        return x

    def backward_subs(U,b):
        rows, cols = U.shape
        if rows != cols:
            raise Exception("La matriz L debe ser cuadrada ...")

        rowb, colb = b.shape
        if colb != 1:
            raise Exception("No implementada la funcion de resolver más de un vector a la vez ...")

        n = rows
        for i in range(1, n  ):
            for j in range(0, n - 2):
                if U[i, j] != 0:
                    raise Exception("La matriz no es de la forma triangular superior ...")

        x = Array.zeros((rowb, colb))

        for m in range(rowb -1 ,-1,-1):
            valor = 0.
            for i in range(m+1, rowb):
                valor += U[m, i] * x[i, 0]
            x[m, 0] = (b[m, 0] - valor) / U[m, m]

        return x

    def pivot(self):
        rows, cols = self.shape
        newArray = Array.eye(rows)

        for j in range(rows):
            row = max(range(j, rows), key=lambda i: abs(self[i,j]))

            if j != row:
                for z in range(cols):
                    newArray[j, z], newArray[row, z] = newArray[row, z], newArray[j, z]
        return newArray

    def LU(self):
        # Doolittle algorithm
        rows, cols = self.shape
        if rows != cols:
            raise Exception("La matriz debe ser cuadrada ....")
        n = rows

        L = Array.zeros(self.shape)
        U = Array.zeros(self.shape)

        P = self.pivot()
        PA = P * self


        for j in range(n):
            L[j,j] = 1.0

            for i in range(j + 1):
                s1 = sum(U[k,j] * L[i,k] for k in range(i))
                U[i,j] = PA[i,j] - s1
            for i in range(j, n):
                s2 = sum(U[k,j] * L[i,k] for k in range(j))
                L[i,j] = (PA[i,j] - s2) / U[j,j]

        return (P, L, U)

    def lu_linsolve(A,b):
        P,L,U = A.LU()
        y = L.forward_subs(P * b)
        x = U.backward_subs(y)
        return x

    def imprimeFormatoWolfram(self):
        # imprime con formato de matrices Wolfram para comprobación
        rows, cols = self.shape
        resultado = ""
        if rows > 1:
            resultado += "{"
        for r in range(rows):
            resultado += "{"
            for c in range(cols):
                if c < cols - 1:
                    resultado += str(self.data[r][c]) + " , "
                else:
                    resultado += str(self.data[r][c]) + "  "
            if r + 1 == rows:
                resultado += "}"
            else:
                resultado += "},\n"
        if rows > 1:
            resultado += "}"
        print(resultado)

# Clase Vector
#
class Vector(Array):
    def __init__(self, list_of_numbers):
        self.vdata = list_of_numbers
        list_of_rows = [[x] for x in list_of_numbers]
        return Array.__init__(self, list_of_rows)

    def __repr__(self):
        return "Vector(" + str(self.vdata) + ")"

    def __str__(self):
        return str(self.vdata)

    def __add__(self, other):
        new_arr = Array.__add__(self, other)
        return Vector([x[0] for x in new_arr.data])

    def __radd__(self, other):
        new_arr = Array.__add__(self, other)
        return Vector([x[0] for x in new_arr.data])

Actividades

1. Un metodo para imprimir mejor...



In [6]:

    
print("Método para imprimir mejor")
A = Array([[-1,0,-5],[3,1,3],[1,1,1]])
print(A)









    



Método para imprimir mejor
[[-1 0 -5 ]
[3 1 3 ]
[1 1 1 ]]



In [7]:

    
A









    Out[7]:





Array([[-1 0 -5 ]
[3 1 3 ]
[1 1 1 ]])

2. Validador

Los valores deben ser numéricos, debe marcar error si se ingresa una letra



In [8]:

    
A = Array([['a',1,2],[1,-2,1],[-2,2,1]])









    



---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-8-f5dae7bfcd32> in <module>()
----> 1 A = Array([['a',1,2],[1,-2,1],[-2,2,1]])

<ipython-input-1-84211876ecfb> in __init__(self, list_of_rows)
     18             for c in range(ncol):
     19                 if not isinstance(self.data[r][c],(int, float, complex)):
---> 20                     raise Exception("Los valores deben ser numéricos")
     21 
     22     def __add__(self, other):

Exception: Los valores deben ser numéricos

3. Indexing and Item assignment



In [11]:

    
print("Método para obtener y asignar un valor de la matriz")
A = Array([[-1,0,-5],[3,1,3],[1,1,1]])
print(A)
print("\n se cambia el valor (1,1) a -10 ")
A[1,1]=-10
print(A)
print("Valor A[1,1] = " + str(A[1,1]))









    



Método para obtener y asignar un valor de la matriz
[[-1 0 -5 ]
[3 1 3 ]
[1 1 1 ]]

 se cambia el valor (1,1) a -10 
[[-1 0 -5 ]
[3 -10 3 ]
[1 1 1 ]]
Valor A[1,1] = -10

4. Iniciar una matriz en ceros y generar la matriz identidad



In [13]:

    
print("Iniciar una matriz en ceros")
Z = Array.zeros((3,3))
print(Z)









    



Iniciar una matriz en ceros
[[0.0 0.0 0.0 ]
[0.0 0.0 0.0 ]
[0.0 0.0 0.0 ]]



In [14]:

    
print("Iniciar una matriz identidad")
Z = Array.eye(3)
print(Z)









    



Iniciar una matriz identidad
[[1.0 0.0 0.0 ]
[0.0 1.0 0.0 ]
[0.0 0.0 1.0 ]]

5. Transpuesta



In [15]:

    
A = Array([[-1,0,-5],[3,1,3],[1,1,1]])
print("A:")
print(A)
At = A.transpose()
print("A Transpuesta ")
print(At)









    



A:
[[-1 0 -5 ]
[3 1 3 ]
[1 1 1 ]]
A Transpuesta 
[[-1 3 1 ]
[0 1 1 ]
[-5 3 1 ]]

6. Suma



In [17]:

    
A = Array([[-1,0,-5],[3,1,3],[1,1,1]])
B = Array([[1,2,6],[3,2,1],[-3,-2,1]])
print("A:")
print(A)
print("B:")
print(B)
C = A + B
print("Suma C = A + B")
print(C)









    



A:
[[-1 0 -5 ]
[3 1 3 ]
[1 1 1 ]]
B:
[[1 2 6 ]
[3 2 1 ]
[-3 -2 1 ]]
Suma C = A + B
[[0 2 1 ]
[6 3 4 ]
[-2 -1 2 ]]

7. Multiplicación matricial



In [18]:

    
A = Array([[-1,0,-5],[3,1,3],[1,1,1]])
B = Array([[1,2,6],[3,2,1],[-3,-2,1]])
valor = -3;

C = A * B
print("Multiplicacion A x B")
print(C)
print("Multiplicacion por escalar de la forma C = z x C")
C = 2 * A
print(C)

print("Multiplicacion por escalar de la forma C = C * z")
C = A * 2
print(C)









    



Multiplicacion A x B
[[14.0 8.0 -11.0 ]
[-3.0 2.0 22.0 ]
[1.0 2.0 8.0 ]]
Multiplicacion por escalar de la forma C = z x C
[[-2 0 -10 ]
[6 2 6 ]
[2 2 2 ]]
Multiplicacion por escalar de la forma C = C * z
[[-2 0 -10 ]
[6 2 6 ]
[2 2 2 ]]

8. Forward substitution, la matriz debe ser de la forma triangular inferior L



In [22]:

    
print("Forward substitution ")
A = Array([[-1,0,0],[-3,1,0],[1,-2,1]])
y = Array([[1],[-2],[-3]])
print("A:")
print(A)
print("y:")
print(y)
print("Resultado x:")
print(A.forward_subs(y))









    



Forward substitution 
A:
[[-1 0 0 ]
[-3 1 0 ]
[1 -2 1 ]]
y:
[[1 ]
[-2 ]
[-3 ]]
Resultado x:
[[-1.0 ]
[-5.0 ]
[-12.0 ]]

9. Back substitution, la matriz debe ser de la forma triangular superior U



In [24]:

    
print("back substitution ")
A = Array([[1,-1,2],[0,1,-3],[0,0,1]])
y = Array([[1],[-2],[-3]])
print("A:")
print(A)
print("y:")
print(y)
print("Resultado x:")
print(A.backward_subs(y))









    



back substitution 
A:
[[1 -1 2 ]
[0 1 -3 ]
[0 0 1 ]]
y:
[[1 ]
[-2 ]
[-3 ]]
Resultado x:
[[-4.0 ]
[-11.0 ]
[-3.0 ]]

10. Funcion LU, regresa tres matrices L, U, P tales que PA = LU, se utiliza el algoritmo Doolitte



In [25]:

    
print("Factorizacion LU ...")
A = Array([[1,3,5],[2,4,7],[1,1,0]])
print(A)
print("\n")
P,L,U = A.LU()
print("P:")
print(P)
print("L:")
print(L)
print("U:")
print(U)









    



Factorizacion LU ...
[[1 3 5 ]
[2 4 7 ]
[1 1 0 ]]


P:
[[0.0 1.0 0.0 ]
[1.0 0.0 0.0 ]
[0.0 0.0 1.0 ]]
L:
[[1.0 0.0 0.0 ]
[0.5 1.0 0.0 ]
[0.5 -1.0 1.0 ]]
U:
[[2.0 4.0 7.0 ]
[0.0 1.0 1.5 ]
[0.0 0.0 -2.0 ]]

Comprobación haciendo PA = LU



In [26]:

    
print("PA")
print(P*A)
print("LU")
print(L*U)









    



PA
[[2.0 4.0 7.0 ]
[1.0 3.0 5.0 ]
[1.0 1.0 0.0 ]]
LU
[[2.0 4.0 7.0 ]
[1.0 3.0 5.0 ]
[1.0 1.0 0.0 ]]

11 Función lu_linsolve, utilizando el resultado de 10.

se resuelve primero Ly = Pb

luego se resuelve

Ux = y



In [27]:

    
print("Resolver sistema de ecuaciones ")
A = Array([[-2,-1,1],[-3,1,1],[2,1,3]])
b = Array([[1],[4],[-1]])
print("A:")
print(A)
print("b:")
print(b)
y = Array.lu_linsolve(A,b)
print("Resultado: ")
print(y)









    



Resolver sistema de ecuaciones 
A:
[[-2 -1 1 ]
[-3 1 1 ]
[2 1 3 ]]
b:
[[1 ]
[4 ]
[-1 ]]
Resultado: 
[[-1.0 ]
[1.0 ]
[0.0 ]]



In [ ]: