import math
import numpy as np
import Seals

sl = Seals.method()

class Algebra:

    def __init__(self, function):
        self.f = function
        self.integral = self.Integral(self.f)
        self.roots = self.Roots(self.f)
        self.edo = self.Edo(self.f)
    
    def d(self, x, e):
        return (self.f(x + e) - self.f(x - e))/(2*e)

    class Integral:

        def __init__(self,function):
            self.f = function
            self.simple = self.Simple(function)
            self.double = self.Double(function)
        
        class Simple:
            def __init__(self, function):
                self.f = function

            def riemann(self,a,b,n=None):
                
                if n is None:
                    n = 10**6

                delta = (b-a)/n

                psi = a
                theta = 0

                while((psi+delta) <= b):
                    
                    theta += (self.f(psi) + self.f(psi + delta))/2
                    psi += delta
                    
                integral = delta*theta

                return integral

            def simpson(self,a,b,n=None):

                if n is None:
                    n = 10**6

                def x(i):
                    return a + i*h

                h = (b-a)/n

                eta = 0
                theta = 0

                psi = 1
                kappa = 1

                while(psi <= (n/2)):

                    eta = eta + self.f(x(2*psi - 1))
                    psi = psi + 1

                while(kappa <= ((n/2)-1)):

                    theta = theta + self.f(x(2*kappa))
                    kappa = kappa + 1

                return (h/3)*( self.f(x(0)) + self.f(x(n)) + 4*eta + 2*theta)


        class Double:
            
            def __init__(self,function):
                self.f = function

            def riemann(self,a,b,c,d,n=None,m=None):
                
                if n is None:
                    n = 10**4

                if m is None:
                    m = n

                dx = (b-a)/n
                dy = (d-c)/m
                kappa = a
                psi = c
                theta = 0

                while((psi + dy) < d):

                    while((kappa + dx) < b):

                        theta = theta + self.f(kappa, psi)
                        kappa = kappa + dx

                    psi = psi + dy
                    kappa = a

                return theta*(dx)*(dy)

            def simpson(self,a,b,c,d,n=None,m=None):
                
                if n is None:
                    n = 10**4

                if m is None:
                    m = n

                dx = (b-a)/n
                dy = (d-c)/m
                
                def x(i):
    
                    x = a + i*dx
                    
                    return x

                def y(i):
                    
                    y = c + i*dy
                    
                    return y

                def g(i):
                    
                    sigma = 0
                    upsilon  = 0
                    
                    zeta = 1
                    csi = 1
                    
                    while(zeta <= (m/2)):
                        
                        sigma += self.f(x(i),y(2*zeta - 1))
                        zeta += 1
                        
                    while(csi <= ((m/2)-1)):
                        
                        upsilon  += self.f(x(i),y(2*csi))
                        csi += 1
                    
                    return (dy/3)*( self.f(x(i),y(0)) + self.f(x(i),y(m)) + 4*sigma + 2*upsilon )

                eta = 0
                theta = 0
                
                psi = 1
                kappa = 1
                
                while(psi <= (n/2)):
                    
                    eta += g(2*psi - 1)
                    psi += 1
                    
                while(kappa <= ((n/2)-1)):
                    
                    theta += g(2*kappa)
                    kappa += 1
                
                return (dx/3)*( g(0) + g(n) + 4*eta + 2*theta)

    class Roots:

        def __init__(self, function=None):
            if function is not None:
                self.f = function

        def bissec(self,a,b,e=None):

            if e is None:
                e = 10**(-6)

            fa = self.f(a)
            
            while abs(a-b) > e:
                
                c = (a+b)/2
                fc = self.f(c)
                
                if (fa*fc) < 0:
                    
                    b = c
                    
                else:
                    
                    a = c
                    fa = fc
                    
                c = (a+b)/2
                
                return c

        def d(self, x, e):
            return (self.f(x + e) - self.f(x - e))/(2*e)

        def newton(self,a,e=None):

            if e is None:
                e = 10**(-6)
            
            fa = self.f(a)
            da = self.d(a,e)
            b = a - fa/da
            
            
            while abs(a-b) > e:
                
                b = a
                a -= (fa/da)
                fa = self.f(a)
                da = self.d(a,e)
                
            return a

        def bissec_newton(self,a,b,e=None):

            if e is None:
                e = 10**(-6)
                
            fa = self.f(a)
            
            c = (a+b)/2                             # 'c' é a raiz calculada
            
            while abs(a-b) > 0.1:
                
                fc = self.f(c)
                
                if fa*fc < 0:
                    
                    b = c
                    
                else:
                    
                    a = c
                    fa = self.f(a)
                    
                c = (a+b)/2
                
            fc = self.f(c)
            dc = self.d(c,e)
            h = c - fc/dc                          # 'h' é uma variável de controle
            
            while abs(c-h) > e:
                
                h = c
                c -= (fc/dc)
                fc = self.f(c)
                dc = self.d(c,e)
                
            return (c)


    class Edo:

        def __init__(self, function):
            self.f = function

        def euler(self,a,y,b,n=None):

            if n is None:
                n = 10**7

            dx = (b-a)/n

            def x(i):
                return a + i*dx

            for i in range(n):

                y = y + (self.f(x(i),y))*dx
                
            return y

        def runge(self,a,y,b,n=None):

            if n is None:
                n = 10**7

            dx = (b-a)/n

            def x(i):
                return (a + i*dx)

            for i in range(n):

                y = y + (dx/2)*(self.f(x(i),y)+self.f(x(i+1),(y+(dx*self.f(x(i),y)))))
                
            return y

class Interpolation:
    """ Data should be organized in two columns: X and Y"""

    def __init__(self, data):

        self.data = data
        self.polinomial = self.Polinomial(self.data)

    def minimus(self,x):

        theta = 0
        # somatorio de x
        for i in range(self.data.shape[0]):

            theta += self.data[i][0]

        eta = 0
        #somatorio de y
        for i in range(self.data.shape[0]):

            eta += self.data[i][1]

        sigma = 0
        #somatorio de xy
        for i in range(self.data.shape[0]):

            sigma += self.data[i][0]*self.data[i][1]

        omega = 0
        #somatorio de x^2
        for i in range(self.data.shape[0]):

            omega += self.data[i][0]**2


        self.a = (self.data.shape[0]*sigma - theta*eta)/(self.data.shape[0]*omega - (theta**2))

        self.b = (theta*sigma - eta*omega)/((theta**2) - self.data.shape[0]*omega)
        
        ym = 0

        for i in range(self.data.shape[0]):

            ym += self.data[i][1]/self.data.shape[0]

        sqreq = 0

        for i in range(self.data.shape[0]):

            sqreq += ((self.a*self.data[i][0] + self.b) - ym)**2

        sqtot = 0

        for i in range(self.data.shape[0]):

            sqtot += (self.data[i][1] - ym)**2

        self.r2 = sqreq/sqtot

        return self.a*x + self.b

    class Polinomial:

        def __init__(self, data):
            self.data = data

        def vandermonde(self, x):

            matrix = np.zeros((self.data.shape[0],self.data.shape[0]))
            
            for k in range(0, self.data.shape[0]):

                matrix[:,k] = self.data[:,0]**k 

            self.A = sl.gauss(np.c_[matrix,self.data[:,1]])
            
            y = 0

            for i in range(0,self.A.shape[0]):

                y += self.A[i]*(x**i)

            return float(y)

        def lagrange(self, x):

            data_x = self.data[:,0]
            data_y = self.data[:,1]

            def L(k,x):
        
                up = down = 1

                for i in [x for x in range(data_x.shape[0]) if x != k]:
                    up = up*(x - data_x[i])

                for i in [x for x in range(data_x.shape[0]) if x != k]:
                    down = down*(data_x[k] - data_x[i])

                return up/down

            y = 0

            for i in range(data_x.shape[0]):

                y += data_y[i]*L(i,x)

            return y

        def newton(self,x):
            
            d = np.array(np.zeros((self.data.shape[0],self.data.shape[0])))

            d[0] = self.data[:,1]

            i = j = 0

            while (i < self.data.shape[0]):

                while (j < (self.data.shape[0]-(i+1))):

                    d[i+1][j] = (d[i][j+1] - d[i][j])/(self.data[(i+1)+j][0]-self.data[j][0])
                    j += 1
                
                i += 1
                j = 0
            
            def f(x):

                y = d[0][0]
                i = 0
                
                while ((i+1) < self.data.shape[0]):

                    mult = 1
                    k = 0
                    while (k <= i):
                        mult = mult*(x - self.data[k][0])
                        k += 1

                    y += d[i+1][0]*mult
                    i += 1

                return y

            self.f = f
            
            return f(x)

        def gregory(self,x):

            h = self.data[0][0] - self.data[1][0]
            
            d = np.array(np.zeros((self.data.shape[0],self.data.shape[0])))

            d[0] = self.data[:,1]

            i = j = 0

            while (i < self.data.shape[0]):

                while (j < (self.data.shape[0]-(i+1))):

                    d[i+1][j] = (d[i][j+1] - d[i][j])/((i+1)*h)
                    j += 1
                
                i += 1
                j = 0
            
            y = d[0][0]
            i = 0
            
            while ((i+1) < self.data.shape[0]):

                mult = 1
                k = 0
                while (k <= i):
                    mult = mult*(x - self.data[k][0])
                    k += 1

                y += d[i+1][0]*mult
                i += 1

            return y