# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

# # Aproximação e Interpolação
#
# Como aproximar um conjunto de dados usando uma sequência de funções simples?
#
#  * Polinômios
#  * Interpolação - Introdução
#  * Interpolação de Lagrange
#  * Método dos mínimos quadrados
#  * Códigos

# # Introdução
#
# Qual o objetivo?
#
#  1. Modelar dados experimentais
#  2. Aproximar uma função complexa usando funções mais simples
#
# Dada uma função $u(x)$ queremos aproximá-la usando um conjunto de funções que são, de alguma, forma convenientes:
#
# $$
# u \approx u^\delta (x) = \sum_{i=1}^N \hat{u}_i \phi_i(x)
# $$
#
#
# Porque queremos fazer isso? 
#
#  * Funções muito complexas 
#  * Funções que não são conhecidas explicitamente (algumas funções especiais)
#  * Queremos desempenho na hora de calcular as funções: praticamente todas as funções básicas são aproximadas utilizando polinômios (muitas vezes polinômios racionais)
#  * Quando conhecemos apenas alguns pontos da função: sempre que você tiver um experimento ou estiver processando dados de uma sequência de simulações numéricas.
#  
#  
#  
#

methods(sin)

# ## Interpolação
#
# É lógico que, em geral, a aproximação é apenas uma aproximação! Então existe um erro:
#
# $$
# u(x) - u^\delta(x) = \varepsilon(x)
# $$
#
# A idéia é minimizar este erro de alguma forma.
#
# Uma possibilidade é escolher $N$ pontos $x_i$ e impor que o erro é zero nestes pontos:
#
# $$
# \varepsilon(x_i) = 0 \qquad i=1, \ldots, N
# $$
#
# Esta abordagem é conhecida na literatura como colocação. Se a função $u(x)$ é conhecida apenas por alguns pontos $(x_i, u_i)$ isto é chamado de *interpolação*.
#
# Isto resulta num sistema de equações lineares:
#
# $$
# \left(\begin{matrix}
# \phi_1(x_1) & \phi_2(x_1) & \cdots & \phi_N(x_1) \\
# \phi_1(x_2) & \phi_2(x_2) & \cdots & \phi_N(x_2) \\
# \vdots & \vdots & \ddots & \vdots \\
# \phi_1(x_N) & \phi_2(x_N) & \cdots & \phi_N(x_N) \\
# \end{matrix}\right)
# \cdot\left(\begin{matrix} \hat{u}_1 \\ \hat{u}_2 \\ \vdots \\ \hat{u}_N \end{matrix}
# \right)
# = \left(\begin{matrix} u(x_1) \\ u(x_2) \\ \vdots \\ u(x_N) \end{matrix}
# \right)
# $$
# Esta é a matriz de Vandermonde.
#
# A escolha de funções $\phi_i(x)$ adequadas simplificam a solução do problema. Se, por exemplo, $\phi_i(x_k) = \delta_{ik}$ onde $\delta$, neste caso, é o delta de Kronecker, a matriz é diagonal.
#

# ## Mínimos quadrados
# Por outro lado, podemos, para uma sequência de pontos minimizar o erro quadrático total:
#
# $$
# R(\hat{u}_1, \ldots, \hat{u}_N) = \sum_{i=1}^Q \left[ u(x_i) - u^\delta(x_i)\right]^2
# $$
# onde $Q$ é o número de pontos. Esta operação se chama método dos mínimos quadrados.
#
# Chamando $u(x_i) = u_i$, para minimizar este resíduo (erro quadrático total), basta derivar e igualar a zero, assim, chega-se ao seguinte sistema de equações lineares:
#
# $$
# \left(
# \begin{matrix}
# \sum_{i=1}^Q \phi_1(x_i)\cdot\phi_1(x_i) & 
# \cdots &
# \sum_{i=1}^Q \phi_1(x_i)\cdot\phi_N(x_i) \\
# \sum_{i=1}^Q \phi_2(x_i)\cdot\phi_1(x_i) & 
# \cdots &
# \sum_{i=1}^Q \phi_2(x_i)\cdot\phi_N(x_i) \\
# \vdots & \ddots & \vdots \\
# \sum_{i=1}^Q \phi_N(x_i)\cdot\phi_1(x_i) & 
# \cdots &
# \sum_{i=1}^Q \phi_N(x_i)\cdot\phi_N(x_i) \\
# \end{matrix}\right)
# \cdot
# \left(\begin{matrix} \hat{u}_1 \\ \hat{u}_2 \\ \vdots \\ \hat{u}_N \end{matrix}\right)
# = 
# \left(\begin{matrix} \sum_{i=1}^Q u_i \phi_1(x_i) \\ \sum_{i=1}^Q u_i \phi_2(x_i) \\  \vdots \\ \sum_{i=1}^Q u_i \phi_N(x_i)\end{matrix}\right)
# $$
#
#
#
#
#

# ## Uma formulação mais geral
#
# Uma formulação mais geral é escolher funções peso $w_k(x)$ de modo que
# $$
# \int_a^b u(x)w_k(x)\:dx = \int_a^b u^\delta (x) w_k(x) \:dx
# $$
#
# Escolhendo $w_k(x)$ = $\delta(x_k)$ (*Delta* de Dirac), recuperamos a interpolação. Escolhendo $w_k(x) = \phi_k(x)$ temos o método de Galerkin, muito usado no método de elementos finitos.
#
# Com o método de elementos finitos, a função não necessariamente vai passar pelos pontos, mas o erro será minimizado e pode-se obter uma aproximação melhor.
#
# Outro ponto a ser considerador é o que ocorre quando as medidas possuem erro e como isso afeta a aproximação.

# # Polinômios
#
# $$
# p_n(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \ldots a_n x^n = \sum_{i=0}^n a_i x^i
# $$
# Um polinômio é dado por seus coeficientes $a_i$.
# Em Julia, a indexação começa em 1, portanto é conveniente reescrever a definição acima como:
#
# $$
# p_n(x) = a_1 + a_2 x + a_3 x^2 + a_4 x^3 + \ldots a_n x^{n-1} + a_{n+1} x^n = \sum_{i=1}^{n+1} a_i x^{i-1}
# $$
#
#

function polyval1(a, x)
    n = length(a)
    
    y = 0.0
    for i in 1:n
        y += a[i] * x^(i-1)
    end
    
    return y
end


# ### Método de Horner
#
# Maneira mais eficiente e mais precisa de calcular o valor de um polinômio:
#
# $$
# p_n(x) = a_1 + x\left(a_2 + x\left(a_3 + x\left(\ldots + a_{n+1} x\right) \right) \right)
# $$

# +
function polyval2(a, x)
    
    y = a[end]
    
    for i = (lastindex(a)-1):-1:1
        y = a[i] + x * y
    end
    
    return y
    
end

# -

a1 = rand(3)
a2 = rand(6)
a3 = rand(15)
a4 = rand(30);


using PyPlot

using BenchmarkTools

x = 0.5
y1 = polyval1(a3, x)
y2 = polyval2(a3, x)
y1≈y2

@btime polyval1(a1, x)
@btime polyval2(a1, x)

@btime polyval1(a2, x)
@btime polyval2(a2, x)

@btime polyval1(a3, x)
@btime polyval2(a3, x)

typeof(a3)
@btime polyval1(a4, x)
@btime polyval2(a4, x)

# ### E os tipos dos dados???

a5 = rand(-5:5, 7)

polyval1(a5, 2)

polyval1(a5, 2.0)

polyval2(a5, 2)

polyval2(a5, 2.0)

a6 = [1//2, 2//3, 3//4, 4//5]

polyval2(a6, 0.5)

polyval2(a6, 1//2)

polyval2(a6, 1)

@code_warntype polyval2(a6, 0.5)

# +
function polyval3(a::Vector{T}, x::S) where {T,S}
    R = promote_type(T,S)
    
    y = convert(R, a[end])
    
    for i = (lastindex(a)-1):-1:1
        y = a[i] + x * y
    end
    
    return y
    
end

# -

@code_warntype polyval3(a6, 0.5)

@btime polyval2(a6, 0.5+0.2im)

@btime polyval3(a6, 0.5+0.2im)

# ### Polynomials.jl

using Polynomials

p = Polynomial([1,2,3], :z)

typeof(p)

p(1//2 + 3//4im)

roots(p)

p2 = Polynomial([1,2]) # Raízes

roots(p2)

integrate(p)

derivative(p)

degree(p)

# ### Macro quando os coeficientes são conhecidos

@evalpoly(0.5, 1, 2, 3)


@macroexpand @evalpoly(0.5, 1, 2, 3)

evalpoly(0.5, [1,2,3])


# # Interpolação
#
# Dado um conjunto de n pontos $(x_i, y_i)$, qual o poliômio que passa por todos?
#
# $$
# y_i = a_0 + a_1 x_i + a_2 x_i^2 + \ldots a_n x_i^n \qquad i=1, \ldots, m
# $$

# ## Vandermonde
#
# Com n+1 pontos distintos, se o polinômio for de grau n, pode-se montar o seguinte sistema linear:
#
# $$
# \begin{bmatrix}
# 1 & x_0 & x_0^2 & \cdots & x_0^n \\
# 1 & x_1 & x_1^2 & \cdots & x_1^n \\
# \vdots & \vdots & \vdots & \ddots & \vdots\\
# 1 & x_n & x_n^2 & \cdots & x_n^n\\
# \end{bmatrix}\cdot
# \left\{ \begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_n\\ \end{matrix}\right\}
# = \left\{\begin{matrix} y_0 \\ y_1 \\ \vdots \\ y_n\\\end{matrix}\right\}
# $$
#
# Mas esta é uma operação cara, $\mathcal{O}(n^3)$!

# ## Outra possibilidade
#
# $$
# y = f(x) \approx a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + \cdots + a_n(x-x_0)(x-x_1)\cdots(x - x_{n-1})
# $$
#
# Com isso chegamos ao seguinte sistema linear triangular:
# $$
# \begin{bmatrix}
# 1 & 0 & 0 &  0 &\cdots & 0\\
# 1 & (x_1-x_0) & 0 & 0 &  \cdots & 0\\
# 1 & (x_2 - x_0) &  (x_2 - x_0)(x_2 - x_1) & 0 & \cdots & 0\\
# \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\
# 1 & (x_n - x_0) &  (x_n - x_0)(x_n - x_1) & (x_n - x_0)(x_n - x_1)(x_n-x_2) & \cdots &(x_n-x_0)(x_n-x_1)\cdots(x_n-x_{n-1})\\
# \end{bmatrix}\cdot\left\{\begin{matrix} a_0\\ a_1 \\ a_2 \\ \vdots \\ a_n\end{matrix}\right\} = 
# \left\{\begin{matrix} y_0\\ y_1 \\ y_2 \\ \vdots \\ y_n\end{matrix}\right\}
# $$
#
# Resolver este sistema é muito mais barato: $\mathcal{O}(n^2)$

# ## Interpolação de Lagrange:
#
# $$
# y(x) = \sum_{i=1}^n y_i h_i(x)
# $$
#
# onde $h_i(x)$ é o interpolador de Lagrange:
#
# $$
# h_k(x) = \prod_{i=1\ldots n,}^n \frac{x - x_i}{x_k - x_i} \qquad i\ne k
# $$
#
# Propriedade:
# $$
# h_i(x_j) = \delta_{ij} \quad \text{onde} \quad \delta_{ij} = \left\{\begin{matrix}1, \: i=j \\ 0, i\ne j\\ \end{matrix}\right.
# $$

# +
function lagrange1(k, z, x)
    
    h = 1.0
    
    n = length(z)
    
    for i = 1:n
        if i != k
            h *= (x - z[i]) / (z[k] - z[i])
        end
    end
    
    return h
end
    
    
# -

function lagrange2(k, z, x)
    
    h = 1.0
    n = length(z)
    for i = 1:(k-1)
        h *= (x - z[i]) / (z[k] - z[i])
    end
    
    for i = (k+1):n
        h *= (x - z[i]) / (z[k] - z[i])
    end
    return h
end

# +
N = 10
x = range(-1.0, 1.0, step=0.2)
#x = [cos(k*π/N) for k in 0:N]
xx = range(-1.0, 1.0, step=0.005)



# -

hh = [lagrange2.(k, Ref(x), xx) for k in 1:length(x)];

for i in 1:length(x)
    plot(xx, hh[i])
end


# +
plot(xx, hh[5])
axhline(y=0, color="black", linestyle = "--")
axhline(y=1, color="black", linestyle = "--")

for xv in x
    axvline(x=xv, color="red", linestyle="--")
end

# -

# ### Vamos organizar a interpolação de Lagrange

# +
struct Lagrange
    x::Vector{Float64}
    y::Vector{Float64}
    Lagrange(x, y) = new(copy(x), copy(y))
end
Base.Broadcast.broadcastable(lgr::Lagrange) = Ref(lgr)

function lagrange(k, z, x)
 h = 1.0
    n = length(z)
    for i = 1:(k-1)
        h *= (x - z[i]) / (z[k] - z[i])
    end
    
    for i = (k+1):n
        h *= (x - z[i]) / (z[k] - z[i])
    end
    return h
end

function interp(lgr::Lagrange, x)
    
    y = lgr.y[1] * lagrange(1, lgr.x, x)
    
    for i = 2:length(lgr.x)
        y += lgr.y[i] * lagrange(i, lgr.x, x)
    end
    
    return y
end

    
(lgr::Lagrange)(x) = interp(lgr, x)    
# -

x1 = range(-1, 1, step=0.2)
y1 = sin.(π.*x1);
xx = range(-1, 1, step=0.005)


lgr = Lagrange(x1, y1)
yy = sin.(π.*xx);
yy1 = lgr.(xx);

# +

plot(x1, y1, "ro")
plot(xx, yy, "r-")
plot(xx, yy1, "b--")
# -

# ## Interpolação Linear
#
#

# +
x2 = range(-1, 1, step=0.2)
y2 = sin.(π.*x2);

plot(x2, y2, "ro")
plot(x2, y2, "b-")

# +
struct LinearInterp
    x::Vector{Float64}
    y::Vector{Float64}
    LinearInterp(x, y) = new(copy(x), copy(y))
end
Base.Broadcast.broadcastable(lin::LinearInterp) = Ref(lin)

function interp1(lin::LinearInterp, x)
    
    if x < lin.x[1] || x > lin.x[end]
        error("Fora do Range")
    end
    
    index = 2
    n = length(lin.x)
    for i = 2:n
        if lin.x[i] >= x
            index = i
            break
        end
    end
    i1 = index-1
    return lin.y[i1] + (lin.y[index] - lin.y[i1]) * (x - lin.x[i1]) / (lin.x[index] - lin.x[i1])
    
end

(lin::LinearInterp)(x) = interp1(lin, x)

# -

lin = LinearInterp(x2, y2)
yy2 = interp1.(lin, xx);
plot(x2, y2, "ro")
plot(xx, yy2, "b-")

@btime yy2 .= interp1.(lin, xx);

# # Mínimos quadrados
#

function linfit(x,y)
    sx = sum(x)
    sx2 = sum(x->x^2, x)
    N = length(x)

    sy = sum(y)
    syx = sum(x[i]*y[i] for i in 1:N)
    
    return [N sx; sx sx2] \ [sy; syx]
end

x = 1.0:20
y = 2 .*x .+ 3.0

linfit(x,y)

# # Pacotes
#
#  * [Interpolations](https://github.com/JuliaMath/Interpolations.jl)
#  * [Dierckx](https://github.com/kbarbary/Dierckx.jl)
#  * [GridInterpolations](https://github.com/sisl/GridInterpolations.jl)

# ## Exercícios
#
# ### Problema 1
#
# Interpole a função de Runge com $-1 \le x \le 1$:
# $$
# f(x) = \frac{1}{1 + 25x^2}
# $$
#
#  1. Use 11 pontos uniformemente distribuídos
#  2. Aumente o número de pontos
#  3. Tente usar os pontos $x_k = \cos\left(\frac{k\pi}{N}\right)$ para $k = 0\ldots N$.
#  4. Brinque com o número de pontos
#
# ### Problema 2
#
# Procure na Net o método de diferenças divididas de Newton a interpole a função anterior nos mesmos pontos. Este método é simplesmente um jeito inteligente de resolver a matriz apresentada lá em cima.
#
# ### Problema 3
#
# Use a biblioteca Interpolations.jl e Dierckx.jl para fazer as interpolações. Compare a interpolação linear com os splines.
#
# ### Problema 4
#
# Crie funções para fazer os seguintes problemas de mínimos quadrados:
#  * $y = a_0 x^ a_1$
#  * $y = a_0 \exp \left( a_1 \cdot x\right)$
#  * Polinômio genérico de ordem n
#

f(x) = 1.0 / (1.0 + 25x^2)


xx = -1:0.01:1
yy = f.(xx);

plot(xx, yy)



StarcoderdataJupyterScriptsDedupFiltered