text:
using HDF5

d = h5read(ARGS[1], "/")
x = d["eigen_vector_0"]
for i in 1:(length(d["eigen_values"])-1)
    x += d["eigen_vector_$i"]
end

println(length(x))
for i in x
    println(i)
end

max_stars_repo_name: ligmg/ligmg

max_issues_repo_name: ligmg/ligmg

max_forks_repo_name: ligmg/ligmg


text:
const
    C♭♭ = Cff
    C♭ = Cf
    C♯ = Cs
    C♯♯ = Css
    D♭♭ = Dff
    D♭ = Df
    D♯ = Ds
    D♯♯ = Dss
    E♭♭ = Eff
    E♭ = Ef
    E♯ = Es
    E♯♯ = Ess
    F♭♭ = Fff
    F♭ = Ff
    F♯ = Fs
    F♯♯ = Fss
    G♭♭ = Gff
    G♭ = Gf
    G♯ = Gs
    G♯♯ = Gss
    A♭♭ = Aff
    A♭ = Af
    A♯ = As
    A♯♯ = Ass
    B♭♭ = Bff
    B♭ = Bf
    B♯ = Bs
    B♯♯ = Bss

export
    C♭♭, C♭, C♯, C♯♯,
    D♭♭, D♭, D♯, D♯♯,
    E♭♭, E♭, E♯, E♯♯,
    F♭♭, F♭, F♯, F♯♯,
    G♭♭, G♭, G♯, G♯♯,
    A♭♭, A♭, A♯, A♯♯,
    B♭♭, B♭, B♯, B♯♯

c♭♭ = cff
c♭ = cf
c♯ = cs
c♯♯ = css
d♭♭ = dff
d♭ = df
d♯ = ds
d♯♯ = dss
e♭♭ = eff
e♭ = ef
e♯ = es
e♯♯ = ess
f♭♭ = fff
f♭ = ff
f♯ = fs
f♯♯ = fss
g♭♭ = gff
g♭ = gf
g♯ = gs
g♯♯ = gss
a♭♭ = aff
a♭ = af
a♯ = as
a♯♯ = ass
b♭♭ = bff
b♭ = bf
b♯ = bs
b♯♯ = bss

export
    c♭♭, c♭, c♯♯, c♯,
    d♭♭, d♭, d♯♯, d♯,
    e♭♭, e♭, e♯♯, e♯,
    f♭♭, f♭, f♯♯, f♯,
    g♭♭, g♭, g♯♯, g♯,
    a♭♭, a♭, a♯♯, a♯,
    b♭♭, b♭, b♯♯, b♯

ext: jl

max_stars_repo_name: basp/Tarana

max_issues_repo_name: basp/Tarana

max_forks_repo_name: basp/Tarana


text:
using Base: add12, mul12, significand_bits
using Base.Math: ldexp

const SysFloat = Union{Float32, Float64}

# N_min^s : The smallest positive subnormal number (=nextfloat(zero(T)))
# N_max^s : The largest positive subnormal number (=prevfloat(floatmin(T)))
# N_min^n : The smallest positive normal number (=floatmin(T))
# N_max^n : The largest positive normal number (=floatmax(T))

# constants
for T in (Float32, Float64)
    # log_2(N_min^s)
    # N_min^s = 2 * 2^{-precision(T)} * N_min^n
    @eval exponent_smallest_subnormal(::Type{$T}) = $(Int(log2(nextfloat(zero(T)))))
    @eval exponent_product_errorfree_threshold(::Type{$T}) = $(exponent_smallest_subnormal(T) + 2 * significand_bits(T))
    @eval product_errorfree_threshold(::Type{$T}) = $(ldexp(one(T), exponent_product_errorfree_threshold(T)))
    @eval exponent_product_underflow_mult(::Type{$T}) = $(ceil(Int, -exponent_smallest_subnormal(T)//2))
    @eval product_underflow_mult(::Type{$T}) = $(ldexp(one(T), exponent_product_underflow_mult(T)))
    @eval exponent_quotient_errorfree_threshold(::Type{$T}) = $(-exponent_smallest_subnormal(T) - 3 * significand_bits(T))
    @eval quotient_errorfree_threshold(::Type{$T}) = $(ldexp(one(T), exponent_quotient_errorfree_threshold(T)))
    @eval exponent_quotient_underflow_mult(::Type{$T}) = $(2 * significand_bits(T) + 1)
    @eval quotient_underflow_mult(::Type{$T}) = $(ldexp(one(T), exponent_quotient_underflow_mult(T)))
    @eval inverse_smallest_normal(::Type{$T}) = $(ldexp(one(T), precision(T)))
end

"""
    add_up(a, b)

Computes `a + b` with the rounding mode
[`Base.Rounding.RoundUp`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundUp).

```jldoctest
julia> add_up(0.1, 0.2)
0.30000000000000004

julia> add_up(10.0^308, 10.0^308)
Inf

julia> add_up(-10.0^308, -10.0^308)
-1.7976931348623157e308

julia> add_up(-0.1, 0.1)
0.0

julia> add_up(0.0, 0.0)
0.0

julia> add_up(0.0, -0.0)
0.0

julia> add_up(-0.0, -0.0)
-0.0
```
"""
function add_up(a::T, b::T) where {T<:SysFloat}
    x, y = add12(a, b) # twosum
    if isinf(x)
        ifelse(x == typemin(x) && isfinite(a) && isfinite(b), -floatmax(x), x)
    else
        y > zero(y) ? nextfloat(x) : x
    end
end

"""
    add_down(a, b)

Computes `a + b` with the rounding mode
[`Base.Rounding.RoundDown`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundDown).

```jldoctest
julia> add_down(0.1, 0.2)
0.3

julia> add_down(10.0^308, 10.0^308)
1.7976931348623157e308

julia> add_down(-10.0^308, -10.0^308)
-Inf

julia> add_down(-0.1, 0.1)
-0.0

julia> add_down(0.0, 0.0)
0.0

julia> add_down(0.0, -0.0)
-0.0

julia> add_down(-0.0, -0.0)
-0.0
```
"""
function add_down(a::T, b::T) where {T<:SysFloat}
    x, y = add12(a, b) # twosum
    if isinf(x)
        ifelse(x == typemax(x) && isfinite(a) && isfinite(b), floatmax(x), x)
    elseif y < zero(y)
        prevfloat(x)
    else
        ifelse(iszero(x) && (signbit(a) || signbit(b)), -zero(x), x)
    end
end

"""
    sub_up(a, b)

Computes `a - b` with the rounding mode
[`Base.Rounding.RoundUp`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundUp).

```jldoctest
julia> sub_up(-0.1, 0.2)
-0.3

julia> sub_up(-10.0^308, 10.0^308)
-1.7976931348623157e308

julia> sub_up(10.0^308, -10.0^308)
Inf

julia> sub_up(0.1, 0.1)
0.0

julia> sub_up(0.0, 0.0)
0.0

julia> sub_up(0.0, -0.0)
0.0

julia> sub_up(-0.0, 0.0)
-0.0

julia> sub_up(-0.0, -0.0)
0.0
```
"""
sub_up(a::T, b::T) where {T<:SysFloat} = add_up(a, -b)

"""
    sub_down(a, b)

Computes `a - b` with the rounding mode
[`Base.Rounding.RoundDown`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundDown).

```jldoctest
julia> sub_down(-0.1, 0.2)
-0.30000000000000004

julia> sub_down(-10.0^308, 10.0^308)
-Inf

julia> sub_down(10.0^308, -10.0^308)
1.7976931348623157e308

julia> sub_down(0.1, 0.1)
-0.0

julia> sub_down(0.0, 0.0)
-0.0

julia> sub_down(0.0, -0.0)
0.0

julia> sub_down(-0.0, 0.0)
-0.0

julia> sub_down(-0.0, -0.0)
-0.0
```
"""
sub_down(a::T, b::T) where {T<:SysFloat} = add_down(a, -b)

"""
    mul_up(a, b)

Computes `a * b` with the rounding mode
[`Base.Rounding.RoundUp`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundUp).

```jldoctest
julia> mul_up(0.1, 0.2)
0.020000000000000004

julia> mul_up(10.0^308, 10.0^308)
Inf

julia> mul_up(10.0^308, -10.0^308)
-1.7976931348623157e308

julia> mul_up(5.0e-324, 5.0e-324)
5.0e-324

julia> mul_up(-0.1, 0.1)
-0.01

julia> mul_up(0.0, 0.0)
0.0

julia> mul_up(0.0, -0.0)
-0.0

julia> mul_up(-0.0, -0.0)
0.0
```
"""
function mul_up(a::T, b::T) where {T<:SysFloat}
    x, y = mul12(a, b)
    if isinf(x)
        ifelse(x == typemin(x) && isfinite(a) && isfinite(b), -floatmax(x), x)
    elseif abs(x) > product_errorfree_threshold(T) # not zero(x): (a, b) = (-2.1634867667116802e-200, 1.6930929484402486e-119) fails
        y > zero(y) ? nextfloat(x) : x
    else
        mult = product_underflow_mult(T)
        s, s2 = mul12(a * mult, b * mult)
        t = (x * mult) * mult
        t < s || (t == s && s2 > zero(s2)) ? nextfloat(x) : x
    end
end

"""
    mul_down(a, b)

Computes `a * b` with the rounding mode
[`Base.Rounding.RoundDown`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundDown).

```jldoctest
julia> mul_down(0.1, 0.2)
0.02

julia> mul_down(10.0^308, 10.0^308)
1.7976931348623157e308

julia> mul_down(10.0^308, -10.0^308)
-Inf

julia> mul_down(5.0e-324, 5.0e-324)
0.0

julia> mul_down(-0.1, 0.1)
-0.010000000000000002

julia> mul_down(0.0, 0.0)
0.0

julia> mul_down(0.0, -0.0)
-0.0

julia> mul_down(-0.0, -0.0)
0.0
```
"""
function mul_down(a::T, b::T) where {T<:SysFloat}
    x, y = mul12(a, b)
    if isinf(x)
        ifelse(x == typemax(x) && isfinite(a) && isfinite(b), floatmax(x), x)
    elseif abs(x) > product_errorfree_threshold(T) # not zero(x): (a, b) = (6.640350825165134e-116, -1.1053488936824272e-202) fails
        y < zero(y) ? prevfloat(x) : x
    else
        mult = product_underflow_mult(T)
        s, s2 = mul12(a * mult, b * mult)
        t = (x * mult) * mult
        t > s || (t == s && s2 < zero(s2)) ? prevfloat(x) : x
    end
end

"""
    div_up(a, b)

Computes `a / b` with the rounding mode
[`Base.Rounding.RoundUp`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundUp).

```jldoctest
julia> div_up(0.1, 0.3)
0.33333333333333337

julia> div_up(2.0^-100, 2.0^1000)
5.0e-324

julia> div_up(-0.0, 1.0)
-0.0
```
"""
function div_up(a::T, b::T) where {T<:SysFloat}
    if iszero(a) || iszero(b) || isinf(a) || isinf(b) || isnan(a) || isnan(b)
        a / b
    else
        # if b < 0, flip sign of a and b
        a = flipsign(a, b)
        b = abs(b)
        if abs(a) < product_errorfree_threshold(T) && abs(b) < quotient_errorfree_threshold(T)
            mult = quotient_underflow_mult(T)
            a *= mult
            b *= mult
        end
        d = a / b
        x, y = mul12(d, b)
        x < a || (x == a && y < zero(y)) ? nextfloat(d) : d
    end
end

"""
    div_down(a, b)

Computes `a / b` with the rounding mode
[`Base.Rounding.RoundDown`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundDown).

```jldoctest
julia> div_down(0.1, 0.3)
0.3333333333333333

julia> div_down(2.0^-100, 2.0^1000)
0.0

julia> div_down(-0.0, 1.0)
-0.0
```
"""
function div_down(a::T, b::T) where {T<:SysFloat}
    if iszero(a) || iszero(b) || isinf(a) || isinf(b) || isnan(a) || isnan(b)
        a / b
    else
        # if b < 0, flip sign of a and b
        a = flipsign(a, b)
        b = abs(b)
        if abs(a) < product_errorfree_threshold(T) && abs(b) < quotient_errorfree_threshold(T)
                mult = quotient_underflow_mult(T)
                a *= mult
                b *= mult
        end
        d = a / b
        x, y = mul12(d, b)
        x > a || (x == a && y > zero(y)) ? prevfloat(d) : d
    end
end

"""
    sqrt_up(a)

Computes `sqrt(a)` with the rounding mode
[`Base.Rounding.RoundUp`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundUp).

```jldoctest
julia> sqrt_up(2.0)
1.4142135623730951

julia> sqrt_up(0.0)
0.0

julia> sqrt_up(-0.0)
-0.0
```
"""
function sqrt_up(a::SysFloat)
    d = sqrt(a)
    if isinf(d)
        typemax(d)
    elseif a < product_errorfree_threshold(typeof(a))
        invn = inverse_smallest_normal(typeof(a))
        a2 = a * invn^2
        d2 = d * invn
        x, y = mul12(d2, d2)
        x < a2 || (x == a2 && y < zero(y)) ? nextfloat(d) : d
    else
        x, y = mul12(d, d)
        x < a || (x == a  && y < zero(y)) ? nextfloat(d) : d
    end
end

"""
    sqrt_down(a)

Computes `sqrt(a)` with the rounding mode
[`Base.Rounding.RoundDown`](https://docs.julialang.org/en/v1/base/math/#Base.Rounding.RoundDown).

```jldoctest
julia> sqrt_down(2.0)
1.414213562373095

julia> sqrt_down(0.0)
0.0

julia> sqrt_down(-0.0)
-0.0
```
"""
function sqrt_down(a::SysFloat)
    d = sqrt(a)
    if isinf(d)
        typemax(d)
    elseif a < product_errorfree_threshold(typeof(a))
        invn = inverse_smallest_normal(typeof(a))
        a2 = a * invn^2
        d2 = d * invn
        x, y = mul12(d2, d2)
        x > a2 || (x == a2 && y > zero(y)) ? prevfloat(d) : d
    else
        x, y = mul12(d, d)
        x > a || (x == a  && y > zero(y)) ? prevfloat(d) : d
    end
end


text:
"""
    fetch_series_as_at(code::AbstractString, as_at_date::AbstractString)

Fetch series with series_id `code`, with whatever the latest vintage of
the series was as at the date specified by `as_at_date`.

as_at_date must be a date in YYYY-MM-DD form.

Returns a dataframe with four colums: `dates`, `as_at_date`, `series_id`, `values`.
"""
function fetch_series_as_at(code::AbstractString, as_at_date::AbstractString = "latest")
    fetch_series_as_at([code => as_at_date])
end

"""
    fetch_series_as_at(codes::AbstractArray{<:AbstractString, 1}, as_at_date::AbstractString)

Fetch a list of series with codes given by the array `codes`, with whatever the latest vintage of
the series was as at the date specified by `as_at_date`.

as_at_date must be a date in YYYY-MM-DD form.

Returns a dataframe with four colums: `dates`, `as_at_date`, `series_id`, `values`.
"""
function fetch_series_as_at(codes::AbstractArray{<:AbstractString, 1}, as_at_date::AbstractString)
    series = map(c -> c => as_at_date, codes)

    fetch_series_as_at(series)
end

"""
    fetch_series_as_at(series::AbstractArray{T, 1}) where T <: Pair{<:AbstractString, <:AbstractString}

Fetch a list of series with specified series codes, with whatever the latest vintage of
the series was as at the date specified by `as_at_date`.

Each element of the list is a pair `series_code` => `as_at_date`.

as_at_date must be a date in YYYY-MM-DD form.

Returns a dataframe with four colums: `dates`, `as_at_date`, `series_id`, `values`.

# Example
```
fetch_series_as_at(["CPI" => "2021-01-01", "CPI_SYD" => "2019-01-01"])
````
"""
function fetch_series_as_at(series::AbstractArray{T, 1}) where T <: Pair{<:AbstractString, <:AbstractString}
    map(p -> begin
        check_valid_code(p[1])
        check_valid_as_at_date(p[2])
    end, series)

    series_req = join(
        map(p -> """{ "id" : "$(p[1])", "as_at_date" : "$(p[2])" }""", series),
        ", "
    )

    res = HTTP.request(
        "POST",
        "https://aeconomica.io/api/v1/fetchseries",
        [("Content-Type", "application/json")],
        """{
        "series": [
            $series_req
        ],
        "apikey" : "$(apikey())"}""",
        status_exception = false)
    if res.status == 200
        response = JSON.parse(String(res.body), null = missing)
        df = reduce(vcat, map(x -> DataFrames.DataFrame(x), response))
        df.dates = Dates.Date.(df.dates)
        df.vintage = Dates.Date.(df.vintage)
        df.values = if any(ismissing.(df.values))
            Vector{Union{Float64, Missing}}(df.values)
        else
            Vector{Float64}(df.values)
        end

        return df
    else
        error = JSON.parse(String(res.body))["error"]
        if error[1:21] == "500 Internal Error - "
            error = error[22:end]
            throw(ErrorException(error))
        else
            throw(ErrorException(error))
        end
    end
end


text:
@testset "PIP" begin
  G = small_group(8, 3) # D_4
  QG = QQ[G]
  ZG = Order(QG, basis(QG))
  I = 1 * ZG
  Hecke._assert_has_refined_wedderburn_decomposition(QG)
  fl, a = Hecke._isprincipal(I, ZG, :right)
  @test fl
  @test a * ZG == I
  fl, a = Hecke.__isprincipal(ZG, I, :right)
  @test fl
  @test a * ZG == I

  G = small_group(8, 5) # Q_8
  QG = QQ[G]
  ZG = Order(QG, basis(QG))
  I = 1 * ZG
  Hecke._assert_has_refined_wedderburn_decomposition(QG)
  fl, a = Hecke.__isprincipal(ZG, I, :right)
  @test fl
  @test a * ZG == I

  G = small_group(16, 9) # Q_16
  QG = QQ[G]
  ZG = Order(QG, basis(QG))
  I = 1 * ZG
  Hecke._assert_has_refined_wedderburn_decomposition(QG)
  fl, a = Hecke.__isprincipal(ZG, I, :right)
  @test fl
  @test a * ZG == I

  N = Hecke.swan_module(ZG, 3)
  fl, a = Hecke.__isprincipal(ZG, N, :right)
  @test !fl
end


text:
using Interpolations
using Test

@testset "IO" begin
    SPACE = " "

    @testset "BSpline" begin
        A = rand(8,20)

        itp = interpolate(A, BSpline(Constant()))
        @test summary(itp) == "8×20 interpolate(::Matrix{Float64}, BSpline(Constant())) with element type Float64" ||
              summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Constant())) with element type Float64"

        itp = interpolate(A, BSpline(Constant()))
        @test summary(itp) == "8×20 interpolate(::Matrix{Float64}, BSpline(Constant())) with element type Float64" ||
              summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Constant())) with element type Float64"

        itp = interpolate(A, BSpline(Linear()))
        @test summary(itp) == "8×20 interpolate(::Matrix{Float64}, BSpline(Linear())) with element type Float64" ||
              summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Linear())) with element type Float64"

        itp = interpolate(A, BSpline(Quadratic(Reflect(OnCell()))))
        @test summary(itp) == "8×20 interpolate(OffsetArray(::Matrix{Float64}, 0:9, 0:21), BSpline(Quadratic(Reflect(OnCell())))) with element type Float64" ||
              summary(itp) == "8×20 interpolate(OffsetArray(::Array{Float64,2}, 0:9, 0:21), BSpline(Quadratic(Reflect(OnCell())))) with element type Float64"

        itp = interpolate(A, (BSpline(Linear()), NoInterp()))
        @test summary(itp) == "8×20 interpolate(::Matrix{Float64}, (BSpline(Linear()), NoInterp())) with element type Float64" ||
              summary(itp) == "8×20 interpolate(::Array{Float64,2}, (BSpline(Linear()), NoInterp())) with element type Float64"

        itp = interpolate!(copy(A), BSpline(Quadratic(InPlace(OnCell()))))
        @test summary(itp) == "8×20 interpolate(::Matrix{Float64}, BSpline(Quadratic(InPlace(OnCell())))) with element type Float64" ||
              summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Quadratic(InPlace(OnCell())))) with element type Float64"
    end

    @testset "Gridded" begin
        A = rand(20)
        A_x = collect(1.0:2.0:40.0)
        knots = (A_x,)
        itp = interpolate(knots, A, Gridded(Linear()))
        @test summary(itp) == "20-element interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())) with element type Float64" ||
              summary(itp) == "20-element interpolate((::Array{Float64,1},), ::Array{Float64,1}, Gridded(Linear())) with element type Float64"

        A = rand(8,20)
        knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
        itp = interpolate(knots, A, Gridded(Linear()))
        @test summary(itp) == "8×20 interpolate((::Vector{Int64},::Vector{Float64}), ::Matrix{Float64}, Gridded(Linear())) with element type Float64" ||
              summary(itp) == "8×20 interpolate((::Array{Int64,1},::Array{Float64,1}), ::Array{Float64,2}, Gridded(Linear())) with element type Float64"

        itp = interpolate(knots, A, (Gridded(Linear()),Gridded(Constant())))
        @test summary(itp) == "8×20 interpolate((::Vector{Int64},::Vector{Float64}), ::Matrix{Float64}, (Gridded(Linear()), Gridded(Constant()))) with element type Float64" ||
              summary(itp) == "8×20 interpolate((::Array{Int64,1},::Array{Float64,1}), ::Array{Float64,2}, (Gridded(Linear()), Gridded(Constant()))) with element type Float64"

        # issue #260
        A = (1:4)/4
        itp = interpolate((range(0.0, stop=0.3, length=4),), A, Gridded(Linear()))
        io = IOBuffer()
        show(io, MIME("text/plain"), itp)
        str1 = String(take!(io))
        str2 = "4-element interpolate((0.0:0.1:0.3,), ::Array{Float64,1}, Gridded(Linear())) with element type Float64:\n 0.25\n 0.5 \n 0.75\n 1.0"
        str3 = "4-elementinterpolate((0.0:0.1:0.3,),::Vector{Float64},Gridded(Linear()))withelementtypeFloat64:\n 0.25 \n 0.5 \n 0.75\n 1.0"
        @test filter(!isspace, str1) == filter(!isspace, str3) ||
              filter(!isspace, str1) == filter(!isspace, str2)
        io2 = IOBuffer()
        show(io2, itp)
        str1 = String(take!(io2))
        str2 = "4-element interpolate((0.0:0.1:0.3,), ::Array{Float64,1}, Gridded(Linear())) with element type Float64:\n 0.25\n 0.5 \n 0.75\n 1.0"
        str3 = "4-elementinterpolate((0.0:0.1:0.3,),::Vector{Float64},Gridded(Linear()))withelementtypeFloat64:\n 0.25\n 0.5\n 0.75\n 1.0"
        @test filter(!isspace, str1) == filter(!isspace, str3) ||
              filter(!isspace, str1) == filter(!isspace, str2)
    end

    @testset "scaled" begin
        itp = interpolate(1:1.0:10, BSpline(Linear()))
        sitp = scale(itp, -3:.5:1.5)
        @test summary(sitp) == "10-element scale(interpolate(::Vector{Float64}, BSpline(Linear())), (-3.0:0.5:1.5,)) with element type Float64" ||
              summary(sitp) == "10-element scale(interpolate(::Array{Float64,1}, BSpline(Linear())), (-3.0:0.5:1.5,)) with element type Float64"
        io = IOBuffer()
        show(io, MIME("text/plain"), sitp)
        str = String(take!(io))
        @test str == "10-element scale(interpolate(::Vector{Float64}, BSpline(Linear())), (-3.0:0.5:1.5,)) with element type Float64:\n  1.0\n  2.0\n  3.0\n  4.0\n  5.0\n  6.0\n  7.0\n  8.0\n  9.0\n 10.0" ||
              str == "10-element scale(interpolate(::Array{Float64,1}, BSpline(Linear())), (-3.0:0.5:1.5,)) with element type Float64:\n  1.0\n  2.0\n  3.0\n  4.0\n  5.0\n  6.0\n  7.0\n  8.0\n  9.0\n 10.0"
        io2 = IOBuffer()
        show(io2, sitp)
        str2 = String(take!(io2))
        @test str2 == "10-element scale(interpolate(::Vector{Float64}, BSpline(Linear())), (-3.0:0.5:1.5,)) with element type Float64:\n  1.0\n  2.0\n  3.0\n  4.0\n  5.0\n  6.0\n  7.0\n  8.0\n  9.0\n 10.0" ||
              str2 == "10-element scale(interpolate(::Array{Float64,1}, BSpline(Linear())), (-3.0:0.5:1.5,)) with element type Float64:\n  1.0\n  2.0\n  3.0\n  4.0\n  5.0\n  6.0\n  7.0\n  8.0\n  9.0\n 10.0"

        gauss(phi, mu, sigma) = exp(-(phi-mu)^2 / (2sigma)^2)
        testfunction(x,y) = gauss(x, 0.5, 4) * gauss(y, -.5, 2)
        xs = -5:.5:5
        ys = -4:.2:4
        zs = Float64[testfunction(x,y) for x in xs, y in ys]
        itp2 = interpolate(zs, BSpline(Quadratic(Flat(OnGrid()))))
        sitp2 = scale(itp2, xs, ys)
        teststring = "21×41 scale(interpolate(OffsetArray(::Array{Float64,2}, 0:22, 0:42), BSpline(Quadratic(Flat(OnGrid())))), (-5.0:0.5:5.0,$SPACE-4.0:0.2:4.0)) with element type Float64"
        teststring2 = "21×41 scale(interpolate(OffsetArray(::Matrix{Float64}, 0:22, 0:42), BSpline(Quadratic(Flat(OnGrid())))), (-5.0:0.5:5.0,$SPACE-4.0:0.2:4.0)) with element type Float64"
        @test summary(sitp2) == teststring2 ||
              summary(sitp2) == teststring
        io3 = IOBuffer()
        show(io3, sitp2)
        str3 = String(take!(io3))
        @test occursin(teststring2, str3) || occursin(teststring, str3)
    end

    @testset "Monotonic" begin
        A_x = collect(1.0:2.0:40.0)
        A = map(x->x^2, A_x)
        itp = interpolate(A_x, A, FritschButlandMonotonicInterpolation())
        teststring = "20-element interpolate(::Array{Float64,1}, ::Array{Float64,1}, FritschButlandMonotonicInterpolation()) with element type Float64"
        teststring2 = "20-element interpolate(::Vector{Float64}, ::Vector{Float64}, FritschButlandMonotonicInterpolation()) with element type Float64"
        @test summary(itp) == teststring2 || summary(itp) == teststring

        io = IOBuffer()
        show(io, itp)
        str = String(take!(io))
        @test occursin(teststring2, str) || occursin(teststring, str)

        io2 = IOBuffer()
        show(io2, MIME("text/plain"), itp)
        str = String(take!(io2))
        @test occursin(teststring2, str) || occursin(teststring, str)
    end

    @testset "Extrapolation" begin
        A = rand(8,20)

        itpg = interpolate(A, BSpline(Linear()))
        etpg = extrapolate(itpg, Flat())
        teststring = "8×20 extrapolate(interpolate(::Array{Float64,2}, BSpline(Linear())), Flat()) with element type Float64"
        teststring16 = "8×20 extrapolate(interpolate(::Matrix{Float64}, BSpline(Linear())), Flat()) with element type Float64"
        @test summary(etpg) == teststring16 || summary(etpg) == teststring
        io = IOBuffer()
        show(io, etpg)
        str = String(take!(io))
        @test occursin(teststring16, str) || occursin(teststring, str)

        etpf = extrapolate(itpg, NaN)
        teststring2 = "8×20 extrapolate(interpolate(::Array{Float64,2}, BSpline(Linear())), NaN) with element type Float64"
        teststring16 = "8×20 extrapolate(interpolate(::Matrix{Float64}, BSpline(Linear())), NaN) with element type Float64"
        @test summary(etpf) == teststring16 || summary(etpf) == teststring2
        io2 = IOBuffer()
        show(io2, etpf)
        str2 = String(take!(io2))
        @test occursin(teststring16, str2) || occursin(teststring2, str2)
    end

    @testset "Combinations" begin
        xs = 1.0:0.5:2.0
        ys = xs.^3
        interp_cubic = CubicSplineInterpolation(xs, ys)
        io = IOBuffer()
        show(io, MIME("text/plain"), interp_cubic)
        str = String(take!(io))
        @test occursin("3.375", str)
    end

    @testset "Issue #260" begin
        itp = scale(interpolate(zeros(3,3,3),BSpline(Cubic(Free(OnGrid())))),2:4,2:4,2:4)
        io = IOBuffer()
        show(io, MIME("text/plain"), itp)
        str = String(take!(io))
        @test str == "3×3×3 scale(interpolate(OffsetArray(::Array{Float64, 3}, 0:4, 0:4, 0:4), BSpline(Cubic(Free(OnGrid())))), (2:4, 2:4, 2:4)) with element type Float64:\n[:, :, 1] =\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n\n[:, :, 2] =\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n\n[:, :, 3] =\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n 0.0  0.0  0.0" ||
              str == "3×3×3 scale(interpolate(OffsetArray(::Array{Float64,3}, 0:4, 0:4, 0:4), BSpline(Cubic(Free(OnGrid())))), (2:4, 2:4, 2:4)) with element type Float64:\n[:, :, 1] =\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n\n[:, :, 2] =\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n\n[:, :, 3] =\n 0.0  0.0  0.0\n 0.0  0.0  0.0\n 0.0  0.0  0.0"
    end
end # testset


text:
struct MaskedArray{T, N, S<:AbstractArray{T, N}} <: AbstractArray{T, N}
    data::S
    mask::BitArray{N}
end

MaskedArray(data::AbstractArray, mask::AbstractArray{Bool}) =
    MaskedArray(data, BitArray(mask))

Base.size(A::MaskedArray) = size(A.data)

Base.IndexStyle(::Type{<:MaskedArray}) = Base.IndexLinear()

Base.getindex(A::MaskedArray, i::Int) = A.data[i] * A.mask[i]

function Base.setindex!(A::MaskedArray, v, i::Int)
    if A.mask[i]
        A.data[i] = v
    end

    return A
end

Base.similar(A::MaskedArray, ::Type{S}, dims::Dims) where S =
    MaskedArray(similar(A.data, S, dims), trues(dims))

struct MaskedSliceArray{T, N, S<:AbstractArray{T, N}, R<:Tuple} <: AbstractArray{T, N}
    data::S
    slices::R

    function MaskedSliceArray(data::S, slices::R) where {T, N, S<:AbstractArray{T, N}, R<:Tuple}
        length(slices) == N ||
            throw(ArgumentError("MaskedSliceArray(data, slices) expects ndims(data) == length(slices) (got $N != $(length(slices)))."))
        slices = to_indices(data, slices)
        for (d, (s, is)) in enumerate(zip(size(data), slices))
            all(i -> i ∈ 1:s, is) ||
                throw(ArgumentError("Attempted to mask out of bounds indices $is ∉ 1:$s at dimension $d."))
        end

        new{T, N, S, typeof(slices)}(data, slices)
    end
end
MaskedSliceArray(data::AbstractArray{<:Any, N}, slices::Vararg{<:Any, N}) where N =
    MaskedSliceArray(data, slices)

Base.size(A::MaskedSliceArray) = size(A.data)

Base.getindex(A::MaskedSliceArray{T, N}, I::Vararg{Int, N}) where {T, N} =
    any(i ∉ s for (i, s) in zip(I, A.slices)) ? zero(T) : A.data[I...]

function Base.setindex!(A::MaskedSliceArray{<:Any, N}, v, I::Vararg{Int, N}) where N
    if all(i ∈ s for (i, s) in zip(I, A.slices))
        A.data[I...] = v
    end

    return A
end


text:
using Colors, ColorVectorSpace, StaticArrays
using GeometryBasics, Interpolations
using ImageView
using Makie: orthographicprojection

@inline function edge_function(a, b, c)
    (c[1] - a[1]) * (b[2] - a[2]) - (c[2] - a[2]) * (b[1] - a[1])
end

@inline function src_alpha(c::T) where T <: Colorant
    a = alpha(c)
    a == 0.0 && return zero(T)
    c ./ a
end

one_minus_alpha(c::T) where {T <: Colorant} = one(T) .- src_alpha(c)
blend(source, dest, src_func, dest_func) = clamp01(src_func(source) .+ dest_func(dest))
ColorTypes.alpha(x::StaticVector) = x[4]

function standard_transparency(source, dest::T) where T
    (alpha(source) .* source) .+ ((one(eltype(T)) - alpha(source)) .* dest)
end

mutable struct FixedGeomView{GeomOut, VT}
    buffer::Vector{GeomOut}
    view::VT
    idx::Int
end

struct TriangleStripView{T} <: AbstractVector{NTuple{3, T}}
    parent::AbstractVector{T}
end

Base.length(x::TriangleStripView) = 2
Base.size(x::TriangleStripView) = (2,)

function Base.getindex(x::TriangleStripView, i)
    if i === 1
        return ntuple(i-> x.parent[i], 3)
    elseif i === 2
        return map(i-> x.parent[i], (3, 2, 4))
    else
        error("Out of bounds")
    end
end

function FixedGeomView(T, max_primitives, primitive_in, primitive_out)
    buffer = Vector{Tuple{Point4f, T}}(undef, max_primitives)
    # TODO implement primitive_in and out correctly
    # this is for triangle_strip and 4 max_primitives
    if max_primitives != 4 || primitive_out != :triangle_strip
        error("Not implemented for max_primitives == $max_primitives and primitive_out == $primitive_out.")
    end
    geometry_view = if primitive_out == :triangle_strip
        TriangleStripView(buffer)
    else
        error("$primitive_out not supported. Only :triangle_strip supported right now")
    end
    return FixedGeomView(buffer, geometry_view, 1)
end

function reset!(A::FixedGeomView)
    A.idx = 1
end

function Base.push!(A::FixedGeomView, element)
    if A.idx > length(A.buffer)
        error("Emit called more often than max_primitives. max_primitives: $(length(A.buffer))")
    end
    A.buffer[A.idx] = element
    A.idx += 1
    return
end

struct JLRasterizer{Vertex, Args, FragN, VS, FS, GS, GV, EF}
    vertexshader::VS
    fragmentshader::FS

    geometryshader::GS
    geometry_view::GV
    emit::EF
end

function JLRasterizer{Vertex, Args, FragN}(
        vertexshader::VS,
        fragmentshader::FS,
        geometryshader::GS,
        geometry_view::GV,
        emit::EF
    ) where {Vertex, Args, FragN, VS, FS, GS, GV, EF}
    JLRasterizer{Vertex, Args, FragN, VS, FS, GS, GV, EF}(
        vertexshader,
        fragmentshader,
        geometryshader,
        geometry_view,
        emit
    )
end

function geometry_return_type(vertex_array, vertexshader, geometryshader, uniforms)
    typ = Any
    emit_t(position, ::T) where {T} = (typ = T)
    face1 = first(vertex_array)
    vertex_stage = map(reverse(face1)) do f
        vertexshader(f, uniforms...)
    end
    geometryshader(emit_t, vertex_stage, uniforms...) # figure out typ
    typ
end

arglength(::Type{T}) where {T <: Tuple} = length(T.parameters)
arglength(::Type{T}) where {T <: AbstractArray} = 1
arglength(::Type{T}) where {T} = nfields(T)

function rasterizer(
        vertexarray::AbstractArray,
        uniforms::Tuple,
        vertexshader::Function,
        fragmentshader::Function;
        geometryshader = nothing,
        max_primitives = 4,
        primitive_in = :points,
        primitive_out = :triangle_strip,
    )

    emit, geometry_view = nothing, nothing
    fragment_in_ndim = if !isnothing(geometryshader)
        T = geometry_return_type(vertexarray, vertexshader, geometryshader, uniforms)
        geometry_view = FixedGeomView(T, max_primitives, primitive_in, primitive_out)
        emit = (position, fragment_args) -> push!(geometry_view, (position, fragment_args))
        arglength(T)
    else
        # when we don't have a geometry shader, vertex shader will feed fragment shader
        T = Base.Core.Inference.return_type(vertexshader, Tuple{eltype(vertexarray), map(typeof, uniforms)...})
        if T <: Tuple
            # TODO error handling
            arglength(T.parameters[2])
        else # if not a tuple, vertex shader doesn't pass any arguments to fragment shader
            0
        end
    end

    raster = JLRasterizer{eltype(vertexarray), typeof(uniforms), fragment_in_ndim}(
        vertexshader,
        fragmentshader,
        geometryshader,
        geometry_view,
        emit
    )
    raster, (vertexarray, uniforms)
end

Base.@pure Next(::Val{N}) where {N} = Val(N - 1)

function interpolate(bary, face::NTuple{N, T}, vn::Val{0}, aggregate) where {N, T}
    if T <: Tuple
        aggregate
    else
        T(aggregate...)
    end
end

function interpolate(bary, face, vn::Val{N}, aggregate = ()) where N
    @inbounds begin
        res = (
            bary[1] * getfield(face[1], N) .+
            bary[2] * getfield(face[2], N) .+
            bary[3] * getfield(face[3], N)
        )
    end
    return interpolate(bary, face, Next(vn), (res, aggregate...))
end

broadcastmin(a, b) = min.(a, b)
broadcastmax(a, b) = max.(a, b)

function clip2pixel_space(position, resolution)
    clipspace = position ./ position[4]
    p = clipspace[Vec(1, 2)]
    (((p .+ 1f0) / 2f0) .* (resolution .- 1f0)) .+ 1f0, clipspace[3]
end

function (r::JLRasterizer{Vert, Args, FragN})(
        canvas, vertex_array::AbstractArray{Vert}, uniforms::Args
    ) where {Vert, Args, FragN}
    framebuffers = canvas.color; depthbuffer = canvas.depth
    resolution = Vec2f(size(framebuffers[1]))
    # hoisting out functions... Seems to help inference a bit. Or not?
    vshader = r.vertexshader
    gshader = r.geometryshader
    fshader = r.fragmentshader
    FragNVal = Val(FragN)
    fragments_drawn = 0
    for face in vertex_array
        vertex_stage = map(reverse(face)) do f
            vshader(f, uniforms...)
        end

        geom_stage = if isnothing(gshader)
            (vertex_stage,)
        else
            reset!(r.geometry_view)
            gshader(r.emit, vertex_stage, uniforms...)
            r.geometry_view.view
        end

        for geom_face in geom_stage
            fdepth = map(geom_face) do vert
                fv = first(vert)
                clip2pixel_space(fv, resolution)
            end
            f = map(first, fdepth)
            depths = map(last, fdepth)
            vertex_out = map(last, geom_face)
            # Bounding rectangle
            mini = max.(reduce(broadcastmin, f), 1f0)
            maxi = min.(reduce(broadcastmax, f), resolution)
            area = edge_function(f[1], f[2], f[3])
            for y = mini[2]:maxi[2]
                for x = mini[1]:maxi[1]
                    p = Vec(x, y)
                    w = Vec(
                        edge_function(f[2], f[3], p),
                        edge_function(f[3], f[1], p),
                        edge_function(f[1], f[2], p)
                    )
                    yi, xi = round(Int, y), round(Int, x)
                    if all(w-> w >= 0f0, w) && checkbounds(Bool, framebuffers[1], yi, xi)

                        bary = w / area
                        depth = bary[1] * depths[1] + bary[2] * depths[2] + bary[3] * depths[3]

                        if depth <= depthbuffer[yi, xi]
                            depthbuffer[yi, xi] = depth
                            fragment_in = interpolate(bary, vertex_out, FragNVal)
                            fragment_out = fshader(fragment_in, uniforms...)
                            for i = eachindex(fragment_out)
                                src_color = framebuffers[i][yi, xi]
                                dest_color = fragment_out[i]
                                fragments_drawn += 1
                                framebuffers[i][yi, xi] = standard_transparency(
                                    src_color,
                                    RGBA{Float32}(dest_color[1], dest_color[2], dest_color[3], dest_color[4])
                                )
                            end
                        end
                    end
                end
            end
        end
    end
    println("fragments drawn: ", fragments_drawn)
    return
end

"""
smoothstep performs smooth Hermite interpolation between 0 and 1 when edge0 < x < edge1. This is useful in cases where a threshold function with a smooth transition is desired. smoothstep is equivalent to:
```
    t = clamp((x - edge0) / (edge1 - edge0), 0.0, 1.0);
    return t * t * (3.0 - 2.0 * t);
```
Results are undefined if edge0 ≥ edge1.
"""
function smoothstep(edge0, edge1, x::T) where T
    t = clamp.((x .- edge0) ./ (edge1 .- edge0), T(0), T(1))
    return t * t * (T(3) - T(2) * t)
end

function aastep(threshold1::T, value) where T
    afwidth = norm(Vec2f(dFdx(value), dFdy(value))) * T(1.05);
    smoothstep(threshold1 - afwidth, threshold1 + afwidth, value)
end

function aastep(threshold1::T, threshold2::T, value::T) where T
    afwidth = norm(Vec2f(dFdx(value), dFdy(value))) * T(1.05);
    return (
        smoothstep(threshold1 - afwidth, threshold1 + afwidth, value) -
        smoothstep(threshold2 - afwidth, threshold2 + afwidth, value)
    )
end

mix(x, y, a::T) where {T} = x .* (T(1) .- a) .+ y .* a
fract(x) = x - floor(x)
fabs(x::AbstractFloat) = abs(x)

"""
Gradient in x direction
This is sadly a bit hard to implement for a pure CPU versions, since it's pretty much backed into the GPU hardware.
How it seems to work is, that it takes the values from neighboring registers, which work in parallel on the pixels
of the triangle, so they actually do hold the neighboring values needed to calculate the gradient.
"""
dFdx(value::T) where {T} = T(0.001) # just default to a small gradient if it's called on the CPU
dFdy(value::T) where {T} = T(0.001) # just default to a small gradient if it's called on the CPU

mutable struct Uniforms{F}
    projection::Mat4f
    strokecolor::Vec4f
    glowcolor::Vec4f
    distance_func::F
end

struct VertexCS{N, T}
    position::Vec{N, T}
    color::Vec4f
    scale::Vec2f
end

struct Vertex2Geom
    uvrect::Vec4f
    color::Vec4f
    rect::Vec4f
end

function vert_particles(vertex, uniforms)
    p = vertex.position
    scale = vertex.scale
    return Vertex2Geom(
        Vec4f(0,0,1,1),
        vertex.color,
        Vec4f(p[1], p[2], scale[1], scale[2])
    )
end

"""
Emits a vertex with
"""
function emit_vertex(emit!, vertex, uv, arg, pos, uniforms)
    datapoint = uniforms.projection * Vec4f(pos[1], pos[2], 0, 1)
    final_position = uniforms.projection * Vec4f(vertex[1], vertex[2], 0, 0)
    emit!(datapoint .+ final_position, (uv, arg.color))
    return
end

function geom_particles(emit!, vertex_out, uniforms, image)
    geom_particles(emit!, vertex_out, uniforms)
    return
end

function geom_particles(emit!, vertex_out, uniforms)
    # get arguments from first face
    # (there is only one in there anywas, since primitive type is point)
    # (position, vertex_out)
    arg = vertex_out[1]
    # emit quad as triangle strip
    # v3. ____ . v4
    #    |\   |
    #    | \  |
    #    |  \ |
    #    |___\|
    # v1*      * v2
    pos_scale = arg.rect
    pos = pos_scale[Vec(1, 2)]
    scale = pos_scale[Vec(3, 4)]
    quad = Vec4f(0f0, 0f0, scale[1], scale[2])
    uv = arg.uvrect
    emit_vertex(emit!, quad[Vec(1, 2)], uv[Vec(1, 4)], arg, pos, uniforms)
    emit_vertex(emit!, quad[Vec(1, 4)], uv[Vec(1, 2)], arg, pos, uniforms)
    emit_vertex(emit!, quad[Vec(3, 2)], uv[Vec(3, 4)], arg, pos, uniforms)
    emit_vertex(emit!, quad[Vec(3, 4)], uv[Vec(3, 2)], arg, pos, uniforms)
    return
end

function sdf2color(dist, bg_color, color)
    inside = aastep(0f0, dist)
    mix(bg_color, color, inside)
end

function frag_particles(geom_out, uniforms, image)
    uv = geom_out[1]; color = geom_out[2]
    dist = -image[uv][1]
    bg_color = Vec4f(0f0, 0f0, 0f0, 0f0)
    (sdf2color(dist, bg_color, color), )
end

function frag_particles(geom_out, uniforms)
    uv = geom_out[1]; color = geom_out[2]
    dist = uniforms.distance_func(uv)
    bg_color = Vec4f(0f0, 0f0, 0f0, 0f0)
    # col = Vec4f(norm(uv .- 0.5), 0, 0, 1)
    (sdf2color(dist, bg_color, color), )
end

resolution = (1024, 1024)

proj = Makie.orthographicprojection(Rect2(0, 0, resolution...), -10_000f0, 10_000f0)

circle(uv::Vec{2, T}) where {T} = -T(norm(uv .- 0.5) - 0.5)

uniforms = Uniforms(
    proj,
    Vec4f(1, 0, 0, 1),
    Vec4f(1, 0, 1, 1),
    circle
)

N = 10
middle = Vec2f(resolution) / 2f0
radius = (min(resolution...) / 2f0) - 50
vertices = [(VertexCS(
    Vec2f((sin(2pi * (i / N)) , cos(2pi * (i / N))) .* radius) .+ middle,
    Vec4f(1, i/N, 0, 1),
    Vec2f(40, 40)
),) for i = 1:N]

raster, rest = rasterizer(
    vertices,
    (uniforms,),
    vert_particles,
    frag_particles;
    geometryshader = geom_particles,
    max_primitives = 4,
    primitive_in = :points,
    primitive_out = :triangle_strip,
)

struct Canvas{N}
    color::NTuple{N, Matrix{RGBA{Float32}}}
    depth::Matrix{Float32}
end

function Canvas(xdim::Integer, ydim::Integer)
    color = fill(RGBA{Float32}(0,0,0,0), xdim, ydim)
    depth = ones(Float32, xdim, ydim)
    return Canvas((color,), depth)
end

c = Canvas(resolution...)
raster(c, vertices, (uniforms,))
save("test.png", c.color[1])


text:
@testset "getGaussPoints" begin
    domain = example_domain(1, 1, 1.0)
    @test getGaussPoints(domain) ≈ [0.21132486540518713 0.21132486540518713; 0.21132486540518713 0.7886751345948129; 0.7886751345948129 0.21132486540518713; 0.7886751345948129 0.7886751345948129]
end

@testset "getStrain" begin
    domain = example_domain(10,10,0.1)
    domain.state = rand(domain.neqs)
    E = getStrain(domain)
    E2 = eval_strain_on_gauss_pts(domain.state, 10,10,0.1)
    @test E≈E2
end

@testset "getStress" begin
    domain = example_domain(10,10,0.1)
    domain.state = rand(domain.neqs)
    E = getStrain(domain)
    S = E*domain.elements[1].mat[1].H
    @test getStress(domain)≈S
end


text:
function basispursuit(A::AbstractMatrix, b::AbstractVector)
    basispursuit(A, b, ones(eltype(A), size(A, 2)))
end
const bp = basispursuit

function basispursuit(A::AbstractMatrix, b::AbstractVector, w::AbstractVector)
    model = JuMP.Model(optimizer_with_attributes(Clp.Optimizer, "LogLevel" => 0))
    n, m = size(A)
    @variable(model, x⁺[1:m] ≥ 0)
    @variable(model, x⁻[1:m] ≥ 0)
    @objective(model, Min, dot(w, x⁺) + dot(w, x⁻))
    @constraint(model, constraint, A * (x⁺ .- x⁻) .== b)
    JuMP.optimize!(model)
    x = @. JuMP.value(x⁺) - JuMP.value.(x⁻)
    sparse(x)
end

function basispursuit_reweighting(A::AbstractMatrix, b::AbstractVector,
                                    reweighting!; maxiter = 8, min_decrease = 1e-8)
    x = basispursuit(A, b)
    w = ones(length(x))
    for i in 2:maxiter
        reweighting!(w, x) # updates w
        xs = basispursuit(A, b, w)
        if norm(xs-x) < min_decrease # stop if we don't make sufficient progress
            return xs
        end
        x .= xs
    end
    x
end

candes_weight(x, ε) = inv(abs(x) + ε)
function candes_weights!(w, x, ε::Real)
    @. w = candes_weight(x, ε)
    if any(isnan, w) || any(isinf, w)
        throw("weights contain NaN or Inf: $w")
    end
end
candes_function(ε::Real) = (w, x) -> candes_weights!(w, x, ε)
function candes_reweighting(A::AbstractMatrix, b::AbstractVector,
                                            ε = 1e-2; maxiter = 8)
    basispursuit_reweighting(A, b, candes_function(ε), maxiter = maxiter)
end
const bp_candes = candes_reweighting

# using LinearAlgebraExtensions: LowRank
# TODO: stop when w has converged, instead of maxiter
function ard_weights!(w, A, x, ε::Real, iter::Int = 8)
    if any(==(0), w)
        error("weights cannot be zero")
    end
    # nzind = x.nzind # this could speed up the weight computation
    # Ai = @view A[:, nzind]
    # wx = @. abs(x.nzval) / w[nzind]
    # K = Woodbury(ε*I(size(A, 1)), Ai, WX, Ai')
    for i in 1:iter
        wx = @. abs(x) / w
        WX = Diagonal(wx)
        K = Woodbury(ε*I(size(A, 1)), A, WX, A')
        K⁻¹ = inverse(factorize(K))
        dK(a) = sqrt(max(dot(a, K⁻¹, a), 0))
        @. w = dK($eachcol(A))
    end
end

ard_function(A::AbstractMatrix, ε::Real) = (w, x) -> ard_weights!(w, A, x, ε)

function ard_reweighting(A::AbstractMatrix, b::AbstractVector, ε::Real = 1e-2;
                                                            maxiter::Int = 8)
    basispursuit_reweighting(A, b, ard_function(A, ε), maxiter = maxiter)
end
const bp_ard = ard_reweighting

############################# Basis Pursuit Denoising ##########################
# noisy variant, where we allow δ-deviations from data
# minimizes l2 norm constrained linear objective (conic program)
# min_x |x|_1 s.t. |Ax-b|_2 < δ
function basis_pursuit_denoising(A::AbstractMatrix, b::AbstractVector, δ::Real,
                                    w::AbstractVector = ones(eltype(A), size(A, 2)))

    model = bpd_ecos_model()
    n, m = size(A)
    @variable(model, x⁺[1:m] ≥ 0)
    @variable(model, x⁻[1:m] ≥ 0)
    @objective(model, Min, dot(w, x⁺) + dot(w, x⁻))
    @constraint(model, [δ; A * (x⁺ .- x⁻) .- b] in SecondOrderCone())
    # return model
    JuMP.optimize!(model)
    try
        x = @. JuMP.value(x⁺) - JuMP.value(x⁻)
        return sparse(x)
    catch e
        println(e)
        x = fill(NaN, size(A, 2))
        return sparse(x)
    end
end
const bpd = basis_pursuit_denoising

function bpd_reweighting(A::AbstractMatrix, b::AbstractVector, δ::Real,
                                    reweighting!; maxiter = 8, min_decrease = 1e-4)
    x = bpd(A, b, δ)
    w = ones(length(x))
    for i in 2:maxiter
        reweighting!(w, x) # updates w
        xs = bpd(A, b, δ, w)
        if norm(xs-x) < min_decrease # stop if we don't make sufficient progress
            return xs
        end
        x .= xs
    end
    x
end

# δ is l2 noise threshold, ε is reweighting constant
# ε below 1e-3 can lead to convergence problems
function bpd_candes(A::AbstractMatrix, b::AbstractVector, δ::Real, ε::Real = δ; maxiter = 8)
    bpd_reweighting(A, b, δ, candes_function(ε), maxiter = maxiter)
end
function bpd_ard(A::AbstractMatrix, b::AbstractVector, δ::Real, ε::Real = δ^2; maxiter = 8)
    bpd_reweighting(A, b, δ, ard_function(A, ε), maxiter = maxiter)
end

# optimizer setup
function bpd_scs_model()
    JuMP.Model(optimizer_with_attributes(SCS.Optimizer,
                "verbose" => true,
                "eps" => 1e-6,
                "alpha" => 1.,
                "cg_rate" => 2.,
                "max_iters" => 5000))
end
function bpd_ecos_model()
    JuMP.Model(optimizer_with_attributes(ECOS.Optimizer, "verbose" => false))
end

################################ (F)ISTA ########################################
struct FISTA end
# TODO: use homotopy to select appropriate λ

# shrinkage and thresholding operator
shrinkage(x::Real, α::Real) = sign(x) * max(abs(x)-α, zero(x))

# TODO fast iterative shrinkage and hard-thresholding algorithm
# for weighted l1-norm minimization
function fista(A::AbstractMatrix, b::AbstractVector, λ::Real,
    x::AbstractVector = spzeros(size(A, 2)); maxiter::Int = 1024, stepsize::Real = 1e-2)
    w = fill(λ, size(x))
    fista(A, b, w, x, maxiter = maxiter, stepsize = stepsize)
end

function l1(x::AbstractVector, w::AbstractVector)
    length(x) == length(w) || throw(DimensionMismatch("length(x) ≠ length(w)"))
    n = zero(eltype(x))
    @simd for i in eachindex(x)
        @inbounds n += w[i] * abs(x[i])
    end
    return n
end

# TODO: stepsize selection
function ista(A, b, λ::Real, x = spzeros(size(A, 2)); maxiter::Int = 1024, stepsize::Real = 1e-2)
    ista(A, b, fill(λ, size(x)), x, maxiter = maxiter, stepsize = stepsize)
end

function ista(A::AbstractMatrix, b::AbstractVector, w::AbstractVector,
                                        x::AbstractVector = spzeros(size(A, 2));
                                        maxiter::Int = 1024, stepsize::Real = 1e-2)
    x = sparse(x)
    r(x) = b-A*x # residual
    f(x) = sum(abs2, r(x)) + l1(x, w)
    g(x) = A'r(x) # negative gradient
    α = stepsize
    fx = f(x)
    for i in 1:maxiter
        ∇ = g(x)
        @. x = shrinkage(x + 2α*∇, w*α)
        dropzeros!(x) # optimize sparse representation
    end
    return x
end

function fista(A::AbstractMatrix, b::AbstractVector, w::AbstractVector,
                                        x::AbstractVector = spzeros(size(A, 2));
                                        maxiter::Int = 1024, stepsize::Real = 1e-2)
    x = sparse(x)
    r(x) = b-A*x # residual
    f(x) = sum(abs2, r(x)) + l1(x, w)
    g(x) = A'r(x) # negative gradient
    α = stepsize
    tk = 1
    for i in 1:maxiter
        ∇ = g(x)
        tkn = (1 + sqrt(1 + 4tk^2)) / 2
        y = xkn + (tk - 1) / tkn * (xkn - xk)
        tk = tkn
        @. x = shrinkage(x + 2*∇, w*α)
        dropzeros!(x) # optimize sparse representation
    end
    return x
end


text:
using NNNCVXF #(https://github.com/PetersBas/NNNCVXF)
using PyPlot
using SetIntersectionProjection #(https://github.com/slimgroup/SetIntersectionProjection.jl)
using InvertibleNetworks #(https://github.com/slimgroup/InvertibleNetworks.jl)
using LinearAlgebra
using Flux
import Flux: gradient
import Flux.Losses.crossentropy

# We use the CamVid data prepared by:
# https://github.com/alexgkendall/SegNet-Tutorial/tree/master/CamVid

use_gpu = true

cd("/home/bpete23/Results/") #working directory
basepath = "/home/bpete23/Data"
#basepath = "/Users/BasP/Downloads/"

nr_output_chan = 3*16
train_data, train_labels, val_data, val_labels = LoadCamvid(basepath,nr_output_chan)

active_channels = [1,2]

n_pix_per_class = 8
output_samples_train = Vector{Array{Float32,2}}(undef,length(train_data))
output_samples_val   = Vector{Array{Float32,2}}(undef,length(val_data))
for i=1:length(output_samples_train)
  output_samples_train[i] = SparseClassSamples(train_labels[i],active_channels,n_pix_per_class)
end
for i=1:length(output_samples_val)
  output_samples_val[i] = SparseClassSamples(val_labels[i],active_channels,n_pix_per_class)
end

mutable struct grid #grid structure for discrete differential operators
    n
    d
end
compgrid = grid(size(train_labels[1])[1:2], (1, 1))

noise_factor = 0.33
TV_list_train = GetTVvalues(grid,train_labels,noise_factor,compgrid)

#set up one projector per example
options=PARSDMM_options()
options.FL=Float32
options=default_PARSDMM_options(options,options.FL)
options.maxit = 1000

#using these (below) more accurate projections onto the anisotropic total-variation constraint set_type
#give slightly better segmentation reguls, but it takes significantly more computational time
options.feas_tol     = 0.0001
options.obj_tol      = 0.0001
options.evol_rel_tol = 0.00001

BLAS.set_num_threads(4)

constraint = Vector{SetIntersectionProjection.set_definitions}()

m_min        = 0.0
m_max        = 1f0#dummy value that will be replaced in the loop below true_TV_list[i]
set_type     = "l1"
TD_OP        = "TV"
app_mode     = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP));

#set up constraints, precompute some things and define projector
(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,compgrid,options.FL)
(TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,compgrid,options)

P_sub_list = Vector{Any}(undef,length(train_data))
P_train    = Vector{Any}(undef,length(train_data))
P_val      = Vector{Any}(undef,length(val_data))
for i=1:length(train_data)
  P_sub_list[i] = Vector{Any}(undef,1)
  P_sub_list[i][1] = x-> project_l1_Duchi!(x,Float32(TV_list_train[i]))
end
for i=1:length(train_data)
  P_train[i] = Vector{Any}(undef,2)
  function prj(x)
    os=size(x)
    (x,dummy1,dummy2,dymmy3) = PARSDMM(deepcopy(vec(x)), AtA, TD_OP, set_Prop, P_sub, compgrid, options)
    return reshape(x,os)
  end
  P_train[i][1] = x -> prj(x)
  P_train[i][2] = x -> x #don't apply any projection to 2nd channel (not necessary for good results in this example)
end

#No projection for validation set
for i=1:length(val_data)
  P_val[i] = Vector{Any}(undef,2)
  P_val[i][1] = x-> x
  P_val[i][2] = x-> x
end

#Identity mapping (in place of the projector in case we don't want to project output for training data)
P_I = Vector{Any}(undef,length(train_data))
for i=1:length(train_data)
  P_I[i] = Vector{Any}(undef,2)
  P_I[i][1] = x-> x
  P_I[i][2] = x-> x
end

architecture = ((0, 48), (0, 48), (0, 48),(0, 48), (0, 48), (0, 48),(0, 48),(0, 48),(0, 48),(0,48),(0, 48),(0,48),(0, 48),(0, 48))

k = 3   # kernel size
s = 1   # stride
p = 1   # padding
n_chan_in = size(train_data[1],3)
#(nx, ny)  = size(train_data[1])[1:2]
α = 0.3^2 #artificial time-step in the nonlinear telegraph equation discretization that is the neural network
#HN = H = NetworkHyperbolic(nx,ny,n_chan_in,1,architecture; α)
HN = H = NetworkHyperbolic(n_chan_in,architecture; α)

# define loss functions for the labels and corresponding gradients
CESM   = (x,y)-> crossentropy(softmax(x,dims=3), y; dims=3, agg=sum,ϵ=eps(x[1]))
wCESM  = (x,y,c)-> crossentropy(softmax(x,dims=3), y; dims=3, agg=x->sum(c .* x),ϵ=eps(x[1]))
gCESM  = (x,y)-> gradient(CESM,x,y)
gwCESM = (x,y,c)-> gradient(wCESM,x,y,c)

########## EXPERIMENT: Point annotations only #################################
TrainOptions = TrOpts()
TrainOptions.eval_every      = 50
TrainOptions.batchsize       = 4
TrainOptions.use_gpu         = use_gpu
TrainOptions.lossf           = wCESM
TrainOptions.lossg           = gwCESM
TrainOptions.active_channels = [1,2]
TrainOptions.flip_dims       = []
TrainOptions.permute_dims    = []

alpha   = [0f0 0f0]
maxiter = [1200 800]
opt     = [Flux.ADAM(1f-3) Flux.ADAM(1f-4)]

  #set all data and labels to GPU
  if TrainOptions.use_gpu==true
    HN         = H |> gpu
    val_data   = val_data|>gpu
    train_data = train_data|>gpu
  end

  logs = Log()

  for i=1:length(maxiter)
    TrainOptions.alpha   = alpha[i]
    TrainOptions.opt     = opt[i]
    TrainOptions.maxiter = maxiter[i]
    logs = Train(HN,logs,TrainOptions,train_data,val_data,train_labels,val_labels,P_I,output_samples_train,output_samples_val)
  end

  if TrainOptions.use_gpu==true
    train_data  = train_data |>cpu
    val_data    = val_data |>cpu
    HN          = HN |>cpu
  end

  if isdir("camvid_unconstrained") == true
    cd("camvid_unconstrained")
  else
    mkdir("camvid_unconstrained")
    cd("camvid_unconstrained")
  end

  PlotCamvidLossUnconstrained(logs,TrainOptions.eval_every)

  #plot training result
  plt_ind = 10
  tag     = "train"
  PlotDataLabelPredictionCamvid(plt_ind,train_data,train_labels,HN,TrainOptions.active_channels,tag)

  #plot validation result
  plt_ind = 3
  tag     = "validation"
  PlotDataLabelPredictionCamvid(plt_ind,val_data,val_labels,HN,TrainOptions.active_channels,tag)

  cd("..")

########## EXPERIMENT: point_annotations_+_TV_constraints #####################
  #set all data and labels to GPU
  if use_gpu==true
    HN         = H |> gpu
    val_data   = val_data|>gpu
    train_data = train_data|>gpu
  end

  TrainOptions = TrOpts()
  TrainOptions.eval_every      = 50
  TrainOptions.batchsize       = 4
  TrainOptions.use_gpu         = true
  TrainOptions.lossf           = wCESM
  TrainOptions.lossg           = gwCESM
  TrainOptions.active_channels = [1,2]
  TrainOptions.flip_dims       = []
  TrainOptions.permute_dims    = []

  alpha   = [0f0 0f0 1f-5 5f-5 1f-4 5f-4 1f-3 5f-3 1f-2]
  maxiter = [1200 800 400 400 400 400 800 800 800]
  opt     = [Flux.ADAM(1f-3) Flux.ADAM(1f-4) Flux.ADAM(1f-4) Flux.ADAM(1f-4) Flux.ADAM(1f-4) Flux.ADAM(1f-4) Flux.ADAM(1f-4) Flux.ADAM(1f-4) Flux.ADAM(1f-4)]

  logs = Log()

  for i=1:2
    TrainOptions.alpha   = alpha[i]
    TrainOptions.opt     = opt[i]
    TrainOptions.maxiter = maxiter[i]
    logs = Train(HN,logs,TrainOptions,train_data,val_data,train_labels,val_labels,P_I,output_samples_train,output_samples_val)
  end
  for i=3:length(maxiter)
    TrainOptions.alpha   = alpha[i]
    TrainOptions.opt     = opt[i]
    TrainOptions.maxiter = maxiter[i]
    logs = Train(HN,logs,TrainOptions,train_data,val_data,train_labels,val_labels,P_train,output_samples_train,output_samples_val)
  end

  if use_gpu==true
    train_data  = train_data |>cpu
    val_data    = val_data |>cpu
    HN          = HN |>cpu
  end

  if isdir("camvid_constrained") == true
    cd("camvid_constrained")
  else
    mkdir("camvid_constrained")
    cd("camvid_constrained")
  end

  PlotCamvidLossConstrained(logs,TrainOptions.eval_every)

  #plot training result
  plt_ind = 10
  tag     = "train"
  PlotDataLabelPredictionCamvid(plt_ind,train_data,train_labels,HN,TrainOptions.active_channels,tag)

  #plot validation result
  plt_ind = 3
  tag     = "validation"
  PlotDataLabelPredictionCamvid(plt_ind,val_data,val_labels,HN,TrainOptions.active_channels,tag)

  cd("..")


text:
export PerfRNN, generate

using StatsBase
using Flux: gate

import ..MusicTransformer: generate

struct PerfRNN
    model
    encoder::PerformanceOneHotEncoding
    num_velocitybins::Int
end

"""
    generate(perfrnn::PerfRNN, performance::Performance)

Generate `PerformanceEvent`s by sampling from a PerformanceRNN model.
"""
function generate(perfrnn::PerfRNN;
                  primer::Performance=Performance(100, velocity_bins=perfrnn.num_velocitybins),
                  numsteps=3000,
                  as_notesequence=false)

    model, encoder = perfrnn.model, perfrnn.encoder
    Flux.reset!(model)

    performance = deepcopy(primer)

    if isempty(performance)
        push!(performance, encoder.defaultevent)
    end

    # Primer is already numsteps or longer
    if performance.numsteps >= numsteps
        return performance
    end

    indices = map(event -> encode_event(event, encoder), performance)
    inputs = map(index -> Flux.onehot(index, encoder.labels), indices)

    outputs = model.(inputs)
    out = wsample(encoder.labels, softmax(outputs[end]))
    push!(performance, decode_event(out, encoder))

    while performance.numsteps < numsteps
        input = Flux.onehot(out, encoder.labels)
        logits = model(input)
        out = wsample(encoder.labels, softmax(logits))
        push!(performance, decode_event(out, encoder))
    end

    ns = getnotesequence(performance)
    as_notesequence == true && return ns
    midifile(ns)
end


text:
using IterativeSolvers, KrylovKit, Arpack, LinearAlgebra

# abstract type for generalised eigenvector
abstract type AbstractGEigenSolver <: AbstractEigenSolver end
# abstract type for Matrix-Free eigensolvers
abstract type AbstractGMFEigenSolver <: AbstractMFEigenSolver end
abstract type AbstractGFloquetSolver <: AbstractFloquetSolver end

convertToGEV(l::AbstractGEigenSolver, B) = l
####################################################################################################
# Solvers for default \ operator (backslash)
####################################################################################################
"""
The struct `Default` is used to  provide the backslash operator to our Package
"""
@with_kw struct DefaultGEig{T, Tb <: AbstractMatrix} <: AbstractGEigenSolver
	which::T = real		# how do we sort the computed eigenvalues
	B::Tb
end

function (l::DefaultGEig)(Jac, nev; kwargs...)
	# I put Array so we can call it on small sparse matrices
	F = eigen(Array(Jac), l.B)
	I = sortperm(F.values, by = l.which, rev = true)
	nev2 = min(nev, length(I))
	J = findall( abs.(F.values[I]) .< 100000)
	# we perform a conversion to Complex numbers here as the type can change from Float to Complex along the branch, this would cause a bug
	return Complex.(F.values[I[J[1:nev2]]]), Complex.(F.vectors[:, I[J[1:nev2]]]), true, 1
end

convertToGEV(l::DefaultEig, B) = DefaultGEig(l.which, Array(B)) # we convert B from sparse to Array
####################################################################################################
# case of sparse matrices or matrix free method
####################################################################################################
"""
$(TYPEDEF)
$(TYPEDFIELDS)

More information is available at [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl). You can pass the following parameters `tol=0.0, maxiter=300, ritzvec=true, v0=zeros((0,))`.

# Constructor

`EigArpack(sigma = nothing, which = :LR; kwargs...)`
"""
struct GEigArpack{T, Tby, Tw, Tb} <: AbstractGEigenSolver
	"Shift for Shift-Invert method with `(J - sigma⋅I)"
	sigma::T

	"Which eigen-element to extract :LR, :LM, ..."
	which::Symbol

	"Sorting function, default to real"
	by::Tby

	"Keyword arguments passed to EigArpack"
	kwargs::Tw

	"Mass matrix"
	B::Tb
end

GEigArpack(sigma = nothing, which = :LR; kwargs...) = EigArpack(sigma, which, real, kwargs)
convertToGEV(l::EigArpack, B) = GEigArpack(l.sigma, l.which, l.by, l.kwargs, B)

function (l::GEigArpack)(J, nev; kwargs...)
	@show typeof(J) typeof(l.B)
	if J isa AbstractMatrix
		λ, ϕ, ncv = Arpack.eigs(J, l.B; nev = nev, which = l.which, sigma = l.sigma, l.kwargs...)
	else
		if !(:v0 in keys(l.kwargs))
			error("The v0 argument must be provided in EigArpack for the matrix-free case")
		end
		N = length(l.kwargs[:v0])
		T = eltype(l.kwargs[:v0])
		Jmap = LinearMap{T}(J, N, N; ismutating = false)
		λ, ϕ, ncv, = Arpack.eigs(Jmap, l.B; nev = nev, which = l.which, sigma = l.sigma, l.kwargs...)
	end
	Ind = sortperm(λ, by = l.by, rev = true)
	ncv < nev && @warn "$ncv eigenvalues have converged using Arpack.eigs, you requested $nev"
	return λ[Ind], ϕ[:, Ind], true, 1
end


text:
using CSV
using DataFrames
using FileIO
using Base.Iterators
using Random
using Flux: onehotbatch

function build_batches(metadatapath::AbstractString, melspectrogramspath::AbstractString, batchsize::Integer; eos='~', pad='_', kwargs...)
   dataset, alphabet = build_dataset(metadatapath, melspectrogramspath; eos=eos, pad=pad)
   batches = build_batches!(dataset, alphabet, batchsize, alphabet[eos], alphabet[pad], kwargs...)
   return batches, alphabet
end

build_batches(dataset::AbstractVector{<:Tuple{AbstractString,DenseMatrix{Float32}}}, alphabet::AbstractDict{Char,<:Integer}, batchsize::Integer, eosindex::Integer, padindex::Integer; kwargs...) =
   build_batches!(copy(dataset), alphabet, batchsize, eosindex, padindex; kwargs...)

function build_batches!(dataset::AbstractVector{<:Tuple{AbstractString,DenseMatrix{Float32}}},
                        alphabet::AbstractDict{Char,<:Integer},
                        batchsize::Integer,
                        eosindex::Integer,
                        padindex::Integer;
                        sorted=false, reorder=true)
   sorted || sort!(dataset; by=length∘first)
   batches = map(partition(dataset, batchsize)) do batch
      pad_batch(batch, alphabet, eosindex, padindex; sorted=true, reorder=reorder)
   end
   reorder && shuffle!(batches)
   return batches
end

function pad_batch(batch::AbstractVector{<:Tuple{String,DenseMatrix{Float32}}}, alphabet::AbstractDict{Char,<:Integer}, eosindex::Integer, padindex::Integer; sorted=false, reorder=true)
   batchsize = length(batch)
   maxlength_x = 1 + (sorted ? (length∘first∘last)(batch) : maximum(length∘first, batch))
   reorder && shuffle!(batch)

   mellengths = map(xy -> size(last(xy), 1), batch)
   maxlength_y = maximum(mellengths)
   nchannels_y = size(last(first(batch)), 2)

   stoptarget = Float32.(permutedims(onehotbatch(mellengths, 1:maxlength_y)))
   textindices = fill(padindex, maxlength_x, batchsize)
   melspectrograms = zeros(Float32, maxlength_y, nchannels_y, batchsize)

   for (n, (text, melspectrogram)) ∈ enumerate(batch)
      melspectrograms[axes(melspectrogram,1),:,n] = melspectrogram
      for (i, char) ∈ enumerate(text)
         textindices[i,n] = alphabet[char]
      end
      textindices[length(text)+1,n] = eosindex
   end
   (time, nmelfeatures, batchsize) = size(melspectrograms)
   melspectrograms⁴ = reshape(melspectrograms, (time, 1, nmelfeatures, batchsize))
   return textindices, melspectrograms⁴, stoptarget
end

function build_dataset(metadatapath::AbstractString, melspectrogramspath::AbstractString; eos='~', pad='_')
   df = DataFrame(CSV.File(metadatapath; delim='|', header=[:file_name, :text, :text_normalized], quotechar='\\', escapechar='\\'))
   ## sanity check that dataframe was parsed correctly
   # lines = readlines(metadatapath)
   # all(i -> collect(df[i,:]) == split(lines[i], '|'), eachindex(lines))
   select!(df, Not(:text))
   melspectrograms = load(melspectrogramspath, "melspectrograms")
   dataset = map(eachrow(df)) do row
      row.text_normalized, melspectrograms[row.file_name]
   end
   alphabet = build_alphabet(df.text_normalized; eos='~', pad='_')
   return dataset, alphabet
end

function build_alphabet(texts::AbstractVector{<:AbstractString}; eos::Char, pad::Char)
   alphabet = sort(unique(reduce(*, texts)))
   # padding character will be used a lot, so make its index to be 2 and make eos char index to be 1
   errormessage(char) = char -> "Cannot add the $char character to the alphabet because it already contains it"
   (eos ∈ alphabet) && error(errormessage("eos"))
   (pad ∈ alphabet) && error(errormessage("pad"))
   pushfirst!(alphabet, eos, pad)
   return Dict(v => k for (k, v) ∈ enumerate(alphabet))
end


text:
using LarSurf
# Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

V,(VV,EV,FV,CV) = LarSurf.Lar.cuboidGrid([10,20,1],true)

GL.VIEW([
      GL.GLGrid(V,FV,GL.Point4d(1,1,1,0.1))
	GL.GLAxis(GL.Point3d(-1,-1,-1),GL.Point3d(1,1,1))
])

GL.VIEW([
	GL.GLGrid(V,VV,GL.COLORS[4])
	GL.GLGrid(V,EV,GL.COLORS[2])
	GL.GLGrid(V,FV,GL.Point4d(0,1,0,0.1))
	GL.GLAxis(GL.Point3d(-1,-1,-1),GL.Point3d(1,1,1))
])


text:
import Dolo

import Dolo: n_nodes, n_inodes, nodes, CachedDecisionRule
import Dolo: invert_jac

model = Dolo.yaml_import("examples/models/rbc_mc.yaml")
dp = Dolo.discretize(model.exogenous)

using StaticArrays
import OrderedCollections
import Dolang
fff = Dolang.FlatFunctionFactory(OrderedCollections.OrderedDict{Symbol,Union{Expr, Number, Symbol}}(:_k__0_=>:((1 - _delta_) * _k_m1_ + _i_m1_)), OrderedCollections.OrderedDict(
:m=>[:_z_m1_, :_z2_m1_],:s=>[:_k_m1_],:x=>[:_n_m1_, :_i_m1_],:M=>[:_z__0_, :_z2__0_],
:p=>[:_beta_, :_sigma_, :_eta_, :_chi_, :_delta_, :_alpha_, :_rho_, :_zbar_, :_sig_z_]),
 Symbol[:_k__0_], OrderedCollections.OrderedDict{Symbol,Union{Expr, Number, Symbol}}(), :transition)

fun_code = Dolang.gen_generated_gufun(fff)
ff = Dolo.eval(fun_code)

s0 = SVector(model.calibration[:states]...)
x0 = SVector(model.calibration[:controls]...)
m0 = SVector(model.calibration[:exogenous]...)
p = SVector(model.calibration[:parameters]...)

println( (s0,x0,m0,p) )

ff(m0, s0, x0, m0, p)

@time sol = Dolo.time_iteration(model, verbose=false, complementarities=false)
@time sol = Dolo.time_iteration(model, sol.dr, verbose=false, complementarities=false)

@time sol = Dolo.time_iteration(model, verbose=false, complementarities=true)

@time sol = Dolo.improved_time_iteration(model, verbose=true, complementarities=false, method=:gmres)
@time sol = Dolo.improved_time_iteration(model, verbose=true, complementarities=true, method=:gmres)

@time sol = Dolo.improved_time_iteration(model, verbose=true, complementarities=false, method=:iti)

@time sol = Dolo.improved_time_iteration(model, verbose=true, complementarities=false, method=:gmres)

@time solv = Dolo.evaluate_policy(model, sol.dr, verbose=false)

@time solv = Dolo.value_iteration(model, verbose=true)

using ProfileView

xl = Profile.init()

yl = (1000000, 0.1)

Profile.init(n=yl[1], delay=yl[2])
Profile.clear()
@profile solv = Dolo.value_iteration(model, verbose=true)
ProfileView.view()

Profile.clear()
@profile sol = Dolo.time_iteration(model, verbose=false, complementarities=true)
ProfileView.view()

Profile.clear()
@profile sol = Dolo.improved_time_iteration(model, verbose=false, complementarities=true, method=:gmres)
ProfileView.view()

Profile.clear()
@profile sol = Dolo.time_iteration_direct(model, verbose=false)
ProfileView.view()

@time sol = Dolo.improved_time_iteration(model, verbose=true, complementarities=true, method=:gmres)

sim = Dolo.simulate(model, sol.dr)
dp = Dolo.discretize(model.exogenous)
grid = model.grid

@time Dolo.evaluate_policy(model, dp, grid, sol.dr, verbose=false)

SVector(model.calibration[:parameters]...)

@time sol = Dolo.time_iteration(model,verbose=true,complementarities=true, verbose=true, dampen=0.1)

Dolo.improved_time_iteration(model, sol.dr,verbose=true)
@time Dolo.improved_time_iteration(model,verbose=true,complementarities=false, verbose=false)

function iidr(m,s)
    k = s[1]
    kss = 9.3549
    n = 0.33-0.00*(k-kss)
    i = 0.234+0.000*(k-kss)
    SVector(n,i)
end

cfdr = InitDR.CFunDR(iidr, dp)

Dolo.improved_time_iteration(model, verbose=true, complementarities=false)

Dolo.improved_time_iteration(model, cfdr, verbose=true, complementarities=false)

Dolo.improved_time_iteration(model, cfdr, verbose=true)

methods(InitDR.CFunDR)

typeof(iidr)
typeof(dp)<:Dolo.AbstractDiscretizedProcess

cfdr(1, SVector(s_ss...))

# model = Dolo.yaml_import("examples/models/sudden_stop.yaml")

@time dr = Dolo.improved_time_iteration(model, verbose=true, method=:gmres, complementarities=true)

function time_it(t=100)
    l = []
    for i=1:t
        sol = Dolo.improved_time_iteration(model, verbose=false, method=:gmres, complementarities=true)
        push!(l, sol)
    end
    return true
end
Profile.clear()
@profile time_it(10)
ProfileView.view()

@time dr = Dolo.improved_time_iteration(model, verbose=false, method=:gmres)

@time dr = Dolo.improved_time_iteration(model, verbose=false, method=:iti)


text:
incy = incx
incz = incx


text:
struct Beta{T<:Real} <: Prior
    pivots::Vector{T}
    weights::Vector{T}
    pdf::Vector{T}
    a::T
    b::T
    low::T
    high::T
end

function Beta(a::Real, b::Real; low::Real = 0, high::Real = 1)
    p, ω = gauss_legendre_25(low, high)
    pdf  = dbeta.(p, convert(Float64, a), convert(Float64, b))
    z    = (high - low)/2 * sum(pdf .* ω)
    pdf  = pdf ./ z # normalize
    return Beta{Float64}(
        convert(Vector{Float64}, p),
        convert(Vector{Float64}, ω),
        convert(Vector{Float64}, pdf),
        convert(Float64, a),
        convert(Float64, b),
        convert(Float64, low),
        convert(Float64, high)
    )
end

# convenience function to fit a and b
function Beta(;mean::Real = .5, sd::Real = sqrt(1/12))
    sd >= .5 ? error("sd must be strictly smaller than .5") : nothing
    # solve for a, b
    a = (-mean^3 + mean^2 - mean*sd^2) / sd^2
    b = (mean - 1) * (mean^2 - mean + sd^2) / sd^2
    Beta(a, b; low = 0, high = 1)
end

function pdf(p::T, prior::Beta{T})::T where {T<:Real}
    if !(prior.low <= p <= prior.high); return 0.0 end
    F(p) = cdf(Distributions.Beta(prior.a, prior.b), p)
    return pdf(Distributions.Beta(prior.a, prior.b), p) / ((F(prior.high) - F(prior.low)))
end

function cdf(p::T, prior::Beta{T})::T where {T<:Real}
    p < prior.low ? (return 0.0) : nothing
    p > prior.high ? (return 1.0) : nothing
    F(p) = cdf(Distributions.Beta(prior.a, prior.b), p)
    return (F(p) - F(prior.low)) / ((F(prior.high) - F(prior.low)))
end

function expectation(f::Function, prior::Beta{T})::T where {T<:Real}

    return (prior.high - prior.low)/2 * sum( prior.pdf .* f.(prior.pivots) .* prior.weights )
end

function mean(prior::Beta{T})::T where {T<:Real}

    return expectation(p -> p, prior)
end

function condition(prior::Beta{T}; low::T = prior.low, high::T = prior.high)::Beta{T} where {T<:Real}

    return Beta(prior.a, prior.b, low = max(low, prior.low), high = min(high, prior.high))
end

function update(prior::Beta{T}, x::TI, n::TI)::Beta{T} where {T<:Real,TI<:Integer}

    !(0 <= x <= n) ? error(@sprintf "invalid x=%i / n=%i " x n) : nothing
    return Beta(prior.a + x, prior.b + n - x; low = prior.low, high = prior.high)
end

function string(prior::Beta)
    if (prior.low != 0) | (prior.high != 1)
        @sprintf "Beta|[%.2f,%.2f](a=%.2f,b=%.2f)" prior.low prior.high prior.a prior.b
    else
        @sprintf "Beta(a=%.2f,b=%.2f)" prior.a prior.b
    end
end

struct BetaMixture{T<:Real} <: Prior
    ω::Vector{T}
    priors::Vector{Beta{T}}
end

BetaMixture(ω::T, priors::Vector{Beta{T}}) where {T<:Real} = BetaMixture{T}(ω, priors)

bounds(prior::BetaMixture) = [minimum(minimum.(bounds.(prior.priors))), maximum(maximum.(bounds.(prior.priors)))]

*(ω::T, prior::Beta{T}) where{T<:Real} = BetaMixture{T}([ω], [prior])
*(ω::T, prior::BetaMixture{T}) where{T<:Real} = BetaMixture{T}(ω.*prior.ω, prior.priors)
+(φ::BetaMixture{T}, η::BetaMixture{T}) where {T<:Real} = BetaMixture{T}(vcat(φ.ω, η.ω), vcat(φ.priors, η.priors))
is_proper(mprior::BetaMixture) = sum(mprior.ω) == 1

function pdf(p::T, mprior::BetaMixture{T})::T where {T<:Real}

    return sum( mprior.ω .* pdf.(p, mprior.priors) )
end

function pdf(x1::TI, x2::TI, mprior::BetaMixture{T}) where {T<:Real,TI<:Integer}

    return sum( mprior.ω .* pdf.(x1, x2, mprior.priors) )
end

function cdf(p::T, mprior::BetaMixture{T})::T where {T<:Real}

    return sum( mprior.ω .* cdf.(p, mprior.priors) )
end

function expectation(f::Function, mprior::BetaMixture{T})::T where {T<:Real}

    return sum( mprior.ω .* expectation.(f::Function, mprior.priors) )
end

function mean(mprior::BetaMixture{T})::T where {T<:Real}

    sum( mprior.ω .* mean.(mprior.priors) )
end

condition(mprior::BetaMixture; low::T1 = 0., high::T1 = 1.) where {T1<:Real, T2<:Real} =
    BetaMixture(mprior.ω,  condition.(mprior.priors; low = low, high = high))

update(mprior::BetaMixture{T}, x::TI, n::TI) where {T<:Real, TI<:Integer} = BetaMixture{T}(mprior.ω,  update.(mprior.priors, x, n))

function string(mprior::BetaMixture)
    n = length(mprior.priors)
    res = ""
    for i in 1:(n - 1)
        res *= (@sprintf "%.2f*" mprior.ω[i]) * string(mprior.priors[i]) * " + "
    end
    res *= (@sprintf "%.2f*" mprior.ω[n]) * string(mprior.priors[n])
    res
end


text:
module Perturbation
using Printf: @printf, @sprintf
#using QuantumLattices:
using ExactDiagonalization.FEDW: matrix, NCBasis
using ExactDiagonalization.FED: FockBasis
using LinearAlgebra: diagm,eigen,Hermitian

import ExactDiagonalization.FEDW: bkron,boplus,dimension
import QuantumLattices: id, Operators, Table
export ProjState,GSpinGenerator,jval,hamiltonianeff,SpinHeff,projstate, projstate!
#abstract type AbstractProjState end

"""
    ProjState{T1<:Real,T2<:Complex}(vals::Vector{T1},vecs::Matrix{T2},basis::NCBasis)
    Projected states contain eigenvalues, eigenvectors, and basis.
"""
struct ProjState{T1 <:Real, T2 <: Complex  } #<: AbstractProjState
    vals :: Vector{T1}
    vecs :: Matrix{T2}
    basis :: NCBasis
end

"""
    bkron(ps1::ProjState{T1,T2},ps2::ProjState{T1,T2}) where {T1<:Real,T2<:Complex} -> ProjState
Direct product of two `ProjState`s.
"""
function bkron(ps1::ProjState{T1,T2},ps2::ProjState{T1,T2}) where {T1<:Real,T2<:Complex}
    vals0 = [ps1.vals[i]+ps2.vals[j] for i =1:length(ps1.vals) for j=1:length(ps2.vals) ]
    vecs0 = kron(ps1.vecs,ps2.vecs)
    basis0 = bkron(ps1.basis,ps2.basis)
    return ProjState{T1,T2}(vals0,vecs0,basis0)
end

"""
    boplus(ps1::ProjState{T1,T2},ps2::ProjState{T1,T2}) where {T1<:Real,T2<:Complex} -> ProjState
Direct sum of two `ProjState`s.
"""
function boplus(ps1::ProjState{T1,T2},ps2::ProjState{T1,T2}) where {T1<:Real,T2<:Complex}
    vals = vcat(ps1.vals,ps2.vals)
    n = size(ps1.vecs)
    m = size(ps2.vecs)
    vecs = zeros(T2,n[1]+m[1],n[2]+m[2])
    vecs[1:n[1],1:n[2]]=ps1.vecs
    vecs[1+n[1]:end,1+n[2]:end]=ps2.vecs
    hcat(ps1.vecs,ps2.vecs)
    basis = boplus(ps1.basis,ps2.basis)
    return ProjState{T1,T2}(vals,vecs,basis)
end

"""
    projstate(h0::Operators,basis::NCBasis,table::Table;pick::Union{UnitRange{Int},Colon}=:) -> ProjState
Construct `ProjState`
"""
function projstate(h0::Operators,basis::NCBasis,table::Table;pick::Union{UnitRange{Int},Colon}=:)
    hm = matrix(h0,basis,table)
    hm2 = Hermitian(Array(hm+hm')/2)
    F = eigen(hm2)
    return ProjState{eltype(F.values),eltype(F.vectors)}(F.values[pick],F.vectors[:,pick],basis)
end
"""
    projstate!(prj::ProjState,U::T) where T<:AbstractMatrix{ComplexF64} -> ProjState
unitary transformation of eigenvectors.
"""
function projstate!(prj::ProjState,U::T) where T<:AbstractMatrix{ComplexF64}
    vecs=prj.vecs*U
    return ProjState(prj.vals,vecs,prj.basis)
end

"""
    GSpinGenerator{T<:Number,N} where N
Struct `GSpinGenerator` contains `names` (general spin generator names, e.g. spin generator σ⁰,σˣ,σʸ,σᶻ) and the corresponding matrix representation `mat` (e.g. Pauli matrix).
"""
struct GSpinGenerator{T<:Number,N}
    names::NTuple{N,Symbol}
    mat:: NTuple{N,Matrix{T}}
end
GSpinGenerator(names::NTuple{N,Symbol},mat::Tuple) where N = GSpinGenerator(names,promote(mat...))
"""
    Base.getindex(gsf::GSpinGenerator,name::Symbol) -> Matrix
"""
Base.getindex(gsf::GSpinGenerator,name::Symbol) = _getinx_(gsf,name)
function _getinx_(gsf::GSpinGenerator,name::Symbol)
    i=findfirst(isequal(name),gsf.names)
    return gsf.mat[i]
end
"""
    Base.getindex(gsf::GSpinGenerator,i::Int) -> Tuple{Symbol,Matrix}
Get the i-th (name::Symbol,matrix::Matrix{T}) of gsf::GSpinGenerator.
"""
Base.getindex(gsf::GSpinGenerator,i::Int) = (gsf.names[i],gsf.mat[i])
"""
    Base.keys(gsf::GSpinGenerator) -> NTuple{N,Symbol}
Get the names of the general spin generators.
"""
Base.keys(gsf::GSpinGenerator) = gsf.names
Base.values(gsf::GSpinGenerator) = gsf.mat
Base.eltype(gsf::GSpinGenerator) = typeof(gsf) |> eltype
Base.eltype(::Type{GSpinGenerator{T,N}}) where {N,T<:Number} = T
"""
    dimension(gsf::GSpinGenerator) -> Tuple{Int,Int}
Get the size of matrix of the generator.
"""
dimension(gsg::GSpinGenerator) = size( gsg.mat[1] )
"""
    Base.length(gsg::GSpinGenerator) -> Int
The number of the general spin generators.
"""
Base.length(gsg::GSpinGenerator) = Base.length(gsg.mat)
"""
    Base.kron(gsg1::GSpinGenerator,gsg2::GSpinGenerator) -> GSpinGenerator
Direct product of two `GSpinGenerator`s.
"""
function Base.kron(gsg1::GSpinGenerator,gsg2::GSpinGenerator)
    temp=tuple([kron(gsg1[i],gsg2[j]) for i in keys(gsg1) for j in keys(gsg2)]... )
    nameS = tuple([Symbol(i,",",j)  for i in keys(gsg1) for j in keys(gsg2)]... )
    return GSpinGenerator(nameS,temp)
end

abstract type AbstractHeff{T1<:Number,N} end

"""
    SpinHeff{T1,N} <: AbstractHeff{T1,N}
    SpinHeff{T1,N}( hterm :: NTuple{N,Matrix{T1}}, jeff :: Vector{T1}, gs :: GSpinGenerator) where T1<:Number,N
To store the matrix representation of the effective Hamiltonian, the exchange coefficients and the general spin generators.
"""
struct SpinHeff{T1,N} <: AbstractHeff{T1,N}
    hterm :: NTuple{N,Matrix{T1}}
    jeff :: Vector{T1}
    gs :: GSpinGenerator
end
"""
    Base.getindex(sph::SpinHeff,i::Int) -> Matrix
Get the i-th term of the effective Hamiltonian, e.g. sph[1] = the zero-order perturbation term P*H_0*P
"""
Base.getindex(sph::SpinHeff,i::Int) = sph.hterm[i]
Base.keys(sph::SpinHeff) = sph.gs.names
Base.values(sph::SpinHeff) = sph.jeff
Base.getindex(sph::SpinHeff,name::Symbol) = _getinx_(sph,name)
function _getinx_(sph::SpinHeff,name::Symbol)
    i=findfirst(isequal(name),sph.gs.names)
    return sph.jeff[i]
end
id(sph::SpinHeff) = ["hterm",length(sph.hterm),"jeff",length(sph.jeff),"gs"]

function Base.show(io::IO,sph::SpinHeff)
    @printf io "%s:\n" repr(sph)
    @assert length(sph.jeff) == length(sph.gs.names) " show error: the number of jeff ($(length(sph.jeff))) == the number of SpinHeff.gs.names ($(length(sph.gs.names)))"
    for i = 1:length(sph.jeff)
        @printf io " \$J_{%s}\$ = %s\n" sph.gs.names[i] sph.jeff[i]
    end
end
@inline Base.repr(sph::SpinHeff) = @sprintf "%s(%s)" nameof(typeof(sph)) join(id(sph), ",")

"""
    jval(heffmat::Matrix{<:Number},gsg::GSpinGenerator) -> Vector
Get the coefficients of exchange terms.
"""
function jval(heffmat::Matrix{<:Number},gsg::GSpinGenerator)
    b=heffmat[:]
    nn = length(gsg)
    @assert nn == length(b) "jval error: length(GSpinGenerator)[$(length(gsg))] == length(heffmat[:])[$(length(b))]."
    @assert nn == dimension(gsg)|>prod "jval error: length(GSpinGenerator)[$(length(gsg))] == the number of element of matrix GSpinGenerator.mat[$(prod(dimension(gsg)))]."
    a=zeros(eltype(gsg),nn,nn)
    for i = 1:nn
        a[:,i] = gsg.mat[i][:]
    end
    return a\b
end

"""
    jval(heff::NTuple{N,Matrix{T1}},gsg::GSpinGenerator) where {N,T1<:Number} -> SpinHeff{T1,N}
Get the struct `SpinHeff`.
"""
function jval(heff::NTuple{N,Matrix{T1}},gsg::GSpinGenerator) where {N,T1<:Number}
    result = zeros(T1,length(gsg))
    for i =1:length(heff)
        result = result+jval(heff[i],gsg)
    end
    return SpinHeff{T1,N}(heff,result,gsg)
end

"""
    hamiltonianeff(psp::ProjState, psq::ProjState,h1::Operators,table1::Table;order::Int=2) ->Tuple{Matrix,Matrix}
Get the effective Hamiltonian, the first and second terms of the result correspond to the zero-th and 2nd perturbations respectively.
"""
function hamiltonianeff(psp::ProjState, psq::ProjState,h1::Operators,table1::Table;order::Int=2)
    @assert order == 2 "hamiltonianeff error: now only support the second order perturbation."
    heff0 = diagm(psp.vals)
    heff01 = psp.vecs'*Array(matrix(h1,psp.basis,table1))*psp.vecs

    tqp = psq.vecs'*Array(matrix(h1,psq.basis,psp.basis,table1))*psp.vecs
    m, n = size(psp.vecs,2), size(psq.vecs,2)
    sqp = zeros(eltype(psp.vecs),n,m)
    for i = 1:m
        for j = 1:n
            sqp[j,i] = -tqp[j,i]/(-psp.vals[i]+psq.vals[j])
        end
    end
    spq = - sqp'
    tpq = tqp'
    heff2 = (tpq*sqp - spq*tqp )/2.0
    return heff0+heff01, heff2
end

end #module


text:
using Test
using ConvDiffMIPDECO
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare

println("\t\t-- test ConvDiffParam --")
domain = [-1 1. -1 1]
n      = [32 32]
M      = getRegularMesh(domain,n)
bc     = (:dir,:neu,:dir,:neu)
h      = M.h
v      = X -> [ones(size(X,1)) zeros(size(X,1))]
gd     = X -> X[:,1].*X[:,2]
gn     = X -> ones(size(X,1))
pFor   = getConvDiffParam(M,v,sig=.5,gd=gd,gn=gn,bc=bc)

m0     = zeros(tuple(n...))
m0[10:20,15:25] .=1.;
isOK, his = checkDerivative(vec(m0),pFor,out=true)
@test isOK
v = randn(M.nc)
w = randn(prod(M.nc))
t1 = dot(w, getSensMatVec(v,vec(m0),pFor))
t2 = dot(v, getSensTMatVec(w,vec(m0),pFor))
@test abs(t1-t2)/abs(t1) < 1e-8

domain = [-1 1. 0 1 -1 1]
n      = [6 8 10]
M      = getRegularMesh(domain,n)
bc     = (:dir,:neu,:dir,:neu,:dir,:neu)
h      = M.h
v      = X -> [ones(size(X,1)) zeros(size(X,1)) zeros(size(X,1))]
gd     = X -> X[:,1].*X[:,2]
gn     = X -> ones(size(X,1))
pFor   = getConvDiffParam(M,v,sig=.5,bc=bc)

m0              = zeros(tuple(n...))
m0[3:4,4:6,4:7] .=1.;
isOK, his = checkDerivative(vec(m0),pFor,out=true)
@test isOK
v = randn(M.nc)
w = randn(prod(M.nc))
t1 = dot(w, getSensMatVec(v,vec(m0),pFor))
t2 = dot(v, getSensTMatVec(w,vec(m0),pFor))
@test abs(t1-t2)/abs(t1) < 1e-8

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
## 3PL
"""
    Parameters3PL <: AbstractParametersBinary

Contains information of a set of item parameters (values, bounds, priors, posteiors, chains, expected Fisher information, calibrated) under the 3-parameter logistic model.
"""
mutable struct Parameters3PL <: AbstractParametersBinary
    a::Float64
    bounds_a::Vector{Float64}
    b::Float64
    bounds_b::Vector{Float64}
    c::Float64
    bounds_c::Vector{Float64}
    prior::Distributions.MultivariateDistribution
    posterior::Distributions.MultivariateDistribution
    chain::Vector{Vector{Float64}}
    expected_information::Matrix{Float64}
    calibrated::Bool

    Parameters3PL(
        a,
        bounds_a,
        b,
        bounds_b,
        c,
        bounds_c,
        prior,
        posterior,
        chain,
        expected_information,
        calibrated,
    ) = new(
        a,
        bounds_a,
        b,
        bounds_b,
        c,
        bounds_c,
        prior,
        posterior,
        chain,
        expected_information,
        calibrated,
    )

    # Random Initializers

    function Parameters3PL(
        trivariate_dist::Distributions.MultivariateDistribution,
        bounds_a,
        bounds_b,
        bounds_c,
    )
        pars = truncate_rand(trivariate_dist, [bounds_a, bounds_b, bounds_c])
        new(
            pars[1][1],
            bounds_a,
            pars[1][2],
            bounds_b,
            pars[1][3],
            bounds_c,
            trivariate_dist,
            trivariate_dist,
            [pars[1]],
            [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
            false,
        )
    end

    function Parameters3PL()
        Parameters3PL(
            Distributions.Product([
                Distributions.LogNormal(0, 0.25),
                Distributions.Normal(0, 1),
                Distributions.Uniform(0, 1),
            ]),
            [1e-5, 5.0],
            [-6.0, 6.0],
            [1e-5, 1.0 - 1e-5],
        )
    end
end

"""
    _empty_chain!(parameters::Parameters3PL)
"""
function _empty_chain!(parameters::Parameters3PL)
    parameters.chain = Vector{Vector{Float64}}(undef, 0)
end

"""
    _chain_append!(parameters::Union{Parameters2PL,Parameters3PL}; sampling = false)
"""
function _chain_append!(parameters::Union{Parameters2PL,Parameters3PL}; sampling = false)
    val = Distributions.rand(parameters.posterior)
    if (sampling && size(parameters.chain, 1) >= 1000)
        parameters.chain[Random.rand(1:1000)] = val
    else
        push!(parameters.chain, val)
    end
    return nothing
end

"""
    _chain_append_and_set_val!(parameters::Parameters3PL; sampling = false)
"""
function _chain_append_and_set_val!(parameters::Parameters3PL; sampling = false)
    val = Distributions.rand(parameters.posterior)
    if (sampling && size(parameters.chain, 1) >= 1000)
        parameters.chain[Random.rand(1:1000)] = val
    else
        push!(parameters.chain, val)
    end
    parameters.a = val[1]
    parameters.b = val[2]
    parameters.c = val[3]
    return nothing
end

"""
    _chain_append_and_set_val!(parameters::Parameters2PL; sampling = false)
"""
function _chain_append_and_set_val!(parameters::Parameters2PL; sampling = false)
    val = Distributions.rand(parameters.posterior)
    if (sampling && size(parameters.chain, 1) >= 1000)
        parameters.chain[Random.rand(1:1000)] = val
    else
        push!(parameters.chain, val)
    end
    parameters.a = val[1]
    parameters.b = val[2]
    return nothing
end


text:
# # Poisson Equation
#
# $$
# \frac{\partial^2 u}{\partial x^2} = b  \qquad x \in [0,1]
# $$
#
# $$
# u(0) = u(1) = 0, \qquad b = \sin(2\pi x)
# $$

using Plots, SparseArrays

Δx = 0.05
x = Δx:Δx:1-Δx ## Solve only interior points: the endpoints are set to zero.
n = length(x)
B = sin.(2π*x) * Δx^2

P = spdiagm( -1 =>    ones(Float64,n-1),
              0 => -2*ones(Float64,n),
              1 =>    ones(Float64,n-1))

u1 = P \ B

# ---

plot([0;x;1],[0;u1;0], label="computed")
scatter!([0;x;1],-sin.(2π*[0;x;1])/(4π^2),label="exact")
savefig("poisson1.png"); nothing #hide

# ![](poisson1.png)

# ---

# # DiffEqOperators.jl

using DiffEqOperators

Δx = 0.05
x = Δx:Δx:1-Δx ## Solve only interior points: the endpoints are set to zero.
n = length(x)
b = sin.(2π*x)

## Second order approximation to the second derivative
order = 2
deriv = 2

Δ = CenteredDifference{Float64}(deriv, order, Δx, n)
bc = Dirichlet0BC(Float64)

u2 = (Δ * bc) \ b

# ---

plot([0;x;1],[0;u2;0], label="computed")
scatter!([0;x;1],-sin.(2π*[0;x;1])/(4π^2),label="exact")
savefig("poisson2.png"); nothing #hide

# ![](poisson2.png)

# ---


text:
# function compose_core(X, B, rdef)
#    nX = size(X,2)
#    r = sumabs2(X, 1) |> sqrt
#    nneigs = zeros(Int, nX)
#    for n = 1:length(B[1])
#       nneigs[B[1][n]] += 1
#       nneigs[B[2][n]] += 1
#    end
#    Icore = find( (nneigs .!= 6) .* (r .<= rdef) )
# end

# ########################## PLOTTING ############################

# # code to plot in a REPL window
# using Gtk
# c = Gtk.@GtkCanvas(400,300);
# w = Gtk.@GtkWindow(c,"data win");
# show(c);
# using Compose
# function sierpinski(n)
#     if n == 0
#         compose(context(), polygon([(1,1), (0,1), (1/2, 0)]))
#     else
#         t = sierpinski(n - 1)
#         compose(context(),
#                 (context(1/4,   0, 1/2, 1/2), t),
#                 (context(  0, 1/2, 1/2, 1/2), t),
#                 (context(1/2, 1/2, 1/2, 1/2), t))
#     end
# end
# co = sierpinski(5);
# Gtk.draw(c) do widget
#     Compose.draw(CAIROSURFACE(c.back),co)
# end
# Gtk.draw(c)


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

using Markdown
using InteractiveUtils

# ╔═╡ 6f915d8a-a7d6-11ec-0039-cbbff20d0352
begin
	using PlutoUI

md"""
# HW 9

Eric Nguyen
"""
end

# ╔═╡ 5a928acb-5f63-4429-9bc6-b03e1ea75911
md"""
### Problem 1

The wave equation is defined as

$\frac{1}{v^2} \frac{∂^2 u}{∂t^2} = \frac{∂^2 u}{∂x^2}$

The general solution to the wave equation is of the form

$u(x,t) = f(x - vt) + g(x + vt)$

for any differentiable functions ``f`` and ``g`` and constant velocity ``v``.
We want to show that

$u(x,t) = \frac{1}{2} p(x - vt) + \frac{1}{2} p(x + vt)$

The initial condition ``u(x,0) = p(x)`` gives

$u(x,0) = f(x) + g(x) = p(x) \tag{1}$

The initial condition ``u̇(x,0) = 0`` gives

$u̇(x,0) = -v \dot{f}(x) + v \dot{g}(x) = 0$

Integrating the expressions in the second initial condition with respect to ``t``, we get

$u(x,0) = -v f(x) + v g(x) = 0$

which reduces to

$u(x,0) = -f(x) + g(x) = 0 \tag{2}$

Subtracting and adding the equations ``(1)`` and ``(2)`` together, we get

$\begin{align*}
2f(x) &= p(x) \implies f(x) = \frac{1}{2} p(x) \implies f(x - vt) = \frac{1}{2} p(x - vt) \\
2g(x) &= p(x) \implies g(x) = \frac{1}{2} p(x) \implies g(x + vt) = \frac{1}{2} p(x + vt)
\end{align*}$

Finally, we can rewrite the solution as

$u(x,t) = \frac{1}{2} p(x - vt) + \frac{1}{2} p(x + vt)$

which indeed corresponds to the "splitting" of ``p(x)`` into two pieces, one moving to the left and the other moving to the right such that each piece is equal to one-half of ``p(x)``.
"""

# ╔═╡ b1b35d3c-8aff-4210-b7e5-26af398f7667
md"""
### Problem 2

**(a).**
We obtain a general solution to the wave equation such that the condition ``u(0,t) = 0`` is satisfied but still is of the form ``u(x,t) = g(x + vt)``.
We get

$u(x,t) = g(x + vt) - g(-x + vt)$

Indeed, this satisfies the condition ``u(0,t) = 0``,

$u(0,t) = g(vt) - g(vt) = 0$

**(b).**
We obtain a general solution to the wave equation such that the condition ``\displaystyle \left.\frac{∂u(x,t)}{∂x}\right|_{x=0} = 0`` is satisfied but is still of the form ``u(x,t) = g(x + vt)``.
We get

$u(x,t) = g(x + vt) + g(-x + vt)$

Indeed, this satisfies the condition ``\displaystyle \left.\frac{∂u(x,t)}{∂x}\right|_{x=0} = 0``,

$\displaystyle \left.\frac{∂u(x,t)}{∂x}\right|_{x=0} = \frac{∂}{∂x} (g(x + vt) + g(-x + vt)) = 1 - 1 = 0$
"""

# ╔═╡ 29ba7a28-2fdd-4e77-96a5-160578539b7a
md"""
### Problem 3

The Fourier Sine decomposition shows a similar motion compared to Lab 8.
"""

# ╔═╡ 244e298f-393a-4a76-b287-f8f1fcd9b94d
PlutoUI.LocalResource("img/hw_9_p3.png")

# ╔═╡ 701dc659-73d1-4264-a179-82d6c61eeaf2
md"""
### Problem 4

$f(x) = \frac{1}{\sigma \sqrt{2\pi}} e^{-x^2 / 2\sigma^2}$

**(a).**
The Fourier Transform ``A(k)`` of ``f(x)`` is

$\begin{align*}
A(k) &= ∫_{-∞}^∞ e^{-2\pi ikx} f(x) \,dx \\
&= \frac{1}{\sigma \sqrt{2\pi}} ∫_{-∞}^∞ e^{-2\pi ikx} e^{-x^2 / 2\sigma^2} \,dx \\
&= \frac{1}{\sigma \sqrt{2\pi}} ∫_{-∞}^∞ e^{-(4 \pi \sigma^2 ikx + x^2) / 2\sigma^2} \,dx \\
&= \frac{1}{\sigma \sqrt{2\pi}} ∫_{-∞}^∞ e^{-[(x + 2\pi \sigma^2 ik)^2 - 4 \pi^2 \sigma^4 i^2 k^2] / 2\sigma^2} \,dx & \text{Complete the square} \\
&= \frac{e^{2 \pi^2 \sigma^2 i^2 k^2}}{\sigma \sqrt{2\pi}} ∫_{-∞}^∞ e^{-(x + 2\pi^2 \sigma^2 ik)^2 / 2\sigma^2} \,dx \\
&= \frac{e^{-2 \pi^2 \sigma^2 k^2}}{\sigma \sqrt{2\pi}} \left[\sqrt{2\pi \sigma^2} \right] & \text{Using } ∫_{-∞}^∞ e^{-a(x+b)^2} \,dx = \sqrt{\frac{\pi}{a}} \\
&= e^{-2\pi^2 \sigma^2 k^2}
\end{align*}$

**(b).**
The width ``Δx`` of ``f(x)`` is found through the following calculations

$\mathcal{A} = ∫_{-∞}^∞ f(x) \,dx = \frac{1}{\sqrt{2 \pi}} ∫_{-∞}^∞ e^{-x^2/2\sigma^2} \,dx = \frac{1}{\sqrt{2\pi}} \left[\sqrt{2\pi\sigma^2}\right] = \sigma$

$\begin{align*}
(Δx)^2 &= \frac{1}{\mathcal{A}} ∫_{-∞}^∞ x^2 f(x) \,dx \\
&= \frac{1}{\sigma^2 \sqrt{2 \pi}} ∫_{-∞}^∞ x^2 e^{-x^2 / 2\sigma^2} \,dx \\
&= \frac{1}{\sigma^2 \sqrt{2 \pi}} \left[\frac{1}{2} \sqrt{\frac{\pi}{(1/2\sigma^2)^3}}\right] \\
&= \frac{1}{\sigma^2 \sqrt{2 \pi}} \left[\frac{1}{2} \sqrt{8\pi \sigma^6}\right] \\
&= \frac{1}{\sigma^2 \sqrt{2 \pi}} \left[\frac{1}{2} \left(2 \sqrt{2 \pi} \sigma^3\right)\right] \\
&= \sigma \\
\end{align*}$

$\implies Δx = \sqrt{\sigma}$

**(c).**
The width ``Δk`` of ``A(k)`` is found through the following calculations

$\mathcal{A} = ∫_{-∞}^∞ A(k) \,dk = ∫_{-∞}^∞ e^{-2\pi^2 \sigma^2 k^2} \,dk = \sqrt{\frac{\pi}{2\pi^2 \sigma^2}}$

$\begin{align*}
(Δk)^2 &= \frac{1}{\mathcal{A}} ∫_{∞}^∞ k^2 A(k) \,dk \\
&= \sqrt{\frac{2\pi^2 \sigma^2}{\pi}} ∫_{∞}^∞ k^2 e^{-2\pi^2 \sigma^2 k^2} \,dk \\
&= \sqrt{\frac{2\pi^2 \sigma^2}{\pi}} \left[\frac{1}{2} \sqrt{\frac{\pi}{(2\pi^2 \sigma^2)^3}}\right] \\
&= \frac{1}{2 (2 \pi^2 \sigma^2)} \\
&= \frac{1}{4 \pi^2 \sigma^2}
\end{align*}$

$\implies Δk = \sqrt{\frac{1}{4\pi^2 \sigma^2}} = \frac{1}{2\pi \sigma}$

**(d).**
The product ``Δx Δk`` is

$Δx Δk = \frac{1}{4\pi^2 \sqrt{\sigma}}$

which implies that it falls to zero even faster than the triangle pulse for ``\displaystyle \sigma ≥ \frac{1}{4 \pi^6}``.
"""

# ╔═╡ 1767eedf-c523-4865-956b-bd9ba83cc9cf
md"""
### Problem 5

**(a).**
Let ``f(x) = x`` for the property

$\int_{-\epsilon}^\epsilon f(x) \delta(x) \,dx = f(0)$

We have ``f(0) = 0``, so we get

$x \delta(x) = 0$

**(b).**

$∫_{-∞}^∞ x \delta'(x) \,dx = x \delta(x)\Big|_{-∞}^∞ - ∫_{-∞}^∞ \delta(x) \,dx = -∫_{-∞}^∞ \delta(x) \,dx$

$\implies x \delta'(x) = -\delta(x)$

**(c).**

$∫_{-∞}^∞ x^2 \delta''(x) \,dx = \left[2x \delta'(x)\right]_{-∞}^∞ - ∫_{-∞}^{∞} 2x \delta'(x) \,dx = ∫_{-∞}^∞ 2\delta(x) \,dx$

$\implies x^2 \delta''(x) \,dx = 2\delta(x)$
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.37"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.8.0-beta1"
manifest_format = "2.0"
project_hash = "e4d36aff991b8dc605f222074120cb6c69173dc8"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.5.0+0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
version = "1.6.0"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
version = "0.6.3"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
version = "7.81.0+0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
version = "1.10.2+0"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
version = "2.28.0+0"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
version = "2022.2.1"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
version = "1.2.0"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
version = "0.3.17+2"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
version = "1.8.0"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "bf0a1121af131d9974241ba53f601211e9303a9e"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.37"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
version = "0.7.0"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
version = "1.0.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
version = "1.10.0"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
version = "1.2.12+1"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
version = "5.0.1+0"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
version = "1.41.0+1"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
version = "16.2.1+1"
"""

# ╔═╡ Cell order:
# ╟─6f915d8a-a7d6-11ec-0039-cbbff20d0352
# ╟─5a928acb-5f63-4429-9bc6-b03e1ea75911
# ╟─b1b35d3c-8aff-4210-b7e5-26af398f7667
# ╟─29ba7a28-2fdd-4e77-96a5-160578539b7a
# ╟─244e298f-393a-4a76-b287-f8f1fcd9b94d
# ╟─701dc659-73d1-4264-a179-82d6c61eeaf2
# ╟─1767eedf-c523-4865-956b-bd9ba83cc9cf
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002