text:
"""
    Epd(μ, σ, p)

The *Exponential power distribution* with location `μ`, scale parameter `σ` and shape
parameter `p` has probability density function

```math
f(x; \\mu, \\sigma, p) = \\frac{1}{\\sigma} \\frac{1}{2 p^{1/p}\\Gamma(1+1/p)} \\exp\\left\\{-\\frac{1}{p}\\Big| \\frac{x-\\mu}{\\sigma} \\Big| \\right\\}
```

for

```
-\\infty < x < \\infty
```

```julia
Epd(μ, σ, p)      # Epd with shape p, scale σ and location μ.
params(d)       # Get the parameters, i.e. (μ, σ, p)
location(d)     # Get the location parameter, i.e. μ
scale(d)        # Get the scale parameter, i.e. σ
shape(d)        # Get the shape parameter, i.e. p (sometimes called θ)
```

External links
* [Epd](https://www.sciencedirect.com/science/article/pii/S0304407608001668?casa_token=kEUNFlIXYBEAAAAA:aCZkSeMVO4y3H7J9cCyR18j5R6QifeKYa_PSkvgDaCBO-xYQKxbeB1YQLpJbTQowZR_4fIKYcOc)
"""
struct Epd{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    p::T
    Epd{T}(µ::T, σ::T, p::T) where {T} = new{T}(µ, σ, p)
end

function Epd(µ::T, σ::T, p::T; check_args=true) where {T <: Real}
    check_args && @check_args(Epd, σ > zero(σ))
    check_args && @check_args(Epd, p > zero(p))
    return Epd{T}(µ, σ, p)
end

Epd(μ::Real, σ::Real, p::Real) = Epd(promote(μ, σ, p)...)
Epd(μ::Integer, σ::Integer, p::Integer) = Laplace(float(μ), float(σ), float(p))
Epd(μ::T, p::T) where {T <: Real} = Epd(μ, one(T), p)
Epd(μ::T) where {T <: Real} = Epd(μ, one(T), 2.)
Epd() = Laplace(0.0, 1.0, 2., check_args=false)

@distr_support Epd -Inf Inf

location(d::Epd) = d.μ
scale(d::Epd) = d.σ
shape(d::Epd) = d.p
params(d::Epd) = (d.μ, d.σ, d.p)
@inline partype(d::Epd{T}) where {T<:Real} = T

#### Statistics

#### Evaluations
function pdf(d::Epd, x::Real)
    μ, σ, p = params(d)
    K = σ * 2 * p^(1/p) * gamma(1 + 1/p)
    exp(-1/p * abs((x - μ)/σ)^p)/K
end

logpdf(d::Epd, x::Real) = log(pdf(d, x))

quantile(d::Epd, x::Real) = sign(2*x-1) * (d.p * quantile(Gamma(1/d.p, d.σ^d.p), abs(2*x-1)) )^(1/d.p)

#### Sampling
rand(rng::AbstractRNG, d::Epd) = quantile(d, rand(rng)) + d.μ


text:
@static if VERSION < v"0.7.0-DEV.2005"
  using Base.Test
else
  using Test
end

@testset "RunTests" begin

include("runsymbolictests.jl")
#include("TestModelsWithError.jl")
include("runsimulationtests.jl")
include("../examples/runexamples.jl")
#include("runlargetests.jl")

end # testset
nothing

ext: jl


text:
"""
Parsing the JSON in `data`
"""
function extract_country_data(country, data)
    deaths = map(e -> e["deaths"], data[country])
    confirmed=  map(e -> e["confirmed"], data[country])
    dates=  map(e -> Date(e["date"],dateformat"y-m-d"), data[country])
    recovered=map(e -> e["recovered"], data[country])
    (deaths=deaths, confirmed=confirmed, dates=dates, recovered=recovered)
end

function days_till_target_number(target_number,d=deaths, v=v)
    c=d[end] # Number of infected today
    x=log(target_number/c)/log(2)
    doubling_time=log(2)/v[2]
    time=ceil(Int,x*doubling_time)
    doubling_time, time
end

function segmented_models(d,model_interval=5)
    segments=[i:i+model_interval-1 for i=1:model_interval:length(d)-model_interval+1]
    confirmed_segments=map(s->d[s], segments)
    v_interval=[[ones(model_interval) [segments[i]...]]\ log.(confirmed_segments[i]) for i=1:length(confirmed_segments)]
    (segments=segments, v_interval=v_interval)
end

function plot_doubling_times(d,v_interval,segments,dates, p_head="Confirmed cases"; st=:bar)
    d_times=[days_till_target_number(20_000,d,v)[1] for v in v_interval] # Doubling times
    plot(d_times, xticks=(1:length(segments),
            string.(dates[[last(s) for s in segments]])),
         xrotation=45, label="", title="$p_head doubling time",
         ylabel="days", xlabel="date", seriestype=st)
end

function plot_whisker_fit(d,v_interval,segments, dates, p_head="Confirmed cases"; model_interval=5, date_density=date_density)
    p=scatter(1:length(dates),d,
        label="data", title="$p_head of Covid-19",
        legend=:topleft, xrotation=45,
        xticks=(1:date_density:length(dates),string.(dates[1:date_density:length(dates)])))
    for i in eachindex(segments)
        plot!(p,first(segments[i]):last(segments[i])+model_interval,
            x->exp(v_interval[i][1]+v_interval[i][2]*x),
            label="")
    end
    p
end


text:
"""
    Orientations

Orientation functions:

 - `vertical_unit_vector`
 - `altitude`
 - `latitude`
 - `longitude`
 - `gravitational_potential`
 - `projection_normal`
 - `projection_tangential`

for orientations:

 - [`NoOrientation`](@ref)
 - [`FlatOrientation`](@ref)
 - [`SphericalOrientation`](@ref)
"""
module Orientations

using CLIMAParameters: AbstractParameterSet
const APS = AbstractParameterSet
using CLIMAParameters.Planet: grav, planet_radius
using StaticArrays
using LinearAlgebra

using ..VariableTemplates
using ..BalanceLaws
using ..MPIStateArrays: MPIStateArray
import ..BalanceLaws: vars_state
using ..DGMethods:
    init_state_auxiliary!, auxiliary_field_gradient!, LocalGeometry

export Orientation, NoOrientation, FlatOrientation, SphericalOrientation

export init_aux!,
    orientation_nodal_init_aux!,
    vertical_unit_vector,
    altitude,
    latitude,
    longitude,
    projection_normal,
    gravitational_potential,
    ∇gravitational_potential,
    projection_tangential,
    sphr_to_cart_vec,
    cart_to_sphr_vec

abstract type Orientation <: BalanceLaw end

#####
##### Fallbacks
#####

function vars_state(m::Orientation, ::Auxiliary, FT)
    @vars begin
        Φ::FT # gravitational potential
        ∇Φ::SVector{3, FT}
    end
end

gravitational_potential(::Orientation, aux::Vars) = aux.orientation.Φ

∇gravitational_potential(::Orientation, aux::Vars) = aux.orientation.∇Φ

function altitude(orientation::Orientation, param_set::APS, aux::Vars)
    FT = eltype(aux)
    return gravitational_potential(orientation, aux) / FT(grav(param_set))
end

function vertical_unit_vector(
    orientation::Orientation,
    param_set::APS,
    aux::Vars,
)
    FT = eltype(aux)
    return ∇gravitational_potential(orientation, aux) / FT(grav(param_set))
end

function projection_normal(
    orientation::Orientation,
    param_set::APS,
    aux::Vars,
    u⃗::AbstractVector,
)
    n̂ = vertical_unit_vector(orientation, param_set, aux)
    return n̂ * (n̂' * u⃗)
end

function projection_tangential(
    orientation::Orientation,
    param_set::APS,
    aux::Vars,
    u⃗::AbstractVector,
)
    return u⃗ .- projection_normal(orientation, param_set, aux, u⃗)
end

function init_aux!(
    m,
    ::Orientation,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
)
    param_set = parameter_set(m)
    init_state_auxiliary!(
        m,
        (m, aux, tmp, geom) ->
            orientation_nodal_init_aux!(m.orientation, param_set, aux, geom),
        state_auxiliary,
        grid,
        direction,
    )

    auxiliary_field_gradient!(
        m,
        state_auxiliary,
        ("orientation.∇Φ",),
        state_auxiliary,
        ("orientation.Φ",),
        grid,
        direction,
    )
end

#####
##### NoOrientation
#####

"""
    NoOrientation

No gravitational force or potential.
"""
struct NoOrientation <: Orientation end

init_aux!(m, ::NoOrientation, state_auxiliary::MPIStateArray, grid, direction) =
    nothing

vars_state(m::NoOrientation, ::Auxiliary, FT) = @vars()

gravitational_potential(::NoOrientation, aux::Vars) = -zero(eltype(aux))
∇gravitational_potential(::NoOrientation, aux::Vars) =
    SVector{3, eltype(aux)}(0, 0, 0)
altitude(orientation::NoOrientation, param_set::APS, aux::Vars) =
    -zero(eltype(aux))

#####
##### SphericalOrientation
#####

"""
    SphericalOrientation <: Orientation

Gravity acts towards the origin `(0,0,0)`, and the gravitational potential is relative
to the surface of the planet.
"""
struct SphericalOrientation <: Orientation end

function orientation_nodal_init_aux!(
    ::SphericalOrientation,
    param_set::APS,
    aux::Vars,
    geom::LocalGeometry,
)
    FT = eltype(aux)
    _grav::FT = grav(param_set)
    _planet_radius::FT = planet_radius(param_set)
    normcoord = norm(geom.coord)
    aux.orientation.Φ = _grav * (normcoord - _planet_radius)
end

# TODO: should we define these for non-spherical orientations?
latitude(orientation::SphericalOrientation, aux::Vars) =
    @inbounds asin(aux.coord[3] / norm(aux.coord, 2))

longitude(orientation::SphericalOrientation, aux::Vars) =
    @inbounds atan(aux.coord[2], aux.coord[1])

"""
    sphr_to_cart_vec(orientation::SphericalOrientation, state::Vars, aux::Vars)

Projects a vector defined based on unit vectors in spherical coordinates to cartesian unit vectors.
"""
function sphr_to_cart_vec(
    orientation::SphericalOrientation,
    vec::AbstractVector,
    aux::Vars,
)
    FT = eltype(aux)
    lat = latitude(orientation, aux)
    lon = longitude(orientation, aux)

    slat, clat = sin(lat), cos(lat)
    slon, clon = sin(lon), cos(lon)

    u = MVector{3, FT}(
        -slon * vec[1] - slat * clon * vec[2] + clat * clon * vec[3],
        clon * vec[1] - slat * slon * vec[2] + clat * slon * vec[3],
        clat * vec[2] + slat * vec[3],
    )

    return u
end

"""
    cart_to_sphr_vec(orientation::SphericalOrientation, state::Vars, aux::Vars)

Projects a vector defined based on unit vectors in cartesian coordinates to a spherical unit vectors.
"""
function cart_to_sphr_vec(
    orientation::SphericalOrientation,
    vec::AbstractVector,
    aux::Vars,
)
    FT = eltype(aux)
    lat = latitude(orientation, aux)
    lon = longitude(orientation, aux)

    slat, clat = sin(lat), cos(lat)
    slon, clon = sin(lon), cos(lon)

    u = MVector{3, FT}(
        -slon * vec[1] + clon * vec[2],
        -slat * clon * vec[1] - slat * slon * vec[2] + clat * vec[3],
        clat * clon * vec[1] + clat * slon * vec[2] + slat * vec[3],
    )

    return u
end

#####
##### FlatOrientation
#####

"""
    FlatOrientation <: Orientation

Gravity acts in the third coordinate, and the gravitational potential is relative to
`coord[3] == 0`.
"""
struct FlatOrientation <: Orientation
    # for Coriolis we could add latitude?
end
function orientation_nodal_init_aux!(
    ::FlatOrientation,
    param_set::APS,
    aux::Vars,
    geom::LocalGeometry,
)
    FT = eltype(aux)
    _grav::FT = grav(param_set)
    @inbounds aux.orientation.Φ = _grav * geom.coord[3]
end

end


text:
"""
    Calculate firm variables and populate struct with steady state values.
"""
function Save_Firm_Solution(l, l_P, l_G, D_I, D_E, D; par::Pars_v3)
    type = [:BS, :LS, :LU, :BU]

    D_S  = sold_data.(D_I, l_G, type; par)
    D_G  = data_gen.(l_G, type; par)
    prof = firm_profits.(l, l_P, l_G, D_S, D_E, type; par)
    K    = knowledge.(l_P, D, type; par)
    res = Firm_sol()
    @set! res.prof_BS_ss = prof[1] # Labor good production big customer base, sophisticated
    @set! res.prof_LS_ss = prof[2] # Labor good production small customer base, sophisticated
    @set! res.prof_LU_ss = prof[3] # Labor good production small customer base, unsophisticated
    @set! res.prof_BU_ss = prof[4] # Labor good production big customer base, unsophisticated

    @set! res.l_BS_ss = l[1] # Labor good production big customer base, sophisticated
    @set! res.l_LS_ss = l[2] # Labor good production small customer base, sophisticated
    @set! res.l_LU_ss = l[3] # Labor good production small customer base, unsophisticated
    @set! res.l_BU_ss = l[4] # Labor good production big customer base, unsophisticated

    @set! res.l_P_BS_ss = l_P[1] # Labor good production big customer base, sophisticated
    @set! res.l_P_LS_ss = l_P[2] # Labor good production small customer base, sophisticated
    @set! res.l_P_LU_ss = l_P[3] # Labor good production small customer base, unsophisticated
    @set! res.l_P_BU_ss = l_P[4] # Labor good production big customer base, unsophisticated

    @set! res.l_G_BS_ss = l_G[1] # Labor good production big customer base, sophisticated
    @set! res.l_G_LS_ss = l_G[2] # Labor good production small customer base, sophisticated
    @set! res.l_G_LU_ss = l_G[3] # Labor good production small customer base, unsophisticated
    @set! res.l_G_BU_ss = l_G[4] # Labor good production big customer base, unsophisticated

    @set! res.D_S_BS_ss = D_S[1] # data sharing big customer base, sophisticated
    @set! res.D_S_LS_ss = D_S[2] # data sharing small customer base, sophisticated
    @set! res.D_S_LU_ss = D_S[3] # data sharing small customer base, unsophisticated
    @set! res.D_S_BU_ss = D_S[4] # data sharing big customer base, unsophisticated

    @set! res.D_G_BS_ss = D_G[1] # data generated big customer base, sophisticated
    @set! res.D_G_LS_ss = D_G[2] # data generated small customer base, sophisticated
    @set! res.D_G_LU_ss = D_G[3] # data generated small customer base, unsophisticated
    @set! res.D_G_BU_ss = D_G[4] # data generated big customer base, unsophisticated

    @set! res.D_BS_ss = D[1] # data generated big customer base, sophisticated
    @set! res.D_LS_ss = D[2] # data generated small customer base, sophisticated
    @set! res.D_LU_ss = D[3] # data generated small customer base, unsophisticated
    @set! res.D_BU_ss = D[4] # data generated big customer base, unsophisticated

    @set! res.D_E_BS_ss = D_E[1] # data generated big customer base, sophisticated
    @set! res.D_E_LS_ss = D_E[2] # data generated small customer base, sophisticated
    @set! res.D_E_LU_ss = D_E[3] # data generated small customer base, unsophisticated
    @set! res.D_E_BU_ss = D_E[4] # data generated big customer base, unsophisticated

    @set! res.D_I_BS_ss = D_I[1] # data generated big customer base, sophisticated
    @set! res.D_I_LS_ss = D_I[2] # data generated small customer base, sophisticated
    @set! res.D_I_LU_ss = D_I[3] # data generated small customer base, unsophisticated
    @set! res.D_I_BU_ss = D_I[4] # data generated big customer base, unsophisticated

    @set! res.K_BS_ss = K[1] # knowledge employed big customer base, sophisticated
    @set! res.K_LS_ss = K[2] # knowledge employed small customer base, sophisticated
    @set! res.K_LU_ss = K[3] # knowledge employed small customer base, unsophisticated
    @set! res.K_BU_ss = K[4] # knowledge employed big customer base, unsophisticated

    return res
end

function weight_of_type(type::Symbol; par::Pars_v3)
    if type == :LU # little, unsophisticated
        (1 - par.w_B) * (1 - par.w_S)
    elseif type == :LS # little, sophisticated
        (1 - par.w_B) * par.w_S
    elseif type == :BU # big, unsophisticated
        par.w_B * (1 - par.w_S)
    elseif type == :BS # big, sophisticated
        par.w_B * par.w_S
    else
        error("Firm type needs to be either LU, LS, BU or BS")
    end
end

"""
    Data Generation Productivity / Customer Base depending on type.
"""
function A_G_of_type(type::Symbol; par::Pars_v3)
    if type == :LU # little, unsophisticated
        par.A_G_S
    elseif type == :LS # little, sophisticated
        par.A_G_S
    elseif type == :BU # big, unsophisticated
        par.A_G_B
    elseif type == :BS # big, sophisticated
        par.A_G_B
    else
        error("Firm type needs to be either LU, LS, BU or BS")
    end
end

"""
    Data Processing Productivity depending on type.
"""
function A_P_of_type(type::Symbol; par::Pars_v3)
    if type == :LU # little, unsophisticated
        par.A_P_U
    elseif type == :LS # little, sophisticated
        par.A_P_S
    elseif type == :BU # big, unsophisticated
        par.A_P_U
    elseif type == :BS # big, sophisticated
        par.A_P_S
    else
        error("Firm type needs to be either LU, LS, BU or BS")
    end
end

"""
    Squares the minimizer and returns the solution values
"""
function get_solution(sol; par::Pars_v3, type::Symbol)
    if length(sol.minimizer) == 4
        l_P = sol.minimizer[1]^2
        l_G = sol.minimizer[2]^2
        D_E = sol.minimizer[3]^2
        D_I = sol.minimizer[4]^2
        D = bundle(D_I, D_E; par, type)
        l = labor_demand(l_P, D; par, type)
    elseif length(sol.minimizer) == 3
        l_P = sol.minimizer[1]^2
        l_G = sol.minimizer[2]^2
        D_E = sol.minimizer[3]^2
        D_I = data_gen(l_G; par, type)
        D = bundle(D_I, D_E; par, type)
        l = labor_demand(l_P, D; par, type)
    else
        error("I should never be here")
    end

    return l_P, l_G, l, D_E, D_I, D
end

function get_solution_nl(sol; par::Pars_v3, type::Symbol)
    if length(sol.zero) == 4
        l_P = sol.zero[1]^2
        l_G = sol.zero[2]^2
        D_E = sol.zero[3]^2
        D_I = sol.zero[4]^2
        D = bundle(D_I, D_E; par, type)
        l = labor_demand(l_P, D; par, type)
    elseif length(sol.zero) == 3
        l_P = sol.zero[1]^2
        l_G = sol.zero[2]^2
        D_E = sol.zero[3]^2
        D_I = data_gen(l_G; par, type)
        D = bundle(D_I, D_E; par, type)
        l = labor_demand(l_P, D; par, type)
    else
        error("I should never be here")
    end

    return l_P, l_G, l, D_E, D_I, D
end

function get_solution_nl2(sol; par::Pars_v3, type::Symbol)
    if length(sol.zero) == 4
        l_P = sol.zero[1]^2
        l_G = sol.zero[2]^2
        D_E = sol.zero[3]^2
        D_I = sol.zero[4]^2
        D = bundle(D_I, D_E; par, type)
        l = labor_demand(l_P, D; par, type)
        return l_P, l_G, l, D_E, D_I, D
    elseif length(sol.zero) == 3
        l_P = sol.zero[1]^2
        l_G = sol.zero[2]^2
        D_E = sol.zero[3]^2
        D_I = data_gen(l_G; par, type)
        D = bundle(D_I, D_E; par, type)
        l = labor_demand(l_P, D; par, type)
        return l_P, l_G, l, D_E, D_I, D
    elseif length(sol.zero) == 8
        l_P = sol.zero[1]^2
        l_G = sol.zero[2]^2
        D_E = sol.zero[3]^2
        D_I = sol.zero[4]^2
        l_P2 = sol.zero[5]^2
        l_G2 = sol.zero[6]^2
        D_E2 = sol.zero[7]^2
        D_I2 = sol.zero[8]^2
        D = bundle(D_I, D_E; par, type)
        D2 = bundle(D_I2, D_E2; par, type)
        l = labor_demand(l_P, D; par, type, sec_period=:no)
        l2 = labor_demand(l_P2, D2; par, type, sec_period=:yes)
        return l_P, l_G, l, D_E, D_I, D, l_P2, l_G2, l2, D_E2, D_I2, D2
    elseif length(sol.zero) == 6 #! Need to take also into account an intermediate case in which firms switch from data sharing to not sharing... but unlikely to arise in our setting
        l_P = sol.zero[1]^2
        l_G = sol.zero[2]^2
        D_E = sol.zero[3]^2
        l_P2 = sol.zero[4]^2
        l_G2 = sol.zero[5]^2
        D_E2 = sol.zero[6]^2
        D_I = data_gen(l_G; par, type)
        D_I2 = data_gen(l_G2; par, type)
        D = bundle(D_I, D_E; par, type)
        D2 = bundle(D_I2, D_E2; par, type)
        l = labor_demand(l_P, D; par, type, sec_period = :no)
        l2 = labor_demand(l_P2, D2; par, type, sec_period = :yes)
    else
        error("I should never be here")
    end
end

function set_a_Y(compl_prod, σ, ζ)
    if compl_prod == :yes
        α_Y = (σ*ζ - σ + 1) / (σ*ζ) # auxiliary paramter
    elseif compl_prod == :no
        α_Y = 0
    else
        error("Field compl_prod needs to be either :no or :yes")
    end
    return α_Y
end

"""
    Return the correct gamma
"""
function gamma_of_type(par::Pars_v3; type::Symbol)
    if type in [:BS, :LS]
        par.γ_S
    elseif type in [:BU, :LU]
        par.γ_U
    else
        error("type needs to be either :BS, :LS, :BU or :LU")
    end
end

function K_last_of_type(par::Pars_v3; type::Symbol)
    if type == :LU # little, unsophisticated
        par.K_last_LU
    elseif type == :LS # little, sophisticated
        par.K_last_LS
    elseif type == :BU # big, unsophisticated
        par.K_last_BU
    elseif type == :BS # big, sophisticated
        par.K_last_BS
    else
        error("Firm type needs to be either LU, LS, BU or BS")
    end
end

function D_last_of_type(par::Pars_v3; type::Symbol)
    if type == :LU # little, unsophisticated
        par.D_last_LU
    elseif type == :LS # little, sophisticated
        par.D_last_LS
    elseif type == :BU # big, unsophisticated
        par.D_last_BU
    elseif type == :BS # big, sophisticated
        par.D_last_BS
    else
        error("Firm type needs to be either LU, LS, BU or BS")
    end
end


text:
using PowerSeries

function series_loop()
  accum = 0.0
  for i = 1:10000
    accum += polyder(sin(series(1.0*i, 1.0)))
  end
  return accum
end

function bare_loop()
  accum = 0.0
  for i = 1:10000
    accum += cos(1.0*i)
  end
  return accum
end

# Warm up JIT
@assert series_loop() == bare_loop()

tseries = @elapsed series_loop()
tbare = @elapsed bare_loop()

println("Series time")
println("Series", ", ", "Bare", ", ", "Ratio")
println(tseries, ", ", tbare, ", ", tseries/tbare)


text:
__precompile__()

module PPT

const codefile = joinpath(dirname(@__FILE__), "literate_org_tangled_code.jl")
if isfile(codefile)
    include(codefile)
else
    error("PPT not properly installed. Please run Pkg.build(\"PPT\") then restart Julia.")
end

end # module


text:
""" Computes the (4) coefficients of the cubic interpolant on the two
    points x using the 4 values for function and derivative in dy
"""
function coefHIN3(x,dy)         # Cubic Hermite interpolation (Newton)
    # vector x: [x_1,x_2]
    # matrix dy contains 2 values, function and derivative each of the 2 points
    n, m = size(dy)
    dd = DividedDifference2(x,dy)   # get the divided difference array

    α = x[2] - x[1]
    a = zeros(typeof(x[1]),4)
    a[1] = dd[1,1]
    a[2] = dd[2,2]
    a[3] = dd[3,3] - dd[4,4]*α
    a[4] = dd[4,4]

    return a
end

""" Computes the (6) coefficients of the quintic interpolant on the two
    points x using the 6 values for function and first and second
    derivative in dy
"""
function coefHIN5(x,dy)         # Quintic Hermite interpolation (Newton)
    # vector x: [x_1,x_2]
    # matrix dy contains 3 values, function, derivative and second derivatives of each of the 2 points
    n, m = size(dy)
    dd = DividedDifference2(x,dy)   # get the divided difference array

    α = x[2] - x[1]
    a = zeros(typeof(x[1]),6)
    a[1] = dd[1,1]
    a[2] = dd[2,2]
    a[3] = dd[3,3]
    a[4] = dd[4,4] - α*dd[5,5] + α^2*dd[6,6]
    a[5] = dd[5,5] - 2*α*dd[6,6]
    a[6] = dd[6,6]

    return a
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export GRULayer
sigmoid(x)  = 1.0f0 / (1.0f0 + exp(-x))
∇sigmoid(x) = 1.0f0 - x * x
∇tanh(x)    = x * (1.0f0 - x)

@neuron type GRUNeuron
    W_z  :: DenseArray{Float32}
    ∇W_z :: DenseArray{Float32}
    U_z  :: DenseArray{Float32}
    ∇U_z :: DenseArray{Float32}
    b_z  :: DenseArray{Float32}
    ∇b_z :: DenseArray{Float32}

    W_r  :: DenseArray{Float32}
    ∇W_r :: DenseArray{Float32}
    U_r  :: DenseArray{Float32}
    ∇U_r :: DenseArray{Float32}
    b_r  :: DenseArray{Float32}
    ∇b_r :: DenseArray{Float32}

    W_h  :: DenseArray{Float32}
    ∇W_h :: DenseArray{Float32}
    U_h  :: DenseArray{Float32}
    ∇U_h :: DenseArray{Float32}
    b_h  :: DenseArray{Float32}
    ∇b_h :: DenseArray{Float32}

    r    :: Batch{Float32}
    z    :: Batch{Float32}
    hh   :: Batch{Float32}
end
function GRUNeuron(
        W_z  :: DenseArray{Float32},
        ∇W_z :: DenseArray{Float32},
        U_z  :: DenseArray{Float32},
        ∇U_z :: DenseArray{Float32},
        b_z  :: DenseArray{Float32},
        ∇b_z :: DenseArray{Float32},

        W_r  :: DenseArray{Float32},
        ∇W_r :: DenseArray{Float32},
        U_r  :: DenseArray{Float32},
        ∇U_r :: DenseArray{Float32},
        b_r  :: DenseArray{Float32},
        ∇b_r :: DenseArray{Float32},

        W_h  :: DenseArray{Float32},
        ∇W_h :: DenseArray{Float32},
        U_h  :: DenseArray{Float32},
        ∇U_h :: DenseArray{Float32},
        b_h  :: DenseArray{Float32},
        ∇b_h :: DenseArray{Float32}
    )
    GRUNeuron(
        W_z, ∇W_z, U_z, ∇U_z, b_z, ∇b_z,
        W_r, ∇W_r, U_r, ∇U_r, b_r, ∇b_r,
        W_h, ∇W_h, U_h, ∇U_h, b_h, ∇b_h,
        Batch(0.0f0),
        Batch(0.0f0),
        Batch(0.0f0))
end

@neuron forward(neuron::GRUNeuron) do
    x_z = neuron.b_z[1]
    x_r = neuron.b_r[1]
    x_h = neuron.b_h[1]
    for i in 1:length(neuron.inputs[1])
        x_z += neuron.W_z[i] * neuron.inputs[1][i]
        x_r += neuron.W_r[i] * neuron.inputs[1][i]
        x_h += neuron.W_h[i] * neuron.inputs[1][i]
    end

    u_z = 0.0f0
    u_r = 0.0f0
    for i = 1:length(neuron.inputs[2])
        u_z += neuron.U_z[i] * neuron.inputs[2][i]
        u_r += neuron.U_r[i] * neuron.inputs[2][i]
    end

    neuron.z = sigmoid(x_z + u_z)
    neuron.r = sigmoid(x_r + u_r)
    u_h = 0.0f0
    for i = 1:length(neuron.inputs[2])
        u_h += neuron.U_h[i] * neuron.r * neuron.inputs[2][i]
    end
    neuron.hh = tanh(x_h + u_h)
    neuron.value = neuron.z * neuron.inputs[2][neuron.index] + (1 - neuron.z) * neuron.hh
end

@neuron backward(neuron::GRUNeuron) do
    ∇_z = neuron.∇ * neuron.hh
    ∇hh = neuron.∇ * (1 - neuron.z)
    ∇z  = ∇sigmoid(neuron.∇ * neuron.inputs[2][neuron.index])
    neuron.∇inputs[2][neuron.index] = neuron.∇ * neuron.z
    ∇h = ∇tanh(∇hh)
    ∇r = 0.0f0
    for i = 1:length(neuron.inputs[2])
        neuron.∇U_h[i] += ∇h * neuron.r * neuron.inputs[2][i]
        neuron.∇inputs[2][i] += ∇h * neuron.r * neuron.U_h[i]
        ∇r += ∇h * neuron.U_h[i] * neuron.inputs[2][i]
    end
    ∇r = ∇sigmoid(∇r)
    for i = 1:length(neuron.inputs[2])
        neuron.∇U_z[i] += ∇z * neuron.inputs[2][i]
        neuron.∇inputs[2][i] += ∇z * neuron.U_z[i]

        neuron.∇U_r[i] += ∇r * neuron.inputs[2][i]
        neuron.∇inputs[2][i] += ∇r * neuron.U_r[i]
    end
    for i = 1:length(neuron.inputs[2])
        neuron.∇W_z[i] += ∇z * neuron.inputs[1][i]
        neuron.∇inputs[1][i] += ∇z * neuron.W_z[i]

        neuron.∇W_r[i] += ∇r * neuron.inputs[1][i]
        neuron.∇inputs[2][i] += ∇r * neuron.W_r[i]

        neuron.∇W_h[i] += ∇h * neuron.inputs[1][i]
        neuron.∇inputs[2][i] += ∇h * neuron.W_h[i]
    end
    neuron.∇b_h[1] = ∇h
    neuron.∇b_r[1] = ∇r
    neuron.∇b_z[1] = ∇z
end

function GRULayer(name::Symbol, net::Net, input_ensemble::AbstractEnsemble, num_outputs::Int)
    weight_init = xavier
    in_out_shape = (length(input_ensemble), num_outputs)
    out_out_shape = (num_outputs, num_outputs)

    W_z = weight_init(Float32, in_out_shape...)
    ∇W_z = zeros(Float32, in_out_shape...)

    U_z = weight_init(Float32, out_out_shape...)
    ∇U_z = zeros(Float32, out_out_shape...)

    b_z = zeros(Float32, 1, num_outputs)
    ∇b_z = zeros(Float32, 1, num_outputs)

    W_r = weight_init(Float32, in_out_shape...)
    ∇W_r = zeros(Float32, in_out_shape...)

    U_r = weight_init(Float32, out_out_shape...)
    ∇U_r = zeros(Float32, out_out_shape...)

    b_r = zeros(Float32, 1, num_outputs)
    ∇b_r = zeros(Float32, 1, num_outputs)

    W_h = weight_init(Float32, in_out_shape...)
    ∇W_h = zeros(Float32, in_out_shape...)

    U_h = weight_init(Float32, out_out_shape...)
    ∇U_h = zeros(Float32, out_out_shape...)

    b_h = zeros(Float32, 1, num_outputs)
    ∇b_h = zeros(Float32, 1, num_outputs)

    neurons = Array(GRUNeuron, num_outputs)
    for i in 1:num_outputs
        neurons[i] = GRUNeuron(
            view(W_z, :, i), view(∇W_z, :, i),
            view(U_z, :, i), view(∇U_z, :, i),
            view(b_z, :, i), view(∇b_z, :, i),

            view(W_r, :, i), view(∇W_r, :, i),
            view(U_r, :, i), view(∇U_r, :, i),
            view(b_r, :, i), view(∇b_r, :, i),

            view(W_h, :, i), view(∇W_h, :, i),
            view(U_h, :, i), view(∇U_h, :, i),
            view(b_h, :, i), view(∇b_h, :, i),
        )
    end

    ens = Ensemble(net, name, neurons, [Param(net, name,:W_z, 1.0),
                                        Param(net, name,:U_z, 1.0),
                                        Param(net, name,:b_z, 2.0),
                                        Param(net, name,:W_r, 1.0),
                                        Param(net, name,:U_r, 1.0),
                                        Param(net, name,:b_r, 2.0),
                                        Param(net, name,:W_h, 1.0),
                                        Param(net, name,:U_h, 1.0),
                                        Param(net, name,:b_h, 2.0)])
    add_connections(net, input_ensemble, ens,
                    (i) -> (tuple([Colon() for d in size(input_ensemble)]... )))
    add_connections(net, ens, ens,
                    (i) -> (tuple([Colon() for d in size(input_ensemble)]... ));
                    recurrent=true)
    ens

end


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ec1fe200-8d73-11eb-0ba3-e593aa79dab2
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ 973abf88-9125-41c7-8372-61efb2976fc2
using SymPy

# ╔═╡ da340240-f845-4591-94c3-9ffce6a33d8a
begin
	using Random
	Random.seed!(123)
	A₄=rand(-4:4,4,4)
end

# ╔═╡ 60d351a5-0ef2-47cd-a5a5-5ecf18c4e898
begin
	using SpecialMatrices
	using Polynomials
end

# ╔═╡ bb8416e3-b682-4733-af6d-36fd48f0a7b8
using Plots

# ╔═╡ 10910df9-240f-4e6b-8076-1a9be4d66ba1
md"""
# Eigenvalue Decomposition - Definitions and Facts

__Prerequisites__

The reader should be familiar with basic linear algebra concepts.

__Competences__

The reader should be able to understand and check the facts about eigenvalue decomposition.

__Selected references__

There are many excellent books on the subject. Here we list a few:

* J.W. Demmel, _Applied Numerical Linear Algebra_, SIAM, Philadelphia, 1997.
* G. H. Golub and C. F. Van Loan, _Matrix Computations_, 4th ed., The John Hopkins University Press, Baltimore, MD, 2013.
* N. Higham, _Accuracy and Stability of Numerical Algorithms_, SIAM, Philadelphia, 2nd ed., 2002.
* L. Hogben, ed., _Handbook of Linear Algebra_, CRC Press, Boca Raton, 2014.
* B. N. Parlett, _The Symmetric Eigenvalue Problem_, Prentice-Hall, Englewood Cliffs, NJ, 1980, also SIAM, Philadelphia, 1998.
* G. W. Stewart, _Matrix Algorithms, Vol. II: Eigensystems_, SIAM, Philadelphia, 2001.
* L. N. Trefethen and D. Bau, III, _Numerical Linear Algebra_, SIAM, Philadelphia, 1997.
* J. H. Wilkinson, _The Algebraic Eigenvalue Problem_, Clarendon Press, Oxford, U.K.,  1965.

"""

# ╔═╡ b87354b9-ee7b-4ce0-bf26-b4e71c6be6a9
md"""
# General matrices

For more details and the proofs of the Facts below, see
[L. M. DeAlba, Determinants and Eigenvalues, pp. 4.1-1.15](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

## Definitions

We state the basic definitions:

Let $\mathbb{F}=\mathbb{R}$ or $F=\mathbb{C}$ and let $A\in \mathbb{F}^{n\times n}$ with elements $a_{ij}\in \mathbb{F}$.

An element $\lambda \in \mathbb{F}$ is an __eigenvalue__ of $A$ if
$\exists x\in \mathbb{F}^n$, $x\neq 0$ such that

$$Ax=\lambda x,$$

and $x$ is an __eigenvector__ of $\lambda$.

__Characteristic polynomial__ of $A$ is $p_A(x)=\det(A-xI)$.

__Algebraic multiplicty__, $\alpha(\lambda)$, is the multiplicity of $\lambda$ as a root of $p_A(x)$.

__Spectrum of $A$__, $\sigma(A)$, is the multiset of all eigenvalues of $A$, with each eigenvalue appearing $\alpha(A)$ times.

__Spectral radius__ of $A$ is

$$\rho(A)=\max\{|\lambda|, \lambda \in \sigma(A)\}.$$

__Eigenspace__ of $\lambda$ is

$$
E_{\lambda}(A)=\ker(A-\lambda I).$$

__Geometric multiplicity__ of $\lambda$ is

$$\gamma(\lambda)=\dim(E_{\lambda}(A)).$$

 $\lambda$ is __simple__ if $\alpha(\lambda)=1$.

 $\lambda$ is __semisimple__ if $\alpha(\lambda)=\gamma(\lambda)$.

 $A$ is __nonderogatory__ if $\gamma(\lambda)=1$ for all $\lambda$.

 $A$ is __nondefective__ if every $\lambda$ is semisimple.

 $A$ is __diagonalizable__ if there exists nonsingular $B$ matrix and diagonal matrix $D$ such that $A=BDB^{-1}$.

__Trace__ of $A$ is

$$\mathop{\mathrm{tr}}(A)=\sum_i a_{ii}.$$

 $Q\in\mathbb{C}^{n\times n}$ is __unitary__ if $Q^*Q=QQ^*=I$,
where $Q^*=(\bar Q)^T$.

__Schur decomposition__ of $A$ is $A=QTQ^*$,
where $Q$ is unitary and $T$ is upper triangular.

 $A$ and $B$ are __similar__ if $B=QAQ^{-1}$ for some nonsingular matrix $Q$.

 $A$ is __normal__ if $AA^*=A^*A$.
"""

# ╔═╡ 31a71b9e-1dba-4a69-b87e-9202292eb3ed
md"""
## Facts

There are many facts related to the eigenvalue problem for general matrices. We state some basic ones:

1.  $\lambda\in\sigma(A) \Leftrightarrow p_A(\lambda)=0$.

1. __Cayley-Hamilton Theorem.__ $p_A(A)=0$.

2. A simple eigenvalue is semisimple.

3.  $\mathop{\mathrm{tr}}(A)=\sum_{i=1}^n \lambda_i$.

4.  $\det(A)=\prod_{i=1}^n \lambda_i$.

5.  $A$ is singular $\Leftrightarrow$ $\det(A)=0$ $\Leftrightarrow$ $0\in\sigma(A)$.

7. If $A$ is triangular, $\sigma(A)=\{a_{11},a_{22},\ldots,a_{nn}\}$.

6. For $A\in\mathbb{C}^{n\times n}$, $\lambda\in\sigma(A)$ $\Leftrightarrow$ $\bar\lambda\in\sigma(A^*)$.

8. __Corollary of the Fundamental theorem of algebra.__ For $A\in\mathbb{R}^{n\times n}$, $\lambda\in\sigma(A)$ $\Leftrightarrow$ $\bar\lambda\in\sigma(A)$.

9. If $(\lambda,x)$ is an eigenpair of a nonsingular $A$, then $(1/\lambda,x)$ is an eigenpair of $A^{-1}$.

10. Eigenvectors corresponding to distinct eigenvalues are linearly independent.

11.  $A$ is diagonalizable $\Leftrightarrow$ $A$ is nondefective $\Leftrightarrow$ $A$ has $n$ linearly independent eigenvectors.

12. Every $A$ has Schur decomposition. Moreover, $T_{ii}=\lambda_i$.  _(For proof see GVL p.351 (375).)_

13. If $A$ is normal, matrix $T$ from its Schur decomposition is normal. Consequently:

    *  $T$ is diagonal, and has eigenvalues of $A$ on diagonal,
    * matrix $Q$ of the Schur decomposition is the unitary matrix of eigenvectors,
    * all eigenvalues of $A$ are semisimple and $A$ is nondefective.

14. Real matrix has a __real Schur decomposition__, $Q^TAQ=T$, where $Q$ is real orthogonal and $T$ is upper block-triangular with $1\times 1$ and $2 \times 2$ blocks on the diagonal. The $1\times 1$ blocks correspond to real eigenvalues, and the $2\times 2$ blocks correspond to pairs of complex conjugate eigenvalues.

15. If $A$ and $B$ are similar, $\sigma(A)=\sigma(B)$. Consequently, $\mathop{\mathrm{tr}}(A)=\mathop{\mathrm{tr}}(B)$ and $\det(A)=\det(B)$.

16. Eigenvalues and eigenvectors are continous and differentiable: if $\lambda$ is a simple eigenvalue of $A$ and $A(\varepsilon)=A+\varepsilon E$ for some $E\in F^{n\times n}$, for small $\varepsilon$ there exist differentiable functions $\lambda(\varepsilon)$ and $x(\varepsilon)$ such that $A(\varepsilon) x(\varepsilon) = \lambda(\varepsilon) x(\varepsilon)$.

17. Classical motivation for the eigenvalue problem is the following: consider the system of linear differential equations with constant coefficients, $\dot y(t)=Ay(t)$. If the solution is $y=e^{\lambda t} x$ for some constant vector $x$, then

$$\lambda e^{\lambda t} x=Ae^{\lambda t} x \quad \textrm{or} \quad Ax=\lambda x.$$

"""

# ╔═╡ 2f5aee03-3c2e-4920-9d10-c509b06be627
md"""
## Examples

We shall illustrate above Definitions and Facts on several small examples, using symbolic computation.
"""

# ╔═╡ 37d24090-8d78-11eb-059d-737014bb6588
md"
### Defective eigenvalue
"

# ╔═╡ 184bc4ec-c566-4f21-99b3-b6e8cd328b47
A=[-3 7 -1; 6 8 -2; 72 -28 19]

# ╔═╡ 5d0aac9e-4f48-4f09-956f-5033a1fb8f11
x=symbols("x")

# ╔═╡ 73ece9ec-97dd-46bc-8e17-10ce7bfc18bc
begin
	using LinearAlgebra
	eye(n)=Matrix{Int}(I,n,n)
	A-x*eye(3)
end

# ╔═╡ ce83f04c-e915-4599-beeb-d4a0a681bc17
# Characteristic polynomial p_A(λ)
p(x)=det(A-x*eye(3))

# ╔═╡ ff5044c0-8d76-11eb-3951-6bb4b5caa7d0
p(x)

# ╔═╡ 55171a69-3557-4ba5-8a4a-813f38c03dd3
# Characteristic polynomial in nicer form
factor(p(x))

# ╔═╡ 6df42190-a29e-4c68-8684-07d8d0c9d200
λ=solve(p(x),x)

# ╔═╡ b94f7fcd-c6f4-4f4f-924a-20f00e82d94e
md"""
The eigenvalues are $\lambda_1=-6$ and $\lambda_2=15$ with algebraic multiplicities
$\alpha(\lambda_1)=1$ and $\alpha(\lambda_2)=2$.
"""

# ╔═╡ 1b90d50a-bf97-4b1a-9580-d729b03c914e
nullspace(float.(A-λ[1]*I))

# ╔═╡ 15b60cd5-eb30-43bb-bbcf-dd533e29c35f
nullspace(float.(A-λ[2]*I))

# ╔═╡ 3619fca4-1374-4adc-b29c-8058b2e147e0
md"""
The geometric multiplicities are $\gamma(\lambda_1)=1$ and $\gamma(\lambda_2)=1$. Thus, $\lambda_2$ is not semisimple, therefore $A$ is defective and not diagonalizable.
"""

# ╔═╡ d46018c0-a5f8-4dd0-8f3d-73322be792e9
# Trace and determinant
tr(A), λ[1]+λ[2]+λ[2]

# ╔═╡ d387b26f-bec1-47a9-b1c6-f20c43f0bde8
det(A), λ[1]*λ[2]*λ[2]

# ╔═╡ 703da5c4-cd09-4832-b64e-5fda75099db5
# Schur decomposition
S=schur(A)

# ╔═╡ e26f113c-dd98-4e97-b495-581de1b431af
fieldnames(typeof(S))

# ╔═╡ 62779bfc-ea5d-4eec-94ad-2641e11fc44d
S.T

# ╔═╡ bdf079a6-beee-44f5-8200-613e54b3a8dc
S.Z

# ╔═╡ e1f5027c-1080-4c95-81cc-a75dde641b88
# Orthogonality
S.Z'*S.Z

# ╔═╡ c71a7760-b4c6-404e-b5ef-f4879cb787cd
begin
	# Similar matrices
	M=rand(-5:5,3,3)
	eigvals(M*A*inv(M)), tr(M*A*inv(M)), det(M*A*inv(M))
end

# ╔═╡ 343c8f45-6226-4361-a149-5313a47c355c
md"
### Real Schur decomposition
"

# ╔═╡ 41a05c50-f6df-4b93-bd8a-50412829105e
begin
	# Fact 14
	Random.seed!(326)
	A₀=rand(-9:9,6,6)
end

# ╔═╡ cd8fb666-d2f2-423c-ac39-405c842e97ba
schur(A₀)

# ╔═╡ 7aa7ac89-837c-4ee8-b9f1-8ad0256f25e1
md"""
### Diagonalizable matrix

This matrix is nondefective and diagonalizable.
"""

# ╔═╡ c3724c4f-142b-4863-8c04-4a0c15411f40
A₁=[57 -21 21; -14 22 -7; -140 70 -55]

# ╔═╡ 792f632e-fdc2-4016-baaa-b366704f02fd
p₁(x)=factor(det(A₁-x*eye(3)))

# ╔═╡ a8664a40-8d78-11eb-325e-41c6497be282
factor(p₁(x))

# ╔═╡ 1d42245a-2ed6-4cc7-8cbd-2b6aa4e64861
λ₁=solve(p₁(x),x)

# ╔═╡ f7aae752-70c4-4a10-b96e-206c13a6ced1
nullspace(float.(A₁-λ₁[2]*I))

# ╔═╡ 60525ba5-2c76-406e-8eb3-63ced2e00bdd
# F=schur(A)
schur(A₁).T

# ╔═╡ 189544ef-122c-4e30-bfd6-cfc126fc0320
md"""
### Symbolic computation for $n=4$

Let us try some random examples of dimension $n=4$ (the largest $n$ for which we can compute eigevalues symbolically).
"""

# ╔═╡ e72bf173-7e55-4476-91ff-789847dc53cf
p₄(x)=factor(det(A₄-x*eye(4)))

# ╔═╡ da4a5d70-92d8-11eb-0635-11c530ff7d21
p₄(x)

# ╔═╡ ec7fc5c0-92d8-11eb-16df-af1133ccb42a
f(x)=x^4-5*x^3+2*x^2-30*x-30*x^0

# ╔═╡ 29de5fd0-92d9-11eb-09c0-0356001caaef
# Cayley-Hamilton teorem
f(A₄)

# ╔═╡ d43480a3-9247-48a4-8a3b-e006ab58539e
λ₄=solve(p₄(x),x)

# ╔═╡ 5995cd31-a75d-4718-8808-4ada01920344
λ₄[1]

# ╔═╡ 57e2742e-8d79-11eb-0a33-f7250ac02840
λ₄[2]

# ╔═╡ 5e069120-8d79-11eb-3c00-0b5249d1f2dc
λ₄[3]

# ╔═╡ 63167eee-8d79-11eb-29c2-2febbd7cbd93
λ₄[4]

# ╔═╡ e051d2eb-3972-4c97-b7eb-8a18629f6171
md"""
Since all eigenvalues are distinct, they are all simple and the matrix is diagonalizable.
With high probability, all eigenvalues of a random matrix are simple.

Do not try to use `nullspace()` here.

Symbolic computation does not work well with floating-point numbers - in the following example, the determinant of $A-xI$ is a rational function:
"""

# ╔═╡ 22f0a627-beb2-4b06-aa94-ee1e90fe8a6c
begin
	A₂=rand(4,4)
    factor(det(A₂-x*eye(4)))
end

# ╔═╡ da6baba0-92d9-11eb-1603-7961ad1f49c8
A₂

# ╔═╡ c6ab93ab-b336-4bb4-9dda-b725108f7827
md"""
Let us try `Rational` numbers:
"""

# ╔═╡ ffcb919f-9754-4650-92f1-7dd163c5df25
Aᵣ=map(Rational,A₂)

# ╔═╡ 32ce58b9-783b-4739-9001-0fe581b2953d
begin
	pᵣ(x)=factor(det(Aᵣ-x*eye(4)))
	pᵣ(x)
end

# ╔═╡ a7e8d711-57a6-477a-a641-b929789dbadf
λᵣ=solve(pᵣ(x),x)

# ╔═╡ 5781508c-388a-4df0-a6f9-0fe312a81335
# For example:
λᵣ[1]

# ╔═╡ 90b4341d-7964-4e41-a268-8a047a058286
md"""
### Circulant matrix

For more details, see
[A. Böttcher and I. Spitkovsky, Special Types of Matrices, pp. 22.1-22.20](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

Given $a_0,a_1,\ldots,a_{n-1}\in \mathbb{C}$, the __circulant matrix__ is

$$C(a_0,a_1,\ldots,a_{n-1})=\begin{bmatrix}
a_0 & a_{n-1} & a_{n-2} & \cdots & a_{1} \\
a_1 & a_0 & a_{n-1} & \cdots & a_{2} \\
a_2 & a_1 & a_{0} & \cdots & a_{3} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
a_{n-1} & a_{n-2} & a_{n-3} & \cdots & a_{0}
\end{bmatrix}.$$

Let $a(z)=a_0+a_1z+a_2z^2+\cdots +a_{n-1}z^{n-1}$ be the associated complex polynomial.

Let $\omega_n=\exp\big(\displaystyle\frac{2\pi i}{n}\big)$. The __Fourier matrix__ of order $n$ is

$$F_n=\frac{1}{\sqrt{n}} \bigg[ \omega_n^{(j-1)(k-1)}\bigg]_{j,k=1}^n=
\frac{1}{\sqrt{n}} \begin{bmatrix}
1 & 1 & \cdots & 1 \\
1& \omega_n & \omega_n^2 & \cdots \omega_n^{n-1} \\
1& \omega_n^2 & \omega_n^4 & \cdots \omega_n^{2(n-1)} \\
\vdots & \vdots & \ddots & \vdots \\
1& \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots \omega_n^{(n-1)(n-1)}
\end{bmatrix}.$$

Fourier matrix is unitary.
Fourier matrix is a Vandermonde matrix,

$$F_n=\displaystyle\frac{1}{\sqrt{n}} V(1,\omega_n,\omega_n^2,\ldots, \omega_n^{n-1}).$$

Circulant matrix is normal and, thus, unitarily diagonalizable, with the eigenvalue decomposition

$$
C(a_0,a_1,\ldots,a_{n-1})=F_n^*\mathop{\mathrm{diag}}[(a(1),a(\omega_n),a(\omega_n^2),\ldots,
a(\omega_n^{n-1})]F_n.$$

We shall use the package [SpecialMatrices.jl](https://github.com/JuliaMatrices/SpecialMatrices.jl).
"""

# ╔═╡ 16795c64-2994-4fca-89f3-f035fb6a2e45
# using Pkg; Pkg.add(PackageSpec(name="SpecialMatrices",rev="master"))

# ╔═╡ 76ffa073-59b8-4118-aa32-522f4ff05979
varinfo(SpecialMatrices)

# ╔═╡ f85ce269-7555-48ed-a2fd-1193d3f33f40
begin
	Random.seed!(124)
	n=6
	a=rand(-9:9,n)
end

# ╔═╡ e6544939-c8c2-464c-83a0-a6bf1751efa1
C=Circulant(a)

# ╔═╡ b559d56b-02bb-47d3-8f0f-8dae09d918d8
# Check for normality
C*C'-C'*C

# ╔═╡ 5e0adf59-29d3-47ac-8f3f-52929cc4e1c1
pᵪ=Polynomial(a)

# ╔═╡ 83ba2d1b-66de-41bf-8ded-7070bed15587
# Polynomials are easy to manipulate
pᵪ*pᵪ

# ╔═╡ b0255c92-3242-442a-a670-1bfef38e955c
# n-th complex root of unity
ω=exp(2*π*im/n)

# ╔═╡ c6116576-5b25-45b2-bc71-a889f302b5ac
begin
	v=[ω^k for k=0:n-1]
	Fₙ=Vandermonde(v)/√n
end

# ╔═╡ 3de2dc38-ba95-479e-89bc-77cb520daf2f
Λᵪ=Fₙ*C*Fₙ'

# ╔═╡ ad32c51c-7cbb-42f9-abc8-561828df1379
# Check
[diag(Λᵪ) pᵪ.(v) eigvals(Matrix(C))]

# ╔═╡ a7999464-7060-4345-b396-e2e15fcbcaac
md"""

# Hermitian and real symmetric matrices

For more details and the proofs of the Facts below, see
[W. Barrett, Hermitian and Positive Definite Matrices, pp. 9.1-9.13](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

## Definitions

Matrix $A\in \mathbb{C}^{n\times n}$ is __Hermitian__ or __self-adjoint__ if $A^*=A$,
or element-wise, $\bar a_{ij}=a_{ji}$. We say $A\in\mathcal{H}_n$.

Matrix $A\in \mathbb{R}^{n\times n}$ is __symmetric__ if $A^T=A$, or element-wise, $a_{ij}=a_{ji}$.
We say $A\in\mathcal{S}_n$.

__Rayleigh qoutient__ of $A\in\mathcal{H}_n$ and nonzero vector $x\in\mathbb{C}^n$ is

$$
R_A(x)=\frac{x^*Ax}{x^*x}.$$

Matrices $A,B \in\mathcal{H}_n$ are __congruent__ if there exists nonsingular matrix $C$ such that
$B=C^*AC$.

__Inertia__ of $A\in\mathcal{H}_n$ is the ordered triple

$$\mathop{\mathrm{in}}(A)=(\pi(A),\nu(A),\zeta(A)),$$

where $\pi(A)$ is the number of positive eigenvalues of $A$,  $\nu(A)$ is the number of negative eigenvalues of $A$,
and $\zeta(A)$ is the number of zero eigenvalues of $A$.

__Gram matrix__ of a set of vectors $x_1,x_2,\ldots,x_k\in\mathbb{C}^{n}$ is the matrix $G$ with entries $G_{ij}=x_i^*x_j$.
"""

# ╔═╡ 01e2e27b-506e-439a-ba9a-bef034328482
md"""
## Facts

Assume $A$ is Hermitian and $x\in\mathbb{C}^n$ is nonzero. Then

1. Real symmetric matrix is Hermitian, and real Hermitian matrix is symmetric.

2. Hermitian and real symmetric matrices are normal.

3.  $A+A^*$, $A^*A$, and $AA^*$ are Hermitian.

4. If $A$ is nonsingular, $A^{-1}$ is Hermitian.

5. Main diagonal entries of $A$ are real.

6. Matrix $T$ from the Schur decomposition of $A$ is Hermitian. Consequently:

    *  $T$ is diagonal and real, and has eigenvalues of $A$ on diagonal,
    * matrix $Q$ of the Schur decomposition is the unitary matrix of eigenvectors,
    * all eigenvalues of $A$ are semisimple and $A$ is nondefective,
    * eigenvectors corresponding to distinct eigenvalues are orthogonal.

7. __Spectral Theorem.__ To summarize:

    * if $A\in\mathcal{H}_n$, there is a unitary matrix $U$ and real diagonal matrix $\Lambda$ such that $A=U\Lambda U^*$. The diagonal entries of $\Lambda$ are the eigenvalues of $A$, and the columns of $U$ are the corresponding eigenvectors.
    * if $A\in\mathcal{S}_n$, the same holds with orthogonal matrix $U$, $A=U\Lambda U^T$.
    * if $A\in\mathcal{H}_n$ with eigenpairs $(\lambda_i,u_i)$, then $A=\lambda_1u_1u_1^*+\lambda_2 u_2u_2^* +\cdots +\lambda_n u_nu_n^*$.

    * similarly, if  $A\in\mathcal{S}_n$, then $A=\lambda_1u_1u_1^T+\lambda_2 u_2u_2^T +\cdots +\lambda_n u_nu_n^T$.

8. Since all eigenvalues of $A$ are real, they can be ordered: $\lambda_1\geq \lambda_2\geq \cdots \geq \lambda_n$.

9. __Rayleigh-Ritz Theorem.__ It holds:

$$\begin{aligned}
\lambda_n &\leq \frac{x^*Ax}{x^*x} \leq \lambda_1, \\
\lambda_1&=\max_x\frac{x^*Ax}{x^*x} =\max_{\|x\|_2=1} x^*Ax, \\
\lambda_n&=\min_x\frac{x^*Ax}{x^*x} =\min_{\|x\|_2=1} x^*Ax.
\end{aligned}$$

10. __Courant-Fischer Theorem.__ It holds:

$$\begin{aligned}
\lambda_k& =\max_{\dim(W)=k}\min_{x\in W} \frac{x^*Ax}{x^*x}\\
& =\min_{\dim(W)=n-k+1}\max_{x\in W} \frac{x^*Ax}{x^*x}.
\end{aligned}$$

11. __Cauchy Interlace Theorem.__ Let $B$ be the principal submatrix of $A$ obtained by deleting the $i$-th row and the $i$-th column of $A$. Let $\mu_1\geq \mu_2\geq \cdots \geq \mu_{n-1}$ be the eignvalues of $B$. Then

$$\lambda_1\geq \mu_1\geq \lambda_2\geq \mu_2\geq \lambda_3\geq\cdots\geq \lambda_{n-1}\geq\mu_{n-1}\geq\lambda_n.$$

12. __Weyl Inequalities.__ For $A,B\in\mathcal{H}_n$, it holds:

$$\begin{aligned}
\lambda_{j+k-1}(A+B)& \leq \lambda_j(A)+\lambda_k(B), & \textrm{for} \ j+k\leq n+1,\\
\lambda_{j+k-n}(A+B)& \geq \lambda_j(A)+\lambda_k(B), & \textrm{for} \ j+k\geq n+1,\\
\lambda_j(A)+\lambda_n(B) & \leq \lambda_j(A+B) \leq \lambda_j(A)+\lambda_1(B), & \textrm{for} \ j=1,2,\ldots,n.
\end{aligned}$$

13.  $\pi(A)+\nu(A)+\zeta(A)=n$.

14.  $\mathop{\mathrm{rank}}(A)=\pi(A)+\nu(A)$.

15. If $A$ is nonsingular, $\mathop{\mathrm{in}}(A)=\mathop{\mathrm{in}}(A^{-1})$.

16. If $A,B \in\mathcal{H}_n$ are similar,  $\mathop{\mathrm{in}}(A)=\mathop{\mathrm{in}}(B)$.

17. __Sylvester's Law of Inertia.__ $A,B \in\mathcal{H}_n$ are congruent if and only if $\mathop{\mathrm{in}}(A)=\mathop{\mathrm{in}}(B)$.

18. __Subadditivity of Inertia.__ For $A,B \in\mathcal{H}_n$,

$$\pi(A+B)\leq \pi(A)+\pi(B), \qquad \nu(A+B)\leq \nu(A)+\nu(B).$$

19. Gram matrix is Hermitian.
"""

# ╔═╡ 40c67883-7c5a-4033-96df-090c58b72f65
md"""
## Examples
### Hermitian matrix
"""

# ╔═╡ 43ee6302-cd9e-472e-8c25-9c079d90bd1c
begin
	# Generating Hermitian matrix
	Random.seed!(432)
	nₕ=6
	Aₕ=rand(ComplexF64,nₕ,nₕ)
	Aₕ=Aₕ+adjoint(Aₕ)
end

# ╔═╡ f2b96941-4a21-4c1a-8c91-cd97786ff285
ishermitian(Aₕ)

# ╔═╡ 7a4e98d5-40ba-4ff1-9644-df99163e250d
# Diagonal entries
Diagonal(Aₕ)

# ╔═╡ 669c86fe-c5ac-4b6e-ac38-be628d2d38f3
# Schur decomposition
schur(Aₕ)

# ╔═╡ c19945f2-1fcc-4a33-acdf-615cc9621b12
λₕ,Uₕ=eigen(Aₕ)

# ╔═╡ e98a08fc-58af-4e07-ad77-80492f28490e
# Spectral theorem
norm(Aₕ-Uₕ*Diagonal(λₕ)*Uₕ')

# ╔═╡ 96f6750c-422d-4c2f-9719-32a3c8042181
# Spectral theorem
Aₕ-sum([λₕ[i]*Uₕ[:,i]*Uₕ[:,i]' for i=1:nₕ])

# ╔═╡ 189b6d8f-a080-473c-a9a7-0507707707c8
λₕ

# ╔═╡ a6339443-0249-4b12-b94b-1f4b2b1b6d10
begin
	# Cauchy Interlace Theorem (repeat several times)
	@show i=rand(1:nₕ)
	eigvals(Aₕ[[1:i-1;i+1:end],[1:i-1;i+1:end]])
end

# ╔═╡ b42042c0-5592-4793-b3ff-72ff269dcf6e
# Inertia
inertia(A)=[sum(eigvals(A).>0), sum(eigvals(A).<0), sum(eigvals(A).==0)]

# ╔═╡ 73eb2ab0-92f7-4d74-8960-0b11e8773f4b
inertia(Aₕ)

# ╔═╡ b6adec3b-d812-4268-9d13-9e82dcf30f22
begin
	# Similar matrices
	Mₕ=rand(ComplexF64,nₕ,nₕ)
	eigvals(Mₕ*Aₕ*inv(Mₕ))
end

# ╔═╡ 3181a27c-8360-4a08-bcb0-f6669a9ce062
md"""
This did not work numerically due to rounding errors!
"""

# ╔═╡ e439d2b6-b581-4b87-9f07-ee0e93338a5f
# Congruent matrices - this does not work either, without some preparation
inertia(A)==inertia(Mₕ'*Aₕ*Mₕ)

# ╔═╡ e7a79e8d-d7d7-4280-88b7-58d6d2f228b8
# However,
inertia(Aₕ)==inertia(Hermitian(Mₕ'*Aₕ*Mₕ))

# ╔═╡ 1549b463-2849-4994-a2a3-45c71efda75e
# or,
inertia((Mₕ'*Aₕ*Mₕ+(Mₕ'*Aₕ*Mₕ)')/2)

# ╔═╡ de537b01-876b-40a3-a941-4b5d8bcd2e0f
begin
	# Weyl's Inequalities
	B=rand(ComplexF64,nₕ,nₕ)
	Bₕ=(B+B')/10
	μₕ=eigvals(Bₕ)
	γₕ=eigvals(Aₕ+Bₕ)
	μₕ,γₕ
end

# ╔═╡ 5e49f18c-d633-469e-af8a-1e3877ebf07a
begin
	# Theorem uses different order!
	j=rand(1:nₕ)
	k=rand(1:nₕ)
	@show j,k
	if j+k<=n+1
	    @show sort(γₕ,rev=true)[j+k-1], sort(λₕ,rev=true)[j]+sort(μₕ,rev=true)[k]
	end
	if j+k>=n+1
	    @show sort(γₕ,rev=true)[j+k-n], sort(λₕ,rev=true)[j]+sort(μₕ,rev=true)[k]
	end
end

# ╔═╡ 1866e2e0-8766-4ebd-b045-1f7dcae30b24
sort(λₕ,rev=true)

# ╔═╡ 80c270bc-e313-4b78-94d0-2fcd9aa29a89
md"""
### Real symmetric matrix
"""

# ╔═╡ da3207e5-1cdb-474f-a5f0-11b7c7ea5640
begin
	# Generating real symmetric matrix
	Random.seed!(531)
	nₛ=5
	Aₛ=rand(-9:9,nₛ,nₛ)
	Aₛ=Aₛ+Aₛ'
end

# ╔═╡ e44e5b0c-83e0-4dec-be0d-c5d9bebc7d55
issymmetric(Aₛ)

# ╔═╡ f804922d-b411-4d6b-a973-0e9486756279
schur(Aₛ)

# ╔═╡ 157fd857-d4c8-4433-9e41-9aba559ce0ea
cond(A)

# ╔═╡ dc690338-6ba8-463d-b35b-e57b4e8b0159
inertia(Aₛ)

# ╔═╡ f8334f36-4e46-423f-8028-8a51caff6fc4
begin
	Cₛ=rand(nₛ,nₛ)
	inertia(Cₛ'*Aₛ*Cₛ)
end

# ╔═╡ 53637040-356d-4310-970e-6a9b9c1260c9
md"""
# Positive definite matrices

These matrices are an important subset of Hermitian or real symmteric matrices.

## Definitions

Matrix $A\in\mathcal{H}_n$ is __positive definite__ (PD) if $x^*Ax>0$ for all nonzero $x\in\mathbb{C}^n$.

Matrix $A\in\mathcal{H}_n$ is __positive semidefinite__ (PSD) if $x^*Ax\geq 0$ for all nonzero $x\in\mathbb{C}^n$.

## Facts

1.  $A\in\mathcal{S}_n$ is PD if $x^TAx>0$ for all nonzero $x\in \mathbb{R}^n$, and is PSD if $x^TAx\geq 0$ for all $x\in \mathbb{R}^n$.

2. If $A,B\in \mathrm{PSD}_n$, then $A+B\in \mathrm{PSD}_n$. If, in addition, $A\in \mathrm{PD}_n$, then $A+B\in \mathrm{PD}_n$.

3. If $A\in \mathrm{PD}_n$, then $\mathop{\mathrm{tr}}(A)>0$ and $\det(A)>0$.

3. If $A\in \mathrm{PSD}_n$, then $\mathop{\mathrm{tr}}(A)\geq 0$ and $\det(A)\geq 0$.

4. Any principal submatrix of a PD matrix is PD. Any principal submatrix of a PSD matrix is PSD. Consequently, all minors are positive or nonnegative, respectively.

5.  $A\in\mathcal{H}_n$ is PD iff __every leading__ principal minor of $A$ is positive. $A\in\mathcal{H}_n$ is PSD iff __every__ principal minor is nonnegative.

6. For $A\in \mathrm{PSD}_n$, there exists unique PSD $k$-th __root__, $A^{1/k}=U\Lambda^{1/k} U^*$.

7. __Cholesky Factorization.__ $A\in\mathcal{H}_n$ if PD iff there is an invertible lower triangular matrix $L$ with positive diagonal entries such that $A=LL^*$.

8. Gram matrix is PSD. If the vectors are linearly independent, Gram matrix is PD.
"""

# ╔═╡ 5082d501-8e0f-45ec-8927-ad1a0be51ae3
md"""
## Examples
### Positive definite matrix
"""

# ╔═╡ e00702a4-2ad5-4dee-8d0c-7d15c4ffdeb8
begin
	# Generating positive definite matrix as a Gram matrix
	nₚ=5
	Aₚ=rand(ComplexF64,nₚ,nₚ)
	Aₚ=Aₚ*Aₚ'
end

# ╔═╡ fb844ff8-1696-45a8-97a9-a08ae0ace55b
ishermitian(Aₚ)

# ╔═╡ c6d47f33-1c78-4568-a0f6-e9d1deb669ca
eigvals(Aₚ)

# ╔═╡ ada85c66-6268-4664-a458-4213d7d7111f
# Positivity of principal leading minors
[det(Aₚ[1:i,1:i]) for i=1:nₚ]

# ╔═╡ a43a9532-f865-4166-b0c6-124045030e64
begin
	# Matrix function - square root
	λₚ,Uₚ=eigen(Aₚ)
	A2=Uₚ*√Diagonal(λₚ)*Uₚ'
	Aₚ-A2^2
end

# ╔═╡ 8a4c37e0-a442-46da-929c-24cdc8280313
# Cholesky factorization - lower triangular factor
Lₚ=cholesky(Aₚ).L

# ╔═╡ 1ad5826e-f849-4cef-82a1-2bd97cf38731
norm(Aₚ-Lₚ*Lₚ')

# ╔═╡ 8b3db768-e5d1-49cf-8320-7cb36a400ce1
md"
### Positive semidefinite matrix
"

# ╔═╡ 18457042-90d2-4180-a999-8d3454134528
begin
	# Generate positive semidefinite matrix as a Gram matrix, try it several times
	A₅=rand(-9:9,6,4)
	A₅=A₅*A₅'
end

# ╔═╡ 853b9c26-a15b-4109-81ee-5ce6a2dd42bd
# There are rounding errors!
eigvals(A₅)

# ╔═╡ ce5c03dd-f8a8-453d-a3b8-d618aa33be9e
rank(A₅)

# ╔═╡ 710b83d0-8d8c-11eb-39c2-e97b1880a16c
md"
How does function `rank()` work?
"

# ╔═╡ 854dfa30-8d8c-11eb-0d82-7d17e7d40d98
@which rank(A₅)

# ╔═╡ 5a229e53-a1e0-4d0f-b479-22bbd4dfaa06
s₅=svdvals(A₅)

# ╔═╡ 7056a7af-f3cf-4fa6-aa4e-b891b6444aa6
s₅[1]*eps()

# ╔═╡ f4b34f27-082b-4bbe-8ee0-c1066de09a9b
# Cholesky factorization can fail
cholesky(A₅)

# ╔═╡ 0528ccb5-04a8-4b8e-a0c8-d845ca716c64
md"
### Covariance and correlation matrices

[Covariance and correlation matrices](https://en.wikipedia.org/wiki/Covariance_matrix) are PSD.

Correlation matrix is diagonally scaled covariance matrix.
"

# ╔═╡ 9daeb277-5cbb-406e-b42f-c21ddaf4366f
begin
	Random.seed!(651)
	y=rand(10,5)
end

# ╔═╡ 69346220-2310-4015-b29b-5b7186791fa8
begin
	using Statistics
	Cov=cov(y)
end

# ╔═╡ 688efeb3-07da-4d75-a35b-a79f9440da76
# Covariance matrix is a Gram matrix
(y.-mean(y,dims=1))'*(y.-mean(y,dims=1))/(size(y,1)-1)-Cov

# ╔═╡ f28e0e9a-23c5-4114-a853-2da285a24641
# Correlation matrix
Cor=cor(y)

# ╔═╡ 867926b1-3d3e-4c7f-95cf-54ed6b54c112
# Diagonal scaling of ovariance matrix
D=inv(√Diagonal(Cov))

# ╔═╡ c5fe8288-4650-4a38-a322-38780978095d
D*Cov*D

# ╔═╡ f6be3df8-6368-4535-8085-537a80dd3ae0
eigvals(Cov)

# ╔═╡ 808904ca-b1d9-4eca-ae40-bcd53451fa83
eigvals(Cor)

# ╔═╡ 34001f70-e2db-4529-a315-caf38ecb87b0
Covₜ=cov(y')

# ╔═╡ 80a7d977-99e7-4d51-9cc9-4e6b7c410a6a
eigvals(Covₜ)

# ╔═╡ c74eb635-2948-4db3-bd53-c7d7a5b2a454
inertia(Covₜ)

# ╔═╡ cc0422bf-2e9b-463a-9b1a-468b09bb370e
rank(Covₜ)

# ╔═╡ 3e29a10b-0e07-4ddf-ad6d-848100804262
md"
# Eigenvalues of random matrices
"

# ╔═╡ 846047e6-56c7-4ce5-b7f6-6d037c9918bb
md"""
k = $(@bind kₘ Slider(10:30,show_value=true))

n = $(@bind nₘ Slider(10:30,show_value=true))

Matrix type = $(@bind mt Select(["Uniform", "Normal", "Uniform symmetric", "Normal Symmetric"]))
"""

# ╔═╡ c36a856f-a726-4ccb-bee9-691f01147b97
begin
	E=Array{Any}(undef,nₘ,kₘ)
	for i=1:kₘ
		if mt=="Uniform"
			# Unsymmetric uniform distribution
	    	A=rand(nₘ,nₘ)
		elseif mt=="Normal"
			# Unsymmetric normal distribution
			A=randn(nₘ,nₘ)
		elseif mt=="Uniform symmetric"
			# Symmetric uniform distribution
			A=Symmetric(rand(nₘ,nₘ))
		else
			# Symmetric normal distribution
			A=Symmetric(randn(nₘ,nₘ))
		end
	    E[:,i]=eigvals(A)
	end
	# We need this since plot cannot handle `Any`
	E=map(eltype(E[1,1]),E)
	scatter(E,legend=false)
end

# ╔═╡ Cell order:
# ╟─10910df9-240f-4e6b-8076-1a9be4d66ba1
# ╟─ec1fe200-8d73-11eb-0ba3-e593aa79dab2
# ╟─b87354b9-ee7b-4ce0-bf26-b4e71c6be6a9
# ╟─31a71b9e-1dba-4a69-b87e-9202292eb3ed
# ╟─2f5aee03-3c2e-4920-9d10-c509b06be627
# ╠═973abf88-9125-41c7-8372-61efb2976fc2
# ╟─37d24090-8d78-11eb-059d-737014bb6588
# ╠═184bc4ec-c566-4f21-99b3-b6e8cd328b47
# ╠═5d0aac9e-4f48-4f09-956f-5033a1fb8f11
# ╠═73ece9ec-97dd-46bc-8e17-10ce7bfc18bc
# ╠═ce83f04c-e915-4599-beeb-d4a0a681bc17
# ╠═ff5044c0-8d76-11eb-3951-6bb4b5caa7d0
# ╠═55171a69-3557-4ba5-8a4a-813f38c03dd3
# ╠═6df42190-a29e-4c68-8684-07d8d0c9d200
# ╟─b94f7fcd-c6f4-4f4f-924a-20f00e82d94e
# ╠═1b90d50a-bf97-4b1a-9580-d729b03c914e
# ╠═15b60cd5-eb30-43bb-bbcf-dd533e29c35f
# ╟─3619fca4-1374-4adc-b29c-8058b2e147e0
# ╠═d46018c0-a5f8-4dd0-8f3d-73322be792e9
# ╠═d387b26f-bec1-47a9-b1c6-f20c43f0bde8
# ╠═703da5c4-cd09-4832-b64e-5fda75099db5
# ╠═e26f113c-dd98-4e97-b495-581de1b431af
# ╠═62779bfc-ea5d-4eec-94ad-2641e11fc44d
# ╠═bdf079a6-beee-44f5-8200-613e54b3a8dc
# ╠═e1f5027c-1080-4c95-81cc-a75dde641b88
# ╠═c71a7760-b4c6-404e-b5ef-f4879cb787cd
# ╟─343c8f45-6226-4361-a149-5313a47c355c
# ╠═41a05c50-f6df-4b93-bd8a-50412829105e
# ╠═cd8fb666-d2f2-423c-ac39-405c842e97ba
# ╟─7aa7ac89-837c-4ee8-b9f1-8ad0256f25e1
# ╠═c3724c4f-142b-4863-8c04-4a0c15411f40
# ╠═792f632e-fdc2-4016-baaa-b366704f02fd
# ╠═a8664a40-8d78-11eb-325e-41c6497be282
# ╠═1d42245a-2ed6-4cc7-8cbd-2b6aa4e64861
# ╠═f7aae752-70c4-4a10-b96e-206c13a6ced1
# ╠═60525ba5-2c76-406e-8eb3-63ced2e00bdd
# ╟─189544ef-122c-4e30-bfd6-cfc126fc0320
# ╠═da340240-f845-4591-94c3-9ffce6a33d8a
# ╠═e72bf173-7e55-4476-91ff-789847dc53cf
# ╠═da4a5d70-92d8-11eb-0635-11c530ff7d21
# ╠═ec7fc5c0-92d8-11eb-16df-af1133ccb42a
# ╠═29de5fd0-92d9-11eb-09c0-0356001caaef
# ╠═d43480a3-9247-48a4-8a3b-e006ab58539e
# ╠═5995cd31-a75d-4718-8808-4ada01920344
# ╠═57e2742e-8d79-11eb-0a33-f7250ac02840
# ╠═5e069120-8d79-11eb-3c00-0b5249d1f2dc
# ╠═63167eee-8d79-11eb-29c2-2febbd7cbd93
# ╟─e051d2eb-3972-4c97-b7eb-8a18629f6171
# ╠═22f0a627-beb2-4b06-aa94-ee1e90fe8a6c
# ╠═da6baba0-92d9-11eb-1603-7961ad1f49c8
# ╟─c6ab93ab-b336-4bb4-9dda-b725108f7827
# ╠═ffcb919f-9754-4650-92f1-7dd163c5df25
# ╠═32ce58b9-783b-4739-9001-0fe581b2953d
# ╠═a7e8d711-57a6-477a-a641-b929789dbadf
# ╠═5781508c-388a-4df0-a6f9-0fe312a81335
# ╟─90b4341d-7964-4e41-a268-8a047a058286
# ╠═16795c64-2994-4fca-89f3-f035fb6a2e45
# ╠═60d351a5-0ef2-47cd-a5a5-5ecf18c4e898
# ╠═76ffa073-59b8-4118-aa32-522f4ff05979
# ╠═f85ce269-7555-48ed-a2fd-1193d3f33f40
# ╠═e6544939-c8c2-464c-83a0-a6bf1751efa1
# ╠═b559d56b-02bb-47d3-8f0f-8dae09d918d8
# ╠═5e0adf59-29d3-47ac-8f3f-52929cc4e1c1
# ╠═83ba2d1b-66de-41bf-8ded-7070bed15587
# ╠═b0255c92-3242-442a-a670-1bfef38e955c
# ╠═c6116576-5b25-45b2-bc71-a889f302b5ac
# ╠═3de2dc38-ba95-479e-89bc-77cb520daf2f
# ╠═ad32c51c-7cbb-42f9-abc8-561828df1379
# ╟─a7999464-7060-4345-b396-e2e15fcbcaac
# ╟─01e2e27b-506e-439a-ba9a-bef034328482
# ╟─40c67883-7c5a-4033-96df-090c58b72f65
# ╠═43ee6302-cd9e-472e-8c25-9c079d90bd1c
# ╠═f2b96941-4a21-4c1a-8c91-cd97786ff285
# ╠═7a4e98d5-40ba-4ff1-9644-df99163e250d
# ╠═669c86fe-c5ac-4b6e-ac38-be628d2d38f3
# ╠═c19945f2-1fcc-4a33-acdf-615cc9621b12
# ╠═e98a08fc-58af-4e07-ad77-80492f28490e
# ╠═96f6750c-422d-4c2f-9719-32a3c8042181
# ╠═189b6d8f-a080-473c-a9a7-0507707707c8
# ╠═a6339443-0249-4b12-b94b-1f4b2b1b6d10
# ╠═b42042c0-5592-4793-b3ff-72ff269dcf6e
# ╠═73eb2ab0-92f7-4d74-8960-0b11e8773f4b
# ╠═b6adec3b-d812-4268-9d13-9e82dcf30f22
# ╟─3181a27c-8360-4a08-bcb0-f6669a9ce062
# ╠═e439d2b6-b581-4b87-9f07-ee0e93338a5f
# ╠═e7a79e8d-d7d7-4280-88b7-58d6d2f228b8
# ╠═1549b463-2849-4994-a2a3-45c71efda75e
# ╠═de537b01-876b-40a3-a941-4b5d8bcd2e0f
# ╠═5e49f18c-d633-469e-af8a-1e3877ebf07a
# ╠═1866e2e0-8766-4ebd-b045-1f7dcae30b24
# ╟─80c270bc-e313-4b78-94d0-2fcd9aa29a89
# ╠═da3207e5-1cdb-474f-a5f0-11b7c7ea5640
# ╠═e44e5b0c-83e0-4dec-be0d-c5d9bebc7d55
# ╠═f804922d-b411-4d6b-a973-0e9486756279
# ╠═157fd857-d4c8-4433-9e41-9aba559ce0ea
# ╠═dc690338-6ba8-463d-b35b-e57b4e8b0159
# ╠═f8334f36-4e46-423f-8028-8a51caff6fc4
# ╟─53637040-356d-4310-970e-6a9b9c1260c9
# ╟─5082d501-8e0f-45ec-8927-ad1a0be51ae3
# ╠═e00702a4-2ad5-4dee-8d0c-7d15c4ffdeb8
# ╠═fb844ff8-1696-45a8-97a9-a08ae0ace55b
# ╠═c6d47f33-1c78-4568-a0f6-e9d1deb669ca
# ╠═ada85c66-6268-4664-a458-4213d7d7111f
# ╠═a43a9532-f865-4166-b0c6-124045030e64
# ╠═8a4c37e0-a442-46da-929c-24cdc8280313
# ╠═1ad5826e-f849-4cef-82a1-2bd97cf38731
# ╟─8b3db768-e5d1-49cf-8320-7cb36a400ce1
# ╠═18457042-90d2-4180-a999-8d3454134528
# ╠═853b9c26-a15b-4109-81ee-5ce6a2dd42bd
# ╠═ce5c03dd-f8a8-453d-a3b8-d618aa33be9e
# ╟─710b83d0-8d8c-11eb-39c2-e97b1880a16c
# ╠═854dfa30-8d8c-11eb-0d82-7d17e7d40d98
# ╠═5a229e53-a1e0-4d0f-b479-22bbd4dfaa06
# ╠═7056a7af-f3cf-4fa6-aa4e-b891b6444aa6
# ╠═f4b34f27-082b-4bbe-8ee0-c1066de09a9b
# ╟─0528ccb5-04a8-4b8e-a0c8-d845ca716c64
# ╠═9daeb277-5cbb-406e-b42f-c21ddaf4366f
# ╠═69346220-2310-4015-b29b-5b7186791fa8
# ╠═688efeb3-07da-4d75-a35b-a79f9440da76
# ╠═f28e0e9a-23c5-4114-a853-2da285a24641
# ╠═867926b1-3d3e-4c7f-95cf-54ed6b54c112
# ╠═c5fe8288-4650-4a38-a322-38780978095d
# ╠═f6be3df8-6368-4535-8085-537a80dd3ae0
# ╠═808904ca-b1d9-4eca-ae40-bcd53451fa83
# ╠═34001f70-e2db-4529-a315-caf38ecb87b0
# ╠═80a7d977-99e7-4d51-9cc9-4e6b7c410a6a
# ╠═c74eb635-2948-4db3-bd53-c7d7a5b2a454
# ╠═cc0422bf-2e9b-463a-9b1a-468b09bb370e
# ╟─3e29a10b-0e07-4ddf-ad6d-848100804262
# ╠═bb8416e3-b682-4733-af6d-36fd48f0a7b8
# ╟─846047e6-56c7-4ce5-b7f6-6d037c9918bb
# ╠═c36a856f-a726-4ccb-bee9-691f01147b97

max_stars_repo_path: src/Module B - Eigenvalue and Singular Value Decompositions/L3a Eigenvalue Decomposition - Definitions and Facts.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
p=50;setprecision(p); setrounding(BigFloat, RoundUp);
L=BigFloat("10"); n=BigInt(100); iN=convert(Int64,n);
dx=BigFloat("2")*L/(n-1); x=convert(Array{Float64},-L:dx:L)
I=zeros(BigFloat, n, n);for i=1:n I[i,i]=BigFloat("1") end

f=sin.(x)+sin.(x./BigFloat("2"))
fi=sin.(xi)+sin.(xi./BigFloat("2"))
f+=convert(Array{BigFloat},(rand(n).-1/2)./1)

ff=3.2
typeof(ff)
ff16=convert(Float16,ff)
ff32=convert(Float32,ff)
ff16-ff32

scatter(x,f)
plot!(x,f)

i_dd=Array{Int64}(undef,0)
for i=2:length(d)-1
    if (f[i+1]-f[i])*(f[i-1]-f[i])>0
        push!(i_dd,i)
    end
end
scatter!(x[i_dd],f[i_dd])

V=TaylorBigF.calc_Vandermonde(x,n,n+BigInt(length(i_dd)))
xi=(x.+dx/BigFloat("2"))[1:end-1]
Vi=TaylorBigF.calc_Vandermonde(xi,n-1,n+BigInt(length(i_dd)))
ddV=TaylorBigF.calc_Vandermonde_dot_dot(x,n,n+BigInt(length(i_dd)))

Eq1=ddV[i_dd,:]
q1=zeros(BigFloat,length(i_dd))

VV=[V;Eq1]
q_all=[f;q1]
a=VV\q_all
plot(a)

r_root=Array{Float64}(undef,0); for i=1:Int64(n-1) push!(r_root,abs(a[i])^(1/(n-1))) end; r_rootF=convert(Array{Float64},r_root);
plot(r_rootF)
error1=maximum(abs.(f-V*a))
plot(x,f)
pred=V*a
scatter!(x[20:80],pred[20:80])
Vi=TaylorBigF.calc_Vandermonde(xi,n-1,BigInt(150));
fi_comp=Vi*a
error2=maximum(abs.(fi_comp-fi))
scatter!(xi,Vi*a)


text:
using Test
using GradientBoost.ML
using GradientBoost.LossFunctions

instances = [
  1.0 1.0;
  1.0 8.0;
  1.0 10.0
]
labels = [
  0.0;
  1.0;
  1.0;
]

@testset "Machine Learning API" begin
  @testset "not implemented functions throw an error" begin
    gbl = GBLearner(GBDT(), :regression)
    instances = 1
    labels = 1

    @test_throws UndefVarError fit!(gbl, instances, labels)
    @test_throws ErrorException predict!(gbl, instances)
  end

  @testset "fit! on Float64 arrays works" begin
    gbl = GBLearner(GBDT(), :regression)
    @test gbl.model === nothing
    GradientBoost.ML.fit!(gbl, instances, labels)
    @test gbl.model != nothing
  end

  @testset "predict! on Float64 arrays works" begin
    gbl = GBLearner(GBDT(;loss_function=BinomialDeviance()), :class)
    GradientBoost.ML.fit!(gbl, instances, labels)
    predictions = GradientBoost.ML.predict!(gbl, instances)
    @test predictions isa Vector{Float64}
  end

  @testset "postprocess_pred works" begin
    predictions = [-Inf, 0.0, Inf]
    expected = [0.0, 1.0, 1.0]
    actual = ML.postprocess_pred(:class, BinomialDeviance(), predictions)
    @test actual ≈ expected

    predictions = [-Inf, 0.0, Inf]
    actual = ML.postprocess_pred(:class_prob, BinomialDeviance(), predictions)
    @test all(i -> (0 <= i <= 1), actual)

    predictions = [-Inf, 0.0, Inf]
    expected = predictions
    actual = ML.postprocess_pred(:regression, LeastSquares(), predictions)
    @test actual == expected

    predictions = [-Inf, 0.0, Inf]
    @test_throws ErrorException ML.postprocess_pred(:class, LeastSquares(), predictions)
  end

  @testset "logistic works" begin
    x = [-Inf, -1, 0, 1, Inf]
    expected = [0.0, 0.0, 1.0, 1.0, 1.0]

    actual = round.(Int, ML.logistic(x), RoundNearestTiesUp)
    @test actual ≈ expected
  end
end


text:
using Pkg
Pkg.activate(".")

using PackageCompiler

create_sysimage(; sysimage_path="Base16.so", incremental=false)


text:
#global tophap_ccf_mask_scale_factor=1.6

global max_spectra_to_use = 200
if max_spectra_to_use < 200
   @warn "param.in setting max_spectra_to_use to " * string(max_spectra_to_use)
end
global fits_target_str
if fits_target_str == "101501"
   global espresso_mask_filename = "G9.espresso.mas"
   global ccf_mid_velocity = -5.0e3
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "10700"
   global espresso_mask_filename = "G8.espresso.mas"
   global ccf_mid_velocity = -16640.0
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "34411"  #G1?
   global espresso_mask_filename = "G2.espresso.mas"
   global ccf_mid_velocity = 66500.0
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "26965" # K0?
   global espresso_mask_filename = "G9.espresso.mas"
   global ccf_mid_velocity = -40320.0
   global df_files
   global df_files_use = df_files |>
      @filter( _.target == fits_target_str ) |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "Sim1"
   global espresso_mask_filename = "G2.espresso.mas"
   global ccf_mid_velocity = 67000.0
   global df_files
   global df_files_use = df_files |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "Sim2"
   global espresso_mask_filename = "G2.espresso.mas"
   global ccf_mid_velocity = 67000.0
   global df_files
   global df_files_use = df_files |>
      @take(max_spectra_to_use) |>
      DataFrame
end

if fits_target_str == "Sim3"
   global espresso_mask_filename = "G2.espresso.mas"
   global ccf_mid_velocity = 67000.0
   global df_files
   global df_files_use = df_files |>
      @take(max_spectra_to_use) |>
      DataFrame
end


text:
module GaussianProcesses

using Optim, PDMats, ElasticPDMats, ElasticArrays, Distances, FastGaussQuadrature, RecipesBase, Distributions
using StaticArrays
using StatsFuns, SpecialFunctions

using LinearAlgebra, Printf, Random, Statistics
import Statistics: mean

# Functions that should be available to package
# users should be explicitly exported here
export GPBase, GP, GPE, GPMC, predict_f, predict_y, Kernel, Likelihood, CompositeKernel, SumKernel, ProdKernel, Masked, FixedKernel, fix, Noise, Const, SE, SEIso, SEArd, Periodic, Poly, RQ, RQIso, RQArd, Lin, LinIso, LinArd, Matern, Mat12Iso, Mat12Ard, Mat32Iso, Mat32Ard, Mat52Iso, Mat52Ard, #kernel functions
    MeanZero, MeanConst, MeanLin, MeanPoly, SumMean, ProdMean, #mean functions
    GaussLik, BernLik, ExpLik, StuTLik, PoisLik, BinLik,       #likelihood functions
    mcmc, optimize!,                                           #inference functions
    set_priors!,set_params!, update_target!

const MatF64 = AbstractMatrix{Float64}
const VecF64 = AbstractVector{Float64}

const φ = normpdf
const Φ = normcdf
const invΦ = norminvcdf

# all package code should be included here
include("means/means.jl")
include("kernels/kernels.jl")
include("likelihoods/likelihoods.jl")
include("common.jl")
include("utils.jl")
include("chol_utils.jl")
include("GP.jl")
include("GPE.jl")
include("GPEelastic.jl")
include("GPMC.jl")
include("mcmc.jl")
include("optimize.jl")
include("plot.jl")

# ScikitLearnBase, which is a skeleton package.
include("ScikitLearn.jl")

end # module


text:
module MarkovTools

export
    #functions
    count_fixed_points,
    count_cycles,
    readbitstream,

    min,
    min2,
    Spritz_Encrypt,
    Spritz_Decrypt

include("permutations.jl")
include("prng/input.jl")
include("spritz/spritz.jl")
include("monotonicites/monotonicities.jl")

end # module

max_stars_repo_name: PawelPamula/MarkovTools

max_issues_repo_name: PawelPamula/MarkovTools

max_forks_repo_name: PawelPamula/MarkovTools


text:
function PhaseShift(x,angle)

	nx = length(x)
	nf = 2*nx
	nw = int(nf/2) + 1
	X = fft([x;zeros(typeof(x[1]),nf-nx)])
	for iw=1:nw
		X[iw] *= exp(-1im*angle)
	end
	# symmetries
	for iw=nw+1:nf
		X[iw] = conj(X[nf-iw+2])
	end
	x = real(ifft(X,1))
	return x[1:nx]

end


text:
# descriptor

mutable struct DropoutDesc
  ptr::Ptr{Nothing}
  states::CuVector{UInt8}
end

Base.unsafe_convert(::Type{Ptr{Nothing}}, dd::DropoutDesc) = dd.ptr

function DropoutDesc(ρ::Real; seed::Integer=0)
  d = [C_NULL]
  s = Csize_t[0]
  cudnnCreateDropoutDescriptor(d)
  cudnnDropoutGetStatesSize(handle(), s)
  states = CuArray{UInt8}(undef, s[]) # TODO: can we drop this when ρ=0?
  desc = DropoutDesc(d[], states)
  cudnnSetDropoutDescriptor(desc, handle(), ρ, states, length(states), seed)
  finalizer(desc) do x
    cudnnDestroyDropoutDescriptor(x)
  end
  return desc
end


text:
export SigmoidalProgram, LinearSP

## problem types
abstract type SigmoidalProgram end

# LinearSP represents the problem
# maximize    sum_i f_i(x_i)
# subject to  A x <= b
#             C x == d
struct LinearSP <: SigmoidalProgram
    fs::Array{Function,1}
    dfs::Array{Function,1}
    z::Array{Float64,1}
    A::Array{Float64,2}
    b::Array{Float64,1}
    C::Array{Float64,2}
    d::Array{Float64,1}
    function LinearSP(fs,dfs,z,A,b,C,d)
        nconstr, nvar = size(A)
        @assert(length(b) == nconstr)
        nconstr, nvar = size(C)
        @assert(length(d) == nconstr)
        @assert(length(fs) == nvar)
        @assert(length(dfs) == nvar)
        @assert(length(z) == nvar)
        new(fs,dfs,z,A,b,C,d)
    end
end
# Default constructor when no equality constraints are present
LinearSP(fs, dfs, #::Array{Function,1},dfs::Array{Function,1},
         z::Array{Float64,1},A::Array{Float64,2},b::Array{Float64,1}) =
         LinearSP(fs,dfs,z,A,b,zeros(0,length(fs)),zeros(0))

function addConstraints!(m::Model, x, p::LinearSP)
    # Ax <= b
    nconstr, nvar = size(p.A)
    for i=1:nconstr
        @constraint(m, sum(p.A[i,j]*x[j] for j=1:nvar) <= p.b[i])
    end
    # Cx == d
    nconstr, nvar = size(p.C)
    for i=1:nconstr
        @constraint(m, sum(p.C[i,j]*x[j] for j=1:nvar) == p.d[i])
    end
end


text:
using PtFEM, Test

data = Dict(
# Plane(ndim, nst, nxe, nye, nip, direction=(:x|:y), finite_element(nod, nodof), axisymmetric=false)
# Plane(ndim, nst, nxe, nye, nip, direction=:z, finite_element(nod, nodof), axisymmetric=true)
# Plane(ndim, nst, nre, nze, nip, direction=:r, finite_element(nod, nodof), axisymmetric=false)
  :struc_el => Plane(2, 6, 1, 5, 4, :r, Quadrilateral(8, 3), false),
  :properties => [1.0e5 0.3;],
  :r_coords => [0.0, 0.5],
  :z_coords => [10.0, 8.0, 6.0, 4.0, 2.0, 0.0],
  :support => [
    ( 1, [1 0 1]),
    ( 4, [1 0 1]),
    ( 6, [1 0 1]),
    ( 9, [1 0 1]),
    (11, [1 0 1]),
    (14, [1 0 1]),
    (16, [1 0 1]),
    (19, [1 0 1]),
    (21, [1 0 1]),
    (24, [1 0 1]),
    (26, [0 0 0]),
    (27, [0 0 0]),
    (28, [0 0 0])
    ],
  :loaded_nodes => [
    (3, [1/pi 0.0 0.0])
    ],
  :lth => 1,
  :iflag => 1,
  :chi => 0.0
)

@time m, dis_df, sigma_df = p52(data)

@test m.loads[0:65] ≈ [0.0,0.0675569552650584,-0.06755339718587051,
  0.0675550611183897,-0.0025282635314929514,-0.06754930096906765,
  0.06755376595399615,-0.005063646782360032,-0.06754378259701828,
  0.05743046037889881,-0.057432490517240836,0.057441944750754405,
  -0.005005739743328588,-0.05742749087610757,0.04753054505535036,
  -0.04752501175531152,0.04752761039720702,-0.0024265026287448132,
  -0.04751826130337188,0.04753200512262878,-0.00485859624896776,
  -0.04750446544560081,0.0380108152341777,-0.03801430582540627,
  0.038043595025430185,-0.004598096664133727,-0.0379939151061887,
  0.029138408322084745,-0.029128707630258533,0.029131687887729216,
  -0.002121143506426845,-0.02911866562356955,0.029141015554487087,
  -0.0042482014668140024,-0.029092503165066133,0.021034656820675384,
  -0.021040857210335804,0.021090965865308207,-0.0037823104185696304,
  -0.021005867287750803,0.014004054458491054,-0.013993277627426427,
  0.013994731459463667,-0.0016135944093911055,-0.013978288565264622,
  0.014005550176560485,-0.00323312039147469,-0.01394313064790418,
  0.008132647707849692,-0.008137157612832228,0.008212440365727661,
  -0.0025592520670215577,-0.008085266130752576,0.003764582745644812,
  -0.0037356644208552323,0.003738861252011861,-0.0009009830493624234,
  -0.003731710717854891,0.0037540460087339984,-0.0018087025205919894,
  -0.0036869554668832335,0.000937846981324757,-0.0009620426947358215,
  0.0010519661074415124,-0.0009218802871890891,-0.0008805978030738389]


text:
@testset "Interval" begin
    @test MOI.Interval(MOI.GreaterThan(1.0)) === MOI.Interval(1.0, Inf)
    @test MOI.Interval(MOI.LessThan(2.0)) === MOI.Interval(-Inf, 2.0)
    @test MOI.Interval(MOI.EqualTo(3.0)) === MOI.Interval(3.0, 3.0)

    @test MOI.Interval(MOI.GreaterThan(1.0f0)) === MOI.Interval(1.0f0, Inf32)
    @test MOI.Interval(MOI.LessThan(2.0f0)) === MOI.Interval(-Inf32, 2.0f0)
    @test MOI.Interval(MOI.EqualTo(3.0f0)) === MOI.Interval(3.0f0, 3.0f0)

    @test_throws MethodError MOI.Interval(MOI.GreaterThan(false))
    @test_throws MethodError MOI.Interval(MOI.LessThan(true))
    @test MOI.Interval(MOI.EqualTo(true)) === MOI.Interval(true, true)

    @test_throws MethodError MOI.Interval(MOI.GreaterThan(1))
    @test_throws MethodError MOI.Interval(MOI.LessThan(2))
    @test MOI.Interval(MOI.EqualTo(3)) === MOI.Interval(3, 3)

    @test MOI.Interval(MOI.Interval(1, 2)) == MOI.Interval(1, 2)
end


text:
#module SpatialRegionTools

using StatsBase:rle
using HDF5
using DataStructures
using NearestNeighbors
include("utils.jl")

#export SpatialRegion, makeVocab!, gps2vocab, saveKNearestVocabs,
#       trip2seq, seq2trip, createTrainVal

const UNK = 3

"""
example:

region = SpatialRegion("porto",
                       -8.735152, 40.953673,
                       -8.156309, 41.307945,
                       cellsize, cellsize,
                       50, # minfreq
                       50_000, # maxvocab_size
                       5, # k
                       4) # vocab_start
"""
mutable struct SpatialRegion
    name::String
    ## bounding box
    minlon::Float64
    minlat::Float64
    maxlon::Float64
    maxlat::Float64
    minx::Float64
    miny::Float64
    maxx::Float64
    maxy::Float64
    xstep::Float64
    ystep::Float64
    numx::Int
    numy::Int
    minfreq::Int
    maxvocab_size::Int
    k::Int
    ## the number of hitting
    cellcount
    ## hot cells
    hotcell::Vector{Int}
    ## hot cell kdtree
    hotcell_kdtree
    ## map a hot cell into the vocabulary id
    hotcell2vocab::Dict{Int, Int}
    ## map the vocabulary id into its hot cell
    vocab2hotcell::Dict{Int, Int}
    ## vocab start
    vocab_start::Int
    ## vocabulary size
    vocab_size::Int
    ## whether indices have been built
    built::Bool
    function SpatialRegion(name::String,
                           minlon::Float64,
                           minlat::Float64,
                           maxlon::Float64,
                           maxlat::Float64,
                           xstep::Float64,
                           ystep::Float64,
                           minfreq::Int,
                           maxvocab_size::Int,
                           k::Int,
                           vocab_start::Int)
        minx, miny = lonlat2meters(minlon, minlat)
        maxx, maxy = lonlat2meters(maxlon, maxlat)
        numx = round(maxx - minx, digits=6) / xstep
        numx = convert(Int, ceil(numx))
        numy = round(maxy - miny, digits=6) / ystep
        numy = convert(Int, ceil(numy))
        new(name,
            minlon, minlat, maxlon, maxlat,
            minx, miny, maxx, maxy,
            xstep, ystep,
            numx, numy, minfreq, maxvocab_size, k,
            Accumulator{Int, Int}(),
            Int[],
            Any,
            Dict(),
            Dict(),
            vocab_start,
            vocab_start,
            false)
    end
end

#function gps2cell(region::SpatialRegion, lon::Float64, lat::Float64)
#    lonoffset = round(lon - region.minlon, 6) / region.hstep
#    latoffset = round(lat - region.minlat, 6) / region.vstep
#    lonoffset = convert(Int64, floor(lonoffset))
#    latoffset = convert(Int64, floor(latoffset))
#    latoffset * region.numlon + lonoffset
#end
#
#function cell2gps(region::SpatialRegion, cell::Int)
#    latoffset = div(cell, region.numlon)
#    lonoffset = mod(cell, region.numlon)
#    lat = region.minlat + (latoffset + 0.5) * region.vstep
#    lon = region.minlon + (lonoffset + 0.5) * region.hstep
#    lon, lat
#end

#cell2gps(region, 0)
#gps2cell(region, (126.48349, 45.6468)...)

"""
Web Mercator coordinate to cell id
"""
function coord2cell(region::SpatialRegion, x::Float64, y::Float64)
    xoffset = round(x - region.minx, digits=6) / region.xstep
    yoffset = round(y - region.miny, digits=6) / region.ystep
    xoffset = convert(Int, floor(xoffset))
    yoffset = convert(Int, floor(yoffset))
    yoffset * region.numx + xoffset
end

"""
Cell id to Web Mercator coordinate
"""
function cell2coord(region::SpatialRegion, cell::Int)
    yoffset = div(cell, region.numx)
    xoffset = mod(cell, region.numx)
    y = region.miny + (yoffset + 0.5) * region.ystep
    x = region.minx + (xoffset + 0.5) * region.xstep
    x, y
end

function gps2cell(region::SpatialRegion, lon::Float64, lat::Float64)
    x, y = lonlat2meters(lon, lat)
    coord2cell(region, x, y)
end

function cell2gps(region::SpatialRegion, cell::Int)
    x, y = cell2coord(region, cell)
    meters2lonlat(x, y)
end

#lon, lat = -8.735152, 40.953673
#x, y = lonlat2meters(lon, lat)
#coord2cell(region, x, y+300)
#cell2coord(region, 10)

function inregion(region::SpatialRegion, lon::Float64, lat::Float64)
    lon >= region.minlon && lon < region.maxlon &&
    lat >= region.minlat && lat < region.maxlat
end

function inregion(region::SpatialRegion, trip::Matrix{Float64})
    for i = 1:size(trip, 2)
        inregion(region, trip[:, i]...) || return false
    end
    true
end

#inregion(region, -8.503, 41.35)
#inregion(region, [-8.503 -8.5; 41.01 41.01])

"""
Build the vocabulary from the raw trajectories stored in the hdf5 file.
For a trip (trajectory), each point lies in a column when reading with Julia
while it lies in a row if reading with Python.
"""
function makeVocab!(region::SpatialRegion, trjfile::String)
    region.cellcount = Accumulator{Int, Int}()
    ## useful for evaluating the size of region bounding box
    num_out_region = 0
    ## scan all trips (trajectories)
    h5open(trjfile, "r") do f
        num = read(attrs(f)["traj_nums"]) - 1
        for i = 0:num-1
            trip = read(f["/trips/$i"])
            for p = 1:size(trip, 2)
                lon, lat = trip[:, p]
                if !inregion(region, lon, lat)
                    num_out_region += 1
                else
                    cell = gps2cell(region, lon, lat)
                    push!(region.cellcount, cell)
                end
            end
            i % 100_000 == 0 && println("Processed $i trips")
            #i >= 8_000 && break
        end
    end
    ## filter out all hot cells
    max_num_hotcells = min(region.maxvocab_size, length(region.cellcount))
    topcellcount = sort(collect(region.cellcount), by=last, rev=true)[1:max_num_hotcells]
    println("Cell count at max_num_hotcells:$(max_num_hotcells) is $(last(topcellcount[end]))")
    #region.hotcell =
    #    filter((k,v)->v >= region.minfreq, Dict(topcellcount))|>keys|>collect
    region.hotcell = filter(p -> last(p) >= region.minfreq, topcellcount) .|> first
    ## build the map between cell and vocab id
    region.hotcell2vocab = Dict([(cell, i-1+region.vocab_start)
        for (i, cell) in enumerate(region.hotcell)])
    #region.vocab2hotcell = map(reverse, region.hotcell2vocab)
    region.vocab2hotcell = Dict(last(p) => first(p) for p in region.hotcell2vocab)
    ## vocabulary size
    region.vocab_size = region.vocab_start + length(region.hotcell)
    ## build the hot cell kdtree to facilitate search
    coord = hcat(map(x->collect(cell2coord(region, x)), region.hotcell)...)
    region.hotcell_kdtree = KDTree(coord)
    region.built = true

    num_out_region
end

function knearestHotcells(region::SpatialRegion, cell::Int, k::Int)
    @assert region.built == true "Build index for region first"
    coord = cell2coord(region, cell) |> collect
    idxs, dists = knn(region.hotcell_kdtree, coord, k)
    region.hotcell[idxs], dists
end

function nearestHotcell(region::SpatialRegion, cell::Int)
    @assert region.built == true "Build index for region first"
    hotcell, _ = knearestHotcells(region, cell, 1)
    hotcell[1]
end

"""
k-nearest vocabs and corresponding distances for each vocab.

This is used in training for KLDiv loss.
"""
function saveKNearestVocabs(region::SpatialRegion)
    V = zeros(Int, region.k, region.vocab_size)
    D = zeros(Float64, region.k, region.vocab_size)
    for vocab in 0:region.vocab_start-1
        V[:, vocab+1] .= vocab
        D[:, vocab+1] .= 0.0
    end
    for vocab in region.vocab_start:region.vocab_size-1
        cell = region.vocab2hotcell[vocab]
        kcells, dists = knearestHotcells(region, cell, region.k)
        kvocabs = map(x->region.hotcell2vocab[x], kcells)
        V[:, vocab+1] .= kvocabs
        D[:, vocab+1] .= dists
    end
    cellsize = Int(region.xstep)
    file = joinpath("../data", region.name * "-vocab-dist-cell$(cellsize).h5")
    h5open(file, "w") do f
        f["V"], f["D"] = V, D
    end
    println("Saved cell distance into $file")
    nothing
end

"""
Return the vocab id for a cell in the region.
If the cell is not hot cell, the function will first search its nearest
hotcell and return the corresponding vocab id.
"""
function cell2vocab(region::SpatialRegion, cell::Int)
    @assert region.built == true "Build index for region first"
    if haskey(region.hotcell2vocab, cell)
        return region.hotcell2vocab[cell]
    else
        hotcell = nearestHotcell(region, cell)
        return region.hotcell2vocab[hotcell]
    end
end

"""
Mapping a gps point to the vocab id in the vocabulary consists of hot cells,
each hot cell has an unique vocab id (hotcell2vocab).

If the point falls out of the region, `UNK` will be returned.
If the point falls into the region, but out of the hot cells its nearest hot cell
will be used.
"""
function gps2vocab(region::SpatialRegion, lon::Float64, lat::Float64)
    inregion(region, lon, lat) || return UNK
    cell2vocab(region, gps2cell(region, lon, lat))
end

function trip2seq(region::SpatialRegion, trip::Matrix{Float64})
    seq = Int[]
    for i in 1:size(trip, 2)
        lon, lat = trip[:, i]
        push!(seq, gps2vocab(region, lon, lat))
    end
    seq |> rle |> first
end

function seq2trip(region::SpatialRegion, seq::Vector{Int})
    trip = zeros(Float64, 2, length(seq))
    for i = 1:length(seq)
        cell = get(region.vocab2hotcell, seq[i], -1)
        cell == -1 && error("i=$i is out of vocabulary")
        lon, lat = cell2gps(region, cell)
        trip[:, i] = [lon, lat]
    end
    trip
end

"""
Create training and validation dataset

createTrainVal(region, "porto.h5", downsampling, 50, 10)
"""
function createTrainVal(region::SpatialRegion,
                        trjfile::String,
                        injectnoise::Function,
                        ntrain::Int,
                        nval::Int;
                        nsplit=5,
                        min_length=20,
                        max_length=100)
    seq2str(seq) = join(map(string, seq), " ") * "\n"
    cityname = region.name
    cellsize = Int(region.xstep)
    h5open(trjfile, "r") do f
        trainsrc, traintrg = open("../data/$cityname-$cellsize-train.src", "w"), open("../data/$cityname-$cellsize-train.trg", "w")
        validsrc, validtrg = open("../data/$cityname-$cellsize-val.src", "w"), open("../data/$cityname-$cellsize-val.trg", "w")
        for i = 1:ntrain+nval
            trip = f["/trips/$i"] |> read
            min_length <= size(trip, 2) <= max_length || continue
            trg = trip2seq(region, trip) |> seq2str
            noisetrips = injectnoise(trip, nsplit)
            srcio, trgio = i <= ntrain ? (trainsrc, traintrg) : (validsrc, validtrg)
            for noisetrip in noisetrips
                ## here: feel weird
                #src = noisetrip |> trip2seq |> seq2str
                src = trip2seq(region, noisetrip) |> seq2str
                write(srcio, src)
                write(trgio, trg)
            end
            i % 100_000 == 0 && println("Scaned $i trips...")
            #i >= 8_000 && break
        end
        close(trainsrc), close(traintrg), close(validsrc), close(validtrg)
    end
    nothing
end

function saveregion(region::SpatialRegion, paramfile::String)
    save(paramfile,
          # JLD cannot handle Accumulator correctly
          #"cellcount", region.cellcount.map,
          "hotcell", region.hotcell,
          "hotcell2vocab", region.hotcell2vocab,
          "vocab2hotcell", region.vocab2hotcell,
          "hotcell_kdtree", region.hotcell_kdtree,
          "vocab_size", region.vocab_size)
    nothing
end

function loadregion!(region::SpatialRegion, paramfile::String)
    jldopen(paramfile, "r") do f
        #region.cellcount = read(f, "cellcount")
        region.hotcell = read(f, "hotcell")
        region.hotcell2vocab = read(f, "hotcell2vocab")
        region.vocab2hotcell = read(f, "vocab2hotcell")
        region.hotcell_kdtree = read(f, "hotcell_kdtree")
        region.vocab_size = read(f, "vocab_size")
        region.built = true
    end
    nothing
end

#end # module
#createTrainVal(region, "porto.h5", downsampling, 50, 10)

max_stars_repo_name: lzzppp/DERT

max_issues_repo_name: lzzppp/DERT

max_forks_repo_name: lzzppp/DERT