text:
module CTImageData

export AbstractCTData, CTData

using ...Monads
using ...CTImages
import ...CTImages: ctimage, ctsinogram, cttomogram
import ...Utils: _atype
import Base: similar, eltype

abstract type AbstractCTData end

# Container type of CTData
_atype(x::AbstractCTData) = _atype(typeof(x))

eltype(::Type{T}) where {T<:AbstractCTData} = eltype(_atype(T))

ctimage(gsd::AbstractCTData) = gsd.image
ctsinogram(gsd::AbstractCTData) = gsd.sinog
cttomogram(gsd::AbstractCTData) = gsd.tomog

"""
    mutable struct CTData{M}

Hold data for reconstruction.
`M` should be a subtype of `AbstractArray`, for example
`Matrix{Float32}` or `AFMatrix{Float64}`.
By default, the library will use CPU objects rather than GPU ones.
"""
mutable struct CTData{T<:Number,M<:AbstractMatrix{T}} <: AbstractCTData
    image::Maybe{CTImage{T,M}}
    sinog::Maybe{CTSinogram{T,M}}
    tomog::Maybe{CTTomogram{T,M}}
end

_atype(::Type{CTData{T,M}}) where {T,M} = M

function CTData{T,M}(gsd::AbstractCTData) where {T,M<:AbstractArray{T}}
    CTData{T,M}(
        mmap(CTImage{T,M}, gsd.image),
        mmap(CTSinogram{T,M}, gsd.sinog),
        mmap(CTTomogram{T,M}, gsd.tomog),
    )
end

CTData{T,M}() where {T,M<:AbstractArray{T}} =
    CTData{T,M}(nothing, nothing, nothing)

"""
    CTData(::Type{M}) where {M<:AbstractArray{<:Number}}

Construct an empty `CTData` object.
"""
CTData(::Type{M}) where {M<:AbstractArray{<:Number}} = CTData{eltype(M),M}()

"""
    CTData([::Type{T} = Float32]) where {T<:Number}

Construct an empty `CTData` object with eltype `T`.
"""
CTData(::Type{T} = Float32) where {T<:Number} = CTData(T,Matrix{T})

"""
    CTData(::AbstractCTImage)

Construct a `CTData` object from its constituents.
"""
CTData(img::CTImage) =
    CTData{eltype(img),_atype(img)}(Some(img), nothing, nothing)
CTData(img::CTSinogram) =
    CTData{eltype(img),_atype(img)}(nothing, Some(img), nothing)
CTData(img::CTTomogram) =
    CTData{eltype(img),_atype(img)}(nothing, nothing, Some(img))

function CTData{T,M}(gsd::AbstractCTData, img::CTImage) where {T,M}
    CTData{T,M}(
        Some(img),
        gsd.sinog,
        gsd.tomog,
    )
end

CTData{T,M}(gsd::AbstractCTData, img::CTSinogram) where {T,M} =
    CTData{T,M}(
        gsd.image,
        Some(img),
        gsd.tomog,
    )

CTData{T,M}(gsd::AbstractCTData, img::CTTomogram) where {T,M} =
    CTData{T,M}(
        gsd.image,
        gsd.sinog,
        Some(img),
    )

function CTData{T,M}(
    gsd::AbstractCTData,
    img::AbstractCTImage,
    rest::AbstractCTImage...
) where {T,M}
    CTData{T,M}(CTData{T,M}(gsd, img), rest...)
end

CTData(gsd::AbstractCTData, imgs::AbstractCTImage...) =
    CTData{eltype(gsd),_atype(gsd)}(gsd, imgs...)

function CTData(img::AbstractCTImage, imgs::AbstractCTImage...)
    CTData{eltype(img),_atype(img)}(CTData(img), imgs...)
end

similar(gsd::AbstractCTData, img::AbstractCTImage, rest::AbstractCTImage...) =
    typeof(gsd)(gsd, img, rest...)
similar(gsd::AbstractCTData) = typeof(gsd)()

end # module


text:
using Test, TOML, StructTypes
using WOML

@testset "WOML" begin

@enum EE Aa Bb Cc

@testset "Union utilities" begin
    @test WOML.union_types(Float32) == (Float32, )
    @test Set(WOML.union_types(Union{Int64, Float32})) ==
          Set([ Int64, Float32 ])
    @test Set(WOML.union_types(Union{Float32, Int64, String})) ==
          Set([ Int64, Float32, String])
    @test Set(WOML.union_types(Union{Int, Float64, String, Vector{Float64}})) ==
          Set([ Int, Float64, String, Vector{Float64} ])
    println("#TODO: Adding a type like `Dict{<:AbstractString, <:Integer} breaks this. Why?")

    @test WOML.union_parse_order(Union{Int, Float64}) ==
            WOML.union_parse_order(Union{Float64, Int}) ==
            [ Int, Float64]
    @test WOML.union_parse_order(Union{Float32, Int64, String}) ==
            [ Int64, Float32, String ]
    @test WOML.union_parse_order(Union{Int, String, Float64}) ==
            WOML.union_parse_order(Union{Float64, Int, String}) ==
            WOML.union_parse_order(Union{String, Float64, Int}) ==
            WOML.union_parse_order(Union{Int, Float64, String}) ==
            [ Int, Float64, String ]
    @test WOML.union_parse_order(Union{Int, Float64, String, Vector, Dict}) ==
          [ Int, Float64, String, Dict, Vector ]

    # Test simple overloads of union_try_parse(), i.e. those which don't need a Converter.
    @test WOML.union_try_parse(nothing, nothing, Nothing) === Some(nothing)
    @test WOML.union_try_parse(nothing, 5, Int) === Some(5)
    @test WOML.union_try_parse(nothing, 5.0, Float64) === Some(5.0)
    @test WOML.union_try_parse(nothing, 5.23, Float64) === Some(5.23)
    @test WOML.union_try_parse(nothing, 5.0, Int8) === Some(Int8(5))
    @test WOML.union_try_parse(nothing, 5.5, Float16) === Some(Float16(5.5))
    @test WOML.union_try_parse(nothing, "hi", Symbol) === Some(:hi)
    @test WOML.union_try_parse(nothing, 5, Symbol) === Some(Symbol(5))
    @test WOML.union_try_parse(nothing, 5.23, Int8) === nothing
end

@testset "Converter Write basic data" begin
    c = WOML.Converter{CM_Write}(
        null_string = "T_NULL",
        numbers_can_be_strings = true
    )

    @test c(5) === 5
    @test c(3.12345) === 3.12345

    s = "hello"
    @test c(s) === s
    @test c(:hello) == s

    @test c(Aa::EE) == "Aa"
    @test c(Aa, Int) == Int(Aa)
    @test c(Bb::EE) == "Bb"
    @test c(Bb, Int) == Int(Bb)
    @test c(Cc::EE) == "Cc"
    @test c(Cc, Int) == Int(Cc)

    @test c((4, 5.5, true, "seven", :eight, Bb)) ==
          [ 4, 5.5, true, "seven", "eight", "Bb" ]

    @test c([ 4, 5, 6 ]) == [ 4, 5, 6 ]
    @test c(Dict(3=>"hi", 5=>:hey_there)) ==
          Dict("3"=>"hi", "5"=>"hey_there")

    @test c([ 3, 4, [ 1, 2, 3 ], 5, 6, 7 ]) ==
            [ 3, 4, [ 1, 2, 3 ], 5, 6, 7 ]

    @test c(nothing) == "T_NULL"
    @test c(NaN) == "NaN"
    @test c(NaN32) == "NaN"
    @test c(+Inf) == "+Inf"
    @test c(-Inf) == "-Inf"
    @test c(+Inf16) == "+Inf"
    @test c(-Inf16) == "-Inf"
    @test c([ nothing, Inf32, NaN64, -20 ]) ==
          [ "T_NULL", "+Inf", "NaN", -20 ]
    @test c(Dict(
                nothing => 5,
                5 => nothing,

                +Inf32 => 70,
                70 => +Inf32,

                NaN => "hi",
                "hi" => NaN,

                Bb => Aa,
                Cc => Bb
          )) ==
          Dict(
              "T_NULL" => 5,
              "5" => "T_NULL",

              "+Inf" => 70,
              "70" => "+Inf",

              "NaN" => "hi",
              "hi" => "NaN",

              "Bb" => "Aa",
              "Cc" => "Bb"
          )

    @test c(5, Union{Int, Float64, EE, Vector, Dict{String, Symbol}}) === 5
    @test c(Aa, Union{Int, Float64, String, EE}) == "Aa"
    @test c([ Aa, Bb, Aa, Cc ], Vector{Int}) == map(Int, [ Aa, Bb, Aa, Cc ])

    @test c(Dict(
                7.5 => :abcd,
                :abc => [ "a", "b", "c" ],
                "def" => [
                    [ 1, 2, 3 ],
                    [ 4, 5, 6 ],
                    [ 7, 8, 9 ]
                ],
                "ghi" => [
                    Dict(:j=>1),
                    Dict(:j=>2),
                    Dict(:j=>3)
                ],
                "jkl" => (
                    4.5, 7, "twenty", :seven_again, Aa
                )
            )) ==
            Dict(
                "7.5"=>"abcd",
                "abc" => [ "a", "b", "c" ],
                "def" => [
                    [ 1, 2, 3 ],
                    [ 4, 5, 6 ],
                    [ 7, 8, 9 ]
                ],
                "ghi" => [
                    Dict("j"=>1),
                    Dict("j"=>2),
                    Dict("j"=>3)
                ],
                "jkl" => [
                    4.5, 7, "twenty", "seven_again", "Aa"
                ]
            )
    @test c(Dict(:a=>:b, :c=>:d)) == Dict("a"=>"b", "c"=>"d")
end
@testset "Converter Read basic data" begin
    c = WOML.Converter{CM_Read}(
        null_string = "T_NULL",
        numbers_can_be_strings = true,
        bools_can_be_strings = true
    )

    @test c(5, Int) === 5
    @test c(5, Float32) === Float32(5)
    @test c(true, Bool) === true

    s = "hello"
    @test c(s, String) === s
    @test c(s, AbstractString) === s
    @test c("hello", String) == s
    @test c("hello", AbstractString) == s

    @test c("Aa", EE) == Aa
    @test c("Bb", EE) == Bb
    @test c("Cc", EE) == Cc
    @test c(Int(Aa), EE) == Aa
    @test c(Int(Bb), EE) == Bb
    @test c(Int(Cc), EE) == Cc

    @test c("4.5", Any) === 4.5
    @test c("false", Any) === false

    @test c("f", Any) === "f"
    @test c("f", Bool) === false
    @test c("True", Any) === "True"
    @test c("True", Bool) === true

    @test c([ 4, 5, 6 ], Vector{Int}) == [ 4, 5, 6 ]
    @test c([ 4, 5, 6 ], Vector{Integer}) == [ 4, 5, 6 ]
    @test c([ 4, 5, 6 ], Vector{Real}) == [ 4, 5, 6 ]

    @test c([ "Aa", "Cc", "Bb", "Aa", "Bb" ], Vector{EE}) ==
          [ Aa, Cc, Bb, Aa, Bb ]

    @test c([1, true, 3.5, "four", :five], Tuple{Int, Bool, Float16, Symbol, String}) ==
          (1, true, Float16(3.5), :four, "five")

    @test c(Dict("3"=>"hi", "5"=>"hey_there"), Dict{Int, Symbol}) ==
          Dict(3=>:hi, 5=>:hey_there)
    @test c(Dict(
                "Aa" => [4, 5],
                "Bb" => [2, 4],
                "Cc" => [10, 11]
            ), Dict{EE, NTuple{2, Int}}) ==
           Dict(
               Aa => (4, 5),
               Bb => (2, 4),
               Cc => (10, 11)
           )

    @test c([ 3, 4, [ 1, 2, 3 ], 5, 6, 7 ], Vector{Any}) ==
            [ 3, 4, [ 1, 2, 3 ], 5, 6, 7 ]

    @test c("T_NULL", Nothing) === nothing
    @test c("NaN", Float32) === NaN32
    @test c("-Inf", Float64) === -Inf64
    @test c("+Inf", Float16) === +Inf16
    @test c("Inf", Float32) === +Inf32

    # Test the overloads of union_try_parse() that do use a converter.
    @test WOML.union_try_parse(c, "abc", Int) === nothing
    @test WOML.union_try_parse(c, "abc", Float64) === nothing
    @test WOML.union_try_parse(c, "4.2", Float64) === Some(4.2)
    @test WOML.union_try_parse(c, "-8.0", Int16) === Some(Int16(-8))
    @test WOML.union_try_parse(c, "T_NULL", Nothing) === Some(nothing)
    @test WOML.union_try_parse(c, "abcd", Nothing) === nothing
    @test WOML.union_try_parse(c, "abcd", Nothing) === nothing

    for (actual, expected) in zip(c([ "T_NULL", "+Inf", "NaN", -20 ], Vector{Any}),
                                    Any[ nothing, Inf64, NaN64, -20 ])
        @test actual === expected
    end

    real_d = c(Dict(
        "+Inf" => 70,
        "70" => "+Inf",

        "NaN" => 30,
        30 => "NaN"
    ), Dict{AbstractFloat, AbstractFloat})
    expected_d = Dict(
        +Inf => 70.0,
        70.0 => +Inf,

        NaN => 30.0,
        30.0 => NaN
    )
    @test length(real_d) == length(expected_d)
    for expected_pair in expected_d
        @test haskey(real_d, expected_pair[1])
        @test real_d[expected_pair[1]] === expected_pair[2]
    end

    in_d = Dict(
        "7.5"=>"abcd",
        "true" => [ "a", "b", "c" ],
        "def" => [
            [ 1, 2, 3 ],
            [ 4, 5, 6 ],
            [ 7, 8, 9 ]
        ],
        "false" => [
            Dict("j"=>1),
            Dict("j"=>2),
            Dict("j"=>3)
        ],
        "Bb" => "Aa"
    )
    expected_d = Dict(
        7.5=>"abcd",
        true => [ "a", "b", "c" ],
        "def" => [
            [ 1, 2, 3 ],
            [ 4, 5, 6 ],
            [ 7, 8, 9 ]
        ],
        false => [
            Dict("j"=>1),
            Dict("j"=>2),
            Dict("j"=>3)
        ],
        "Bb" => "Aa"
    )
    @test c(in_d, Dict{Any, Any}) == expected_d

    @test c(Dict("a"=>"b", "c"=>"d"), Dict{Symbol, Symbol}) ==
          Dict(:a=>:b, :c=>:d)
    @test c(Dict("a" => [ 4, 5, 6 ], "b" => [-3, -1, 0 ]),
                 Dict{Symbol, NTuple{3, Int}}) ==
          Dict(:a => (4, 5, 6), :b => (-3, -1, 0))

    @test c(5, Union{Int, String}) === 5
    @test c(5.6, Union{Int, Float64}) === 5.6
    @test c(5.0, Union{Int, Float64}) === 5
    @test c("hey", Union{Symbol, String}) === :hey
    @test c("hey", Union{String, Int}) == "hey"
    @test c("5", Union{String, Int}) === 5
    @test c("-Inf", Union{String, Int, Float32}) === -Inf32
    @test c("abcd", Union{String, Int, EE}) == "abcd"
    @test c("AA", Union{String, EE}) == "AA"
    @test c("Aa", Union{String, EE}) == Aa

    @test c([ 2, 3, 4 ], Union{String, Vector{Int}, Vector{Float64}}) ==
          [ 2, 3, 4 ]
end

@testset "Converter Write struct data" begin
    c = WOML.Converter{CM_Write}(
        default_struct_type = StructTypes.Struct(),
        write_null_fields = false,
        numbers_can_be_strings = true
    )

    struct Sw
        x
    end
    StructTypes.StructType(::Type{Sw}) = StructTypes.Struct()

    @test c(Sw(5)) == Dict("x"=>5)
    @test c(Sw(:a)) == Dict("x"=>"a")
    @test c(Sw("1\"2\"3")) == Dict("x"=>"1\"2\"3")

    @test c(Sw([ 1, 2, 3 ])) == Dict("x" => [ 1, 2, 3 ])
    @test c(Sw(Dict(:y => NaN))) == Dict("x" => Dict("y" => "NaN"))

    @test c(Sw(Sw(-1.5))) == Dict("x" => Dict("x" => -1.5))

    struct SSw
        a::Real
        b::Sw
        c::Union{Nothing, Vector{Int}}
    end
    # Don't explicitly give the struct-type, it should default to `Struct()`.

    @test c(SSw(true, Sw(5), [ 4, 5, 6 ])) ==
          Dict("a" => true,
               "b" => Dict("x" => 5),
               "c" => [ 4, 5, 6 ])
    # Note that the null field will be omitted in writing.
    @test c(SSw(-Inf16, Sw(Symbol("x\"y\"z-3")), nothing)) ==
          Dict("a" => "-Inf",
               "b" => Dict("x" => "x\"y\"z-3"))
end
@testset "Converter Read struct data" begin
    c = WOML.Converter{CM_Read}(
        null_string = "T_NULL",
        default_struct_type = StructTypes.Struct(),
        numbers_can_be_strings = true,
        bools_can_be_strings = true
    )

    struct Sr
        x
    end
    StructTypes.StructType(::Type{Sr}) = StructTypes.Struct()
    Base.:(==)(a::Sr, b::Sr) = isequal(a.x, b.x)  # To handle things like NaN fields

    @test c(Dict("x"=>5), Sr) === Sr(5)
    @test c(Dict("x"=>"a"), Sr) === Sr("a")
    @test c(Dict("x"=>"1\"2\"3"), Sr) === Sr("1\"2\"3")

    @test c(Dict("x" => [ 1, 2, 3 ]), Sr) == Sr([ 1, 2, 3 ])
    @test c(Dict("x" => Dict("y" => "NaN")), Sr) == Sr(Dict("y" => NaN))

    @test c(Dict("x" => -1.5), Sr) === Sr(-1.5)
    @test c(Dict("x" => "-1.5"), Sr) === Sr(-1.5)

    struct SSr
        a::Real
        b::Sr
        c::Vector{Int}
        d::Nothing
    end
    # Don't explicitly give the struct-type, it should default to `Struct()`.
    Base.:(==)(a::SSr, b::SSr) = all(isequal(getfield(a, n), getfield(b, n))
                                       for n in fieldnames(SSr))

    @test c(Dict("a" => true,
                 "b" => Dict("x" => 5),
                 "c" => [ 4, 5, 6 ],
                 "d" => "T_NULL"),
            SSr) ==
          SSr(true, Sr(5), [ 4, 5, 6 ], nothing)
    # Try again, omitting the field that is null.
    @test c(Dict("a" => "-Inf",
                 "b" => Dict("x" => "x\"y\"z-3"),
                 "c" => [ ]),
            SSr) ==
          SSr(-Inf64, Sr("x\"y\"z-3"), Int[ ], nothing)
end

end # Main test-set

println("#TODO: Test enums in struct data")
println("#TODO: Test non-default converter settings")
println("#TODO: Test abstract types")
println("#TODO: Test custom types")
println("#TODO: Add thorough union tests everywhere")

ext: jl


text:
###############################################################################
#
# Calculate SIF fluxes
#
###############################################################################
"""
    SIF_fluxes!(leaves::Array{LeafBios{FT},1},
                can_opt::CanopyOpticals{FT},
                can_rad::CanopyRads{FT},
                can::Canopy4RT{FT},
                soil::SoilOpticals{FT},
                wls::WaveLengths{FT},
                rt_con::RTCache{FT},
                rt_dim::RTDimensions;
                photon::Bool = true
    ) where {FT<:AbstractFloat}

Computes 2-stream diffusive radiation transport for SIF radiation (calls
    [`diffusive_S!`] internally). Layer reflectance and transmission is
    computed from SW optical properties, layer sources from absorbed light and
    SIF efficiencies. Boundary conditions are zero SIF incoming from atmosphere
    or soil.
- `leaves` Array of [`LeafBios`](@ref) type struct
- `can_opt` [`CanopyOpticals`](@ref) type struct
- `can_rad` [`CanopyRads`](@ref) type struct
- `can` [`Canopy4RT`](@ref) type struct
- `soil` [`SoilOpticals`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
- `rt_con` [`RTCache`](@ref) type cache
- `rt_dim` [`RTDimensions`](@ref) type struct
- `photon` If true, use photon unit in the matrix conversion

"""
function SIF_fluxes!(
            leaves::Array{LeafBios{FT},1},
            can_opt::CanopyOpticals{FT},
            can_rad::CanopyRads{FT},
            can::Canopy4RT{FT},
            soil::SoilOpticals{FT},
            wls::WaveLengths{FT},
            rt_con::RTCache{FT},
            rt_dim::RTDimensions;
            photon::Bool = true
) where {FT<:AbstractFloat}
    # unpack variables from structures
    @unpack LAI, lidf, nLayer, Ω = can;
    @unpack a, absfo, absfs, absfsfo, cosΘ_l, cos2Θ_l, fo, fs, fsfo, Po, Ps, Pso, sigb, vb, vf = can_opt;
    @unpack E_down, E_up, ϕ_shade, ϕ_sun = can_rad;
    @unpack ρ_SW_SIF = soil;
    @unpack dWL_iWLE, iWLE, iWLF, WLE, WLF = wls;
    sf_con = rt_con.sf_con;

    # 1. define some useful parameters
    iLAI = LAI * Ω / nLayer;

    # 2. calculate some useful parameters
    sf_con.τ_dd .= 1 .- view(a, iWLF, :) .* iLAI;
    sf_con.ρ_dd .= view(sigb, iWLF, :) .* iLAI;

    # 3. Compute mid layer Ps,Po,Pso
    #    Qso = (Pso[1:nLayer] + Pso[2:nLayer+1]) / 2;
    #    Qs  = ( Ps[1:nLayer] +  Ps[2:nLayer+1]) / 2;
    #    Qo  = ( Po[1:nLayer] +  Po[2:nLayer+1]) / 2;
    Qso = view(Pso, 1:nLayer);
    Qs  = view(Ps , 1:nLayer);
    Qo  = view(Po , 1:nLayer);

    # 4.  reflectance, transmittance factors in a thin layer the following are
    #     vectors with length [nl,nWL]
    sf_con.sun_dwl_iWlE .= view(can_opt.Es_, iWLE, 1) .* dWL_iWLE;
    if photon sf_con.sun_dwl_iWlE .*= WLE .* _FAC(FT) end;
    @inbounds for i=1:nLayer
        if length(leaves)>1
            Mb = leaves[i].Mb;
            Mf = leaves[i].Mf;
        else
            Mb = leaves[1].Mb;
            Mf = leaves[1].Mf;
        end
        sf_con.M⁺ .= (Mb .+ Mf) ./ 2;
        sf_con.M⁻ .= (Mb .- Mf) ./ 2;

        # Need to normalize incoming radiation bin
        # change from mSCOPE to enable different WL grids!
        mul!(sf_con.M⁻_sun, sf_con.M⁻, sf_con.sun_dwl_iWlE);
        mul!(sf_con.M⁺_sun, sf_con.M⁺, sf_con.sun_dwl_iWlE);

        sf_con.tmp_dwl_iWlE .= view(E_down, iWLE, i) .* dWL_iWLE;
        if photon sf_con.tmp_dwl_iWlE .*= WLE .* _FAC(FT) end;
        mul!(sf_con.M⁺⁻, sf_con.M⁺, sf_con.tmp_dwl_iWlE);

        sf_con.tmp_dwl_iWlE .= view(E_up, iWLE, i+1) .* dWL_iWLE;
        if photon sf_con.tmp_dwl_iWlE .*= WLE .* _FAC(FT) end;
        mul!(sf_con.M⁺⁺, sf_con.M⁺, sf_con.tmp_dwl_iWlE);

        sf_con.tmp_dwl_iWlE .= view(E_up, iWLE, i+1) .* dWL_iWLE;
        if photon sf_con.tmp_dwl_iWlE .*= WLE .* _FAC(FT) end;
        mul!(sf_con.M⁻⁺, sf_con.M⁻, sf_con.tmp_dwl_iWlE);

        sf_con.tmp_dwl_iWlE .= view(E_down, iWLE, i) .* dWL_iWLE;
        if photon sf_con.tmp_dwl_iWlE .*= WLE .* _FAC(FT) end;
        mul!(sf_con.M⁻⁻, sf_con.M⁻, sf_con.tmp_dwl_iWlE);

        # Here comes the tedious part:
        # TODO move them to a seprate function
        sf_con.ϕ_cosΘ .= view(ϕ_sun, :, :, i) .* cosΘ_l;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        sunCos = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.ϕ_cosΘ .= ϕ_shade[i] .* cosΘ_l;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        shadeCos = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.ϕ_cosΘ .= view(ϕ_sun, :, :, i) .* cos2Θ_l;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        sunCos2 = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.ϕ_cosΘ .= ϕ_shade[i] .* cos2Θ_l;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        shadeCos2 = mean(sf_con.ϕ_cosΘ_lidf);

        mul!(sf_con.ϕ_cosΘ_lidf, view(ϕ_sun, :, :, i)', lidf);
        sunLidf = mean(sf_con.ϕ_cosΘ_lidf);
        shadeLidf = mean(lidf) * ϕ_shade[i];

        sf_con.ϕ_cosΘ .= view(ϕ_sun, :, :, i) .* absfsfo;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        _mean_absfsfo = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.ϕ_cosΘ .= view(ϕ_sun, :, :, i) .* fsfo;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        _mean_fsfo = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.wfEs .= _mean_absfsfo .* sf_con.M⁺_sun .+
                       _mean_fsfo    .* sf_con.M⁻_sun;

        sf_con.ϕ_cosΘ .= view(ϕ_sun, :, :, i) .* absfs;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        _mean_absfs = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.ϕ_cosΘ .= view(ϕ_sun, :, :, i) .* fs .* cosΘ_l;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        _mean_fs = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.sfEs .= _mean_absfs .* sf_con.M⁺_sun .- _mean_fs .* sf_con.M⁻_sun;
        sf_con.sbEs .= _mean_absfs .* sf_con.M⁺_sun .+ _mean_fs .* sf_con.M⁻_sun;

        sf_con.ϕ_cosΘ .= ϕ_shade[i] .* absfo;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        _mean_absfo = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.ϕ_cosΘ .= ϕ_shade[i] .* fo .* cosΘ_l;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        _mean_fo = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.vfEplu_shade .= _mean_absfo .* sf_con.M⁺⁺ .- _mean_fo .* sf_con.M⁻⁺;
        sf_con.vbEmin_shade .= _mean_absfo .* sf_con.M⁺⁻ .+ _mean_fo .* sf_con.M⁻⁻;

        sf_con.ϕ_cosΘ .= view(ϕ_sun, :, :, i) .* absfo;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        _mean_absfo = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.ϕ_cosΘ .= view(ϕ_sun, :, :, i) .* fo .* cosΘ_l;
        mul!(sf_con.ϕ_cosΘ_lidf, adjoint(sf_con.ϕ_cosΘ), lidf);
        _mean_fo = mean(sf_con.ϕ_cosΘ_lidf);

        sf_con.vfEplu_sun .= _mean_absfo .* sf_con.M⁺⁺ .- _mean_fo .* sf_con.M⁻⁺;
        sf_con.vbEmin_sun .= _mean_absfo .* sf_con.M⁺⁻ .+ _mean_fo .* sf_con.M⁻⁻;

        sf_con.sigfEmin_shade .= shadeLidf .* sf_con.M⁺⁻ .- shadeCos2 .* sf_con.M⁻⁻;
        sf_con.sigbEmin_shade .= shadeLidf .* sf_con.M⁺⁻ .+ shadeCos2 .* sf_con.M⁻⁻;
        sf_con.sigfEmin_sun   .= sunLidf   .* sf_con.M⁺⁻ .- sunCos2   .* sf_con.M⁻⁻;
        sf_con.sigbEmin_sun   .= sunLidf   .* sf_con.M⁺⁻ .+ sunCos2   .* sf_con.M⁻⁻;
        sf_con.sigfEplu_shade .= shadeLidf .* sf_con.M⁺⁺ .- shadeCos2 .* sf_con.M⁻⁺;
        sf_con.sigbEplu_shade .= shadeLidf .* sf_con.M⁺⁺ .+ shadeCos2 .* sf_con.M⁻⁺;
        sf_con.sigfEplu_sun   .= sunLidf   .* sf_con.M⁺⁺ .- sunCos2   .* sf_con.M⁻⁺;
        sf_con.sigbEplu_sun   .= sunLidf   .* sf_con.M⁺⁺ .+ sunCos2   .* sf_con.M⁻⁺;

        # Fluxes:
        #    sunlit for each layer
        #    shade leaf for each layer
        #    Eq. 29a for sunlit leaf
        #    Eq. 29b for sunlit leaf
        #    Eq. 29a for shade leaf
        #    Eq. 29b for shade leaf
        sf_con.piLs[:,i]  .= sf_con.wfEs .+ sf_con.vfEplu_sun .+ sf_con.vbEmin_sun;
        sf_con.piLd[:,i]  .= sf_con.vbEmin_shade .+ sf_con.vfEplu_shade;
        sf_con.Fsmin[:,i] .= sf_con.sfEs .+ sf_con.sigfEmin_sun .+ sf_con.sigbEplu_sun;
        sf_con.Fsplu[:,i] .= sf_con.sbEs .+ sf_con.sigbEmin_sun .+ sf_con.sigfEplu_sun;
        sf_con.Fdmin[:,i] .= sf_con.sigfEmin_shade .+ sf_con.sigbEplu_shade;
        sf_con.Fdplu[:,i] .= sf_con.sigbEmin_shade .+ sf_con.sigfEplu_shade;

        # Total weighted fluxes
        _qs_iLAI   = Qs[i] * iLAI;
        _1_qs_iLAI = (1 - Qs[i]) * iLAI;
        sf_con.S⁻[:,i]   .= _qs_iLAI .* view(sf_con.Fsmin, :, i) .+ _1_qs_iLAI .* view(sf_con.Fdmin, :, i);
        sf_con.S⁺[:,i]   .= _qs_iLAI .* view(sf_con.Fsplu, :, i) .+ _1_qs_iLAI .* view(sf_con.Fdplu, :, i);
        sf_con.Femo[:,i] .= _qs_iLAI .* view(sf_con.piLs , :, i) .+ _1_qs_iLAI .* view(sf_con.piLd , :, i);
    end

    # 5. Compute diffusive fluxes within canopy
    #    Use Zero SIF fluxes as top and bottom boundary:
    # TODO pre-allocate these!
    diffusive_S!(sf_con, soil, rt_dim);

    # 6. Save in output structures!
    #    direct Sunlit leaves
    #    direct shaded leaves
    _iLAI_pi = iLAI / FT(pi);

    # why this step is so slow?
    sf_con.tmp_1d_nLayer .= Qso .* _iLAI_pi;
    mul!(can_rad.SIF_obs_sunlit, sf_con.piLs, sf_con.tmp_1d_nLayer);
    sf_con.tmp_1d_nLayer .= (Qo .- Qso) .* _iLAI_pi;
    mul!(can_rad.SIF_obs_shaded, sf_con.piLd, sf_con.tmp_1d_nLayer);

    # 7. SIF from scattered internally and soil contribution
    sf_con.tmp_2d_nWlF_nLayer .= view(vb, iWLF, :) .* view(sf_con.F⁻, :, 1:nLayer) .+ view(vf, iWLF, :) .* view(sf_con.F⁺, :, 1:nLayer);
    mul!(can_rad.SIF_obs_scattered, sf_con.tmp_2d_nWlF_nLayer, Qo);
    can_rad.SIF_obs_scattered .= can_rad.SIF_obs_scattered .* _iLAI_pi;
    can_rad.SIF_obs_soil .= ( ρ_SW_SIF .* view(sf_con.F⁻, :, nLayer+1) ) .* Po[end] ./ FT(pi);

    can_rad.SIF_hemi .= view(sf_con.F⁺, :, 1);
    can_rad.SIF_obs  .= can_rad.SIF_obs_sunlit .+ can_rad.SIF_obs_shaded .+ can_rad.SIF_obs_scattered .+ can_rad.SIF_obs_soil;
    # can_rad.SIF_sum[:]  = sum(sf_con.S⁻ + sf_con.S⁺, dims=2);
    @inbounds for j in eachindex(can_rad.SIF_sum)
        can_rad.SIF_sum[j] = sum( view(sf_con.S⁻, j, :) ) + sum( view(sf_con.S⁺, j, :) );
    end

    if photon
        can_rad.SIF_obs ./= WLF .* _FAC(FT);
    end;

    return nothing
end

"""
    SIF_fluxes!(leaf::LeafBios{FT},
                in_rad::IncomingRadiation{FT},
                wls::WaveLengths{FT},
                rt_con::RTCache{FT},
                fqe::FT = FT(0.01);
                photon::Bool = true
    ) where {FT<:AbstractFloat}

Leaf level SIF, given
- `leaf` [`LeafBios`](@ref) type struct
- `in_rad` [`IncomingRadiation`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
- `rt_con` [`RTCache`](@ref) type cache
- `fqe` Fluorescence quantum yield (default at 1%)
- `photon` If true, use photon unit in the matrix conversion

Note that `in_rad` assumes direct light with zenith angle of 0, and a zenith
    angle correction needs to be made before passing it to this function. The
    up- and down-ward SIF are stored in `sf_con` as `M⁻_sun` and `M⁺_sun`.
"""
function SIF_fluxes!(
            leaf::LeafBios{FT},
            in_rad::IncomingRadiation{FT},
            wls::WaveLengths{FT},
            rt_con::RTCache{FT},
            fqe::FT = FT(0.01);
            photon::Bool = true
) where {FT<:AbstractFloat}
    # unpack the values
    @unpack Mb, Mf = leaf;
    @unpack dWL_iWLE, iWLE, WLE, WLF = wls;
    sf_con = rt_con.sf_con;
    sf_con.tmp_dwl_iWlE  .= (view(in_rad.E_direct , iWLE, 1) .+ view(in_rad.E_diffuse, iWLE, 1)) .* dWL_iWLE;
    if photon
        sf_con.tmp_dwl_iWlE .*= WLE .* _FAC(FT);
    end;

    # calculate the SIF spectra for direct light
    # sf_con.M⁺ .= (Mb .+ Mf) ./ 2;
    # sf_con.M⁻ .= (Mb .- Mf) ./ 2;
    # mul!(sf_con.M⁺_sun, sf_con.M⁺, sf_con.tmp_dwl_iWlE);
    # mul!(sf_con.M⁻_sun, sf_con.M⁻, sf_con.tmp_dwl_iWlE);
    mul!(sf_con.M⁺_sun, Mb, sf_con.tmp_dwl_iWlE);
    mul!(sf_con.M⁻_sun, Mf, sf_con.tmp_dwl_iWlE);
    if photon
        sf_con.M⁺_sun ./= WLF .* _FAC(FT);
        sf_con.M⁻_sun ./= WLF .* _FAC(FT);
    end;

    # divide by pi to account for scattering
    sf_con.M⁻_sun .*= fqe / pi;
    sf_con.M⁺_sun .*= fqe / pi;

    return nothing
end

max_stars_repo_name: climate-machine/Land

max_issues_repo_name: climate-machine/Land

max_forks_repo_name: climate-machine/Land


text:
module StipplePlotly

using Genie, Stipple, Stipple.Reexport

#===#

const assets_config = Genie.Assets.AssetsConfig(package = "StipplePlotly.jl")

#===#

function deps() :: String
  if ! Genie.Assets.external_assets(Stipple.assets_config)

    Genie.Router.route(Genie.Assets.asset_path(assets_config, :js, file="plotly.min")) do
      Genie.Renderer.WebRenderable(
        Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), file="plotly.min.js")),
        :javascript) |> Genie.Renderer.respond
    end

    Genie.Router.route(Genie.Assets.asset_path(assets_config, :js, file="resizesensor.min")) do
      Genie.Renderer.WebRenderable(
        Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), file="resizesensor.min.js")),
        :javascript) |> Genie.Renderer.respond
    end

    Genie.Router.route(Genie.Assets.asset_path(assets_config, :js, file="lodash.min")) do
      Genie.Renderer.WebRenderable(
        Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), file="lodash.min.js")),
        :javascript) |> Genie.Renderer.respond
    end

    Genie.Router.route(Genie.Assets.asset_path(assets_config, :js, file="vueresize.min")) do
      Genie.Renderer.WebRenderable(
        Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), file="vueresize.min.js")),
        :javascript) |> Genie.Renderer.respond
    end

    Genie.Router.route(Genie.Assets.asset_path(assets_config, :js, file="vueplotly.min")) do
      Genie.Renderer.WebRenderable(
        Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..")), file="vueplotly.min.js")),
        :javascript) |> Genie.Renderer.respond
    end

  end

  string(
    Genie.Renderer.Html.script(src=Genie.Assets.asset_path(assets_config, :js, file="plotly.min")),
    Genie.Renderer.Html.script(src=Genie.Assets.asset_path(assets_config, :js, file="resizesensor.min")),
    Genie.Renderer.Html.script(src=Genie.Assets.asset_path(assets_config, :js, file="lodash.min")),
    Genie.Renderer.Html.script(src=Genie.Assets.asset_path(assets_config, :js, file="vueresize.min")),
    Genie.Renderer.Html.script(src=Genie.Assets.asset_path(assets_config, :js, file="vueplotly.min"))
  )
end

#===#

include("Charts.jl")
@reexport using .Charts

function __init__()
  push!(Stipple.DEPS, deps)
end

end # module


text:
using Distributions, Random, Plots; pyplot()
Random.seed!(0)

beta, delta, gamma = 0.25, 0.4, 0.1
initialInfect = 0.025
M = 1000
I0 = Int(floor(initialInfect*M))
N = 30

function simulateSIRDoobGillespie(beta,delta,gamma,I0,M,T)
    t, S, E, I, R = 0.0, M-I0, 0, I0, 0
    tValues, sValues, eValues, iValues, rValues = [0.0], [S], [E], [I], [R]
    while t<T
        infectionRate = beta*I*S
        symptomRate = delta*E
        removalRate = gamma*I
        totalRate = infectionRate + symptomRate + removalRate
        probs = [infectionRate, symptomRate, removalRate]/totalRate
        t += rand(Exponential(1/(totalRate)))
        u = rand()
        if u < probs[1]
            S -= 1; E += 1
        elseif u < probs[1] + probs[2]
            E -=1; I+=1
        else
            I -= 1; R += 1
        end
        push!(tValues,t)
        push!(sValues,S);push!(eValues,E);push!(iValues,I);push!(rValues,R)
        I == 0 && break
    end
    return [tValues, sValues, eValues, iValues, rValues]
end

tV,sV,eV,iV,rV = simulateSIRDoobGillespie(beta/M,delta,gamma,I0,M,Inf)
lastT = tV[end]

finals = [simulateSIRDoobGillespie(beta/M,delta,gamma,I0,M,Inf)[5][end]
                for _ in 1:N]/M

p1 = plot(tV,sV/M,label = "Susceptible", c=:green)
plot!(tV,eV/M,label = "Exposed", c=:blue)
plot!(tV,iV/M,label = "Infected",c=:red)
plot!(tV,rV/M,label = "Removed", c=:yellow,
    xlabel = "Time", ylabel = "Proportion",
    legend = :topleft, xlim = (0,lastT*1.05))
scatter!(lastT*1.025*ones(N),finals, c = :yellow,label= "Final Infected")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@connector function Pin(;name)
    @variables v(t) i(t)
    ODESystem(Equation[], t, [v, i], [], name=name, defaults=Dict(v=>1.0, i=>1.0))
end

@connector function DigitalPin(; name)
    @variables val(t) v(t) i(t)
    eqs = [
        val ~ IfElse.ifelse((0.0 <= v) & (v <= 0.8) | (2.0 <= v) & (v <= 5.0),
                                IfElse.ifelse(v > 2.0, 1, 0), X)
    ]
    ODESystem(Equation[], t, [val, v, i], [], defaults=Dict(val=>0, i=>0), name=name)
end

abstract type ElectricalPin end
ModelingToolkit.promote_connect_rule(::Type{DigitalPin}, ::Type{Pin}) = ElectricalPin
ModelingToolkit.promote_connect_rule(::Type{Pin}, ::Type{DigitalPin}) = ElectricalPin
ModelingToolkit.promote_connect_rule(::Type{ElectricalPin}, ::Type{DigitalPin}) = ElectricalPin
ModelingToolkit.promote_connect_rule(::Type{ElectricalPin}, ::Type{Pin}) = ElectricalPin

function ModelingToolkit.connect(::Type{<:Pin}, ps...)
    eqs = [
           0 ~ sum(p->p.i, ps) # KCL
          ]
    # KVL
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end

    return eqs
end

function ModelingToolkit.connect(::Type{DigitalPin}, ps...)
    eqs = [
           0 ~ sum(p->p.i, ps) # KCL
          ]
    # KVL
    for i in 1:length(ps)-1
        push!(eqs, ps[i].val ~ ps[i+1].val)
    end
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end
    return eqs
end

function ModelingToolkit.connect(::Type{ElectricalPin}, ps...)
    eqs = [
           0 ~ sum(p->p.i, ps) # KCL
          ]

    # KVL
    digpins = ModelingToolkit.ODESystem[]
    for p in ps
        ModelingToolkit.get_connection_type(p) == DigitalPin && push!(digpins, p)
    end
    for i in 1:length(digpins)-1
        push!(eqs, digpins[i].val ~ digpins[i+1].val)
    end
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end
    return eqs
end


text:
# (1) 利用するデータを持つライブラリのインストール（初回のみ）
## ] + addで
# using RCall
# @rlibrary devtools
# R"devtools::install_github('itamarcaspi/experimentdatar')"

# (2) ライブラリとデータの読み込み
using DataFrames
using DataFramesMeta
using GLM
using RCall
@rimport tidyr
@rlibrary experimentdatar
vouchers = rcopy(R"experimentdatar::vouchers")

## 可視化用。@rlibraryでは都度ライブラリの指定が必要なのでR"library()"でインポートする
R"library(ggplot2)"
R"library(magrittr)"

## broom::tidyと同じ働きをする関数を自作
function broom_tidy_julia(model)
  table = coeftable(model)
  return DataFrame(
    term = coefnames(model),
    estimate = table.cols[1],
    stderror = table.cols[2],
    statistic = table.cols[3],
    pvalue = table.cols[4],
  )
end

# (3) Angrist(2002)のTable 3. bogota 1995の再現
## 回帰式の準備
### 回帰式で使う文字列の準備
formula_x_base = :VOUCH0
formula_x_covariates = [
  :SVY,
  :HSVISIT,
  :AGE,
  :STRATA1,
  :STRATA2,
  :STRATA3,
  :STRATA4,
  :STRATA5,
  :STRATA6,
  :STRATAMS,
  :D1993,
  :D1995,
  :D1997,
  :DMONTH1,
  :DMONTH2,
  :DMONTH3,
  :DMONTH4,
  :DMONTH5,
  :DMONTH6,
  :DMONTH7,
  :DMONTH8,
  :DMONTH9,
  :DMONTH10,
  :DMONTH11,
  :DMONTH12,
  :SEX2,
]
formula_y = [
  :TOTSCYRS,
  :INSCHL,
  :PRSCH_C,
  :USNGSCH,
  :PRSCHA_1,
  :FINISH6,
  :FINISH7,
  :FINISH8,
  :REPT6,
  :REPT,
  :NREPT,
  :MARRIED,
  :HASCHILD,
  :HOURSUM,
  :WORKING3,
]

### formula_yの各要素に対して共変量を含まない回帰式を作成
base_reg_formula =
  [FormulaTerm(Term(y), Term(formula_x_base)) for y in formula_y]
base_reg_formula_names = String.(formula_y) .* "_base"

### formula_yの各要素に対して共変量を含む回帰式を作成
covariate_reg_formula = [
  FormulaTerm(
    Term(y),
    Tuple(Term.(vcat(formula_x_base, formula_x_covariates))),
  ) for y in formula_y
]
covariate_reg_formula_names = String.(formula_y) .* "_covariate"

### モデル式のベクトルを作成
### 直接データフレーム化するのでスキップ

### モデル式のベクトルをデータフレーム化する
models = DataFrame(
  model_index = vcat(base_reg_formula_names, covariate_reg_formula_names),
  formula = vcat(base_reg_formula, covariate_reg_formula),
)

## 回帰分析を実行
### bogota 1995のデータを抽出する
@linq regression_data = vouchers |> where(:TAB3SMPL .== 1, :BOG95SMP .== 1)

### まとめて回帰分析を実行
@linq df_models =
  models |> transform(
    model = map(formula -> lm(formula, regression_data, true), :formula),
  ) |> transform(lm_result = map(broom_tidy_julia, :model))

### モデルの結果を整形
@linq df_results =
  df_models |> select(:formula, :model_index, :lm_result) |>
  tidyr.unnest(cols = :lm_result) |> rcopy()

# ------------------
## 通学率と奨学金の利用傾向の可視化(ch2_plot2.png)
### PRSCHA_1, USNGSCHに対するVOUCH0の効果を取り出す
@linq using_voucher_results =
  df_results |> where(
    :term .== "VOUCH0",
    match.(r"PRSCHA_1|USNGSCH", :model_index) .!= nothing,
  ) |> select(:model_index, :term, :estimate, :stderror, :pvalue) |>
  orderby(:model_index)

### 取り出した効果をggplotで可視化
### RCallでggplot2にObjectを渡してやればいい
R"$(using_voucher_results) %>%
  ggplot(aes(y = estimate, x = model_index)) +
  geom_point() +
  geom_errorbar(aes(ymax = estimate + stderror*1.96,
                    ymin = estimate - stderror*1.96,
                    width = 0.1)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5),
        legend.position = 'bottom',
        plot.margin = margin(0.5,1,0.5,1, 'cm'))"

## 留年の傾向を可視化
### PRSCH_C,INSCHL,FINISH6-8,REPTに対するVOUCH0の効果を取り出す
@linq going_private_results =
  df_results |> where(
    :term .== "VOUCH0",
    match.(r"PRSCH_C|INSCHL|FINISH|REPT", :model_index) .!= nothing,
  ) |> select(:model_index, :term, :estimate, :stderror, :pvalue) |>
  orderby(:model_index)

### 取り出した効果をggplotで可視化
going_private_results_gg =
  @where(going_private_results, match.(r"covariate", :model_index) .!= nothing)

R"$(going_private_results_gg) %>%
  ggplot(aes(y = estimate, x = model_index)) +
  geom_point() +
  geom_errorbar(aes(ymax = estimate + stderror*1.96,
                    ymin = estimate - stderror*1.96,
                    width = 0.1)) +
  geom_hline(yintercept = 0, linetype = 2) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5),
        legend.position = 'bottom',
        plot.margin = margin(0.5,1,0.5,1, 'cm'))"

# (4) Angrist(2002)のTable.4 & 6 bogota 1995の再現
## table4に使うデータを抜き出す
@linq data_tbl4_bog95 =
  vouchers |> where(
    :BOG95SMP .== 1,
    :TAB3SMPL .== 1,
    :SCYFNSH .!== missing,
    :FINISH6 .!== missing,
    :PRSCHA_1 .!== missing,
    :REPT6 .!== missing,
    :NREPT .!== missing,
    :INSCHL .!== missing,
    :FINISH7 .!== missing,
    :PRSCH_C .!== missing,
    :FINISH8 .!== missing,
    :PRSCHA_2 .!== missing,
    :TOTSCYRS .!== missing,
    :REPT .!== missing,
  ) |> select(
    :VOUCH0,
    :SVY,
    :HSVISIT,
    :DJAMUNDI,
    :PHONE,
    :AGE,
    :STRATA1,
    :STRATA2,
    :STRATA3,
    :STRATA4,
    :STRATA5,
    :STRATA6,
    :STRATAMS,
    :DBOGOTA,
    :D1993,
    :D1995,
    :D1997,
    :DMONTH1,
    :DMONTH2,
    :DMONTH3,
    :DMONTH4,
    :DMONTH5,
    :DMONTH6,
    :DMONTH7,
    :DMONTH8,
    :DMONTH9,
    :DMONTH10,
    :DMONTH11,
    :DMONTH12,
    :SEX_MISS,
    :FINISH6,
    :FINISH7,
    :FINISH8,
    :REPT6,
    :REPT,
    :NREPT,
    :SEX2,
    :TOTSCYRS,
    :MARRIED,
    :HASCHILD,
    :HOURSUM,
    :WORKING3,
    :INSCHL,
    :PRSCH_C,
    :USNGSCH,
    :PRSCHA_1,
  )

## 女子生徒のみのデータでの回帰分析
### 女子生徒のデータだけ取り出す
@linq regression_data = data_tbl4_bog95 |> where(:SEX2 .== 0)

### まとめて回帰分析を実行
@linq df_models =
  models |> transform(
    model = map(formula -> lm(formula, regression_data, true), :formula),
  ) |> transform(lm_result = map(broom_tidy_julia, :model))

### モデルの結果を整形
@linq df_results_female =
  df_models |> transform(gender = fill("female", length(:model_index))) |>
  select(:formula, :model_index, :lm_result, :gender) |>
  tidyr.unnest(cols = :lm_result) |> rcopy()

## 男子生徒のみのデータでの回帰分析
@linq regression_data = data_tbl4_bog95 |> where(:SEX2 .== 1)

### まとめて回帰分析を実行
@linq df_models =
  models |> transform(
    model = map(formula -> lm(formula, regression_data, true), :formula),
  ) |> transform(lm_result = map(broom_tidy_julia, :model))

### モデルの結果を整形
@linq df_results_male =
  df_models |> transform(gender = fill("male", length(:model_index))) |>
  select(:formula, :model_index, :lm_result, :gender) |>
  tidyr.unnest(cols = :lm_result) |> rcopy()

## 通学傾向への分析結果の可視化(ch2_plot3.png)
### PRSCHA_1,USNGSCHに対する分析結果を抜き出す
@linq using_voucher_results_gender =
  vcat(df_results_male, df_results_female) |> where(
    :term .== "VOUCH0",
    match.(r"PRSCHA_1|USNGSCH", :model_index) .!= nothing,
  ) |> select(:gender, :model_index, :term, :estimate, :stderror, :pvalue) |>
  orderby(:gender, :model_index) |>
  where(match.(r"covariate", :model_index) .!= nothing)

### ggplotによる可視化
R"$(using_voucher_results_gender) %>%
  ggplot(aes(y = estimate, x = model_index)) +
  geom_point() +
  geom_errorbar(aes(ymax = estimate + stderror*1.96,
                    ymin = estimate - stderror*1.96,
                    width = 0.1)) +
  geom_hline(yintercept = 0, linetype = 2) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5),
        legend.position ='bottom',
        plot.margin = margin(0.5,1,0.5,1, 'cm')) +
  facet_grid(gender ~ .)"

## 留年と通学年数への分析結果の可視化(ch2_plot4.png)
### PRSCH_C,INSCHL,REPT,TOTSCYRS,FINISHに対する分析結果を抜き出す
@linq going_private_results_gender =
  vcat(df_results_male, df_results_female) |> where(
    :term .== "VOUCH0",
    match.(r"PRSCH_C|INSCHL|REPT|TOTSCYRS|FINISH", :model_index) .!= nothing,
  ) |> select(:gender, :model_index, :term, :estimate, :stderror, :pvalue) |>
  orderby(:model_index) |> where(match.(r"covariate", :model_index) .!= nothing)

### ggplotによる可視化
R"$(going_private_results_gender) %>%
  ggplot(aes(y = estimate, x = model_index)) +
  geom_point() +
  geom_errorbar(aes(ymax = estimate + stderror*1.96,
                    ymin = estimate - stderror*1.96,
                    width = 0.1)) +
  geom_hline(yintercept = 0, linetype = 2) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5),
        legend.position = 'bottom',
        plot.margin = margin(0.5,1,0.5,1, 'cm')) +
  facet_grid(gender ~ .)"

## 労働時間に対する分析結果の可視化(ch2_plot5.png)
### HOURに対する分析結果を抜き出す
@linq working_hour_results_gender =
  vcat(df_results_male, df_results_female) |>
  where(:term .== "VOUCH0", match.(r"HOUR", :model_index) .!= nothing) |>
  select(:gender, :model_index, :term, :estimate, :stderror, :pvalue) |>
  orderby(:gender, :model_index) |>
  where(match.(r"covariate", :model_index) .!= nothing)

### ggplotによる可視化
R"$(working_hour_results_gender) %>%
  ggplot(aes(y = estimate, x = model_index)) +
  geom_point() +
  geom_errorbar(aes(ymax = estimate + stderror*1.96,
                    ymin = estimate - stderror*1.96,
                    width = 0.1)) +
  geom_hline(yintercept = 0, linetype = 2) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5),
        legend.position = 'bottom',
        plot.margin = margin(0.5,1,0.5,1, 'cm')) +
  facet_grid(. ~ gender)"


text:
using Random
using Colors
using DataFrames
using DataFramesMeta
using NearestNeighbors
using Statistics

import Base64
import MultivariateStats
import Plots
import AbstractPlotting
import CairoMakie
import Base.show
MK = CairoMakie

plot_molecules!(args...; kwargs...) = plot_molecules(args...; append=true, kwargs...)

function plot_molecules(df_spatial::DataFrame, polygons::Array{Matrix{Float64}, 1}=Matrix{Float64}[]; markersize=2,
        color::Union{Vector, Symbol, String}=:gene, size=(800, 800), poly_strokewidth=1, xlims=nothing, ylims=nothing, offset=(0, 0),
        is_noise::Union{Vector, BitArray, Symbol, Nothing}=nothing, annotation::Union{<:AbstractVector, Nothing} = nothing,
        ann_colors::Union{Nothing, Dict} = nothing, legend=(annotation !== nothing), fontsize=8,
        noise_ann = nothing, shuffle_colors::Bool=false, append::Bool=false,
        polygon_kwargs::KWArgT=nothing, axis_kwargs::KWArgT=nothing, noise_kwargs::KWArgT=nothing, legend_kwargs::KWArgT=nothing, kwargs...)

    noise_args_default = (marker=:xcross, markersize=(0.75 * markersize), strokewidth=0, color="black");
    axis_args_default = (xticklabelsize=12, yticklabelsize=12, :xticksvisible => false, :xticklabelsvisible => false, :yticksvisible => false, :yticklabelsvisible => false);
    legend_args_default = (bgcolor=Colors.RGBA(1, 1, 1, 0.85),);
    polygon_args_default = (strokecolor="black", color="transparent", strokewidth=poly_strokewidth, label="")

    noise_kwargs = update_args(update_args(noise_args_default, Dict(kwargs...)), noise_kwargs)
    axis_kwargs = update_args(axis_args_default, axis_kwargs)
    legend_kwargs = update_args(legend_args_default, legend_kwargs)
    polygon_kwargs = update_args(polygon_args_default, polygon_kwargs)

    if typeof(color) === Symbol
        color = df_spatial[!,color]
    end

    if annotation !== nothing
        annotation = ["$a" for a in annotation]
        if noise_ann !== nothing
            noise_ann = "$noise_ann"
        end
    end

    xlims = something(xlims, val_range(df_spatial.x))
    ylims = something(ylims, val_range(df_spatial.y))

    if !append
        fig = MK.Figure(resolution=size)
        fig[1, 1] = MK.Axis(fig; axis_kwargs...);
    end

    df_noise = nothing

    if typeof(is_noise) === Symbol
        is_noise = df_spatial[!,is_noise]
    end

    if is_noise !== nothing
        df_noise = df_spatial[is_noise,:]
        df_spatial = df_spatial[.!is_noise,:]
        color = color[.!is_noise]
    end

    if is_noise !== nothing
        MK.scatter!(df_noise.x .+ offset[1], df_noise.y .+ offset[2]; noise_kwargs...)
    end

    if annotation === nothing
        MK.scatter!(df_spatial.x .+ offset[1], df_spatial.y .+ offset[2]; color=color,
            strokewidth=0, markersize=markersize, kwargs...)
    else
        ann_vals = annotation[annotation .!= noise_ann] |> unique |> sort
        c_map = Colors.distinguishable_colors(length(ann_vals), Colors.colorant"#007a10", lchoices=range(20, stop=70, length=15))
        if shuffle_colors
            Random.shuffle!(c_map)
        end

        for (color, ann) in zip(c_map, ann_vals)
            style_dict = (ann_colors === nothing) ? Dict() : Dict(:color => ann_colors[ann])
            MK.scatter!(df_spatial.x[annotation .== ann] .+ offset[1], df_spatial.y[annotation .== ann] .+ offset[2];
                strokewidth=0, markersize=markersize, label=ann, color=color, style_dict..., kwargs...)
        end

        if noise_ann in annotation
            MK.scatter!(df_spatial.x[annotation .== noise_ann] .+ offset[1], df_spatial.y[annotation .== noise_ann] .+ offset[2];
                label=noise_ann, noise_kwargs...)
        end

        if legend
            MK.axislegend(;legend_kwargs...)
        end
    end

    if length(polygons) > 0
        MK.poly!([MK.Point2.(eachrow(p .+ [offset[1] offset[2]])) for p in polygons]; polygon_kwargs...)
    end

    MK.xlims!(MK.current_axis(), xlims .+ offset[1])
    MK.ylims!(MK.current_axis(), ylims .+ offset[2])

    return MK.current_figure()
end

function shuffle_labels(labels::Array{Int})
    new_labs = deepcopy(labels)
    mask = (new_labs .!= 0)
    new_labs[mask] = shuffle(1:maximum(labels))[new_labs[mask]]
    return new_labs
end

function shuffle_colors(colors::Vector)
    uniq_cols = unique(colors);
    col_ord = Dict(Pair.(uniq_cols, shuffle(1:length(uniq_cols))));
    return [uniq_cols[col_ord[tc]] for tc in colors]
end

function plot_expression_vectors(vecs...; gene_names::Vector{String}, min_expr_frac::Float64=0.05, alpha::Float64=0.5, fontsize::Int=5, text_offset::Float64=0.005,
        labels::Vector{String}=["y$i" for i in 1:length(vecs)], xrotation=60, xticks::Bool=false, kwargs...)
    y_vals = maximum(hcat(vecs...), dims=2) |> vec
    scale = sum(y_vals)

    p = Plots.plot(;widen=false, ylims=(0, maximum(y_vals) + text_offset * scale), kwargs...)
    for (v,l) in zip(vecs, labels)
        p = Plots.bar!(v, alpha=alpha, label=l)
    end

    ann_genes = findall(y_vals .>= min_expr_frac * scale)
    p = Plots.annotate!(collect(zip(ann_genes, y_vals[ann_genes] .+ text_offset * scale, Plots.text.(gene_names[ann_genes], fontsize))))
    if xticks
        Plots.xticks!(1:length(gene_names), gene_names, xrotation=xrotation, xgrid=false)
    end
    return p
end

function clustermap(mtx::T where T <: AbstractMatrix{Float64}, gene_names::Vector{String}; gene_ord::Union{<:AbstractVector{<:Integer}, Bool}=true,
        cell_ord::Union{<:AbstractVector{<:Integer},Bool}=true, diag_genes::Bool=false, kwargs...)
    if typeof(cell_ord) === Bool
        if cell_ord
            cell_dists = 1 .- cor(mtx);
            cell_ord = Clustering.hclust(cell_dists, linkage=:ward).order;
        else
            cell_ord = 1:size(mtx, 2)
        end
    end

    if typeof(gene_ord) === Bool
        if gene_ord
            if diag_genes
                gene_ord = sortperm(vec(mapslices(x -> findmax(x)[2], mtx[:, cell_ord], dims=2)), rev=true);
            else
                gene_dists = 1 .- cor(mtx');
                gene_ord = Clustering.hclust(gene_dists, linkage=:ward).order;
            end
        else
            gene_ord = 1:size(mtx, 1)
        end
    end

    Plots.heatmap(mtx[gene_ord, cell_ord]; yticks=(1:length(gene_names), gene_names[gene_ord]), kwargs...), cell_ord, gene_ord
end

### Tracing

function plot_num_of_cells_per_iterarion(tracer::Dict{Symbol, Any}; kwargs...)
    if !(:n_components in keys(tracer)) || length(tracer[:n_components]) == 0
        error("No data about #components per iteration was stored")
    end

    n_components_per_iter = hcat(collect.(values.(tracer[:n_components]))...);
    labels = collect(keys(tracer[:n_components][1]));

    n_components_per_iter = n_components_per_iter[sortperm(labels),:]
    labels = sort(labels)

    p = Plots.plot(; legendtitle="Min #molecules", title="Convergence", xlabel="Iteration", ylabel="#Cells",
                   background_color_legend=Colors.RGBA(1.0, 1.0, 1.0, 0.5), legend=:topleft, kwargs...)
    for i in 1:size(n_components_per_iter, 1)
        p = Plots.plot!(n_components_per_iter[i,:], label="$(labels[i])")
    end

    p = Plots.plot!(xlabel="Iteration", ylabel="#Cells", inset = (1, Plots.bbox(0.05,0.05,0.35,0.35,:top,:right)), subplot=2,
        bg_inside=nothing, legend=nothing, kwargs...)

    subplot_start = ceil(Int, 0.75 * size(n_components_per_iter, 2))
    iter_vals = subplot_start:size(n_components_per_iter, 2)
    for i in 1:size(n_components_per_iter, 1)
        p = Plots.plot!(iter_vals, n_components_per_iter[i, iter_vals], subplot=2, ylims=[0, maximum(n_components_per_iter[:, iter_vals]) * 1.05])
    end

    return p
end

### Colormaps

function map_to_colors(vals::Array{T, 1} where T; lims=nothing, palette=Colors.sequential_palette(0, 11))
    offset = (lims === nothing) ? minimum(vals) : lims[1]
    scale = (lims === nothing) ? maximum(vals) - offset : (lims[2] - lims[1])

    if lims !== nothing
        vals = min.(max.(vals, lims[1]), lims[2])
    end

    color_ticks = collect(range(0.0, 1.0, length=length(palette))) .* scale .+ offset
    colors = palette[floor.(Int, ((vals .- offset) ./ scale) .* (length(palette) - 1) .+ 1)]

    return Dict(:colors=>colors, :ticks=>color_ticks, :palette=>palette)
end

function distinguishable_colors(vals::T where T <: AbstractVector, args...; kwargs...)
    id_per_type = Dict(c => i for (i,c) in enumerate(unique(vals)));
    colors_uniq = Colors.distinguishable_colors(length(id_per_type), args...; kwargs...);
    colors = colors_uniq[get.(Ref(id_per_type), vals, 0)];

    return Dict(:colors=>colors, :ticks=>unique(vals), :palette=>colors_uniq)
end

plot_colorbar(colors; kwargs...) = plot_colorbar(colors[:ticks], colors[:palette]; kwargs...)

function plot_colorbar(color_ticks, palette; size=(500, 60), rotation=0, kwargs...)
    p = Plots.bar(color_ticks, ones(length(palette)); color=palette, size=size, legend=false, yticks=false, kwargs...)
    p.subplots[1][:xaxis][:rotation] = rotation;

    return p
end

### Utils

function Base.show(io::IO, m::MIME"text/html", fig::AbstractPlotting.Scene)
    io64 = Base64.IOBuffer();
    iob64 = Base64.Base64EncodePipe(io64)
    show(iob64, MIME("image/png"), fig);
    close(iob64)
    str = String(take!(io64))
    close(io64)

    print(io, "")
    # show(io, MIME("text/plain"), "")
end

max_stars_repo_name: YinuoJin/Baysor

max_issues_repo_name: YinuoJin/Baysor

max_forks_repo_name: YinuoJin/Baysor


text:
module KnapsackLib

export KnapData, KnapItem, solveKnapModel, solveKnapNaive, solveKnapExpCore, solveKnapHeur, genKnap

using Random, JuMP

const EPS = 1e-7

include("Data.jl")
include("Model.jl")
include("Naive.jl")
include("ExpandingCore.jl")
include("Heuristic.jl")
include("Generator.jl")

end


text:
__precompile__()
module Ska

export MSID, MSIDset, interpolate
export CxoTime
#export XijaModel, make, calc
export CmdStates, Commands, export_to_dataframe

include("cheta.jl")
include("time.jl")
include("kadi.jl")
#include("maude.jl")
#include("xija.jl")

import .cheta: MSID, MSIDset, interpolate
import .time: CxoTime
import .kadi: CmdStates, Commands, export_to_dataframe

#import .xija: XijaModel, make, calc

end


text:
using WAV

# create a CMD command that twill extract information from a WAV file given as a full path string.
function command(file; TIME_STEP = 0.01)
    # read wav file
    x,fs = WAV.wavread(file)
    x = x / maximum(x)

    # save wav file with WAVE_FORMAT_PCM/16bits
    WAV.wavwrite(x, raw"extraFiles\Temp.wav",
                 Fs=fs,
                 compression=WAVE_FORMAT_PCM,
                 nbits = 16) # this doesn't work with other files than test.wav in WavFiles

    # .exe praat file
    path_PraatExecutable = raw".\extraFiles\praatcon.exe"
    # praat script filemyF
    path_PraatScript = raw".\extraFiles\AutoRhythmPitch.praat"

    temp = raw" Temp.wav "
    features = raw"Temp.features "

    cmd = path_PraatExecutable * " " * path_PraatScript * temp  * features * raw".\ " * string(TIME_STEP)

    return cmd
end


text:
using StateSpaceRoutines, JLD2, HDF5, Random
path = dirname(@__FILE__)

# Initialize arguments to function
@load "$path/reference/kalman_filter_args.jld2" y T R C Q Z D E z0 P0

file  = h5open("$path/reference/kalman_filter_out.h5", "r")
pred  = read(file, "pred")
vpred = read(file, "vpred")
filt  = read(file, "filt")
vfilt = read(file, "vfilt")
close(file)

# Run smoothers
states = Dict{Symbol, Matrix{Float64}}()
shocks = Dict{Symbol, Matrix{Float64}}()

states[:hamilton], shocks[:hamilton] = hamilton_smoother(y, T, R, C, Q, Z, D, E, z0, P0)
states[:koopman], shocks[:koopman] = koopman_smoother(y, T, R, C, Q, Z, D, E, z0, P0, pred, vpred)
states[:carter_kohn], shocks[:carter_kohn] = carter_kohn_smoother(y, T, R, C, Q, Z, D, E, z0, P0; draw_states = false)
states[:durbin_koopman], shocks[:durbin_koopman] = durbin_koopman_smoother(y, T, R, C, Q, Z, D, E, z0, P0; draw_states = false)

# Check that last-period smoothed states equal last-period filtered states
@testset "Ensure last smoothed state = the last filtered state" begin
    for smoother in [:hamilton, :koopman, :carter_kohn, :durbin_koopman]
        @test filt[:, end] ≈ states[smoother][:, end]
    end
end

# Compare to expected output
exp_states, exp_shocks = h5open("$path/reference/smoothers_out.h5", "r") do file
    read(file, "exp_states"), read(file, "exp_shocks")
end

@testset "Check expected output for all smoothers" begin
    for smoother in [:hamilton, :koopman, :carter_kohn, :durbin_koopman]
        @test exp_states ≈ states[smoother]
        @test exp_shocks ≈ shocks[smoother]
    end
end

# Make sure that simulation smoothers run with `draw_states` on
carter_kohn_smoother(y, T, R, C, Q, Z, D, E, z0, P0; draw_states = true)
durbin_koopman_smoother(y, T, R, C, Q, Z, D, E, z0, P0; draw_states = true)

# TODO: tests for regime-switching and make sure they match on states across smoothers too, including where the data is missing
# Furthermore, you should test that the implied observables are matching.
nreg = 7
TTTs = Vector{Matrix{Float64}}(undef, nreg)
RRRs = Vector{Matrix{Float64}}(undef, nreg)
CCCs = Vector{Vector{Float64}}(undef, nreg)
ZZs = Vector{Matrix{Float64}}(undef, nreg)
DDs = Vector{Vector{Float64}}(undef, nreg)
QQs = Vector{Matrix{Float64}}(undef, nreg)
EEs = Vector{Matrix{Float64}}(undef, nreg)
regime_inds0 = Vector{Vector{Int64}}(undef, nreg)
regime_inds = Vector{UnitRange{Int64}}(undef, nreg)

for i in 1:nreg
    TTTs[i] = h5read("reference/time_varying_system.h5", "T$(i)")
    RRRs[i] = h5read("reference/time_varying_system.h5", "R$(i)")
    CCCs[i] = h5read("reference/time_varying_system.h5", "C$(i)")
    ZZs[i] = h5read("reference/time_varying_system.h5", "Z$(i)")
    DDs[i] = h5read("reference/time_varying_system.h5", "D$(i)")
    QQs[i] = h5read("reference/time_varying_system.h5", "Q$(i)")
    EEs[i] = h5read("reference/time_varying_system.h5", "E$(i)")
    regime_inds0[i] = h5read("reference/time_varying_system.h5", "regime_inds$(i)")
    regime_inds[i] = regime_inds0[i][1]:regime_inds0[i][end]
end

y = h5read("reference/time_varying_system.h5", "data")

# Run smoothers
states = Dict{Symbol, Matrix{Float64}}()
shocks = Dict{Symbol, Matrix{Float64}}()

z0 = zeros(size(TTTs[1], 1))
P0 = StateSpaceRoutines.solve_discrete_lyapunov(TTTs[1], RRRs[1] * QQs[1] * RRRs[1]')
_, pred, vpred, ~ = kalman_filter(regime_inds, y, TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs, z0, P0)

states[:hamilton], shocks[:hamilton] =
    hamilton_smoother(regime_inds, y, TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs, z0, P0)
states[:koopman], shocks[:koopman] =
    koopman_smoother(regime_inds, y, TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs, z0, P0, pred, vpred)
states[:carter_kohn], shocks[:carter_kohn] =
    carter_kohn_smoother(regime_inds, y, TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs, z0, P0; draw_states = false)
states[:durbin_koopman], shocks[:durbin_koopman] =
    durbin_koopman_smoother(regime_inds, y, TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs, z0, P0; draw_states = false)

Random.seed!(1793)
states[:carter_kohn_draw], shocks[:carter_kohn_draw] =
    carter_kohn_smoother(regime_inds, y, TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs, z0, P0; draw_states = true)
states[:durbin_koopman_draw], shocks[:durbin_koopman_draw] =
    durbin_koopman_smoother(regime_inds, y, TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs, z0, P0; draw_states = true)

obs = Dict{Symbol, Matrix{Float64}}()
for k in [:hamilton, :koopman, :carter_kohn, :durbin_koopman, :carter_kohn_draw, :durbin_koopman_draw]
    obs[k] = zeros(size(y, 1), size(states[:hamilton], 2))
    for (i, inds) in enumerate(regime_inds)
        obs[k][:, inds] = ZZs[i] * states[k][:, inds] .+ DDs[i]
    end
end

@testset "Regime-switching in T, R, etc., with smoothers" begin
    @test states[:hamilton] ≈ states[:carter_kohn]
    @test states[:koopman] ≈ states[:durbin_koopman]
    @test shocks[:hamilton] ≈ shocks[:carter_kohn]
    @test shocks[:koopman] ≈ shocks[:durbin_koopman]
    @test maximum(abs.(states[:hamilton] - states[:koopman])) < 5e-3
    @test maximum(abs.(shocks[:hamilton] - shocks[:koopman])) < 1.5e-3
    for i in 1:size(y, 1)
        not_nan = findall(.!isnan.(y[i, :]))
        if !isempty(not_nan)
            for k in [:hamilton, :koopman, :carter_kohn, :durbin_koopman]
                @test obs[k][i, not_nan] ≈ y[i, not_nan] atol=6e-6
            end
            if i in [1, 3, 7, 8, 13]
                @test obs[:carter_kohn_draw][i, not_nan] ≈ y[i, not_nan] atol=1e-2
            elseif i in [4, 5]
                @test obs[:carter_kohn_draw][i, not_nan] ≈ y[i, not_nan] atol=5e-4
            else
                @test obs[:carter_kohn_draw][i, not_nan] ≈ y[i, not_nan] atol=5e-4
            end
            @test obs[:durbin_koopman_draw][i, not_nan] ≈ y[i, not_nan] atol=5e-4
        end
    end
end

nothing


text:
# helper functions to use adm along with POMDPs.jl

"""
    state_to_scene(pomdp::SingleOCPOMDP, s::SingleOCState)
Convert an SingleOCState object to a scene
"""
function state_to_scene(pomdp::SingleOCPOMDP, s::SingleOCState)
    scene = Scene()
    ego = Vehicle(s.ego, pomdp.ego_type, EGO_ID)
    ped = Vehicle(s.ped, pomdp.ped_type, 2)
    push!(scene, ego)
    push!(scene, ped)
    return scene
end

"""
    rand_scene(rng::AbstractRNG, b::Dict{Int64, SparseCat}, env::CrosswalkEnv)
Sample a scene from a belief representation
"""
function rand_scene(rng::AbstractRNG, b::Dict{Int64, SparseCat}, env::CrosswalkEnv)
    scene = Scene()
    for id in keys(b)
        s = rand(rng, b[id])
        if scene.n == 0
            ego = Vehicle(s.ego, pomdp.ego_type, EGO_ID)
            push!(scene, ego)
        end
        ped = Vehicle(s.ped, PEDESTRIAN_DEF, id)
        push!(scene, ped)
    end
    return scene
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
"""
# JuBasicMD.jl

A collection of basic routines for Molecular Dynamics simulations written in
Julia.  These include:
* Euler–Maruyama (EM)
* Random Walk Metropolis (RWM)
* Metropolis Adjusted Langevin (MALA)
* Brünger, Brooks and Karplus (BBK)
* ABOBA, BAOAB
* Verlet
* Hamiltonian/Hybrid Monte Carlo (HMC)

## REPL help
`?` followed by an sampler name (`?RWM`)to obtain infromation about individual
samplers.  Use `?` on `sample_trajectory!` and `sample_trajectory` to obtain
information about running the samplers.  Finally, `?` on `Options` will provide
information about how to set additional options, such as the number of
information.

## Performing Sampling

Having constructed a `sampler` structure, samplers are executed with the the
commands
```
sample_trajectory!(X, sampler, options=opts);
```
for an in place transformation of `X` or
```
Xvals = sample_trajectory(X₀, sampler);
```
to obtain a full trajectory.  Metropolis methods will also return the running
acceptance rate.

"""
module JuBasicMD

using LinearAlgebra

export sample_trajectory, sample_trajectory!,
    Options,
    RWM, MALA, EM, BBK, ABOBA, BAOAB, HMC, Verlet, GJF

include("types.jl")
include("sample.jl")
include("utils.jl")
# RWM methods
include("metropolis/zeroth_order/rwm.jl")
# MALA methods
include("metropolis/first_order/mala.jl")
# EM methods
include("nonmetropolis/first_order/em.jl")
# BBK methods
include("nonmetropolis/second_order/bbk.jl")
# ABOBA methods
include("nonmetropolis/second_order/aboba.jl")
# BAOAB methods
include("nonmetropolis/second_order/baoab.jl")
# HMC methods
include("metropolis/second_order/hmc.jl")
# Verlet methods
include("nonmetropolis/second_order/verlet.jl")
# G-JF methods
include("nonmetropolis/second_order/g_jf.jl")
end

max_stars_repo_name: cortner/JuBasicMD

max_issues_repo_name: cortner/JuBasicMD

max_forks_repo_name: cortner/JuBasicMD


text:
using Flux
using Flux: Tracker, mse, throttle
using Base.Iterators: repeated

function norm_input(x, min, max)
     # res = []
     res = zeros(size(x))
     for i=1:length(res)
          res[i] = (x[i]-min)/(max-min)
          # push!(res,[(x[i]-min)/(max-min)])
     end
     return res
end

function extract_output(y, min, max)
     res = deepcopy(y)
     for i=1:length(y)
          res[i] = y[i]*((max-min)+min)
     end
     return res
end

# global parameters
input_dim = 1
hidden_dim = 2
output_dim = 2
α = 0.1                 # training rate
iterations = 5000

#### dataset ####
# input (control value)
raw_input = [1 1 1 1 1 1 2 2 2 2 2 1 1 1 1]
# raw_input = [1 1 1 1 1 1 2 2 2 2 2 1 1 1 1; 1 1 1 1 1 1 2 2 2 2 2 1 1 1 1]

# output (rotation speed)
raw_labels = [15 15 15 15 15 15 30 30 30 30 30 15 15 15 15]
# raw_labels = [15 15 15 15 15 15 30 30 30 30 30 15 15 15 15; 15 15 15 15 15 15 30 30 30 30 30 15 15 15 15]

x = norm_input(raw_input,0,10)
y = norm_input(raw_labels,0,50)

dataset = collect(repeated((x,y),iterations))
# dataset = (x,y)

#### model ####

# >>> without using Flux
# W1 = param(randn(hidden_dim, input_dim))
# b1 = param(randn(hidden_dim))
# layer1(x) = W1 * x .+ b1
#     # hidden neurons in between the two layers
# W2 = param(randn(output_dim, hidden_dim))
# b2 = param(randn(hidden_dim))
# layer2(x) = W2 * x .+ b2
# model(x) = layer2(layer1(x))

# using CuArrays (GPU)
# W, b, x, y = cu.((W, b, x, y))

# >>> with Flux
model = Chain(
    Dense(input_dim, hidden_dim),
    Dense(hidden_dim, output_dim))

#### loss function ####
loss(x, y) = mse(model(x), y)

# dataset = zip(x, y)
# dataset = repeated((x, y),100)
# opt = SGD(params(m), 0.1)          # Gradient descent with learning rate
# evalcb = () -> @show(loss(x, y))   # output the progress

#### training ####
function update!(ps, η = .1)
  for w in ps
    w.data .-= w.grad .* η
    w.grad .= 0
  end
end

for i = 1:length(dataset)
    xs = dataset[i][1]
    ys = dataset[i][2]
    back!(loss(xs, ys))             # automatic backpropagation
    update!(params(model), α)
    # update!((W1, b1, W2, b2), α)

    if i%1000 == 0  # do not print every step
        @show loss(x, y)
    end
end

# Flux.train!(loss, dataset, opt, cb = throttle(evalcb, 0.5))

# output final error
# @show(loss(x, y))

# testing the results
# @show m([1,0,0])    # should be close to 1
# @show m([0,1,0])    # should be close to 0


text: collect(i % 15 == 0 ? "FizzBuzz" : i % 5 == 0 ? "Buzz" : i % 3 == 0 ? "Fizz" : i for i in 1:100) |> println


text:
using LinearAlgebra, DiffEqOperators, Random, Test, BandedMatrices, SparseArrays

# 4th derivative, 4th order
function fourth_deriv_approx_stencil(N)
    A = zeros(N,N+2)
    A[1,1:8] = [3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3]
    A[2,1:8] = [2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6]

    A[N-1,N-5:end] = reverse([2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6], dims=2)
    A[N,N-5:end] = reverse([3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3], dims=2)

    for i in 3:N-2
        A[i,i-2:i+4] = [-1/6 2.0 -13/2 28/3 -13/2 2.0 -1/6]
    end
    return A
end

function second_deriv_fourth_approx_stencil(N)
    A = zeros(N,N+2)
    A[1,1:6] = [5/6 -15/12 -1/3 7/6 -6/12 5/60]
    A[N,N-3:end] = [1/12 -6/12 14/12 -4/12 -15/12 10/12]
    for i in 2:N-1
        A[i,i-1:i+3] = [-1/12 4/3 -5/2 4/3 -1/12]
    end
    return A
end

function second_derivative_stencil(N)
  A = zeros(N,N+2)
  for i in 1:N, j in 1:N+2
      (j-i==0 || j-i==2) && (A[i,j]=1)
      j-i==1 && (A[i,j]=-2)
  end
  A
end

@test_skip begin
@testset "Test Constructor, Multiplication, and Concretization" begin
    # Generate random parameters
    al = rand()
    bl = rand()
    cl = rand()

    ar = rand()
    br = rand()
    cr = rand()
    dx = rand()

    Q = RobinBC((al, bl, cl), (ar, br, cr), dx)
    N = 20
    L = CenteredDifference(4, 4, dx, N)
    L2 = CenteredDifference(2,4, dx, N)
    L1 = UpwindDifference(1, 1, dx, N, 1.)

    function coeff_func(du,u,p,t)
        du .= u
    end

    cL = coeff_func*L
    coeffs = rand(N)
    DiffEqOperators.update_coefficients!(cL,coeffs,nothing,0.0)
    cA = cL*Q
    DiffEqOperators.update_coefficients!(cA,coeffs,nothing,0.0)
    @test cL.coefficients == coeffs == cA.L.coefficients

    cL1 = coeff_func*L1
    coeffs = rand(N)
    DiffEqOperators.update_coefficients!(cL1,coeffs,nothing,0.0)
    cA1 = cL1*Q
    DiffEqOperators.update_coefficients!(cA1,coeffs,nothing,0.0)
    @test cL1.coefficients == coeffs == cA1.L.coefficients

    A = L*Q
    c2A = coeff_func*A
    DiffEqOperators.update_coefficients!(c2A,coeffs,nothing,0.0)
    @test coeffs == c2A.L.coefficients

    # Test GhostDerivativeOperator constructor by *
    u = rand(N)
    A = L*Q
    A1 = L1*Q

    # Test for consistency of GhostDerivativeOperator*u with L*(Q*u)
    @test A == L * Q
    @test A*u ≈ L*(Q*u)
    @test A1 == L1 * Q
    @test A1*u ≈ L1*(Q*u)

    # Test for consistency of c*GhostDerivativeOperator*u with alternative methods
    c = 2.1
    @test c * A == (c * L) * Q == c * (L * Q)
    @test c * A * u ≈ (c * L) * (Q * u) ≈ c * (L * Q) * u
    @test c * A1 == (c * L1) * Q == c * (L1 * Q)
    @test c * A1 * u ≈ (c * L1) * (Q * u) ≈ c * (L1 * Q) * u

    # check A + B, where A and B are GhostDerivativeOperators
    B = c * A
    B1 = c * A1
    @test (A + B) * u == (A + c * A) * u == (B + A) * u == (c * A + A) * u
    @test (A + B) * u == A * u + B * u
    @test (A1 + B1) * u == (A1 + c * A1) * u == (B1 + A1) * u == (c * A1 + A1) * u
    @test (A1 + B1) * u == A1 * u + B1 * u

    # Check (L + L) * Q works
    LLQ = (L + L) * Q
    LLQ1 = (L1 + L1) * Q
    @test LLQ * u == A * u + A * u == (A + A) * u
    @test LLQ1 * u == A1 * u + A1 * u == (A1 + A1) * u

    # Test for consistency of c*GhostDerivativeOperator*u when c is a vector
    c = rand(N)
    L = CenteredDifference(4, 4, dx, N)
    L1 = UpwindDifference(1, 1, 1., N, 1.)
    A1 = L1 * Q
    cA = c * A
    cL = c * L
    cLQ = (c * L) * Q
    cA1 = c * A1
    cL1 = c * L1
    cLQ1 = (c * L1) * Q
    @test cA.L == cL
    @test cA.L.coefficients == cL.coefficients == cLQ.L.coefficients
    @test c * A == (c * L) * Q == c * (L * Q)
    @test c * A * u ≈ (c * L) * (Q * u) ≈ c * (L * Q) * u
    @test cA1.L.coefficients == cL1.coefficients == cLQ1.L.coefficients
    @test c * A1 == (c * L1) * Q == c * (L1 * Q)
    @test c * A1 * u ≈ (c * L1) * (Q * u) ≈ c * (L1 * Q) * u

    # check A + B, where A and B are GhostDerivativeOperators
    B = c * A
    B1 = c * A1
    # @test (A + B) == A + c * A == B + A == c * A + A # uncomment if we want to implement and test for equality of linear combinations of operators
    @test (A + B) * u == (A + c * A) * u == (B + A) * u == (c * A + A) * u
    @test (A + B) * u == A * u + B * u
    # @test (A1 + B1) == A + c * A1 == B1 + A1 == c * A1 + A1
    @test (A1 + B1) * u == (A1 + c * A1) * u == (B1 + A1) * u == (c * A1 + A1) * u
    @test (A1 + B1) * u == A1 * u + B1 * u

    # Test for consistency of GhostDerivativeOperator*M with L*(Q*M)

    M = rand(N,10)
    Qx = MultiDimBC{1}(Q, size(M))
    Am = L*Qx
    LQM = zeros(N,10)
    for i in 1:10
        mul!(view(LQM,:,i), L, Q*M[:,i])
    end
    ghost_LQM = Am*M
    @test ghost_LQM ≈ LQM

    u = rand(N + 2)
    @test (L + L2) * u ≈ convert(AbstractMatrix,L + L2) * u ≈ (BandedMatrix(L) + BandedMatrix(L2)) * u

    # Test concretization
    @test Array(A)[1] ≈ (Array(L)*Array(Q,N)[1], Array(L)*Array(Q,N)[2])[1]
    @test Array(A)[2] ≈ (Array(L)*Array(Q,N)[1], Array(L)*Array(Q,N)[2])[2]
    @test SparseMatrixCSC(A)[1] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[2])[1]
    @test SparseMatrixCSC(A)[2] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[2])[2]
    @test sparse(A)[1] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[1]
    @test sparse(A)[2] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[2]
    # BandedMatrix not implemented for boundary operator
    @test_broken BandedMatrix(A)[1] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[1]
    @test_broken BandedMatrix(A)[2] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[2]

    # Test that concretization works with multiplication
    u = rand(N)
    @test Array(A)[1]*u + Array(A)[2] ≈ L*(Q*u) ≈ A*u
    @test sparse(A)[1]*u + sparse(A)[2] ≈ L*(Q*u) ≈ A*u
end

@testset "Test Constructor, Multiplication, and Concretization (Non-uniform grid)" begin
    # Generate random parameters
    al = rand()
    bl = rand()
    cl = rand()

    ar = rand()
    br = rand()
    cr = rand()
    N = 20
    dx = rand(N + 2)

    Q = RobinBC((al, bl, cl), (ar, br, cr), dx)
    L = UpwindDifference(1, 1, dx, N + 1, 1.)
    L2 = CenteredDifference(2, 2, dx, N)
    L4 = CenteredDifference(4, 4, dx, N - 2)

    function coeff_func(du,u,p,t)
        du .= u
    end

    cL = coeff_func*L
    coeffs = rand(N + 1)
    DiffEqOperators.update_coefficients!(cL,coeffs,nothing,0.0)
    cA = cL*Q
    DiffEqOperators.update_coefficients!(cA,coeffs,nothing,0.0)
    @test cL.coefficients == coeffs == cA.L.coefficients

    cL2 = coeff_func*L2
    coeffs = rand(N)
    DiffEqOperators.update_coefficients!(cL2,coeffs,nothing,0.0)
    cA2 = cL2*Q
    DiffEqOperators.update_coefficients!(cA2,coeffs,nothing,0.0)
    @test cL2.coefficients == coeffs == cA2.L.coefficients

    cL4 = coeff_func*L4
    coeffs = rand(N - 2)
    DiffEqOperators.update_coefficients!(cL4,coeffs,nothing,0.0)
    cA4 = cL4*Q
    DiffEqOperators.update_coefficients!(cA4,coeffs,nothing,0.0)
    @test cL4.coefficients == coeffs == cA4.L.coefficients

    # Test GhostDerivativeOperator constructor by *
    u = rand(N + 1)
    u2 = rand(N)
    u4 = rand(N - 2)
    A = L*Q
    A2 = L2*Q
    A4 = L4*Q

    # Test for consistency of GhostDerivativeOperator*u with L*(Q*u)
    @test A == L * Q
    @test A*u ≈ L*(Q*u)
    @test A2 == L2 * Q
    @test A2*u2 ≈ L2*(Q*u2)
    @test A4 == L4 * Q
    @test A4*u4 ≈ L4*(Q*u4)

    # Test for consistency of c*GhostDerivativeOperator*u with alternative methods
    c = 2.1
    cA = c * A
    cL = c * L
    @test c * A == (c * L) * Q == c * (L * Q)
    @test c * A * u ≈ (c * L) * (Q * u) ≈ c * (L * Q) * u
    @test c * A2 == (c * L2) * Q == c * (L2 * Q)
    @test c * A2 * u2 ≈ (c * L2) * (Q * u2) ≈ c * (L2 * Q) * u2
    @test c * A4 == (c * L4) * Q == c * (L4 * Q)
    @test c * A4 * u4 ≈ (c * L4) * (Q * u4) ≈ c * (L4 * Q) * u4

    # check A + B, where A and B are GhostDerivativeOperators
    B = c * A
    B2 = c * A2
    B4 = c * A4
    @test (A + B) * u == (A + c * A) * u == (B + A) * u == (c * A + A) * u
    @test (A + B) * u == A * u + B * u
    @test (A2 + B2) * u2 == (A2 + c * A2) * u2 == (B2 + A2) * u2 == (c * A2 + A2) * u2
    @test (A2 + B2) * u2 == A2 * u2 + B2 * u2
    @test (A4 + B4) * u4 == (A4 + c * A4) * u4 == (B4 + A4) * u4 == (c * A4 + A4) * u4
    @test (A4 + B4) * u4 == A4 * u4 + B4 * u4

    # Check (L + L) * Q works
    LLQ = (L + L) * Q
    LLQ2 = (L2 + L2) * Q
    LLQ4 = (L4 + L4) * Q
    @test LLQ * u == A * u + A * u == (A + A) * u
    @test LLQ2 * u2 == A2 * u2 + A2 * u2 == (A2 + A2) * u2 # this comparison fails, even though the operators seem alike
    @test LLQ4 * u4 == A4 * u4 + A4 * u4 == (A4 + A4) * u4

    # Test for consistency of c*GhostDerivativeOperator*u when c is a vector
    c = rand(N + 1)
    L1 = L
    A1 = L1 * Q
    cA = c * A
    cL = c * L
    cLQ = (c * L) * Q
    cA1 = c * A1
    cL1 = c * L1
    cLQ1 = (c * L1) * Q
    @test cA.L == cL
    @test cA.L.coefficients == cL.coefficients == cLQ.L.coefficients
    @test c * A == (c * L) * Q == c * (L * Q)
    @test c * A * u ≈ (c * L) * (Q * u) ≈ c * (L * Q) * u
    @test cA1.L.coefficients == cL1.coefficients == cLQ1.L.coefficients
    @test c * A1 == (c * L1) * Q == c * (L1 * Q)
    @test c * A1 * u ≈ (c * L1) * (Q * u) ≈ c * (L1 * Q) * u

    # check A + B, where A and B are GhostDerivativeOperators
    B = c * A
    # @test (A + B) == A + c * A == B + A == c * A + A # uncomment if implement equality of combinations of linear operators
    @test (A + B) * u == (A + c * A) * u == (B + A) * u == (c * A + A) * u
    @test (A + B) * u == A * u + B * u

    # Test for consistency of GhostDerivativeOperator*M with L*(Q*M)
    M = rand(N+1,10)
    LQM = zeros(N+1,10)
    for i in 1:10
        mul!(view(LQM,:,i), L, Q*M[:,i])
    end
    ghost_LQM = A*M
    @test ghost_LQM ≈ LQM

    # Test concretization for UpwindDifference
    @test Array(A)[1] ≈ (Array(L)*Array(Q,N + 1)[1], Array(L)*Array(Q,N + 1)[2])[1]
    @test Array(A)[2] ≈ (Array(L)*Array(Q,N + 1)[1], Array(L)*Array(Q,N + 1)[2])[2]
    @test SparseMatrixCSC(A)[1] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N + 1)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N + 1)[2])[1]
    @test SparseMatrixCSC(A)[2] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N + 1)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N + 1)[2])[2]
    @test sparse(A)[1] ≈ (sparse(L)*sparse(Q,N + 1)[1], sparse(L)*sparse(Q,N + 1)[2])[1]
    @test sparse(A)[2] ≈ (sparse(L)*sparse(Q,N + 1)[1], sparse(L)*sparse(Q,N + 1)[2])[2]
    # BandedMatrix not implemeted for boundary operator
    @test_broken BandedMatrix(A)[1] ≈ (BandedMatrix(L)*BandedMatrix(Q,N + 1)[1], BandedMatrix(L)*BandedMatrix(Q,N + 1)[2])[1]
    @test_broken BandedMatrix(A)[2] ≈ (BandedMatrix(L)*BandedMatrix(Q,N + 1)[1], BandedMatrix(L)*BandedMatrix(Q,N + 1)[2])[2]

    # Test concretization for CenteredDifference
    A2 = L2 * Q
    @test Array(A2)[1] ≈ (Array(L2)*Array(Q,N)[1], Array(L2)*Array(Q,N)[2])[1]
    @test Array(A2)[2] ≈ (Array(L2)*Array(Q,N)[1], Array(L2)*Array(Q,N)[2])[2]
    @test SparseMatrixCSC(A2)[1] ≈ (SparseMatrixCSC(L2)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L2)*SparseMatrixCSC(Q,N)[2])[1]
    @test SparseMatrixCSC(A2)[2] ≈ (SparseMatrixCSC(L2)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L2)*SparseMatrixCSC(Q,N)[2])[2]
    @test sparse(A2)[1] ≈ (sparse(L2)*sparse(Q,N)[1], sparse(L2)*sparse(Q,N)[2])[1]
    @test sparse(A2)[2] ≈ (sparse(L2)*sparse(Q,N)[1], sparse(L2)*sparse(Q,N)[2])[2]
    # BandedMatrix not implemeted for boundary operator
    @test_broken BandedMatrix(A2)[1] ≈ (BandedMatrix(L2)*BandedMatrix(Q,N)[1], BandedMatrix(L2)*BandedMatrix(Q,N)[2])[1]
    @test_broken BandedMatrix(A2)[2] ≈ (BandedMatrix(L2)*BandedMatrix(Q,N)[1], BandedMatrix(L2)*BandedMatrix(Q,N)[2])[2]

    # Test that concretization works with multiplication, UpwindDifference
    u = rand(N + 1)
    @test Array(A)[1]*u + Array(A)[2] ≈ L*(Q*u) ≈ A*u
    @test sparse(A)[1]*u + sparse(A)[2] ≈ L*(Q*u) ≈ A*u

    # Test that concretization works with multiplication, CenteredDifference
    u = rand(N)
    @test Array(A2)[1]*u + Array(A2)[2] ≈ L2*(Q*u) ≈ A2*u
    @test sparse(A2)[1]*u + sparse(A2)[2] ≈ L2*(Q*u) ≈ A2*u
end

@testset "Test Left Division L2 (second order)" begin
    # Test \ Homogenous BC
    # f(x) = -x^2 + x
    # f''(x) = -2.0
    # f'(0) = f'(1) = 0
    f(x) = -x^2 + x
    f2(x) = -2.0

    dx = 0.01
    x = dx:dx:1.0-dx
    N = length(x)

    L = CenteredDifference(2, 2, dx, N)
    Q = RobinBC((1.0, 0.0, 0.0), (1.0, 0.0, 0.0), dx)
    A = L*Q

    analytic_L = second_derivative_stencil(N) ./ dx^2
    analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
    analytic_AL = analytic_L*analytic_QL

    analytic_QM = zeros((length(x)+2)*2, length(x)*2)
    interior = CartesianIndices(Tuple([2:length(x)+1, 1:2]))
    I1 = CartesianIndex(1,0)
    for I in interior
        i = DiffEqOperators.cartesian_to_linear(I, (length(x)+2, 2)) #helper function, see utils.jl
        j = DiffEqOperators.cartesian_to_linear(I-I1, (length(x), 2))
        analytic_QM[i,j] = 1.0
    end
    analytic_Am = kron(Diagonal(ones(2)), analytic_L)*analytic_QM

    # No affine component to the this system
    analytic_f = analytic_AL \ f2.(x)
    ghost_f = A \ f2.(x)

    # Check that A\f2.(x) is consistent with analytic_AL \ f2.(x)
    @test analytic_f ≈ ghost_f

    # Additionally test that A\f2.(x) ≈ f.(x)
    @test f.(x) ≈ ghost_f ≈ analytic_f

    # Check that left division with matrices works
    M = [f2.(x) f2.(x)]
    Qx = MultiDimBC{1}(Q, size(M))

    Am = L*Qx
    ghost_fM = Am \ M
    s = size(M)
    analytic_fM = analytic_Am \ reshape(M, prod(s))
    @test ghost_fM ≈ reshape(analytic_fM, s)

    fM_temp = zeros(N,2)

    # Test \ Inhomogenous BC
    # f(x) = -x^2 + x + 4.0
    # f''(x) = -2.0
    # f'(0) = f'(1) = 0
    f(x) = -x^2 + x + 4.0
    f2(x) = -2.0

    dx = 0.01
    x = dx:dx:1.0-dx
    N = length(x)

    L = CenteredDifference(2, 2, dx, N)
    Q = RobinBC((1.0, 0.0, 4.0), (1.0, 0.0, 4.0), dx)
    A = L*Q

    analytic_L = second_derivative_stencil(N) ./ dx^2
    analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
    analytic_Qb = [4.0; zeros(N); 4.0]
    analytic_AL = analytic_L*analytic_QL
    analytic_Ab = analytic_L*analytic_Qb
    analytic_Lm = kron(Diagonal(ones(3)), analytic_L)

    analytic_f = analytic_AL \ (f2.(x) - analytic_Ab)
    ghost_f = A \ f2.(x)

    # Check that A\f2.(x) is consistent with analytic_AL \ f2.(x)
    @test analytic_f ≈ ghost_f

    # Additionally test that A\f2.(x) ≈ f.(x)
    @test f.(x) ≈ ghost_f ≈ analytic_f

    # Check \ for Matrix
    M2 = [f2.(x) 2.0*f2.(x) 10.0*f2.(x)]

    s = size(M2)
    Qx = MultiDimBC{1}(Q, size(M2))

    analytic_QLm, analytic_Qbm = Array(Qx, s)

    analytic_ALm = analytic_Lm*analytic_QLm
    analytic_Abm = analytic_Lm*analytic_Qbm

    Am = L*Qx
    analytic_M = analytic_ALm \ (reshape(M2 , prod(s)).- analytic_Abm)
    ghost_M = Am \ M2
    @test reshape(analytic_M, s) ≈ ghost_M

    # Additionally test that A\M2 ≈ [f, 2.0(f-4.0)+4.0, 10.0(f-4.0)+4.0]
    M = [f.(x) 2.0*(f.(x) .- 4.0).+4.0 10.0*(f.(x) .- 4.0).+4.0]
    @test M ≈ reshape(analytic_M, s) ≈ ghost_M
end

@testset "Test Left Division L4 (fourth order)" begin

    # Test \ homogenous and inhomogenous BC
    dx = 0.01
    x = 0.01:dx:0.2
    N = length(x)
    u = sin.(x)

    L = CenteredDifference(4, 4, dx, N)
    Q = RobinBC((1.0, 0.0, 0.0), (1.0, 0.0, sin(0.2+dx)), dx)
    A = L*Q

    analytic_L = fourth_deriv_approx_stencil(N) ./ dx^4
    analytic_QL = [transpose(zeros(N)); Diagonal(ones(N)); transpose(zeros(N))]
    analytic_AL = analytic_L*analytic_QL
    analytic_Qb = [zeros(N+1); sin(0.2+dx)]
    analytic_Ab = analytic_L*analytic_Qb

    analytic_QM = zeros((length(x)+2)*3, length(x)*3)
    interior = CartesianIndices(Tuple([2:length(x)+1, 1:3]))
    I1 = CartesianIndex(1,0)
    for I in interior
        i = DiffEqOperators.cartesian_to_linear(I, (length(x)+2, 3)) #helper function, see utils.jl
        j = DiffEqOperators.cartesian_to_linear(I-I1, (length(x), 3))
        analytic_QM[i,j] = 1.0
    end
    analytic_Am = kron(Diagonal(ones(3)), analytic_L)*analytic_QM

    analytic_u = analytic_AL \ Vector(u .- analytic_Ab)
    ghost_u = A \ u

    As_l, As_b = sparse(A, length(u))

    @test As_l ≈ analytic_AL #This is true
    @test As_b ≈  analytic_Ab #This is true
    @test Vector(u .- analytic_Ab) ≈ Vector(u .- As_b) #this is true
    u_translated = u .- analytic_Ab |> Vector
    # TODO validate analytic_AL as having reasonable condition (what is reasonable condition?)
    @test norm(analytic_AL \ u_translated - As_l \ u_translated) / norm(analytic_AL \ u_translated) < cond(analytic_AL)*eps()

    sp_u = As_l\Vector(u .- As_b)
    # Check that A\u.(x) is consistent with analytic_AL \ u.(x)
    @test ghost_u ≈ sp_u
    # TODO make cond function in ghost_derivative_operator?
    # however, note: ArgumentError: 2-norm condition number is not implemented for sparse matrices, try cond(Array(A), 2) instead
    # so the actual function would probably just propagate that error
    _cond(A::GhostDerivativeOperator, u) = cond(sparse(A,length(u))[1] |> Array)
    @test norm(analytic_u - ghost_u) / norm(ghost_u) < _cond(A,u) * eps() #

    M2 = [u 2.0*u 10.0*u]
    s = size(M2)
    Qx = MultiDimBC{1}(Q, size(M2))
    Am = L*Qx
    #Somehow the operator is singular
    @test_broken analytic_M = analytic_Am \ (reshape(M2, prod(s)) .- repeat(analytic_Ab, 3))
    @test_broken ghost_M = Am \ M2
    @test_broken reshape(analytic_M, s) ≈ ghost_M

end

@testset "Test Operator and BC combinations" begin
    N = 40
    x = range(-pi, stop = pi, length=N)
    Δx = x[2]-x[1]
    u₀=1.0
    Γ=1.0
    Dx=u₀*CenteredDifference{1}(1,4,Δx,N)
    Dxx=Γ*CenteredDifference{1}(2,4,Δx,N)
    Q=NeumannBC((2x[1]+2, 2x[end]+2), Δx, 6)

    A = Dx*Q + Dxx*Q

    y = A*(x.^2)
    analytic_y = 2x.+2

    @test_broken y ≈ analytic_y
end
end


text:
using FractionalDiffEq

function ϕ(x)
    if x == 0
        return 19.00001
    else
        return 19.0
    end
end

f(t, y, ϕ) = 3.5*y*(1-ϕ/19)

h = 0.05
α = 0.97
τ = 0.8
T = 56
fddeprob = FDDEProblem(f, ϕ, α, τ, T)
V, y = solve(fddeprob, h, DelayPECE())

using Plots
plot(y, V, xlabel="y(t)", ylabel="y(t-τ)")


text:
using DiffEqBase, DiffEqJump
using Test
const DJ = DiffEqJump

# test data
minpriority = 2.0^exponent(1e-12)
maxpriority = 2.0^exponent(1e12)
priorities = [1e-13, .99*minpriority, minpriority,1.01e-4, 1e-4, 5., 0., 1e10]

mingid = exponent(minpriority)   # = -40
ptog = priority -> DJ.priortogid(priority, mingid)
pt = DJ.PriorityTable(ptog, priorities, minpriority, maxpriority)

#display(priorities)
#display(pt)

# test insert
grpcnt = DJ.numgroups(pt)
push!(priorities, maxpriority*.99)
DJ.insert!(pt, length(priorities), priorities[end])
@test grpcnt == DJ.numgroups(pt)
@test pt.groups[end].pids[1] == length(priorities)

push!(priorities, maxpriority*.99999)
DJ.insert!(pt, length(priorities), priorities[end])
@test grpcnt == DJ.numgroups(pt)
@test pt.groups[end].pids[2] == length(priorities)

numsmall = length(pt.groups[2].pids)
push!(priorities, minpriority*.6)
DJ.insert!(pt, length(priorities), priorities[end])
@test grpcnt == DJ.numgroups(pt)
@test pt.groups[2].pids[end] == length(priorities)

push!(priorities, maxpriority)
DJ.insert!(pt, length(priorities), priorities[end])
@test grpcnt == DJ.numgroups(pt)-1
@test pt.groups[end].pids[1] == length(priorities)

# test updating
DJ.update!(pt, 5, priorities[5], 2*priorities[5])   # group 29
priorities[5] *= 2
@test pt.groups[29].numpids == 1
@test pt.groups[30].numpids == 1

DJ.update!(pt, 9, priorities[9], maxpriority*1.01)
priorities[9] = maxpriority*1.01
@test pt.groups[end].numpids == 2
@test pt.groups[end-1].numpids == 1

DJ.update!(pt, 10, priorities[10], 0.)
priorities[10] = 0.
@test pt.groups[1].numpids == 2

# test sampling
cnt = 0
Nsamps = Int(1e7)
for i = 1:Nsamps
    global cnt; pid = DJ.sample(pt, priorities)
    (pid == 8) && (cnt += 1)
end
@test abs(cnt // Nsamps -  0.008968535978248484) / 0.008968535978248484 < .05


text:
#=
Core functions

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones, and others (see Julia contributors)
Licensed under MIT License, see LICENSE.md

Inspired by / derived from code in Julia
=#

_lastindex(::CodeUnitSingle, str) = (@_inline_meta(); ncodeunits(str))

@propagate_inbounds function _getindex(::CodeUnitSingle, T, str, pos::Int)
    @_inline_meta()
    @boundscheck checkbounds(str, pos)
    T(get_codeunit(pointer(str), pos))
end

@propagate_inbounds function _next(::CodeUnitSingle, T, str, pos)
    @_inline_meta()
    @boundscheck 0 < pos <= ncodeunits(str) || boundserr(str, pos)
    T(get_codeunit(str, pos)), pos + 1
end

_nextcpfun(::CodeUnitSingle, ::Type{S}, pnt::Ptr{T}) where {S,T<:CodeUnitTypes} =
    get_codeunit(pnt), pnt + sizeof(T)
_nextcp(::Type{T}, pnt) where {T} = _nextcpfun(CodePointStyle(T), T, pnt)

@propagate_inbounds _getindex(::CodeUnitMulti, T, str, pos::Int) =
    first(_next(CodeUnitMulti(), T, str, pos))

@inline _length(::CodeUnitSingle, str) = ncodeunits(str)

# Use more generic length check
@inline _length_check(str::SubString{<:Str{C}}, cnt) where {C<:CSE} =
    _length(CodeUnitMulti(), C, pointer(str), cnt)

# Go directly to aligned length check
@inline _length_check(str::Str{C}, cnt) where {C<:CSE} =
    _length_al(CodeUnitMulti(), C, pointer(str), cnt)

@inline _length(::CodeUnitMulti, str::MaybeSub{T}) where {T<:Str} =
    (cnt = ncodeunits(str); cnt < 2 ? Int(cnt > 0) : @preserve str _length_check(str, cnt))

@inline _length(::CodeUnitSingle, ::Type{<:CSE}, ::Ptr{<:CodeUnitTypes}, cnt::Int) = cnt

@propagate_inbounds function _length(cs::CodePointStyle, str, i::Int, j::Int)
    @boundscheck begin
        # I think the bounds of these should be 1:siz
        lim = ncodeunits(str)+1
        0 <  i <= lim || boundserr(str, i)
        0 <= j <  lim || boundserr(str, j)
    end
    (cnt = j - i + 1) <= 1 ? Int(cnt > 0) :
        @preserve str _length(cs, cse(str), bytoff(pointer(str), i - 1), cnt)
end

@inline _thisind(::CodeUnitSingle, str, len, pnt, pos) = Int(pos)

@propagate_inbounds function _thisind(cs::CS, str, pos) where {CS<:CodePointStyle}
    @_inline_meta()
    # I do think thisind should not return anything outside of the valid range
    # but for now, to make it compatible with the current String API, do this:
    pos == 0 && return 0
    len = ncodeunits(str)
    pos == len + 1 && return pos
    @boundscheck 0 < pos <= len || boundserr(str, pos)
    @preserve str _thisind(cs, str, len, pointer(str), pos)
end

@propagate_inbounds function _prevind(::CodeUnitSingle, str, i)
    @_inline_meta()
    @boundscheck 0 < i <= ncodeunits(str)+1 || boundserr(str, i)
    Int(i) - 1
end

@propagate_inbounds function _prevind(::CodeUnitSingle, str, i, nchar)
    @_inline_meta()
    nchar < 0 && ncharerr(nchar)
    @boundscheck 0 < i <= ncodeunits(str)+1 || boundserr(str, i)
    max(Int(i) - nchar, 0)
end

@propagate_inbounds function _nextind(::CodeUnitSingle, str, i)
    @_inline_meta()
    @boundscheck 0 <= i <= ncodeunits(str) || boundserr(str, i)
    Int(i) + 1
end

@propagate_inbounds function _nextind(::CodeUnitSingle, str, i, nchar)
    @_inline_meta()
    nchar < 0 && ncharerr(nchar)
    @boundscheck 0 <= i <= ncodeunits(str) || boundserr(str, i)
    min(Int(i) + nchar, ncodeunits(str)+1)
end

_index(cs::CodePointStyle, str, i)               = _thisind(cs, str, i)
_index(cs::CodePointStyle, ::Fwd, str, i)        = _nextind(cs, str, i)
_index(cs::CodePointStyle, ::Fwd, str, i, nchar) = _nextind(cs, str, i, nchar)
_index(cs::CodePointStyle, ::Rev, str, i)        = _prevind(cs, str, i)
_index(cs::CodePointStyle, ::Rev, str, i, nchar) = _prevind(cs, str, i, nchar)

#  Call to specialized version via trait
@propagate_inbounds lastindex(str::MaybeSub{T}) where {T<:Str} =
    (@_inline_meta(); _lastindex(CodePointStyle(T), str))
@propagate_inbounds getindex(str::MaybeSub{T}, i::Int) where {T<:Str} =
    (@_inline_meta(); R = eltype(T) ; _getindex(CodePointStyle(T), R, str, i)::R)
@propagate_inbounds next(str::T, i::Int) where {T<:Str} =
    (@_inline_meta(); R = eltype(T) ; _next(CodePointStyle(T), R, str, i)::Tuple{R,Int})
@propagate_inbounds next(str::SubString{T}, i::Int) where {T<:Str} =
    (@_inline_meta(); R = eltype(T) ; _next(CodePointStyle(T), R, str, i)::Tuple{R,Int})
@propagate_inbounds length(str::MaybeSub{T}) where {T<:Str} =
    (@_inline_meta(); _length(CodePointStyle(T), str))
@propagate_inbounds length(str::MaybeSub{T}, i::Int, j::Int) where {T<:Str} =
    (@_inline_meta(); _length(CodePointStyle(T), str, i, j))
@propagate_inbounds thisind(str::MaybeSub{T}, i::Int) where {T<:Str} =
    (@_inline_meta(); _thisind(CodePointStyle(T), str, i))
@propagate_inbounds prevind(str::MaybeSub{T}, i::Int) where {T<:Str} =
    (@_inline_meta(); _prevind(CodePointStyle(T), str, i))
@propagate_inbounds nextind(str::MaybeSub{T}, i::Int) where {T<:Str} =
    (@_inline_meta(); _nextind(CodePointStyle(T), str, i))
@propagate_inbounds prevind(str::MaybeSub{T}, i::Int, nchar::Int) where {T<:Str} =
    (@_inline_meta(); _prevind(CodePointStyle(T), str, i, nchar))
@propagate_inbounds nextind(str::MaybeSub{T}, i::Int, nchar::Int) where {T<:Str} =
    (@_inline_meta(); _nextind(CodePointStyle(T), str, i, nchar))

@propagate_inbounds index(str::MaybeSub{T}, i::Integer) where {T<:Str} =
    (@_inline_meta(); _index(CodePointStyle(T), str, Int(i)))
@propagate_inbounds index(::D, str::MaybeSub{T}, i::Integer) where {T<:Str,D<:Direction} =
    (@_inline_meta(); _index(CodePointStyle(T), D(), str, Int(i)))
@propagate_inbounds index(::D, str::MaybeSub{T}, i::Integer,
                          nchar::Integer) where {T<:Str,D<:Direction} =
    (@_inline_meta(); _index(CodePointStyle(T), D(), str, Int(i), Int(nchar)))

@propagate_inbounds reverseind(str::MaybeSub{T}, i::Integer) where {T<:Str} =
    (@_inline_meta(); _index(CodePointStyle(T), str, Int(ncodeunits(str) - i + 1)))

@static if V6_COMPAT
    @propagate_inbounds ind2chr(str::MaybeSub{T}, i::Int) where {T<:Str} = length(str, 1, i)
    @propagate_inbounds chr2ind(str::MaybeSub{T}, i::Int) where {T<:Str} = nextind(str, 0, i)
end

@propagate_inbounds function is_valid(str::MaybeSub{T}, i::Integer) where {T<:Str}
    @_inline_meta()
    @boundscheck 1 <= i <= ncodeunits(str) || return false
    _isvalid_char_pos(CodePointStyle(T), cse(T), str, i)
end

_isvalid_char_pos(::CodeUnitSingle, C, str, i) = true

@propagate_inbounds function _collectstr(::CodeUnitMulti, ::Type{S},
                                         str::MaybeSub{T}) where {S,T<:Str}
    len = _length(CodeUnitMulti(), str)
    vec = create_vector(S, len)
    pos = 1
    @inbounds for i = 1:len
        vec[i], pos = _next(CodeUnitMulti(), S, str, pos)
    end
    vec
end

@propagate_inbounds function _collectstr(::CodeUnitSingle, ::Type{S},
                                         str::MaybeSub{T}) where {S,T<:Str}
    len = ncodeunits(str)
    vec = create_vector(S, len)
    cpt = eltype(T)
    @preserve str begin
        pnt = pointer(str)
        if S == cpt
            unsafe_copyto!(reinterpret(Ptr{basetype(cpt)}, pointer(vec)), pnt, len)
        else
            @inbounds for i = 1:len
                vec[i] = T(get_codeunit(pnt, i))
            end
        end
    end
    vec
end

function map(fun, str::MaybeSub{T}) where {C<:CSE, T<:Str{C}}
    out = get_iobuffer(sizeof(str))
    CP = eltype(T)
    for ch in str
        retc = fun(ch)
        isa(retc, AbstractChar) || throw(ArgumentError(
            "map($fun, str::AbstractString) requires $fun to return AbstractChar; " *
            "try map($fun, collect(str)) or a comprehension instead"))
        is_valid(CP, retc) || codepoint_error(CP, retc)
        _write(C, out, retc)
    end
    Str{C}(String(take!(out)))
end

@propagate_inbounds collect(str::MaybeSub{T}) where {T<:Str} =
    @preserve str _collectstr(CodePointStyle(T), eltype(T), str)

# An optimization here would be to check just if they are the same type, but
# rather if they are the same size with a compatible encoding, i.e. like
# UTF32Chr and UInt32, but not Char and UInt32.

@propagate_inbounds Base._collect(::Type{S}, str::T, isz::Base.HasLength) where {S,T<:Str} =
    _collectstr(CodePointStyle(T), S, str)

@inline function check_valid(ch, pos)
    is_surrogate_codeunit(ch) && unierror(UTF_ERR_SURROGATE, pos, ch)
    ch <= 0x10ffff || unierror(UTF_ERR_INVALID, pos, ch)
    ch
end

# Convert single characters to strings
convert(::Type{T}, ch::Char) where {T<:Str} = convert(T, UInt32(ch))

function _convert(::Type{C}, ch::T) where {C<:CSE,T<:CodeUnitTypes}
    buf, pnt = _allocate(T, 1)
    set_codeunit!(pnt, ch)
    Str(C, buf)
end

# Todo: These should be made more generic, work for all CodeUnitSingle types

convert(::Type{<:Str{ASCIICSE}}, ch::Unsigned) =
    is_ascii(ch) ? _convert(ASCIICSE, ch%UInt8) : unierror(UTF_ERR_INVALID_ASCII, 0, ch)
convert(::Type{<:Str{LatinCSE}}, ch::Unsigned) =
    is_latin(ch) ? _convert(LatinCSE, ch%UInt8) : unierror(UTF_ERR_INVALID_LATIN1, 0, ch)
convert(::Type{<:Str{UCS2CSE}}, ch::Unsigned) =
    is_bmp(ch) ? _convert(UCS2CSE, ch%UInt16) : unierror(UTF_ERR_INVALID, 0, ch)
convert(::Type{<:Str{UTF32CSE}}, ch::Unsigned) =
    is_unicode(ch) ? _convert(UTF32CSE, ch%UInt32) : unierror(UTF_ERR_INVALID, 0, ch)

convert(::Type{T}, ch::Signed) where {T<:Str} = ch < 0 ? ncharerr(ch) : convert(T, ch%Unsigned)

Str(str::SubString{<:Str{C}}) where {C<:Byte_CSEs} =
    Str(C, unsafe_string(pointer(str.string, str.offset+1), str.ncodeunits))

# don't make unnecessary copies when passing substrings to C functions
cconvert(::Type{Ptr{UInt8}}, str::SubString{<:Str{<:Byte_CSEs}}) = str
cconvert(::Type{Ptr{Int8}},  str::SubString{<:Str{<:Byte_CSEs}}) = str

unsafe_convert(::Type{Ptr{UInt8}}, s::SubString{<:Str{Byte_CSEs}}) =
    convert(Ptr{UInt8}, pointer(s.string)) + s.offset
unsafe_convert(::Type{Ptr{Int8}}, s::SubString{<:Str{Byte_CSEs}}) =
    convert(Ptr{Int8}, pointer(s.string)) + s.offset

function _reverse(::CodeUnitSingle, ::Type{C}, len, str::Str{C}) where {C<:CSE}
    len < 2 && return str
    @preserve str begin
        pnt = pointer(str)
        T = codeunit(C)
        buf, beg = _allocate(T, len)
        out = bytoff(beg, len)
        while out > beg
            set_codeunit!(out -= sizeof(T), get_codeunit(pnt))
            pnt += sizeof(T)
        end
        Str(C, buf)
    end
end

function _reverse(::CodeUnitMulti, ::Type{C}, len, str) where {C<:CSE}
    @inbounds ((t = nextind(str, 0)) > len || nextind(str, t) > len) && return str
    @preserve str _reverse(CodeUnitMulti(), C, len, pointer(str))
end

reverse(str::MaybeSub{T}) where {C<:CSE,T<:Str{C}} =
    _reverse(CodePointStyle(T), C, ncodeunits(str), str)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/ModelReduction.jl/blob/master/LICENSE

using ModelReduction

@testset "Perform Guyan Reduction" begin
    K = [
         1.0 -1.0  0.0  0.0  0.0
        -1.0  2.0 -1.0  0.0  0.0
         0.0 -1.0  2.0 -1.0  0.0
         0.0  0.0 -1.0  2.0 -1.0
         0.0  0.0  0.0 -1.0  1.0]
    m = [1, 5]
    s = [2, 3, 4]
    expected = [0.25 -0.25; -0.25 0.25]
    result = ModelReduction.guyan_reduction(K, m, s)
    @test isapprox(result, expected)
end


text:
 export VmoveAllChildren, MoveAll, MoveNodeToLeft, MoveNodeToRight

#======================================================================================#
#    Similar to MoveAll() in LayoutBuild.jl
#======================================================================================#
# may be: TextLine -or- Element
function VmoveAllChildren(node, y::Float64, moveNode::Bool)

  shape = node.shape
  if shape.flags[Fixed]
    return
  end
  if moveNode
      shape.top  += y
  end
    if isdefined(node, :rows) # ..it has rows of children so let's move them!
      for i in 1:length(node.rows)
        row = node.rows[i]
        row.top += y
        for j in 1:length(row.nodes)
            VmoveAllChildren(row.nodes[j],y, true) # do the same for each child
        end
      end
    end
end
#======================================================================================#
#    This is only to translate a shape with all children by x,y
#    It may also be nesesary to move the node to another location in its row!
#======================================================================================#
# may be: TextLine -or- Element
function MoveAll(node, x::Float64, y::Float64)
  shape = node.shape
  if shape.flags[Fixed]
    return
  end
  shape.left += x # Move this object!
  shape.top  += y
    if isdefined(node, :rows) # ..it has rows of children so let's move them!
      for i in 1:length(node.rows)
        row = node.rows[i]
        row.left += x # ...don't forget to move the actual row
        row.top += y
        for j in 1:length(row.nodes)
            MoveAll(row.nodes[j],x,y) # do the same for each child
        end
      end
    end
end
#======================================================================================#
#======================================================================================#
function MoveNodeToLeft(row::Row, index::Int64)
  node = row.nodes[index]
  for n in index:-1:2
    row.nodes[n] = row.nodes[n-1]
  end
  row.nodes[1] = node
end
#======================================================================================#
function MoveNodeToRight(row::Row, index::Int64)

  node = row.nodes[index]
  for n in index:length(row.nodes)-1
    row.nodes[n] = row.nodes[n+1]
  end
  row.nodes[end] = node
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

module MyDict

using Actors
import Actors: spawn

struct DictSrv{L}
    lk::L
end
(ds::DictSrv)() = call(ds.lk)
(ds::DictSrv)(f::Function, args...) = call(ds.lk, f, args...)
# indexing interface
Base.getindex(d::DictSrv, key) = call(d.lk, getindex, key)
Base.setindex!(d::DictSrv, value, key) = call(d.lk, setindex!, value, key)

# dict server behavior
ds(d::Dict, f::Function, args...) = f(d, args...)
ds(d::Dict) = copy(d)
# start dict server
dictsrv(d::Dict; remote=false) = DictSrv(spawn(ds, d, remote=remote))

export DictSrv, dictsrv

end

using .MyDict, .Threads

d = dictsrv(Dict{Int,Int}())

@threads for i in 1:1000
    d[i] = threadid()
end

d()