text:
# Autogenerated wrapper script for LibUV_jll for armv7l-linux-gnueabihf
export libuv

JLLWrappers.@generate_wrapper_header("LibUV")
JLLWrappers.@declare_library_product(libuv, "libuv.so.2")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libuv,
        "lib/libuv.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function compute_ldb_vectors(ldb::LocalDiscriminantBasis, N::Integer)
    n = ldb.sz[1]                               # Signal length
    topN_subspace = ldb.order[1:N]              # Position of top N subspace
    vectors = zeros(n,N)                        # Build matrix to containing top N subspace
    for (j, i) in enumerate(topN_subspace)
        vectors[i,j] = 1
    end
    vectors = iwptall(vectors, wt, ldb.tree)    # Inverse transform matrix
    return vectors
end

function plot_coefficients(data)
    p = scatter(data[1,1:33], data[2,1:33], shape=:circle, color=:black, label="Class 1")
    scatter!(p, data[1,34:66], data[2,34:66], shape=:star5, color=:yellow, label="Class 2")
    scatter!(p, data[1,67:99], data[2,67:99], shape=:cross, color=:green, label="Class 3")
    return p
end

"""

"""
function extract_feature!(ldb::Union{LocalDiscriminantBasis, Nothing},
                          X::AbstractArray{T₁}, y::AbstractVector{T₂};
                          return_plots::Bool = false,
                          N::Integer = 3) where
                         {T₁<:AbstractFloat, T₂}
    # If no available LDB object, then plots shouldn't be returned
    @assert !(isnothing(ldb) && return_plots)
    # Fit LDB
    train_features = isnothing(ldb) ? X : WaveletsExt.fit_transform(ldb, X, y)
    # Build plots if necessary
    if return_plots
        ldbvec = compute_ldb_vectors(ldb, N)
        ldbvec_plot = plot(ldbvec, label="")
        tfbdry_plot = plot_tfbdry(ldb.tree)
        coef_plot = plot_coefficients(train_features)
        return (ldbvec_plot = ldbvec_plot, tfbdry_plot = tfbdry_plot, coef_plot = coef_plot)
    end
end

function fit_model(ldb::Union{LocalDiscriminantBasis, Nothing},
                   classifier::Supervised,
                   X::AbstractArray{T₁}, y::AbstractVector{T₂}) where
                  {T₁<:AbstractFloat, T₂}
    # Get LDB features
    train_features = isnothing(ldb) ? X : WaveletsExt.transform(ldb, X)
    # Model fitting
    Xₜ = DataFrame(train_features', :auto)
    yₜ = coerce(y, Multiclass)
    model = machine(classifier, Xₜ, yₜ)
    MLJ.fit!(model)
    return model
end

function evaluate_model(ldb::Union{LocalDiscriminantBasis, Nothing},
                        model::Machine,
                        X::AbstractArray{T₁}, y::AbstractVector{T₂},
                        measures::Vector{T₃}) where
                       {T₁<:AbstractFloat, T₂, T₃<:MLJ.MLJBase.Measure}
    # Get LDB features
    train_features = isnothing(ldb) ? X : WaveletsExt.transform(ldb, X)
    # Model prediction
    Xₜ = DataFrame(train_features', :auto)
    yₜ = coerce(y, Multiclass)
    ŷ = predict_mode(model, Xₜ)
    # Model evaluation
    result = Dict{String, Float64}()
    for measure in measures
        metric = typeof(measure)
        result["$metric"] = measure(ŷ, yₜ)
    end
    return result
end

function vector2dict(vec::Vector{T₁}, suffix::String = "method") where T₁
    return Dict("$(suffix)_$i" => data for (i, data) in enumerate(vec))
end

function run_experiment(ldbs::Union{Dict{String, T₁}, Vector{T₁}},
                        classifiers::Union{Dict{String, T₂}, Vector{T₂}},
                        X_train::AbstractArray{T₃}, y_train::AbstractVector{T₄},
                        X_test::AbstractArray{T₃}, y_test::AbstractVector{T₄},
                        measures::Vector{T₅}; kwargs...) where
                       {T₁<:Union{LocalDiscriminantBasis, Nothing},
                        T₂<:Supervised,
                        T₃<:AbstractFloat, T₄,
                        T₅<:MLJ.MLJBase.Measure}
    ldbs = isa(ldbs, Vector) ? vector2dict(ldbs, "ldb") : ldbs
    classifiers = isa(classifiers, Vector) ? vector2dict(classifiers, "clf") : classifiers
    models = Dict(name => Dict{String, Machine}() for name in keys(ldbs))
    result = Dict(name => Dict{String, Dict}() for name in keys(ldbs))
    for (ldb_name, ldb) in ldbs
        for (clf_name, classifier) in classifiers
            # Extract LDB features
            extract_feature!(ldb, X_train, y_train)
            # Fit model
            models[ldb_name][clf_name] = fit_model(ldb, classifier, X_train, y_train)
            # Evaluate model
            train_result = evaluate_model(ldb, models[ldb_name][clf_name], X_train, y_train, measures)
            test_result = evaluate_model(ldb, models[ldb_name][clf_name], X_test, y_test, measures)
            result[ldb_name][clf_name] = Dict("train" => train_result, "test" => test_result)
        end
    end
    return (model = models, result = result)
end

function repeat_experiment(ldbs::Union{Dict{String, T₁}, Vector{T₁}},
                           classifiers::Union{Dict{String, T₂}, Vector{T₂}},
                           measures::Vector{T₃};
                           config_train::ClassData = ClassData(:cbf, 33, 33, 33),
                           config_test::ClassData = ClassData(:cbf, 333, 333, 333),
                           repeats::Integer = 10, kwargs...) where
                          {T₁<:Union{LocalDiscriminantBasis, Nothing},
                           T₂<:Supervised,
                           T₃<:MLJ.MLJBase.Measure}
    ldbs = isa(ldbs, Vector) ? vector2dict(ldbs, "ldb") : ldbs
    classifiers = isa(classifiers, Vector) ? vector2dict(classifiers, "clf") : classifiers
    results = Dict{String, Dict}()
    for i in 1:repeats
        X_train, y_train = generateclassdata(config_train, false)
        X_test, y_test = generateclassdata(config_test, true)
        _, resultᵢ = run_experiment(ldbs, classifiers, X_train, y_train, X_test, y_test, measures, kwargs...)
        results["experiment $i"] = resultᵢ
    end
    return results
end

h₁(i::Int) = max(6 - abs(i-7), 0)
h₂(i::Int) = h₁(i - 8)
h₃(i::Int) = h₁(i - 4)

"""
  triangular_test_function(c1::Int, c2::Int, c3::Int, L::Int=32)

Generates a set of triangluar test functions with 3 classes.
"""
function triangular_test_functions(c1::Int, c2::Int, c3::Int; L::Int=32, shuffle::Bool=false)
  @assert c1 >= 0
  @assert c2 >= 0
  @assert c3 >= 0

  u = rand(Uniform(0,1),1)[1]
  ϵ = rand(Normal(0,1),(L,c1+c2+c3))

  y = string.(vcat(ones(c1), ones(c2) .+ 1, ones(c3) .+ 2))

  H₁ = Array{Float64,2}(undef,L,c1)
  H₂ = Array{Float64,2}(undef,L,c2)
  H₃ = Array{Float64,2}(undef,L,c3)
  for i in 1:L
    H₁[i,:] .= u * h₁(i) + (1 - u) * h₂(i)
    H₂[i,:] .= u * h₁(i) + (1 - u) * h₃(i)
    H₃[i,:] .= u * h₂(i) + (1 - u) * h₃(i)
  end

  H = hcat(H₁, H₂, H₃) + ϵ

  if shuffle
    idx = [1:(c1+c2+c3)...]
    shuffle!(idx)
    return H[:,idx], y[idx]
  end

  return H, y
end

ext: jl


text:
parse_projection(m::Number) = m
function parse_projection(m::AbstractString)
    n = tryparse(Float64, m)
    isnothing(n) || return n
    if m == "α"
        half(1)
    elseif m == "β"
        -half(1)
    elseif occursin('/', m)
        n,d = split(m, "/")
        parse(Int,n)//parse(Int,d)
    else
        throw(ArgumentError("Don't know how to parse projection quantum number $(m)"))
    end
end

"""
    struct SpinOrbital{O<:Orbital} <: AbstractOrbital

Spin orbitals are fully characterized orbitals, i.e. the projections
of all angular momenta are specified.
"""
struct SpinOrbital{O<:AbstractOrbital,M<:Tuple} <: AbstractOrbital
    orb::O
    m::M
    function SpinOrbital(orb::O, m::Tuple) where O
        O <: SpinOrbital &&
            throw(ArgumentError("Cannot create SpinOrbital from $O"))

        am = angular_momentum_ranges(orb)
        aml = angular_momentum_labels(orb)
        nam = length(am)
        nm = length(m)
        nm == nam ||
            throw(ArgumentError("$(nam) projection quantum numbers required, got $(nm)"))
        newm = ()
        for (mi,ar,l) in zip(m,am,aml)
            mm = parse_projection(mi)
            mm ∈ ar ||
                throw(ArgumentError("Projection $mi of quantum number $l not in valid set $ar"))
            newm = (newm..., convert(eltype(ar), mm))
        end

        new{O,typeof(newm)}(orb, newm)
    end
end
SpinOrbital(orb, m...) = SpinOrbital(orb, m)

Base.:(==)(a::SpinOrbital, b::SpinOrbital) =
    a.orb == b.orb && a.m == b.m

Base.hash(o::SpinOrbital, h::UInt) = hash(o.orb, hash(o.m, h))

function Base.show(io::IO, so::SpinOrbital)
    show(io, so.orb)
    projections = map(((l,m),) -> "m_$l = $m", zip(angular_momentum_labels(so.orb),so.m))
    write(io, "(", join(projections, ", "), ")")
end

function Base.show(io::IO, so::SpinOrbital{<:Orbital})
    show(io, so.orb)
    mℓ,ms = so.m
    write(io, to_subscript(mℓ))
    write(io, ms == half(1) ? "α" : "β")
end

function Base.show(io::IO, so::SpinOrbital{<:RelativisticOrbital})
    show(io, so.orb)
    write(io, "($(so.m[1]))")
end

degeneracy(::SpinOrbital) = 1

# We cannot order spin-orbitals of differing orbital types
Base.isless(a::SpinOrbital{O,M}, b::SpinOrbital{O,N}) where {O<:AbstractOrbital,M,N} = false

function Base.isless(a::SpinOrbital{<:O,M}, b::SpinOrbital{<:O,M}) where {O,M}
    a.orb < b.orb && return true
    a.orb > b.orb && return false

    for (ma,mb) in zip(a.m,b.m)
        ma < mb && return true
        ma > mb && return false
    end
    # All projections were equal
    return false
end

function Base.isless(a::SpinOrbital{<:Orbital}, b::SpinOrbital{<:Orbital})
    a.orb < b.orb && return true
    a.orb > b.orb && return false

    a.m[1] < b.m[1] ||
        a.m[1] == b.m[1] && a.m[2] > b.m[2] # We prefer α to appear before β
end

parity(so::SpinOrbital) = parity(so.orb)
symmetry(so::SpinOrbital) = (symmetry(so.orb), so.m...)

isbound(so::SpinOrbital) = isbound(so.orb)

Base.promote_type(::Type{SO}, ::Type{SO}) where {SO<:SpinOrbital} = SO

Base.promote_type(::Type{SpinOrbital{O}}, ::Type{SpinOrbital}) where O = SpinOrbital
Base.promote_type(::Type{SpinOrbital}, ::Type{SpinOrbital{O}}) where O = SpinOrbital

Base.promote_type(::Type{SpinOrbital{A,M}}, ::Type{SpinOrbital{B,M}}) where {A,B,M} =
    SpinOrbital{<:promote_type(A,B),M}

# ** Saving/loading

Base.write(io::IO, o::SpinOrbital{<:Orbital}) = write(io, 'n', o.orb, o.m[1], o.m[2].twice)
Base.write(io::IO, o::SpinOrbital{<:RelativisticOrbital}) = write(io, 'r', o.orb, o.m[1].twice)

function Base.read(io::IO, ::Type{SpinOrbital})
    kind = read(io, Char)
    if kind == 'n'
        orb = read(io, Orbital)
        mℓ = read(io, Int)
        ms = half(read(io, Int))
        SpinOrbital(orb, (mℓ, ms))
    elseif kind == 'r'
        orb = read(io, RelativisticOrbital)
        mj = half(read(io, Int))
        SpinOrbital(orb, mj)
    else
        error("Unknown SpinOrbital type $(kind)")
    end
end

# * Orbital construction from strings

function Base.parse(::Type{O}, orb_str) where {OO<:AbstractOrbital,O<:SpinOrbital{OO}}
    m = match(r"^(.*)\((.*)\)$", orb_str)
    # For non-relativistic spin-orbitals, we also support specifying
    # the m_ℓ quantum number using Unicode subscripts and the spin
    # label using α/β
    m2 = match(r"^(.*?)([₋]{0,1}[₁₂₃₄₅₆₇₈₉₀]+)([αβ])$", orb_str)
    if !isnothing(m)
        o = parse(OO, m[1])
        SpinOrbital(o, (split(m[2], ",")...,))
    elseif !isnothing(m2) && OO <: Orbital
        o = parse(OO, m2[1])
        SpinOrbital(o, from_subscript(m2[2]), m2[3])
    else
        throw(ArgumentError("Invalid spin-orbital string: $(orb_str)"))
    end
end

"""
    @so_str -> SpinOrbital{<:Orbital}

String macro to quickly construct a non-relativistic spin-orbital; it
is similar to [`@o_str`](@ref), with the added specification of
the magnetic quantum numbers ``m_ℓ`` and ``m_s``.

# Examples

```jldoctest
julia> so"1s(0,α)"
1s₀α

julia> so"kd(2,β)"
kd₂β

julia> so"2p(1,0.5)"
2p₁α

julia> so"2p(1,-1/2)"
2p₁β
```
"""
macro so_str(orb_str)
    parse(SpinOrbital{Orbital}, orb_str)
end

"""
    @rso_str -> SpinOrbital{<:RelativisticOrbital}

String macro to quickly construct a relativistic spin-orbital; it
is similar to [`@o_str`](@ref), with the added specification of
the magnetic quantum number ``m_j``.

# Examples

```jldoctest
julia> rso"2p-(1/2)"
2p-(1/2)

julia> rso"2p(-3/2)"
2p(-3/2)

julia> rso"3d(2.5)"
3d(5/2)
```
"""
macro rso_str(orb_str)
    parse(SpinOrbital{RelativisticOrbital}, orb_str)
end

"""
    spin_orbitals(orbital)

Generate all permissible spin-orbitals for a given `orbital`, e.g. 2p
-> 2p ⊗ mℓ = {-1,0,1} ⊗ ms = {α,β}

# Examples

```jldoctest
julia> spin_orbitals(o"2p")
6-element Array{SpinOrbital{Orbital{Int64},Tuple{Int64,HalfIntegers.Half{Int64}}},1}:
 2p₋₁α
 2p₋₁β
 2p₀α
 2p₀β
 2p₁α
 2p₁β

julia> spin_orbitals(ro"2p-")
2-element Array{SpinOrbital{RelativisticOrbital{Int64},Tuple{HalfIntegers.Half{Int64}}},1}:
 2p-(-1/2)
 2p-(1/2)

julia> spin_orbitals(ro"2p")
4-element Array{SpinOrbital{RelativisticOrbital{Int64},Tuple{HalfIntegers.Half{Int64}}},1}:
 2p(-3/2)
 2p(-1/2)
 2p(1/2)
 2p(3/2)
```

"""
function spin_orbitals(orb::O) where {O<:AbstractOrbital}
    map(reduce(vcat, Iterators.product(angular_momentum_ranges(orb)...))) do m
        SpinOrbital(orb, m...)
    end |> sort
end

"""
    @sos_str -> Vector{<:SpinOrbital{<:Orbital}}

Can be used to easily construct a list of [`SpinOrbital`](@ref)s.

# Examples

```jldoctest
julia> sos"3[s-p]"
8-element Array{SpinOrbital{Orbital{Int64},Tuple{Int64,HalfIntegers.Half{Int64}}},1}:
 3s₀α
 3s₀β
 3p₋₁α
 3p₋₁β
 3p₀α
 3p₀β
 3p₁α
 3p₁β
```
"""
macro sos_str(orbs_str)
    reduce(vcat, map(spin_orbitals, orbitals_from_string(Orbital, orbs_str)))
end

"""
    @rsos_str -> Vector{<:SpinOrbital{<:RelativisticOrbital}}

Can be used to easily construct a list of [`SpinOrbital`](@ref)s.

# Examples

```jldoctest
julia> rsos"3[s-p]"
8-element Array{SpinOrbital{RelativisticOrbital{Int64},Tuple{HalfIntegers.Half{Int64}}},1}:
 3s(-1/2)
 3s(1/2)
 3p-(-1/2)
 3p-(1/2)
 3p(-3/2)
 3p(-1/2)
 3p(1/2)
 3p(3/2)
```
"""
macro rsos_str(orbs_str)
    reduce(vcat, map(spin_orbitals, orbitals_from_string(RelativisticOrbital, orbs_str)))
end

export SpinOrbital, spin_orbitals, @so_str, @rso_str, @sos_str, @rsos_str


text:
const ISO4217 = Dict{Symbol, Tuple{Int, String, Int}}(
    :AED => (2,"UAE Dirham",784),
    :AFN => (2,"Afghani",971),
    :ALL => (2,"Lek",8),
    :AMD => (2,"Armenian Dram",51),
    :ANG => (2,"Netherlands Antillean Guilder",532),
    :AOA => (2,"Kwanza",973),
    :ARS => (2,"Argentine Peso",32),
    :AUD => (2,"Australian Dollar",36),
    :AWG => (2,"Aruban Florin",533),
    :AZN => (2,"Azerbaijanian Manat",944),
    :BAM => (2,"Convertible Mark",977),
    :BBD => (2,"Barbados Dollar",52),
    :BDT => (2,"Taka",50),
    :BGN => (2,"Bulgarian Lev",975),
    :BHD => (3,"Bahraini Dinar",48),
    :BIF => (0,"Burundi Franc",108),
    :BMD => (2,"Bermudian Dollar",60),
    :BND => (2,"Brunei Dollar",96),
    :BOB => (2,"Boliviano",68),
    :BOV => (2,"Mvdol",984),
    :BRL => (2,"Brazilian Real",986),
    :BSD => (2,"Bahamian Dollar",44),
    :BTN => (2,"Ngultrum",64),
    :BWP => (2,"Pula",72),
    :BYR => (0,"Belarussian Ruble",974),
    :BZD => (2,"Belize Dollar",84),
    :CAD => (2,"Canadian Dollar",124),
    :CDF => (2,"Congolese Franc",976),
    :CHE => (2,"WIR Euro",947),
    :CHF => (2,"Swiss Franc",756),
    :CHW => (2,"WIR Franc",948),
    :CLF => (4,"Unidad de Fomento",990),
    :CLP => (0,"Chilean Peso",152),
    :CNY => (2,"Yuan Renminbi",156),
    :COP => (2,"Colombian Peso",170),
    :COU => (2,"Unidad de Valor Real",970),
    :CRC => (2,"Costa Rican Colon",188),
    :CUC => (2,"Peso Convertible",931),
    :CUP => (2,"Cuban Peso",192),
    :CVE => (2,"Cabo Verde Escudo",132),
    :CZK => (2,"Czech Koruna",203),
    :DJF => (0,"Djibouti Franc",262),
    :DKK => (2,"Danish Krone",208),
    :DOP => (2,"Dominican Peso",214),
    :DZD => (2,"Algerian Dinar",12),
    :EGP => (2,"Egyptian Pound",818),
    :ERN => (2,"Nakfa",232),
    :ETB => (2,"Ethiopian Birr",230),
    :EUR => (2,"Euro",978),
    :FJD => (2,"Fiji Dollar",242),
    :FKP => (2,"Falkland Islands Pound",238),
    :GBP => (2,"Pound Sterling",826),
    :GEL => (2,"Lari",981),
    :GHS => (2,"Ghana Cedi",936),
    :GIP => (2,"Gibraltar Pound",292),
    :GMD => (2,"Dalasi",270),
    :GNF => (0,"Guinea Franc",324),
    :GTQ => (2,"Quetzal",320),
    :GYD => (2,"Guyana Dollar",328),
    :HKD => (2,"Hong Kong Dollar",344),
    :HNL => (2,"Lempira",340),
    :HRK => (2,"Kuna",191),
    :HTG => (2,"Gourde",332),
    :HUF => (2,"Forint",348),
    :IDR => (2,"Rupiah",360),
    :ILS => (2,"New Israeli Sheqel",376),
    :INR => (2,"Indian Rupee",356),
    :IQD => (3,"Iraqi Dinar",368),
    :IRR => (2,"Iranian Rial",364),
    :ISK => (0,"Iceland Krona",352),
    :JMD => (2,"Jamaican Dollar",388),
    :JOD => (3,"Jordanian Dinar",400),
    :JPY => (0,"Yen",392),
    :KES => (2,"Kenyan Shilling",404),
    :KGS => (2,"Som",417),
    :KHR => (2,"Riel",116),
    :KMF => (0,"Comoro Franc",174),
    :KPW => (2,"North Korean Won",408),
    :KRW => (0,"Won",410),
    :KWD => (3,"Kuwaiti Dinar",414),
    :KYD => (2,"Cayman Islands Dollar",136),
    :KZT => (2,"Tenge",398),
    :LAK => (2,"Kip",418),
    :LBP => (2,"Lebanese Pound",422),
    :LKR => (2,"Sri Lanka Rupee",144),
    :LRD => (2,"Liberian Dollar",430),
    :LSL => (2,"Loti",426),
    :LYD => (3,"Libyan Dinar",434),
    :MAD => (2,"Moroccan Dirham",504),
    :MDL => (2,"Moldovan Leu",498),
    :MGA => (2,"Malagasy Ariary",969),
    :MKD => (2,"Denar",807),
    :MMK => (2,"Kyat",104),
    :MNT => (2,"Tugrik",496),
    :MOP => (2,"Pataca",446),
    :MRO => (2,"Ouguiya",478),
    :MUR => (2,"Mauritius Rupee",480),
    :MVR => (2,"Rufiyaa",462),
    :MWK => (2,"Kwacha",454),
    :MXN => (2,"Mexican Peso",484),
    :MXV => (2,"Mexican Unidad de Inversion (UDI)",979),
    :MYR => (2,"Malaysian Ringgit",458),
    :MZN => (2,"Mozambique Metical",943),
    :NAD => (2,"Namibia Dollar",516),
    :NGN => (2,"Naira",566),
    :NIO => (2,"Cordoba Oro",558),
    :NOK => (2,"Norwegian Krone",578),
    :NPR => (2,"Nepalese Rupee",524),
    :NZD => (2,"New Zealand Dollar",554),
    :OMR => (3,"Rial Omani",512),
    :PAB => (2,"Balboa",590),
    :PEN => (2,"Nuevo Sol",604),
    :PGK => (2,"Kina",598),
    :PHP => (2,"Philippine Peso",608),
    :PKR => (2,"Pakistan Rupee",586),
    :PLN => (2,"Zloty",985),
    :PYG => (0,"Guarani",600),
    :QAR => (2,"Qatari Rial",634),
    :RON => (2,"Romanian Leu",946),
    :RSD => (2,"Serbian Dinar",941),
    :RUB => (2,"Russian Ruble",643),
    :RWF => (0,"Rwanda Franc",646),
    :SAR => (2,"Saudi Riyal",682),
    :SBD => (2,"Solomon Islands Dollar",90),
    :SCR => (2,"Seychelles Rupee",690),
    :SDG => (2,"Sudanese Pound",938),
    :SEK => (2,"Swedish Krona",752),
    :SGD => (2,"Singapore Dollar",702),
    :SHP => (2,"Saint Helena Pound",654),
    :SLL => (2,"Leone",694),
    :SOS => (2,"Somali Shilling",706),
    :SRD => (2,"Surinam Dollar",968),
    :SSP => (2,"South Sudanese Pound",728),
    :STD => (2,"Dobra",678),
    :SVC => (2,"El Salvador Colon",222),
    :SYP => (2,"Syrian Pound",760),
    :SZL => (2,"Lilangeni",748),
    :THB => (2,"Baht",764),
    :TJS => (2,"Somoni",972),
    :TMT => (2,"Turkmenistan New Manat",934),
    :TND => (3,"Tunisian Dinar",788),
    :TOP => (2,"Pa’anga",776),
    :TRY => (2,"Turkish Lira",949),
    :TTD => (2,"Trinidad and Tobago Dollar",780),
    :TWD => (2,"New Taiwan Dollar",901),
    :TZS => (2,"Tanzanian Shilling",834),
    :UAH => (2,"Hryvnia",980),
    :UGX => (0,"Uganda Shilling",800),
    :USD => (2,"US Dollar",840),
    :USN => (2,"US Dollar (Next day)",997),
    :UYI => (0,"Uruguay Peso en Unidades Indexadas (URUIURUI)",940),
    :UYU => (2,"Peso Uruguayo",858),
    :UZS => (2,"Uzbekistan Sum",860),
    :VEF => (2,"Bolivar",937),
    :VND => (0,"Dong",704),
    :VUV => (0,"Vatu",548),
    :WST => (2,"Tala",882),
    :XAF => (0,"CFA Franc BEAC",950),
    :XAG => (-1,"Silver",961),
    :XAU => (-1,"Gold",959),
    :XBA => (-1,"Bond Markets Unit European Composite Unit (EURCO)",955),
    :XBB => (-1,"Bond Markets Unit European Monetary Unit (E.M.U.-6)",956),
    :XBC => (-1,"Bond Markets Unit European Unit of Account 9 (E.U.A.-9)",957),
    :XBD => (-1,"Bond Markets Unit European Unit of Account 17 (E.U.A.-17)",958),
    :XCD => (2,"East Caribbean Dollar",951),
    :XDR => (-1,"SDR (Special Drawing Right)",960),
    :XOF => (0,"CFA Franc BCEAO",952),
    :XPD => (-1,"Palladium",964),
    :XPF => (0,"CFP Franc",953),
    :XPT => (-1,"Platinum",962),
    :XSU => (-1,"Sucre",994),
    :XTS => (-1,"TESTING",963),
    :XUA => (-1,"ADB Unit of Account",965),
    :XXX => (-1,"NO CURRENCY",999),
    :YER => (2,"Yemeni Rial",886),
    :ZAR => (2,"Rand",710),
    :ZMW => (2,"Zambian Kwacha",967),
    :ZWL => (2,"Zimbabwe Dollar",932))


text:
module day10

using DelimitedFiles, UnicodePlots, REPL

readInput() = permutedims(reshape([line[x] == '#' for line in readdlm(joinpath(@__DIR__, "input.txt")) for x in 1:length(line)], 24, 24), (2,1))

function part1(data)
    maxCount = 0
    for i = 1:size(data, 1), j= 1:size(data, 2)
        if data[i,j]
            maxCount = max(maxCount, length(countSight(data, i, j)))
        end
    end
    maxCount
end

function countSight(data, x, y)
    blocked = Set{Tuple{Int, Int}}()
    countAsteroid = 0
    counted = Vector{Tuple{Int, Int}}()
    for rayon = 1:max(size(data)...) - 1
        dx = -rayon
        for dy = -rayon:rayon
            if checkAsteroid(data, blocked, x, y, dx, dy)
                blockView!(data, blocked, x, y, dx, dy)
                countAsteroid += 1
                push!(counted, (x+dx, y+dy))
            end
        end
        dx = rayon
        for dy = -rayon:rayon
            if checkAsteroid(data, blocked, x, y, dx, dy)
                blockView!(data, blocked, x, y, dx, dy)
                countAsteroid += 1
                push!(counted, (x+dx, y+dy))
            end
        end
        dy = -rayon
        for dx = (-rayon+1):(rayon-1)
            if checkAsteroid(data, blocked, x, y, dx, dy)
                blockView!(data, blocked, x, y, dx, dy)
                countAsteroid += 1
                push!(counted, (x+dx, y+dy))
            end
        end
        dy = rayon
        for dx = (-rayon+1):(rayon-1)
            if checkAsteroid(data, blocked, x, y, dx, dy)
                blockView!(data, blocked, x, y, dx, dy)
                countAsteroid += 1
                push!(counted, (x+dx, y+dy))
            end
        end
    end

    # p = scatterplot([y], [size(data, 1) - x + 1], width = 2 * size(data, 1), height = 2 * size(data, 2), xlim=(1, size(data, 1)), ylim=(1, size(data, 2)), color = :blue, canvas = BlockCanvas)
    # p = scatterplot!(p, last.(counted), size(data, 1) .- first.(counted) .+ 1, color=:green)
    # p = scatterplot!(p, last.(blocked), size(data, 1) .- first.(blocked) .+ 1, color=:red)
    # display(p)

    counted
end

@inline checkAsteroid(data, blocked, x, y, dx, dy) = !((x+dx, y+dy) in blocked) && (1 <= x + dx <= size(data, 1)) && (1 <= y + dy <= size(data, 2)) && data[x+dx, y+dy]

function blockView!(data, blocked, x, y, dx, dy)
    xx = x + dx
    yy = y + dy
    step = gcd(abs(dx), abs(dy))
    dx, dy = dx ÷ step, dy ÷ step
    xx += dx
    yy += dy
    while 1 <= xx <= size(data, 1) && 1 <= yy <= size(data, 2)
        push!(blocked, (xx, yy))
        xx += dx
        yy += dy
    end
end

function part2(data, showplot = false, target = 200)
    data = copy(data)
    maxCount, maxI, maxJ = 0, 0, 0
    for i = 1:size(data, 1), j= 1:size(data, 2)
        if data[i,j]
            cpt = length(countSight(data, i, j))
            if cpt > maxCount
                maxCount, maxI, maxJ = cpt, i, j
            end
        end
    end

    vaporized = 0
    inSight = countSight(data, maxI, maxJ)
    while vaporized + length(inSight) < target
        for (i, j) in inSight
            data[i, j] = false
        end
        vaporized += length(inSight)
        inSight = countSight(data, maxI, maxJ)
    end

    sortRad(x) = begin
        d = atand((x .- (maxI, maxJ))...) .+ (3 * 90)
        d > 180 ? d - 360 : d
    end
    sort!(inSight, by = sortRad)

    if showplot
        p = scatterplot([maxJ], [size(data, 1) - maxI + 1], width = 2 * size(data, 1), height = 2 * size(data, 2),
            xlim=(1, size(data, 1)), ylim=(1, size(data, 2)), color = :blue, canvas = BrailleCanvas)
        display(p)
        terminal = REPL.Terminals.TTYTerminal("", stdin, stdout, stderr)
        for asteroid in inSight
            sleep(0.1)
            REPL.Terminals.clear(terminal)
            p = scatterplot!(p, [last(asteroid)], [size(data, 1) - first(asteroid) + 1], color=:green)
            display(p)
        end
    end

    x, y = inSight[target - 1] .- 1
    100y+x
end

end


text:
function AssetPricingTest3dPs(Re1,Re2,f,h,bandwidth=0,q=0,CrossRegType="beta",vvR=0)
#AssetPricingTest3dPs   Testing difference of two linear factor models by testing the
#                       difference of cross-sectional pricing errors. Allows Re1 ~= Re2
#
#
#
#
#  Usage: fnOutput = AssetPricingTest3dPs(Re1,Re2,f,h,[,bandwidth[,q[,CrossRegType[,vvR]]]])
#
#  Input:   Re1            Txn1, test assets
#           Re2            Txn2, test assets     (in many applications, Re2 = Re1)
#           f              TxK, 1st set of factors, not necessarily excess returns
#           h              TxL, 2nd set of factors
#           bandwidth      (optional) scalar of Newey-West band width,  [0]
#           q              (optional) vector, expected values of R*theta: [alpha delta]
#           CrossRegType   (optional) ('beta' or 'GLS') or {Kxn1;Lxn2} cell arrray, ['beta']
#           vvR            (optional) heterogeneous array of vectors, indices of
#                          moment conditions to test with R*g-q
#
#
#  Output:  fnOutput      8x1 heterogeneous array with the following contents:
#             [1] Joint         3x1, in each row WaldStat for alpha, delta, alpha-delta (if n1=n2)
#             [2] Indiv         3xmax(n1,n2), t-stat of alfa,delta and alpha-delta (if n1=n2)
#             [3] bf
#             [4] bh
#             [5] EReFit        1x2 cell array, each cell is a nx2 matrix with [actual,fitted ERe]
#             [6] theta         [n1(1+K)+n2(1+L)]x1, estimated parameters, [alpha;vec(beta);delta;vec(gamma)]
#                              Notice the order of parameters: [for const; for factor 1;...]
#             [7] V             matrix, Cov[sqrt(T)*theta]
#
#
#
#  Uses:    excise, HDirProdPs, NewEst3Ps
#
#
#
#
#  Paul.Soderlind@unisg.ch, June 2015, to Julia Dec 2015
#------------------------------------------------------------------------------

  n1 = size(Re1,2)
  n2 = size(Re2,2)
  if q == 0                                  #scalar 0
    q = Any[zeros(n1,1),zeros(n2,1)]         #default values for options
  end
  if vvR == 0
    vvR = Any[1:n1,1:n2]
  end
  #----------------------------------------------------------

  K  = size(f,2)
  L  = size(h,2)

  (Re1,Re2,f,h,) = excise4mPs(Re1,Re2,f,h)
  T              = size(Re1,1)                      #no. obs

  Ef   = mean(f,1)'
  Eh   = mean(h,1)'
  ERe1 = mean(Re1,1)'
  ERe2 = mean(Re2,1)'

  f1 = [ones(T,1) f]
  K1 = 1 + K
  h1 = [ones(T,1) h]
  L1 = 1 + L

  bf         = f1\Re1
  epsf       = Re1 - f1*bf
  bf         = bf'                       #n1xK1
  alfaM      = bf[:,1]
  betaM      = bf[:,2:end]               #n1xK, Re1(t) = alfa * beta*f(t) + e(t)
  gf         = HDirProdPs(f1,epsf)       #moment conditions

  bh         = h1\Re2
  epsh       = Re2 - h1*bh
  bh         = bh'                       #n2xL1
  deltaM     = bh[:,1]
  gammaM     = bh[:,2:end]               #n2xL
  gh         = HDirProdPs(h1,epsh)
  #----------------------------------------------------------

  if CrossRegType == "beta"
    B = betaM'                               #traditional cross-sectional regression
    G = gammaM'
  elseif CrossRegType == "GLS"
    B = betaM'* inv(cov(epsf))               #GLS
    G = gammaM'*inv(cov(epsh))
  else
    B = CrossRegType[1]                      #predefined matrix
    G = CrossRegType[2]
  end
                                             #cross-sectional estimate of price of factor risk
  lambda  = (B*betaM)\(B*ERe1)               #solve (B*betaM)*lambda = B*ERe
  psi     = (G*gammaM)\(G*ERe2)
  theta   = [vec(bf);lambda;vec(bh);psi]     #parameters (combined)

  EReHatf = betaM*lambda                     #nxK1 * K1x1
  EReHath = gammaM*psi                       #nxL1 * L1x1
  EReFit  = Any[[ERe1 EReHatf],[ERe2 EReHath]]  #actual and fitted ERe

  gfb     = Re1 .- EReHatf'                  #moment condition ERe - beta*lambda
  ghb     = Re2 .- EReHath'                  #moment condition ERe - gamma*psi

  g       = [gf gfb gh ghb]                  #moment conditions (combined)
  S       = NewEst3Ps(g,bandwidth)           #Cov[sqrt(T)*gbar]
  gf   = nothing                               #can save some memory
  gfb  = nothing
  gh   = nothing
  ghb  = nothing

  gbar = mean(g,1)'                          #moment conditions, sample average
  gNum = size(g,2)                           #no. moment conditions

  Sff  = f1'f1/T                             #Jacobian
  Df11 = kron(Sff,eye(n1))                   #upper left corner
  Shh  = h1'h1/T
  Dh11 = kron(Shh,eye(n2))
  Df12 = kron([0.0 vec(lambda)'],eye(n1))    #lower left
  Dh12 = kron([0.0 vec(psi)'],   eye(n2))
  Df   = - [ Df11  zeros(n1*(1+K),K);
             Df12  betaM           ]
  Dh   = - [ Dh11  zeros(n2*(1+L),L);
             Dh12  gammaM          ]
  D    = cat([1,2],Df,Dh)                    #block diagonal cat

  A    = cat([1,2],eye(n1*K1),B,eye(n2*L1),G)

  AD_1 = inv(A*D)
  V    = AD_1*A*S*A'AD_1'                    #Cov[sqrt(T)*theta]
  Psia = eye(gNum) - D*AD_1*A
  Psi  = Psia*S*Psia'                        #Cov[sqrt(T) gbar], rank gNum-k
  #----------------------------------------------------------

  RM1 = [zeros(n1,n1*K1) eye(n1)      zeros(n1,n2*L1) zeros(n1,n2)]    #Re1 - beta*lambda
  RM2 = [zeros(n2,n1*K1) zeros(n2,n1) zeros(n2,n2*L1) eye(n2)     ]    #Re2 - gamma*psi
  RM  = Any[RM1[vvR[1],:],RM2[vvR[2],:]]                               #pick out rows

  if size(RM[1],1) == size(RM[2],1)
    RM = Any[RM[1],RM[2],RM[1] - RM[2]]
    q  = Any[ q[1], q[2], q[1] -  q[2]]
  end

  Indiv = fill(NaN,(length(RM),max(n1,n2)))     #t-stat
  Joint = fill(NaN,length(RM))                  #WaldStat
  for i = 1:length(RM)                          #loop over (alpha,delta,alpha-delta)
    R          = RM[i]
    q_i        = q[i]
    Rstat      = R*gbar - vec(q_i)              #eg. gbar[2:3]
    Lambda     = R*Psi*R'                       #eg. Cov[sqrt(T)*gbar[2:3]
    stat_std   = sqrt(diag(Lambda/T))
    tstat      = Rstat./stat_std
    Joint_i    =  T * Rstat'pinv(Lambda)*Rstat  #reduced rank
    Joint[i]   = Joint_i[1]                     #[1] to make it a scalar
    Indiv[i,1:length(tstat)] = tstat'
  end
  #----------------------------------------------------------

              #    1     2   3   4   5      6   7
  fnOutput = Any[Joint,Indiv,bf,bh,EReFit,theta,V]

  return fnOutput

end
#------------------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
import PredictMD
import PredictMD.Cleaning

Test.@test(
    PredictMD.Cleaning.x_contains_y("abc", ["xyz", "abc", "123",])
    )

Test.@test(
    !PredictMD.Cleaning.x_contains_y("abc", ["xyz", "opqrst", "123",])
    )

Test.@test(
    PredictMD.Cleaning.symbol_begins_with(:abcdefg, "abc")
    )

Test.@test(
    !PredictMD.Cleaning.symbol_begins_with(:abcdefg, "xyz")
    )


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>2.0, :alg=>BOX(δ=0.01), :disjointness_method=>BoxEnclosure(), :intersect_source_invariant=>false),
    "X0"   => Hyperrectangle(low=[0.5, 0.5, 0.0, 0.0], high=[0.55, 0.55, 0, 0]),
    "info" => Dict(
        "name" => "navigation_system",
        "dim" => 4,
        "linear" => true,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
#=
 Tutorial 2
=#
# Start repl Julia to open a command line in Julia

my_integer = 5
my_float = .5
## To run section of the code in VSCode
my_range = 1:7
my_vector = [1,3,5]
my_string = "hello"
my_character = 'o'
my_boolean = true

typeof(my_integer)
typeof(my_float)

# Alt + enter to run the code line by line

# \alpha
α = 5;

@show α + my_integer

print("α = $α \n", "My string is: $(my_string)")

## Array can contain any type of element
A = [1, 2, 3]
A = [1 2 3]
ones(3)
zeros(3)
C = [1 3; 5 9]

A = 1:9
B = 1:0.5:9

# Julia starts indexing at 1
B[3]

typeof(C)
size(C)
length(C)
ndims(C)

C[end]
C[1,2]

@show C
@show C[1,2]
@show C[2, :]

#=
Create a new environment
] to active the package manager in the command line
active .
activate test_env

add LinearAlgebra to add packages
status to see what packages are installed

using LinearAlgebra to use the package

backspace to go back to the Julia Command Line
=#

using LinearAlgebra
C*I
C

sin.([1,2,3])
inv(C)*C
C\C # to inverse a matrix

using Random
rand(1)
randn(3)

Random.seed!(1234)

# If Else statements
a = 3
if a == 3
    print("a is 3")
end

c = true
if c
    print("True")
end

# Condition ? If true : Else
c ? print("c is true") : nothing

if c
    print("c is true")
else
    print("c is not true")
end

# Logical operators
b = true
if c | b # or, a & b and, && check if the first term is true then goes to the second term
    print("a or b are true")
end

# For loops
vals = ["a", "b", "c"]
for val in vals
    println(val)
end

e = 3
for val in vals
    e = e + 1
end

for i in eachindex(vals)
    println(vals[i])
end

for (index, value) in enumerate(vals)
    print("Element $index of vals is $value\n")
end

# Functions

function cobbdouglas(l, k, α, β = 1-α; A = 1) # To have default value A = 1
    y = A*l^α*k^β
    return y
end

cobbdouglas(1, 3, 0.5)
cobbdouglas(1, 3, 0.5, 0.6)

# To store the returned output of the function
@show output = cobbdouglas(1, 3, 0.5)

# To evaluate the function at the two points
cobbdouglas.([1,2], 3, 0.5)

cobbdouglasshort(l,k,α) = l^α*k^(1-α)
cobbdouglasshort(1,3,0.5)


text:
# Demo using DiscreteHMM and Baum-Welch algorithm
# to estimate a simple model of natural language.
include("models.jl")

module demo_discrete
using DiscreteHMM
can_save = (Pkg.installed("JLD")!=nothing)  # check if JLD is installed
if can_save; using JLD; end

# Settings
m = 100  # number of hidden states to use in estimated model
n_max = 100000  # maximum number of words/symbols to use
from_file = false  # load estimated parameters from file
run_profiler = false  # run profiler to assess performance
tolerance = 1e-2  # convergence tolerance for Baum-Welch
n_gen = 1000  # length of sequence to generate using estimated parameters
filename = "simple.txt"  # file containing a sample of simple natural language

# Read and preprocess data
function preprocess(filename,n_max)
	# Read file
	f = open(filename)
	s = readall(f)
	close(f)
	# Remove whitespace and remove some punctuation
	s = replace(s,r"\s+"," ")
	s = replace(s,r"[\"]","")
	s = replace(s,r"' | '"," ")
	s = replace(s,r"(\W) | (\W)",s" \1 ")
	s = strip(s)
	words = split(s,' ')
	# Consider only the first n_max words/symbols
	words = words[1:min(n_max,length(words))]
	# Extract vocabulary
	vocabulary = sort(unique(words))
	V = length(vocabulary)
	code = Dict(zip(vocabulary,1:V))  # mapping from words to indices
	return words,vocabulary,code
end

# Produce the numerical sequence corresponding to a sequence of words.
function encode(words,code)
	return Int[code[w] for w in words]
end

# Produce the sequence of text corresponding to a numerical sequence.
function compose(x,vocabulary)
	s = join(vocabulary[x]," ")
	s = replace(s,r" (\W) ",s"\1 ")
	return s
end

# Preprocess and encode the text in the given file
words,vocabulary,code = preprocess(filename,n_max)
V = length(vocabulary)
x = encode(words,code)
n = length(x)

# Save encoding to file
#writedlm("x.txt",x)
#writedlm("code.txt",[(1:V) vocabulary],' ')

# Sanity check of forward algorithm implementation
if true
	log_pi = log(ones(m)/m)
	log_T = log(ones(m,m)/m)
	phi = [ones(V)/V for r = 1:m]
	G,log_m = DiscreteHMM.forward(x,log_pi,log_T,phi)
	log_m_exact = -n*log(V)
	@assert(abs(log_m - log_m_exact) < 1e-10, "Log marginal does not match")
end

# Estimate parameters using Baum-Welch
if from_file
	(log_pi,log_T,phi,log_m) = load("estimated_params.jld","params")
else
    if run_profiler
        @profile log_pi,log_T,phi,log_m = DiscreteHMM.estimate(x,m,tolerance)
        Profile.print() #format = :flat)
        Profile.clear()
    else
		log_pi,log_T,phi,log_m = DiscreteHMM.estimate(x,m,tolerance)
	end
	if can_save
		save("estimated_params.jld","params",(log_pi,log_T,phi,log_m))
	end
end

# Generate a new sequence from the estimated model
println("Generating a sample of text using the estimated model:")
y,zy = DiscreteHMM.generate(n_gen,log_pi,log_T,phi)
println(compose(y,vocabulary))

end # module

max_stars_repo_name: jwmi/HMM

max_issues_repo_name: jwmi/HMM

max_forks_repo_name: jwmi/HMM


text:
"""
    readphase(filename; rescale=true, keyargs...)

Reads the NIfTI phase with sanity checking and optional rescaling to [-π;π].

# Examples
```julia-repl
julia> phase = readphase("phase.nii")
```

### Optional keyargs are forwarded to `niread`:
$(@doc niread)
"""
function readphase(filename; rescale=true, keyargs...)
    phase = niread(filename; keyargs...)
    if phase.header.scl_slope == 0 # slope of 0 is always wrong
        phase.header.scl_slope = 1
    end
    if rescale
        minp, maxp = Float32.(approxextrema(phase))
        if isapprox(maxp - minp, 2π; atol=0.1) # no rescaling required
            return phase
        end
        minp, maxp = Float32.(approxextrema(phase.raw))
        if isapprox(maxp - minp, 2π; atol=0.1) # no rescaling required, but header wrong
            phase.header.scl_slope = 1
            phase.header.scl_inter = 0
        else # rescaling
            phase.header.scl_slope = 2pi / (maxp - minp)
            phase.header.scl_inter = -pi - minp * phase.header.scl_slope
        end
    end
    return phase
end

"""
    readmag(filename; rescale=false, keyargs...)

Reads the NIfTI magnitude with sanity checking and optional rescaling to [0;1].

# Examples
```julia-repl
julia> mag = readmag("mag.nii")
```

### Optional keyargs are forwarded to `niread`:
$(@doc niread)
"""
function readmag(fn; rescale=false, keyargs...)
    mag = niread(fn; keyargs...)
    if mag.header.scl_slope == 0
        mag.header.scl_slope = 1
    end
    if rescale
        mini, maxi = Float32.(approxextrema(mag.raw))
        mag.header.scl_slope = 1 / (maxi - mini)
        mag.header.scl_inter = - mini * mag.header.scl_slope
    end
    return mag
end

Base.copy(x::NIfTI.NIfTI1Header) = NIfTI.NIfTI1Header([getfield(x, k) for k ∈ fieldnames(NIfTI.NIfTI1Header)]...)

function Base.similar(header::NIfTI.NIfTI1Header)
    hdr = copy(header)
    hdr.scl_inter = 0
    hdr.scl_slope = 1
    return hdr
end

"""
    header(v::NIVolume)

Returns a copy of the header with the orientation information.

# Examples
```julia-repl
julia> vol = readmag("image.nii")
julia> hdr = header(vol)
julia> savenii(vol .+ 10, "vol10.nii"; header=hdr)
```
"""
header(v::NIVolume) = similar(v.header)

function savenii(image, name, writedir, header=nothing)
    if isnothing(writedir) return end
    if splitext(name)[2] != ".nii"
        name = name * ".nii"
    end
    savenii(image, joinpath(writedir, name); header)
end
"""
    savenii(image::AbstractArray, filepath; header=nothing)

    savenii(image::AbstractArray, name, writedir, header=nothing)

Warning: MRIcro can only open images with types Int32, Int64, Float32, Float64

# Examples
```julia-repl
julia> savenii(ones(64,64,5), "image.nii")

julia> savenii(ones(64,64,5), "image2", "folder")
```
"""
function savenii(image::AbstractArray, filepath; header=nothing)
    vol = NIVolume([h for h in [header] if h !== nothing]..., image)
    niwrite(filepath, vol)
    return filepath
end

ConvertTypes = Union{BitArray, AbstractArray{UInt8}} #TODO debug NIfTI
MriResearchTools.savenii(image::ConvertTypes, args...;kwargs...) = savenii(Float32.(image), args...;kwargs...)

"""
    write_emptynii(size, path; datatype=Float32, header=NIVolume(zeros(datatype, 1)).header)

Writes an empty NIfTI image to disk that can be used for memory-mapped access.

# Examples
```julia-repl
julia> vol = write_emptynii((64,64,64), "empty.nii")
julia> vol.raw[:,:,1] .= ones(64,64) # synchronizes mmapped file on disk
```

Warning: MRIcro can only open images with types Int32, Int64, Float32, Float64
"""
function write_emptynii(sz, path; datatype=Float32, header=NIVolume(zeros(datatype, 1)).header)
    header = copy(header)
    header.dim = Int16.((length(sz), sz..., ones(8-1-length(sz))...))
    header.datatype = NIfTI.eltype_to_int16(datatype)
    header.bitpix = NIfTI.nibitpix(datatype)

    if isfile(path) rm(path) end
    open(path, "w") do file
        write(file, header)
        write(file, Int32(0)) # offset of 4 bytes
    end
    return niread(path; mmap=true, mode="r+")
end

mmtovoxel(sizemm, nii::NIVolume) = mmtovoxel(sizemm, nii.header)
mmtovoxel(sizemm, header::NIfTI.NIfTI1Header) = mmtovoxel(sizemm, header.pixdim)
mmtovoxel(sizemm, pixdim) = sizemm ./ pixdim

getcomplex(mag::NIVolume, phase::NIVolume) = getcomplex(mag.raw, phase.raw)

function Base.setindex!(vol::NIVolume{<:AbstractFloat}, v, i...)
    scaled = v / vol.header.scl_slope + vol.header.scl_inter
    setindex!(vol.raw, scaled, i...)
end


text:
"""
Low pass filter (LPF) -like FDI.
"""
struct LPFFDI <: AbstractFDI
    τ
    function LPFFDI(τ)
        @assert τ > 0
        new(τ)
    end
end

function State(fdi::LPFFDI)
    return function (dim_input::Int)
        Λ̂ = Diagonal(ones(dim_input))
    end
end

"""
Λ: effectiveness matrix
"""
function Dynamics!(fdi::LPFFDI)
    @unpack τ = fdi
    return function (dΛ̂, Λ̂, p, t; Λ)
        dΛ̂ .= (Λ - Λ̂) / τ
    end
end


text:
mutable struct ExecFlag{T}
    symbol      :: Symbol
    name        :: String
    description :: String
    value       :: T
    minus_count :: Int
end
function ExecFlag(e::ExecFlag, value)
    return ExecFlag(e.symbol, e.name, e.description, value, e.minus_count)
end
function ExecFlag(p::Pair)
    return ExecFlag(first(p), String(first(p)), "", last(p), 1)
end
function ExecFlag(p::Pair, count::Int)
    return ExecFlag(first(p), String(first(p)), "", last(p), count)
end

StructTypes.StructType(::Type{<:ExecFlag}) = StructTypes.Struct()

Base.hash(e::ExecFlag, h::UInt)        = hash(e.symbol, hash(e.value, h))
Base.:(==)(e1::ExecFlag, e2::ExecFlag) = e1.symbol == e2.symbol && e1.value == e2.value
Base.eltype(::ExecFlag{T}) where {T}   = T

"""
    Exec(;exec::String = "", dir::String = "", flags::Vector{ExecFlag} = ExecFlag[])

Representation of an `executable` that will run the [`Calculation`](@ref Calculation).
Basically `dir/exec --<flags>` inside a job script.

    Exec(exec::String, dir::String, flags::Pair{Symbol}...)

Will first transform `flags` into a `Vector{ExecFlag}`, and construct the [`Exec`](@ref).
"""
@with_kw mutable struct Exec
    name::String = ""
    exec::String = ""
    dir::String = ""
    flags::Vector{ExecFlag} = ExecFlag[]
    modules::Vector{AbstractString} = String[]
    parallel::Bool = true
end

function Exec(name::String, exec::String, dir::String, flags::Pair...; kwargs...)
    _flags = ExecFlag[]
    ismpi = occursin("mpi", exec)
    for (f, v) in flags
        if ismpi
            mflag = mpi_flag(f)
            @assert mflag !== nothing "$f is not a recognized mpirun flag."
        end
        push!(_flags, ExecFlag(f => v))
    end
    return Exec(name=name,exec=exec, dir=dir, flags=_flags; kwargs...)
end

StructTypes.StructType(::Type{Exec}) = StructTypes.Mutable()

const EXEC_DIR = DFC.config_path("execs")

function load_exec(name)
    return JSON3.read(read(joinpath(EXEC_DIR, "$name.json"), String), Exec)
end

function load_execs()
    return JSON3.read.(read.(joinpath.((EXEC_DIR,), readdir(EXEC_DIR)), String), Exec)
end

save(e::Exec) =
    JSON3.write(joinpath(EXEC_DIR, "$(e.name).json"), e)

function isrunnable(e::Exec)
    # To find the path to then ldd on
    fullpath = !isempty(e.dir) ? joinpath(e.dir, e.exec) : Sys.which(e.exec)
    if fullpath !== nothing && ispath(fullpath)
        out = Pipe()
        err = Pipe()
        #by definition by submitting something with modules this means the server has them
        if !isempty(e.modules)
            cmd = `echo "source /etc/profile && module load $(join(e.modules, " ")) && ldd $fullpath"`
        else
            cmd = `echo "source /etc/profile && ldd $fullpath"`
        end
        run(pipeline(pipeline(cmd, stdout=ignorestatus(`bash`)),stdout = out, stderr=err))
        close(out.in)
        close(err.in)

        stderr = String(read(err))
        stdout = String(read(out))
        # This basically means that the executable would run
        return !occursin("not found", stdout) && !occursin("not found", stderr)
    else
        return false
    end
end

function known_exec_names()
    d = readdir(EXEC_DIR)
    return isempty(d) ? String[] : map(x -> splitext(x)[1], d)
end

"Verifies the validity of an executable and also registers it if it wasn't known before."
function verify_exec(e::Exec)
    known_execs = known_exec_names()
    if e.name ∈ known_execs && e == load_exec(e.name)
        return true
    else
        valid = isrunnable(e)
        if valid
            save(e)
        end
        return valid
    end
end

"Finds all executables that are known with the same exec name."
function known_execs(exec::AbstractString, dir="")
    return filter(x -> isempty(dir) ? x.exec == exec : x.exec == exec && x.exec == dir, load_execs())
end

const RUN_EXECS = ["mpirun", "mpiexec", "srun"]

hasflag(exec::Exec, s::Symbol) = findfirst(x -> x.symbol == s, exec.flags) != nothing

path(exec::Exec) = joinpath(exec.dir, exec.exec)

function Base.:(==)(e1::Exec, e2::Exec)
    if e1.exec != e2.exec || e1.dir != e2.dir
        return false
    else
        for f in e1.flags
            f2 = getfirst(isequal(f), e2.flags)
            f2 === nothing && return false
        end
        return all(x->x ∈ e1.modules, e2.modules)
    end
end

Base.hash(e::Exec, h::UInt) = hash(e.exec, hash(e.dir, hash(e.modules, hash(e.flags, h))))

function Base.string(e::Exec)
    direxec = joinpath(e.dir, e.exec)
    str = "$direxec"
    for flag in e.flags
        str *= " $(join(fill('-', flag.minus_count)))$(flag.symbol)"
        if !isa(flag.value, AbstractString)
            for v in flag.value
                str *= " $v"
            end
        else
            str *= " $(flag.value)"
        end
    end
    return str
end

Base.setindex!(e::Exec, v, s::Symbol) =
    set_flags!(e, s => v)

function Base.getindex(e::Exec, s::Symbol)
    out = getfirst(x -> x.symbol == s, e.flags)
    out === nothing && throw(KeyError(s))
    return out
end

function set_flags!(exec::Exec, flags...)
    for (f, val) in flags
        flag = isa(f, String) ? getfirst(x -> x.name == f, exec.flags) :
               getfirst(x -> x.symbol == f, exec.flags)
        if flag != nothing
            flag.value = convert(eltype(flag), val)
        else
            found = false
            #TODO generalize this
            for (f1, f2) in zip((is_qe_exec, is_wannier_exec, is_mpi_exec),
                    (qe_execflag, wan_execflag, mpi_flag))
                if f1(exec)
                    def_flag = f2(f)
                    if def_flag !== nothing
                        push!(exec.flags, ExecFlag(def_flag, val))
                        found = true
                        break
                    end
                end
            end
            if !found
                error("Flag $f was not found in allowed executable flags.")
            end
        end
    end
    return exec.flags
end

function rm_flags!(exec::Exec, flags...)
    for f in flags
        if isa(f, String)
            filter!(x -> x.name != f, exec.flags)
        elseif isa(f, Symbol)
            filter!(x -> x.symbol != f, exec.flags)
        end
    end
    return exec.flags
end

#### MPI Exec Functionality ###
include(joinpath(DFC.DEPS_DIR, "mpirunflags.jl"))
const MPI_FLAGS = _MPI_FLAGS()

mpi_flag(flag::AbstractString) = getfirst(x -> x.name == flag, MPI_FLAGS)
mpi_flag(flag::Symbol) = getfirst(x -> x.symbol == flag, MPI_FLAGS)

function mpi_flag_val(::Type{String}, line, i)
    v = line[i+1]
    return v, i + 2
end

function mpi_flag_val(::Type{Vector{String}}, line, i)
    tval = String[]
    while i + 1 <= length(line) && !occursin('-', line[i+1])
        push!(tval, line[i+1])
        i += 1
    end
    return tval, i + 1
end

function mpi_flag_val(::Type{Int}, line, i)
    if line[i+1][1] == '$'
        tval = line[i+1]
    else
        tval = parse(Int, line[i+1])
    end
    return tval, i + 2
end

function mpi_flag_val(::Type{Vector{Int}}, line, i)
    tval = Union{Int,String}[]
    while i + 1 <= length(line) && !occursin('-', line[i+1])
        if line[i+1][1] == '$'
            push!(tval, line[i+1])
        else
            push!(tval, parse(Int, line[i+1]))
        end
        i += 1
    end
    return tval, i + 1
end

function mpi_flag_val(::Type{Pair{Int,String}}, line, i)
    tval = Union{Int,String}[]
    while i + 1 <= length(line) && !occursin('-', line[i+1])
        if line[i+1][1] == '$'
            push!(tval, line[i+1])
        else
            tparse = tryparse(Int, line[i+1])
            if tparse != nothing
                push!(tval, tparse)
            else
                push!(tval, string(line[i+1]))
            end
        end
        i += 1
    end
    return tval, i + 1
end
function parse_mpi_flags(line::Vector{<:SubString})
    eflags = ExecFlag[]
    i = 1
    while i <= length(line)
        s = line[i]
        if s[1] != '-'
            break
        end
        if s[2] == '-'
            mflag = mpi_flag(strip(s, '-'))
        else
            mflag = mpi_flag(Symbol(strip(s, '-')))
        end
        @assert mflag !== nothing "$(strip(s, '-')) is not a recognized mpi_flag"
        val, i = mpi_flag_val(eltype(mflag), line, i)
        push!(eflags, ExecFlag(mflag, val))
    end
    return eflags
end

is_mpi_exec(exec::Exec) = occursin("mpi", exec.exec)

#Wannier
const WAN_EXECS = ["wannier90.x"]

const WAN_EXECFLAGS = ExecFlag[ExecFlag(:pp, "preprocess",
                                        "Whether or not to preprocess the wannier calculation",
                                        false, 1),]

wan_execflag(flag::AbstractString) = getfirst(x -> x.name == flag, WAN_EXECFLAGS)
wan_execflag(flag::Symbol) = getfirst(x -> x.symbol == flag, WAN_EXECFLAGS)
function parse_wan_execflags(line::Vector{<:AbstractString})
    flags = ExecFlag[]
    i = 1
    while i <= length(line)
        s = strip(line[i], '-')
        push!(flags, ExecFlag(wan_execflag(Symbol(s)), nothing))
        i += 1
    end
    return flags
end

is_wannier_exec(exec::Exec) = exec.exec ∈ WAN_EXECS

#QE
const QE_EXECS = ["pw.x", "projwfc.x", "pp.x", "ld1.x", "ph.x", "pw2wannier90.x", "hp.x"]

const QE_EXECFLAGS = ExecFlag[ExecFlag(:nk, "kpoint-pools",
                                       "groups k-point parallelization into nk processor pools",
                                       0, 1),
                              ExecFlag(:ntg, "task-groups", "FFT task groups", 0, 1),
                              ExecFlag(:ndiag, "diag",
                                       "Number of processes for linear algebra", 0, 1),
                              ExecFlag(:ni, "images",
                                       "Number of processes used for the images", 0, 1),
                               ExecFlag(:npool, "pools",
                                       "Number of processes used for the pools", 0, 1)]

qe_execflag(flag::AbstractString) = getfirst(x -> x.name == flag, QE_EXECFLAGS)
qe_execflag(flag::Symbol) = getfirst(x -> x.symbol == flag, QE_EXECFLAGS)

function parse_qe_execflags(line::Vector{<:AbstractString})
    flags = ExecFlag[]
    i = 1
    while i <= length(line)
        s = strip(line[i], '-')
        push!(flags, ExecFlag(qe_execflag(Symbol(s)), parse(Int, line[i+1])))
        i += 2
    end
    return flags
end

is_qe_exec(exec::Exec) = exec.exec ∈ QE_EXECS

# Elk
const ELK_EXECS = ["elk", "elk-omp"]

is_elk_exec(exec::Exec) = exec.exec ∈ ELK_EXECS

allexecs() = vcat(RUN_EXECS, QE_EXECS, WAN_EXECS, ELK_EXECS)
parseable_execs() = vcat(QE_EXECS, WAN_EXECS, ELK_EXECS)
has_parseable_exec(l::String) = occursin(">", l) && any(occursin.(parseable_execs(), (l,)))

isparseable(exec::Exec) = exec.exec ∈ parseable_execs()

function parse_generic_flags(flags::Vector{<:SubString})
    out = ExecFlag[]
    f = :none
    c = 1
    for s in flags
        if s[1] == '-'
            c = count(isequal('-'), s)
            f = Symbol(strip(s, '-'))
        else
            v = tryparse(Int, s)
            if v === nothing
                v = tryparse(Float64, s)
                if v === nothing
                    v = s
                end
            end
            push!(out, ExecFlag(f => v, c))
        end
    end
    return out
end


text:
"""
    mutable struct Container
        ID::String
        Image::String
        Command::Union{Missing, String} = missing
        CreatedAt::Union{Missing, String} = missing
        Names::Union{Missing, String} = missing
        Ports::Union{Missing, String} = missing
        exists::Bool = true
    end

Represents a docker container on the local environment. Should not be instantiated directly. If
`exists` is `true`, then the container is assumed to exit and so should be visible from utilities
such as `docker container ls --all`.
"""
@with_kw mutable struct Container
  ID::String
  Image::String
  Command::Union{Missing, String} = missing
  CreatedAt::Union{Missing, String} = missing
  Names::Union{Missing, String} = missing
  Ports::Union{Missing, String} = missing
  exists::Bool = true
end
Base.:(==)(a::Container, b::Container) = a.ID[1:docker_hash_limit] == b.ID[1:docker_hash_limit]

"""
    is_container_running(con::Container)::Bool

Returns `true` if the given docker container is currently running (not stopped).
"""
function is_container_running(con::Container)::Bool
  running_containers = get_all_containers()
  @debug running_containers
  @debug con
  con in running_containers
end

"""
    stop_container(con::Container)

Stops the given docker container, if currently running.
"""
function stop_container(con::Container)
  if is_container_running(con)
    run(`docker stop $(con.ID)`)
  end
end

"""
    delete!(con::Container)

Deletes the passed container from the local docker system. The `Container` instance continues to
exist, but has its `exists` attribute set to `false`.

"""
function delete!(con::Container)
  con.exists || error("Container does not exist")
  run(`docker container rm $(con.ID)`)
  con.exists = false
end

"""
    get_all_containers(args::Vector{String} = Vector{String}())::Vector{Container}

Returns a list of containers currently available on the local machine.

`args` are additional arguments passed to the `docker ps` call that this function wraps.
"""
function get_all_containers(args::Vector{String} = Vector{String}())::Vector{Container}
  docker_output = readchomp(`docker ps $args --format '{{json .}}'`)
  parse_docker_ls_output(Container, docker_output)
end

function get_all_containers(image::LocalImage; args::Vector{String}=Vector{String}())::Vector{Container}
  @debug image.ID
  get_all_containers([
                  ["--filter", "ancestor=$(image.ID[1:docker_hash_limit])"]
                  args
                 ])
end


text:
@testset "hyper_to_matrix_tests" begin
    # Trivial graphs
    @test to_adjacency_matrix(Hypergraphs(Array{Int64,1}[], 1, 0))[1,:] == [0]
    @test to_adjacency_matrix(Hypergraphs([[1]], 1, 1))[1,:] == [0]
    @test to_adjacency_matrix(Hypergraphs([[1,1]], 1, 1))[1,:] == [2]

    # Only self-loops
    @test to_adjacency_matrix(Hypergraphs([[1,1,1], [2,2], [3]], 3, 3)) == [6 0 0; 0 2 0; 0 0 0]

    # Multi-edges and self-loops
    @test to_adjacency_matrix(Hypergraphs([[1,4], [2,3], [1,2,3]], 4, 3)) == [0 1 1 1; 1 0 2 0; 1 2 0 0; 1 0 0 0]
    @test to_adjacency_matrix(Hypergraphs([[1,2], [1,3], [3,3]], 3, 3)) == [0 1 1; 1 0 0; 1 0 2]

    # Conversion to simple
    @test to_adjacency_matrix(Hypergraphs([[1,4], [2,3], [1,2,3]], 4, 3); simple=true) == [0 1 1 1; 1 0 1 0; 1 1 0 0; 1 0 0 0]
    @test to_adjacency_matrix(Hypergraphs([[1,2,3], [2,3,4], [1,4], [5,6]], 6, 4); simple=true) == [0 1 1 1 0 0; 1 0 1 1 0 0; 1 1 0 1 0 0; 1 1 1 0 0 0; 0 0 0 0 0 1; 0 0 0 0 1 0]
end

@testset "matrix_to_hyper_tests" begin # STUB
end


text:
using Plots
# using Random
using LaTeXStrings
using JLD
# include("./cmd.jl")
# include("./simmd.jl")
# include("./extramd.jl")
# include("./pextramd.jl")
# include("./AutoDiffUtils.jl")
# include("./Potentials.jl")

lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"
brown = colorant"rgb(57%,34%,16%)"

α = 100
β = 1

#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(2.00)

#Extracting stored varibles
ld = load("./out/jld/regression_alpha_$(α)_beta_$(β).jld")
for key in keys(ld)
  symb = Symbol(key)
  data = ld[key]
  @eval $(symb) = $(data)
end

iterplot = plot(evals_cmd, label=L"\mathrm{CMW}", color=silver, title=latexstring("\\alpha = $(α), \\beta = $(β)"), linewidth=5, ylims=(0,40))
plot!(iterplot, evals_extragd, label=L"\mathrm{PX}", color=orange, linestyle=:dot, linewidth=5)
plot!(iterplot, evals_extramd, label=L"\mathrm{PXM}", color=rust, linestyle=:dash, linewidth=5)

savefig(iterplot, "./out/figures/regression_iterplot_alpha_$(α)_beta_$(β).pdf")

gradplot = plot(Int.(iters_cmd),  evals_cmd, label=L"\mathrm{CMW}", color=silver, title=latexstring("\\alpha = $(α), \\beta = $(β)"), linewidth=5, ylims=(0,40))
plot!(gradplot, Int.(iters_extragd), evals_extragd, label=L"\mathrm{PX}", color=orange, linestyle=:dot, linewidth=5)

plot!(gradplot, Int.(iters_extramd), evals_extramd, label=L"\mathrm{PXM}", color=rust, linestyle=:dash, linewidth=5)

savefig(gradplot, "./out/figures/regression_gradplot_alpha_$(α)_beta_$(β).pdf")

max_stars_repo_name: f-t-s/CMD

max_issues_repo_name: f-t-s/CMD

max_forks_repo_name: f-t-s/CMD


text:
using Kahuna
using Test
using SafeTestsets
using Libdl

push!(Libdl.DL_LOAD_PATH, "/usr/lib")

@safetestset "File I/O" begin include("test_fileio.jl") end


text:
export local_dimension, spacedimension
export nsites, toint, index, shape
export state_i

"""
    nsites(hilbert) -> Int

Returns the number of lattice sites in the Hilbert space.
"""
function nsites end

"""
    shape(hilbert) -> Vector{Int}

Returns a vector containing the local hilbert space dimensions of every
mode in the Hilbert space.

In the case of homogeneous spaces, it is usually more efficient to call
`local_dim`.
"""
function shape end

"""
    spacedimension(hilbert) = prod(shape(hilbert)) -> Int

Returns the total dimension of the vector space `hilbert`.
This is only valid if `indexable(hilbert) == true`, otherwise it's 0.
"""
function spacedimension end

"""
    indexable(hilbert) -> Bool

Returns true if the space is can be indexed with an Int64, which means that htis
tests true when

```
    spacedimension(hilbert) <= typemax(Int64)
```
"""
function indexable end

"""
    is_homogeneous(hilbert) -> Bool

Returns true if the space is homogeneous, that is, if all the modes have the
same local hilbert space dimension.
"""
function is_homogeneous end

# Global somewhat slow methods
"""
    local_dim(hilbert [, i]) -> Int

Returns the local dimension of the `hilbert` space on site `i`.
If the space is homogeneous then specifying `i` is not needed.
"""
@inline local_dim(h::AbstractHilbert, i) = shape(h)[i]

"""
    index(hilb, state)

Converts to an int the `state` of the hilbert space `hilb`
assuming it's indexable.
"""
@inline index(h::AbstractHilbert, s) = toint(s, h)

function index(h::AbstractHilbert, s::Union{AStateBatch, AStateBatchVec,ADoubleStateBatch, ADoubleStateBatchVec})
    sz = state_size(s)
    out = zeros(Int, sz...)
    for (i,sᵢ)=enumerate(states(s))
        @inbounds out[i] = index(h, sᵢ)
    end

    return out
end

"""
    state([arrT=Vector], [T=STD_REAL_PREC], hilbert, [i=0])

Constructs a state for the `hilbert` space with precision `T` and array type
`arrT`. By default that's the lowest state, otherwise if the hilbert space it's
indexable you can specify with a second argument.
"""
state(arrT::AbstractArray,  T::Type, h::AbstractHilbert, dims::Integer...) = state(arrT, T, h, dims)
state(arrT::AbstractArray,           h::AbstractHilbert, dims::Integer...) = state(arrT, STD_REAL_PREC, h, dims)
state(                      T::Type, h::AbstractHilbert, dims::Integer...) = state(zeros(0), T, h, dims)
state(                               h::AbstractHilbert, dims::Integer...) = state(STD_REAL_PREC, h, dims)

state(arrT::AbstractArray,  h::AbstractHilbert, dims::Dims) = state(arrT, STD_REAL_PREC, h, dims)
state(   T::Type{<:Number}, h::AbstractHilbert, dims::Dims) = state(zeros(0), T, h, dims)
state(                      h::AbstractHilbert, dims::Dims) = state(STD_REAL_PREC, h, dims)

Base.@propagate_inbounds state_i(h::AbstractHilbert, i::Int) = state_i(STD_REAL_PREC, h, i)
Base.@propagate_inbounds state_i(T::Type{<:Number}, h::AbstractHilbert, i::Int) = set!(state(T, h), h, i)

Base.@propagate_inbounds function apply!(σ::ADoubleState, h::AbstractHilbert, cngs_l, cngs_r)
    apply!(row(σ), h, cngs_l)
    apply!(col(σ), h, cngs_r)
    return σ
end

function apply!(σ::AbstractVector, h::AbstractHilbert, cngs)
    @boundscheck checkbounds(σ, site)

    for (site, val)=cngs
        @inbounds σ[site] = val
    end
    return σ
end

"""
    flipat!(state, hilb, i) -> (old_val, new_val)

Randomly flips `state[i]` to another available state. Returns the
old value of `state[i]` and the new value. state is changed in-place.
"""
Base.@propagate_inbounds flipat!(σ, hilb::AbstractHilbert, i::Int) = flipat!(GLOBAL_RNG, σ, hilb, i)

"""
    rand!([rng=GLOBAL_RNG], state, hilb)

Generates a random state of hilbert space and stores on the preallocated
state `state`. It can also be a batch of states.
Optionally you can pass the rng.
"""
@inline Random.rand!(σ::AbstractArray, h::AbstractHilbert) = rand!(GLOBAL_RNG, σ, h)
@inline Random.rand!(σ::NTuple{2,<:AbstractArray}, h::AbstractHilbert) = rand!(GLOBAL_RNG, σ, h)

# AbstractHilbert implementation for batches/batchvecs of states
function Random.rand!(rng::AbstractRNG, σ::AbstractArray, h::AbstractHilbert)
    for σᵢ=states(σ)
        rand!(rng, σᵢ, h)
    end
    return σ
end

"""
    rand([rng=GLOBAL_RNG], hilb, [batch_size, chain_length])

Generates a random state of hilbert space on the state.
Optionally you can pass the rng or the batch size/chain length
"""
@inline Random.rand(rng::AbstractRNG, h::AbstractHilbert, dims::Int...) = rand!(rng, state(h, dims...), h)
@inline Random.rand(rng::AbstractRNG, h::AbstractHilbert, dims::Dims)   = rand!(rng, state(h, dims), h)

@inline Random.rand(                  h::AbstractHilbert, dims::Int...) = rand(GLOBAL_RNG, h, dims...)
@inline Random.rand(                  h::AbstractHilbert, dims::Dims)   = rand(GLOBAL_RNG, h, dims)

@inline function checkbounds_hilbert(hilb::AbstractHilbert, i)
    (i >= 1 && i<= spacedimension(hilb)) || Base.throw_boundserror(hilb, i)
    return nothing
end


text:
struct PropaneRef <: EmpiricHelmholtzModel end

const PropaneRef_consts = (
    R = 8.314472 #J·mol-1·K-1
    ,Mw = 44.09562 #g·mol-1
    ,T_c = 369.89    #K
    ,P_c = 4.2512e6 #Pa
    ,rho_c= 5000.0 # mol·m-3
    ,T_triple = 85.525 #K
    ,P_Triple = 0.00017 #Pa
    ,rho_v_triple = 2.4e-07 # mol·m-3
    ,rho_l_triple = 16626.0 # mol·m-3
    ,T_0 = 231.036 #K
    ,rho_v_0 = 54.8 # mol·m-3
    ,rho_l_0 = 13173.0 # mol·m-3
    ,ref_T_ideal = 273.15
    ,ref_P_0 = 100000.0
    ,ref_T_ideal_H = 26148.48 ## J·mol-1
    ,ref_T_ideal_S = 157.9105 ##J·mol-1·K-1
    ,components = ["propane"]
    ,lengthcomponents=1
    )
const PropaneRef_a_consts = (
    v = (0.0,0.0,3.043,5.874,9.337,7.922)
    ,u = (0.0,0.0, 393.0, 1237.0,1984.0,4351)
    ,b = (0.0,0.0,1.062478, 3.344237,5.363757,11.762957)
    ,N = [0.042910051,
        1.7313671,
        -2.4516524,
        0.34157466,
        -0.46047898,
        -0.66847295,
        0.20889705,
        0.19421381,
        -0.22917851,
        -0.60405866,
        0.066680654,
        0.017534618,
        0.33874242,
        0.22228777,
        -0.23219062,
        -0.092206940,
        -0.47575718,
        -0.017486824]
    ,t = [1.00,
        0.33,
        0.80,
        0.43,
        0.90,
        2.46,
        2.09,
        0.88,
        1.09,
        3.25,
        4.62,
        0.76,
        2.50,
        2.75,
        3.05,
        2.55,
        8.40,
        6.75]
    ,d = [4,1,1,2,2,1,3,6,6,2,3,1,1,1,2,2,4,1]
    ,l = [1,1,1,1,2,2]
    ,η = [0.963,1.977,1.917,2.307,2.546,3.28,14.6]
    ,β = [2.33,3.47,3.15,3.19,0.92,18.8,547.8]
    ,γ = [0.684,0.829,1.419,0.817,1.500,1.426,1.093]
    ,ε = [1.283,0.6936,0.788,0.473,0.8577,0.271,0.948]

)

function _propane_ref_a0(δ,τ)
    δ,τ = promote(δ,τ)
    _1 = one(δ)
    a₁ = -4.970583
    a₂ = 4.29352
    b = (1.062478, 3.344237,5.363757,11.762957)
    v = (3.043,5.874,9.337,7.922)
    α₀ = log(δ) + 3*log(τ) + a₁ + a₂*τ +
     v[1]*log(_1-exp(-b[1]*τ)) +
     v[2]*log(_1-exp(-b[2]*τ)) +
     v[3]*log(_1-exp(-b[3]*τ)) +
     v[4]*log(_1-exp(-b[4]*τ))
     return α₀
end

function _propane_ref_ar(δ,τ)
    δ,τ = promote(δ,τ)
    N = PropaneRef_a_consts.N
    t = PropaneRef_a_consts.t
    d = PropaneRef_a_consts.d
    l = PropaneRef_a_consts.l
    η = PropaneRef_a_consts.η
    β = PropaneRef_a_consts.β
    γ = PropaneRef_a_consts.γ
    ε = PropaneRef_a_consts.ε
    αᵣ = zero(δ)

    for k in 1:5
        αᵣ += N[k]*(δ^d[k])*(τ^t[k])
    end

    for k in 6:11
        αᵣ += N[k]*(δ^d[k])*(τ^t[k])*exp(-δ^l[k-5])
    end

    for k in 12:18
        _k = k - 11
        mul = exp(-η[_k]*(δ-ε[_k])^2  - β[_k]*(τ-γ[_k])^2)
        αᵣ += N[k]*(δ^d[k])*(τ^t[k])*mul
    end
    return αᵣ
end
#ancillary equations for calculation of P_sat, rhovsat y rholsat
function _propaneref_tsat(T)
    T_c = PropaneRef_consts.T_c
    P_c = PropaneRef_consts.P_c
    Tr = T/T_c
    θ = 1.0-Tr
    lnPsatPc = (-6.7722*θ + 1.6938*θ^1.5 -1.3341*θ^2.2 -3.1876*θ^4.8 + 0.94937*θ^6.2)/Tr
    Psat = exp(lnPsatPc)*P_c
    return Psat
end

function _propaneref_rholsat(T)
    Tc = PropaneRef_consts.T_c
    ρ_c =PropaneRef_consts.rho_c
    Tr = T/Tc
    θ = 1.0-Tr
    #
    ρ_l = (1.0 + 1.82205*θ^0.345 + 0.65802*θ^0.74 + 0.21109*θ^2.6 + 0.083973*θ^7.2)*ρ_c
    return ρ_l
end

function _propaneref_rhovsat(T)
    Tc = PropaneRef_consts.T_c
    ρ_c =PropaneRef_consts.rho_c
    Tr = T/Tc
    θ = 1.0 - Tr
    log_ρ_v_ρ_c = (-2.4887*θ^0.3785 -5.1069*θ^1.07 -12.174*θ^2.7 -30.495*θ^5.5 -52.192*θ^10 -134.89*θ^20)
    ρ_v = exp(log_ρ_v_ρ_c)*ρ_c
    return ρ_v
end

function a_scaled(model::PropaneRef,δ,τ)
    return  _propane_ref_a0(δ,τ)+_propane_ref_ar(δ,τ)
end

function eos(model::PropaneRef, V, T, z=SA[1.0];phase=:unknown)
    R =PropaneRef_consts.R
    T_c = PropaneRef_consts.T_c
    rho_c = PropaneRef_consts.rho_c
    N = only(z)
    rho = (N/V)
    δ = rho/rho_c
    τ = T_c/T
    return N*R*T*a_scaled(model::PropaneRef,δ,τ)
end

function eos_res(model::PropaneRef,V,T,z=SA[1.0];phase=:unknown)
    R =PropaneRef_consts.R
    T_c = PropaneRef_consts.T_c
    rho_c = PropaneRef_consts.rho_c
    N = only(z)
    rho = (N/V)
    δ = rho/rho_c
    τ = T_c/T
    return N*R*T*_propane_ref_ar(δ,τ)
end

mw(model::PropaneRef) = SA[PropaneRef_consts.Mw]
molecular_weight(model::PropaneRef,z = @SVector [1.]) = PropaneRef_consts.Mw*0.001
T_scale(model::PropaneRef,z=SA[1.0]) = PropaneRef_consts.T_c
p_scale(model::PropaneRef,z=SA[1.0]) = PropaneRef_consts.P_c
lb_volume(model::PropaneRef,z=SA[1.0]; phase=:l) = 6.0647250138479785e-5 #calculated at 1000 MPa and 650 K

function Base.show(io::IO,mime::MIME"text/plain",model::PropaneRef)
    print(io,"Propane Reference Equation of State")
end
function x0_sat_pure(model::PropaneRef,T,z=SA[1.0])
    log10vv = log10(1.0/_propaneref_rhovsat(T))
    log10vl = log10(1.0/_propaneref_rholsat(T))
    return [log10vl,log10vv]
end

function Base.getproperty(model::PropaneRef,sym::Symbol)
    return PropaneRef_consts[sym]
end

export PropaneRef


text:
module Qux

using Hyperspecialize

struct Wobble <: Real
  w::Int
end

f(::Real) = false
g(::Real) = false
@concretize TypicalTag []
@replicable f(::@hyperspecialize(TypicalTag)) = true
@replicable g(::@hyperspecialize(TypicalTag)) = true

h(::Real, ::Wobble) = false

export h, Wobble

end #module


text:
# BinaryLayer and BinaryMERA types, and methods thereof.
# To be `included` in MERAKit.jl.

"""
    BinaryLayer{ST, ET, Tan} <: SimpleLayer

The type for layers of a binary MERA.

Each layer consists of two tensors, a 2-to-2 disentangler, often called `u`, and a 2-to-1
isometry, often called `w`.

The type parameters are `ST` for space type, e.g. `ComplexSpace` or `SU2Space`; `ET` for
element type, e.g. `Complex{Float64}`; and `Tan` for whether this layer is a tangent layer
or not.  If `Tan = false`, the layer in question is an actual MERA layer. If `Tan = true` it
consists, instead of the actual tensors, of Stiefel/Grassmann tangent vectors of these
tensors.

`BinaryLayer` implements the iteration interface, returning first the disentangler, then the
isometry.

Index numbering convention is as follows, where the physical indices are at the bottom:
Disentangler:
```
 3│    4│
 ┌┴─────┴┐
 │   u   │
 └┬─────┬┘
 1│    2│
```

Isometry:
```
    3│
 ┌───┴───┐
 │   w   │
 └┬─────┬┘
 1│    2│
```
"""
struct BinaryLayer{ST, ET, Tan} <: SimpleLayer
    # Even though the types are here marked as any, they are restricted to be specific,
    # concrete types dependent on ST, ET and Tan, both in the constructor and in
    # getproperty.
    disentangler::Any
    isometry::Any

    function BinaryLayer{ST, ET, Tan}(disentangler, isometry) where {ST, ET, Tan}
        DisType = disentangler_type(ST, ET, Tan)
        IsoType = binaryisometry_type(ST, ET, Tan)
        disconv = convert(DisType, disentangler)::DisType
        isoconv = convert(IsoType, isometry)::IsoType
        return new{ST, ET, Tan}(disconv, isoconv)
    end
end

function BinaryLayer(disentangler::DisType, isometry::IsoType) where {
    ST,
    DisType <: AbstractTensorMap{ST, 2, 2},
    IsoType <: AbstractTensorMap{ST, 2, 1}
}
    ET = eltype(DisType)
    @assert eltype(IsoType) === ET
    return BinaryLayer{ST, ET, false}(disentangler, isometry)
end

function BinaryLayer(disentangler::DisTanType, isometry::IsoTanType) where {
    ST,
    DisType <: AbstractTensorMap{ST, 2, 2},
    IsoType <: AbstractTensorMap{ST, 2, 1},
    DisTanType <: Stiefel.StiefelTangent{DisType},
    IsoTanType <: Grassmann.GrassmannTangent{IsoType},
}
    ET = eltype(DisType)
    @assert eltype(IsoType) === ET
    return BinaryLayer{ST, ET, true}(disentangler, isometry)
end

function Base.getproperty(l::BinaryLayer{ST, ET, Tan}, sym::Symbol) where {ST, ET, Tan}
    if sym === :disentangler
        T = disentangler_type(ST, ET, Tan)
    elseif sym === :isometry
        T = binaryisometry_type(ST, ET, Tan)
    else
        T = Any
    end
    return getfield(l, sym)::T
end

"""
    BinaryMERA{N}

A binary MERA is a MERA consisting of `BinaryLayer`s.
"""
BinaryMERA{N} = GenericMERA{N, T, O} where {T <: BinaryLayer, O}
layertype(::Type{BinaryMERA}) = BinaryLayer
#Base.show(io::IO, ::Type{BinaryMERA}) = print(io, "BinaryMERA")
#Base.show(io::IO, ::Type{BinaryMERA{N}}) where {N} = print(io, "BinaryMERA{($N)}")

# Implement the iteration and indexing interfaces. Allows things like `u, w = layer`.
# See simplelayer.jl for details.
_tuple(layer::BinaryLayer) = (layer.disentangler, layer.isometry)

function operatortype(::Type{BinaryLayer{ST, ET, false}}) where {ST, ET}
    return tensormaptype(ST, 3, 3, ET)
end
operatortype(::Type{BinaryLayer{ST, ET, true}}) where {ST, ET} = Nothing

TensorKit.spacetype(::Type{BinaryLayer{ST, ET, Tan}}) where {ST, ET, Tan} = ST

scalefactor(::Type{<:BinaryLayer}) = 2

causal_cone_width(::Type{<:BinaryLayer}) = 3

Base.eltype(::Type{BinaryLayer{ST, ET, Tan}}) where {ST, ET, Tan} = ET

outputspace(layer::BinaryLayer) = space(layer.disentangler, 1)
inputspace(layer::BinaryLayer) = space(layer.isometry, 3)'
internalspace(layer::BinaryLayer) = space(layer.isometry, 1)
internalspace(m::BinaryMERA, depth) = internalspace(getlayer(m, depth))

function expand_inputspace(layer::BinaryLayer, V_new)
    u, w = layer
    w = pad_with_zeros_to(w, 3 => V_new')
    return BinaryLayer(u, w)
end

function expand_outputspace(layer::BinaryLayer, V_new)
    u, w = layer
    u = pad_with_zeros_to(u, 1 => V_new, 2 => V_new)
    return BinaryLayer(u, w)
end

function expand_internalspace(layer::BinaryLayer, V_new)
    u, w = layer
    u = pad_with_zeros_to(u, 3 => V_new', 4 => V_new')
    w = pad_with_zeros_to(w, 1 => V_new, 2 => V_new)
    return BinaryLayer(u, w)
end

function randomlayer(
    ::Type{BinaryLayer}, ::Type{T}, Vin, Vout, Vint = Vout; random_disentangler = false
) where {T}
    w = randomisometry(T, Vint ⊗ Vint, Vin)
    u = initialize_disentangler(T, Vout, Vint, random_disentangler)
    return BinaryLayer(u, w)
end

function ascending_fixedpoint(layer::BinaryLayer)
    width = causal_cone_width(layer)
    Vtotal = ⊗(ntuple(n->inputspace(layer), Val(width))...)
    return id(storagetype(operatortype(layer)), Vtotal)
end

function scalingoperator_initialguess(l::BinaryLayer, irreps...)
    width = causal_cone_width(l)
    V = inputspace(l)
    interlayer_space = ⊗(ntuple(n->V, Val(width))...)
    outspace = interlayer_space
    inspace = interlayer_space
    for irrep in irreps
        inspace = inspace ⊗ spacetype(V)(irrep => 1)
    end
    typ = eltype(l)
    t = TensorMap(randn, typ, outspace ← inspace)
    return t
end

function gradient(layer::BinaryLayer, env::BinaryLayer; metric = :euclidean)
    u, w = layer
    uenv, wenv = env
    # The environment is the partial derivative. We need to turn that into a tangent vector
    # of the Stiefel manifold point u or w.
    # The factor of two is from the partial_x + i partial_y derivative of the cost function,
    # and how it depends on both v and v^dagger.
    ugrad = Stiefel.project!(2*uenv, u; metric = metric)
    wgrad = Grassmann.project!(2*wenv, w)
    return BinaryLayer(ugrad, wgrad)
end

function precondition_tangent(layer::BinaryLayer, tan::BinaryLayer, rho)
    u, w = layer
    utan, wtan = tan
    @planar rho_wl[-1; -11] := rho[-1 1 2; -11 1 2]
    @planar rho_wm[-1; -11] := rho[1 -1 2; 1 -11 2]
    @planar rho_wr[-1; -11] := rho[1 2 -1; 1 2 -11]
    rho_w = (rho_wl + rho_wm + rho_wr) / 3.0
    @planar rho_twosite_l[-1 -2; -11 -12] := rho[-1 -2 1; -11 -12 1]
    @planar rho_twosite_r[-1 -2; -11 -12] := rho[1 -1 -2; 1 -11 -12]
    rho_twosite = (rho_twosite_l + rho_twosite_r) / 2.0
    @planar(
        rho_u[-1 -2; -11 -12] :=
        w'[12; 1 -11] * w'[22; -12 2] *
        rho_twosite[11 21; 12 22] *
        w[1 -1; 11] * w[-2 2; 21]
    )
    utan_prec = precondition_tangent(utan, rho_u)
    wtan_prec = precondition_tangent(wtan, rho_w)
    return BinaryLayer(utan_prec, wtan_prec)
end

# # # Invariants

function space_invar_intralayer(layer::BinaryLayer)
    u, w = layer
    matching_bonds = (
        (space(u, 3)', space(w, 2)),
        (space(u, 4)', space(w, 1))
    )
    allmatch = all(pair -> ==(pair...), matching_bonds)
    # Check that the dimensions are such that isometricity can hold.
    allmatch &= all((u, w)) do v
        codom, dom = fuse(codomain(v)), fuse(domain(v))
        infimum(dom, codom) == dom
    end
    return allmatch
end

function space_invar_interlayer(layer::BinaryLayer, next_layer::BinaryLayer)
    u, w = layer.disentangler, layer.isometry
    unext, wnext = next_layer.disentangler, next_layer.isometry
    matching_bonds = (
        (space(w, 3)', space(unext, 1)),
        (space(w, 3)', space(unext, 2))
    )
    allmatch = all(pair -> ==(pair...), matching_bonds)
    return allmatch
end

# # # Ascending and descending superoperators
const BinaryOperator{S} = AbstractTensorMap{S, 3, 3}
const ChargedBinaryOperator{S} = AbstractTensorMap{S, 3, 4}
const DoubleChargedBinaryOperator{S} = AbstractTensorMap{S, 3, 5}

function ascend(
    op::Union{BinaryOperator, ChargedBinaryOperator, DoubleChargedBinaryOperator},
    layer::BinaryLayer
)
    l = ascend_left(op, layer)
    r = ascend_right(op, layer)
    scaled_op = (l+r)/2
    return scaled_op
end

function ascend(
    op::Union{SquareTensorMap{1}, SquareTensorMap{2}}, layer::BinaryLayer
)
    return ascend(expand_support(op, causal_cone_width(BinaryLayer)), layer)
end

# TODO Think about how to best remove the code duplication of having the separate methods
# for ordinary, charged, and double charged operators. Note also that there's a mix of using
# @tensor and @planar.
function ascend_left(op::ChargedBinaryOperator, layer::BinaryLayer)
    u, w = layer
    @tensor(
        scaled_op[-100 -200 -300; -400 -500 -600 -1000] :=
        w[5 6; -400] * w[9 8; -500] * w[16 15; -600] *
        u[1 2; 6 9] * u[10 12; 8 16] *
        op[3 4 14; 1 2 10 -1000] *
        u'[7 13; 3 4] * u'[11 17; 14 12] *
        w'[-100; 5 7] * w'[-200; 13 11] * w'[-300; 17 15]
    )
    return scaled_op
end

function ascend_right(op::ChargedBinaryOperator, layer::BinaryLayer)
    u, w = layer
    @tensor(
        scaled_op[-100 -200 -300; -400 -500 -600 -1000] :=
        w[15 16; -400] * w[8 9; -500] * w[6 5; -600] *
        u[12 10; 16 8] * u[1 2; 9 6] *
        op[14 3 4; 10 1 2 -1000] *
        u'[17 11; 12 14] * u'[13 7; 3 4] *
        w'[-100; 15 17] * w'[-200; 11 13] * w'[-300; 7 5]
    )
    return scaled_op
end

function ascend_left(op::BinaryOperator, layer::BinaryLayer)
    u, w = layer
    @planar(
        scaled_op[-100 -200 -300; -400 -500 -600] :=
        w[5 6; -400] * w[9 8; -500] * w[16 15; -600] *
        u[1 2; 6 9] * u[10 12; 8 16] *
        op[3 4 14; 1 2 10] *
        u'[7 13; 3 4] * u'[11 17; 14 12] *
        w'[-100; 5 7] * w'[-200; 13 11] * w'[-300; 17 15]
    )
    return scaled_op
end

function ascend_right(op::BinaryOperator, layer::BinaryLayer)
    u, w = layer
    @planar(
        scaled_op[-100 -200 -300; -400 -500 -600] :=
        w[15 16; -400] * w[8 9; -500] * w[6 5; -600] *
        u[12 10; 16 8] * u[1 2; 9 6] *
        op[14 3 4; 10 1 2] *
        u'[17 11; 12 14] * u'[13 7; 3 4] *
        w'[-100; 15 17] * w'[-200; 11 13] * w'[-300; 7 5]
    )
    return scaled_op
end

function ascend_left(op::DoubleChargedBinaryOperator, layer::BinaryLayer)
    u, w = layer
    @planar temp1[-1 -2; -3 -4 -5 -6 -7 -8] := op[-1 -2 -8; 1 2 -5 -6 -7] * u[1 2; -3 -4]
    temp2 = braid(temp1, (11, 12, 13, 14, 15, 16, 1, 100), (1, 2), (3, 6, 7, 4, 5, 8))
    @planar(
        temp3[-100 -200 -300; -400 -1000 -2000 -500 -600] :=
        w[5 6; -400] * w[9 8; -500] * w[16 15; -600] *
        u[10 12; 8 16] *
        temp2[3 4; 6 -1000 -2000 9 10 14] *
        u'[7 13; 3 4] * u'[11 17; 14 12] *
        w'[-100; 5 7] * w'[-200; 13 11] * w'[-300; 17 15]
    )
    scaled_op = braid(temp3, (11, 12, 13, 14, 1, 100, 15, 16), (1, 2, 3), (4, 7, 8, 5, 6))
    return scaled_op
end

function ascend_right(op::DoubleChargedBinaryOperator, layer::BinaryLayer)
    u, w = layer
    @planar temp1[-1 -2; -3 -4 -5 -6 -7 -8] := op[-3 -1 -2; -4 1 2 -7 -8] * u[1 2; -5 -6]
    temp2 = braid(temp1, (11, 12, 13, 14, 15, 16, 1, 100), (1, 2), (3, 4, 5, 7, 8, 6))
    @planar(
        temp3[-100 -200 -300; -400 -500 -1000 -2000 -600] :=
        w[15 16; -400] * w[8 9; -500] * w[6 5; -600] *
        u[12 10; 16 8] *
        temp2[3 4; 14 10 9 -1000 -2000 6] *
        u'[17 11; 12 14] * u'[13 7; 3 4] *
        w'[-100; 15 17] * w'[-200; 11 13] * w'[-300; 7 5]
    )
    scaled_op = braid(temp3, (11, 12, 13, 14, 15, 1, 100, 16), (1, 2, 3), (4, 5, 8, 6, 7))
    return scaled_op
end

function descend_left(rho::BinaryOperator, layer::BinaryLayer)
    u, w = layer
    @planar(
        scaled_rho[-100 -200 -300; -400 -500 -600] :=
        u'[16 17; -400 -500] * u'[2 10; -600 11] *
        w'[12; 1 16] * w'[9; 17 2] * w'[5; 10 4] *
        rho[13 7 6; 12 9 5] *
        w[1 14; 13] * w[15 3; 7] * w[8 4; 6] *
        u[-100 -200; 14 15] * u[-300 11; 3 8]
    )
    return scaled_rho
end

function descend_right(rho::BinaryOperator, layer::BinaryLayer)
    u, w = layer
    @planar(
        scaled_rho[-100 -200 -300; -400 -500 -600] :=
        u'[10 2; 11 -400] * u'[17 16; -500 -600] *
        w'[5; 4 10] * w'[9; 2 17] * w'[12; 16 1] *
        rho[6 7 13; 5 9 12] *
        w[4 8; 6] * w[3 15; 7] * w[14 1; 13] *
        u[11 -100; 8 3] * u[-200 -300; 15 14]
    )
    return scaled_rho
end

function descend(rho::BinaryOperator, layer::BinaryLayer)
    l = descend_left(rho, layer)
    r = descend_right(rho, layer)
    scaled_rho = (l+r)/2
    return scaled_rho
end

# # # Optimization

function environment(op, layer::BinaryLayer, rho; vary_disentanglers = true)
    if vary_disentanglers
        env_u = environment_disentangler(op, layer, rho)
    else
        # The adjoint is just for type stability.
        env_u = zero(layer.disentangler')'
    end
    env_w = environment_isometry(op, layer, rho)
    return BinaryLayer(env_u, env_w)
end

function minimize_expectation_ev(
    layer::BinaryLayer, env::BinaryLayer; vary_disentanglers = true
)
    u = if vary_disentanglers
        projectisometric(env.disentangler; alg = Polar())
    else
        layer.disentangler
    end
    w = projectisometric(env.isometry; alg = Polar())
    return BinaryLayer(u, w)
end

function environment_disentangler(h::BinaryOperator, layer::BinaryLayer, rho)
    u, w = layer
    @planar(
        env1[-1 -2; -3 -4] :=
        rho[17 18 10; 15 14 9] *
        w[5 7; 17] * w[12 11; 18] * w[19 8; 10] *
        u[3 4; 7 12] * u[13 -2; 11 19] *
        h[1 2 -1; 3 4 13] *
        u'[6 16; 1 2] *
        w'[15; 5 6] * w'[14; 16 -3] * w'[9; -4 8]
    )

    @planar(
        env2[-1 -2; -3 -4] :=
        rho[4 15 6; 3 10 5] *
        w[1 13; 4] * w[14 16; 15] * w[17 2; 6] *
        u[12 18; 13 14] * u[7 8; 16 17] *
        h[19 -1 -2; 18 7 8] *
        u'[11 9; 12 19] *
        w'[3; 1 11] * w'[10; 9 -3] * w'[5; -4 2]
    )

    @planar(
        env3[-1 -2; -3 -4] :=
        rho[6 15 4; 5 10 3] *
        w[2 17; 6] * w[16 14; 15] * w[13 1; 4] *
        u[8 7; 17 16] * u[18 12; 14 13] *
        h[-1 -2 19; 8 7 18] *
        u'[9 11; 19 12] *
        w'[5; 2 -3] * w'[10; -4 9] * w'[3; 11 1]
    )

    @planar(
        env4[-1 -2; -3 -4] :=
        rho[10 18 17; 9 14 15] *
        w[8 19; 10] * w[11 12; 18] * w[7 5; 17] *
        u[-1 13; 19 11] * u[4 3; 12 7] *
        h[-2 2 1; 13 4 3] *
        u'[16 6; 2 1] *
        w'[9; 8 -3] * w'[14; -4 16] * w'[15; 6 5]
    )

    env = (env1 + env2 + env3 + env4)/2
    return env
end

function environment_disentangler(
    h::Union{SquareTensorMap{1}, SquareTensorMap{2}}, layer::BinaryLayer, rho
)
    h = expand_support(h, causal_cone_width(BinaryLayer))
    return environment_disentangler(h, layer, rho)
end

function environment_isometry(h::BinaryOperator, layer, rho)
    u, w = layer
    @planar(
        env1[-1 -2; -3] :=
        rho[16 15 19; 18 17 -3] *
        w[5 7; 16] * w[14 13; 15] * w[20 -2; 19] *
        u[4 3; 7 14] * u[12 11; 13 20] *
        h[2 1 10; 4 3 12] *
        u'[6 9; 2 1] * u'[8 -1; 10 11] *
        w'[18; 5 6] * w'[17; 9 8]
    )

    @planar(
        env2[-1 -2; -3] :=
        rho[18 17 19; 16 15 -3] *
        w[12 14; 18] * w[11 10; 17] * w[20 -2; 19] *
        u[9 8; 14 11] * u[4 3; 10 20] *
        h[7 2 1; 8 4 3] *
        u'[13 5; 9 7] * u'[6 -1; 2 1] *
        w'[16; 12 13] * w'[15; 5 6]
    )

    @planar(
        env3[-1 -2; -3] :=
        rho[19 20 15; 18 -3 14] *
        w[5 7; 19] * w[10 8; 20] * w[16 13; 15] *
        u[4 3; 7 10] * u[9 11; 8 16] *
        h[2 1 12; 4 3 9] *
        u'[6 -1; 2 1] * u'[-2 17; 12 11] *
        w'[18; 5 6] * w'[14; 17 13]
    )

    @planar(
        env4[-1 -2; -3] :=
        rho[15 20 19; 14 -3 18] *
        w[13 16; 15] * w[8 10; 20] * w[7 5; 19] *
        u[11 9; 16 8] * u[3 4; 10 7] *
        h[12 1 2; 9 3 4] *
        u'[17 -1; 11 12] * u'[-2 6; 1 2] *
        w'[14; 13 17] * w'[18; 6 5]
    )

    @planar(
        env5[-1 -2; -3] :=
        rho[19 17 18; -3 15 16] *
        w[-1 20; 19] * w[10 11; 17] * w[14 12; 18] *
        u[3 4; 20 10] * u[8 9; 11 14] *
        h[1 2 7; 3 4 8] *
        u'[-2 6; 1 2] * u'[5 13; 7 9] *
        w'[15; 6 5] * w'[16; 13 12]
    )

    @planar(
        env6[-1 -2; -3] :=
        rho[19 15 16; -3 17 18] *
        w[-1 20; 19] * w[13 14; 15] * w[7 5; 16] *
        u[11 12; 20 13] * u[3 4; 14 7] *
        h[10 1 2; 12 3 4] *
        u'[-2 8; 11 10] * u'[9 6; 1 2] *
        w'[17; 8 9] * w'[18; 6 5]
    )

    env = (env1 + env2 + env3 + env4 + env5 + env6)/2
    return env
end

function environment_isometry(
    h::Union{SquareTensorMap{1}, SquareTensorMap{2}}, layer::BinaryLayer, rho
)
    h = expand_support(h, causal_cone_width(BinaryLayer))
    return environment_isometry(h, layer, rho)
end


text:
function write_input(	file_path::String,
						output_generation_path::String,
						number_of_samples::Int64,
						d2::Array{Float64, 1},
						d3::Array{Float64, 1},
						theta3::Array{Float64, 1},
						number_of_cells_x::Int64,
						number_of_cells_y::Int64,
						number_of_cells_z::Int64,
						output_file_path::String)

	file_stream::IOStream = open(file_path, "w")

	@printf(file_stream, "%s", "<input_characterization>\n")

	write_key(file_stream, "output_generation_path", output_generation_path)
	write_key(file_stream, "number_of_samples", number_of_samples)
	write_key(file_stream, "d2", d2)
	write_key(file_stream, "d3", d3)
	write_key(file_stream, "theta3", theta3)
	write_key(file_stream, "number_of_cells_x", number_of_cells_x)
	write_key(file_stream, "number_of_cells_y", number_of_cells_y)
	write_key(file_stream, "number_of_cells_z", number_of_cells_z)
	write_key(file_stream, "output_file_path", output_file_path)

	@printf(file_stream, "%s", "</input_characterization>")

	close(file_stream)

	nothing
end


text:
using Gen, Genify
using Distributions: Normal

##  Define model of noisy weighing scale in Julia

function scale(n::Int)
    mass = rand(Normal(5, 1))
    obs = measure(mass, n)
end

function measure(mass::Real, n::Int)
    obs = Float64[]
    for i in 1:n
        m = rand(Normal(mass, 2))
        push!(obs, m)
    end
    return obs
end

## Genify the model

genscale = genify(scale, Int)

## Construct choicemap of noisy measurements

observations = choicemap(((:obs => :m => i, 4) for i in 1:10)...)

## Perform importance sampling to estimate the true mass

trs, ws, _ = importance_sampling(genscale, (10,), observations, 100)
mass_est = sum(tr[:mass] * exp(w) for (tr, w) in zip(trs, ws))


text:
struct PNInput
    pInterval::Interval
    nInterval::Interval
    nMaterial::SemiconductorData
    pMaterial::SemiconductorData

end

PNInput(T,N,L) = PNInput(Domain(Interval{:close,:open}(0,L/2.0)),Domain(Interval{:close,:close}(L/2.0,L)),SemiconductorData(T,BoltzmanDist(),NSilicon(N,0)),SemiconductorData(T,BoltzmanDist(),PSilicon(N,0)))

function PoissonEq(PN::PNInput)
    F(x) = DistFun(PN.pMaterial,x)
    @variables x ψ(..)
    @parameters δ λ Cfun(x)
    Dxx = Differential(x)^2
    eqs = [
        λ*Dxx(ψ(x)) ~ δ*F(ψ(x))-p(ψ(x))-Cfun(x*lscal)
    ]
    return eqs
end

function Poisson(Vg,PN::PNInput,Nres = 21.0)
    HeavisideC(x) = (1+sign(x+eps(Float64)))*0.5
    Heaviside(x) = (1+sign(x-eps(Float64)))*0.5
    Il(x,i::Interval) = HeavisideC(x-i.left)-Heaviside(x-i.right)
    Ir(x,i::Interval) = HeavisideC(x-i.left)-Heaviside(x-i.right)

    Vscal,lscalp,CN,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i = scalling(PN.pMaterial,PN.pInterval.right)
    Vscal,lscaln,CP,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i = scalling(PN.nMaterial,PN.nInterval.right)
    Cscal = max(abs.(CN),abs.(CP))
    lscal = max(lscaln,lscalp)
    f(x,a,b,Nn,Np) = C(Nn)*Il(x,a) + C(Np)*Ir(x,b)
    g(x) = f(x,PN.pInterval,PN.nInterval,PN.pMaterial,PN.nMaterial)

    ul = -Ei(PN.pMaterial)/Vscal
    ur = -Ei(PN.nMaterial)/Vscal

    #C0 = Cdoping./Cscal

    dx = Ld0/Nres/lscal

    F(x) = DistFun(PN.pMaterial,x)
    n(ψ) = δ*F(ψ)
    p(ψ) = δ*F(-ψ)
    @variables x ψ(..)

    Dxx = Differential(x)^2
    Dx = Differential(x)

    Cfun(x) = g(x)/Cscal
    eqs = [
        λ*Dxx(ψ(x)) ~ n(ψ(x))-p(ψ(x))-Cfun(x*lscal)
    ]
    bcs = [ψ(0) ~ ul+Vg/Vscal,
            ψ(1.0) ~ ur]
    domain = [x ∈ (0.0,1.0)]
    pde_system = PDESystem(eqs,bcs,domain,[x],[ψ(x)])
    discretization = MOLFiniteDifference([x=>dx],nothing;centered_order=2)
    prob = discretize(pde_system,discretization)
    L = length(prob.u0)
    L ÷ 2
    prob.u0[1:L÷2] .= ul
    prob.u0[1] = ul+ Vg/Vscal
    prob.u0[L÷2:L] .= ur
    prob.u0
    sol = solve(prob,jac=true)*Vscal
    if iseven(L)
        dxP = 0
    else
        dxP = dx*lscal
    end
    xP = PN.pInterval.left:dx*lscal:(PN.pInterval.right-dxP)
    LP = length(xP)
    SDP = Semiconductor1D(PN.pMaterial,xP,sol[1:LP].+Ei(PN.pMaterial))
    xN = PN.nInterval.left:dx*lscal:PN.nInterval.right
    LN = length(xN)
    SDN = Semiconductor1D(PN.nMaterial,xN,sol[LP+1:LP+LN].+Ei(PN.nMaterial))
    return Semiconductor1D(vcat(SDP.points,SDN.points))

end

function DD1D(Vg,PN::PNInput,Nres=21.0)
    HeavisideC(x) = (1+sign(x+eps(Float64)))*0.5
    Heaviside(x) = (1+sign(x-eps(Float64)))*0.5
    Il(x,i::Interval) = HeavisideC(x-i.left)-Heaviside(x-i.right)
    Ir(x,i::Interval) = HeavisideC(x-i.left)-Heaviside(x-i.right)

    Vscal,lscalp,CN,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i = scalling(PN.pMaterial,PN.pInterval.right)
    Vscal,lscaln,CP,μscal,Dscal,Jscal,Rscal,tscal,δ,λ,Ld0,n_i = scalling(PN.nMaterial,PN.nInterval.right)
    Cscal = max(abs.(CN),abs.(CP))
    lscal = max(lscaln,lscalp)
    f(x,a,b,Nn,Np) = C(Nn)*Il(x,a) + C(Np)*Ir(x,b)
    g(x) = f(x,PN.pInterval,PN.nInterval,PN.pMaterial,PN.nMaterial)

    τ = 1e-6
    τn = τ/tscal
    τp = τ/tscal

    un = 1400/μscal
    up = 450/μscal

    Rfunc(n,p) = δ/(τn*(p/δ+1)+τp*(n/δ+1))

    dx = Ld0/Nres/lscal

    F(x) = DistFun(PN.pMaterial,x)
    nf(ψ,ϕ) = δ*F(ψ-ϕ)
    pf(ψ,ϕ) = δ*F(-ψ+ϕ)

    @variables x u(..) ϕn(..) ϕp(..)

    Dxx = Differential(x)^2
    Dx = Differential(x)

    Cfun(x) = g(x)/Cscal
    ul = -Ei(PN.pMaterial)/Vscal
    ur = -Ei(PN.nMaterial)/Vscal

    #nD0 = (n(PN.pMaterial)+Nd(PN.pMaterial))/Cscal
    nD0 = (n(PN.pMaterial))/Cscal
    ψn0 = find_zero(x->nf(ul,x/Vscal)-nD0,(-100.0,10.0))*Vscal +  Vg/Vscal

    #pD0 = (p(PN.pMaterial)+Na(PN.pMaterial))/Cscal
    pD0 = (p(PN.pMaterial))/Cscal
    pD0 = ni(PN.pMaterial)
    ψp0 = find_zero(x->pf(ul,x/Vscal)-nD0,(-100.0,10.0))*Vscal + Vg/Vscal

    #nDL = (n(PN.nMaterial)+Nd(PN.nMaterial))/Cscal
    nDL = (n(PN.nMaterial))/Cscal
    ψnL = find_zero(x->nf(ul,x/Vscal)-nDL,(-100.0,100.0))*Vscal

    #pDL = (p(PN.nMaterial)+Na(PN.nMaterial))/Cscal
    pDL = (p(PN.nMaterial))/Cscal
    ψpL = find_zero(x->pf(ul,x/Vscal)-nDL,(-100.0,100.0))*Vscal

    eqs = [
        λ*Dxx(u(x)) ~ nf(u(x),0)-pf(u(x),0)-Cfun(x*lscal),
        Dxx(ϕn(x)) ~ -Rfunc(nf(u(x),ϕn(x)),pf(u(x),0.0))/nf(u(x),ϕn(x))/un+Dx(ϕn(x))*(Dx(ϕn(x)) - Dx(u(x))),
        Dxx(ϕp(x)) ~  Rfunc(nf(u(x),ϕn(x)),pf(u(x),ϕp(x)))/pf(u(x),ϕp(x))/up-Dx(ϕp(x))*(Dx(ϕp(x)) - Dx(u(x)))
        #Dx(up*p(ψ(x),ϕp(x)))*Dx(ϕp(x))+(up*p(ψ(x),ϕp(x)))*Dxx(ϕp(x)) ~ -R(n(ψ(x),ϕn(x)),p(ψ(x),ϕp(x)))
    ]

    bcs = [u(0.0) ~ ul,
           u(1.0) ~ ur,
           ϕn(0.0) ~ ψn0,
           ϕn(1.0) ~ ψnL,
           ϕp(0.0) ~ ψp0,
           ϕp(1.0) ~ ψpL,
            ]
        domain = [x ∈ (0.0,1.0)]
        pde_system = PDESystem(eqs,bcs,domain,[x],[u(x),ϕn(x),ϕp(x)])
        discretization = MOLFiniteDifference([x=>dx],nothing;centered_order=2)
        prob = discretize(pde_system,discretization)
        L = length(prob.u0)÷3
        prob.u0[1:L÷2] .= ul+ Vg/Vscal
        prob.u0[L÷2+1:L] .= ur
        sol = solve(prob,jac=true)*Vscal
        if iseven(L)
            dxP = 0
        else
            dxP = dx*lscal
        end
        ψsol = sol[1:L]
        ϕnsol = sol[L+1:2L]
        ϕpsol = sol[2L+1:3L]

        xP = PN.pInterval.left:dx*lscal:(PN.pInterval.right-dxP)
        LP = length(xP)
        SDP = Semiconductor1D(PN.pMaterial,xP,ψsol[1:LP].+Ei(PN.pMaterial),ϕnsol[1:LP],ϕpsol[1:LP])
        xN = PN.nInterval.left:dx*lscal:PN.nInterval.right
        LN = length(xN)
        SDN = Semiconductor1D(PN.nMaterial,xN,ψsol[LP+1:LP+LN].+Ei(PN.nMaterial),ϕnsol[LP+1:LP+LN],ϕpsol[LP+1:LP+LN])

        return Semiconductor1D(vcat(SDP.points,SDN.points))

end