# -*- coding: utf-8 -*-
"""
Created on Wed May 12 17:56:25 2021

@author: Gabri and Vinícius
"""

import numpy as np
from scipy.linalg import lu_factor, lu_solve
from scipy.interpolate import interp1d


# =============================================================================
# Projeto Métodos Numéricos: Calculadora de métodos
# =============================================================================


class raizes:

    def __init__(self, f, n=0, Emax=0.001, maxit=10):
        """
        Encontra a raiz de um polinomio.
        :param n: Palpite inicial.
        :param f: Função a ser usada.
        :param f_linha: Derivada da função a ser usada.
        :param sigma: Função atualizadora de x para o metodo do ponto fixo.
        :param Emax: Erro maximo.
        :param maxit: Numero maximo de iterações.
        :param a, b: Limites inferior e superior de um intervalo no qual deverá conter a raiz.
        :param x0, x1: Chutes iniciais para o metodo da secante.
        :param d: Valor de pertubação da raiz, caso não especificado será 0.0001.
        :param x: Raiz da função.
        """
        self.f = f
        self.Emax = Emax
        self.maxit = maxit
        self.n = n
        self.x = None
        self.erro = None
        self.iter = None

    def bissecao(self, a, b):

        xrold = 0
        cont = 0

        while True:
            self.x = (a + b) / 2
            cont += 1

            if self.f(a) * self.f(self.x) < 0:
                b = self.x
                self.erro = abs(self.x - xrold)
            elif self.f(b) * self.f(self.x) < 0:
                a = self.x
                self.erro = abs(self.x - xrold)
            else:
                self.erro = 0

            xrold = self.x
            self.iter = cont

            if self.erro <= self.Emax or cont >= self.maxit:
                break

    def posicao_falsa(self, a, b):

        xrold = 0
        cont = 0

        while True:
            self.x = ((a * self.f(b)) - (b * self.f(a))) / (self.f(b) - self.f(a))
            cont += 1

            if self.f(a) * self.f(self.x) < 0:
                b = self.x
                self.erro = abs(self.x - xrold)
            elif self.f(b) * self.f(self.x) < 0:
                a = self.x
                self.erro = abs(self.x - xrold)
            else:
                self.erro = 0

            xrold = self.x
            self.iter = cont

            if self.erro <= self.Emax or cont >= self.maxit:
                break

    def ponto_fixo(self, sigma):

        cont = 0
        xrold = 0
        self.x = self.f(self.n)

        while True:
            cont += 1
            self.iter = cont
            if abs(self.x - xrold) < self.Emax or cont >= self.maxit:
                self.erro = abs(self.x - xrold)
                break
            xrold = self.x
            self.x = sigma(self.x)

    def secante(self, x0, x1):

        cont = 0
        xrold1 = x1
        xrold2 = x0

        while True:
            self.x = xrold1 - (self.f(xrold1) * (xrold2 - xrold1)) / (self.f(xrold2) - self.f(xrold1))
            cont += 1

            self.erro = abs((self.x - xrold1) / self.x)

            xrold2 = xrold1
            xrold1 = self.x
            self.iter = cont

            if self.erro <= self.Emax or cont >= self.maxit:
                break

    def secante_mod(self, d=0.001):

        cont = 0
        xrold = self.n

        while True:
            self.x = xrold - (d * xrold * self.f(xrold)) / (self.f(xrold + d * xrold) - self.f(xrold))
            cont += 1

            self.erro = abs((self.x - xrold) / self.x)

            xrold = self.x
            self.iter = cont

            if self.erro <= self.Emax or cont >= self.maxit:
                break

    def Newton_Raphson(self, f_linha):
        contador = 0
        self.x = [self.n]  # x(n), sendo inicialmente nosso palpite
        self.iter = contador
        while True:
            contador += 1
            # a cada iteração vamos adicionar um novo valor de x à lista
            # x[-1] é o ultimo elemento de uma lista
            if f_linha(self.x[-1]) == 0:
                print(f'Não foi possivel encontrar a raiz pois f_linha(x) = 0 \nEscolher outro palpite')
                return 0
            self.x.append(self.x[-1] - self.f(self.x[-1]) / f_linha(self.x[-1]))
            # x(n+1) pela formula de Newton-Raphson
            self.erro = abs(self.x[-1] - self.x[-2])  # erro absoluto
            if self.erro < self.Emax:
                break
            if self.iter > self.maxit:
                break
        self.x = self.x[-1]


class interpol:

    def __init__(self, x, y, po, pf):
        """
        Resolve problemas de interpolação.
        :param x: Vetor com os valores de X.
        :param y: Vetor com os valores de Y.
        :param po: O valor inicial de X onde se pretende calcular a aproximação.
        :param pf: O valor final de X onde se pretende calcular a aproximação.
        :param tipo: O tipo de spline a ser utilizado, caso não especificado será spline cubico.
        :param grau: O grau do polinomio MMQ gerado, caso não especificado será grau 1.
        :param resultado: Dicionario com os valores entre os ponto pedidos e suas respectivas aproximações de F(x).
        """
        self.x = x
        self.y = y
        self.po = po
        self.pf = pf
        self.resultado = {None}

    def Newton(self):
        tam = len(self.x)
        res = np.array([], dtype='float32')
        pontos = np.linspace(self.po, self.pf, 100)  # Intervalo a ser calculado
        for p in pontos:
            FDD = np.ones(tam ** 2).reshape(
                (tam, tam)) * -1  # Tabela de diferenças divididas(começa em -1 para visualização)
            for i in range(tam):
                FDD[i][0] = self.y[i]  # Primeira coluna
            for i in range(1, tam):
                for j in range(tam - i):  # Demais colunas
                    FDD[j][i] = (FDD[j + 1][i - 1] - FDD[j][i - 1]) / (
                            self.x[j + i] - self.x[j])  # Faz todas as colunas antes de iterar a linha
            x_termo = 1  # Auxiliar para a multiplicaçãos pelos termos de X
            pol = FDD[0][0]  # Polinomio de Newton
            for ordem in range(1, tam):
                x_termo = x_termo * (p - self.x[ordem - 1])
                pol += FDD[0][ordem] * x_termo  # Operador(primeiro termo de cada coluna)*multiplicação dos x
            res = np.append(res, pol)
        self.resultado = {"X": pontos,
                          "Y": res}

    def Lagrange(self):
        res = np.array([], dtype='float32')
        r = 0
        tam = len(self.x)
        pontos = np.linspace(self.po, self.pf, 100)  # Intervalo a ser calculado
        for ponto in pontos:
            for i in range(0, tam):  # Cada ordem de y(i)L(i)
                produtorio = 1  # L(i)
                for j in range(0, tam):
                    if j != i:
                        produtorio *= (ponto - self.x[j]) / (self.x[i] - self.x[j])
                r += produtorio * self.y[i]
            res = np.append(res, r)
            r = 0
        self.resultado = {"X": pontos,
                          "Y": res}

    def Spline(self, tipo=3):
        if tipo == 1:
            s = interp1d(self.x, self.y, kind='linear')
            pontos = np.linspace(self.po, self.pf, 100)  # Intervalo a ser calculado
            res = np.array([s(i) for i in pontos])
            self.resultado = {"X": pontos,
                              "Y": res}
        elif tipo == 2:
            s = interp1d(self.x, self.y, kind='quadratic')
            pontos = np.linspace(self.po, self.pf, 100)  # Intervalo a ser calculado
            res = np.array([s(i) for i in pontos])
            self.resultado = {"X": pontos,
                              "Y": res}
        else:
            s = interp1d(self.x, self.y, kind='cubic')
            pontos = np.linspace(self.po, self.pf, 100)  # Intervalo a ser calculado
            res = np.array([s(i) for i in pontos])
            self.resultado = {"X": pontos,
                              "Y": res}

    def MMQ(self, grau=1):

        def mmq(x, th, grau):  # Devolve a soma dos termos do polinomio
            aux = grau
            soma = 0
            for c in range(0, grau):
                soma += th[c] * (x ** aux)
                aux -= 1
            soma += th[-1]
            return soma

        M_x = np.ones(len(self.x) * (grau + 1)).reshape(len(self.x), grau + 1)  # Cria a matrix a ser usada
        aux = grau
        for c in range(0, grau):
            for i in range(0, len(self.x)):
                M_x[i][c] = self.x[i] ** aux
            aux -= 1
        th = np.linalg.inv(M_x.transpose().dot(M_x)).dot(M_x.transpose().dot(self.y))  # Aplica a formula do MMQ
        pontos = np.linspace(self.po, self.pf, 100)  # Intervalo a ser calculado
        res = np.array([mmq(i, th, grau) for i in pontos])
        self.resultado = {"X": pontos,
                          "Y": res}


class Sistemas_Lineares():
    def __init__(self, A, P, max_it=10):
        """
        Resolve problemas de sistemas lineares do tipo Ax=P.
        :param A: Matriz com os coeficientes de x.
        :param P: Matriz com os valores independentes.
        :param max_it: Numero maximo de iterações.
        :param resultado: Vetor com as soluções.
        """
        self.A = A
        self.P = P
        self.max_it = max_it
        self.resultado = None

    def LU(self):
        lu, piv = lu_factor(self.A)
        self.resultado = lu_solve((lu, piv), self.P)

    def Gauss_Siedel(self):
        self.resultado = np.copy(self.P) * 0
        for c in range(0, self.max_it):  # Repete quantas vezes forem passadas
            for i in range(len(self.resultado)):
                s = self.P[i].copy()  # s = posição de P
                for j in range(len(self.resultado)):
                    if i != j:
                        s -= self.A[i, j] * self.resultado[j]  # Subtraimos tds A*x exceto onde no x[i]
                s /= self.A[i, i]  # Dividimos pelo A em questão
                self.resultado[i] = s  # Atualizamos vetor resultado


class EDO():
    def __init__(self, f_linha, xi, fi, h=0.1):
        """
        Encontra a aproximação da EDO por um metodo numerico.
        :param f_linha: Derivada isolada da função.
        :param xi: Posição da condição conhecida.
        :param fi: Condição conhecida.
        :param h: Valor do passo, caso não especificado será 0.1.
        :param x_anterior: Variavel auxiliar para aplicação do metodo.
        :param y_anterior: Valor passado da aproximação da solução da EDO no ponto t.
        :param y_novo: Valor da aproximação da solução da EDO no ponto t.
        :param t: Posição do valor no qual deseja ser encontrado.
        """
        self.f_linha = f_linha
        self.xi = xi
        self.fi = fi
        self.h = h
        self.x_anterior = None
        self.y_anterior = None
        self.y_novo = None

    def Euler(self, t):
        self.t = t
        if self.xi == self.t:  # Retorna valor conhecido caso busquemos sua posição.
            self.y_novo = self.fi

        n = int((self.t - self.xi) / self.h)
        self.y_anterior = self.fi
        self.x_anterior = self.xi
        for p in range(0, n):  # Aplica a formula n vezes.
            self.y_novo = self.y_anterior + self.h * self.f_linha(self.x_anterior, self.y_anterior)
            self.x_anterior += self.h  # Atualiza o valor de x.
            self.y_anterior = self.y_novo  # Atualiza o valor de y.

    def Heun(self, t):
        self.t = t
        if self.xi == self.t:  # Retorna valor conhecido caso busquemos sua posição.
            self.y_novo = self.fi

        n = int((self.t - self.xi) / self.h)
        for p in range(0, n):  # Aplica a formula n vezes.
            if p == 0:
                self.y_anterior = self.fi
                self.x_anterior = self.xi
            self.y_novo = self.y_anterior + (self.h / 2) * (self.f_linha(self.x_anterior, self.y_anterior) +
                                                            self.f_linha(self.x_anterior + self.h,
                                                                         self.y_anterior + self.h * self.f_linha(
                                                                             self.x_anterior,
                                                                             self.y_anterior)))
            self.x_anterior += self.h  # Atualiza o valor de x.
            self.y_anterior = self.y_novo  # Atualiza o valor de y.

    def Runge_Kutta_4(self, t):
        self.t = t
        if self.xi == self.t:  # Retorna valor conhecido caso busquemos sua posição.
            self.y_novo = self.fi

        n = int((self.t - self.xi) / self.h)
        for p in range(0, n):  # Aplica a formula n vezes.
            if p == 0:
                self.y_anterior = self.fi
                self.x_anterior = self.xi
            k1 = self.f_linha(self.x_anterior, self.y_anterior)
            k2 = self.f_linha(self.x_anterior + 0.5 * self.h, self.y_anterior + k1 * 0.5 * self.h)
            k3 = self.f_linha(self.x_anterior + 0.5 * self.h, self.y_anterior + k2 * 0.5 * self.h)
            k4 = self.f_linha(self.x_anterior + self.h, self.y_anterior + k3 * self.h)
            self.y_novo = self.y_anterior + (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4) * self.h
            self.x_anterior += self.h  # Atualiza o valor de x.
            self.y_anterior = self.y_novo  # Atualiza o valor de y.


class Integral():
    def __init__(self, f, a, b, n=1):
        """
        Encontra a aproximação da integral por um metodo numerico.
        :param f: Função a ser integrada.
        :param a: Limite de integração inferior.
        :param b: Limite de integração superior.
        :param n: Quantidade de subdivições do metodo, será 1 caso não especificado.
        :param resultado: Valor da aproximação.
        """
        self.f = f
        self.a = a
        self.b = b
        self.n = n
        self.resultado = None

    def Trapezio(self):
        h = (self.b - self.a) / self.n
        x = np.arange(self.a, self.b + h, h)
        y = np.array([self.f(i) for i in x])
        soma = y.sum()
        self.resultado = (h / 2) * ((2 * soma) - y[0] - y[-1])

    def Simpson(self):
        h = (self.b - self.a) / (2 * self.n)
        x = np.arange(self.a, self.b + h, h)
        y = np.array([self.f(i) for i in x])
        for i in range(len(y)-1):
            if i == 0:
                par = 0
                impar = 0
            elif i % 2 == 0:
                par += y[i]
            else:
                impar += y[i]
        self.resultado = (h / 3) * (y[0] + y[-1] + 2 * par + 4 * impar)
