text:
# ============================================================
# Solver Initializer
# ============================================================

"""
    initialize(configfilename::T) where {T<:AbstractString}
    initialize(config::T) where {T<:AbstractDict}

Initialize solver from input file or dictionary.
This can also be done from a Julia script directly.

"""
function initialize(configfilename::T) where {T<:AbstractString}
    println("==============================================================")
    println("Kinetic.jl")
    println("Portable Kinetic Simulation and Scientific Machine Learning")
    println("==============================================================")
    println("")
    @info "initializing solver"
    println("")

    if configfilename[end-3:end] == "jld2"
        D = load(configfilename)
        ks = D["set"]
        ctr = D["ctr"]
        t = D["t"]
        face = init_fvm(ks, ks.pSpace)[2]

        return ks, ctr, face, t
    else
        ks = SolverSet(configfilename)
        cftuple = init_fvm(ks)

        return (ks, cftuple..., 0.0)
    end
end

function initialize(config::T) where {T<:AbstractDict}
    println("==============================================================")
    println("Kinetic.jl")
    println("A Lightweight Toolbox for Kinetic Modeling and Simulation")
    println("==============================================================")
    println("")
    @info "initializing solver"
    println("")

    ks = SolverSet(config)
    cftuple = init_fvm(ks)

    return (ks, cftuple..., 0.0)
end

"""
    init_fvm(KS::T, ps::T1) where {T<:AbstractSolverSet,T1<:AbstractPhysicalSpace}

Initialize finite volume method

"""
function init_fvm(
    KS::T,
    array = :dynamic_array;
    structarray = false,
) where {T<:AbstractSolverSet}
    return init_fvm(KS, KS.ps, array; structarray = structarray)
end

function init_fvm(
    KS::T,
    ps::T1,
    array = :dynamic_array;
    structarray = false,
) where {T<:AbstractSolverSet,T1<:AbstractPhysicalSpace1D}

    funcar = eval(array)
    funcst = ifelse(structarray, StructArray, dynamic_array)

    funcprim = ifelse(KS.set.nSpecies == 1, conserve_prim, mixture_conserve_prim)
    γ = begin
        if KS.gas isa Scalar
            KS.gas.a
        else
            KS.gas.γ
        end
    end

    if KS.set.space[3:4] == "0f"

        ctr = OffsetArray{ControlVolume1D}(undef, eachindex(KS.pSpace.x))
        face = Array{Interface1D}(undef, KS.pSpace.nx + 1)

        for i in eachindex(ctr)
            w = KS.ib.fw(KS.ps.x[i])
            prim = begin
                if γ == 0
                    funcprim(w)
                else
                    funcprim(w, γ)
                end
            end

            ctr[i] = ControlVolume1D(funcar(w), funcar(prim))
        end

        for i = 1:KS.pSpace.nx+1
            fw = deepcopy(KS.ib.fw(KS.ps.x[1])) |> funcar
            face[i] = Interface1D(fw)
        end

    elseif KS.set.space[3:4] == "1f"

        ctr = OffsetArray{ControlVolume1D1F}(undef, eachindex(KS.pSpace.x))
        face = Array{Interface1D1F}(undef, KS.pSpace.nx + 1)

        for i in eachindex(ctr)
            w = KS.ib.fw(KS.ps.x[i])
            prim = funcprim(w, γ)
            f = KS.ib.ff(KS.ps.x[i])

            ctr[i] = ControlVolume1D1F(funcar(w), funcar(prim), funcar(f))
        end

        for i = 1:KS.pSpace.nx+1
            fw = deepcopy(KS.ib.fw(KS.ps.x[1])) |> funcar
            ff = deepcopy(KS.ib.ff(KS.ps.x[1])) |> funcar
            face[i] = Interface1D1F(fw, ff)
        end

    elseif KS.set.space[3:4] == "2f"

        ctr = OffsetArray{ControlVolume1D2F}(undef, eachindex(KS.pSpace.x))
        face = Array{Interface1D2F}(undef, KS.pSpace.nx + 1)

        for i in eachindex(ctr)
            w = KS.ib.fw(KS.ps.x[i])
            prim = funcprim(w, γ)
            h, b = KS.ib.ff(KS.ps.x[i])

            ctr[i] = ControlVolume1D2F(funcar(w), funcar(prim), funcar(h), funcar(b))
        end

        for i = 1:KS.pSpace.nx+1
            fw = deepcopy(KS.ib.fw(KS.ps.x[1])) |> funcar
            ff = deepcopy(KS.ib.ff(KS.ps.x[1])[1]) |> funcar
            face[i] = Interface1D2F(fw, ff)
        end

    elseif KS.set.space[3:4] == "3f"

        ctr = OffsetArray{ControlVolume1D3F}(undef, eachindex(KS.pSpace.x))
        face = Array{Interface1D3F}(undef, KS.pSpace.nx + 1)

        for i in eachindex(ctr)
            w = KS.ib.fw(KS.ps.x[i])
            prim = funcprim(w, γ)
            h0, h1, h2 = KS.ib.ff(KS.ps.x[i])
            E = KS.ib.fE(KS.ps.x[i])
            B = KS.ib.fB(KS.ps.x[i])
            L = KS.ib.fL(KS.ps.x[i])

            ctr[i] = ControlVolume1D3F(
                funcar(w),
                funcar(prim),
                funcar(h0),
                funcar(h1),
                funcar(h2),
                funcar(E),
                funcar(B),
                funcar(L),
            )
        end

        for i = 1:KS.pSpace.nx+1
            fw = deepcopy(KS.ib.fw(KS.ps.x[1])) |> funcar
            ff = deepcopy(KS.ib.ff(KS.ps.x[1])[1]) |> funcar
            fe = deepcopy(KS.ib.fE(KS.ps.x[1])) |> funcar
            face[i] = Interface1D3F(fw, ff, fe)
        end

    elseif KS.set.space[3:4] == "4f"

        ctr = OffsetArray{ControlVolume1D4F}(undef, eachindex(KS.pSpace.x))
        face = Array{Interface1D4F}(undef, KS.pSpace.nx + 1)

        for i in eachindex(ctr)
            w = KS.ib.fw(KS.ps.x[i])
            prim = funcprim(w, γ)
            h0, h1, h2, h3 = KS.ib.ff(KS.ps.x[i])
            E = KS.ib.fE(KS.ps.x[i])
            B = KS.ib.fB(KS.ps.x[i])
            L = KS.ib.fL(KS.ps.x[i])

            ctr[i] = ControlVolume1D4F(
                funcar(w),
                funcar(prim),
                funcar(h0),
                funcar(h1),
                funcar(h2),
                funcar(h3),
                funcar(E),
                funcar(B),
                funcar(L),
            )
        end

        for i = 1:KS.pSpace.nx+1
            fw = deepcopy(KS.ib.fw(KS.ps.x[1])) |> funcar
            ff = deepcopy(KS.ib.ff(KS.ps.x[1])[1]) |> funcar
            fe = deepcopy(KS.ib.fE(KS.ps.x[1])) |> funcar
            face[i] = Interface1D4F(fw, ff, fe)
        end

    end

    return ctr |> funcst, face |> funcst
end

function init_fvm(
    KS::T,
    ps::T1,
    array = :dynamic_array;
    structarray = false,
) where {T<:AbstractSolverSet,T1<:AbstractPhysicalSpace2D}

    funcar = eval(array)
    funcst = ifelse(structarray, StructArray, dynamic_array)
    funcprim = ifelse(KS.set.nSpecies == 1, conserve_prim, mixture_conserve_prim)

    nx, ny, dx, dy = begin
        if ps isa CSpace2D
            ps.nr, ps.nθ, ps.dr, ps.darc
        else
            ps.nx, ps.ny, ps.dx, ps.dy
        end
    end

    if KS.set.space[3:4] == "0f"

        ctr =
            OffsetArray{ControlVolume2D}(undef, axes(KS.pSpace.x, 1), axes(KS.pSpace.y, 2))
        a1face = Array{Interface2D}(undef, nx + 1, ny)
        a2face = Array{Interface2D}(undef, nx, ny + 1)

        for j in axes(ctr, 2), i in axes(ctr, 1)
            w = KS.ib.fw(KS.ps.x[i, j], KS.ps.y[i, j])
            prim = funcprim(w, KS.gas.γ)

            ctr[i, j] = ControlVolume2D(funcar(w), funcar(prim))
        end

        for j = 1:ny
            for i = 1:nx
                n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :])
                n .= ifelse(
                    dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0,
                    n,
                    -n,
                )

                a1face[i, j] = Interface2D1F(
                    point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :]),
                    n[1],
                    n[2],
                    funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                )
            end
            n = unit_normal(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :])
            n .= ifelse(
                dot(n, ps.vertices[nx, j, 2, :] .- [ps.x[nx, j], ps.y[nx, j]]) >= 0,
                n,
                -n,
            )

            a1face[nx+1, j] = Interface2D1F(
                point_distance(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]),
                n[1],
                n[2],
                funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
            )
        end
        for i = 1:nx
            for j = 1:ny
                n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :])
                n .= ifelse(
                    dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0,
                    n,
                    -n,
                )

                a2face[i, j] = Interface2D1F(
                    point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :]),
                    n[1],
                    n[2],
                    funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                )
            end
            n = unit_normal(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :])
            n .= ifelse(
                dot(n, ps.vertices[i, ny, 3, :] .- [ps.x[i, ny], ps.y[i, ny]]) >= 0,
                n,
                -n,
            )

            a2face[i, ny+1] = Interface2D1F(
                point_distance(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]),
                n[1],
                n[2],
                funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
            )
        end

    elseif KS.set.space[3:4] == "1f"

        ctr = OffsetArray{ControlVolume2D1F}(
            undef,
            axes(KS.pSpace.x, 1),
            axes(KS.pSpace.y, 2),
        )
        a1face = Array{Interface2D1F}(undef, nx + 1, ny)
        a2face = Array{Interface2D1F}(undef, nx, ny + 1)

        for j in axes(ctr, 2), i in axes(ctr, 1)
            w = KS.ib.fw(KS.ps.x[i, j], KS.ps.y[i, j])
            prim = funcprim(w, KS.gas.γ)
            h = KS.ib.ff(KS.ps.x[i, j], KS.ps.y[i, j])

            ctr[i, j] = ControlVolume2D1F(funcar(w), funcar(prim), funcar(h))
        end

        for j = 1:ny
            for i = 1:nx
                n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :])
                n .= ifelse(
                    dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0,
                    n,
                    -n,
                )

                a1face[i, j] = Interface2D1F(
                    point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :]),
                    n[1],
                    n[2],
                    funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                    funcar(KS.ib.ff(KS.ps.x[1], KS.ps.y[1])),
                )
            end
            n = unit_normal(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :])
            n .= ifelse(
                dot(n, ps.vertices[nx, j, 2, :] .- [ps.x[nx, j], ps.y[nx, j]]) >= 0,
                n,
                -n,
            )

            a1face[nx+1, j] = Interface2D1F(
                point_distance(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]),
                n[1],
                n[2],
                funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                funcar(KS.ib.ff(KS.ps.x[1], KS.ps.y[1])),
            )
        end
        for i = 1:nx
            for j = 1:ny
                n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :])
                n .= ifelse(
                    dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0,
                    n,
                    -n,
                )

                a2face[i, j] = Interface2D1F(
                    point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :]),
                    n[1],
                    n[2],
                    funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                    funcar(KS.ib.ff(KS.ps.x[1], KS.ps.y[1])),
                )
            end
            n = unit_normal(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :])
            n .= ifelse(
                dot(n, ps.vertices[i, ny, 3, :] .- [ps.x[i, ny], ps.y[i, ny]]) >= 0,
                n,
                -n,
            )

            a2face[i, ny+1] = Interface2D1F(
                point_distance(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]),
                n[1],
                n[2],
                funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                funcar(KS.ib.ff(KS.ps.x[1], KS.ps.y[1])),
            )
        end

    elseif KS.set.space[3:4] == "2f"

        ctr = OffsetArray{ControlVolume2D2F}(
            undef,
            axes(KS.pSpace.x, 1),
            axes(KS.pSpace.y, 2),
        )
        a1face = Array{Interface2D2F}(undef, nx + 1, ny)
        a2face = Array{Interface2D2F}(undef, nx, ny + 1)

        for j in axes(ctr, 2), i in axes(ctr, 1)
            w = KS.ib.fw(KS.ps.x[i, j], KS.ps.y[i, j])
            prim = funcprim(w, KS.gas.γ)
            h, b = KS.ib.ff(KS.ps.x[i, j], KS.ps.y[i, j])

            ctr[i, j] = ControlVolume2D2F(funcar(w), funcar(prim), funcar(h), funcar(b))
        end

        for j = 1:ny
            for i = 1:nx
                n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :])
                n .= ifelse(
                    dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0,
                    n,
                    -n,
                )

                a1face[i, j] = Interface2D2F(
                    point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :]),
                    n[1],
                    n[2],
                    funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                    funcar(KS.ib.ff(KS.ps.x[1], KS.ps.y[1])[1]),
                )
            end
            n = unit_normal(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :])
            n .= ifelse(
                dot(n, ps.vertices[nx, j, 2, :] .- [ps.x[nx, j], ps.y[nx, j]]) >= 0,
                n,
                -n,
            )

            a1face[nx+1, j] = Interface2D2F(
                point_distance(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]),
                n[1],
                n[2],
                funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                funcar(KS.ib.ff(KS.ps.x[1], KS.ps.y[1])[1]),
            )
        end
        for i = 1:nx
            for j = 1:ny
                n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :])
                n .= ifelse(
                    dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0,
                    n,
                    -n,
                )

                a2face[i, j] = Interface2D2F(
                    point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :]),
                    n[1],
                    n[2],
                    funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                    funcar(KS.ib.ff(KS.ps.x[1], KS.ps.y[1])[1]),
                )
            end
            n = unit_normal(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :])
            n .= ifelse(
                dot(n, ps.vertices[i, ny, 3, :] .- [ps.x[i, ny], ps.y[i, ny]]) >= 0,
                n,
                -n,
            )

            a2face[i, ny+1] = Interface2D2F(
                point_distance(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]),
                n[1],
                n[2],
                funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])),
                funcar(KS.ib.ff(KS.ps.x[1], KS.ps.y[1])[1]),
            )
        end

    end

    return ctr |> funcst, a1face |> funcst, a2face |> funcst

end

function init_fvm(
    KS::T,
    ps::UnstructPSpace,
    array = :dynamic_array;
    structarray = false,
) where {T<:AbstractSolverSet}

    funcar = eval(array)
    funcst = ifelse(structarray, StructArray, dynamic_array)
    funcprim = ifelse(KS.set.nSpecies == 1, conserve_prim, mixture_conserve_prim)

    if KS.set.space[3:4] == "0f"

        ctr = Array{KitBase.ControlVolumeUS}(undef, size(ps.cellid, 1))
        for i in eachindex(ctr)
            n = Vector{Float64}[]
            for j = 1:3
                push!(
                    n,
                    KitBase.unit_normal(
                        ps.points[ps.facePoints[ps.cellFaces[i, j], 1], :],
                        ps.points[ps.facePoints[ps.cellFaces[i, j], 2], :],
                    ) |> funcar,
                )

                if dot(
                    ps.faceCenter[ps.cellFaces[i, j], 1:2] .- ps.cellCenter[i, 1:2],
                    n[j],
                ) < 0
                    n[j] .= -n[j]
                end
            end

            dx =
                [
                    KitBase.point_distance(
                        ps.cellCenter[i, :],
                        ps.points[ps.cellid[i, 1], :],
                        ps.points[ps.cellid[i, 2], :],
                    ),
                    KitBase.point_distance(
                        ps.cellCenter[i, :],
                        ps.points[ps.cellid[i, 2], :],
                        ps.points[ps.cellid[i, 3], :],
                    ),
                    KitBase.point_distance(
                        ps.cellCenter[i, :],
                        ps.points[ps.cellid[i, 3], :],
                        ps.points[ps.cellid[i, 1], :],
                    ),
                ] |> funcar

            w = KS.ib.fw(ps.cellCenter[i, 1], ps.cellCenter[i, 2])
            prim = funcprim(w, KS.gas.γ)

            ctr[i] = KitBase.ControlVolumeUS(
                n,
                funcar(ps.cellCenter[i, :]),
                dx,
                funcar(w),
                funcar(prim),
            )
        end

        face = Array{KitBase.Interface2D}(undef, size(ps.facePoints, 1))
        for i in eachindex(face)
            len =
                norm(ps.points[ps.facePoints[i, 1], :] .- ps.points[ps.facePoints[i, 2], :])
            n = KitBase.unit_normal(
                ps.points[ps.facePoints[i, 1], :],
                ps.points[ps.facePoints[i, 2], :],
            )

            if !(-1 in ps.faceCells[i, :])
                n0 =
                    ps.cellCenter[ps.faceCells[i, 2], :] .-
                    ps.cellCenter[ps.faceCells[i, 1], :]
            else
                idx =
                    ifelse(ps.faceCells[i, 1] != -1, ps.faceCells[i, 1], ps.faceCells[i, 2])
                n0 = ps.cellCenter[idx, :] .- ps.faceCenter[i, :]
            end
            if dot(n, n0[1:2]) < 0
                n .= -n
            end

            fw = zero(KS.ib.fw(ctr[1].x[1], ctr[1].x[2])) |> funcar

            face[i] = KitBase.Interface2D(len, n[1], n[2], fw)
        end

    elseif KS.set.space[3:4] == "1f"

        ctr = Array{KitBase.ControlVolumeUS1F}(undef, size(ps.cellid, 1))
        for i in eachindex(ctr)
            n = Vector{Float64}[]
            for j = 1:3
                push!(
                    n,
                    KitBase.unit_normal(
                        ps.points[ps.facePoints[ps.cellFaces[i, j], 1], :],
                        ps.points[ps.facePoints[ps.cellFaces[i, j], 2], :],
                    ) |> funcar,
                )

                if dot(
                    ps.faceCenter[ps.cellFaces[i, j], 1:2] .- ps.cellCenter[i, 1:2],
                    n[j],
                ) < 0
                    n[j] .= -n[j]
                end
            end

            dx =
                [
                    KitBase.point_distance(
                        ps.cellCenter[i, :],
                        ps.points[ps.cellid[i, 1], :],
                        ps.points[ps.cellid[i, 2], :],
                    ),
                    KitBase.point_distance(
                        ps.cellCenter[i, :],
                        ps.points[ps.cellid[i, 2], :],
                        ps.points[ps.cellid[i, 3], :],
                    ),
                    KitBase.point_distance(
                        ps.cellCenter[i, :],
                        ps.points[ps.cellid[i, 3], :],
                        ps.points[ps.cellid[i, 1], :],
                    ),
                ] |> funcar

            w = KS.ib.fw(ps.cellCenter[i, 1], ps.cellCenter[i, 2])
            prim = funcprim(w, KS.gas.γ)
            h = KS.ib.ff(ps.cellCenter[i, 1], ps.cellCenter[i, 2])

            ctr[i] = KitBase.ControlVolumeUS1F(
                n,
                ps.cellCenter[i, :],
                dx,
                funcar(w),
                funcar(prim),
                funcar(h),
            )

        end

        face = Array{KitBase.Interface2D1F}(undef, size(ps.facePoints, 1))
        for i in eachindex(face)
            len =
                norm(ps.points[ps.facePoints[i, 1], :] .- ps.points[ps.facePoints[i, 2], :])
            n = KitBase.unit_normal(
                ps.points[ps.facePoints[i, 1], :],
                ps.points[ps.facePoints[i, 2], :],
            )

            if !(-1 in ps.faceCells[i, :])
                n0 =
                    ps.cellCenter[ps.faceCells[i, 2], :] .-
                    ps.cellCenter[ps.faceCells[i, 1], :]
            else
                idx =
                    ifelse(ps.faceCells[i, 1] != -1, ps.faceCells[i, 1], ps.faceCells[i, 2])
                n0 = ps.cellCenter[idx, :] .- ps.faceCenter[i, :]
            end
            if dot(n, n0[1:2]) < 0
                n .= -n
            end

            fw = zero(KS.ib.fw(ctr[1].x[1], ctr[1].x[2]))
            ff = zero(KS.ib.ff(ctr[1].x[1], ctr[1].x[2]))

            face[i] = KitBase.Interface2D1F(len, n[1], n[2], funcar(fw), funcar(ff))
        end

    elseif KS.set.space[3:4] == "2f"

        ctr = Array{KitBase.ControlVolumeUS2F}(undef, size(ps.cellid, 1))
        for i in eachindex(ctr)
            n = Vector{Float64}[]
            for j = 1:3
                push!(
                    n,
                    KitBase.unit_normal(
                        ps.points[ps.facePoints[ps.cellFaces[i, j], 1], :],
                        ps.points[ps.facePoints[ps.cellFaces[i, j], 2], :],
                    ) |> funcar,
                )

                if dot(
                    ps.faceCenter[ps.cellFaces[i, j], 1:2] .- ps.cellCenter[i, 1:2],
                    n[j],
                ) < 0
                    n[j] .= -n[j]
                end
            end

            dx =
                [
                    KitBase.point_distance(
                        ps.cellCenter[i, :],
                        ps.points[ps.cellid[i, 1], :],
                        ps.points[ps.cellid[i, 2], :],
                    ),
                    KitBase.point_distance(
                        ps.cellCenter[i, :],
                        ps.points[ps.cellid[i, 2], :],
                        ps.points[ps.cellid[i, 3], :],
                    ),
                    KitBase.point_distance(
                        ps.cellCenter[i, :],
                        ps.points[ps.cellid[i, 3], :],
                        ps.points[ps.cellid[i, 1], :],
                    ),
                ] |> funcar

            w = KS.ib.fw(ps.cellCenter[i, 1], ps.cellCenter[i, 2])
            prim = funcprim(w, KS.gas.γ)
            h, b = KS.ib.ff(ps.cellCenter[i, 1], ps.cellCenter[i, 2])

            ctr[i] = KitBase.ControlVolumeUS2F(
                n,
                funcar(ps.cellCenter[i, :]),
                dx,
                funcar(w),
                funcar(prim),
                funcar(h),
                funcar(b),
            )
        end

        face = Array{KitBase.Interface2D2F}(undef, size(ps.facePoints, 1))
        for i in eachindex(face)
            len =
                norm(ps.points[ps.facePoints[i, 1], :] .- ps.points[ps.facePoints[i, 2], :])
            n = KitBase.unit_normal(
                ps.points[ps.facePoints[i, 1], :],
                ps.points[ps.facePoints[i, 2], :],
            )

            if !(-1 in ps.faceCells[i, :])
                n0 =
                    ps.cellCenter[ps.faceCells[i, 2], :] .-
                    ps.cellCenter[ps.faceCells[i, 1], :]
            else
                idx =
                    ifelse(ps.faceCells[i, 1] != -1, ps.faceCells[i, 1], ps.faceCells[i, 2])
                n0 = ps.cellCenter[idx, :] .- ps.faceCenter[i, :]
            end
            if dot(n, n0[1:2]) < 0
                n .= -n
            end

            fw = zero(KS.ib.fw(ctr[1].x[1], ctr[1].x[2]))
            fh = zero(KS.ib.ff(ctr[1].x[1], ctr[1].x[2])[1])

            face[i] = KitBase.Interface2D2F(len, n[1], n[2], funcar(fw), funcar(fh))
        end

    end

    return ctr |> funcst, face |> funcst

end

"""
    init_ptc!(
        KS::SolverSet,
        ctr::T;
        mode = :soa::Symbol,
        factor = 1::Real,
    ) where {T<:AbstractArray{<:AbstractControlVolume1D,1}}

Initialize particles based on flow conditions

"""
function init_ptc!(
    KS::SolverSet,
    ctr::T;
    mode = :soa::Symbol,
    factor = 1::Real,
) where {T<:AbstractArray{<:AbstractControlVolume1D,1}}
    if mode == :soa
        init_ptc_soa!(KS, ctr, factor)
    elseif mode == :aos
        init_ptc_aos!(KS, ctr, factor)
    end
end

"""
    init_ptc_aos!(
        KS::SolverSet,
        ctr::T,
        factor = 1,
    ) where {T<:AbstractArray{<:AbstractControlVolume1D,1}}

Initialize particles with array of structs

"""
function init_ptc_aos!(
    KS::SolverSet,
    ctr::T,
    factor = 1,
) where {T<:AbstractArray{<:AbstractControlVolume1D,1}}

    np = 0
    for i in eachindex(ctr)
        np += Int(round(ctr[i].w[1] * KS.ps.dx[i] / KS.gas.m))
    end
    KS.gas.np = np
    np *= round(factor) |> Int

    ptc = Array{Particle1D}(undef, np)
    for i in eachindex(ptc)
        m = KS.gas.m
        x = 0.0
        v = zeros(3)
        e = 0.0
        idx = -7
        flag = 0
        tc = 0.0

        ptc[i] = Particle1D(m, x, v, e, idx, flag, tc)
    end

    np_tmp = 0
    for i in eachindex(ctr)
        npl = Int(round(ctr[i].w[1] * KS.ps.dx[i] / KS.gas.m))
        for j = 1:npl
            np_tmp += 1
            sample_particle!(ptc[np_tmp], KS, ctr[i], i, KS.ps.x[i], KS.ps.dx[i])
        end
    end

    return ptc |> StructArray

end

"""
    init_ptc_soa!(
        KS::SolverSet,
        ctr::T,
        factor = 1,
    ) where {T<:AbstractArray{<:AbstractControlVolume1D,1}}

Initialize particles with struct of arrays

"""
function init_ptc_soa!(
    KS::SolverSet,
    ctr::T,
    factor = 1,
) where {T<:AbstractArray{<:AbstractControlVolume1D,1}}

    np = 0
    for i in eachindex(ctr)
        np += round(ctr[i].prim[1] * KS.ps.dx[i] / KS.gas.m) |> Int
    end
    KS.gas.np = np
    np *= round(factor) |> Int
    np_tmp = 0

    m = zeros(np)
    x = zeros(np)
    v = zeros(np, 3)
    e = zeros(np)
    idx = zeros(Int, np)
    ref = zeros(Int, np) # default
    flag = zeros(Int, np) # default
    tc = zeros(np)

    # in-cell particles
    for i in eachindex(ctr)
        npl = Int(round(ctr[i].w[1] * KS.ps.dx[i] / KS.gas.m))
        for j = 1:npl
            np_tmp += 1

            m[np_tmp] = KS.gas.m
            x[np_tmp] = KS.ps.x[i] + (rand() - 0.5) * KS.ps.dx[i]
            v[np_tmp, :] .= sample_maxwell(ctr[i].prim)
            e[np_tmp] = 0.5 / ctr[i].prim[end]
            idx[np_tmp] = i
            flag[np_tmp] = 0
            if i < 1
                flag[np_tmp] = 1
            elseif i > KS.pSpace.nx
                flag[np_tmp] = 2
            end
            τ = vhs_collision_time(ctr[i].prim, KS.gas.μᵣ, KS.gas.ω)
            tc[np_tmp] = next_collision_time(τ)
        end
    end

    # placeholder particles
    for i = KS.gas.np+1:np
        m[i] = KS.gas.m
        x[i] = 0.0
        idx[i] = -7
        v[i, :] .= 0.0
        e[i] = 0.0
        τ = -1e8
        tc[i] = 1e8
    end

    ptc = Particle(m, x, v, e, idx, ref, flag, tc)

    return ptc

end


text:
"""
$(TYPEDEF)

Pose3 is currently a Euler angle mechanization of three Euclidean translations and three Circular rotation.

Future:
------
- Work in progress on AMP3D for proper non-Euler angle on-manifold operations.
"""
struct Pose3 <: IncrementalInference.InferenceVariable end
getDimension(::Pose3) = 6
# TODO the AMP upgrade is aimed at resolving 3D to Quat/SE3/SP3 -- current Euler angles will be replaced
getManifolds(::Pose3) = (:Euclid,:Euclid,:Euclid,:Circular,:Circular,:Circular)

ext: jl


text:
using SimilaritySearch, Random

function main(res, n=30_000_000)
    Random.seed!(0)
    @info "===========", typeof(res), n
    st = initialstate(res)
    st = push!(res, st, 1, rand())
    st = push!(res, st, 2, rand())

    for i in 3:n
        if length(res, st) > 0 && rand() < 0.001
            _, st = pop!(res, st)
        end

        if length(res, st) > 0 && rand() < 0.001
            _, st = popfirst!(res, st)
        end
        st = push!(res, st, i, rand())
    end

    res
end

ksearch = 10
n = 1000
main(KnnResultMatrix(ksearch), n)
main(KnnResult(ksearch), n)

n = 30_000_000
@timev main(KnnResultMatrix(ksearch), n)
@timev main(KnnResult(ksearch), n)

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

@userplot ExcessPlot

@recipe function f(p::ExcessPlot)
  # get user input
  xs = p.args[1]

  ks = 2:length(xs)
  ξs = meanexcess(xs, ks)

  seriestype --> :path
  xlabel --> "number of maxima"
  ylabel --> "extreme value index"

  ks, ξs
end


text:
using DrWatson, Test

@testset "Naming" begin include("naming_tests.jl"); end
@testset "Project Setup" begin include("project_tests.jl"); end
@testset "Saving tools" begin include("stools_tests.jl"); end
@testset "Produce or Save" begin include("savefiles_tests.jl"); end
@testset "Collect Results" begin include("update_results_tests.jl"); end
@testset "Parameter Customization" begin include("customize_savename.jl"); end


text:
@testset "package issues" begin
    @test isempty(detect_ambiguities(Optim))
end


text:
#!/usr/bin/env julia
using Gtk, Gtk.ShortNames

ui = Builder(filename=joinpath(@__DIR__, "ui.glade"))
showall(ui["win"])

function preview()
    get_gtk_property(ui["text"], :buffer, GtkTextBuffer) |> x -> get_gtk_property(x, :text, String) |> x -> set_gtk_property!(ui["label"], :label, x)

    return nothing
end

function open_about()
    showall(ui["about"])

    return nothing
end

signal_connect(x-> preview(), ui["preview_button"], "clicked")
signal_connect((x,y) -> open_about(), ui["menu_about"], :activate, Nothing, (), false)

signal_connect(ui["file_quit"], :activate) do w
  if !isinteractive()
    Gtk.gtk_quit()
  else
    exit()
  end
end

if !isinteractive()
    c = Condition()
    signal_connect(ui["win"], :destroy) do widget
        notify(c)
    end
    wait(c)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
using GeometryProcessingPlayground
using Test

@testset "GeometryProcessingPlayground.jl" begin
    # Write your tests here.
end


text:
#=
Loads a single image from a filename and converts to the expected 256x256 size
=#
function load_single_image(filename::String,imgs,ii,x_off,y_off)

    this_img=load(filename)

    x_off[ii] = 256 / size(this_img,2)
    y_off[ii] = 256 / size(this_img,1)

    new_img=imresize(this_img,(256,256))
    for i=1:size(imgs,2)
        for j=1:size(imgs,1)
            imgs[i,j,1,ii]=convert(Float32,new_img[j,i].r)
            imgs[i,j,2,ii]=convert(Float32,new_img[j,i].g)
            imgs[i,j,3,ii]=convert(Float32,new_img[j,i].b)
        end
    end

    nothing
end

function make_all_labels(ll,xo,yo,j_a)
    for i=1:size(ll,4)
        add_single_label(ll,i,xo,yo,j_a);
    end
end

function add_single_label(ll,ii,x_off,y_off,j_a,ds=4)

    for i=1:size(ll,3)
        xo=round(Int64,x_off[ii] * j_a[1,i,ii]/ds)
        yo=round(Int64,y_off[ii] * j_a[2,i,ii]/ds)

        ll[:,:,i,ii] = gaussian_2d(1.0:1.0:size(ll,1),1.0:1.0:size(ll,2),xo,yo)

        if maximum(ll[:,:,i,ii])>0
            ll[:,:,i,ii] = ll[:,:,i,ii] ./ maximum(ll[:,:,i,ii])
        end

    end
end

#=
Saving Hourglass Model
=#
function save_hourglass(name,x)
    count=[1];
    file=matopen(name,"w")
        save_nn(x,file,count)
    close(file)
end

function save_nn(x,file,count)

    for f in fieldnames(typeof(x))
        f1=getfield(x,f)
        if typeof(f1) <: NN
            save_nn(f1,file,count)
        elseif eltype(f1) <: NN
            for i=1:length(f1)
                save_nn(f1[i],file,count)
            end
        else
            if f == :w
                write(file,string("w_",count[1]),convert(Array,f1.value))
                count[1]+=1
            elseif f == :b
                write(file,string("b_",count[1]),convert(Array,f1.value))
                count[1]+=1
            elseif f == :ms
                write(file,string("ms_mo_",count[1]),f1.momentum)
                write(file,string("ms_mean_",count[1]),convert(Array,f1.mean))
                write(file,string("ms_var_",count[1]),convert(Array,f1.var))
                count[1]+=1
            elseif f == :bn_p
                write(file,string("bn_p_",count[1]),convert(Array,f1.value))
                count[1]+=1
            end
        end
    end
end

#=
Loading Hourglass Model
=#
function load_hourglass(name,x,atype=KnetArray)
    count=[1];
    file=matopen(name,"r")
        load_nn(x,file,count,atype)
    close(file)
end

function load_nn(x,file,count,atype=KnetArray)

    for f in fieldnames(typeof(x))
        f1=getfield(x,f)
        if typeof(f1) <: NN
            load_nn(f1,file,count,atype)
        elseif eltype(f1) <: NN
            for i=1:length(f1)
                load_nn(f1[i],file,count,atype)
            end
        else
            if f == :w
                setfield!(x,f,Param(convert(atype,read(file,string("w_",count[1])))))
                count[1]+=1
            elseif f == :b
                xx = read(file,string("b_",count[1]))
                if typeof(xx) == Float32
                    xxx = zeros(Float32,1,1,1,1)
                    xxx[1] = xx
                    setfield!(x,f,Param(convert(atype,xxx)))
                else
                    setfield!(x,f,Param(convert(atype,xx)))
                end

                count[1]+=1
            elseif f == :ms
                x.ms.momentum=read(file,string("ms_mo_",count[1]))

                xx=read(file,string("ms_mean_",count[1]))
                if typeof(xx) == Float32
                    x.ms.mean = convert(atype,[xx])
                else
                    x.ms.mean = convert(atype,xx)
                end

                xx=read(file,string("ms_var_",count[1]))
                if typeof(xx) == Float32
                    x.ms.var = convert(atype,[xx])
                else
                    x.ms.var = convert(atype,xx)
                end

                count[1]+=1
            elseif f == :bn_p
                xx = read(file,string("bn_p_",count[1]))
                if typeof(xx) == Float32
                    x.bn_p = Param(convert(atype,[xx]))
                else
                    x.bn_p = Param(convert(atype,xx))
                end
                count[1]+=1
            end
        end
    end
end

#=
Change the number of input or output features or dimensions for hourglass model
=#
function change_hourglass(hg::HG2,feature_num::Int,input_dim::Int,output_dim::Int,atype=KnetArray)

    change_hourglass_input(hg,feature_num,input_dim,atype)
    change_hourglass_output(hg,feature_num,output_dim,atype)

    nothing
end

function change_hourglass_input(hg::HG2,feature_num::Int,input_dim::Int,atype=KnetArray)

    #Input transform
    hg.fb.c1.w = Param(convert(atype,xavier_normal(Float32,7,7,input_dim,64)))
    hg.fb.c1.b = Param(convert(atype,xavier_normal(Float32,1,1,64,1)))
    hg.fb.c1.bn_p = Param(convert(atype{Float32,1},bnparams(input_dim)))
    hg.fb.c1.ms = bnmoments()

    nothing
end

function change_hourglass_output(hg::HG2,feature_num::Int,output_dim::Int,atype=KnetArray)

    for i=1:length(hg.c1)
        hg.c1[i].w = Param(convert(atype,xavier_normal(Float32,1,1,feature_num,output_dim)))
        hg.c1[i].b = Param(convert(atype,xavier_normal(Float32,1,1,output_dim,1)))
    end

    for i=1:length(hg.merge_preds)
        hg.merge_preds[i].w = Param(convert(atype,xavier_normal(Float32,1,1,output_dim,feature_num)))
        hg.merge_preds[i].b = Param(convert(atype,xavier_normal(Float32,1,1,feature_num,1)))
    end

    nothing
end


text:
using Test
using Cumulants
using Distributed
using SymmetricTensors
using CumulantsUpdates
using FileIO
using JLD2
using Random

import CumulantsUpdates: rep, cnorms

include("updatetest.jl")
include("operationstest.jl")


text:
outsize = 2 * length(data.IRF) .- 1
nffts = nextfastfft(outsize)

upad = zeros(Float64, nffts)
vpad = zeros(Float64, nffts)
convpad = zeros(Float64, nffts)

copyto!(upad, data.IRF)
pl = plan_rfft(upad)

upadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
upadRfft = pl * upad

RtDT = Array{Float64}(undef, length(data.IRF))

vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
u_vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))

#copyto!(vpad, RtDT)

function testing1(t, pl, upadRfft, vpad, RtDT, u_vpadRfft, vpadRfft)

    RtDT = DT_model(t, [0.1,10.], 1., 1., 1.)

    copyto!(vpad, RtDT)

    vpadRfft = pl * vpad

    u_vpadRfft = upadRfft .* vpadRfft

    irfft(u_vpadRfft, length(upad))
end

function conv_DT(t, β::Array{Float64,1}, data::input_data,
                 vpad, RtDT, pl, vpadRfft, upadRfft, u_vpadRfft,
                 convpad
                 )

    RtDT = DT_model(data.t, β, data.ρ, data.nmed, data.ndet)
    copyto!(vpad, RtDT)

    vpadRfft = pl * vpad
    u_vpadRfft = upadRfft .* vpadRfft
    convpad = irfft(u_vpadRfft, length(upad))
	convpad = convpad./maximum(convpad)

	tidx = findfirst(x -> x == t[1], data.t)

	convpad = convpad[tidx:tidx+length(t)-1]
	return log.(convpad)
end

function getfit(input_data, model_params; alpha=0.1)

    ## get FFT padding and plans
    outsize = 2 * length(data.IRF) .- 1
    nffts = nextfastfft(outsize)

    upad = zeros(Float64, nffts)
    vpad = zeros(Float64, nffts)
    convpad = zeros(Float64, nffts)

    copyto!(upad, data.IRF)
    pl = plan_rfft(upad)

    upadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
    upadRfft = pl * upad

    RtDT = Array{Float64}(undef, length(data.IRF))

    vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))
    u_vpadRfft = zeros(Complex{Float64}, 1 .+ div(length(upad), 2, RoundDown))

    ind1, ind2 = get_fit_window(input_data.DTOF, model_params.risefactor, model_params.tailfactor)

	weights = 1 ./ input_data.DTOFerrors

    fit = curve_fit((t, β) -> conv_DT(t, β, input_data, vpad, RtDT, pl, vpadRfft, upadRfft, u_vpadRfft, convpad), input_data.t[ind1:ind2], log.(input_data.DTOF[ind1:ind2]),
     model_params.initparams, lower=model_params.lb, upper=model_params.ub)

	perrors = estimate_errors(fit,alpha)
	marginerror = margin_error(fit,alpha)

	chisq = sum(fit.resid.^2)
	chisqdist = Distributions.Chisq(2)
	probability = Distributions.ccdf(chisqdist,chisq)
	return fitresult(input_data.t[ind1:ind2],log.(input_data.DTOF[ind1:ind2]),input_data.DTOFerrors, DT_model, fit.resid,fit.param,marginerror,perrors,2,chisq,probability)
end

function test3(input_data, model_params)

## get FFT padding and plans
outsize = 2 * length(input_data.IRF) .- 1
nffts = nextfastfft(outsize)

upad = zeros(Float64, nffts)
vpad = zeros(Float64, nffts)
convpad = zeros(Float64, nffts)

copyto!(upad, input_data.IRF)

end[[[[[[[[[]]]]]]]]]

function get_fit_window(counts, risefactor, tailfactor)

    maxval = maximum(counts)
    ind1 = findfirst(x -> x > maxval*risefactor, counts)
    ind2 =  findlast(x -> x >maxval*tailfactor, counts)

    return ind1, ind2
end

function pads(u::Array{Float64,1})
    outsize = 2 .* length(u) .- 1
    nffts = nextfastfft(outsize)
    out = zeros(eltype(u), nffts)
end

function get_fftplan(u::Array{Float64,1}, out)
    upad = zeros(eltype(u), length(out))
    copyto!(upad,u)
    pl = plan_rfft(upad)
    upadRfft = pl * upad
    return pl, upadRfft
end

out = pads(u)
pl, upadRfft = get_fftplan(u, out)

function conv_DT(t, β::Array{Float64,1},
    data::input_data,
    out,
    pl,
    upadRfft;
    vpad = zeros(eltype(data.DTOF), length(out)),
    RtDT = Array{Float64}(undef, length(data.IRF)),
    vpadRfft = zeros(eltype(upadRfft), length(upadRfft))
    )

        RtDT = DT_model(data.t, β, data.ρ, data.nmed, data.ndet)
        copyto!(vpad, RtDT)

        vpadRfft = pl * vpad
        vpadRfft .*= upadRfft
        out = irfft(vpadRfft, length(out))

        out = out./maximum(out)

	    tidx = findfirst(x -> x == t[1], data.t)

	    out = out[tidx:tidx+length(t)-1]
	return log.(out)
end

function getfit(input_data, model_params; alpha=0.1)
    ## get FFT padding and plans
    out = pads(input_data.IRF)
    pl, upadRfft = get_fftplan(input_data.IRF, out)

    #get fit windows
    ind1, ind2 = get_fit_window(input_data.DTOF, model_params.risefactor, model_params.tailfactor)

    fit = curve_fit((t, β) -> conv_DT(t, β, input_data, out, pl, upadRfft), input_data.t[ind1:ind2], log.(input_data.DTOF[ind1:ind2]),
     model_params.initparams, lower=model_params.lb, upper=model_params.ub)

end

function DT_model(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)
	#optimzed semi-infinite model

    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed

    Rt = Array{Float64}(undef, length(t))

	if n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	elseif n < 1.0
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	else
		A = 1.0
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

	z3m::Float64 = - zs
    z4m::Float64 = 2ze +zs

	Threads.@threads for n in eachindex(t)

   		Rt[n] = -exp(-(ρ^2/(4D*ν*t[n])) - μa*ν*t[n])
   		Rt[n] = Rt[n]/(2*(4π*D*ν)^(3/2)*t[n]^(5/2))

   		Rt[n] = Rt[n]*(z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n]))))

	end

   replace!(Rt, NaN => 0)
   return Rt./maximum(Rt)

end

function DT_model(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)
	#optimzed semi-infinite model

    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed

    Rt = Array{Float64}(undef, length(t))

	if n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	elseif n < 1.0
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	else
		A = 1.0
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

	z3m::Float64 = - zs
    z4m::Float64 = 2ze +zs

    Threads.@threads for n in eachindex(t)

        Rt[n] = -exp(-(ρ^2/(4D*ν*t[n])) - μa*ν*t[n])
        Rt[n] = Rt[n]/(2*(4π*D*ν)^(3/2)*t[n]^(5/2))

        Rt[n] = Rt[n]*(z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n]))))

        if Rt[n] == NaN
            Rt[n] = 0
        end

 end

     return Rt./maximum(Rt)

end

t = 0:0.01:10

ydata =  DT_model3(0:0.01:10, [0.1,10.], 1.0, 1.0, 1.0) + 0.00001*randn(length(t))

function curvetest(ydata, t)

    fit = curve_fit((t, β) -> DT_model3(t, β, 1.0, 1.0, 1.0), t, ydata, [0.5, 15.])
end

function curvetestj(ydata, t)
    fit = curve_fit((t, β) -> DT_model3(t, β, 1.0, 1.0, 1.0), (t, β) -> jacob_semiinf(t, β, 1.0, 1.0, 1.0), t, ydata, [0.5, 15.])
end

function curvetestff(ydata, t)

    fit = curve_fit((t, β) -> DT_model3(t, β, 1.0, 1.0, 1.0), t, ydata, [0.5, 15.]; autodiff=:finiteforward)
end

m(t, p) = p[1] * exp.(p[2] * t)
model_inplace(F, x, p) = (@. F = p[1] * exp.(p[2] * t))

t = 0:0.001:10
ydata = m(t, [1.0,2.0]) + 0.01*randn(length(t))

function ctest(ydata, t)
    fit = curve_fit(m, t, ydata, [0.5, 0.5])
end

function ctest_in(ydata,t)
    fit = curve_fit(model_inplace, t, ydata, [0.5, 0.5]; inplace = true)
end

function j_m(t,p)
    J = Array{Float64}(undef, length(t),length(p))
    J[:,1] = exp.(p[2] .* t)       #df/dp[1]
    J[:,2] = t .* p[1] .* J[:,1]   #df/dp[2]
    J
end

function jtest(ydata, t)

    fit = curve_fit(m, j_m, t, ydata, [0.5, 0.5])
end

function jacob_semiinf(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)

    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed

    J = Array{Float64}(undef, length(t),length(β))

	if n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	elseif n < 1.0
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	else
		A = 1.0
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

	z3m::Float64 = - zs
    z4m::Float64 = 2ze +zs

    Threads.@threads for n in eachindex(t)

        J[n,1] = ν*t[n]*exp(-(ρ^2/(4D*ν*t[n])) - μa*ν*t[n])
        J[n,1] = J[n,1]/(2*(4π*D*ν)^(3/2)*t[n]^(5/2))

        J[n,1] = J[n,1]*(z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n]))))

        if isnan(J[n,1])
            J[n,1] = 0
        end

        fx = -exp(-(ρ^2/(4D*ν*t[n])) - μa*ν*t[n])
        fx = fx/(2*(4π*D*ν)^(3/2)*t[n]^(5/2))

        dfx = 18*sqrt(3μsp)*ν*t[n]*exp((-3*ρ^2*μsp - 4μa*ν^2*t[n]^2))
        dfx = dfx - 9*sqrt(3)*μsp^(3/2)*ρ^2*exp((-3*ρ^2*μsp - 4*μa*ν^2*t[n]^2)/(4*ν*t[n]))
        dfx = dfx/(64*π^(3/2)*ν^(5/2)*t[n]^(7/2))

        gx = z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n])))

        dgx = exp(-3/(4*μsp*ν*t[n]))/μsp^2 - 3*exp(-3/(4*μsp*ν*t[n]))/(4*μsp^3*ν*t[n])
        dgx = dgx - (4A+3)*(-12*μsp*ν*t[n]*exp(-(4A + 3)^2/(12*μsp*ν*t[n])) + (4A + 3)^2*(exp(-(4A + 3)^2/(12*μsp*ν*t[n]))))/(36*μsp^3*ν*t[n])

        J[n,2] = dfx*gx + fx*dgx

        if isnan(J[n,2])
            J[n,2] = 0
        end

    end

    return J

 end

 fx[n] = -exp(-(ρ^2/(4D*ν*t[n])) - μa*ν*t[n])
 fx[n] = fx[n]/(2*(4π*ν)^(3/2)*t[n]^(5/2))

 dfx[n] = -fx[n]*(3*ρ^2)/(4ν*t[n])

 gx = 1/(D^(3/2))
 dgx = 3^(5/2)*μsp/2

 hx[n] = z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n])))

 dhx[n] = (4*ν*t[n]*μsp - 3)/(4*ν*t[n]*μsp^3)*exp(-3/(4*ν*t[n]*μsp))
 dhx[n] = dhx[n] + (4A + 3)(12*t*ν*μsp - 16*A^2 - 24*A - 9)*exp(-(4A+3)^2 / (12*t*ν*μsp))

 J[n,2] = dfx[n]*gx*hx[n] + fx[n]*dgx*hx[n] + fx[n]*gx*dhx[n]

function get_fit_window(counts, risefactor, tailfactor)

    maxval = maximum(counts)
    ind1::Int = findfirst(x -> x > maxval*risefactor, counts)
    ind2::Int =  findlast(x -> x >maxval*tailfactor, counts)

    return ind1, ind2
end

function pads(u::Array{Float64,1})
    outsize = 2 .* length(u) .- 1
    out = zeros(eltype(u), nextfastfft(outsize))
    copyto!(out, u)
end

function compute_plans(IRFpad)
    pl = plan_rfft(IRFpad)
    fftIRF = pl * IRFpad
    pli = plan_irfft(fftIRF, length(IRFpad))

    return pl, fftIRF, pli
end

function conv_DT(t, β::Array{Float64,1}, data::fitDTOF, out, fftIRF, pl, pli;
    RtDT = Array{Float64}(undef, length(data.IRF))
    )

        RtDT = TPSF_DA_semiinf_refl(data.t, β, data.ρ, data.nmed, data.ndet)
        RtDT = RtDT./maximum(RtDT)
        copyto!(out, RtDT)

        fftRtDT = pl * out
        fftRtDT .*= fftIRF
        out = pli * fftRtDT

        out = out./maximum(out)

	    tidx = findfirst(x -> x == t[1], data.t)

	    out = out[tidx:tidx+length(t)-1]
	return log.(out)
end

function _getfit(data::fitDTOF, model_params::DTOF_fitparams, TPSF_DA_semiinf_refl::Function)
     ## get FFT padding and plans
     out = pads(data.IRF)

     pl, fftIRF, pli = compute_plans(out)

    #get fit windows
    ind1, ind2 = get_fit_window(data.DTOF, model_params.risefactor, model_params.tailfactor)

    curve_fit((t, β) -> conv_DT(t, β, data, out, fftIRF, pl, pli), data.t[ind1:ind2], log.(data.DTOF[ind1:ind2]),
    model_params.initparams, lower=model_params.lb, upper=model_params.ub; autodiff=:finiteforward)

end

function _getfit(data::fitDTOF, model_params::DTOF_fitparams, TPSF_DA_slab_refl::Function)

end

function _getfit(data::fitDTOF, model_params::DTOF_fitparams, TPSF_DA_paralpip_refl::Function)

end

function getfit(data::fitDTOF, model_params::DTOF_fitparams)
    _getfit(data::fitDTOF, model_params::DTOF_fitparams, model_params.model)
end

function DT_model1(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)

    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
  	ν::Float64 = 29.9792345/nmed

    Rt = Array{Float64}(undef, length(t))

	  if n > 1.0
		  A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		  5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	  elseif n < 1.0
		  A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	  else
		  A = 1.0
	  end

	  zs::Float64 = 1/μsp
	  ze::Float64 = 2A*D

	  z3m::Float64 = - zs
    z4m::Float64 = 2ze +zs

    Threads.@threads for n in eachindex(t)

        Rt[n] = -(ρ^2/(4D*ν*t[n]))
        Rt[n] = Rt[n] - μa*ν*t[n]
        Rt[n] = -exp(Rt[n])/(2*(4π*D*ν)^(3/2)*sqrt(t[n]*t[n]*t[n]*t[n]*t[n]))

        Rt[n] = Rt[n]*(z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n]))))

        if isnan(Rt[n])
            Rt[n] = 0
        end

    end

     return Rt

end

function DT_model2(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)

    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
  	ν::Float64 = 29.9792345/nmed

    Rt = Array{Float64}(undef, length(t))

	  if n > 1.0
		  A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		  5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	  elseif n < 1.0
		  A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	  else
		  A = 1.0
	  end

	  zs::Float64 = 1/μsp
	  ze::Float64 = 2A*D

	  z3m::Float64 = - zs
    z4m::Float64 = 2ze +zs

    for n in eachindex(t)

        Rt[n] = -(ρ^2/(4D*ν*t[n]))
        Rt[n] = Rt[n] - μa*ν*t[n]

        Rt[n] = Rt[n]*(z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n]))))

    end

     return Rt

end

# steady-state solution Kienle 2 layer 1998

@. D = 1/3musp

function fluence_spatial_freq(z, s, μa, μsp, l)

    D = @. 1/3*μsp
    α = @. sqrt((D*s^2 + μa)/D)
    A = 1
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])

    ϕ = BigFloat(0.0)

    if z < l && z >= 0

        ϕ = 1/(D[1]*α[1]*cosh(α[1]*(l + zb)) + D[2]*α[2]*sinh(α[1]*(l + zb)))

        ϕ = ϕ*sinh(α[1]*(zb + z0))/(D[1]* α[1])
        ϕ = ϕ*(D[1]*α[1]*cosh(α[1]*(l - z)) + D[2]*α[2]*sinh(α[1]*(l - z)))

        if z < z0
            ϕ  = ϕ - sinh(α[1]*(z0 - z))/(D[1]*α[1])
        end

    else
        ϕ = sinh(α[1]*(zb + z0))*exp(α[2]*(l - z))
        ϕ = ϕ/(D[1]*α[1]*cosh(α[1]*(l + zb)) + D[2]*α[2]*sinh(α[1]*(l + zb)))
    end

    return ϕ
end

function fluence_steadystate(ρ, z, μa, μsp, l)

   ϕ = quadgk((s) -> compute_integral(z, s, μa, μsp, l, ρ), 0, 400)
   return ϕ[1]/2π

end

quadgk((s) -> compute_integral(z, s, μa, μsp, l, ρ), 0, 50)
# 1/2π * integral from 0 to inf ϕ*s*besselj0(ρ*s) ds

sinh(2x)*(0.1*x*cosh(2x) + 0.1*x*sinh(2x))/(0.1*x*(0.1*x*cosh(3x) + 0.1*x*sinh(3x)))

function compute_integral(z, s, μa, μsp, l, ρ)
    return fluence_spatial_freq(z, s, μa, μsp, l)*s*besselj0(ρ*s)
end

function fluence_spatial_freq(z, s, μa, μsp, l)
    D = @. 1/3*μsp
    α = @. sqrt((D*s^2 + μa)/D)
    A = 1
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])

    if isnan(fluence_spatial_freq_small(z, α, zb, z0, l, D))
        ϕ = fluence_spatial_freq_large(z, α, z0, l, D)

    else
        ϕ = fluence_spatial_freq_small(z, α, zb, z0, l, D)
    end

end

function fluence_spatial_freq_small(z, α, zb, z0, l, D)

    if z < l && z >= 0

        ϕ = sinh(α[1]*(zb + z0))/(D[1]*α[1])
        ϕ = ϕ*(D[1]*α[1]*cosh(α[1]*(l - z)) + D[2]*α[2]*sinh(α[1]*(l - z)))
        ϕ = ϕ/(D[1]*α[1]*cosh(α[1]*(l + zb)) + D[2]*α[2]*sinh(α[1]*(l + zb)))

        if z < z0
            ϕ  = ϕ - sinh(α[1]*(z0 - z))/(D[1]*α[1])
        end

    else
        ϕ = sinh(α[1]*(zb + z0))*exp(α[2]*(l - z))
        ϕ = ϕ/(D[1]*α[1]*cosh(α[1]*(l + zb)) + D[2]*α[2]*sinh(α[1]*(l + zb)))
    end

    if isinf(ϕ)
        return NaN
    end

    return ϕ
end

function fluence_spatial_freq_large(z, α, z0, l, D)

    if z < z0
        ϕ = 0.0

    elseif z < l && z >= z0
        ϕ = exp(α[1]*z0 - α[1]*z)/(2*D[1]*α[1])

    else
        ϕ = exp(α[1]*z0 + α[2]*l - α[2]*z - α[1]*l)/(D[1]*α[1] + D[2]*α[2])
    end

end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
module PkgFix

# package code goes here
include("pkg_adder.jl")

end # module


text:
"
    gather(w::AbstractMatrix{T}, xs::OneHotArray) where

getting vector at the given onehot encoding.
"
gather(w::AbstractMatrix{T}, xs::OneHotArray) where T = gather(w, onehot2indices(xs))

"
    gather(w::AbstractMatrix{T}, xs) where

getting vector at the given indices, `xs` is a array of indices(`Int` type).
"
function gather(w::AbstractMatrix{T}, xs::AbstractArray{Int}) where T
    ys = similar(w, size(w, 1), size(xs)...)

    Threads.@threads for i = 1:length(xs)
        @inbounds ind = Tuple(CartesianIndices(xs)[i])
        @inbounds ys[:, ind...] .= w[:, xs[i]]
    end
    return ys
end

"
    gather(w::AbstractArray{T}, xs) where

getting vector at the given indices, `xs` is a array of cartesian indices(`Tuple{Int}` type).
"
function gather(w::AbstractArray{T}, xs::AbstractArray{<:Tuple}) where T
    ys = similar(w, size(w, 1), size(xs)...)

    Threads.@threads for i = 1:length(xs)
        @inbounds ind = Tuple(CartesianIndices(xs)[i])
        @inbounds ys[:, ind...] .= w[:, xs[i]...]
    end
    return ys
end

# cpu gather back
function ∇gather(Δ::AbstractArray{T}, w::AbstractMatrix{T}, xs::AbstractArray{Int}) where T
    ys = fill!(similar(w), zero(T))
    scatter_add!(ys, Δ, xs)
    return ys
end

function ∇gather(Δ::AbstractArray{T}, w::AbstractArray{T}, xs::AbstractArray{<:Tuple}) where T
    ys = fill!(similar(w), zero(T))
    scatter_add!(ys, Δ, xs)
    return ys
end

using ZygoteRules: @adjoint, AContext, pullback
import ZygoteRules: _pullback

_pullback(cx::AContext, ::typeof(gather), w, xs::OneHotArray) = _pullback(cx, gather, w, onehot2indices(xs))
@adjoint gather(w, xs::AbstractArray{<: Union{Int, <:Tuple}}) = gather(w, xs), Δ->(∇gather((similar(w, size(Δ)) .= Δ), w, xs),nothing)


text:
using Sigma
using Base.Test

X = uniform(0,10)
Y = uniform(0,10)
Z = uniform(0,10)

formula = 3X + 2Y - Z >= 4
x,y,z = rand((X, Y, Z), formula)
@test 3x + 2y - z >= 4
@test 0 < x < 10
@test 0 < y < 10
@test 0 < z < 10


text:
boson = @theory begin :c * :cdag => :cdag * :c + 1 end

distr_assoc = Library.distrib(:*, :+) ∪ Library.associativity(:*)

ident = @theory begin
    1 * x => x
end

function boson_expand(ex)
    G = EGraph(ex)
    saturate!(G, boson ∪ distr_assoc, timeout=4)
    extractor = addanalysis!(G, ExtractionAnalysis, astsize_inv)
    ex = extract!(G, extractor)

    G = EGraph(ex)
    saturate!(G, distr_assoc ∪ ident, timeout=4)
    extractor = addanalysis!(G, ExtractionAnalysis, astsize)
    ex = extract!(G, extractor)
    # ex = rewrite(ex, ident)
    println(ex)
    return ex
end

e1 = :(c * c * cdag * cdag)
println(boson_expand(e1))

println(normalize_nocycle(boson_expand, e1))


text:
# Named Tuples
x = (name="abc", title="def")
println(x.name)
println(x.title)
println("$(x[1]), $(x[2])")

# Functions
f(x, y) = x + y
intro = x -> x * " is my name"
println(f(2, 3))
println(intro("John"))

using Plots
plot(x -> x^2 + 2x - 1, 0, 10, xlabel="x", ylabel="y")

function named_args(x, y; a="john", size=3, height=100)
    println("This function works!")
end

# Closure, i.e. return type is a function
foo(x) = ()->x
bar = foo(1)
println(bar())
println(foo(3)())

# Blocks
a, b, c = 1, 2, 3
let x = a, y = b, z = c
    @show x y z
end
# The following is not valid because there are no global x,y,z variables:
# println(x, y, z)

# Error occurs for the following:
# let a = 1, b
#     @show a b c
# end
# a = 1 is local, b is undefined local, c = 3 is global.

function dotest(seq)
    map(seq) do x
        if x < 0 && iseven(x)
            return 0
        elseif x == 0
            return 1
        else
            return x
        end
    end
end

println(dotest([-3,-2,-1,0,1,2,3]))
println(dotest(1:10))

function dotest2(seq1, seq2)
    map(seq1, seq2) do x, y
        x + y
    end
end

println(dotest2([1, 2], [3, 4]))
println(dotest2(1:10, 11:20))

# Ternary Operator: condition ? iftrue: else
a = 5
println(a < 10 ? println("Is < 10") : println("Is >= 10"))

# Short-Circuit Evaluation
# && and || only evaluate subsequent arguments / statements when it is required
function fact(n::Integer)
    n >= 0 || error("n must be non-negative")
    n == 0 && return 1
    n * fact(n-1)
end

# Tasks / coroutines
function fib(c::Channel)
    a = 0
    b = 1
    put!(c, a)
    while true
        put!(c, b)
        (a, b) = (b, a+b)
    end
end

fib_gen = Channel(fib);
for i in 1:10
    println(take!(fib_gen))
end

for val in Channel(fib)
    print(val, " ")
    val > 20 && break
end

# Object-like Functions
struct Point{T<:Real}
    x::T
    y::T
end

function (p::Point)(x)
    return (p.x + p.y) * x
end

p = Point(1, 2)
println(p(1.5))

arr = [1, missing, 2]
println(sum(arr))
println(typeof(missing))
println(sum(skipmissing(arr)))

function squaresum(a::Float64, b::Float64)
    a^2 + b^2
end

using InteractiveUtils

println(@code_lowered squaresum(3.0, 4.0))
println(@code_typed squaresum(3.0, 4.0))
println(@code_llvm squaresum(3.0, 4.0))
println(@code_native squaresum(3.0, 4.0))

@warn "This is an example of a warning."

# Turn on debugging for target module
# ENV["JULIA_DEBUG"] = "all"
ENV["JULIA_DEBUG"] = Main
@debug "The sum of some values $(sum(rand(100)))"


text:
function set_volume(nodes,cellxmax,cellymax)
    volume = zeros(cellxmax,cellymax)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            vec_r1x = nodes[i+1,j+1,1] - nodes[i,j,1]
            vec_r1y = nodes[i+1,j+1,2] - nodes[i,j,2]
            vec_r2x = nodes[i,j+1,1] - nodes[i+1,j,1]
            vec_r2y = nodes[i,j+1,2] - nodes[i+1,j,2]

            volume[i,j] = abs(vec_r1x*vec_r2y - vec_r1y*vec_r2x) /2
        end
    end
    return volume
end

function setup_Qbase_hat(Qbase,cellxmax,cellymax,volume)
    Qbase_hat = zeros(cellxmax,cellymax)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            Qbase_hat[i,j] = Qbase[i,j] * volume[i,j]
        end
    end
    return Qbase_hat
end

function Qhat_to_Q(Qbase_hat,cellxmax,cellymax,volume)
    Qbase = zeros(cellxmax,cellymax)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            Qbase[i,j] = Qbase_hat[i,j] / volume[i,j]
        end
    end
    return Qbase
end


text:
using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
using Plots
gr()

# Create a test problem
function simple(u, p, t)
    return [(2.0u[2]^2 - 3.0)/(1.0 + u[1]^2); -u[1]^2/(2.0 + u[2]^2); (1-u[2])/(1+u[3]^2)]
end

u0 = [2.37; 1.58; -3.10]
tspan = (0.0, 10.0)
prob = ODEProblem(simple, u0, tspan)
sol = solve(prob, Tsit5(), saveat = 0.1)
plot(sol)

# Create the differential data
X = sol[:,:]
DX = similar(X)
for (i, xi) in enumerate(eachcol(X))
    DX[:, i] = simple(xi, [], 0.0)
end

# Create a basis
@variables u[1:3]
polys = Operation[]
# Lots of basis functions
for i ∈ 0:6
    if i == 0
        push!(polys, u[1]^0)
    end
    for ui in u
        if i > 0
            push!(polys, ui^i)
        end
    end
end

basis= Basis(polys, u)
Ψ = ISINDy(X, DX, basis, ADM(1e-2), maxiter = 10, rtol = 0.1)
println(Ψ)
print_equations(Ψ)

# Lets use SINDy pi for a parallel implicit implementation
# Create a basis
@variables u[1:3]
polys = Operation[]
# Lots of basis functions -> SINDy pi can handle more than ADM()
for i ∈ 0:10
    if i == 0
        push!(polys, u[1]^0)
    end
    for ui in u
        if i > 0
            push!(polys, ui^i)
        end
    end
end

basis= Basis(polys, u)

# Simply use any optimizer you would use for SINDy
Ψ = ISINDy(X[:, :], DX[:, :], basis, STRRidge(1e-1), maxiter = 100, normalize = true)

# Transform into ODE System
sys = ODESystem(Ψ)
dudt = ODEFunction(sys)
ps = parameters(Ψ)

# Simulate
estimator = ODEProblem(dudt, u0, tspan, ps)
sol_ = solve(estimator, Tsit5(), saveat = 0.1)

# Yeah! We got it right
plot(sol.t[:], sol[:,:]', color = :red, label = nothing)
plot!(sol_.t, sol_[:, :]', color = :green, label = "Estimation")

plot(sol.t, abs.(sol-sol_)')
norm(sol[:,:]-sol_[:,:], 2) # approx 9e-7


text:
using Morton
using Test
using Random

@test cartesian2morton([5,2]) == 19
@test cartesian3morton([5,2,1]) == 67
@test morton2cartesian(19) == [5,2]
@test morton3cartesian(67) == [5,2,1]
@test tree2morton([2,1,3]) == 19
@test tree3morton([2,1,3]) == 67
@test morton2tree(19) == [2,1,3]
@test morton3tree(67) == [2,1,3]
@test tree2cartesian([2,1,3]) == [5,2]
@test tree3cartesian([2,1,3]) == [5,2,1]
@test cartesian2tree([5,2]) == [2,1,3]
@test cartesian3tree([5,2,1]) ==  [2,1,3]

Random.seed!(1234)

m = rand(0x01:0xff,10)
@test map(cartesian2morton ∘ morton2cartesian, m) == m
@test map(cartesian3morton ∘ morton3cartesian, m) == m
@test map(tree2morton ∘ morton2tree, m) == m
@test map(tree3morton ∘ morton3tree, m) == m

for i=1:10
    c = rand(0x01:0xff,2)
    @test (tree2cartesian ∘ cartesian2tree)(c) == c
    @test (morton2cartesian ∘ cartesian2morton)(c) == c
    c = rand(0x01:0xff,3)
    @test (tree3cartesian ∘ cartesian3tree)(c) == c
    @test (morton3cartesian ∘ cartesian3morton)(c) == c
end

for i=1:10
    t = rand(1:4,10)
    while t[1]==1; popfirst!(t); end
    @test (morton2tree ∘ tree2morton)(t) == t
    @test (cartesian2tree ∘ tree2cartesian)(t) == t
    t = rand(1:8,10)
    while t[1]==1; popfirst!(t); end
    @test (morton3tree ∘ tree3morton)(t) == t
    @test (cartesian3tree ∘ tree3cartesian)(t) == t
end


text:
# ------------ GFI interface implementations for Jaynes ------------ #

# Utilities.
is_choice(::Choice) = true
is_choice(a) = false
is_value(::Value) = true
is_value(a) = false

# ------------ Selections ------------ #

struct JSelection{K <: AddressMap{Select}} <: Selection
    sel::K
end
unwrap(js::JSelection) = js.sel

Base.in(addr, selection::JSelection) = haskey(selection.sel, addr)
Base.getindex(selection::JSelection, addr) = get_sub(selection.sel, addr)
Base.isempty(selection::JSelection, addr) = isempty(selection.sel, addr)

select(v::Vector{T}) where T <: Tuple = JSelection(target(v))
display(js::JSelection) = display(unwrap(js))

# ------------ Choice map ------------ #

struct JChoiceMap{K <: AddressMap} <: ChoiceMap
    chm::K
end
unwrap(jcm::JChoiceMap) = jcm.chm

@inline has_value(choices::JChoiceMap, addr) = has_value(unwrap(choices), addr)
@inline get_value(choices::JChoiceMap, addr) = getindex(unwrap(choices), addr)
@inline get_submap(choices::JChoiceMap, addr) = get_sub(unwrap(choices), addr)
@inline get_sub(choices::JChoiceMap, addr) = get_sub(unwrap(choices), addr)
@inline projection(choices::JChoiceMap, sel) = projection(unwrap(choices), sel)
@inline filter(fn, choices::JChoiceMap) = filter(fn, unwrap(choices))
@inline filter(fn, tup, choices::JChoiceMap) = filter(fn, tup, unwrap(choices))

# TODO: fix.
function get_values_shallow(choices::JChoiceMap)
    shallow_iterator(unwrap(choices))
end
function get_submaps_shallow(choices::JChoiceMap)
    shallow_iterator(unwrap(choices))
end

function merge(chm1::JChoiceMap, chm2::JChoiceMap)
    new, check = merge(unwrap(chm1), unwrap(chm2))
    JChoiceMap(new)
end

to_array(choices::JChoiceMap, ::Type{T}) where T = array(unwrap(choices), T)
from_array(choices::JChoiceMap, arr::Vector) = target(unwrap(choices), arr)
display(jcm::JChoiceMap) = Jaynes.display(unwrap(jcm))

choicemap(c::Vector{Pair{T, K}}) where {T <: Tuple, K} = JChoiceMap(target(c))

@inline collect!(par::T, addrs::Vector, chd::Dict, chm::JChoiceMap, meta) where T <: Tuple = collect!(par, addrs, chd, unwrap(chm), meta)

function convert(::Type{DynamicMap{Value}}, chm::DynamicChoiceMap)
    dm = DynamicMap{Value}()
    for (k, v) in get_values_shallow(chm)
        set_sub!(dm, k, Value(v))
    end
    for (k, v) in get_submaps_shallow(chm)
        sub = convert(DynamicMap{Value}, v)
        set_sub!(dm, k, sub)
    end
    dm
end

@inline convert(::Type{DynamicMap{Value}}, d::JChoiceMap) = convert(DynamicMap{Value}, unwrap(d))

function convert(::Type{DynamicMap{Select}}, sel::Gen.DynamicSelection)
    dm = DynamicMap{Select}()
    for (k, v) in Gen.get_subselections(sel)
        if v isa Gen.AllSelection
            set_sub!(dm, k, SelectAll())
        else
            sub = convert(DynamicMap{Select}, v)
            set_sub!(dm, k, sub)
        end
    end
    dm
end

static(chm::DynamicChoiceMap) = static(convert(DynamicMap{Value}, chm))

# ------------ Trace ------------ #

mutable struct JTrace{A, R} <: Trace
    chm::JChoiceMap
    score::Float64
    jfn::GenerativeFunction
    args::A
    ret::R
    isempty::Bool
end

@inline get_record(trace::JTrace) = trace.chm
@inline set_retval!(trace::JTrace, retval) = (trace.retval = retval)
@inline has_choice(trace::JTrace, addr) = haskey(trace.record, addr) && is_choice(get_sub(trace.record, addr))
@inline has_value(trace::JTrace, addr) = has_value(get_record(trace), addr)
@inline get_value(trace::JTrace, addr) = getindex(get_record(trace), addr)
@inline filter(fn::Function, tr::JTrace) = JChoiceMap(filter(fn, convert(DynamicMap{Value}, unwrap(get_record(tr)))))

function get_choice(trace::JTrace, addr)
    ch = get_sub(trace.record, addr)
    !is_choice(ch) && throw(KeyError(addr))
    ch
end

get_choices(trace::JTrace) = get_record(trace)

Base.display(jtr::JTrace) = Base.display(get_record(jtr))

# Trace GFI methods.
get_args(trace::JTrace) = trace.args
get_retval(trace::JTrace) = trace.ret
get_score(trace::JTrace) = trace.score
get_gen_fn(trace::JTrace) = trace.jfn

# ------------ Generative function ------------ #

struct JFunction{C <: CompilationOptions, N, R, T} <: TypedGenerativeFunction{N, R, JTrace, T}
    fn::Function
    params::DynamicMap{Value}
    params_grad::DynamicMap{Value}
    arg_types::NTuple{N, Type}
    has_argument_grads::NTuple{N, Bool}
    accepts_output_grad::Bool
    ir::IR
    reachability::FlowAnalysis
    trace_type::T
end

function JFunction(opt::J,
                   func::Function,
                   arg_types::NTuple{N, Type},
                   has_argument_grads::NTuple{N, Bool},
                   accepts_output_grad::Bool,
                   ::Type{R}) where {J <: CompilationOptions, N, R}
    (tt, ir) = instantiation_pipeline(func, arg_types, R, opt)
    return JFunction{J, N, R, typeof(tt)}(func,
                                          DynamicMap{Value}(),
                                          DynamicMap{Value}(),
                                          arg_types,
                                          has_argument_grads,
                                          accepts_output_grad,
                                          ir,
                                          flow_analysis(ir),
                                          tt)
end

@inline (jfn::JFunction)(args...) = jfn.fn(args...)
@inline has_argument_grads(jfn::JFunction) = jfn.has_argument_grads
@inline get_params(jfn::JFunction) = jfn.params
@inline get_param(jfn::JFunction, name) = getindex(jfn.params, name)
@inline set_param!(jfn::JFunction, name, v) = set_sub!(jfn.params, name, Value(v))
@inline get_params_grad(jfn::JFunction) = jfn.params_grad
@inline get_param_grad(jfn::JFunction, name) = getindex(jfn.params_grad, name)
@inline zero_param_grad!(jfn::JFunction, name) = set_sub!(jfn.params_grad, name, Value(zero(jfn.params[name])))
@inline set_param_grad!(jfn::JFunction, name::Symbol, v) = set_sub!(jfn.params_grad, name, Value(v))

init_param!(jfn, addr, val) = set_sub!(jfn.params, addr, Value(val))
init_param!(jfn, v::Vector{Pair{T, K}}) where {T <: Tuple, K} = begin
    for (addr, val) in v
        init_param!(jfn, addr, val)
    end
end

# Jaynes extensions.
@inline get_fn(jfn::JFunction) = jfn.fn
@inline get_analysis(jfn::JFunction) = jfn.reachability
@inline get_ir(jfn::JFunction) = jfn.ir
@inline get_trace_type(jfn::JFunction{C, N, R, T}) where {C, N, R, T} = T
@inline get_opt_type(::JFunction{C}) where C = C

# Typed JFunction instances have a defined notion of AC (absolute continuity).
@inline Base.:(<<)(jfn1::JFunction{N1, R1, T1}, jfn2::JFunction{N2, R2, T2}) where {N1, N2, R1, R2, T1, T2} = T1 << T2

# ------------ Model GFI interface ------------ #

function simulate(jfn::JFunction{C}, args::Tuple) where C
    ret, cl = simulate(C(),
                       get_params(jfn),
                       jfn.fn,
                       args...)
    JTrace(get_trace(cl) |> JChoiceMap, cl.score, jfn, args, ret, false)
end

function generate(jfn::JFunction{C}, args::Tuple, chm::JChoiceMap) where C
    ret, cl, w = generate(C(),
                          unwrap(chm),
                          get_params(jfn),
                          jfn.fn,
                          args...)
    JTrace(get_trace(cl) |> JChoiceMap, cl.score, jfn, args, ret, false), w
end
@inline generate(jfn::JFunction, args::Tuple, choices::DynamicChoiceMap) = generate(jfn, args, JChoiceMap(convert(DynamicMap{Value}, choices)))
@inline generate(jfn::JFunction, args::Tuple, choices::Gen.EmptyChoiceMap) = generate(jfn, args, JChoiceMap(Empty()))
@inline generate(jfn::JFunction, args::Tuple, choices::DynamicMap) = generate(jfn, args, JChoiceMap(choices))

function assess(jfn::JFunction{C}, args::Tuple, choices::JChoiceMap) where C
    ret, w = assess(C(),
                    unwrap(choices),
                    get_params(jfn),
                    jfn.fn,
                    args...)
    w, ret
end
@inline assess(jfn::JFunction, args::Tuple, choices::DynamicChoiceMap) = assess(jfn, args, JChoiceMap(convert(DynamicMap{Value}, choices)))

# TODO: must accept any inheritor of Selection.
function project(jtr::JTrace, sel::JSelection)
    w, proj = projection(get_record(jtr), unwrap(sel))
    w
end

function project(jtr::JTrace, sel::EmptySelection)
    w, proj = projection(get_record(jtr), Empty())
    w
end

function propose(jfn::JFunction, args::Tuple)
    ret, chm, w = propose(get_params(jfn),
                          jfn.fn,
                          args...)
    JChoiceMap(chm), w, ret
end

function update(trace::JTrace, args::Tuple, arg_diffs::Tuple, constraints::JChoiceMap)
    ret, cl, w, rd, d = update(get_opt_type(get_gen_fn(trace))(),
                               unwrap(constraints),
                               get_params(get_gen_fn(trace)),
                               DynamicCallSite(unwrap(get_record(trace)),
                                               get_score(trace),
                                               get_gen_fn(trace).fn,
                                               get_args(trace),
                                               get_retval(trace)),
                               map(zip(args, arg_diffs)) do (a, d)
                                   Diffed(a, d)
                               end...)
    JTrace(get_trace(cl) |> JChoiceMap, cl.score, get_gen_fn(trace), args, ret, false), w, rd, JChoiceMap(d)
end
@inline update(trace::JTrace, args::Tuple, arg_diffs::Tuple, constraints::DynamicChoiceMap) = update(trace, args, arg_diffs, JChoiceMap(convert(DynamicMap{Value}, constraints)))
@inline update(trace::JTrace, args::Tuple, arg_diffs::Tuple, constraints::StaticMap) = update(trace, args, arg_diffs, JChoiceMap(constraints))

function regenerate(trace::JTrace, args::Tuple, arg_diffs::Tuple, selection::JSelection)
    ret, cl, w, rd, d = regenerate(get_opt_type(get_gen_fn(trace))(),
                                   unwrap(selection),
                                   get_params(get_gen_fn(trace)),
                                   DynamicCallSite(unwrap(get_record(trace)),
                                                   get_score(trace),
                                                   get_gen_fn(trace).fn,
                                                   get_args(trace),
                                                   get_retval(trace)),
                                   map(zip(args, arg_diffs)) do (a, d)
                                       Diffed(a, d)
                                   end...)
    JTrace(get_trace(cl) |> JChoiceMap, cl.score, get_gen_fn(trace), args, ret, false), w, rd, JChoiceMap(d)
end
@inline regenerate(trace::JTrace, args::Tuple, arg_diffs::Tuple, selection::Gen.DynamicSelection) = regenerate(trace, args, arg_diffs, JSelection(convert(DynamicMap{Select}, selection)))
@inline regenerate(trace::JTrace, args::Tuple, arg_diffs::Tuple, selection::StaticMap) = regenerate(trace, args, arg_diffs, JSelection(selection))

# ------------ Gradients ------------ #

function choice_gradients(tr::JTrace,
                          selection::JSelection = JSelection(SelectAll()),
                          retgrad = nothing)
    vals, as, cgs = get_choice_gradients(unwrap(selection),
                                         get_params(get_gen_fn(tr)),
                                         get_record(tr),
                                         retgrad)
    as, vals, cgs
end

function get_learnable_gradients(ps::P, tr::JTrace, ret_grad...; scaler::Float64 = 1.0) where P <: AddressMap
    param_grads = Gradients()
    cl = DynamicCallSite(get_record(tr), get_score(tr), get_gen_fn(tr).fn, get_args(tr), get_retval(tr))
    arg_grads = accumulate_learnable_gradients!(target(), ps, param_grads, cl, ret_grad...; scaler = scaler)
    return arg_grads, param_grads
end

function accumulate_param_gradients!(tr::JTrace, retgrad = nothing, scale_factor = 1.0)
    gen_fn = get_gen_fn(tr)
    ps = get_params(gen_fn)
    ps isa Empty && return
    as, pgs = get_learnable_gradients(ps, tr, retgrad; scaler = scale_factor)
    accumulate!(gen_fn.params_grad, pgs)
end

# ------------ Optimization ------------ #

# Fixed step size gradient descent.
mutable struct FixedStepGradientDescentJState
    step_size::Float64
    jfn::JFunction
    param_list::Vector
end

function init_update_state(conf::Gen.FixedStepGradientDescent,
                           jfn::JFunction,
                           param_list::Vector{T}) where T <: Tuple
    FixedStepGradientDescentJState(conf.step_size, jfn, param_list)
end

function apply_update!(state::FixedStepGradientDescentJState)
    for param_name in state.param_list
        value = get_param(state.jfn, param_name)
        grad = get_param_grad(state.jfn, param_name)
        set_param!(state.jfn, param_name, value + grad * state.step_size)
        zero_param_grad!(state.jfn, param_name)
    end
end

# Gradient descent.
mutable struct GradientDescentJState
    step_size_init::Float64
    step_size_beta::Float64
    jfn::JFunction
    param_list::Vector
    t::Int
end

function init_update_state(conf::Gen.GradientDescent,
                           jfn::JFunction,
                           param_list::Vector{T}) where T <: Tuple
    GradientDescentJState(conf.step_size_init, conf.step_size_beta, jfn, param_list, 1)
end

function apply_update!(state::GradientDescentJState)
    step_size = state.step_size_init * (state.step_size_beta + 1) / (state.step_size_beta + state.t)
    for param_name in state.param_list
        value = get_param(state.jfn, param_name)
        grad = get_param_grad(state.jfn, param_name)
        set_param!(state.jfn, param_name, value + grad * step_size)
        zero_param_grad!(state.jfn, param_name)
    end
    state.t += 1
end

# ------------ Utilities ------------ #

function display(jfn::JFunction{C, N, R, T}; show_all = false) where {C, N, R, T}
    println(" ___________________________________\n")
    println("             JFunction\n")
    println(" fn : $(jfn.fn)")
    println(" arg_types : $(jfn.arg_types)")
    println(" ret_type : $(R)")
    println(" trace_type: $(T)")
    println(" has_argument_grads : $(jfn.has_argument_grads)")
    println(" accepts_output_grad : $(jfn.accepts_output_grad)")
    println("\n compilation options: $C")
    if show_all
        println(" ___________________________________\n")
        display(get_analysis(jfn))
    else
        println(" ___________________________________\n")
    end
end


text:
# Define and register smooth functions
_cos_wave(x, f, A, st, ϕ) = A*cos(2*π*f*(x-st) + ϕ)
_damped_sine_wave(x, f, A, st, ϕ, d) = exp((st-x)*d)*A*sin(2*π*f*(x-st) + ϕ)
_ramp(x, δ, st, et, h) = h/(et-st)*(_xH(x, δ, st) - _xH(x, δ, et))
_square_wave(x, δ, f, A, st) = A*2atan(sin(2π*(x-st)*f)/δ)/π
_step(x, δ, h, a) = h*(atan((x-a)/δ)/π + 0.5)
_triangular_wave(x, δ, f, A, st) = A*(1-2acos((1 - δ)sin(2π*(x-st)*f))/π)
_xH(x, δ, tₒ) = 0.5*(x-tₒ)*(1+((x-tₒ)/sqrt((x-tₒ)^2+δ^2)))

@register _cos_wave(x, f, A, st, ϕ)
@register _damped_sine_wave(x, f, A, st, ϕ, damping)
@register _ramp(x, δ, st, et, h)
@register _square_wave(x, δ, f, A, st)
@register _step(x, δ, h, a)
@register _triangular_wave(x, δ, f, A, st)

# Voltage sources
function ConstantVoltage(;name, V=1.0)
    val = V

    @named p = Pin()
    @named n = Pin()
    @parameters V
    @variables v(t)

    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           v ~ V
          ]
    ODESystem(eqs, t, [v], [V], systems=[p, n], defaults=Dict(V => val), name=name)
end

function CosineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0)
    o, A, f, st, ϕ = offset, amplitude, frequency, starttime, phase
    δ = 0.00001

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime phase
    @variables v(t)

    eqs = [
           v ~ p.v - n.v
           v ~ _cos_wave(t, f, A, st, ϕ) * _step(t, δ, 1.0, st) + offset
           0 ~ p.i + n.i
          ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime, phase), (o, A, f, st, ϕ)))
    ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime, phase], systems=[p, n], defaults=defaults, name=name)
end

function DampedSineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0, damping_coef=0.0)
    o, A, f, st, ϕ, d = offset, amplitude, frequency, starttime, phase, damping_coef
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime phase damping_coef
    @variables v(t)

    eqs = [
           v ~ p.v - n.v
           v ~ _step(t, δ, o, 0.0) + _damped_sine_wave(t, f, A, st, ϕ, d) * _step(t, δ, 1.0, st)
           0 ~ p.i + n.i
          ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime, phase, damping_coef), (o, A, f, st, ϕ, d)))
    ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime, phase, damping_coef], systems=[p, n], defaults=defaults, name=name)
end

function RampVoltage(;name, offset=0.0, starttime=0.0, endtime=1.0, height=1.0)
    o, st, et, h = offset, starttime, endtime, height
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset starttime endtime height
    @variables v(t)

    eqs = [
           v ~ p.v - n.v
           v ~ offset + _ramp(t, δ, st, et, h)
           0 ~ p.i + n.i
          ]
    defaults = Dict(zip((offset, starttime, endtime, height), (o, st, et, h)))
    ODESystem(eqs, t, [v], [offset, starttime, endtime, height], systems=[p, n], defaults=defaults, name=name)
end

function SineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0)
    o, A, f, st, ϕ = offset, amplitude, frequency, starttime, phase

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime phase
    @variables v(t)

    eqs = [
           v ~ p.v - n.v
           v ~ offset + (t > st) * (A*sin(2*π*f*(t - st) + ϕ))
           0 ~ p.i + n.i
          ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime, phase), (o, A, f, st, ϕ)))
    ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime, phase], systems=[p, n], defaults=defaults, name=name)
end

function SquareVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0)
    o, A, f, st  = offset, amplitude, frequency, starttime
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime
    @variables v(t)

    eqs = [
        v ~ p.v - n.v
        0 ~ p.i + n.i
        v ~ o + _square_wave(t, δ, f, A, st) * (t > st)
        ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime), (o, A, f, st)))
    ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime], systems=[p, n], defaults=defaults, name=name)
end

function StepVoltage(;name, offset=0.0, starttime=0.0, height=1.0)
    o, st, h = offset, starttime, height
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset starttime height
    @variables v(t)

    eqs = [
           v ~ p.v - n.v
           v ~ offset + _step(t, δ, h, st)
           0 ~ p.i + n.i
          ]
    defaults = Dict(zip((offset, starttime, height), (o, st, h)))
    ODESystem(eqs, t, [v], [offset, starttime, height], systems=[p, n], defaults=defaults, name=name)
end

function TriangularVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0)
    o, A, f, st  = offset, amplitude, frequency, starttime
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime
    @variables v(t)

    eqs = [
        v ~ p.v - n.v
        0 ~ p.i + n.i
        v ~ offset + (t>st) * _triangular_wave(t, δ, f, A, st)
    ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime), (o, A, f, st)))
    ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime], systems=[p, n], defaults=defaults, name=name)
end

# Current Sources
function ConstantCurrent(;name, I=1.0)
    val = I

    @named p = Pin()
    @named n = Pin()
    @parameters I
    @variables i(t)

    eqs = [
           0 ~ p.i + n.i
           i ~ p.i
           i ~ I
          ]
    ODESystem(eqs, t, [i], [I], systems=[p, n], defaults=Dict(I => val), name=name)
end

function CosineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0)
    o, A, f, st, ϕ = offset, amplitude, frequency, starttime, phase
    δ = 0.00001

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime phase
    @variables i(t)

    eqs = [
           i ~ _cos_wave(t, f, A, st, ϕ) * _step(t, δ, 1.0, st) + offset
           0 ~ p.i + n.i
           i ~ p.i
          ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime, phase), (o, A, f, st, ϕ)))
    ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime, phase], systems=[p, n], defaults=defaults, name=name)
end

function DampedSineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0, damping_coef=1.0)
    o, A, f, st, ϕ, d = offset, amplitude, frequency, starttime, phase, damping_coef
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime phase damping_coef
    @variables i(t)

    eqs = [
           i ~ _step(t, δ, o, 0.0) + _damped_sine_wave(t, f, A, st, ϕ, d) * _step(t, δ, 1.0, st)
           0 ~ p.i + n.i
           i ~ p.i
          ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime, phase, damping_coef), (o, A, f, st, ϕ, d)))
    ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime, phase, damping_coef], systems=[p, n], defaults=defaults, name=name)
end

function RampCurrent(;name, offset=0.0, starttime=0.0, endtime=1.0, height=1.0)
    o, st, et, h = offset, starttime, endtime, height
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset starttime endtime height
    @variables i(t)

    eqs = [
           i ~ _step(t, δ, o, 0.0) + _ramp(t, δ, st, et, h)
           0 ~ p.i + n.i
           i ~ p.i
          ]
    defaults = Dict(zip((offset, starttime, endtime, height), (o, st, et, h)))
    ODESystem(eqs, t, [i], [offset, starttime, endtime, height], systems=[p, n], defaults=defaults, name=name)
end

function SineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0)
    o, A, f, st, ϕ = offset, amplitude, frequency, starttime, phase

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime phase
    @variables i(t)

    eqs = [
           i ~ offset + (t > st) * (A*sin(2*π*f*(t - st) + ϕ))
           0 ~ p.i + n.i
           i ~ p.i
          ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime, phase), (o, A, f, st, ϕ)))
    ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime, phase], systems=[p, n], defaults=defaults, name=name)
end

function SquareCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0)
    o, A, f, st  = offset, amplitude, frequency, starttime
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime
    @variables i(t)

    eqs = [
        0 ~ p.i + n.i
        i ~ o + _square_wave(t, δ, f, A, st) * (t > st)
        i ~ p.i
        ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime), (o, A, f, st)))
    ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime], systems=[p, n], defaults=defaults, name=name)
end

function StepCurrent(;name, offset=0.0, starttime=0.0, height=1.0)
    o, st, h = offset, starttime, height
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset starttime height
    @variables i(t)

    eqs = [
           i ~ offset + _step(t, δ, h, st)
           0 ~ p.i + n.i
           i ~ p.i
          ]
    defaults = Dict(zip((offset, starttime, height), (o, st, h)))
    ODESystem(eqs, t, [i], [offset, starttime, height], systems=[p, n], defaults=defaults, name=name)
end

function TriangularCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0)
    o, A, f, st  = offset, amplitude, frequency, starttime
    δ = 0.0001

    @named p = Pin()
    @named n = Pin()
    @parameters offset amplitude frequency starttime
    @variables i(t)

    eqs = [
        0 ~ p.i + n.i
        i ~ offset + _step(t, δ, 1, st) * _triangular_wave(t, δ, f, A, st)
        i ~ p.i
    ]
    defaults = Dict(zip((offset, amplitude, frequency, starttime), (o, A, f, st)))
    ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime], systems=[p, n], defaults=defaults, name=name)
end


text:
using InverseLaplace
import InverseLaplace: talbot, gwr
using Test

include("weeks_test.jl")
include("interface_test.jl")

#### plain function interface

@test isapprox( talbot(s -> 1/s,1.0) , 1.000000000007737)
@test isapprox( talbot(s -> 1/s,1) , 1.0)

# ilt wraps talbot
@test isapprox( talbot(s -> 1/s,1), ilt(s -> 1/s,1))

@test isapprox( gwr(s -> 1/s,1) , 1.0)
@test isapprox( gwr(s -> 1/s,1.0) , 1.0000000002465594)
@test isapprox( talbot(s -> s/(1+s^2),1) , cos(1))
@test isapprox( talbot(s -> 1/(1+s),1) , exp(-1))
@test isapprox( gwr(s -> 1/(1+s),1) , exp(-1))

if Int != Int32 && VERSION >= v"0.5.0"
    @test isapprox( talbot(s -> 1/s^4,1//10) * 6 , 1e-3)
    @test isapprox( gwr(s -> 1/s^4,1//10) * 6 , 0.0010000012595365085; atol = 1e-9)

    # At one point, talbot returned rational with rational input. Now we convert to BigFloat
    @test isapprox( talbot.(s -> 1/s^3,[1//10,2//10], 64) , [5e-3, 2e-2]; atol = 1e-18)
    @test (talbot.(s -> 1/s^3,[[.1,.2] ; [.4,.5]]); true)
    @test isapprox( talbot.(s -> 1/s^3, [.1 .2; .3 .4] ),  [0.005 0.02; 0.045 0.08], atol = 1e-13)
else
    @test isapprox( talbot(s -> 1/s^4,1//10) * 6 , 1e-3; atol = 1e-9)
    @test isapprox( gwr(s -> 1/s^4,1//10) * 6 , 0.0010000012595365085; atol = 1e-9)
    @test isapprox( talbot.(s -> 1/s^3, [.1 .2; .3 .4] ),  [0.005 0.02; 0.045 0.08]; atol = 1e-9)
end

@test isapprox( talbot.(s -> 1/s^3,[.1,.2]) , [0.005,0.02])

# gwr is often less accurate.

if VERSION >= v"0.5.0"
    @test isapprox( gwr.(s -> 1/s^3,[.1,.2]) , [0.005,0.02]; atol = 1e-10)
else
    @test isapprox( gwr.(s -> 1/s^3,[.1,.2]) , [0.005,0.02]; atol = 1e-7)
end

@test isapprox( InverseLaplace.talbotarr(s -> 1/s^3,[.1,.2]) , [0.005,0.02])

@test isapprox( map(x -> convert(Float64,x) , talbot.(s -> 1/s^3,[BigFloat(1//10) , BigFloat(2//10)])), [0.005,0.02])

@test typeof(talbot(s -> 1/s^3,1//10)) == BigFloat
@test typeof(talbot(s -> 1/s^3,1//10, 64)) == BigFloat