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

# ## IRLS

using Plots
using LinearAlgebra: norm, opnorm, pinv
using Interact
using Random: seed!

m, n = 100, 50
seed!(0)
A = randn(m,n)
b = randn(m)
σ1 = opnorm(A); # largest singular value is the operator 2-norm of A

# +
xls = pinv(A)*b # SVD based solution
iters = 1000

@manipulate for μ_factor in range(0.0, 2.5, 51)

    μ = μ_factor / σ1^2

    x = zeros(size(A,2))
    ils_err = zeros(iters+1)

    ils_err[1] = norm(x - xls)
    for k = 1:iters
        x = x - μ * (A'*(A*x-b))
        ils_err[k+1] = norm(x - xls) / norm(xls)
    end
    scatter(0:iters, log10.(ils_err),
        xlabel="Iteration", xtick=[0,iters],
        ylabel="log10(NRMSD)", ylim=(-15,5),
        label = "μ_factor = $(round(μ_factor,digits=3))",
        title="Iterative LS convergence plot")

end

plot!() # for saving as html

# -