text:
"""
    mutable struct SolidStateDetector{T <: SSDFloat, CS} <: AbstractConfig{T}

CS: Coordinate System: -> :cartesian / :cylindrical
"""
mutable struct SolidStateDetector{T <: SSDFloat, CS} <: AbstractConfig{T}
    name::String  # optional
    inputunits::Dict{String, Unitful.Units}
    world::World{T, 3}

    config_dict::Dict

    medium::NamedTuple # this should become a struct at some point

    semiconductors::Vector{Semiconductor{T}}
    contacts::Vector{Contact{T}}
    passives::Vector{Passive{T}}

    virtual_drift_volumes::Vector{AbstractVirtualVolume{T}}
end

get_precision_type(d::SolidStateDetector{T}) where {T} = T
get_coordinate_system(d::SolidStateDetector{T, CS}) where {T, CS} = CS

function SolidStateDetector{T, S}()::SolidStateDetector{T} where {T <: SSDFloat, S}
    semiconductors::Vector{Semiconductor{T}}, contacts::Vector{Contact{T}}, passives::Vector{Passive{T}} = [], [], []
    virtual_drift_volumes::Vector{AbstractVirtualVolume{T}} = []
    world_limits = get_world_limits_from_objects(Val(S), semiconductors, contacts, passives )
    world = World(Val(S), world_limits)

    return SolidStateDetector{T, S}(
        "EmptyDetector",
        default_unit_dict(),
        world,
        Dict(),
        material_properties[materials["vacuum"]],
        semiconductors,
        contacts,
        passives,
        virtual_drift_volumes
    )
end

function SolidStateDetector{T}()::SolidStateDetector{T} where {T <: SSDFloat}
    S::Symbol = :cartesian
    return SolidStateDetector{T, S}()
end
function SolidStateDetector()::SolidStateDetector{Float32, :cartesian}
    return SolidStateDetector{Float32, :cartesian}()
end

function default_unit_dict()::Dict{String, Unitful.Units}
    return Dict{String, Unitful.Units}(
        "length" => u"m", # change this to u"m" ? SI Units
        "potential" => u"V",
        "angle" => u"°",
        "temperature" => u"K"
    )
end

function construct_units(config_file_dict::Dict)
    dunits::Dict{String, Unitful.Units} = default_unit_dict()
    if haskey(config_file_dict, "units")
        d = config_file_dict["units"]
        if haskey(d, "length") dunits["length"] = unit_conversion[d["length"]] end
        if haskey(d, "angle") dunits["angle"] = unit_conversion[d["angle"]] end
        if haskey(d, "potential") dunits["potential"] = unit_conversion[d["potential"]] end
        if haskey(d, "temperature") dunits["temperature"] = unit_conversion[d["temperature"]] end
    end
    dunits
end

function construct_semiconductor(T, sc::Dict, inputunit_dict::Dict{String, Unitful.Units})
    Semiconductor{T}(sc, inputunit_dict)
end

function construct_passive(T, pass::Dict, inputunit_dict::Dict{String, Unitful.Units})
    Passive{T}(pass, inputunit_dict)
end

function construct_contact(T, contact::Dict, inputunit_dict::Dict{String, Unitful.Units})
    Contact{T}(contact, inputunit_dict)
end

function construct_virtual_volume(T, pass::Dict, inputunit_dict::Dict{String, Unitful.Units})
    construct_virtual_volume(T, pass, inputunit_dict, Val{Symbol(pass["model"])} )
end
function construct_virtual_volume(T, pass::Dict, inputunit_dict::Dict{String, Unitful.Units}, ::Type{Val{:dead}})
    DeadVolume{T}(pass, inputunit_dict)
end
function construct_virtual_volume(T, pass::Dict, inputunit_dict::Dict{String, Unitful.Units}, ::Type{Val{:arbitrary}})
    ArbitraryDriftModificationVolume{T}(pass, inputunit_dict)
end

function construct_objects(T, objects::Vector, semiconductors, contacts, passives, virtual_drift_volumes, inputunit_dict)::Nothing
    for obj in objects
        if obj["type"] == "semiconductor"
            push!(semiconductors, construct_semiconductor(T, obj, inputunit_dict))
        elseif obj["type"] == "contact"
            push!(contacts, construct_contact(T, obj, inputunit_dict))
        elseif obj["type"] == "passive"
            push!(passives, construct_passive(T, obj, inputunit_dict))
        elseif obj["type"] == "virtual_drift_volume"
            push!(virtual_drift_volumes, construct_virtual_volume(T, obj, inputunit_dict))
        else
            @warn "please specify the class to be either a \"semiconductor\", a \"contact\", or \"passive\""
        end
    end
    nothing
end

function get_world_limits_from_objects(S::Val{:cylindrical}, s::Vector{Semiconductor{T}}, c::Vector{Contact{T}}, p::Vector{Passive{T}}) where {T <: SSDFloat}
    ax1l::T, ax1r::T, ax2l::T, ax2r::T, ax3l::T, ax3r::T = 0, 1, 0, 1, 0, 1
    imps_1::Vector{T} = []
    imps_3::Vector{T} = []
    for objects in [s, c, p]
        for object in objects
            for posgeo in object.geometry_positive
                append!(imps_1, get_important_points( posgeo, Val{:r}()))
                append!(imps_3, get_important_points( posgeo, Val{:z}()))
            end
        end
    end
    imps_1 = uniq(sort(imps_1))
    imps_3 = uniq(sort(imps_3))
    if length(imps_1) > 1
        ax1l = minimum(imps_1)
        ax1r = maximum(imps_1)
    elseif length(imps_1) == 1
        ax1l = minimum(imps_1)
        ax1r = maximum(imps_1) + 1
    end
    if length(imps_3) > 1
        ax3l = minimum(imps_3)
        ax3r = maximum(imps_3)
    elseif length(imps_3) == 1
        ax3l = minimum(imps_3)
        ax3r = maximum(imps_3) + 1
    end
    return ax1l, ax1r, ax2l, ax2r, ax3l, ax3r
end
function get_world_limits_from_objects(S::Val{:cartesian}, s::Vector{Semiconductor{T}}, c::Vector{Contact{T}}, p::Vector{Passive{T}}) where {T <: SSDFloat}
    ax1l::T, ax1r::T, ax2l::T, ax2r::T, ax3l::T, ax3r::T = 0, 1, 0, 1, 0, 1
    imps_1::Vector{T} = []
    imps_2::Vector{T} = []
    imps_3::Vector{T} = []
    for objects in [s, c, p]
        for object in objects
            for posgeo in object.geometry_positive
                append!(imps_1, get_important_points( posgeo, Val{:x}()))
                append!(imps_2, get_important_points( posgeo, Val{:y}()))
                append!(imps_3, get_important_points( posgeo, Val{:z}()))
            end
        end
    end
    imps_1 = uniq(sort(imps_1))
    imps_2 = uniq(sort(imps_2))
    imps_3 = uniq(sort(imps_3))
    if length(imps_1) > 1
        ax1l = minimum(imps_1)
        ax1r = maximum(imps_1)
    elseif length(imps_1) == 1
        ax1l = minimum(imps_1)
        ax1r = maximum(imps_1) + 1
    end
    if length(imps_2) > 1
        ax2l = minimum(imps_2)
        ax2r = maximum(imps_2)
    elseif length(imps_2) == 1
        ax2l = minimum(imps_2)
        ax2r = maximum(imps_2) + 1
    end
    if length(imps_3) > 1
        ax3l = minimum(imps_3)
        ax3r = maximum(imps_3)
    elseif length(imps_3) == 1
        ax3l = minimum(imps_3)
        ax3r = maximum(imps_3) + 1
    end
    return ax1l, ax1r, ax2l, ax2r, ax3l, ax3r
end

function SolidStateDetector{T}(config_file::Dict)::SolidStateDetector{T} where{T <: SSDFloat}
    grid_type::Symbol = :cartesian
    semiconductors::Vector{Semiconductor{T}}, contacts::Vector{Contact{T}}, passives::Vector{Passive{T}} = [], [], []
    virtual_drift_volumes::Vector{AbstractVirtualVolume{T}} = []
    medium::NamedTuple = material_properties[materials["vacuum"]]
    inputunits = dunits::Dict{String, Unitful.Units} = Dict{String, Unitful.Units}(
        "length" => u"m", # change this to u"m" ? SI Units
        "potential" => u"V",
        "angle" => u"°",
        "temperature" => u"K"
    )
    inputunits = construct_units(config_file)
    if haskey(config_file, "medium")
        medium = material_properties[materials[config_file["medium"]]]
    end

    if haskey(config_file, "objects")
        construct_objects(T, config_file["objects"], semiconductors, contacts, passives, virtual_drift_volumes, inputunits)
    end

    world = if haskey(config_file, "grid")
        if isa(config_file["grid"], Dict)
            grid_type = Symbol(config_file["grid"]["coordinates"])
            World(T, config_file["grid"], inputunits)
        elseif isa(config_file["grid"], String)
            grid_type = Symbol(config_file["grid"])
            world_limits = get_world_limits_from_objects(Val(grid_type), semiconductors, contacts, passives)
            World(Val(grid_type), world_limits)
        end
    else
        world_limits = get_world_limits_from_objects(Val(grid_type), semiconductors, contacts, passives )
        World(Val(grid_type), world_limits)
    end

    c = SolidStateDetector{T, grid_type}()
    c.name = haskey(config_file, "name") ? config_file["name"] : "NoNameDetector"
    c.config_dict = config_file
    c.semiconductors = semiconductors
    c.contacts = contacts
    c.passives = passives
    c.inputunits = inputunits
    c.medium = medium
    c.world = world
    c.virtual_drift_volumes = virtual_drift_volumes
    return c
end

function SolidStateDetector(parsed_dict::Dict)
    SolidStateDetector{Float32}(parsed_dict)
end

function Base.sort!(v::AbstractVector{<:AbstractGeometry})
    hierarchies::Vector{Int} = map(x->x.hierarchy,v)
    v_result::typeof(v) = []
    for idx in sort!(unique!(hierarchies))
        push!(v_result,filter(x->x.hierarchy == hierarchies[idx],v)...)
    end
    return v_result
end

function SolidStateDetector{T}(parsed_dict::Dict) where T
    SolidStateDetector{T}(parsed_dict)
end

function contains(c::SolidStateDetector, point::AbstractCoordinatePoint{T,3})::Bool where T
    for contact in c.contacts
        if point in contact
            return true
        end
    end
    for sc in c.semiconductors
        if point in sc
            return true
        end
    end
    return false
end

function println(io::IO, d::SolidStateDetector{T, CS}) where {T <: SSDFloat, CS}
    println("________"*d.name*"________\n")
    # println("Class: ",d.class)
    println("---General Properties---")
    println("- Precision type: $(T)")
    println("- Environment Material: $(d.medium.name)")
    println("- Grid Type: $(CS)")
    println()
    println("# Semiconductors: $(length(d.semiconductors))")
    for (isc, sc)  in enumerate(d.semiconductors)
        println("\t_____Semiconductor $(isc)_____\n")
        println(sc)
    end
    println()
    println("# Contacts: $(length(d.contacts))")
    if length(d.contacts)<=5
        for c in d.contacts
            println(c)
        end
    end
    println()
    println("# Passives: $(length(d.passives))")
    if length(d.passives)<=5
        for p in d.passives
            # println(c)
        end
    end
end

function show(io::IO, d::SolidStateDetector{T}) where {T <: SSDFloat} println(io, d) end
function print(io::IO, d::SolidStateDetector{T}) where {T <: SSDFloat} println(io, d) end
function show(io::IO,::MIME"text/plain", d::SolidStateDetector) where {T <: SSDFloat}
    show(io, d)
end

# ToDo: Test it
function generate_random_startpositions(d::SolidStateDetector{T}, n::Int, Volume::NamedTuple=bounding_box(d), rng::AbstractRNG = MersenneTwister(), min_dist_from_boundary = 0.0001) where T
    delta = T(min_dist_from_boundary)
    n_filled::Int = 0
    positions = Vector{CartesianPoint{T}}(undef,n)
    while n_filled < n
        sample=CylindricalPoint{T}(rand(rng,Volume[:r_range].left:0.00001:Volume[:r_range].right),rand(rng,Volume[:φ_range].left:0.00001:Volume[:φ_range].right),rand(rng,Volume[:z_range].left:0.00001:Volume[:z_range].right))
        if !(sample in d.contacts) && contains(d,sample) && contains(d,CylindricalPoint{T}(sample.r+delta,sample.φ,sample.z))&& contains(d,CylindricalPoint{T}(sample.r-delta,sample.φ,sample.z))&& contains(d,CylindricalPoint{T}(sample.r,sample.φ,sample.z+delta))&& contains(d,CylindricalPoint{T}(sample.r,sample.φ,sample.z-delta))
            n_filled += 1
            positions[n_filled]=CartesianPoint(sample)
        end
    end
    positions
end

function paint_object(det::SolidStateDetector{T}, object::AbstractObject, grid::Grid{T, 3, S}) where {T <: SSDFloat, S}
    samples = []
    axes_syms = S == :cylindrical ? [:r, :φ, :z] : [:x, :y, :z]
    all_imps = [get_important_points(det, s) for s in axes_syms]
    for g in object.geometry_positive
        stepsizes = T[]
        for (iax, sax) in enumerate(axes_syms)
            imps_ax = all_imps[iax]
            ax::Vector{T} = grid.axes[iax].ticks
            imin::Int = searchsortednearest(ax, minimum(imps_ax))
            imax::Int = searchsortednearest(ax, maximum(imps_ax))

            ax = ax[imin:imax]
            delete_inds::Vector{Int} = Int[]
            for imp in imps_ax
                push!(delete_inds, searchsortednearest(ax, imp))
            end
            unique!(sort!(delete_inds))
            deleteat!(ax, delete_inds)
            stepsize::T = length(ax) <= 1 ? T(1) : geom_round((minimum(diff(ax)) / 4))
            imps_g = get_important_points(g, Val(sax))
            unique!(sort!(imps_g))
            if length(imps_g) > 1
                min_imps_g::T = minimum(diff(imps_g)) / 4
                if min_imps_g < stepsize
                    stepsize = geom_round(min_imps_g)
                end
            end
            imps_g_min::T, imps_g_max::T = if length(imps_g) > 0
                minimum(imps_g), maximum(imps_g)
            else
                0, 0
            end
            Δg::T = (imps_g_max - imps_g_min)
            if Δg > 0
                g_imin::Int = searchsortednearest(ax, imps_g_min)
                g_imax::Int = searchsortednearest(ax, imps_g_max)
                n_grid_points::Int = g_imax - g_imin + 1
                n::Int = Int(round(Δg / stepsize))
                if n > 2 * n_grid_points
                    stepsize = Δg / (4 * n_grid_points)
                end
            else
                stepsize = 1
            end
            if iszero(stepsize) stepsize = 1 end
            push!(stepsizes, stepsize)
        end
        append!(samples, filter( x-> x in object.geometry, sample(g, stepsizes)) )
    end
    object_gridpoints = unique!([find_closest_gridpoint(sample_point, grid) for sample_point in samples])
    return object_gridpoints
end

function paint_object(det::SolidStateDetector{T}, object::AbstractObject{T}, grid::CylindricalGrid{T}, ::Val{:φ}, φ::T )  where {T <: SSDFloat}
    closest_φ_idx=searchsortednearest(grid.axes[2].ticks, φ)
    stepsize::Vector{T}= [minimum(diff(grid.axes[1].ticks)), IntervalSets.width(grid.axes[2].interval) == 0.0 ? 0.05236 : minimum(diff(grid.axes[2].ticks)), minimum(diff(grid.axes[3].ticks))]
    stepsize /= 2
    samples = filter(x-> x in object.geometry, vcat([sample(g, stepsize) for g in object.geometry_positive]...))
    object_gridpoints = unique!([find_closest_gridpoint(sample_point,grid) for sample_point in samples])
    return filter(x -> x[2]==closest_φ_idx, object_gridpoints)
end
function paint_object(det::SolidStateDetector{T}, object::AbstractObject{T}, grid::CylindricalGrid{T}, ::Val{:r}, r::T )  where {T <: SSDFloat}
    return CartesianPoint{T}[]
end
function paint_object(det::SolidStateDetector{T}, object::AbstractObject{T}, grid::Grid{T}, ::Val{:z}, z::T )  where {T <: SSDFloat}
    return CartesianPoint{T}[]
end
function paint_object(det::SolidStateDetector{T}, object::AbstractObject{T}, grid::CartesianGrid{T}, ::Val{:x}, x::T )  where {T <: SSDFloat}
    return CartesianPoint{T}[]
end
function paint_object(det::SolidStateDetector{T}, object::AbstractObject{T}, grid::CartesianGrid{T}, ::Val{:y}, y::T )  where {T <: SSDFloat}
    return CartesianPoint{T}[]
end


text:
@testset "Comparison" begin
    #@test SymPoint([1, 2]) ≈ SymPoint([-1, -2])
    #@test !(SymPoint([1, 2]) ≈ SymPoint([-2, -1]))
    @test Polyhedra._lt(1, 1.5)
    @test !Polyhedra._lt(2, 1.5)
end

ext: jl


text:
using Knet
module CopySeq
using Main, Knet, ArgParse
using Knet: copysync!
@useifgpu CUDArt
@useifgpu CUSPARSE

function main(args=ARGS)
    info("Learning to copy sequences to test the S2S model.")
    s = ArgParseSettings()
    @add_arg_table s begin
        ("datafiles"; nargs='+'; required=true; help="First file used for training")
        ("--dictfile"; help="Dictionary file, first datafile used if not specified")
        ("--epochs"; arg_type=Int; default=1)
        ("--hidden"; arg_type=Int; default=100)
        ("--batchsize"; arg_type=Int; default=128)
        ("--lossreport"; arg_type=Int; default=0)
        ("--gclip"; arg_type=Float64; default=10.0)
        ("--lr"; arg_type=Float64; default=2.0)
        ("--fbias"; arg_type=Float64; default=0.0)
        ("--ftype"; default="Float32")
        ("--winit"; default="Gaussian(0,0.01)")
        ("--dense"; action=:store_true)
        ("--fast"; help="skip norm and loss calculations."; action=:store_true)
        ("--gcheck"; arg_type=Int; default=0)
        ("--seed"; arg_type=Int; default=42)
        ("--nosharing"; action = :store_true)
    end
    isa(args, AbstractString) && (args=split(args))
    opts = parse_args(args, s)
    println(opts)
    for (k,v) in opts; @eval ($(symbol(k))=$v); end
    seed > 0 && setseed(seed)
    dict = (dictfile == nothing ? datafiles[1] : dictfile)
    global data = Any[]
    for f in datafiles
        push!(data, S2SData(f; batchsize=batchsize, ftype=eval(parse(ftype)), dense=dense, dict=dict))
    end
    vocab = maxtoken(data[1],2)
    # global model = S2S(lstm; fbias=fbias, hidden=hidden, vocab=vocab, winit=eval(parse(winit)))
    global model = compile(:copyseq; fbias=fbias, out=hidden, vocab=vocab, winit=eval(parse(winit)))
    setp(model; lr=lr)
    if nosharing
        set!(model, :forwoverwrite, false)
        set!(model, :backoverwrite, false)
    end

    perp = zeros(length(data))
    (maxnorm,losscnt) = fast ? (nothing,nothing) : (zeros(2),zeros(2))
    t0 = time_ns()
    println("epoch  secs    ptrain  ptest.. wnorm  gnorm")
    myprint(a...)=(for x in a; @printf("%-6g ",x); end; println(); flush(STDOUT))
    for epoch=1:epochs
        fast || (fill!(maxnorm,0); fill!(losscnt,0))
        train(model, data[1], softloss; gclip=gclip, maxnorm=maxnorm, losscnt=losscnt, lossreport=lossreport)
        fast || (perp[1] = exp(losscnt[1]/losscnt[2]))
        for d=2:length(data)
            loss = test(model, data[d], softloss)
            perp[d] = exp(loss)
        end
        myprint(epoch, (time_ns()-t0)/1e9, perp..., (fast ? [] : maxnorm)...)
        gcheck > 0 && gradcheck(model,
                                f->(train(f,data[1],softloss;losscnt=fill!(losscnt,0),gcheck=true);losscnt[1]),
                                f->(test(f,data[1],softloss;losscnt=fill!(losscnt,0),gcheck=true);losscnt[1]);
                                gcheck=gcheck)
    end
    return (fast ? (perp...) :  (perp..., maxnorm...))
end

# This copies lstm exactly for replicatability:
@knet function copyseq(word; fbias=0, vocab=0, o...)
    if !decoding
        x = wdot(word; o...)
        input  = wbf2(x,h; o..., f=:sigm)
        forget = wbf2(x,h; o..., f=:sigm, binit=Constant(fbias))
        output = wbf2(x,h; o..., f=:sigm)
        newmem = wbf2(x,h; o..., f=:tanh)
    else
        x = wdot(word; o...)
        input  = wbf2(x,h; o..., f=:sigm)
        forget = wbf2(x,h; o..., f=:sigm, binit=Constant(fbias))
        output = wbf2(x,h; o..., f=:sigm)
        newmem = wbf2(x,h; o..., f=:tanh)
    end
    cell = input .* newmem + cell .* forget
    h  = tanh(cell) .* output
    if decoding
        tvec = wdot(h; out=vocab)
        return soft(tvec)
    end
end

@knet function copyseq1(word; fbias=0, vocab=0, o...)
    if decoding
        x = wdot(word; o...)
        input  = sigm(aff2(x,h; o...))
        forget = sigm(aff2(x,h; o..., binit=Constant(fbias)))
        output = sigm(aff2(x,h; o...))
        newmem = tanh(aff2(x,h; o...))
    else
        x = wdot(word; o...)
        input  = sigm(aff2(x,h; o...))
        forget = sigm(aff2(x,h; o..., binit=Constant(fbias)))
        output = sigm(aff2(x,h; o...))
        newmem = tanh(aff2(x,h; o...))
    end
    c = input .* newmem + c .* forget
    h  = tanh(c) .* output
    if decoding
        return soft(wdot(h; out=vocab))
    end
end

@knet function aff2(x,y; out=0, o...)
    a = par(; o..., dims=(out,0))
    b = par(; o..., dims=(out,0))
    c = par(; o..., dims=(0,))
    # return a*x+b*y+c
    xy = a*x+b*y
    return c+xy
end

function train(m, data, loss; o...)
    s2s_loop(m, data, loss; trn=true, ystack=Any[], o...)
end

function test(m, data, loss; losscnt=zeros(2), o...)
    s2s_loop(m, data, loss; losscnt=losscnt, o...)
    losscnt[1]/losscnt[2]
end

# Persistent storage for ygold and mask
s2s_ygold = nothing
s2s_mask = nothing
@gpu copytogpu(y,x::Array)=CudaArray(x)
@gpu copytogpu(y,x::SparseMatrixCSC)=CudaSparseMatrixCSC(x)
@gpu copytogpu{T}(y::CudaArray{T},x::Array{T})=(size(x)==size(y) ? copysync!(y,x) : copytogpu(nothing,x))
@gpu copytogpu{T}(y::CudaSparseMatrixCSC{T},x::SparseMatrixCSC{T})=(size(x)==size(y) ? copysync!(y,x) : copytogpu(nothing,x))

function s2s_loop(m, data, loss; gcheck=false, o...)
    global s2s_ygold, s2s_mask
    s2s_lossreport()
    decoding = false
    reset!(m)
    for (x,ygold,mask) in data
        nwords = (mask == nothing ? size(x,2) : sum(mask))
        # x,ygold,mask are cpu arrays; x gets copied to gpu by forw; we should do the other two here
        if ygold != nothing && gpu()
            ygold = s2s_ygold = copytogpu(s2s_ygold,ygold)
            mask != nothing && (mask = s2s_mask  = copytogpu(s2s_mask,mask)) # mask not used when ygold=nothing
        end
        if decoding && ygold == nothing # the next sentence started
            gcheck && break
            s2s_eos(m, data, loss; gcheck=gcheck, o...)
            reset!(m)
            decoding = false
        end
        if !decoding && ygold != nothing # source ended, target sequence started
            # s2s_copyforw!(m)
            decoding = true
        end
        if decoding && ygold != nothing # keep decoding target
            s2s_decode(m, x, ygold, mask, nwords, loss; o...)
        end
        if !decoding && ygold == nothing # keep encoding source
            s2s_encode(m, x; o...)
        end
    end
    s2s_eos(m, data, loss; gcheck=gcheck, o...)
end

function s2s_encode(m, x; trn=false, o...)
    # forw(m.encoder, x; trn=trn, seq=true, o...)
    (trn?sforw:forw)(m, x; decoding=false)
end

function s2s_decode(m, x, ygold, mask, nwords, loss; trn=false, ystack=nothing, losscnt=nothing, o...)
    # ypred = forw(m.decoder, x; trn=trn, seq=true, o...)
    ypred = (trn?sforw:forw)(m, x; decoding=true)
    ystack != nothing  && push!(ystack, (copy(ygold),copy(mask))) # TODO: get rid of alloc
    losscnt != nothing && s2s_loss(m, ypred, ygold, mask, nwords, loss; losscnt=losscnt, o...)
end

function s2s_loss(m, ypred, ygold, mask, nwords, loss; losscnt=nothing, lossreport=0, o...)
    (yrows, ycols) = size2(ygold)  # TODO: loss should handle mask, currently only softloss does.
    losscnt[1] += loss(ypred,ygold;mask=mask) # loss divides total loss by minibatch size ycols.  at the end the total loss will be equal to
    losscnt[2] += nwords/ycols                # losscnt[1]*ycols.  losscnt[1]/losscnt[2] will equal totalloss/totalwords.
    # we could scale losscnt with ycols so losscnt[1] is total loss and losscnt[2] is total words, but I think that breaks gradcheck since the scaled versions are what gets used for parameter updates in order to prevent batch size from effecting step size.
    lossreport > 0 && s2s_lossreport(losscnt,ycols,lossreport)
end

function s2s_eos(m, data, loss; trn=false, gcheck=false, ystack=nothing, maxnorm=nothing, gclip=0, o...)
    if trn
        s2s_bptt(m, ystack, loss; o...)
        g = (gclip > 0 || maxnorm!=nothing ? gnorm(m) : 0)
        if !gcheck
            gscale=(g > gclip > 0 ? gclip/g : 1)
            update!(m; gscale=gscale, o...)
        end
    end
    if maxnorm != nothing
        w=wnorm(m)
        w > maxnorm[1] && (maxnorm[1]=w)
        g > maxnorm[2] && (maxnorm[2]=g)
    end
end

function s2s_bptt(m, ystack, loss; o...)
    while !isempty(ystack)
        (ygold,mask) = pop!(ystack)
        # back(m.decoder, ygold, loss; seq=true, mask=mask, o...)
        sback(m, ygold, loss; mask=mask, o...) # back passes mask on to loss
    end
    # @assert m.decoder.sp == 0
    # s2s_copyback!(m)
    # while m.encoder.sp > 0
    while m.sp > 0
        # back(m.encoder; seq=true, o...)
        sback(m)                # TODO: what about mask here?
        # error(:ok)
    end
end

s2s_time0 = s2s_time1 = s2s_inst = 0

function s2s_lossreport()
    global s2s_time0, s2s_time1, s2s_inst
    s2s_inst = 0
    s2s_time0 = s2s_time1 = time_ns()
    # println("time inst speed perp")
end

s2s_print(a...)=(for x in a; @printf("%.2f ",x); end; println(); flush(STDOUT))

function s2s_lossreport(losscnt,batchsize,lossreport)
    global s2s_time0, s2s_time1, s2s_inst
    s2s_time0 == 0 && s2s_lossreport()
    losscnt == nothing && return
    losscnt[2]*batchsize < lossreport && return
    curr_time = time_ns()
    batch_time = Int(curr_time - s2s_time1)/10^9
    total_time = Int(curr_time - s2s_time0)/10^9
    s2s_time1 = curr_time
    batch_inst = losscnt[2]*batchsize
    batch_loss = losscnt[1]*batchsize
    s2s_inst += batch_inst
    s2s_print(total_time, s2s_inst, batch_inst/batch_time, exp(batch_loss/batch_inst))
    losscnt[1] = losscnt[2] = 0
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

end # module

### DEAD CODE

    # info("Warm-up epoch")
    # f=datafiles[1]; mini = S2SData(f, f; batch=batchsize, ftype=eval(parse(ftype)), dense=dense, dict1=dict[1], dict2=dict[2], stop=3200)
    # @date train(model, mini, softloss; gcheck=gcheck, gclip=gclip, getnorm=getnorm, getloss=getloss) # pretrain to compile for timing
    # info("Starting profile")


text:
module MyExample

using ForwardDiff

include("extra_file.jl")

export my_f, derivative_of_my_f

end


text:
"""
## DataUtils
`DataUtils` provides the following new constructors for `Dict`:

    Dict(keys, values)

One can provide `Dict` with (equal length) vector arguments.  The first
vector provides a list of keys, while the second provides a list of values.
If the vectors constain nulls, only key, value pairs with *both* their
elements non-null will be added.
"""
function Dict{K, V}(keys::Vector{K}, values::Vector{V})::Dict
    @assert length(keys) == length(values) ("Vectors for constructing
                                             Dict must be of equal length.")
    Dict(k=>v for (k, v) ∈ zip(keys, values))
end

function Dict(keys::Vector{Union{K,Null}}, values::Vector{Union{V,Null}})::Dict where {K,V}
    @assert length(keys) == length(values) ("Vectors for constructing
                                             Dict must be of equal length.")
    dict = Dict{K, V}()
    size_ = sum(!isnull(k) && !isnull(v) ? 1 : 0 for (k, v) ∈ zip(keys, values))::Integer
    sizehint!(dict, size_)
    # we only insert pairs if both values are not null
    for (k, v) ∈ zip(keys, values)
        if !isnull(k) && !isnull(v)
            dict[get(k)] = get(v)
        end
    end
    return dict
end

"""
    serialize(filename::AbstractString, object)

Serializes an object and stores on the local file system to file `filename`.
This is quite similar to the functionality in `Base`, except that the default
`serialize` method requires an `IOStream` object instead of a file name, so
this eliminates an extra line of code.
"""
function serialize(filename::AbstractString, object)
    f = open(filename, "w")
    serialize(f, object)
    close(f)
    return nothing
end
export serialize

"""
    deserialize(filename::AbstractString)

Opens a file from the local file system and deserializes what it finds there.
This is quite similar to the functionality in `Base` except that the default
`deserialize` method requires an `IOStream` object instead of a file name so this
eliminates an extra line of code.
"""
function deserialize(filename::AbstractString)
    f = open(filename)
    o = deserialize(f)
    close(f)
    return o
end
export deserialize

"""
    getNormedHistogramData(X)

Very annoyingly Gadfly does not yet support normed histograms.

This function returns an ordered pair of vectors which can be
fed to gadfly to create a normed histogram.  If the output
is `m, w` one can do
`plot(x=m, y=w, Geom.bar)` to create a histogram.
"""
function getNormedHistogramData{T <: Real}(X::Vector{T};
                                           nbins::Integer=StatsBase.sturges(length(X)))
    h = fit(Histogram, X, nbins=nbins)
    edges = h.edges[1]
    midpoints = Vector{Float64}(length(edges)-1)
    for i in 1:length(midpoints)
        midpoints[i] = (edges[i+1] + edges[i])/2.0
    end
    width = midpoints[2] - midpoints[1]
    weights = h.weights./(sum(h.weights)*width)
    midpoints, weights
end
export getNormedHistogramData

"""
    infast(x, collection)

Checks whether the object `x` is in `collection`. This is done efficiently by creating
hashes for the objects in `collection`.  This should only be used if `collection` is
large, as there is overhead in hashing and allocating.
"""
function infast{T}(x::T, collection::Vector{T})
    # not sure if this is the most efficient way to create the hash
    dict = Dict{T, Void}(c=>nothing for c ∈ collection)
    x ∈ keys(dict)
end
export infast

# TODO these really should be kept in another package
"""
    outer(A, B)

Performs the outer product of two tensors A_{i₁…iₙ}B_{j₁…jₙ}.

**TODO** Currently only implemented for A and B as vectors.
"""
function outer{T}(A::Vector{T}, B::Vector{T})
    C = Matrix{T}(length(A), length(B))
    for j ∈ 1:size(C, 2), i ∈ 1:size(C, 1)
        C[i, j] = A[i]*B[j]
    end
    C
end

# there's no way around iterating over every element if A and B are dense
function outer{T<:AbstractSparseArray}(::Type{T}, A::Vector, B::Vector)
    C = spzeros(length(A), length(B))  # type?
    for j ∈ 1:size(C, 2), i ∈ 1:size(C, 1)
        C[i, j] = A[i]*B[j] # note that Julia keeps this element sparse if prod is 0.0
    end
    convert(T, C)
end
export outer

"""
    _info_expr(message, code)

Private method used by `@info`.
"""
function _info_expr(message::String, code::Union{Expr,Symbol})
    quote
        info(string("Executing ", $message, " ..."))
        $code
        info(string("Done executing ", $message, " ."))
    end
end

"""
    _infotime_expr(message, code)

Private method used by `@infotime`.
"""
function _infotime_expr(message::String, code::Union{Expr,Symbol})
    quote
        info(string("Executing ", $message, " ..."))
        @time $code
        info(string("Done executing ", $message, " ."))
    end
end

"""
    @info code

Executes code sandwhiched between informative info messages telling the user that the
code is being executed.
"""
macro info(message::String, code)
    esc(_info_expr(message, code))
end
macro info(code)
    message = string('`', code, '`')
    esc(_info_expr(message, code))
end
export @info

"""
    @infotime code

Executes code sandwhiched between informative info messages telling the user that the code
is being executed, while applying the `@time` macro to the code.
"""
macro infotime(message::String, code)
    esc(_infotime_expr(message, code))
end
macro infotime(code)
    message = string('`', code, '`')
    esc(_infotime_expr(message, code))
end
export @infotime

#=========================================================================================
This section is all for breaking a matrix into submatrices based on a class label.
=========================================================================================#
getDefaultComparator{T}(::Type{T}) = (==)
getDefaultComparator{T<:AbstractFloat}(::Type{T}) = (≈)
export getDefaultComparator

"""
    findboundaries(X[, ncol; comparator=DEFAULT, check_sort=true])

Finds the boundaries between distinct values in the sorted array `X`.  For example, if
`X = [1, 1, 1, 2, 2, 3, 4]`, `boundaries(X) == [1, 4, 6, 7]`.  This algorithm requires
the input array to be sorted in the column in which boundaries are determined.
Do `sort_array=false` in cases where the input array is already sorted

This algorithm only works properly if the input array is already sorted.  By default,
this function checks if the input array is sorted.  To skip this check, set
`check_sort=false`.
"""
function findBoundaries{T}(X::Vector{T}; check_sort::Bool=true,
                           comparator=getDefaultComparator(T))
    if check_sort && !issorted(X)
        throw(ArgumentError("Trying to find boundaries on unsorted array."))
    end
    boundaries = Int64[1]
    for i ∈ 2:length(X)
        !comparator(X[i-1], X[i]) && push!(boundaries, i)
    end
    boundaries
end

function findBoundaries{T}(X::Matrix{T}, ncol::Integer; check_sort::Bool=true,
                           comparator=getDefaultComparator(T))
    findBoundaries(X[:,ncol], check_sort=check_sort, comparator=comparator)
end

"""
    findboundaries!(v[; comparator=DEFAULT])

This is the same as `findboundaries` except that it sorts the vector `v` first.
See documentation for `findboundaries`.
"""
function findBoundaries!{T}(X::Vector{T}; comparator=getDefaultComparator(T))
    if sort_array
        X = sort!(X)
    end
    findBoundaries(X, comparator=comparator)
end

"""
    boundaryClassValues(X[, ncol,] boundaries)

Given `boundaries` between distinct values in column `ncol` of `Matrix` `X` (can also be a
`Vector`) (see documentation for `findBoundaries`), find the values of the distinct values.

This is essentially the same functionality as unique, except that it is faster since
boundaries have already been determined.
"""
function boundaryClassValues{T}(X::Matrix{T}, ncol::Integer, boundaries::Vector)
    [X[boundaries[i], ncol] for i ∈ 1:length(boundaries)]
end
function boundaryClassValues{T}(v::Vector{T}, boundaries::Vector)
    [X[boundaries[i]] for i ∈ 1:length(boundaries)]
end

"""
    findBoundariesVlaues(X[, ncol; check_sort=true, comparator=DEFAULT])

Returns `boundaries, classes` where `boundaries` is the output of `findBoundaries` and
`classes` is the output of `boundaryClassValues`.  Note that `X` should be sorted before
inputting it into this function.
"""
function findBoundariesValues{T}(X::Vector{T}; check_sort::Bool=true,
                                 comparator=getDefaultComparator(T))
    boundaries = findBoundaries(X, check_sort=check_sort, comparator=comparator)
    classes = boundaryClassValues(X, boundaries)
    boundaries, classes
end

function findBoundariesValues{T}(X::Matrix{T}, ncol::Integer; check_sort::Bool=true,
                                 comparator=getDefaultComparator(T))
    boundaries = findBoundaries(X, ncol, check_sort=check_sort, comparator=comparator)
    classes = boundaryClassValues(X, ncol, boundaries)
    boundaries, classes
end

"""
    findBoundaryDict(X[, ncol; check_sort=true, comparator=DEFAULT])

For a sorted array `X`, find the boundaries between distinct values of column number
`ncol`.
"""
function findBoundaryDict{T}(X::Vector{T}; check_sort::Bool=true,
                             comparator=getDefaultComparator(T))
    boundaries, classes = findBoundariesValues(X, check_sort=check_sort,
                                               comparator=comparator)
    dict = Dict{T,Tuple{Int64,Int64}}();  sizehint!(dict, length(boundaries))
    for i ∈ 1:length(boundaries)
        dict[classes[i]] = (boundaries[i], get(boundaries, i+1, length(X)+1)-1)
    end
    dict
end

function findBoundaryDict{T}(X::Matrix{T}, ncol::Integer; check_sort::Bool=true,
                             comparator=getDefaultComparator(T))
    boundaries, classes = findBoundariesValues(X, ncol, check_sort=check_sort,
                                               comparator=comparator)
    dict = Dict{T,Tuple{Int64,Int64}}();  sizehint!(dict, length(boundaries))
    for i ∈ 1:length(boundaries)
        dict[classes[i]] = (boundaries[i], get(boundaries, i+1, size(X,1)+1)-1)
    end
    dict
end
export findBoundaryDict

_order_y(y::Vector, order::Vector) = y[order]
_order_y(y::Matrix, order::Vector) = y[order, :]

_get_y_interval(y::Vector, r::UnitRange) = y[r]
_get_y_interval(y::Matrix, r::UnitRange) = y[r, :]

"""
    subMatricesByClass(X[, y], ncol)

Breaks the arrays `X` and `y` (optional) into dictionaries of `class=>submatrix`
pairs where `class` is one of the distinct values of the column `ncol` of `X` and
`submatrix` is the range of rows of `X` or `y` where `X[:,ncol]` has the value `class`.

This function is intended for breaking up training and test sets to be used with sets of
different models.
"""
function subMatricesByClass{T,U}(X::Matrix{T}, y::U, ncol::Integer)
    if size(X,1) ≠ size(y,1)
        throw(ArgumentError("X and y must have same number of rows."))
    end
    order = sortperm(X[:, ncol])
    X = X[order, :]
    y = _order_y(y, order)
    bdict = findBoundaryDict(X, ncol, check_sort=false)
    Xdict = Dict{T,Matrix{T}}();  sizehint!(Xdict, length(bdict))
    ydict = Dict{T,U}();  sizehint!(ydict, length(bdict))
    for (k, v) ∈ bdict
        Xdict[k] = X[v[1]:v[2], :]
        ydict[k] = _get_y_interval(y, v[1]:v[2])
    end
    Xdict, ydict
end

function subMatricesByClass{T}(X::Matrix{T}, ncol::Integer)
    order = sortperm(X[:, ncol])
    X = X[order, :]
    bdict = findBoundaryDict(X, ncol, check_sort=false)
    Xdict = Dict{T,Matrix{T}}();  sizehint!(Xdict, length(bdict))
    for (k, v) ∈ bdict
        Xdict[k] = X[v[1]:v[2], :]
    end
    Xdict
end
export subMatricesByClass

"""
    invert(dict)

Inverts a dictionary or other `Associative`.  Output is always `Dict`.
"""
invert(dict::Associative) = Dict(v=>k for (k,v) ∈ dict)

# helper function for freqdict
function _check_insert(dict::Dict, a)
    if a ∈ keys(dict)
        dict[a] += 1
    else
        dict[a] = 1
    end
end
function _check_insert(dict::Dict, ::Null)
    if any(isnull(d) for d ∈ keys(dict))
        dict[null] += 1
    else
        dict[null] = 1
    end
end

"""
    freqdict(A::AbstractVector)
    freqdict(A::AbstractVector...)

Returns a `Dict` the keys of which are the distinct elements of `A` and the
values of which are the number of times that value occurs in `A`.

If multiple vectors are passed, a `Dict` will be created where the keys are
tuples containing a combination of elements from the vectors that were passed,
with values equal to the number of times that each combination occurred.
"""
function freqdict(A::AbstractVector{T}) where T
    dict = Dict{T,Int}()
    for a ∈ A
        _check_insert(dict, a)
    end
    dict
end
function freqdict(As, dtypes)
    dict = Dict{Tuple{dtypes...},Int}()
    for i ∈ 1:length(As[1])
        _check_insert(dict, tuple((a[i] for a ∈ As)...))
    end
    dict
end
freqdict(As::AbstractVector...) = freqdict(As, eltype.(As))
export freqdict


text:
"""
DuckDB type
"""
mutable struct LogicalType
    handle::duckdb_logical_type

    function LogicalType(type::DUCKDB_TYPE)
        handle = duckdb_create_logical_type(type)
        result = new(handle)
        finalizer(_destroy_type, result)
        return result
    end
    function LogicalType(handle::duckdb_logical_type)
        result = new(handle)
        finalizer(_destroy_type, result)
        return result
    end
end

function _destroy_type(type::LogicalType)
    if type.handle != C_NULL
        duckdb_destroy_logical_type(type.handle)
    end
    type.handle = C_NULL
    return
end

create_logical_type(::Type{T}) where {T <: String} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_VARCHAR)
create_logical_type(::Type{T}) where {T <: Bool} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_BOOLEAN)
create_logical_type(::Type{T}) where {T <: Int8} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_TINYINT)
create_logical_type(::Type{T}) where {T <: Int16} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_SMALLINT)
create_logical_type(::Type{T}) where {T <: Int32} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_INTEGER)
create_logical_type(::Type{T}) where {T <: Int64} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_BIGINT)
create_logical_type(::Type{T}) where {T <: UInt8} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_UTINYINT)
create_logical_type(::Type{T}) where {T <: UInt16} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_USMALLINT)
create_logical_type(::Type{T}) where {T <: UInt32} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_UINTEGER)
create_logical_type(::Type{T}) where {T <: UInt64} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_UBIGINT)
create_logical_type(::Type{T}) where {T <: Float32} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_FLOAT)
create_logical_type(::Type{T}) where {T <: Float64} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_DOUBLE)
create_logical_type(::Type{T}) where {T <: Date} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_DATE)
create_logical_type(::Type{T}) where {T <: Time} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_TIME)
create_logical_type(::Type{T}) where {T <: DateTime} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_TIMESTAMP)
create_logical_type(::Type{T}) where {T <: AbstractString} = DuckDB.LogicalType(DuckDB.DUCKDB_TYPE_VARCHAR)
function create_logical_type(::Type{T}) where {T <: FixedDecimal}
    int_type = T.parameters[1]
    width = 0
    scale = T.parameters[2]
    if int_type == Int16
        width = 4
    elseif int_type == Int32
        width = 9
    elseif int_type == Int64
        width = 18
    elseif int_type == Int128
        width = 38
    else
        throw(NotImplementedException("Unsupported internal type for decimal"))
    end
    return DuckDB.LogicalType(duckdb_create_decimal_type(width, scale))
end

function create_logical_type(::Type{T}) where {T}
    throw(NotImplementedException("Unsupported type for create_logical_type"))
end

function get_type_id(type::LogicalType)
    return duckdb_get_type_id(type.handle)
end

function get_internal_type_id(type::LogicalType)
    type_id = get_type_id(type)
    if type_id == DUCKDB_TYPE_DECIMAL
        type_id = duckdb_decimal_internal_type(type.handle)
    elseif type_id == DUCKDB_TYPE_ENUM
        type_id = duckdb_enum_internal_type(type.handle)
    end
    return type_id
end

function get_decimal_scale(type::LogicalType)
    return duckdb_decimal_scale(type.handle)
end

function get_enum_dictionary(type::LogicalType)
    dict::Vector{String} = Vector{String}()
    dict_size = duckdb_enum_dictionary_size(type.handle)
    for i in 1:dict_size
        val = duckdb_enum_dictionary_value(type.handle, i)
        str_val = String(unsafe_string(val))
        push!(dict, str_val)
        duckdb_free(val)
    end
    return dict
end

function get_list_child_type(type::LogicalType)
    return LogicalType(duckdb_list_type_child_type(type.handle))
end

function get_struct_child_count(type::LogicalType)
    return duckdb_struct_type_child_count(type.handle)
end

function get_struct_child_name(type::LogicalType, index::UInt64)
    val = duckdb_struct_type_child_name(type.handle, index)
    result = unsafe_string(val)
    duckdb_free(val)
    return result
end

function get_struct_child_type(type::LogicalType, index::UInt64)
    return LogicalType(duckdb_struct_type_child_type(type.handle, index))
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
"""
For this code to work the 3-dim representation should
be unpacked as a vector

This is row major unpacked so
first 1,...,n entries are susceptibles in urban areas 1,...,n
then n+1,...,2n are exposeds in urban areas 1,...,n
.
.
.
then n*num_states + 1, ..., n*num_states + n entries are susceptibles in rural areas 1,...,n
.
.
.

States:
1 -> S
2 -> E
3 -> I_subclinical
4 -> I_diseased
5 -> H(ospitalised)
6 -> Recovered
7 -> Cumulative I_sub
8 -> Cumulative I_dis
9 -> Cumulative Dead

Events for each location/county:

1-> Urban transmission
2-> Rural transmission
3-> #urban E->A
4-> rural E->A
5-> urban E->D
6-> rural E->D
7-> urban D->H
8-> rural D->H
9-> urban H->R
10-> rural H->R
11-> urban D->R
12-> rural D->R
13-> urban A->R
14-> rural A->R
15-> urban H->death
16-> rural H->death

"""
dc = zeros(Int64,n*n_s*2,n_t*n)

function import_rate_mom(t,into_mom,global_prev)
    if t+1>min(length(into_mom),length(global_prev))
        t_int=min(length(into_mom),length(global_prev))
    else
        t_int=Int(floor(t))+1
    end
    return into_mom[t_int]*global_prev[t_int]
end

function import_rate_nai(t,into_nai,global_prev)
    if t+1>min(length(into_nai),length(global_prev))
        t_int=min(length(into_nai),length(global_prev))
    else
        t_int=Int(floor(t))+1
    end
    return into_nai[t_int]*global_prev[t_int]
end

function calculate_infection_rates!(u,p::CoVParameters,t)
    I_urb_A = @view u[((3-1)*n + 1):((3-1)*n + n)]
    I_urb_D = @view u[((4-1)*n + 1):((4-1)*n + n)]
    I_rur_A = @view u[((3-1)*n + n_s*n + 1):((3-1)*n + n_s*n + n)]
    I_rur_D = @view u[((4-1)*n + n_s*n + 1):((4-1)*n + n_s*n + n)]
    mul!(p.Î,p.T,p.ϵ*I_urb_A .+ I_urb_D  )
    p.Î .+=  p.ϵ*I_rur_A .+ I_rur_D
    mul!(p.λ_urb,p.T',p.β .*(p.Î ./p.N̂))
    p.λ_urb[ind_mombasa] += p.ext_inf_rate*import_rate_mom(t,p.into_mom,p.global_prev)
    p.λ_urb[ind_nairobi] += p.ext_inf_rate*import_rate_nai(t,p.into_nai,p.global_prev)
    p.λ_rur .= p.β .*(p.Î ./p.N̂)
    return nothing
end

function rates(out,u,p::CoVParameters,t)
    @unpack λ_urb,λ_rur,β,γ,σ,δ,τ,μ₁,into_mom,into_nai,global_prev = p
    calculate_infection_rates!(u,p,t)
    for i = 1:n
        out[(i-1)*n_t+1] = λ_urb[i]*u[(1-1)*n + i] #urban transmission
        out[(i-1)*n_t+2] =  λ_rur[i]*u[(1-1)*n + n_s*n + i] #rural transmission
        out[(i-1)*n_t+3] = (1-δ)*σ*u[(2-1)*n + i] #urban E->A
        out[(i-1)*n_t+4] =  (1-δ)*σ*u[(2-1)*n + n_s*n + i] #rural E->A
        out[(i-1)*n_t+5] = δ*σ*u[(2-1)*n + i] #urban E->D
        out[(i-1)*n_t+6] =  δ*σ*u[(2-1)*n + n_s*n + i] #rural E->D
        out[(i-1)*n_t+7] = τ*u[(4-1)*n + i] #urban D->H
        out[(i-1)*n_t+8] =  τ*u[(4-1)*n + n_s*n + i] #rural D->H
        out[(i-1)*n_t+9] = γ*u[(5-1)*n + i] #urban H->R
        out[(i-1)*n_t+10] =  γ*u[(5-1)*n + n_s*n + i] #rural H->R
        out[(i-1)*n_t+11] = γ*u[(4-1)*n + i] #urban D->R
        out[(i-1)*n_t+12] =  γ*u[(4-1)*n + n_s*n + i] #rural D->R
        out[(i-1)*n_t+13] = γ*u[(3-1)*n + i] #urban A->R
        out[(i-1)*n_t+14] =  γ*u[(3-1)*n + n_s*n + i] #rural A->R
        out[(i-1)*n_t+15] = μ₁*u[(5-1)*n + i] #urban H->death
        out[(i-1)*n_t+16] =  μ₁*u[(5-1)*n + n_s*n + i] #rural H->death
    end
end

function change_matrix(dc,u,p,t,mark)
    for i = 1:n
        dc[(1-1)*n + i,(i-1)*n_t+1] = -1
        dc[(2-1)*n + i,(i-1)*n_t+1] = 1 #change due to urban transmission
        dc[(1-1)*n + n_s*n + i,(i-1)*n_t+2] = -1
        dc[(2-1)*n + n_s*n + i,(i-1)*n_t+2] = 1 #change due to rural transmission
        dc[(2-1)*n + i,(i-1)*n_t+3] = -1
        dc[(3-1)*n + i,(i-1)*n_t+3] = 1
        dc[(7-1)*n + i,(i-1)*n_t+3] = 1#change due to urban E->A
        dc[(2-1)*n + n_s*n + i,(i-1)*n_t+4] = -1
        dc[(3-1)*n + n_s*n + i,(i-1)*n_t+4] = 1
        dc[(7-1)*n + n_s*n + i,(i-1)*n_t+4] = 1 #change due to rural E->A
        dc[(2-1)*n + i,(i-1)*n_t+5] = -1
        dc[(4-1)*n + i,(i-1)*n_t+5] = 1
        dc[(8-1)*n + i,(i-1)*n_t+5] = 1 #change due to urban E->D
        dc[(2-1)*n + n_s*n + i,(i-1)*n_t+6] = -1
        dc[(4-1)*n + n_s*n + i,(i-1)*n_t+6] = 1
        dc[(8-1)*n + n_s*n + i,(i-1)*n_t+6] = 1 #change due to rural E->D
        dc[(4-1)*n + i,(i-1)*n_t+7] = -1
        dc[(5-1)*n + i,(i-1)*n_t+7] = 1 #change due to urban D->H
        dc[(4-1)*n + n_s*n + i,(i-1)*n_t+8] = -1
        dc[(5-1)*n + n_s*n + i,(i-1)*n_t+8] = 1 #change due to rural D->H
        dc[(5-1)*n + i,(i-1)*n_t+9] = -1
        dc[(6-1)*n + i,(i-1)*n_t+9] = 1 #change due to urban H->R
        dc[(5-1)*n + n_s*n + i,(i-1)*n_t+10] = -1
        dc[(6-1)*n + n_s*n + i,(i-1)*n_t+10] = 1 #change due to rural H->R
        dc[(4-1)*n + i,(i-1)*n_t+11] = -1
        dc[(6-1)*n + i,(i-1)*n_t+11] = 1 #change due to urban D->R
        dc[(4-1)*n + n_s*n + i,(i-1)*n_t+12] = -1
        dc[(6-1)*n + n_s*n + i,(i-1)*n_t+12] = 1 #change due to rural D->R
        dc[(3-1)*n + i,(i-1)*n_t+13] = -1
        dc[(6-1)*n + i,(i-1)*n_t+13] = 1 #change due to urban A->R
        dc[(3-1)*n + n_s*n + i,(i-1)*n_t+14] = -1
        dc[(6-1)*n + n_s*n + i,(i-1)*n_t+14] = 1 #change due to rural A->R
        dc[(5-1)*n + i,(i-1)*n_t+15] = -1
        dc[(9-1)*n + i,(i-1)*n_t+15] = 1#change due to urban H->death
        dc[(5-1)*n + n_s*n + i,(i-1)*n_t+16] = -1
        dc[(9-1)*n + n_s*n + i,(i-1)*n_t+16] = 1#change due to rural H->death
    end
end

reg_jumps_forKenyaCoV = RegularJump(rates,change_matrix,dc;constant_c=true)

#This generates the underlying Poisson drivers according to the rate function

function PP_drivers(dN::Vector{Int64},rates,p)
    for i = 1:length(dN)
        if rates[i] >= 0.
            dN[i] = rand(Poisson(p.dt*rates[i]))
        else
            dN[i] = 0
        end
    end
end

#This caps the Poisson processes at causing no more transitions than the group being effected
function max_change(out,u,p::CoVParameters)
    for i = 1:n
        out[(i-1)*n_t+1] = min(out[(i-1)*n_t+1],u[(1-1)*n + i] )#urban transmission
        out[(i-1)*n_t+2] =  min(out[(i-1)*n_t+2],u[(1-1)*n + n_s*n + i] )#rural transmission
        # out[(i-1)*n_t+3] = min(out[(i-1)*n_t+3],u[(2-1)*n + i] )#urban E->A
        # out[(i-1)*n_t+4] =  min(out[(i-1)*n_t+4],u[(2-1)*n + n_s*n + i] )#rural E->A
        # out[(i-1)*n_t+5] = min(out[(i-1)*n_t+5],u[(2-1)*n + i] )#urban E->D
        # out[(i-1)*n_t+6] =  min(out[(i-1)*n_t+6],u[(2-1)*n + n_s*n + i] )#rural E->D
        # out[(i-1)*n_t+7] = min(out[(i-1)*n_t+7],u[(4-1)*n + i] )#urban D->H
        # out[(i-1)*n_t+8] =  min(out[(i-1)*n_t+8],u[(4-1)*n + n_s*n + i] )#rural D->H
        out[(i-1)*n_t+9] = min(out[(i-1)*n_t+9],u[(5-1)*n + i] )#urban H->R
        out[(i-1)*n_t+10] = min(out[(i-1)*n_t+10],u[(5-1)*n + n_s*n + i] )#rural H->R
        # out[(i-1)*n_t+11] = min(out[(i-1)*n_t+11],u[(4-1)*n + i] )#urban D->R
        # out[(i-1)*n_t+12] =  min(out[(i-1)*n_t+12],u[(4-1)*n + n_s*n + i] )#rural D->R
        out[(i-1)*n_t+13] = min(out[(i-1)*n_t+13],u[(3-1)*n + i] )#urban A->R
        out[(i-1)*n_t+14] =  min(out[(i-1)*n_t+14],u[(3-1)*n + n_s*n + i] )#rural A->R
        out[(i-1)*n_t+15] = min(out[(i-1)*n_t+15],u[(5-1)*n + i] )#urban H->death
        out[(i-1)*n_t+16] =  min(out[(i-1)*n_t+16],u[(5-1)*n + n_s*n + i] )#rural H->death
        if out[(i-1)*n_t+3] + out[(i-1)*n_t+5] > u[(2-1)*n + i] #More incubations than actual urban E population
            out[(i-1)*n_t+3] = rand(Binomial(u[(2-1)*n + i],p.δ)) #Binomially distributed the urban incubations between Asympotomatic and symptomatic
            out[(i-1)*n_t+5] = u[(2-1)*n + i] - out[(i-1)*n_t+3]
        end
        if out[(i-1)*n_t+4] + out[(i-1)*n_t+6] > u[(2-1)*n + n_s*n + i]  #More incubations than actual rural E population
            out[(i-1)*n_t+4] = rand(Binomial(u[(2-1)*n + n_s*n + i],p.δ)) #Binomially distributed the rural incubations between Asympotomatic and symptomatic
            out[(i-1)*n_t+6] = u[(2-1)*n + n_s*n + i] - out[(i-1)*n_t+4]
        end
        if out[(i-1)*n_t+7] + out[(i-1)*n_t+11] > u[(4-1)*n + i] # More end of infection events than urban symptomatic infecteds
            out[(i-1)*n_t+7] = rand(Binomial(u[(4-1)*n + i],p.τ/(p.τ + p.γ))) #Binomially distributed the urban end of infections between hospitalisation and recovery
            out[(i-1)*n_t+11] = u[(4-1)*n + i] - out[(i-1)*n_t+7]
        end
        if out[(i-1)*n_t+8] + out[(i-1)*n_t+12] > u[(4-1)*n + n_s*n + i] # More end of infection events than rural symptomatic infecteds
            out[(i-1)*n_t+8] = rand(Binomial(u[(4-1)*n + n_s*n + i],p.τ/(p.τ + p.γ))) #Binomially distributed the rural end of infections between hospitalisation and recovery
            out[(i-1)*n_t+12] = u[(4-1)*n + n_s*n + i] - out[(i-1)*n_t+7]
        end
    end
end

function nonneg_tauleap(du,u,p::CoVParameters,t)
    @unpack dc,dN,poi_rates = p
    rates(poi_rates,u,p,t) #calculate rates of underlying Poisson processes
    PP_drivers(dN,poi_rates,p)#Generate Poisson rvs with rates scaled by time step dt
    max_change(dN,u,p)#Cap the size of the Poisson rvs to maintain non-negativity
    mul!(du,dc,dN)#Calculates the effect on the state in the inplace du vector
    du .+= u #Calculates how the state should change
end

function ode_model(du,u,p::CoVParameters,t)
    @unpack λ_urb,λ_rur,β,γ,σ,δ,τ,μ₁,into_mom,into_nai,global_prev = p
    calculate_infection_rates!(u,p,t)
    for i = 1:n
        du[i,1,1] = (-1)*λ_urb[i]*u[i,1,1]
        du[i,2,1] = λ_urb[i]*u[i,1,1] - σ*u[i,2,1]
        du[i,3,1] = (1-δ)*σ*u[i,2,1] - γ*u[i,3,1]
        du[i,4,1] = δ*σ*u[i,2,1] - γ*u[i,4,1] - τ*u[i,4,1]
        du[i,5,1] = τ*u[i,4,1] - γ*u[i,5,1] - μ₁*u[i,5,1]
        du[i,6,1] = γ*(u[i,3,1] + u[i,4,1] + u[i,5,1])
        du[i,7,1] = (1-δ)*σ*u[i,2,1]
        du[i,8,1] = δ*σ*u[i,2,1]
        du[i,9,1] = μ₁*u[i,5,1]
        du[i,1,2] = (-1)*λ_rur[i]*u[i,1,2]
        du[i,2,2] = λ_rur[i]*u[i,1,2] - σ*u[i,2,2]
        du[i,3,2] = (1-δ)*σ*u[i,2,2] - γ*u[i,3,2]
        du[i,4,2] = δ*σ*u[i,2,2] - γ*u[i,4,2] - τ*u[i,4,2]
        du[i,5,2] = τ*u[i,4,2] - γ*u[i,5,2] - μ₁*u[i,5,2]
        du[i,6,2] = γ*(u[i,3,2] + u[i,4,2] + u[i,5,2])
        du[i,7,2] = (1-δ)*σ*u[i,2,2]
        du[i,8,2] = δ*σ*u[i,2,2]
        du[i,9,2] = μ₁*u[i,5,2]
     end
end


text:
"""
	colorbar(cmd0::String="", arg1=[]; kwargs...)

Plots gray scales or color scales on maps.

Full option list at [`psscale`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html)

- **D** : **position** : -- Str --

    Defines the reference point on the map for the color scale using one of four coordinate systems.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#d)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **F** : **box** : -- Str --

    Draws a rectangular border around the scale.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#f)
- **G** : **truncate** : -- Str --

    Truncate the incoming CPT so that the lowest and highest z-levels are to zlo and zhi.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#g)
- **I** : **shade** : -- Number or [] --

    Add illumination effects.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#i)
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- **L** : **equal_size** : -- Str or [] --

    Gives equal-sized color rectangles. Default scales rectangles according to the z-range in the CPT.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#l)
- **M** : **monochrome** : -- Bool or [] --

    Force conversion to monochrome image using the (television) YIQ transformation.
    [`-M`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#m)
- **N** : **dpi** : -- Str or number --

    Controls how the color scale is represented by the PostScript language.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#n)
- **Q** : **log** : -- Str --

    Selects a logarithmic interpolation scheme [Default is linear].
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#q)
- $(GMT.opt_R)
- **S** : **nolines** : -- Bool or [] --

    Do not separate different color intervals with black grid lines.
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** : **zscale** : -- Number --

    Multiply all z-values in the CPT by the provided scale.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#w)
- **Z** : **zfile** : -- Str --

    File with colorbar-width per color entry.
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/psscale.html#z)
"""
function colorbar(cmd0::String="", arg1=[]; K=false, O=false, first=true, kwargs...)

    length(kwargs) == 0 && isempty(data) && return monolitic("psscale", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

    cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, "")
	cmd = parse_UVXY(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd  = parse_params(cmd, d)

	cmd, arg1, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', 0, arg1, [])

	cmd = add_opt(cmd, 'D', d, [:D :position])
	cmd = add_opt(cmd, 'F', d, [:F :box])
	cmd = add_opt(cmd, 'G', d, [:G :truncate])
    cmd = add_opt(cmd, 'I', d, [:I :shade])
	cmd = add_opt(cmd, 'L', d, [:L :equal_size])
	cmd = add_opt(cmd, 'M', d, [:M :monochrome])
	cmd = add_opt(cmd, 'N', d, [:N :dpi])
    cmd = add_opt(cmd, 'Q', d, [:Q :log])
	cmd = add_opt(cmd, 'S', d, [:S :nolines])
	cmd = add_opt(cmd, 'W', d, [:W :zscale])
	cmd = add_opt(cmd, 'Z', d, [:Z :zfile])

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	cmd = finish_PS(d, cmd, output, K, O)
    return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "psscale", arg1)
end

# ---------------------------------------------------------------------------------------------------
colorbar!(cmd0::String="", arg1=[]; K=false, O=false, first=false, kw...) =
    colorbar(cmd0, arg1; K=K, O=O, first=first, kw...)

psscale  = colorbar         # Alias
psscale! = colorbar!        # Alias


text:
using Documenter, MakieLayout, CairoMakie, AbstractPlotting, Animations, ColorSchemes, Colors

CairoMakie.activate!()

makedocs(
    sitename="MakieLayout.jl",
    pages = [
        "index.md",
        "GridLayout" => "grids.md",
        "LAxis" => "laxis.md",
        "LLegend" => "llegend.md",
        "Layoutables Examples" => "layoutables_examples.md",
        "Theming Layoutables" => "theming.md",
        "How Layouting Works" => "layouting.md",
        "Frequently Asked Questions" => "faq.md",
        "API Reference" => "api_reference.md",
    ],
    format = Documenter.HTML(
            prettyurls = get(ENV, "CI", nothing) == "true"
        )
    )

deploydocs(
    repo = "github.com/jkrumbiegel/MakieLayout.jl.git",
)

max_stars_repo_name: jkrumbiegel/MakieLayout


text:
# Uncertain population consisting of CertainValues (scalars get promoted to CertainValue)s
# theoretical distributions and KDE distributions
p1 = ConstrainedUncertainScalarPopulation(
        [3.0, UncertainValue(Normal, 0, 1),
        UncertainValue(Gamma, 2, 3),
        UncertainValue(Uniform, rand(1000))],
        [0.5, 0.5, 0.5, 0.5]
    )

# Uncertain population consisting of scalar values
p2 = ConstrainedUncertainScalarPopulation([1, 2, 3], rand(3))
p3 = ConstrainedUncertainScalarPopulation([1, 2, 3], Weights(rand(3)))

# Uncertain population consisting of uncertain populations
p4 = ConstrainedUncertainScalarPopulation([p1, p2], [0.1, 0.5])
p5 = ConstrainedUncertainScalarPopulation([p1, p2], Weights([0.1, 0.5]));

@test p1 isa ConstrainedUncertainScalarPopulation{T1, T2} where {T1 <: AbstractUncertainValue, T2 <: AbstractWeights}
@test p2 isa ConstrainedUncertainScalarPopulation{T1, T2} where {T1 <: Number, T2 <: AbstractWeights}
@test p3 isa ConstrainedUncertainScalarPopulation{T1, T2} where {T1 <: Number, T2 <: AbstractWeights}
@test p4 isa ConstrainedUncertainScalarPopulation{T1, T2} where {T1 <: AbstractUncertainValue, T2 <: AbstractWeights}
@test p5 isa ConstrainedUncertainScalarPopulation{T1, T2} where {T1 <: AbstractUncertainValue, T2 <: AbstractWeights}


text:
function hailstonelength(n::Integer)
    len = 1
    while n > 1
        n = ifelse(iseven(n), n ÷ 2, 3n + 1)
        len += 1
    end
    return len
end

@show hailstonelength(27); nothing
@show findmax([hailstonelength(i) for i in 1:100_000]); nothing

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    AutocorLenAlgorithm

Abstract type for integrated autocorrelation length estimation algorithms.
"""
abstract type AutocorLenAlgorithm end
export AutocorLenAlgorithm

"""
    bat_integrated_autocorr_len(
        v::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}},
        algorithm::AutocorLenAlgorithm = GeyerAutocorLen()
    )

Estimate the integrated autocorrelation length of variate series `v`,
separately for each degree of freedom.

Returns a NamedTuple of the shape

```julia
(result = integrated_autocorr_len, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable BAT API.

!!! note

    Do not add add methods to `bat_integrated_autocorr_len`, add methods to
    `bat_integrated_autocorr_len_impl` instead (same arguments).
"""
function bat_integrated_autocorr_len end
export bat_integrated_autocorr_len

function bat_integrated_autocorr_len_impl end

function bat_integrated_autocorr_len(
    v::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}},
    algorithm = bat_default_withdebug(bat_integrated_autocorr_len, Val(:algorithm), v)
)
    r = bat_integrated_autocorr_len_impl(v, algorithm)
    result_with_args(r, (algorithm = algorithm,))
end

function argchoice_msg(::typeof(bat_integrated_autocorr_len), ::Val{:algorithm}, x::AutocorLenAlgorithm)
    "Using integrated autocorrelation length estimator $x"
end

@doc doc"""
    bat_eff_sample_size(
        v::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}},
        [algorithm::AutocorLenAlgorithm]
    )

    bat_eff_sample_size(
        smpls::DensitySampleVector,
        [algorithm::AutocorLenAlgorithm]
    )

Estimate effective sample size estimation for variate series `v`, resp.
density samples `smpls`, separately for each degree of freedom.

Returns a NamedTuple of the shape

```julia
(result = X::AbstractVector{<:Real}, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable BAT API.

!!! note

    Do not add add methods to `bat_eff_sample_size`, add methods to
    `bat_eff_sample_size_impl` instead (same arguments).
"""
function bat_eff_sample_size end
export bat_eff_sample_size

function bat_eff_sample_size_impl end

function bat_eff_sample_size(
    v::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}},
    algorithm = bat_default_withdebug(bat_eff_sample_size, Val(:algorithm), v)
)
    r = bat_eff_sample_size_impl(v, algorithm)
    result_with_args(r, (algorithm = algorithm,))
end

function bat_eff_sample_size(
    smpls::DensitySampleVector,
    algorithm = bat_default_withdebug(bat_eff_sample_size, Val(:algorithm), smpls);
    kwargs...
)
    r = bat_eff_sample_size_impl(smpls, algorithm; kwargs...)
    result_with_args(r, (algorithm = algorithm,))
end

function argchoice_msg(::typeof(bat_eff_sample_size), ::Val{:algorithm}, x::AutocorLenAlgorithm)
    "Using integrated autocorrelation length estimator $x"
end


text:
"""
    $(TYPEDEF)

Make a two-dimensional bipedal robot walk forward as fast as possible. The robot model is
based on the work T. Erez, Y. Tassa, and E. Todorov: ["Infinite Horizon Model Predictive
Control for Nonlinear Periodic Tasks"](https://homes.cs.washington.edu/~todorov/papers/ErezRSS11.pdf), 2011.

* **State: (29, )**
* **Action: (6, )**
* **Observation: (7, )**
"""
mutable struct Walker2DV2{Sim, SSpace, OSpace} <: AbstractMuJoCoEnvironment
    sim::Sim
    statespace::SSpace
    observationspace::OSpace
    last_torso_x::Float64
    randreset_distribution::Uniform{Float64}
    function Walker2DV2(sim::MJSim)
        sspace = MultiShape(
            simstate=statespace(sim),
            last_torso_x=ScalarShape(Float64)
        )
        ospace = MultiShape(
            cropped_qpos = VectorShape(Float64, sim.m.nq - 1),
            qvel = VectorShape(Float64, sim.m.nv)
        )
        env = new{typeof(sim), typeof(sspace), typeof(ospace)}(
            sim, sspace, ospace, 0, Uniform(-0.005, 0.005)
        )
        reset!(env)
    end
end

Walker2DV2() = first(tconstruct(Walker2DV2, 1))

function LyceumBase.tconstruct(::Type{Walker2DV2}, n::Integer)
    modelpath = joinpath(@__DIR__, "walker2d.xml")
    Tuple(Walker2DV2(s) for s in tconstruct(MJSim, n, modelpath, skip=4))
end

@inline getsim(env::Walker2DV2) = env.sim

@inline LyceumBase.statespace(env::Walker2DV2) = env.statespace

function LyceumBase.getstate!(state, env::Walker2DV2)
    checkaxes(statespace(env), state)
    @uviews state begin
        shaped = statespace(env)(state)
        getstate!(shaped.simstate, env.sim)
        shaped.last_torso_x = env.last_torso_x
    end
    state
end

function LyceumBase.setstate!(env::Walker2DV2, state)
    checkaxes(statespace(env), state)
    @uviews state begin
        shaped = statespace(env)(state)
        setstate!(env.sim, shaped.simstate)
        env.last_torso_x = shaped.last_torso_x
    end
    env
end

@inline LyceumBase.observationspace(env::Walker2DV2) = env.observationspace

function LyceumBase.getobservation!(obs, env::Walker2DV2)
    checkaxes(observationspace(env), obs)
    qpos = env.sim.d.qpos
    @views @uviews qpos obs begin
        shaped = observationspace(env)(obs)
        copyto!(shaped.cropped_qpos, qpos[2:end])
        copyto!(shaped.qvel, env.sim.d.qvel)
        clamp!(shaped.qvel, -10, 10)
    end
    obs
end

function LyceumBase.getreward(state, action, ::Any, env::Walker2DV2)
    checkaxes(statespace(env), state)
    checkaxes(actionspace(env), action)
    @uviews state begin
        shapedstate = statespace(env)(state)
        alive_bonus = 1.0
        reward = (_torso_x(shapedstate, env) - shapedstate.last_torso_x) / timestep(env)
        reward += alive_bonus
        reward -= 1e-3 * sum(x->x^2, action)
        reward
    end
end

function LyceumBase.reset!(env::Walker2DV2)
    reset!(env.sim)
    env.last_torso_x = _torso_x(env)
    env
end

function LyceumBase.randreset!(rng::AbstractRNG, env::Walker2DV2)
    reset_nofwd!(env.sim)
    perturb!(rng, env.randreset_distribution, env.sim.d.qpos)
    perturb!(rng, env.randreset_distribution, env.sim.d.qvel)
    forward!(env.sim)
    env.last_torso_x = _torso_x(env)
    env
end

function LyceumBase.step!(env::Walker2DV2)
    env.last_torso_x = _torso_x(env)
    step!(env.sim)
    env
end

function LyceumBase.isdone(state, o, env::Walker2DV2)
    checkaxes(statespace(env), state)
    @uviews state begin
        shapedstate = statespace(env)(state)
        torso_x = _torso_x(shapedstate, env)
        height = _torso_height(shapedstate, env)
        torso_ang = _torso_ang(shapedstate, env)
        qpos = shapedstate.simstate.qpos
        qvel = shapedstate.simstate.qvel

        done = !(
            height > 0.8
            && height < 2.0
            && abs(torso_ang) < 1.0
        )
        done
    end
end

@inline _torso_x(shapedstate::ShapedView, ::Walker2DV2) = shapedstate.simstate.qpos[1]
@inline _torso_x(env::Walker2DV2) = env.sim.d.qpos[1]
@inline _torso_height(shapedstate::ShapedView, ::Walker2DV2) = shapedstate.simstate.qpos[2]
@inline _torso_ang(shapedstate::ShapedView, ::Walker2DV2) = shapedstate.simstate.qpos[3]


text:
@test isequal(Float64(Float128(42.0)), 42.0)
@test isequal(Float64(Float128(42.0f0)), 42.0)
@test isequal(Float64(Float128(42)), 42.0)
@test isequal(Float64(Float128(0x2a)), 42.0)

@test Float64(eps(Float128)) <= 1.0e-30

@test isequal(big(zero(Float256)), big(0))
@test isequal(big(one(Float256)), big(1))
@test abs(big(Float128(π)) - big(π)) <= eps(Float128)
@test abs(big(Float128(ℯ)) - big(ℯ)) <= eps(Float128)

pi128 = f128"3.141592653589793238462643383279502884197169399375105820974944592307816406286198"
@test isequal(big(pi128), Float128(π))
@test isequal(parse(Float128, string(pi128)), pi128)

for n in 1:1000
    x64 = rand(Float64)
    y64 = rand(Float64)

    @test isequal(Float64(Float128(x64)), x64)

    xbig = rand(BigFloat)
    @test isequal(big(Float128(x64)), big(x64))
    @test abs(big(Float128(xbig)) - xbig) <= eps(Float128)

    x = Float128(rand(BigFloat))
    y = Float128(rand(BigFloat))

    @test isequal(big(+ x), + big(x))
    @test isequal(big(- x), - big(x))
    @test isequal(big(abs(x)), abs(big(x)))
    @test approxeq(Float128, 2, [x], big(inv(x)),  inv(big(x)))

    @test abs(big(x + y) - (big(x) + big(y))) <= eps(Float128)
    @test abs(big(x64 + y) - (big(x64) + big(y))) <= eps(Float128)
    @test abs(big(x + y64) - (big(x) + big(y64))) <= eps(Float128)

    @test abs(big(x - y) - (big(x) - big(y))) <= eps(Float128)
    @test abs(big(x64 - y) - (big(x64) - big(y))) <= eps(Float128)
    @test abs(big(x - y64) - (big(x) - big(y64))) <= eps(Float128)

    @test abs(big(x * y) - (big(x) * big(y))) <= eps(Float128)
    @test abs(big(x64 * y) - (big(x64) * big(y))) <= eps(Float128)
    @test abs(big(x * y64) - (big(x) * big(y64))) <= eps(Float128)

    @test approxeq(Float128, 2, [x, y], big(x / y), big(x) / big(y))
    @test approxeq(Float128, 2, [x, y], big(x64 / y), big(x64) / big(y))
    @test approxeq(Float128, 2, [x, y], big(x / y64), big(x) / big(y64))

    @test approxeq(Float128, 2, [x, y], big(x \ y), big(x) \ big(y))
    @test approxeq(Float128, 2, [x, y], big(x64 \ y), big(x64) \ big(y))
    @test approxeq(Float128, 2, [x, y], big(x \ y64), big(x) \ big(y64))

    @test isequal(x == x, big(x) == big(x))
    @test isequal(x != x, big(x) != big(x))
    @test isequal(x < x, big(x) < big(x))
    @test isequal(x <= x, big(x) <= big(x))
    @test isequal(x > x, big(x) > big(x))
    @test isequal(x >= x, big(x) >= big(x))

    @test isequal(x == y, big(x) == big(y))
    @test isequal(x != y, big(x) != big(y))
    @test isequal(x < y, big(x) < big(y))
    @test isequal(x <= y, big(x) <= big(y))
    @test isequal(x > y, big(x) > big(y))
    @test isequal(x >= y, big(x) >= big(y))

    @test isequal(x64 == y, big(x64) == big(y))
    @test isequal(x64 != y, big(x64) != big(y))
    @test isequal(x64 < y, big(x64) < big(y))
    @test isequal(x64 <= y, big(x64) <= big(y))
    @test isequal(x64 > y, big(x64) > big(y))
    @test isequal(x64 >= y, big(x64) >= big(y))

    @test isequal(x == y64, big(x) == big(y64))
    @test isequal(x != y64, big(x) != big(y64))
    @test isequal(x < y64, big(x) < big(y64))
    @test isequal(x <= y64, big(x) <= big(y64))
    @test isequal(x > y64, big(x) > big(y64))
    @test isequal(x >= y64, big(x) >= big(y64))

    @test isequal(big(floor(1.0e10 * x)), floor(big(1.0e10 * x)))
    @test isequal(big(ceil(1.0e10 * x)), ceil(big(1.0e10 * x)))
    @test isequal(big(round(1.0e10 * x)), round(big(1.0e10 * x)))
    @test isequal(big(trunc(1.0e10 * x)), trunc(big(1.0e10 * x)))

    @test abs(big(x ^ 2) - big(x) ^ 2) <= eps(Float128)
    @test approxeq(Float128, 10, [x], big(x ^ 5), big(x) ^ 5)
    @test approxeq(Float128, 2, [], big(sqrt(abs(x))), sqrt(big(abs(x))))
    @test approxeq(Float128, 10, [abs(x)],
                   big(nroot(abs(x), 5)), big(abs(x)) ^ (big(1)/5))

    @test approxeq(Float128, 10, [x], big(exp(x)), exp(big(x)))
    @test approxeq(Float128, 10, [x], big(exp10(x)), exp10(big(x)))
    @test approxeq(Float128, 10, [x], big(exp2(x)), exp2(big(x)))
    @test approxeq(Float128, 10,
                   [abs(x)], big(log(abs(x))), log(big(abs(x))))
    @test approxeq(Float128, 10,
                   [abs(x)], big(log10(abs(x))), log10(big(abs(x))))
    @test approxeq(Float128, 10,
                   [abs(x)], big(log2(abs(x))), log2(big(abs(x))))

    @test approxeq(Float128, 10, [x], big(sin(x)), sin(big(x)))
    @test approxeq(Float128, 10, [x], big(cos(x)), cos(big(x)))
    @test approxeq(Float128, 10, [x], big(tan(x)), tan(big(x)))
    @test approxeq(Float128, 10, [x], big(sincos(x)[1]), sin(big(x)))
    @test approxeq(Float128, 10, [x], big(sincos(x)[2]), cos(big(x)))

    @test approxeq(Float128, 10, [x], big(asin(x)), asin(big(x)))
    @test approxeq(Float128, 100, [x], big(acos(x)), acos(big(x)))
    @test approxeq(Float128, 10, [x], big(atan(x)), atan(big(x)))
    @test approxeq(Float128, 10, [x, y],
                   big(atan(x, y)), atan(big(x), big(y)))
    @test approxeq(Float128, 10, [x, y],
                   big(hypot(x, y)), hypot(big(x), big(y)))

    @test approxeq(Float128, 10, [x], big(sinh(x)), sinh(big(x)))
    @test approxeq(Float128, 10, [x], big(cosh(x)), cosh(big(x)))
    @test approxeq(Float128, 10, [x], big(tanh(x)), tanh(big(x)))
    @test approxeq(Float128, 10, [x], big(sincosh(x)[1]), sinh(big(x)))
    @test approxeq(Float128, 10, [x], big(sincosh(x)[2]), cosh(big(x)))

    @test approxeq(Float128, 10, [1, x], big(asinh(x)), asinh(big(x)))
    @test approxeq(Float128, 10, [x + 2], big(acosh(x + 2)), acosh(big(x + 2)))
    @test approxeq(Float128, 10, [1, x], big(atanh(x)), atanh(big(x)))
end

n = 100

A64 = rand(Float64, n, n)
B64 = inv(A64)
@test norm(B64 * A64 - I) <= n * norm(A64) * norm(B64) * eps(Float64)

A128 = rand(Float128, n, n)
B128 = inv(A128)
@test norm(B128 * A128 - I) <= n * norm(A128) * norm(B128) * eps(Float128)


text:
# Modelo de trocador de calor, dada a carga termica
type Trocador
	Trocador()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Potencia(),
			Dif_Pres(),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}"
			)),
			[
				:(Fout.F = Fin.F),
				:(Fout.P = Fin.P - DP),
				:(Fout.F * (Fout.H - Fin.H) = Q),
				:([Fout.S,Fout.T] = PP2.propPH(Fout.P,Fout.H)),
			],
			[
				"","","","",
			],
			[:PP2,],
			[:Q,:DP,:Fin,:Fout,]
		)
	end
	PP2::DanaPlugin
	Q::Potencia
	DP::Dif_Pres
	Fin::Corrente
	Fout::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
end
export Trocador
function setEquationFlow(in::Trocador)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
end


text:
# This file was generated, do not modify it. # hide
feature2 = coerce(feature, Multiclass)
elscitype(feature2)


text:
# ---------------------------------------------------------------------------- #
#
#   meshF.jl
#
#   Abstract meshF type
#   This allows for writing an n-dimensional methods for lattice optimization
#
#   intrico
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

using GroupSlices

"""
    MeshF

MeshF abstract type:
Overarching abstract type for mesh types (2D and 3D).
"""
abstract type MeshF

end

include("meshF2D.jl")
include("meshF3D.jl")


text:
function iterate_interpolations(source::String)
    result = Union{Expr, JSString, Symbol}[]
    lastidx = 1; i = 1; lindex = lastindex(source)
    while true
        c = source[i]
        if c == '$'
            # add elements before $
            if !isempty(lastidx:(i - 1))
                push!(result, JSString(source[lastidx:(i - 1)]))
            end
            # parse the $ expression
            expr, i2 = Meta.parse(source, i + 1, greedy = false, raise = false)
            if i2 >= lindex && expr === nothing
                error("Invalid interpolation at index $(i)-$(lindex): $(source[i:lindex])")
            end
            i = i2
            push!(result, esc(expr))
            lastidx = i
            i > lindex && break
        else
            if i == lindex
                if !isempty(lastidx:lindex)
                    push!(result, JSString(source[lastidx:lindex]))
                end
                break
            end
            i = Base.nextind(source, i)
        end
    end
    return result
end

macro js_str(js_source)
    value_array = :([])
    append!(value_array.args, iterate_interpolations(js_source))
    return :(JSCode($value_array))
end

function Base.show(io::IO, jsc::JSCode)
    print_js_code(io, jsc, SerializationContext(nothing))
end

function print_js_code(io::IO, @nospecialize(object), context)
    serialized = serialize_js(context, object)
    if serialized isa Union{Number, String}
        return print_js_code(io, serialized, context)
    end
    if isnothing(context.interpolated)
        json = JSON3.write(serialized)
        print(io, "JSServe.deserialize_js($(json))")
    else
        index = length(context.interpolated) # 0 indexed
        push!(context.interpolated, serialized)
        print(io, "__eval_context__[$(index)]")
    end
end

function print_js_code(io::IO, x::Number, context)
    print(io, x)
    return context
end

function print_js_code(io::IO, x::String, context)
    print(io, "'", x, "'")
    return context
end

function print_js_code(io::IO, jss::JSString, context)
    print(io, jss.source)
    return context
end

function print_js_code(io::IO, dep::Dependency, context)
    print(io, dep.name)
    return context
end

function print_js_code(io::IO, jsc::JSCode, context)
    for elem in jsc.source
        print_js_code(io, elem, context)
    end
    return context
end

function print_js_code(io::IO, jsss::AbstractVector{JSCode}, context)
    for jss in jsss
        print_js_code(io, jss, context)
        println(io)
    end
    return context
end

function jsrender(session::Session, js::JSCode)
    register_resource!(session, js)
    source = sprint() do io
        println(io)
        println(io, js)
    end
    return DOM.script(source, type="text/javascript")
end


text:
using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
@info "Loading DiffEqSensitivity"
using DiffEqSensitivity
@info "Loading Optim"
using Optim
@info "Loading DiffEqFlux"
using DiffEqFlux
@info "Loading Flux"
using Flux
@info "Loading OrdinaryDiffEq"
using OrdinaryDiffEq
using Test
using JLD2
@info "Finished loading packages"

const losses = []

@info "Started Lotka Volterra UODE Testset"
@info "Generate data"

function lotka_volterra(du, u, p, t)
    α, β, γ, δ = p
    du[1] = α*u[1] - β*u[2]*u[1]
    du[2] = γ*u[1]*u[2]  - δ*u[2]
end

tspan = (0.0f0,3.0f0)
u0 = Float32[0.44249296,4.6280594]
p_ = Float32[1.3, 0.9, 0.8, 1.8]
prob = ODEProblem(lotka_volterra, u0,tspan, p_)
solution = solve(prob, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 0.1)

# Ideal data
tsdata = Array(solution)
# Add noise to the data
noisy_data = tsdata + Float32(1e-5)*randn(eltype(tsdata), size(tsdata))

@info "Setup neural network and auxillary functions"
ann = FastChain(FastDense(2, 32, tanh),FastDense(32, 64, tanh),FastDense(64, 32, tanh), FastDense(32, 2))
p = initial_params(ann)

function dudt_(u, p,t)
    x, y = u
    z = ann(u,p)
    [p_[1]*x + z[1],
    -p_[4]*y + z[2]]
end

prob_nn = ODEProblem(dudt_,u0, tspan, p)
s = solve(prob_nn, Tsit5(), u0 = u0, p = p, saveat = solution.t)

function predict(θ)
    Array(solve(prob_nn, Vern7(), u0 = u0, p = θ, saveat = solution.t,
                         abstol=1e-6, reltol=1e-6))
end

# No regularisation right now
function loss(θ)
    pred = predict(θ)
    sum(abs2, noisy_data .- pred), pred # + 1e-5*sum(sum.(abs, params(ann)))
end

callback(θ,l,pred) = begin
    push!(losses, l)
    if length(losses)%10 == 0
        println("Loss after $(length(losses)) iterations $(losses[end])")
    end
    false
end

if !isfile(joinpath(dirname(@__FILE__), "partial_lotka_volterra.jld2"))
    @info "Initial loss $(loss(p)[1])"
    @info "Train neural network until converged"
    res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters = 100)
    @info "Finished initial training with loss $(losses[end])"
    res2 = DiffEqFlux.sciml_train(loss, res1.minimizer, BFGS(initial_stepnorm=0.01), cb=callback, maxiters = 10000)
    @info "Finished extended training with loss $(losses[end])"
    p_trained = res2.minimizer
    @save "$(joinpath(dirname(@__FILE__), "partial_lotka_volterra.jld2"))" p_trained

    # Check for the convergence of loss
    # Necessary for result
    @test losses[end] <= 1e-2
else
    @info "Loading pretrained parameters"
    @load "$(joinpath(dirname(@__FILE__), "partial_lotka_volterra.jld2"))" p_trained
end

# Plot the data and the approximation
NNsolution = predict(p_trained)

@test norm(NNsolution .- solution, 2) < 0.1

Y = ann(noisy_data, p_trained)
X = noisy_data
# Create a Basis
@variables u[1:2]
# Lots of polynomials
polys = Any[1]
for i ∈ 1:5
    push!(polys, u[1]^i)
    push!(polys, u[2]^i)
    for j ∈ i:5
        if i != j
            push!(polys, (u[1]^i)*(u[2]^j))
            push!(polys, u[2]^i*u[1]^i)
        end
    end
end

# And some other stuff
h = [cos.(u)...; sin.(u)...; polys...]
basis = Basis(h, u)

# Create an optimizer for the SINDy problem
opt = SR3()
# Create the thresholds which should be used in the search process
λ = exp10.(-10:0.05:-0.5)
# Target function to choose the results from; x = L0 of coefficients and L2-Error of the model
g(x) = x[1] < 1 ? Inf : norm(x, 2)
@info "Start SINDy regression with unknown threshold"
# Test on uode derivative data
Ψ = SINDy(X[:, 2:end], Y[:, 2:end], basis, λ, opt, g = g, maxiter = 10000, normalize = true, denoise = true) # Succeed
print_equations(Ψ)
p̂ = parameters(Ψ)
@info "Build initial guess system"
# The parameters are a bit off, so we reiterate another SINDy term to get closer to the ground truth
# Create function
unknown_sys = ODESystem(Ψ)
unknown_eq = ODEFunction(unknown_sys)
# Just the equations
b = Basis((u, p, t)->unknown_eq(u, ones(size(p̂)), t), u)
# Test on uode derivative data
@info "Refine the guess"
Ψ = SINDy(X[:, 2:end], Y[:, 2:end],b, SR3(0.1), maxiter = 1000) # Succeed
p̂ = parameters(Ψ)

@info "Checking equations"
found_basis = map(x->simplify(x.rhs), equations(Ψ.equations))
pps = parameters(Ψ.equations)

expected_eqs = Vector{Any}()
for _p in pps
    push!(expected_eqs, _p*u[1]*u[2])
end
@test all(isequal.(found_basis, expected_eqs))
@test isapprox(abs.(p̂), p_[2:3], atol = 9e-2)

# The parameters are a bit off, so we reiterate another SINDy term to get closer to the ground truth
# Create function
unknown_sys = ODESystem(Ψ)
unknown_eq = ODEFunction(unknown_sys)

# Build a ODE for the estimated system
function approx(du, u, p, t)
    # Add SINDy Term
    α, δ, β, γ = p
    z = unknown_eq(u, [β; γ], t)
    du[1] = α*u[1] + z[1]
    du[2] = -δ*u[2] + z[2]
end
@info "Simulate system"
# Create the approximated problem and solution
ps = [p_[[1,4]]; p̂]
a_prob = ODEProblem(approx, u0, tspan, ps)
a_solution = solve(a_prob, Vern7(), abstol=1e-12, reltol=1e-12, saveat = solution.t)
@test norm(a_solution-solution, 2) < 5e-1
@test norm(a_solution-solution, Inf) < 1.5e-1

@info "Finished"


text:
module Configurations

export @option, from_dict, from_kwargs, from_toml, to_dict, to_toml

using TOML
using ExproniconLite
using OrderedCollections

using Crayons.Box

export no_default,
    Maybe,
    # reflection
    field_aliases,
    field_default,
    field_alias,
    type_alias,
    # traits
    is_option,
    # parse
    from_dict,
    from_kwargs,
    from_toml,
    from_toml_if_exists,
    # serialize
    to_dict

include("types.jl")
include("convert.jl")
include("reflection.jl")
include("utils.jl")
include("parse.jl")
include("codegen.jl")
include("printing.jl")

include("serialize.jl")

end


text:
using LinearAlgebra

# Define a generic linear operator with preallocation
function preallocated_LinearOperator(A)
  (n, m) = size(A)
  Ap = zeros(n)
  Atq = zeros(m)
  return LinearOperator(n, m, false, false, p -> mul!(Ap, A, p),
                        q -> mul!(Atq, transpose(A), q), q -> mul!(Atq, transpose(A), q))
end

# Variants where A is a matrix without specific properties
for fn in (:cgls, :cgne, :craig, :craigmr, :crls, :crmr, :lslq, :lsmr, :lsqr, :dqgmres, :diom, :cgs)
  @eval begin
    $fn(A :: AbstractMatrix{TA}, b :: AbstractVector{Tb}, args...; kwargs...) where {TA, Tb <: Number} =
      $fn(preallocated_LinearOperator(A), convert(Vector{Float64}, b), args...; kwargs...)
  end
end

# Define a symmetric linear operator with preallocation
function preallocated_symmetric_LinearOperator(A)
  (n, m) = size(A)
  Ap = zeros(n)
  return LinearOperator(n, m, true, true, p -> mul!(Ap, A, p))
end

# Variants where A must be symmetric
for fn in (:cg_lanczos, :cg_lanczos_shift_seq, :cg, :cr, :minres, :symmlq)
  @eval begin
    $fn(A :: AbstractMatrix{TA}, b :: AbstractVector{Tb}, args...; kwargs...) where {TA, Tb <: Number} =
      $fn(preallocated_symmetric_LinearOperator(A), convert(Vector{Float64}, b), args...; kwargs...)
  end
end


text:
using Test
import GeoInterface, GDAL, ArchGDAL; const AG = ArchGDAL

@testset "Incomplete GeoInterface geometries" begin
    @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(AG.creategeom(GDAL.wkbCircularString))
    @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(AG.creategeom(GDAL.wkbCompoundCurve))
    @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(AG.creategeom(GDAL.wkbCurvePolygon))
    @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(AG.creategeom(GDAL.wkbMultiSurface))
    @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(AG.creategeom(GDAL.wkbPolyhedralSurface))
    @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(AG.creategeom(GDAL.wkbTIN))
    @test_logs (:warn, "unknown geometry type") GeoInterface.geotype(AG.creategeom(GDAL.wkbTriangle))
end

@testset "Create a Point" begin
    # Method 1
    AG.createpoint(100, 70) do point
        @test GeoInterface.geotype(point) == :Point
        @test isapprox(GeoInterface.coordinates(point), [100,70], atol=1e-6)
        @test AG.getdim(point) == 0
        @test AG.getcoorddim(point) == 2
        AG.setcoorddim!(point, 3)
        @test AG.getcoorddim(point) == 3
        @test AG.isvalid(point) == true
        @test AG.issimple(point) == true
        @test AG.isring(point) == false
        @test AG.getz(point, 0) == 0

        @test sprint(print, AG.getenvelope(point)) == "GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)"
        @test sprint(print, AG.getenvelope3d(point)) == "GDAL.OGREnvelope3D(100.0, 100.0, 70.0, 70.0, 0.0, 0.0)"
        @test AG.toISOWKB(point, GDAL.wkbNDR) == UInt8[0x01,0xe9,0x03,0x00,0x00,
        0x00,0x00,0x00,0x00,0x00,0x00,0x59,0x40,0x00,0x00,0x00,0x00,0x00,0x80,
        0x51,0x40,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00]
        @test AG.toISOWKB(point, GDAL.wkbXDR) == UInt8[0x00,0x00,0x00,0x03,0xe9,
        0x40,0x59,0x00,0x00,0x00,0x00,0x00,0x00,0x40,0x51,0x80,0x00,0x00,0x00,
        0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00]
        @test AG.toKML(point, "relativeToGround") == "<Point><altitudeMode>relativeToGround</altitudeMode><coordinates>100,70,0</coordinates></Point>"
        @test AG.toKML(point, "clampToGround") == "<Point><altitudeMode>clampToGround</altitudeMode><coordinates>100,70,0</coordinates></Point>"
        @test AG.toKML(point) == "<Point><coordinates>100,70,0</coordinates></Point>"
        @test AG.toJSON(point) == "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0, 0.0 ] }"
        AG.createpoint(100,70,0) do point2
            @test isapprox(GeoInterface.coordinates(point2), [100,70,0], atol=1e-6)
            @test AG.equals(point, point2) == true
        end
        AG.createpoint((100,70,0)) do point3
            @test AG.equals(point, point3) == true
        end
        AG.createpoint([100,70,0]) do point4
            @test AG.equals(point, point4) == true
        end
        point5 = AG.createpoint([100,70,0])
        @test AG.equals(point, point5) == true
        AG.flattento2d!(point)
        @test AG.getcoorddim(point) == 2
        @test AG.getnonlineargeomflag() == true
        AG.setnonlineargeomflag!(false)
        @test AG.getnonlineargeomflag() == false
        AG.setnonlineargeomflag!(true)
        @test AG.getnonlineargeomflag() == true
        AG.closerings!(point)
        @test AG.toJSON(point) == "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0 ] }"
    end

    # Method 2
    point = AG.createpoint(100, 70)
    @test AG.getdim(point) == 0
    @test AG.getcoorddim(point) == 2
    AG.setcoorddim!(point, 3)
    @test AG.getcoorddim(point) == 3
    @test AG.isvalid(point) == true
    @test AG.issimple(point) == true
    @test AG.isring(point) == false
    @test AG.getz(point, 0) == 0
    @test sprint(print, AG.getenvelope(point)) == "GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)"
    @test sprint(print, AG.getenvelope3d(point)) == "GDAL.OGREnvelope3D(100.0, 100.0, 70.0, 70.0, 0.0, 0.0)"
    @test AG.toISOWKB(point, GDAL.wkbNDR) == UInt8[0x01,0xe9,0x03,0x00,0x00,
    0x00,0x00,0x00,0x00,0x00,0x00,0x59,0x40,0x00,0x00,0x00,0x00,0x00,0x80,
    0x51,0x40,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00]
    @test AG.toISOWKB(point, GDAL.wkbXDR) == UInt8[0x00,0x00,0x00,0x03,0xe9,
    0x40,0x59,0x00,0x00,0x00,0x00,0x00,0x00,0x40,0x51,0x80,0x00,0x00,0x00,
    0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00]
    @test AG.toKML(point, "relativeToGround") == "<Point><altitudeMode>relativeToGround</altitudeMode><coordinates>100,70,0</coordinates></Point>"
    @test AG.toKML(point, "clampToGround") == "<Point><altitudeMode>clampToGround</altitudeMode><coordinates>100,70,0</coordinates></Point>"
    @test AG.toKML(point) == "<Point><coordinates>100,70,0</coordinates></Point>"
    @test AG.toJSON(point) == "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0, 0.0 ] }"
    @test AG.equals(point, AG.createpoint(100,70,0)) == true
    @test AG.equals(point, AG.createpoint((100,70,0))) == true
    AG.flattento2d!(point)
    @test AG.getcoorddim(point) == 2
    @test AG.getnonlineargeomflag() == true
    AG.setnonlineargeomflag!(false)
    @test AG.getnonlineargeomflag() == false
    AG.setnonlineargeomflag!(true)
    @test AG.getnonlineargeomflag() == true
    AG.closerings!(point)
    @test AG.toJSON(point) == "{ \"type\": \"Point\", \"coordinates\": [ 100.0, 70.0 ] }"
end

@testset "Testing construction of complex geometries" begin
    @test AG.toWKT(AG.createlinestring([1.,2.,3.], [4.,5.,6.])) == "LINESTRING (1 4,2 5,3 6)"
    AG.createlinestring([1.,2.,3.], [4.,5.,6.]) do geom
        @test GeoInterface.geotype(geom) == :LineString
        @test isapprox(GeoInterface.coordinates(geom), [[1,4],[2,5],[3,6]], atol=1e-6)
        @test AG.toWKT(geom) == "LINESTRING (1 4,2 5,3 6)"
        AG.closerings!(geom)
        @test AG.toWKT(geom) == "LINESTRING (1 4,2 5,3 6)"
        AG.setpoint!(geom, 1, 10, 10)
        @test AG.toWKT(geom) == "LINESTRING (1 4,10 10,3 6)"
    end
    AG.createlinestring([1.,2.,3.], [4.,5.,6.], [7.,8.,9.]) do geom
        @test AG.toWKT(geom) == "LINESTRING (1 4 7,2 5 8,3 6 9)"
        AG.setpoint!(geom, 1, 10, 10, 10)
        @test AG.toWKT(geom) == "LINESTRING (1 4 7,10 10 10,3 6 9)"
        AG.addpoint!(geom, 11, 11, 11)
        @test AG.toWKT(geom) == "LINESTRING (1 4 7,10 10 10,3 6 9,11 11 11)"
    end

    @test AG.toWKT(AG.createlinearring([1.,2.,3.], [4.,5.,6.])) == "LINEARRING (1 4,2 5,3 6)"
    AG.createlinearring([1.,2.,3.], [4.,5.,6.]) do geom
        # @test GeoInterface.geotype(geom) == :LinearRing
        @test isapprox(GeoInterface.coordinates(geom), [[1,4],[2,5],[3,6]], atol=1e-6)
        @test AG.toWKT(geom) == "LINEARRING (1 4,2 5,3 6)"
        AG.setpointcount!(geom, 5)
        @test AG.toWKT(geom) == "LINEARRING (1 4,2 5,3 6,0 0,0 0)"
        AG.empty!(geom)
        @test AG.toWKT(geom) == "LINEARRING EMPTY"
    end
    AG.createlinearring([1.,2.,3.], [4.,5.,6.], [7.,8.,9.]) do geom
        @test AG.toWKT(geom) == "LINEARRING (1 4 7,2 5 8,3 6 9)"
        AG.closerings!(geom)
        @test AG.toWKT(geom) == "LINEARRING (1 4 7,2 5 8,3 6 9,1 4 7)"
    end

    @test AG.toWKT(AG.createpolygon([1.,2.,3.], [4.,5.,6.])) == "POLYGON ((1 4,2 5,3 6))"
    AG.createpolygon([1.,2.,3.], [4.,5.,6.]) do geom
        @test GeoInterface.geotype(geom) == :Polygon
        @test isapprox(GeoInterface.coordinates(geom), [[[1,4],[2,5],[3,6]]], atol=1e-6)
        @test AG.toWKT(geom) == "POLYGON ((1 4,2 5,3 6))"
    end
    AG.createpolygon([1.,2.,3.], [4.,5.,6.], [7.,8.,9.]) do geom
        @test AG.toWKT(geom) == "POLYGON ((1 4 7,2 5 8,3 6 9))"
        AG.closerings!(geom)
        @test AG.toWKT(geom) == "POLYGON ((1 4 7,2 5 8,3 6 9,1 4 7))"
    end

    @test AG.toWKT(AG.createmultipoint([1.,2.,3.], [4.,5.,6.])) == "MULTIPOINT (1 4,2 5,3 6)"
    AG.createmultipoint([1.,2.,3.], [4.,5.,6.]) do geom
        @test GeoInterface.geotype(geom) == :MultiPoint
        @test isapprox(GeoInterface.coordinates(geom), [[1,4],[2,5],[3,6]], atol=1e-6)
        @test AG.toWKT(geom) == "MULTIPOINT (1 4,2 5,3 6)"
    end
    AG.createmultipoint([1.,2.,3.], [4.,5.,6.], [7.,8.,9.]) do geom
        @test AG.toWKT(geom) == "MULTIPOINT (1 4 7,2 5 8,3 6 9)"
    end

    @test AG.toWKT(AG.createmultipolygon(Vector{Vector{Tuple{Cdouble,Cdouble}}}[
        Vector{Tuple{Cdouble,Cdouble}}[
            [(0,0),(0,4),(4,4),(4,0)],
            [(1,1),(1,3),(3,3),(3,1)]],
        Vector{Tuple{Cdouble,Cdouble}}[
            [(10,0),(10,4),(14,4),(14,0)],
            [(11,1),(11,3),(13,3),(13,1)]]]
    )) == "MULTIPOLYGON (((0 0,0 4,4 4,4 0),(1 1,1 3,3 3,3 1)),((10 0,10 4,14 4,14 0),(11 1,11 3,13 3,13 1)))"
    AG.createmultipolygon(Vector{Vector{Tuple{Cdouble,Cdouble}}}[
                            Vector{Tuple{Cdouble,Cdouble}}[
                                [(0,0),(0,4),(4,4),(4,0)],
                                [(1,1),(1,3),(3,3),(3,1)]],
                            Vector{Tuple{Cdouble,Cdouble}}[
                                [(10,0),(10,4),(14,4),(14,0)],
                                [(11,1),(11,3),(13,3),(13,1)]]]) do geom
        @test GeoInterface.geotype(geom) == :MultiPolygon
        @test isapprox(
            GeoInterface.coordinates(geom), [
                [[[0,0],[0,4],[4,4],[4,0]], [[1,1],[1,3],[3,3],[3,1]]],
                [[[10,0],[10,4],[14,4],[14,0]], [[11,1],[11,3],[13,3],[13,1]]]
            ],
            atol=1e-6
        )
        @test AG.toWKT(geom) == "MULTIPOLYGON (((0 0,0 4,4 4,4 0),(1 1,1 3,3 3,3 1)),((10 0,10 4,14 4,14 0),(11 1,11 3,13 3,13 1)))"
    end

    AG.fromWKT("CURVEPOLYGON(CIRCULARSTRING(-2 0,-1 -1,0 0,1 -1,2 0,0 2,-2 0),(-1 0,0 0.5,1 0,0 1,-1 0))") do geom
        @test AG.toWKT(AG.getcurvegeom(AG.getlineargeom(geom, 0.5))) == "CURVEPOLYGON (CIRCULARSTRING (-2 0,-1 -1,0 0,1 -1,2 0,0 2,-2 0),(-1 0,0.0 0.5,1 0,0 1,-1 0))"
        AG.getlineargeom(geom, 0.5) do lgeom
            AG.getcurvegeom(lgeom) do clgeom
                @test AG.toWKT(clgeom) == "CURVEPOLYGON (CIRCULARSTRING (-2 0,-1 -1,0 0,1 -1,2 0,0 2,-2 0),(-1 0,0.0 0.5,1 0,0 1,-1 0))"
            end
            @test AG.ngeom(AG.polygonize(AG.forceto(lgeom, GDAL.wkbMultiLineString))) == 2
            AG.forceto(lgeom, GDAL.wkbMultiLineString) do mlsgeom
                AG.polygonize(mlsgeom) do plgeom
                    @test AG.ngeom(plgeom) == 2
                end
            end
        end
    end
end

# TODO: fix #47
# @testset "Testing remaining methods for geometries" begin
#     AG.createmultipolygon(Vector{Vector{Tuple{Cdouble,Cdouble}}}[
#                             Vector{Tuple{Cdouble,Cdouble}}[
#                                 [(0,0),(0,4),(4,4),(4,0)],
#                                 [(1,1),(1,3),(3,3),(3,1)]],
#                             Vector{Tuple{Cdouble,Cdouble}}[
#                                 [(10,0),(10,4),(14,4),(14,0)],
#                                 [(11,1),(11,3),(13,3),(13,1)]]]) do geom1
#     AG.createmultipoint([1.,2.,3.], [4.,5.,6.], [7.,8.,9.]) do geom2
#         AG.closerings!(geom1)
#         @test AG.disjoint(geom1, geom2) == false
#         @test AG.touches(geom1, geom2) == true
#         @test AG.crosses(geom1, geom2) == false
#         @test AG.overlaps(geom1, geom2) == false

#         @test AG.toWKT(AG.boundary(geom2)) == "GEOMETRYCOLLECTION EMPTY"
#         AG.boundary(geom2) do result
#             @test AG.toWKT(result) == "GEOMETRYCOLLECTION EMPTY"
#         end

#         @test AG.toWKT(AG.union(geom1, geom2)) == "GEOMETRYCOLLECTION (POINT (2 5 8),POINT (3 6 9),POLYGON ((0 0 8,0 4 8,4 4 8,4 0 8,0 0 8),(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8)),POLYGON ((10 0 8,10 4 8,14 4 8,14 0 8,10 0 8),(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8)))"
#         AG.union(geom1, geom2) do result
#             @test AG.toWKT(result) == "GEOMETRYCOLLECTION (POINT (2 5 8),POINT (3 6 9),POLYGON ((0 0 8,0 4 8,4 4 8,4 0 8,0 0 8),(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8)),POLYGON ((10 0 8,10 4 8,14 4 8,14 0 8,10 0 8),(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8)))"
#             @test AG.hascurvegeom(result, true) == false
#             @test AG.hascurvegeom(result, false) == false
#         end

#         @test AG.toWKT(AG.difference(geom1, geom2)) == "MULTIPOLYGON (((0 0 8,0 4 8,4 4 8,4 0 8,0 0 8),(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8)),((10 0 8,10 4 8,14 4 8,14 0 8,10 0 8),(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8)))"
#         AG.difference(geom1, geom2) do result
#             @test AG.toWKT(result) == "MULTIPOLYGON (((0 0 8,0 4 8,4 4 8,4 0 8,0 0 8),(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8)),((10 0 8,10 4 8,14 4 8,14 0 8,10 0 8),(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8)))"
#             AG.segmentize!(result, 20)
#             @test AG.toWKT(result) == "MULTIPOLYGON (((0 0 8,0 4 8,4 4 8,4 0 8,0 0 8),(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8)),((10 0 8,10 4 8,14 4 8,14 0 8,10 0 8),(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8)))"
#             AG.segmentize!(result, 2)
#             @test AG.toWKT(result) == "MULTIPOLYGON (((0 0 8,0.0 1.33333333333333 8,0.0 2.66666666666667 8,0 4 8,1.33333333333333 4.0 8,2.66666666666667 4.0 8,4 4 8,4.0 2.66666666666667 8,4.0 1.33333333333333 8,4 0 8,2.66666666666667 0.0 8,1.33333333333333 0.0 8,0 0 8),(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8)),((10 0 8,10.0 1.33333333333333 8,10.0 2.66666666666667 8,10 4 8,11.3333333333333 4.0 8,12.6666666666667 4.0 8,14 4 8,14.0 2.66666666666667 8,14.0 1.33333333333333 8,14 0 8,12.6666666666667 0.0 8,11.3333333333333 0.0 8,10 0 8),(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8)))"
#         end

#         @test AG.toWKT(AG.symdifference(geom1, geom2)) == "GEOMETRYCOLLECTION (POINT (2 5 8),POINT (3 6 9),POLYGON ((0 0 8,0 4 8,4 4 8,4 0 8,0 0 8),(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8)),POLYGON ((10 0 8,10 4 8,14 4 8,14 0 8,10 0 8),(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8)))"
#         AG.symdifference(geom1, geom2) do result
#             @test GeoInterface.geotype(result) == :GeometryCollection
#             @test AG.toWKT(result) == "GEOMETRYCOLLECTION (POINT (2 5 8),POINT (3 6 9),POLYGON ((0 0 8,0 4 8,4 4 8,4 0 8,0 0 8),(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8)),POLYGON ((10 0 8,10 4 8,14 4 8,14 0 8,10 0 8),(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8)))"
#             AG.removegeom!(result, 1)
#             @test AG.toWKT(result) == "GEOMETRYCOLLECTION (POINT (2 5 8),POLYGON ((0 0 8,0 4 8,4 4 8,4 0 8,0 0 8),(1 1 8,3 1 8,3 3 8,1 3 8,1 1 8)),POLYGON ((10 0 8,10 4 8,14 4 8,14 0 8,10 0 8),(11 1 8,13 1 8,13 3 8,11 3 8,11 1 8)))"
#             AG.removeallgeoms!(result)
#             @test AG.toWKT(result) == "GEOMETRYCOLLECTION EMPTY"
#         end
#     end
#     end
# end

# Untested
# toISOWKT(geom::Geometry)
# unsafe_polygonfromedges(lines::Geometry, besteffort::Bool,autoclose::Bool, tol::Real)
# setspatialref!(geom::Geometry, spatialref::SpatialRef)
# getspatialref(geom::Geometry) = GDAL.getspatialreference(geom)
# transform!(geom::Geometry, coordtransform::CoordTransform)
# transform!(geom::Geometry, spatialref::SpatialRef)


text:
        # union
        r1 = AABox2D(Point2D{T}(0, 0), Point2D{T}(2, 2))
        r2 = AABox2D(Point2D{T}(1, 1), Point2D{T}(4, 4))
        r3 = r1 ∪ r2
        @test r3.xmin ≈ 0
        @test r3.ymin ≈ 0
        @test r3.xmax ≈ 4
        @test r3.ymax ≈ 4
        r1 = AABox2D(Point2D{T}(0, 0), Point2D{T}(2, 2))
        r2 = AABox2D(Point2D{T}(1, 1), Point2D{T}(3, 3))
        r3 = AABox2D(Point2D{T}(2, 2), Point2D{T}(4, 4))
        r4 = union([r1,r2,r3])
        @test r4.xmin ≈ 0
        @test r4.ymin ≈ 0
        @test r4.xmax ≈ 4
        @test r4.ymax ≈ 4

#        # union
#        r1 = AABox3D(Point3D{T}(0, 0), Point3D{T}(2, 2))
#        r2 = AABox3D(Point3D{T}(1, 1), Point3D{T}(4, 4))
#        r3 = r1 ∪ r2
#        @test r3.xmin ≈ 0
#        @test r3.ymin ≈ 0
#        @test r3.xmax ≈ 4
#        @test r3.ymax ≈ 4

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport