using LinearAlgebra, PyPlot, Interact

A = [1 -1
     3  2
     1  4]

A'A # AᵀA

b = [1,2,3]
A'A # Aᵀb

x̂ = (A'A) \ (A'b) # the least-squares solution

p = A*x̂

norm(b - A*x̂)

x̂ = A \ b # = least-square solution when A is non-square!

year = [1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022]
ΔT = [0.28, -0.19, 0.11, -0.02, 0.13, 0.16, 0.15, 0.33, 0.51, 0.14, 0.53, 0.3, 0.22, 0.31, 0.32, 0.56, 0.17, 0.36, 0.43, 0.46, 0.36, 0.27, 0.56, 0.25, 0.34, 0.6, 0.51, 0.34, 0.47, 0.71, 0.72, 0.61, 0.65, 0.5, 0.92, 0.27, 0.6, 0.73, 0.46, 0.44, 0.62, 0.69, 0.83, 1.12, 0.98, 0.75, 0.94, 1.14, 0.78, 0.89];

using PyPlot
plot(year, ΔT, "r.-")
xlabel("year")
ylabel("ΔT (°C) vs. 1901–2000 baseline")
title("January global average temperature change")

b = ΔT
A = [year.^0 year.-1973]

x̂ = A \ b

# Alternatively, the normal equations solution:
(A'A) \ A'b

using PyPlot
plot(year, ΔT, "r.-")
xlabel("year")
ylabel("ΔT (°C) vs. 1901–2000 baseline")
title("January global average temperature change: Linear fit")

plot(year, A * x̂, "k-")
legend(["data", "fit: slope $(round(x̂[2], sigdigits=2)) °C/year"])

#=
let fig = figure(), d = range(0,3,length=100)
    @manipulate for d₁=0:0.1:1, d₂=1.1:0.1:2, R₁=1:0.1:2, R₂=1:0.1:2
        withfig(fig) do
            x = [1 d₁^2; 1 d₂^2] \ [R₁, R₂]
            plot([d₁, d₂], [R₁, R₂], "ro")
            plot(d, x[1] .+ x[2]*d.^2, "k-")
            xlim(0,3)
            ylim(0,3)
            title("fit \$x\$ = $x")
            xlabel(L"d")
            ylabel(L"R")
        end
    end
end
=#

let d = range(0,3,length=100), d₁=0.5, d₂=1.5, R₁=1, leg=[]
    for R₂=1:0.2:2
        x = [1 d₁^2; 1 d₂^2] \ [R₁, R₂]
        plot([d₁, d₂], [R₁, R₂], "ro", label="_nolegend_")
        plot(d, x[1] .+ x[2]*d.^2, "k-")
        xlim(0,3)
        ylim(0,3)
        push!(leg, L"$R_1 = %$R₁$, $x = %$(round.(x, digits=2))$")
        xlabel(L"d")
        ylabel(L"R")
    end
    title(L"fits through various $R_1$")
    legend(leg)
end

d = range(0,2,length=200) # 20 points from 1 to 4
b = 1 .+ 2*d.^2 + randn(200)*0.1 # measurements with Gaussian random noise
plot(d, b, "r.")
xlabel(L"d")
ylabel(L"R")
title("noisy \"measurements\" of \$R = 1 + d^2\$")

A = [d.^0 d.^2]

rank(A)

inv(A)

x̂ = A \ b

d = range(0,2,length=200) # 20 points from 1 to 4
b = 1 .+ 2*d.^2 + randn(200)*0.1 # measurements with Gaussian random noise
plot(d, b, "r.")
plot(d, x̂[1] .+ x̂[2] * d.^2, "k-")
xlabel(L"d")
ylabel(L"R")
title("Fit A\\b to noisy \"measurements\" of \$R = 1 + d^2\$")

norm(b - A*x̂) / norm(b)  # the length of the "residual" b - Ax̂ is not zero!

let a = range(0,1.5,length=50),
    afine = range(0,1.5,length=1000),
    b = 1 .+ 2a + 3a.^2 + 4a.^3 + randn(length(a)),
    fig = figure()
    #@manipulate for n=slider(1:40, value=2)
    for n=1:40
        display(
        withfig(fig) do
            plot(a, b, "r.")
            A = a .^ (0:n-1)'
            x̂ = A \ b
            plot(afine, (afine .^ (0:n-1)') * x̂, "k-")
            xlabel(L"a")
            ylabel(L"b")
            xlim(0,1.6)
            ylim(-5,30)
            title("noisy cubic: least-square fit of degree $(n-1)")
        end
        )
    end
end

let a = range(-1,1,length=50),
    afine = range(-1,1,length=1000),
    b = 1 ./ (1 .+ 25 * a.^2),
    fig = figure()
    # @manipulate for n=slider(1:50, value=3)
    for n=1:50
        display(
        withfig(fig) do
            plot(a, b, "r.")
            A = a .^ (0:n-1)'
            x̂ = A \ b
            plot(afine, (afine .^ (0:n-1)') * x̂, "k-")
            xlabel(L"a")
            ylabel(L"b")
            xlim(-1,1)
            ylim(0,1)
            title("smooth function: polynomial fit of degree $(n-1)")
        end
        )
    end
end

n = 50
a = @. cos(((1:n)-0.5) * pi/n) # = Chebyshev nodes of order n
afine = range(-1,1,length=1000)
b = 1 ./ (1 .+ 25 * a.^2)
plot(a, b, "r.")
A = a .^ (0:n-1)'
x̂ = A \ b
plot(afine, (afine .^ (0:n-1)') * x̂, "k-")
xlabel(L"a")
ylabel(L"b")
xlim(-1,1)
ylim(0,1)
title("smooth function: Chebyshev polynomial interpolation of degree $(n-1)")