text:
function utils_unit_tests()
    generate_symmetric_distance_matrix_tests()
    set_indices_tests()
    build_computed_data_tests()
    create_singleton_location_groups_tests()
    preprocess_instance_unit_tests()
    create_intersection_tws_tests()
end

function create_intersection_tws_tests()
    tws_1 = [RVRP.Range(10.0, 20.0), RVRP.Range(35.0, 40.0)]
    tws_2 = [RVRP.Range(15.0, 25.0), RVRP.Range(28.0, 36.0)]
    intersection = RVRP.ranges_intersection(tws_1, tws_2)
    @test intersection == [RVRP.Range(15.0, 20.0), RVRP.Range(35.0, 36.0)]

    tws_1 = [RVRP.Range(10.0, 100.0), RVRP.Range(110.0, 130.0)]
    tws_2 = [RVRP.Range(9.0, 25.0), RVRP.Range(28.0, 36.0), RVRP.Range(88.0, 120.0)]
    intersection = RVRP.ranges_intersection(tws_1, tws_2)
    @test intersection == [RVRP.Range(10.0, 25.0), RVRP.Range(28.0, 36.0), RVRP.Range(88.0, 100.0), RVRP.Range(110.0, 120.0)]
end

function generate_symmetric_distance_matrix_tests()
    xs = [rand(1:20) for i in 1:5]
    ys = [rand(1:20) for i in 1:5]
    matrix = RVRP.generate_symmetric_distance_matrix(xs, ys)
    @test size(matrix) == (5,5)
    for i in 1:5
        @test matrix[i,i] == 0.0
        for j in i+1:5
            @test matrix[i,j] == matrix[j,i]
        end
    end
end

function set_indices_tests()
    data = RVRP.RvrpInstance()
    data.locations = [RVRP.Location(
        id = "", index = rand(-100:-1)
    ) for i in 1:10]
    RVRP.set_indices(data)
    for loc_idx in 1:length(data.locations)
        @test data.locations[loc_idx].index == loc_idx
    end
end

function build_computed_data_tests()
    data = RVRP.RvrpInstance()
    data.locations = [RVRP.Location(
        id = string("loc_", i), index = rand(-100:-1)
    ) for i in 1:9]
    data.vehicle_categories = [RVRP.VehicleCategory(
        id = string("vc_", i)
    ) for i in 1:9]
    data.vehicle_sets = [RVRP.HomogeneousVehicleSet(
        id = string("vs_", i)
    ) for i in 1:12]
    data.product_specification_classes = [RVRP.ProductSpecificationClass(
        id = string("psc_", i)
    ) for i in 1:5]
    data.requests = [RVRP.Request(
        id = string("req_", i)
    ) for i in 10:13]
    RVRP.set_indices(data)
    data.location_groups = RVRP.create_singleton_location_groups(data.locations)
    data.vehicle_categories[1].capacities.of_vehicle = Dict{String,Float64}(
        "weird_name_1" => 10, "wolow_34" => 12, "volume" => 15
    )
    data.vehicle_categories[2].capacities.of_vehicle = Dict{String,Float64}(
        "weird_name_1" => 10, "wolow_34" => 12, "new_name_1" => 15
    )
    data.product_specification_classes[1].capacity_consumptions = Dict{String,Tuple{Float64,Float64}}(
        "weird_name_1" => (1,10), "new_name_2" => (1,12), "wolow_34" => (1,15)
    )
    data.vehicle_categories[1].properties.of_vehicle = Dict{String,Float64}(
        "prop_1" => 10, "prop_34" => 12, "prop_3" => 15
    )
    data.vehicle_categories[2].properties.of_vehicle = Dict{String,Float64}(
        "prop_1" => 10, "prop_34" => 12, "new_prop_1" => 15
    )
    data.product_specification_classes[1].property_requirements = Dict{String,Float64}(
        "prop_1" => 10, "new_prop_2" => 12, "prop_34" => 15
    )
    computed_data = RVRP.build_computed_data(data)
    for (k,v) in computed_data.location_id_2_index
        @test v >= 1
        @test v <= 9
    end
    for (k,v) in computed_data.location_group_id_2_index
        @test v >= 1
        @test v <= 9
    end
    for (k,v) in computed_data.vehicle_category_id_2_index
        @test v >= 1
        @test v <= 9
    end
    for (k,v) in computed_data.vehicle_set_id_2_index
        @test v >= 1
        @test v <= 12
    end
    for (k,v) in computed_data.product_specification_class_id_2_index
        @test v >= 1
        @test v <= 5
    end
    for (k,v) in computed_data.request_id_2_index
        @test v >= 1
        @test v <= 4
    end
    for (k,v) in computed_data.capacity_id_2_index
        @test v >= 1
        @test v <= 5
    end
    for (k,v) in computed_data.property_id_2_index
        @test v >= 1
        @test v <= 5
    end
end

function create_singleton_location_groups_tests()
    data = RVRP.RvrpInstance()
    data.locations = [RVRP.Location(
        index = rand(-100:-1)
    ) for i in 1:10]
    def_loc_groups = RVRP.create_singleton_location_groups(data.locations)
    for i in 1:length(data.locations)
        @test string(data.locations[i].id, "_loc_group") == def_loc_groups[i].id
    end
end

function preprocess_instance_unit_tests()
    data = RVRP.RvrpInstance()
    RVRP.preprocess_instance(data)
    @test length(data.locations) == 0
    @test length(data.location_groups) == 0
    @test length(data.product_compatibility_classes) == 1
    @test data.product_compatibility_classes[1].id == "default_id"
    @test length(data.product_sharing_classes) == 1
    @test data.product_sharing_classes[1].id == "default_id"
    @test length(data.product_specification_classes) == 1
    @test data.product_specification_classes[1].id == "default_id"
    @test length(data.requests) == 0
    @test length(data.vehicle_categories) == 1
    @test data.vehicle_categories[1].id == "default_id"
    @test length(data.vehicle_sets) == 0
end


text:
using IntervalArithmetic
using .DynamicDefinition

"""
Return Branches for a given dynamic, in an iterable
"""
function branches(D::PwMap)
    return [Branch(D.Ts[k], (D.endpoints[k], D.endpoints[k+1]), D.is_full[k] ? (Interval(0),Interval(1)) : (D.Ts[k](Interval(D.endpoints[k]), D.Ts[k](Interval(D.endpoints[k+1])))), D.orientations[k]==1) for k in 1:length(D.Ts)]
end

"""
Return preimages of a certain sequence on all branches of a dynamic
"""
function preimages(seq, D::Dynamic, ϵ = 0.0)
    return [preimages(seq, branch, ϵ) for branch in branches(D)]
end
# TODO: overload for IterateDynamic

"""
Constructs associated "Ulam duals" for a branch b,
i.e., the sequence (k, (T⁻¹(p[k]), T⁻¹(p[k+1]))),
handling endpoints correctly. A callback `f(k, (a, b))` is called on each dual.

This should eventually replace `DualComposedWithDynamic`.
"""
function callback_duals(f, B::Ulam, branch::Branch, preims=nothing)
    if preims === nothing
        preims = preimages(PointSequence(B.p), b)
    end
    n = length(preims.v)

    # duals = Tuple{Int, Tuple{eltype(preims), eltype(preims)}}[]
    # sizehint!(duals, n+2)

    if n==0 # special case: only one preimage
        first_endpoint = preims.increasing ? branch.X[1] : branch.X[2]
        last_endpoint = preims.increasing ? branch.X[2] : branch.X[1]

        f(preims.skip, (first_endpoint, last_endpoint))
        return nothing
    end

    if preims.skip > 0
        first_endpoint = preims.increasing ? branch.X[1] : branch.X[2]
        if first_endpoint != preims.v[1]
            f(preims.skip, (first_endpoint, preims.v[1]))
        end
    end
    for k in 1:n-1
        f(preims.skip+k, (preims.v[k], preims.v[k+1]))
    end
    if preims.skip + n < length(B.p)  # if there are skipped entries at the end of preims
        last_endpoint = preims.increasing ? b.X[2] : b.X[1]
        if last_endpoint != preims.v[n]
            f(preims.skip+n, (preims.v[n], last_endpoint))
        end
    end
    return nothing
end

function callback_duals(f, B::Hat, branch::Branch, preims=nothing)
    if preims === nothing
        preims = preimages(PointSequence(B.p), b)
    end
    @assert preims.skip == 0
    @assert length(preims.v) == length(B.p) # The Hat basis estimation works only for full-branch maps
    for i = 1:length(B) # this skips p[n+1]==1
        x = preims.v[i]
        absT′ = abs(derivative(branch.f, x))
        f(i, (x, absT′))
    end
end

function callback_duals(f, B::Basis, D::Dynamic)
    for (branch, preim) in zip(branches(D), preimages(PointSequence(B.p), D))
        callback_duals(f, B, branch, preim)
    end
end

function assemble2(B::Basis, D::Dynamic, ϵ=0.0; T = Float64)
    # putting types here in hope to improve callback inference
	I::Vector{Int64} = Int64[]
	J::Vector{Int64} = Int64[]
	nzvals::Vector{Interval{T}} = Interval{T}[]
	n::Int64 = length(B)

	# TODO: reasonable size hint?

    function assemble_callback(i, dual_element)
		if !is_dual_element_empty(B, dual_element)
			for (j, x) in ProjectDualElement(B, dual_element)
				push!(I, i)
				push!(J, mod(j,1:n))
				push!(nzvals, x)
			end
		end
    end
    callback_duals(assemble_callback, B, D)

	return sparse(I, J, nzvals, n, n)
end

ext: jl


text:
module test

using VegaLite
include("graph_functions.jl")

end


text:
module Static

export render
import HttpServer: Response, mimetypes

function mime(s)
    if in(s, mimetypes |> keys)
        return mimetypes[s]
    end
    "text/plain"
end

"""
# public serve some static files?
"""
function public(public_directory)
    function (app, ctx)
        ctx.data[:public_dir] = public_directory
        resourcefile = joinpath(public_directory,
                                join(ctx.data[:request][:path], "/")) |> abspath

        if isfile(resourcefile)
            ext = split(resourcefile |> basename, ".")[end]
            f = open(resourcefile)
            res = (join(readall(resourcefile), ""))
            close(f)
            ctx.response.headers["Content-Type"] = mime(ext) * "; charset=utf-8"
            return res
        end

        app(ctx)
    end
end

"""
# favicon

for this to work, you need to have the public middleware installed
"""
function favicon(app, ctx)
    if join(ctx.data[:request][:path], "/") == "favicon.ico"
        return "awesome" # this should be the location of a favicon
    end

    app(ctx)
end

end


text:
module Fields

export Field
export interior, zero_fields!
export fill_halo_nut!, fill_halo_Gcs!, fill_halo_vel!
export fill_halo_u!, fill_halo_v!, fill_halo_w!
export apply_bcs!, getbc
export default_bcs
export set_bc!, validate_bcs
export nut_names

using KernelAbstractions
using CUDA

using PlanktonIndividuals.Grids
using PlanktonIndividuals.Architectures: device, Architecture, array_type

include("halo_regions.jl")
include("boundary_conditions.jl")
include("apply_bcs.jl")

struct Field
    data::AbstractArray{Float64,3}
    bc::BoundaryConditions
end

"""
    Field(arch::Architecture, grid::AbstractGrid; bcs = default_bcs())
Construct a `Field` on `grid` with data and boundary conditions on architecture `arch`
"""
function Field(arch::Architecture, grid::AbstractGrid; bcs = default_bcs())
    total_size = (grid.Nx+grid.Hx*2, grid.Ny+grid.Hy*2, grid.Nz+grid.Hz*2)
    data = zeros(total_size) |> array_type(arch)
    return Field(data,bcs)
end

@inline interior(c, grid) = c[grid.Hx+1:grid.Hx+grid.Nx, grid.Hy+1:grid.Hy+grid.Ny, grid.Hz+1:grid.Hz+grid.Nz]

function zero_fields!(a)
    for i in 1:length(a)
        @inbounds a[i].data .= 0.0
    end
end

const nut_names=(:DIC,:NH4,:NO3,:PO4,:DOC,:DON,:DOP,:POC,:PON,:POP)

end


text:
module NewtonsMethod

using LinearAlgebra, Statistics, Compat, ForwardDiff

function newtonroot(f, f′; x₀, tolerance = 1E-7, maxiter = 1000)
    # setup the algorithm
    x_old = x₀
    normdiff = Inf
    iter = 1
    while normdiff > tolerance && iter <= maxiter
        x_new = x_old - f(x_old)/f′(x_old) # use the passed in map
        normdiff = norm(x_new - x_old)
        x_old = x_new
        iter = iter + 1
    end

    if iter < maxiter+1
        return (root = x_old, normdiff = normdiff, iter = iter) # A named tuple
    else
        return nothing
    end
end

# for auto-differentiation, can just forward to the existing one
D(f) = x -> ForwardDiff.derivative(f, x)

# same name! different parameters
newtonroot(f; x₀, tolerance = 1E-7, maxiter = 1000) = newtonroot(f, D(f), x₀=x₀, tolerance = tolerance, maxiter = maxiter)

export newtonroot
end # module


text:
_is_incremental(::Type{Cbc.Optimizer}) = false

function set_default_config!(ext::ExtensionType, d::GlobalOptimizer, m::Incremental{Cbc.Optimizer}, local_solver::Bool)
    #MOI.set(m, MOI.RawOptimizerAttribute("allowableGap"), _absolute_tol(d)*1E-2)
    #MOI.set(m, MOI.RawOptimizerAttribute("ratioGap"),     _absolute_tol(d)*1E-2)
    #MOI.set(m, MOI.RawOptimizerAttribute("logLevel"),     0)
    #MOI.set(m, MOI.RawOptimizerAttribute("threads"), Threads.nthreads())
    return
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
products = [
    LibraryProduct(prefix, ["libcrypto", "libcrypto-1_1", "libcrypto-1_1-x64"], :libcrypto),
    LibraryProduct(prefix, ["libssl", "libssl-1_1", "libssl-1_1-x64"], :libssl),
    ExecutableProduct(prefix, "openssl", :openssl),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaBinaryWrappers/OpenSSL_jll.jl/releases/download/OpenSSL-v1.1.1+2"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/OpenSSL.v1.1.1.aarch64-linux-gnu.tar.gz", "deae2e855536f8ca336870450c555335db94647e93d7150134269ee234b2f004"),
    Linux(:aarch64, libc=:musl) => ("$bin_prefix/OpenSSL.v1.1.1.aarch64-linux-musl.tar.gz", "385eedaafd409abf4114ee0b33c70a89145695984282c20ca7c4cede3ed62260"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf) => ("$bin_prefix/OpenSSL.v1.1.1.armv7l-linux-gnueabihf.tar.gz", "00e1413c4d5f022b20d8824eeb6e1d50b90cc938a2c5e179de64f815b2980e66"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf) => ("$bin_prefix/OpenSSL.v1.1.1.armv7l-linux-musleabihf.tar.gz", "c0ba637c8301834a0357aac7ef3cfe6c60254b95d62d06de7802dc2441185c66"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/OpenSSL.v1.1.1.i686-linux-gnu.tar.gz", "844a452a6a5dea681bb60dde94c67a880ae96aec97bbc5e111310a5be5bdbf89"),
    Linux(:i686, libc=:musl) => ("$bin_prefix/OpenSSL.v1.1.1.i686-linux-musl.tar.gz", "fd4beae0cb0d1cd17923985a219c1569346b7b32f5c579746919cb81971e37fb"),
    Windows(:i686) => ("$bin_prefix/OpenSSL.v1.1.1.i686-w64-mingw32.tar.gz", "3c942cd0a69ab5a4cef72042d2ba1562b92ceb58112a1d62ae8d324ff603da52"),
    Linux(:powerpc64le, libc=:glibc) => ("$bin_prefix/OpenSSL.v1.1.1.powerpc64le-linux-gnu.tar.gz", "07f1161572ccf1beaf1f7affb9996d6387f7417a4a7d8b128e4efaf7812f1927"),
    MacOS(:x86_64) => ("$bin_prefix/OpenSSL.v1.1.1.x86_64-apple-darwin14.tar.gz", "72a73b0015dd3116186197b3adeb6dfd1715974aff7cccc9cc0742f80f07678e"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/OpenSSL.v1.1.1.x86_64-linux-gnu.tar.gz", "77be5b30230b67396a691ba8be144824b515e8ea5070df7bd5b06026fdc5ca56"),
    Linux(:x86_64, libc=:musl) => ("$bin_prefix/OpenSSL.v1.1.1.x86_64-linux-musl.tar.gz", "5c40cef594c4126885e5409fcb5c18a354dd437d86c77f6be631ac704958e254"),
    FreeBSD(:x86_64) => ("$bin_prefix/OpenSSL.v1.1.1.x86_64-unknown-freebsd11.1.tar.gz", "6736ce1e2daa69ff443675ae3f784433b58bcebac3d084e510e44bba7f0fc026"),
    Windows(:x86_64) => ("$bin_prefix/OpenSSL.v1.1.1.x86_64-w64-mingw32.tar.gz", "c38348da07a449a96b767f782b053b150d7837aa0c7bd94a3f6f632886b26d93"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps_OpenSSL.jl"), products, verbose=verbose)


text:
# ____________________________________________________________________
# Cutoff powerlaw
#
mutable struct cutoff_powerlaw <: AbstractComponent
    norm::Parameter
    x0::Parameter
    alpha::Parameter
    beta::Parameter # when beta is positive the cutoff affects x > x0 (and viceversa)

    function cutoff_powerlaw(x0::Number)
        out = new(Parameter(1),
                  Parameter(x0),
                  Parameter(-1),
                  Parameter(1))
        out.norm.low = 0
        out.x0.low = 0
        out.alpha.low = -5
        out.alpha.high = 5
        out.beta.low = 0.1
        out.beta.high = 10.
        return out
    end
end

function evaluate!(buffer, comp::cutoff_powerlaw, x::Domain{1},
                   norm, x0, alpha, beta)
    buffer .= norm .* (x ./ x0).^alpha .* exp.(1 .- ((x ./ x0) .^beta))
end


text:
using ReachabilityAnalysis, SparseArrays

LazySets.set_ztol(Float64, 1e-14)

function platoon_connected(; deterministic_switching::Bool=true,
                             c1=5.0)  # clock constraints
    n = 10 # 9 dimensions + time
    # x' = Ax + Bu + c
    A = Matrix{Float64}(undef, n, n)
    A[1, :] = [0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0]
    A[2, :] = [0, 0, -1.0, 0, 0, 0, 0, 0, 0, 0]
    A[3, :] = [1.6050, 4.8680, -3.5754, -0.8198, 0.4270, -0.0450, -0.1942,  0.3626, -0.0946, 0.]
    A[4, :] = [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0,]
    A[5, :] = [0, 0, 1.0, 0, 0, -1.0, 0, 0, 0, 0]
    A[6, :] = [0.8718, 3.8140, -0.0754,  1.1936, 3.6258, -3.2396, -0.5950,  0.1294, -0.0796, 0.]
    A[7, :] = [0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0]
    A[8, :] = [0, 0, 0, 0, 0, 1.0, 0, 0, -1.0, 0]
    A[9, :] = [0.7132, 3.5730, -0.0964,  0.8472, 3.2568, -0.0876,  1.2726,  3.0720, -3.1356, 0.]
    A[10, :] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0]; # t' = 1

    if deterministic_switching
        invariant = HalfSpace(sparsevec([n], [1.], n), c1) # t <= c1
    else
        invariant = Universe(n)
    end

    # acceleration of the lead vehicle + time
    B = sparse([2], [1], [1.0], n, 1)
    U = Hyperrectangle(low=[-9.], high=[1.])
    c = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0]
    @system(x' = Ax + Bu + c, x ∈ invariant, u ∈ U)
end

function platoon_disconnected(; deterministic_switching::Bool=true,
                                c2=5.0)  # clock constraints
    n = 10 # 9 dimensions + time
    # x' = Ax + Bu + c
    A = Matrix{Float64}(undef, n, n)
    A[1, :] = [0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0]
    A[2, :] = [0, 0, -1.0, 0, 0, 0, 0, 0, 0, 0]
    A[3, :] = [1.6050, 4.8680, -3.5754, 0, 0, 0, 0, 0, 0, 0]
    A[4, :] = [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0,]
    A[5, :] = [0, 0, 1.0, 0, 0, -1.0, 0, 0, 0, 0]
    A[6, :] = [0, 0, 0,  1.1936, 3.6258, -3.2396, 0, 0, 0, 0.]
    A[7, :] = [0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0]
    A[8, :] = [0, 0, 0, 0, 0, 1.0, 0, 0, -1.0, 0]
    A[9, :] = [0.7132, 3.5730, -0.0964,  0.8472, 3.2568, -0.0876,  1.2726,  3.0720, -3.1356, 0.]
    A[10, :] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0]; # t' = 1

    if deterministic_switching
        invariant = HalfSpace(sparsevec([n], [1.], n), c2) # t <= c2
    else
        invariant = Universe(n)
    end

    # acceleration of the lead vehicle + time
    B = sparse([2], [1], [1.0], n, 1)
    U = Hyperrectangle(low=[-9.], high=[1.])
    c = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0]
    @system(x' = Ax + Bu + c, x ∈ invariant, u ∈ U)
end

function platoon(; deterministic_switching::Bool=true,
                   c1=5.0,  # clock constraints
                   c2=5.0,  # clock constraints
                   tb=10.0,  # lower bound for loss of communication
                   tc=20.0, tr=20.0) # upper bound for loss of communication (tc) and reset time (tr)

    # three variables for each vehicle, (ei, d(et)/dt, ai) for
    # (spacing error, relative velocity, speed), and the last dimension is time
    n = 9 + 1

    # transition graph
    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 2, 1)
    add_transition!(automaton, 2, 1, 2)

    # modes
    mode1 = platoon_connected(deterministic_switching=deterministic_switching, c1=c1)
    mode2 = platoon_disconnected(deterministic_switching=deterministic_switching, c2=c2)
    modes = [mode1, mode2]

    # common reset
    reset = Dict(n => 0.)

    # transition l1 -> l2
    if deterministic_switching
        guard = Hyperplane(sparsevec([n], [1.], n), c1) # t == c1
    else
        # tb <= t <= tc
        guard = HPolyhedron([HalfSpace(sparsevec([n], [-1.], n), -tb),
                             HalfSpace(sparsevec([n], [1.], n), tc)])
    end
    t1 = ConstrainedResetMap(n, guard, reset)

    # transition l2 -> l1
    if deterministic_switching
        guard = Hyperplane(sparsevec([n], [1.], n), c2) # t == c2
    else
        guard = HalfSpace(sparsevec([n], [1.], n), tr) # t <= tr
    end
    t2 = ConstrainedResetMap(n, guard, reset)
    resetmaps = [t1, t2]

    H = HybridSystem(automaton, modes, resetmaps, [AutonomousSwitching()])

    # initial condition is at the orgin in mode 1
    X0 = BallInf(zeros(n), 0.0)
    initial_condition = [(1, X0)]

    return IVP(H, initial_condition)
end

function dmin_specification(sol, dmin)
    return (-ρ(sparsevec([1], [-1.0], 10), sol) > dmin) &&
           (-ρ(sparsevec([4], [-1.0], 10), sol) > dmin) &&
           (-ρ(sparsevec([7], [-1.0], 10), sol) > dmin)
end


text:
module TestDTM
    using Base.Test
    using Languages
    using TextAnalysis

    doc1 = "a a a sample text text"
    doc2 = "another example example text text"
    doc3 = "another another text text text text"

    # TODO: this should work!
    # crps = Corpus(map(StringDocument, [doc1 doc2 doc3]))

    crps = Corpus(Any[StringDocument(doc1), StringDocument(doc2), StringDocument(doc3)])

    update_lexicon!(crps)
    m = DocumentTermMatrix(crps)

    # Terms are in alphabetical ordering
    correctweights = [0.5493061443340548 0.0 0.0 0.18310204811135158 0.0
                      0.0 0.08109302162163289 0.43944491546724385 0.0 0.0
                      0.0 0.13515503603605478 0.0 0.0 0.0]

    myweights = tf_idf(m)
    @test_approx_eq myweights correctweights

    myweights = tf_idf(dtm(m))
    @test_approx_eq myweights correctweights
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = tf_idf(dtm(m, :dense))
    @test_approx_eq myweights correctweights
    @test typeof(myweights) <: Matrix

    myweights = float(dtm(m));
    tf_idf!(myweights)
    @test_approx_eq myweights correctweights
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = float(dtm(m, :dense));
    tf_idf!(myweights)
    @test_approx_eq myweights correctweights
    @test typeof(myweights) <: Matrix
end


text:
export @UNION, @DEFAULT, @ALIGN, @STRUCT, @with_kw

const __module__ = 0

function indexof(needle, haystack)
	for (i, v) in enumerate(haystack)
		v == needle && return i-1
	end
	return -1
end

macro UNION(T, TT)
	typeof(T) == Symbol || throw(ArgumentError("1st argument must be a symbol to represent a Union type"))
	TT.head == :tuple || throw(ArgumentError("2nd argument must be a tuple of types like `(T1,T2,...)`"))
	return esc(quote
		const $T = $(Expr(:curly, :Union, TT.args...))
		FlatBuffers.typeorder(::Type{$T}, ::Type{TT}) where {TT} = FlatBuffers.indexof(TT, $TT)
		FlatBuffers.typeorder(::Type{$T}, i::Integer) = ($TT)[i+1]
		FlatBuffers.isunionwithnothing(::Type{$T}) = false
	end)
end

macro ALIGN(T, sz)
	return esc(quote
		FlatBuffers.alignment(::Type{$T}) = $sz
	end)
end

macro enumtype(T, typ)
	return esc(quote
		FlatBuffers.enumtype(::Type{$T}) = $typ
	end)
end

# recursively finds largest field of a STRUCT
fbsizeof(::Type{T}) where {T<:Enum} = sizeof(enumtype(T))
fbsizeof(::Type{T}) where {T} = sizeof(T)

maxsizeof(::Type{T}) where {T<:Enum} = sizeof(enumtype(T))
maxsizeof(::Type{T}) where {T} = isbitstype(T) ? sizeof(T) : maximum(map(x->maxsizeof(x), T.types))

nextsizeof(::Type{T}) where {T} = isbitstype(T) ? sizeof(T) : nextsizeof(T.types[1])

function fieldlayout(mod, typ, exprs...)
	fields = Expr[]
	values = []
	largest_field = maximum(map(x->maxsizeof(Core.eval(mod, x.args[2])), exprs))
	sz = cur_sz = 0
	x = 0
	for (i,expr) in enumerate(exprs)
		T = Core.eval(mod, expr.args[2])
		if !isbitstype(T)
			exprs2 = [Expr(:(::), nm, typ) for (nm,typ) in zip(fieldnames(T),T.types)]
			fields2, values2 = fieldlayout(mod, T, exprs2...)
			append!(fields, map(x->Expr(:(::), Symbol(string(expr.args[1],'_',x.args[1])), x.args[2]), fields2))
			append!(values, map(x->x == 0 ? 0 : Expr(:call, :getfield, expr.args[1], QuoteNode(x)), values2))
		else
			push!(fields, expr)
			push!(values, expr.args[1])
		end
		sz += cur_sz = fbsizeof(T)
		if sz % largest_field == 0
			sz = cur_sz = 0
			continue
		end
		nextsz = i == length(exprs) ? 0 : nextsizeof(Core.eval(mod, exprs[i+1].args[2]))
		if i == length(exprs) || cur_sz < nextsz || (sz + nextsz) > largest_field
			# this is the last field and we're not `sz % largest_field`
			# potential diffs = 7, 6, 5, 4, 3, 2, 1
			sym = expr.args[1]
			diff = cur_sz < nextsz ? nextsz - cur_sz : largest_field - sz
			if diff == 7
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt8)); x += 1
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt16)); x += 1
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt32)); x += 1
				push!(values, 0); push!(values, 0); push!(values, 0)
			elseif diff == 6
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt16)); x += 1
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt32)); x += 1
				push!(values, 0); push!(values, 0)
			elseif diff == 5
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt8)); x += 1
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt32)); x += 1
				push!(values, 0); push!(values, 0)
			elseif diff == 4
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt32)); x += 1
				push!(values, 0)
			elseif diff == 3
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt8)); x += 1
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt16)); x += 1
				push!(values, 0); push!(values, 0)
			elseif diff == 2
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt16)); x += 1
				push!(values, 0)
			elseif diff == 1
				push!(fields, Expr(:(::), Symbol("_pad_$(sym)_$(typ)_$x"), :UInt8)); x += 1
				push!(values, 0)
			end
			sz = (sz + diff) % largest_field == 0 ? 0 : (cur_sz < nextsz ? sz + diff : 0)
			cur_sz = 0
		end
	end
	return fields, values
end

linefilter = x->typeof(x) != LineNumberNode

macro STRUCT(expr)
	!expr.args[1] || throw(ArgumentError("@STRUCT is only applicable for immutable types"))
	exprs = filter(linefilter, expr.args[3].args)
	fields, values = FlatBuffers.fieldlayout(__module__, expr.args[2], exprs...)
	expr.args[3].args = fields
	# generate convenience outer constructors if necessary
	# if there are nested structs or padding:
	# build an outer constructor that takes all direct, original fields
	# recursively flatten/splat all nested structs into one big args tuple
	# adding zeros for padded arguments
	# pass big, flat, args tuple to inner constructor
	T = expr.args[2]
	if any(x->!FlatBuffers.isbitstype(Core.eval(__module__, x.args[2])), exprs) ||
		length(fields) > length(exprs)
		exprs2 = map(x->FlatBuffers.isbitstype(Core.eval(__module__, x.args[2])) ? x.args[1] : x, exprs)
		sig = Expr(:call, T, exprs2...)
		body = Expr(:call, T, values...)
		outer = Expr(:function, sig, body)
	else
		outer = :(nothing)
	end
	return esc(quote
		$expr
		$outer
	end)
end

macro DEFAULT(T, kwargs...)
	ifblock = quote end
	if length(kwargs) > 0 && isa(kwargs[1], Expr) && length(kwargs[1].args) > 1
		for kw in kwargs
			push!(ifblock.args, :(if sym == $(QuoteNode(kw.args[1]))
				return $(kw.args[2])
			end))
		end
	end
	esc(quote
		if $T <: Enum
			FlatBuffers.default(::Type{$T}) = FlatBuffers.enumtype($T)($(kwargs[1]))
		else
			function FlatBuffers.default(::Type{$T}, TT, sym)
				$ifblock
				return FlatBuffers.default(TT)
			end
		end
	end)
end

import Parameters

function getdef(typedef::Expr)
	isstructexpr(x) = x isa Expr && x.head == :struct
	i = findfirst(x -> x isa Expr && any(isstructexpr.(x.args)), typedef.args)
	if i == nothing
		throw(ArgumentError("malformed @with_kw expression"))
	end
	wrapper = typedef.args[i]
	i = findfirst(isstructexpr, wrapper.args)
	def = wrapper.args[i]
	return def
end

function getfielddefs(typedef::Expr)
	[a for a in getdef(typedef).args[end].args if a isa Expr && a.head == :(::)]
end

function getconstructor(typedef::Expr)
	cons = getdef(typedef).args[end].args[end-1].args[1]
	typevars = []
	if :head in propertynames(cons) && cons.head == :where
		typevars = cons.args[2:end]
		cons = cons.args[1]
	end
	cons, typevars
end

function getkwtype(defs, name)
	for d in defs
		if :args in propertynames(d) && d.args[1] == name
			return d.args[end]
		end
	end
	return nothing
end

function createdefaultfns(typedef::Expr)
	cons, typevars = getconstructor(typedef)
	T = cons.args[1]
	params = cons.args[2]

	@assert params.head == :parameters

	kwargs = []
	defs = getfielddefs(typedef)
	kwdict = Dict{Any, Any}()
	for p in params.args
		name = p.args[1]
		t = getkwtype(defs, name)
		if t == nothing
			continue
		end
		value = p.args[end]
		ifblock = get(kwdict, t, quote end)
		push!(ifblock.args, :(if sym == $(QuoteNode(name))
			return convert($t, $value)
		end))
		kwdict[t] = ifblock
	end

	[:(function FlatBuffers.default(::Type{$T}, ::Type{$TT}, sym) where {$(typevars...)}
		$(kwdict[TT])
		return FlatBuffers.default($TT)
	end) for TT in keys(kwdict)]
end

macro with_kw(typedef)
	body = Parameters.with_kw(typedef, __module__, true)
	defaults = createdefaultfns(body)
	defaultsblock = Expr(:block, body, defaults...)
	esc(defaultsblock)
end

#TODO:
# handle id?
# nested_flatbuffer


text:
module Visualizations

using Plots

end # module


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# Stream.jl: functions to control stream interface

"""
	buffermode(::Camera) -> String

	The current stream buffer mode of the camera.
"""
function buffermode(cam::Camera)
	get(SpinEnumNode(cam::Camera, "StreamBufferHandlingMode", CameraTLStreamNodeMap()))
end

"""
	buffermode!(::Camera, mode::String) -> String

	Set the stream buffer mode of the camera, options include:
	 - OldestFirst
	 - OldestFirstOverwrite
	 - NewestFirst
	 - NewestFirstOverwrite
	 - NewestOnly

	See https://www.ptgrey.com/tan/11174 for more details.
"""
function buffermode!(cam::Camera, mode::String)
	set!(SpinEnumNode(cam, "StreamBufferHandlingMode", CameraTLStreamNodeMap()), mode)
end

"""
	buffercount(::Camera) -> (Int, String)

	Return the buffer count mode and specified number of buffers.
"""
function buffercount(cam::Camera)
	count = get(SpinIntegerNode(cam, "StreamBufferCountResult", CameraTLStreamNodeMap()))
	mode = get(SpinEnumNode(cam, "StreamBufferCountMode", CameraTLStreamNodeMap()))
	return Int(count), mode
end

"""
	buffercount!(::Camera) -> (Int, String)

	Set buffer count mode to automatic. Returns the buffer count and mode.
"""
function buffercount!(cam::Camera)
	set!(SpinEnumNode(cam, "StreamBufferCountMode", CameraTLStreamNodeMap()), "Auto")
	buffercount(cam)
end

"""
	buffercount!(::Camera, count) -> (Int, String)

	Set buffer count mode to manual, and specify the number of buffers. Returns the
	buffer count and mode.
"""
function buffercount!(cam::Camera, count)
	set!(SpinEnumNode(cam, "StreamBufferCountMode", CameraTLStreamNodeMap()), "Manual")
	set!(SpinIntegerNode(cam, "StreamBufferCountManual", CameraTLStreamNodeMap()), count)
	buffercount(cam)
end

"""
	bufferunderrun(::Camera) -> Int

	The number of buffer underruns on the camera stream.
"""
function bufferunderrun(cam::Camera)
	get(SpinIntegerNode(cam, "StreamBufferUnderrunCount", CameraTLStreamNodeMap()))
end

"""
	bufferfailed(::Camera) -> Int

	The number of failed (invalid) buffers on the camera stream.
"""
function bufferfailed(cam::Camera)
	get(SpinIntegerNode(cam, "StreamFailedBufferCount", CameraTLStreamNodeMap()))
end


text:
function wat2wasm(wat::AbstractString)
    out = WasmByteVec()
    @wt_check wasmtime_wat2wasm(
        wat,
        length(wat),
        out,
    )
    out
end

macro wat_str(wat::String)
    :(wat2wasm($wat))
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ThermalStandard <: ThermalGen
        name::String
        available::Bool
        bus::Bus
        activepower::Float64
        reactivepower::Float64
        tech::Union{Nothing, TechThermal}
        op_cost::ThreePartCost
        services::Vector{Service}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Data Structure for thermal generation technologies.

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `activepower::Float64`, validation range: tech.activepowerlimits, action if invalid: warn
- `reactivepower::Float64`, validation range: tech.reactivepowerlimits, action if invalid: warn
- `tech::Union{Nothing, TechThermal}`
- `op_cost::ThreePartCost`
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct ThermalStandard <: ThermalGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    tech::Union{Nothing, TechThermal}
    op_cost::ThreePartCost
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function ThermalStandard(name, available, bus, activepower, reactivepower, tech, op_cost, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    ThermalStandard(name, available, bus, activepower, reactivepower, tech, op_cost, services, ext, forecasts, InfrastructureSystemsInternal(), )
end

function ThermalStandard(; name, available, bus, activepower, reactivepower, tech, op_cost, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    ThermalStandard(name, available, bus, activepower, reactivepower, tech, op_cost, services, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function ThermalStandard(::Nothing)
    ThermalStandard(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        tech=TechThermal(nothing),
        op_cost=ThreePartCost(nothing),
        services=Device[],
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::ThermalStandard) = value.name
"""Get ThermalStandard available."""
get_available(value::ThermalStandard) = value.available
"""Get ThermalStandard bus."""
get_bus(value::ThermalStandard) = value.bus
"""Get ThermalStandard activepower."""
get_activepower(value::ThermalStandard) = value.activepower
"""Get ThermalStandard reactivepower."""
get_reactivepower(value::ThermalStandard) = value.reactivepower
"""Get ThermalStandard tech."""
get_tech(value::ThermalStandard) = value.tech
"""Get ThermalStandard op_cost."""
get_op_cost(value::ThermalStandard) = value.op_cost
"""Get ThermalStandard services."""
get_services(value::ThermalStandard) = value.services
"""Get ThermalStandard ext."""
get_ext(value::ThermalStandard) = value.ext

InfrastructureSystems.get_forecasts(value::ThermalStandard) = value.forecasts
"""Get ThermalStandard internal."""
get_internal(value::ThermalStandard) = value.internal


text:
## foldl & foldr

# folds -- reduce.jl
@test foldl(-,[1:5]) == -13
@test foldl(-,10,[1:5]) == foldl(-,[10,1:5])

@test foldr(-,[1:5]) == 3
@test foldr(-,10,[1:5]) == foldr(-,[1:5,10])

# reduce -- reduce.jl
@test reduce((x,y)->"($x+$y)", [9:11]) == "((9+10)+11)"
@test reduce(max, [8 6 7 5 3 0 9]) == 9
@test reduce(+, 1000, [1:5]) == (1000 + 1 + 2 + 3 + 4 + 5)

# mapreduce -- reduce.jl
@test mapreduce(-, +, [-10 -9 -3]) == ((10 + 9) + 3)
@test mapreduce((x)->x[1:3], (x,y)->"($x+$y)", ["abcd", "efgh", "01234"]) == "((abc+efg)+012)"

z = zeros(2,2,2,2)
for i=1:16
    z[i] = i
end

@test sum(z) == sum(z,(1,2,3,4))[1] == 136

# check variants of summation for type-stability and other issues (#6069)
sum2(itr) = invoke(sum, (Any,), itr)
plus(x,y) = x + y
sum3(A) = reduce(plus, A)
sum4(itr) = invoke(reduce, (Function, Any), plus, itr)
sum5(A) = reduce(plus, 0, A)
sum6(itr) = invoke(reduce, (Function, Int, Any), plus, 0, itr)
sum7(A) = mapreduce(x->x, plus, A)
sum8(itr) = invoke(mapreduce, (Function, Function, Any), x->x, plus, itr)
sum9(A) = mapreduce(x->x, plus, 0, A)
sum10(itr) = invoke(mapreduce, (Function, Function, Int, Any), x->x,plus,0,itr)
for f in (sum2, sum5, sum6, sum9, sum10)
    @test sum(z) == f(z)
    @test sum(Int[]) == f(Int[]) == 0
    @test sum(Int[7]) == f(Int[7]) == 7
    @test typeof(f(Int8[])) == typeof(f(Int8[1])) == typeof(f(Int8[1 7]))
end
for f in (sum3, sum4, sum7, sum8)
    @test sum(z) == f(z)
    @test_throws ErrorException f(Int[])
    @test sum(Int[7]) == f(Int[7]) == 7
end
@test typeof(sum(Int8[])) == typeof(sum(Int8[1])) == typeof(sum(Int8[1 7]))

prod2(itr) = invoke(prod, (Any,), itr)
@test prod(Int[]) == prod2(Int[]) == 1
@test prod(Int[7]) == prod2(Int[7]) == 7
@test typeof(prod(Int8[])) == typeof(prod(Int8[1])) == typeof(prod(Int8[1 7])) == typeof(prod2(Int8[])) == typeof(prod2(Int8[1])) == typeof(prod2(Int8[1 7]))

v = cell(2,2,1,1)
v[1,1,1,1] = 28.0
v[1,2,1,1] = 36.0
v[2,1,1,1] = 32.0
v[2,2,1,1] = 40.0

@test isequal(v,sum(z,(3,4)))

@test sum_kbn([1,1e100,1,-1e100]) == 2

z = rand(10^6)
let es = sum_kbn(z), es2 = sum_kbn(z[1:10^5])
    @test (es - sum(z)) < es * 1e-13
    cs = cumsum(z)
    @test (es - cs[end]) < es * 1e-13
    @test (es2 - cs[10^5]) < es2 * 1e-13
end

@test_throws ErrorException sum(sin, Int[])
@test Base.sumabs(Float64[]) === 0.0
@test Base.sumabs2(Float64[]) === 0.0

@test sum(sin, [1]) == sin(1)
@test Base.sumabs([int8(-2)]) === 2
@test Base.sumabs2([int8(-2)]) === 4

x = -2:3
@test sum(sin, x) == sum(sin(x))
@test Base.sumabs(x) === 9
@test Base.sumabs2(x) === 19

@test_approx_eq sum(sin, z) sum(sin(z))
@test_approx_eq Base.sumabs(z) sum(abs(z))
@test_approx_eq Base.sumabs2(z) sum(abs2(z))

@test maximum(5) == 5
@test minimum(5) == 5

@test maximum([4, 3, 5, 2]) == 5
@test minimum([4, 3, 5, 2]) == 2
@test extrema([4, 3, 5, 2]) == (2, 5)

@test isnan(maximum([NaN]))
@test isnan(minimum([NaN]))
@test isequal(extrema([NaN]), (NaN, NaN))

@test maximum([4., 3., NaN, 5., 2.]) == 5.
@test minimum([4., 3., NaN, 5., 2.]) == 2.
@test extrema([4., 3., NaN, 5., 2.]) == (2., 5.)

@test extrema(1:5) == (1,5)

@test Base.maxabs(-2) == 2
@test Base.minabs(-2) == 2
@test Base.maxabs([1, -2, 3, -4]) == 4
@test Base.minabs([-1, 2, -3, 4]) == 1

@test maximum(abs2, 3:7) == 49
@test minimum(abs2, 3:7) == 9

@test any([true false; false false], 2) == [true false]'
@test any([true false; false false], 1) == [true false]

@test all([true true; false true], 2) == [true false]'
@test all([true false; false true], 1) == [false false]

## cumsum, cummin, cummax

@test isequal(cummin([1, 2, 5, -1, 3, -2]), [1, 1, 1, -1, -1, -2])
@test isequal(cummax([1, 2, 5, -1, 3, -2]), [1, 2, 5, 5, 5, 5])

@test isequal(cummax([1 0; 0 1], 1), [1 0; 1 1])
@test isequal(cummax([1 0; 0 1], 2), [1 1; 0 1])
@test isequal(cummin([1 0; 0 1], 1), [1 0; 0 0])
@test isequal(cummin([1 0; 0 1], 2), [1 0; 0 0])

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
@testset "Active Information" begin
    @test_throws MethodError ActiveInfo(0.5; k=1.0)
    for b in -2:1
        @test_throws ArgumentError ActiveInfo(b; k=2)
    end
    for k in -2:0
        @test_throws ArgumentError ActiveInfo(2; k)
    end

    let ai = ActiveInfo(2; k=1)
        @test ai.joint == zeros(Int, 2, 2)
        @test ai.future == zeros(Int, 2)
        @test ai.history == zeros(Int, 2)
        @test ai.bs == (2,)
        @test ai.k == 1
        @test ai.N == zero(Int)
    end

    let ai = ActiveInfo(2; k=3)
        @test ai.joint == zeros(Int, 2, 8)
        @test ai.future == zeros(Int, 2)
        @test ai.history == zeros(Int, 8)
        @test ai.bs == (2,)
        @test ai.k == 3
        @test ai.N == zero(Int)
    end

    let ai = ActiveInfo(2; k=2)
        observe!(ai, [2,1,2,1])
        @test ai.joint == [0 1 0 0; 0 0 1 0]
        @test ai.future == [1,1]
        @test ai.history == [0,1,1,0]
        @test ai.bs == (2,)
        @test ai.k == 2
        @test ai.N == 2
    end

    let ai = ActiveInfo(2; k=2)
        observe!(ai, [1,1,1,1])
        observe!(ai, [1,2,2,2])
        @test ai.joint == [2 0 0 0; 0 1 0 1]
        @test ai.future == [2,2]
        @test ai.history == [2,1,0,1]
        @test ai.N == 4
    end

    @test_throws ArgumentError ActiveInfo(Int[]; k=2)
    @test_throws ArgumentError ActiveInfo([1]; k=2)
    @test_throws ArgumentError ActiveInfo([1,2]; k=2)
    @test_throws BoundsError ActiveInfo([0,1,1]; k=2)

    @test_throws ArgumentError observe!(ActiveInfo(2; k=2), [1,2])
    @test_throws BoundsError observe!(ActiveInfo(2; k=2), [1,2,3])
    @test_throws BoundsError observe!(ActiveInfo(2; k=2), [0,0,1,1])
    @test_throws BoundsError observe!(ActiveInfo(2; k=2), [2,3,2,3])

    let as = [2,2,1,1,2,1,1,2]
        @test activeinfo(as; k=2) ≈ 0.918296 atol=1e-6
    end

    let as = [1,2,2,2,2,2,2,2,2]
        @test activeinfo(as; k=2) ≈ 0.0 atol=1e-6
    end

    let as = [1,1,2,2,2,2,1,1,1]
        @test activeinfo(as; k=2) ≈ 0.305958 atol=1e-6
    end

    let as = [2,1,1,1,1,1,1,2,2]
        @test activeinfo(as; k=2) ≈ 0.347458 atol=1e-6
    end

    let as = [1,1,1,1,1,2,2,1,1]
        @test activeinfo(as; k=2) ≈ 0.399533 atol=1e-6
    end

    let as = [1,1,1,1,2,2,1,1,1]
        @test activeinfo(as; k=2) ≈ 0.399533 atol=1e-6
    end

    let as = [2,2,2,1,1,1,1,2,2]
        @test activeinfo(as; k=2) ≈ 0.305958 atol=1e-6
    end

    let as = [1,1,1,2,2,2,2,1,1]
        @test activeinfo(as; k=2) ≈ 0.305958 atol=1e-6
    end

    let as = [1,1,1,1,1,1,2,2,1]
        @test activeinfo(as; k=2) ≈ 0.347458 atol=1e-6
    end

    let as = [4,4,4,3,2,1,1,1,2]
        @test activeinfo(as; k=2) ≈ 1.270942 atol=1e-6
    end

    let as = [3,3,4,4,4,4,3,2,1]
        @test activeinfo(as; k=2) ≈ 1.270942 atol=1e-6
    end

    let as = [3,3,3,3,3,3,2,2,2]
        @test activeinfo(as; k=2) ≈ 0.469566 atol=1e-6
    end

    let as = [2 2 2 1 2 1 1 1 2;
              2 2 2 2 1 1 1 1 1]
        @test activeinfo(as; k=2) ≈ 1.270942 atol=1e-6
    end

    let as = [1 1 2 2 2 2 1 2 1;
              2 2 2 2 2 2 2 1 1]
        @test activeinfo(as; k=2) ≈ 1.270942 atol=1e-6
    end

    let as = [1 1 1 1 1 1 2 2 2;
              2 2 2 2 2 2 1 1 1]
        @test activeinfo(as; k=2) ≈ 0.469566 atol=1e-6
    end

    let as = zeros(Int, 2, 9, 2)
        as[:,:,1] = [1 1 2 2 2 2 1 2 1;
                     2 2 2 2 2 2 2 1 1]
        as[:,:,2] = [1 1 1 1 1 1 2 2 2;
                     2 2 2 2 2 2 1 1 1]

        @test activeinfo(as; k=2) ≈ 0.959190 atol=1e-6
    end

    @testset "AI is a Mutual Information" for k in 1:9
        let xs = rand(1:2, 1000)
            ai = activeinfo(xs; k)
            @test ai ≈ mutualinfo(encodehistories(xs, k; warn=false)[1:end-1], xs[1+k:end]) atol=1e-6
            @test ai ≈ mutualinfo(encodehistories(xs, k; warn=false, recode=true)[1:end-1], xs[1+k:end]) atol=1e-6
        end
        let xs = rand(1:3, 1000)
            ai = activeinfo(xs; k)
            @test ai ≈ mutualinfo(encodehistories(xs, k; warn=false)[1:end-1], xs[1+k:end]) atol=1e-6
            @test ai ≈ mutualinfo(encodehistories(xs, k; warn=false, recode=true)[1:end-1], xs[1+k:end]) atol=1e-6
        end
        let xs = rand(1:4, 1000)
            ai = activeinfo(xs; k)
            @test ai ≈ mutualinfo(encodehistories(xs, k; warn=false)[1:end-1], xs[1+k:end]) atol=1e-6
            @test ai ≈ mutualinfo(encodehistories(xs, k; warn=false, recode=true)[1:end-1], xs[1+k:end]) atol=1e-6
        end
    end
end


text:
using PyPlot, MRIReco, FFTW, LinearAlgebra, BenchmarkTools

FFTW.set_num_threads(1);BLAS.set_num_threads(1)

N = 256
numCoils = 8
I = shepp_logan(N)
I = circularShutterFreq!(I,1)

coilsens = birdcageSensitivity(N, numCoils, 1.5)
I = circularShutterFreq!(I,1)
cmap = 1im*quadraticFieldmap(N,N,125*2pi)

# simulation parameters
params = Dict{Symbol, Any}()
params[:simulation] = "fast"
params[:trajName] = "Spiral"
params[:numProfiles] = 6
params[:numSamplingPerProfile] = div(N*N,16)
params[:windings] = div(N,16)
params[:AQ] = 2.0e-2
params[:senseMaps] = coilsens
params[:verbose] = false
#params[:correctionMap] = cmap #[:,:,1]

# do simulation
acqData = simulation(I, params)

raw = RawAcquisitionData(acqData)
fout = ISMRMRDFile(joinpath(@__DIR__(),"simparallel.h5"))
save(fout, raw)


text:
##
## A cowering cow
##
"""
    function cower()

A cowering cow

# Example

```jldoctest
julia> cowsay("Is it safe to come out yet?", cow=Cowsay.cower)
 _____________________________
< Is it safe to come out yet? >
 -----------------------------
     \\
      \\
        ,__, |    |
        (oo)\\|    |___
        (__)\\|    |   )\\_
             |    |_w |  \\
             |    |  ||   *

             Cower....

```
"""
function cower(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
     $thoughts
      $thoughts
        ,__, |    |
        (oo)\\|    |___
        (__)\\|    |   )\\_
             |    |_w |  \\
             |    |  ||   *

             Cower....
"""
return the_cow
end


text:
using Nextion
using Test

@testset "waveform" begin

    function get_values()
        values = Vector{UInt8}(undef, 4)
        delta_v = UInt8(51)
        v = 0
        v = v + delta_v
        values[1] = v
        v = v + delta_v
        values[2] = v
        v = v + delta_v
        values[3] = v
        v = v + delta_v
        values[4] = v
        values
    end

    nexSerial = NexSerial("/dev/ttyUSB0")
    nexWaveform = NexWaveform(nexSerial, Name("s0"), cid=ComponentID(2))

    send(nexSerial, "cls WHITE")
    send(nexSerial, "page pg_waveform")

    nexWaveform.grid.width = 20
    nexWaveform.grid.height = 20
    nexWaveform.grid.color = Colour.GRAY

    sleep(1)

    #send(nexSerial, "s0.pco0=65535")

    #nb_channels = length(nexWaveform.channels)
    #@test nb_channels == 4
    nb_channels = 4

    channel = nexWaveform.channels[2] # 0, 1, 2, 3
    channel.color = Colour.RED

    #push!(channel, UInt8(123))  # values in 0:255 (add)
    #push!(channel, UInt8[123, 123, 100])  # (addt)

    cid = 2
    values = get_values()
    while true
        for ch in 0:nb_channels - 1
            channel = nexWaveform.channels[ch]
            value = values[ch + 1]
            push!(channel, value)
            values[ch + 1] = values[ch + 1] + 2 * rand(-1:1)
        end
        sleep(0.01)
    end

    close(nexSerial)
    @test true

end


text:
using DelimitedFiles

function readSMAT(filename::AbstractString)
    f = open(filename)
    header = readline(f)
    headerparts = split(header)
    nedges = parse(Int,headerparts[3])
    ei = zeros(Int64,nedges)
    ej = zeros(Int64, nedges)
    ev = zeros(Float64, nedges)
    @inbounds for i = 1:nedges
        curline = readline(f)
        parts = split(curline)
        ei[i] = parse(Int, parts[1])+1
        ej[i] = parse(Int, parts[2])+1
        ev[i] = parse(Float64, parts[3])
    end
    close(f)
    A = sparse(ei, ej, ev,
               parse(Int,headerparts[1]),
               parse(Int,headerparts[2])
               )
    return A
end

mutable struct MatrixNetworkMetadata
    A::SparseMatrixCSC{Int64,Int64}
    labels::Vector{AbstractString}
    xy::Array{Float64,2}
    source::AbstractString
end

function load_matrix_network_all(name::AbstractString)
    A = load_matrix_network(name)
    pathname = joinpath(dirname(dirname(@__FILE__)),"data")

    meta_source = joinpath(pathname,"$(name).source")
    if isfile(meta_source)
        # source = open(readstring, meta_source)
        source = open(s->read(s,String),meta_source)
    else
        source = "(None given)"
    end

    meta_xy = joinpath(pathname,"$(name).xy")
    if isfile(meta_xy)
        xy = readdlm(meta_xy)
    else
        xy = zeros(0,2)
    end

    meta_labels = joinpath(pathname,"$(name).labels")
    if isfile(meta_labels)
        labels = open(readlines, meta_labels)
    else
        labels = map(x -> @sprintf("%i",x), 1:size(A,1))
    end

    return MatrixNetworkMetadata(A,labels,xy,source)
end

function load_matrix_network(name::AbstractString)
    pathname = joinpath(dirname(dirname(@__FILE__)),"data")
    smatfile = joinpath(pathname,"$(name).smat")
    if isfile(smatfile)
        return readSMAT(smatfile)
    else
        error(@sprintf "The example datafile '%s' does not seem to exist where it should\n" name)
    end
end

function load_matrix_network_metadata(name::AbstractString)
    pathname = joinpath(dirname(dirname(@__FILE__)),"data")
    smatfile = joinpath(pathname,"$(name).smat")
    meta_xy = joinpath(pathname,"$(name).xy")
    meta_labels = joinpath(pathname,"$(name).labels")
    if isfile(smatfile)
        if isfile(meta_xy) && isfile(meta_labels)
            xy = readdlm(meta_xy)
            labels = readdlm(meta_labels)
            return (readSMAT(smatfile),xy,labels)
        end
    else
        error(@sprintf "The example datafile '%s' does not seem to exist where it should\n" name)
    end
end

function matrix_network_datasets()
    datasets_location = joinpath(dirname(dirname(@__FILE__)),"data")
    content = readdir(datasets_location)
    cc = map(i->match(r".smat",content[i]),1:length(content))
    smat_files = content[findall(typeof.(cc).!=Nothing)]
    for i = 1:length(smat_files)
        smat_files[i] = smat_files[i][1:end-5]
    end
    return smat_files
end


text:
include("downsample.jl")
include("slice.jl")


text:
using ACE, StaticArrays

𝒓 = EuclideanVectorState("𝒓")
μ = SpeciesState("μ")
X1 = μ ⊗ 𝒓
X = AtomState()

X1.rr
X1.mu

@code_llvm X.rr
@code_llvm X1.rr

length(X)  # wait and see what we want `length` to mean...

using BenchmarkTools

function test_perf(X)
   rr = X.rr + rand(SVector{3, Float64})
   return rr
end

@code_warntype test_perf(X1)
@code_warntype getproperty(X1, :rr)

struct S3{T}
   fields::T
end

Base.getproperty(x::S3, sym::Symbol) =
      sym === :fields ? getfield(x, :fields) : getproperty(x.fields, sym)

X3 = S3(X2)

X2 = (mu = AtomicNumber(0), rr = zero(SVector{3, Float64}))

@code_warntype test_perf(X3)

#----

using StaticArrays

struct S{T}
   fields::T
end

Base.getproperty(x::S, sym::Symbol) =
      sym === :fields ? getfield(x, :fields) : getindex(x.fields, sym)

T = SVector{3, Float64}
X1 = (a = 0, b = zero(T))
X2 = S(X1)

function test(X, T)
   rand(T) + X.b
end

@code_warntype test(X1, T)

@code_warntype test(X2, T)


text:
"""
    bijector(model::Model[, sym2ranges = Val(false)])

Returns a `Stacked <: Bijector` which maps from the support of the posterior to ℝᵈ with `d`
denoting the dimensionality of the latent variables.
"""
function Bijectors.bijector(
    model::DynamicPPL.Model,
    ::Val{sym2ranges} = Val(false),
) where {sym2ranges}
    varinfo = DynamicPPL.VarInfo(model)
    num_params = sum([size(varinfo.metadata[sym].vals, 1)
                      for sym ∈ keys(varinfo.metadata)])

    dists = vcat([varinfo.metadata[sym].dists for sym ∈ keys(varinfo.metadata)]...)

    num_ranges = sum([length(varinfo.metadata[sym].ranges)
                      for sym ∈ keys(varinfo.metadata)])
    ranges = Vector{UnitRange{Int}}(undef, num_ranges)
    idx = 0
    range_idx = 1

    # ranges might be discontinuous => values are vectors of ranges rather than just ranges
    sym_lookup = Dict{Symbol, Vector{UnitRange{Int}}}()
    for sym ∈ keys(varinfo.metadata)
        sym_lookup[sym] = Vector{UnitRange{Int}}()
        for r ∈ varinfo.metadata[sym].ranges
            ranges[range_idx] = idx .+ r
            push!(sym_lookup[sym], ranges[range_idx])
            range_idx += 1
        end

        idx += varinfo.metadata[sym].ranges[end][end]
    end

    bs = Bijectors.bijector.(tuple(dists...))

    if sym2ranges
        return (
            Bijectors.Stacked(bs, ranges),
            (; collect(zip(keys(sym_lookup), values(sym_lookup)))...),
        )
    else
        return Bijectors.Stacked(bs, ranges)
    end
end

"""
    meanfield([rng, ]model::Model)

Creates a mean-field approximation with multivariate normal as underlying distribution.
"""
meanfield(model::DynamicPPL.Model) = meanfield(Random.GLOBAL_RNG, model)
function meanfield(rng::Random.AbstractRNG, model::DynamicPPL.Model)
    # setup
    varinfo = DynamicPPL.VarInfo(model)
    num_params = sum([size(varinfo.metadata[sym].vals, 1)
                      for sym ∈ keys(varinfo.metadata)])

    dists = vcat([varinfo.metadata[sym].dists for sym ∈ keys(varinfo.metadata)]...)

    num_ranges = sum([length(varinfo.metadata[sym].ranges)
                      for sym ∈ keys(varinfo.metadata)])
    ranges = Vector{UnitRange{Int}}(undef, num_ranges)
    idx = 0
    range_idx = 1
    for sym ∈ keys(varinfo.metadata)
        for r ∈ varinfo.metadata[sym].ranges
            ranges[range_idx] = idx .+ r
            range_idx += 1
        end

        # append!(ranges, [idx .+ r for r ∈ varinfo.metadata[sym].ranges])
        idx += varinfo.metadata[sym].ranges[end][end]
    end

    # initial params
    μ = randn(rng, num_params)
    σ = StatsFuns.softplus.(randn(rng, num_params))

    # construct variational posterior
    d = DistributionsAD.TuringDiagMvNormal(μ, σ)
    bs = inv.(Bijectors.bijector.(tuple(dists...)))
    b = Bijectors.Stacked(bs, ranges)

    return Bijectors.transformed(d, b)
end

# Overloading stuff from `AdvancedVI` to specialize for Turing
function AdvancedVI.update(d::DistributionsAD.TuringDiagMvNormal, μ, σ)
    return DistributionsAD.TuringDiagMvNormal(μ, σ)
end
function AdvancedVI.update(td::Bijectors.TransformedDistribution, θ...)
    return Bijectors.transformed(AdvancedVI.update(td.dist, θ...), td.transform)
end
function AdvancedVI.update(
    td::Bijectors.TransformedDistribution{<:DistributionsAD.TuringDiagMvNormal},
    θ::AbstractArray,
)
    μ, ω = θ[1:length(td)], θ[length(td) + 1:end]
    return AdvancedVI.update(td, μ, StatsFuns.softplus.(ω))
end

function AdvancedVI.vi(
    model::DynamicPPL.Model,
    alg::AdvancedVI.ADVI;
    optimizer = AdvancedVI.TruncatedADAGrad(),
)
    q = meanfield(model)
    return AdvancedVI.vi(model, alg, q; optimizer = optimizer)
end

function AdvancedVI.vi(
    model::DynamicPPL.Model,
    alg::AdvancedVI.ADVI,
    q::Bijectors.TransformedDistribution{<:DistributionsAD.TuringDiagMvNormal};
    optimizer = AdvancedVI.TruncatedADAGrad(),
)
    # Initial parameters for mean-field approx
    μ, σs = StatsBase.params(q)
    θ = vcat(μ, StatsFuns.invsoftplus.(σs))

    # Optimize
    AdvancedVI.optimize!(elbo, alg, q, make_logjoint(model), θ; optimizer = optimizer)

    # Return updated `Distribution`
    return AdvancedVI.update(q, θ)
end


text:
using PinGo
using Test


text:
module Day09

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day09PuzzleInfo <: PuzzleInfo end

"""
Day 9 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day09PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 9 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day09PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
    $(FUNCTIONNAME)(input)

Read two-dimensional integer (1:9) heightmap data from input stream `input`.
"""
function readheightmap(input::IO)
    reduce(vcat, [parse(Int, c) for c in line]' for line in eachline(input))
end

"""
    $(FUNCTIONNAME)(index, shape)

Get iterator over valid offset indices at index `index` for 2D array with size `shape`.
"""
function offsetsat(index::CartesianIndex{2}, shape::Tuple{Integer, Integer})
    i_offsets(i) = i == 1 ? (1,) : (i == shape[1] ? (-1,) : (-1, 1))
    j_offsets(j) = j == 1 ? (1,) : (j == shape[2] ? (-1,) : (-1, 1))
    Iterators.flatten((
        (CartesianIndex(i, 0) for i in i_offsets(index[1])),
        (CartesianIndex(0, j) for j in j_offsets(index[2]))
    ))
end

"""
    $(FUNCTIONNAME)(heightmap)

Get bitarray of same size as `heightmap` indicating lowpoints.
"""
function getlowpoints(heightmap::AbstractArray{T, 2}) where T
    islowpoint = trues(size(heightmap))
    for index in CartesianIndices(heightmap)
        for offset in offsetsat(index, size(heightmap))
            islowpoint[index] &= heightmap[index] < heightmap[index + offset]
        end
    end
    islowpoint
end

"""
    $(FUNCTIONNAME)(visited, index, heightmap)

Recursively visit points in basin of `heightmap` from `index` recording in `visited`.
"""
function expandbasin!(
    visited::BitArray{2},
    index::CartesianIndex{2},
    heightmap::AbstractArray{T, 2}
) where T
    visited[index] = true
    for offset in offsetsat(index, size(heightmap))
        if heightmap[index + offset] < 9 && !visited[index + offset]
            expandbasin!(visited, index + offset, heightmap)
        end
    end
end

function solve(info::Part1PuzzleInfo)
    heightmap = open(readheightmap, info.datapath)
    islowpoint = getlowpoints(heightmap)
    sum(heightmap[islowpoint] .+ 1)
end

function solve(info::Part2PuzzleInfo)
    heightmap = open(readheightmap, info.datapath)
    islowpoint = getlowpoints(heightmap)
    basinsizes = Array{Int}(undef, count(islowpoint))
    for (i, lowpointindex) in enumerate(findall(islowpoint))
        visited = falses(size(islowpoint))
        expandbasin!(visited, lowpointindex, heightmap)
        basinsizes[i] = count(visited)
    end
    prod(partialsort(basinsizes, 1:3; rev=true))
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day09.txt"))) == 15
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day09.txt"))) == 1134
end

end # module


text:
using RoadRunner
using Test

ant_str = """
    S1 -> S2; k1*S1;
    k1 = 0.1; S1 = 10; S2 = 2.5
"""

rr = RoadRunner.loada(ant_str)

@testset "compartment" begin
    @test RoadRunner.getNumberOfCompartments(rr) == 1
end

@testset "reaction" begin
    @test RoadRunner.getNumberOfReactions(rr) == 1
    @test RoadRunner.getNumberOfRules(rr) == 0
    @test RoadRunner.getReactionIds(rr) == ["_J0"]
end

@testset "parameters" begin
    @test RoadRunner.getNumberOfGlobalParameters(rr) == 1
    @test RoadRunner.getGlobalParameterIds(rr) == ["k1"]
    @test RoadRunner.getGlobalParameterValues(rr) == [0.1]
end

@testset "species" begin
    @test RoadRunner.getNumberOfFloatingSpecies(rr) == 2
    @test RoadRunner.getFloatingSpeciesIds(rr) == ["S1", "S2"]
    @test RoadRunner.getNumberOfBoundarySpecies(rr) == 0
    @test RoadRunner.getBoundarySpeciesIds(rr) == []
    @test RoadRunner.getFloatingSpeciesInitialConcentrationByIndex(rr, 0) == 10
    @test RoadRunner.getFloatingSpeciesInitialConcentrationByIndex(rr, 1) == 2.5
end


text:
function mark_zebra!(robot, line_side, ortogonal_line_side, num_passes, num_start_passes=0)
    start_side = get_start_side(line_side,ortogonal_line_side)
    #УТВ: start_side - кортеж типа NTuple{2,HorizonSide}, например: (Sud,Ost)

    nun_steps = [get_num_movements!(robot,start_side[i]) for i in 1:2]
    #УТВ: Робот - в стартовом углу

    movements_if_posible!(robot, ortogonal_line_side, num_start_passes) || return
    line_mark!(robot,line_side)
    #ИНВАРИАНТ: линия с Роботом и все предыдущие (по ходу движения в направлении ortogonal_line_side) замаркированы
    while movements_if_posible!(robot,ortogonal_line_side, num_passes) == true
        line_side = inverse(line_side)
        line_mark!(robot,line_side)
    end
    #УТВ: все линии замаркированы

    for s in start_side
        movements!(r,s)
    end
    #УТВ: Робот - снова в стартовом углу

    back_side=inverse(start_side)
    for (i,num) in arange(num_steps)
        movements!(robot,back_side[i], num)
    end
    #УТВ: Робот - в исходном положении
end

function movements_if_posible!(robot, side, max_num_steps)
    for _ in 1:max_num_steps
        isborder(robot,side) && (return false)
        move!(robot,side)
    end
    return true
end

function line_mark!(robot,side)
    putmarker!(robot)
    putmarkers!(robot,side)
end


text:
module KLU

using SparseArrays
using SparseArrays: SparseMatrixCSC
import SparseArrays: nnz

export KLUFactorization, klu

include("../lib/klu_lib.jl")
using .LibKLU

import Base: (\), size, getproperty, setproperty!, propertynames, show

# Convert from 1-based to 0-based indices
function decrement!(A::AbstractArray{T}) where T<:Integer
    for i in eachindex(A); A[i] -= oneunit(T) end
    A
end
decrement(A::AbstractArray{<:Integer}) = decrement!(copy(A))

# Convert from 0-based to 1-based indices
function increment!(A::AbstractArray{T}) where T<:Integer
    for i in eachindex(A); A[i] += oneunit(T) end
    A
end
increment(A::AbstractArray{<:Integer}) = increment!(copy(A))

using .LibKLU:
    SuiteSparse_long,
    klu_l_common,
    klu_common,
    klu_l_defaults,
    klu_defaults,
    klu_l_symbolic,
    klu_symbolic,
    klu_l_free_symbolic,
    klu_free_symbolic,
    klu_l_numeric,
    klu_numeric,
    klu_l_free_numeric,
    klu_zl_free_numeric,
    klu_free_numeric,
    klu_z_free_numeric,
    KLU_OK,
    KLU_SINGULAR,
    KLU_OUT_OF_MEMORY,
    KLU_INVALID,
    KLU_TOO_LARGE,
    klu_l_analyze,
    klu_analyze,
    klu_l_factor,
    klu_factor,
    klu_zl_factor,
    klu_z_factor,
    klu_l_solve,
    klu_zl_solve,
    klu_solve,
    klu_z_solve,
    klu_l_tsolve,
    klu_zl_tsolve,
    klu_tsolve,
    klu_z_tsolve,
    klu_l_extract,
    klu_zl_extract,
    klu_extract,
    klu_z_extract,
    klu_l_sort,
    klu_zl_sort,
    klu_sort,
    klu_z_sort,
    klu_refactor,
    klu_z_refactor,
    klu_l_refactor,
    klu_zl_refactor

using LinearAlgebra

const KLUTypes = Union{Float64, ComplexF64}
const KLUValueTypes = (:Float64, :ComplexF64)

if sizeof(SuiteSparse_long) == 4
    const KLUITypes = Int32
    const KLUIndexTypes = (:Int32,)
else
    const KLUITypes = Union{Int32, Int64}
    const KLUIndexTypes = (:Int32, :Int64)
end

function kluerror(status::Integer)
    if status == KLU_OK
        return
    elseif status == KLU_SINGULAR
        throw(LinearAlgebra.SingularException(0))
    elseif status == KLU_OUT_OF_MEMORY
        throw(OutOfMemoryError())
    elseif status == KLU_INVALID
        throw(ArgumentError("Invalid Status"))
    elseif status == KLU_TOO_LARGE
        throw(OverflowError("Integer overflow has occured"))
    else
        throw(ErrorException("Unknown KLU error code: $status"))
    end
end

kluerror(common::Union{klu_l_common, klu_common}) = kluerror(common.status)

"""
Data structure for parameters of and statistics generated by KLU functions.

# Fields
- `tol::Float64`: Partial pivoting tolerance for diagonal preference
- `btf::Int64`: If `btf != 0` use BTF pre-ordering
- `ordering::Int64`: If `ordering == 0` use AMD to permute, if `ordering == 1` use COLAMD,
if `ordering == 3` use the user provided ordering function.
- `scale::Int64`: If `scale == 1` then `A[:,i] ./= sum(abs.(A[:,i]))`, if `scale == 2` then
`A[:,i] ./= maximum(abs.(A[:,i]))`. If `scale == 0` no scaling is done, and the input is
checked for errors if `scale >= 0`.

See the [KLU User Guide](https://github.com/DrTimothyAldenDavis/SuiteSparse/raw/master/KLU/Doc/KLU_UserGuide.pdf)
for more information.
"""
klu_l_common

"""
Data structure for parameters of and statistics generated by KLU functions.

This is the `Int32` version of [`klu_l_common`](@ref).
"""
klu_common

macro isok(A)
    :(kluerror($(esc(A))))
end

function _klu_name(name, Tv, Ti)
    outname = "klu_" * (Tv === :Float64 ? "" : "z") * (Ti === :Int64 ? "l_" : "_") * name
    return Symbol(replace(outname, "__"=>"_"))
end
function _common(T)
    if T == Int64
        common = klu_l_common()
        ok = klu_l_defaults(Ref(common))
    elseif T == Int32
        common = klu_common()
        ok = klu_defaults(Ref(common))
    else
        throw(ArgumentError("Index type must be Int64 or Int32"))
    end
    if ok == 1
        return common
    else
        throw(ErrorException("Could not initialize common struct."))
    end
end

"""
    KLUFactorization <: Factorization

Matrix factorization type of the KLU factorization of a sparse matrix `A`.
This is the return type of [`klu`](@ref), the corresponding matrix factorization function.

The factors can be obtained from `K::KLUFactorization` via `K.L`, `K.U` and `K.F`
See the [`klu`](@ref) docs for more information.

You typically should not construct this directly, instead use [`klu`](@ref).
"""
mutable struct KLUFactorization{Tv<:KLUTypes, Ti<:KLUITypes} <: LinearAlgebra.Factorization{Tv}
    common::Union{klu_l_common, klu_common}
    _symbolic::Ptr{Cvoid}
    _numeric::Ptr{Cvoid}
    n::Int
    colptr::Vector{Ti}
    rowval::Vector{Ti}
    nzval::Vector{Tv}
    function KLUFactorization(n, colptr, rowval, nzval)
        Ti = eltype(colptr)
        common = _common(Ti)
        obj = new{eltype(nzval), Ti}(common, C_NULL, C_NULL, n, colptr, rowval, nzval)
        function f(klu)
            _free_symbolic(klu)
            _free_numeric(klu)
        end
        return finalizer(f, obj)
    end
end

function _free_symbolic(K::KLUFactorization{Tv, Ti}) where {Ti<:KLUITypes, Tv}
    if Ti == Int64
        klu_l_free_symbolic(Ref(Ptr{klu_l_symbolic}(K._symbolic)), Ref(K.common))
    elseif Ti == Int32
        klu_free_symbolic(Ref(Ptr{klu_symbolic}(K._symbolic)), Ref(K.common))
    end
    K._symbolic = C_NULL
end

for Ti ∈ KLUIndexTypes, Tv ∈ KLUValueTypes
    klufree = _klu_name("free_numeric", Tv, Ti)
    ptr = _klu_name("numeric", :Float64, Ti)
    @eval begin
        function _free_numeric(K::KLUFactorization{$Tv, $Ti})
            $klufree(Ref(Ptr{$ptr}(K._numeric)), Ref(K.common))
            K._numeric = C_NULL
        end
    end
end

function KLUFactorization(A::SparseMatrixCSC{Tv, Ti}) where {Tv<:KLUTypes, Ti<:KLUITypes}
    n = size(A, 1)
    n == size(A, 2) || throw(ArgumentError("KLU only accepts square matrices."))
    # Copying here to match UMFPACK
    # Seems overly defensive, it's probably rare that the user will modify A before the
    # numeric factorization is done.
    return KLUFactorization(n, decrement(A.colptr), decrement(A.rowval), copy(A.nzval))
end

size(K::KLUFactorization) = (K.n, K.n)
function size(K::KLUFactorization, dim::Integer)
    if dim < 1
        throw(ArgumentError("size: dimension $dim out of range"))
    elseif dim == 1 || dim == 2
        return Int(K.n)
    else
        return 1
    end
end

nnz(K::KLUFactorization) = K.lnz + K.unz + K.nzoff

Base.adjoint(K::KLUFactorization) = Adjoint(K)
Base.transpose(K::KLUFactorization) = Transpose(K)

function setproperty!(klu::KLUFactorization, ::Val{:(_symbolic)}, x)
    _free_symbolic(klu)
    setfield!(klu, :(_symbolic), x)
end

function setproperty!(klu::KLUFactorization, ::Val{:(_numeric)}, x)
    _free_numeric(klu)
    setfield!(klu, :(_numeric), x)
end

# Certain sets of inputs must be non-null *together*:
# [Lp, Li, Lx], [Up, Ui, Ux], [Fp, Fi, Fx]
for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
    extract = _klu_name("extract", Tv, Ti)
    sort = _klu_name("sort", Tv, Ti)
    if Tv === :ComplexF64
        call = :($extract(klu._numeric, klu._symbolic, Lp, Li, Lx, Lz, Up, Ui, Ux, Uz, Fp, Fi, Fx, Fz, P, Q, Rs, R, Ref(klu.common)))
    else
        call = :($extract(klu._numeric, klu._symbolic, Lp, Li, Lx, Up, Ui, Ux, Fp, Fi, Fx, P, Q, Rs, R, Ref(klu.common)))
    end
    @eval begin
        function _extract!(
            klu::KLUFactorization{$Tv, $Ti};
            Lp = C_NULL, Li = C_NULL, Up = C_NULL, Ui = C_NULL, Fp = C_NULL, Fi = C_NULL,
            P = C_NULL, Q = C_NULL, R = C_NULL, Lx = C_NULL, Lz = C_NULL, Ux = C_NULL, Uz = C_NULL,
            Fx = C_NULL, Fz = C_NULL, Rs = C_NULL
        )
            $sort(klu._symbolic, klu._numeric, Ref(klu.common))
            ok = $call
            if ok == 1
                return nothing
            else
                kluerror(klu.common)
            end
        end
    end
end

function Base.propertynames(::KLUFactorization, private::Bool=false)
    publicnames = (:lnz, :unz, :nzoff, :L, :U, :F, :q, :p, :Rs, :symbolic, :numeric,)
    privatenames = (:nblocks, :maxblock,  :(_L), :(_U), :(_F))
    if private
        return (publicnames..., privatenames...)
    else
        return publicnames
    end
end

function getproperty(klu::KLUFactorization{Tv, Ti}, s::Symbol) where {Tv<:KLUTypes, Ti<:KLUITypes}
    # Forwards to the numeric struct:
    if s ∈ [:lnz, :unz, :nzoff]
        klu._numeric == C_NULL && throw(ArgumentError("This KLUFactorization has not yet been factored. Try `klu_factor!`."))
        return getproperty(klu.numeric, s)
    end
    if s ∈ [:nblocks, :maxblock]
        klu._symbolic == C_NULL && throw(ArgumentError("This KLUFactorization has not yet been analyzed. Try `klu_analyze!`."))
        return getproperty(klu.symbolic, s)
    end
    if s === :symbolic
        klu._symbolic == C_NULL && throw(ArgumentError("This KLUFactorization has not yet been analyzed. Try `klu_analyze!`."))
        if Ti == Int64
            return unsafe_load(Ptr{klu_l_symbolic}(klu._symbolic))
        else
            return unsafe_load(Ptr{klu_symbolic}(klu._symbolic))
        end
    end
    if s === :numeric
        klu._numeric == C_NULL && throw(ArgumentError("This KLUFactorization has not yet been factored. Try `klu_factor!`."))
        if Ti == Int64
            return unsafe_load(Ptr{klu_l_numeric}(klu._numeric))
        else
            return unsafe_load(Ptr{klu_numeric}(klu._numeric))
        end
    end
    # Non-overloaded parts:
    if s ∉ [:L, :U, :F, :p, :q, :R, :Rs, :(_L), :(_U), :(_F)]
        return getfield(klu, s)
    end
    # Factor parts:
    if s === :(_L)
        lnz = klu.lnz
        Lp = Vector{Ti}(undef, klu.n + 1)
        Li = Vector{Ti}(undef, lnz)
        Lx = Vector{Float64}(undef, lnz)
        Lz = Tv == Float64 ? C_NULL : Vector{Float64}(undef, lnz)
        _extract!(klu; Lp, Li, Lx, Lz)
        return Lp, Li, Lx, Lz
    elseif s === :(_U)
        unz = klu.unz
        Up = Vector{Ti}(undef, klu.n + 1)
        Ui = Vector{Ti}(undef, unz)
        Ux = Vector{Float64}(undef, unz)
        Uz = Tv == Float64 ? C_NULL : Vector{Float64}(undef, unz)
        _extract!(klu; Up, Ui, Ux, Uz)
        return Up, Ui, Ux, Uz
    elseif s === :(_F)
        fnz = klu.nzoff
        # We often don't have an F, so create the right vectors for an empty SparseMatrixCSC
        if fnz == 0
            Fp = zeros(Ti, klu.n + 1)
            Fi = Vector{Ti}()
            Fx = Vector{Float64}()
            Fz = Tv == Float64 ? C_NULL : Vector{Float64}()
        else
            Fp = Vector{Ti}(undef, klu.n + 1)
            Fi = Vector{Ti}(undef, fnz)
            Fx = Vector{Float64}(undef, fnz)
            Fz = Tv == Float64 ? C_NULL : Vector{Float64}(undef, fnz)
            _extract!(klu; Fp, Fi, Fx, Fz)
            # F is *not* sorted on output, so we'll have to do it here:
            for i ∈ 1:length(Fp) - 1
                Fp[end] -= 1 # deal with last index being last + 1
                # find each segment
                first = Fp[i] + 1
                last = Fp[i+1] + 1
                # sort each column of rowval, nzval, and Fz for complex numbers if necessary
                #by the ascending permutation of rowval.
                Fiview = view(Fi, first:last)
                Fxview = view(Fx, first:last)
                P = sortperm(Fiview)
                Fiview .= Fiview[P]
                Fxview .= Fxview[P]
                if length(Fz) == length(Fx)
                    Fzview = view(Fz, first:last)
                    Fzview .= Fzview[P]
                end
                # reset the last entry
                Fp[end] += 1
            end
        end
        return Fp, Fi, Fx, Fz
    end
    if s ∈ [:q, :p, :R, :Rs]
        if s === :Rs
            out = Vector{Float64}(undef, klu.n)
        elseif s === :R
            out = Vector{Ti}(undef, klu.nblocks + 1)
        else
            out = Vector{Ti}(undef, klu.n)
        end
        # This tuple construction feels hacky, there's a better way I'm sure.
        s === :q && (s = :Q)
        s === :p && (s = :P)
        _extract!(klu; NamedTuple{(s,)}((out,))...)
        if s ∈ [:Q, :P, :R]
            increment!(out)
        end
        return out
    end
    if s ∈ [:L, :U, :F]
        if s === :L
            p, i, x, z = klu._L
        elseif s === :U
            p, i, x, z = klu._U
        elseif s === :F
            p, i, x, z = klu._F
        end
        if Tv == Float64
            return SparseMatrixCSC(klu.n, klu.n, increment!(p), increment!(i), x)
        else
            return SparseMatrixCSC(klu.n, klu.n, increment!(p), increment!(i), Complex.(x, z))
        end
    end
end

function LinearAlgebra.issuccess(K::KLUFactorization)
    return K.common.status == KLU_OK && K._numeric != C_NULL
end
function show(io::IO, mime::MIME{Symbol("text/plain")}, K::KLUFactorization)
    if issuccess(K)
        summary(io, K); println(io)
        println(io, "L factor:")
        show(io, mime, K.L)
        println(io, "\nU factor:")
        show(io, mime, K.U)
        F = K.F
        if F !== nothing
            println(io, "\nF factor:")
            show(io, mime, K.F)
        end
    else
        throw(ArgumentError("Failed factorization of type $(typeof(K)). Try `klu_factor!(K)`."))
    end
end

function klu_analyze!(K::KLUFactorization{Tv, Ti}) where {Tv, Ti<:KLUITypes}
    if K._symbolic != C_NULL return K end
    if Ti == Int64
        sym = klu_l_analyze(K.n, K.colptr, K.rowval, Ref(K.common))
    else
        sym = klu_analyze(K.n, K.colptr, K.rowval, Ref(K.common))
    end
    if sym == C_NULL
        kluerror(K.common)
    else
        K._symbolic = sym
    end
    return K
end

# User provided permutation vectors:
function klu_analyze!(K::KLUFactorization{Tv, Ti}, P::Vector{Ti}, Q::Vector{Ti}) where {Tv, Ti<:KLUITypes}
    if K._symbolic != C_NULL return K end
    if Ti == Int64
        sym = klu_l_analyze_given(K.n, K.colptr, K.rowval, P, Q, Ref(K.common))
    else
        sym = klu_analyze_given(K.n, K.colptr, K.rowval, P, Q, Ref(K.common))
    end
    if sym == C_NULL
        kluerror(K.common)
    else
        K._symbolic = sym
    end
    return K
end

for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
    factor = _klu_name("factor", Tv, Ti)
    @eval begin
        function klu_factor!(K::KLUFactorization{$Tv, $Ti})
            K._symbolic == C_NULL  && klu_analyze!(K)
            num = $factor(K.colptr, K.rowval, K.nzval, K._symbolic, Ref(K.common))
            if num == C_NULL
                kluerror(K.common)
            else
                K._numeric = num
            end
            return K
        end
    end
end

"""
    klu_factor!(K::KLUFactorization)

Factor `K` into components `K.L`, `K.U`, and `K.F`.
This function will perform both the symbolic and numeric steps of factoriation on an
existing `KLUFactorization` instance.

The `K.common` struct can be used to modify certain options and parameters, see the
[KLU documentation](https://github.com/DrTimothyAldenDavis/SuiteSparse/raw/master/KLU/Doc/KLU_UserGuide.pdf)
or [`klu_common`](@ref) for more information.
"""
klu_factor!

"""
    klu(A::SparseMatrixCSC) -> K::KLUFactorization
    klu(n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) -> K::KLUFactorization

Compute the LU factorization of a sparse matrix `A` using KLU.

For sparse `A` with real or complex element type, the return type of `K` is
`KLUFactorization{Tv, Ti}`, with `Tv` = `Float64` or `ComplexF64`
respectively and `Ti` is an integer type (`Int32` or `Int64`).

The individual components of the factorization `K` can be accessed by indexing:

| Component | Description                                                      |
|:----------|:-----------------------------------------------------------------|
| `L`       | `L` (lower triangular) part of `LU` of the diagonal blocks       |
| `U`       | `U` (upper triangular) part of `LU` of the diagonal blocks       |
| `F`       | `F` (upper triangular) part of `LU + F`, the off-diagonal blocks |
| `p`       | right permutation `Vector`                                       |
| `q`       | left permutation `Vector`                                        |
| `Rs`      | `Vector` of scaling factors                                      |

The relation between `K` and `A` is

`K.L * K.U + K.F  == K.Rs `\\` A[K.p, K.q]`

`K` further supports the following functions:

- `LinearAlgebra.\\`

!!! note
    `klu(A::SparseMatrixCSC)` uses the KLU library that is part of
    SuiteSparse. As this library only supports sparse matrices with [`Float64`](@ref) or
    `ComplexF64` elements, `lu` converts `A` into a copy that is of type
    `SparseMatrixCSC{Float64}` or `SparseMatrixCSC{ComplexF64}` as appropriate.
"""
function klu(n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) where {Ti<:KLUITypes, Tv<:AbstractFloat}
    if Tv != Float64
        nzval = convert(Vector{Float64}, nzval)
    end
    K = KLUFactorization(n, colptr, rowval, nzval)
    return klu_factor!(K)
end

function klu(n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) where {Ti<:KLUITypes, Tv<:Complex}
    if Tv != ComplexF64
        nzval = convert(Vector{ComplexF64}, nzval)
    end
    K = KLUFactorization(n, colptr, rowval, nzval)
    return klu_factor!(K)
end

function klu(A::SparseMatrixCSC{Tv, Ti}) where {Tv<:Union{AbstractFloat, Complex}, Ti<:KLUITypes}
    n = size(A, 1)
    n == size(A, 2) || throw(DimensionMismatch())
    return klu(n, decrement(A.colptr), decrement(A.rowval), A.nzval)
end

for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
    refactor = _klu_name("refactor", Tv, Ti)
    @eval begin
        function klu!(K::KLUFactorization{$Tv, $Ti}, nzval::Vector{$Tv})
            length(nzval) != length(K.nzval)  && throw(DimensionMismatch())
            K.nzval = nzval
            ok = $refactor(K.colptr, K.rowval, K.nzval, K._symbolic, K._numeric, Ref(K.common))
            if ok == 1
                return K
            else
                kluerror(K.common)
            end
        end
    end
end

function klu!(K::KLUFactorization{ComplexF64}, nzval::Vector{U}) where {U<:Complex}
    return klu!(K, convert(Vector{ComplexF64}, nzval))
end

function klu!(K::KLUFactorization{Float64}, nzval::Vector{U}) where {U<:AbstractFloat}
    return klu!(K, convert(Vector{Float64}, nzval))
end

function klu!(K::KLUFactorization{U}, S::SparseMatrixCSC{U}) where {U}
    size(K) == size(S) || throw(ArgumentError("Sizes of K and S must match."))
    return klu!(K, S.nzval)
end
#B is the modified argument here. To match with the math it should be (klu, B). But convention is
# modified first. Thoughts?
for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
    solve = _klu_name("solve", Tv, Ti)
    @eval begin
        function solve!(klu::KLUFactorization{$Tv, $Ti}, B::StridedVecOrMat{$Tv})
            stride(B, 1) == 1 || throw(ArgumentError("B must have unit strides"))
            klu._numeric == C_NULL && klu_factor!(klu)
            size(B, 1) == size(klu, 1) || throw(DimensionMismatch())
            isok = $solve(klu._symbolic, klu._numeric, size(B, 1), size(B, 2), B, Ref(klu.common))
            isok == 0 && kluerror(klu.common)
            return B
        end
    end
end

for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
    tsolve = _klu_name("tsolve", Tv, Ti)
    if Tv === :ComplexF64
        call = :($tsolve(klu._symbolic, klu._numeric, size(B, 1), size(B, 2), B, conj, Ref(klu.common)))
    else
        call = :($tsolve(klu._symbolic, klu._numeric, size(B, 1), size(B, 2), B, Ref(klu.common)))
    end
    @eval begin
        function solve!(klu::Adjoint{$Tv, KLUFactorization{$Tv, $Ti}}, B::StridedVecOrMat{$Tv})
            conj = 1
            klu = parent(klu)
            stride(B, 1) == 1 || throw(ArgumentError("B must have unit strides"))
            klu._numeric == C_NULL && klu_factor!(klu)
            size(B, 1) == size(klu, 1) || throw(DimensionMismatch())
            isok = $call
            isok == 0 && kluerror(klu.common)
            return B
        end
        function solve!(klu::Transpose{$Tv, KLUFactorization{$Tv, $Ti}}, B::StridedVecOrMat{$Tv})
            conj = 0
            klu = parent(klu)
            stride(B, 1) == 1 || throw(ArgumentError("B must have unit strides"))
            klu._numeric == C_NULL && klu_factor!(klu)
            size(B, 1) == size(klu, 1) || throw(DimensionMismatch())
            isok = $call
            isok == 0 && kluerror(klu.common)
            return B
        end
    end
end

function solve(klu, B)
    X = copy(B)
    return solve!(klu, X)
end
LinearAlgebra.ldiv!(klu::KLUFactorization{Tv}, B::StridedVecOrMat{Tv}) where {Tv<:KLUTypes} =
    solve!(klu, B)
LinearAlgebra.ldiv!(klu::LinearAlgebra.AdjOrTrans{Tv, KLUFactorization{Tv, Ti}}, B::StridedVecOrMat{Tv}) where {Tv, Ti} =
    solve!(klu, B)
function LinearAlgebra.ldiv!(klu::KLUFactorization{<:AbstractFloat}, B::StridedVecOrMat{<:Complex})
    imagX = solve(klu, imag(B))
    realX = solve(klu, real(B))
    map!(complex, B, realX, imagX)
end

function LinearAlgebra.ldiv!(klu::LinearAlgebra.AdjOrTrans{Tv, KLUFactorization{Tv, Ti}}, B::StridedVecOrMat{<:Complex}) where {Tv<:AbstractFloat, Ti}
    imagX = solve(klu, imag(B))
    realX = solve(klu, real(B))
    map!(complex, B, realX, imagX)
end

end


text:
# KDE kernels
# unordered categorical default kernel: aitchison-aitken
function aitchison_aitken(bandwidth::Real, observations::Vector, x::Real)
    num_levels = length(unique(observations))
    kernel_value = ones(length(observations)) * bandwidth / (num_levels - 1)
    idx = observations .== x
    kernel_value[idx] .= (idx * (1 - bandwidth))[idx]
    return kernel_value
end

# ordered categorical default kernel: wang-ryzin
function wang_ryzin(bandwidth::Real, observations::Vector, x::Real)
    kernel_value = 0.5 * (1 - bandwidth) * (bandwidth .^ abs.(observations .- x))
    idx = observations .== x
    kernel_value[idx] = (idx * (1-bandwidth))[idx]
    return kernel_value
end

# continuous default kernel: gaussian
function gaussian(bandwidth::Real, observations::Vector, x::Real)
    (1 / sqrt(2*π)) * exp.(-(observations.-x).^2 / (bandwidth^2*2))
end


text:
#Written By Johnathan Bizzano
using Combinatorics
using BenchmarkTools

using ..Utils
using ..HSMath
import ..Fun

export hditer, uhditer, flatten_uhditer, flatten_hditer, total_count,
      sum_comb, unflatten_uhditer, unflatten_hditer, partition, partition_test, partition, HIterator

"iteration_lambda(index_vector::AbstractArray{T}, total_iteration::Integer)::Bool. Return true when it should stop"
const HIterator{T} = Fun{Bool, Tuple{AbstractArray{T}, Int}}

""" High Dimensional Iteration
      Performs an iteration over length(param_1) dimensions, with the value at param_1 being number of times in that dimension
      Leftward direction means that the direction of iteration is going left
            For Example [1, 1, 1] -> [1, 1, 2] if forward direction
      Forward direction means incrementing positively
            For Example [1, 1, 2] -> [1, 1, 1] if left direction
"""
function hditer(length_vector::AbstractArray{T}, iteration_lambda::HIterator{T}; start_vec=nothing, leftward=true, forward=true) where T
      start_vector::Vector{T} = forward ? (@nullc(start_vec, ones(T, length(length_vector)))) : (@nullc(start_vec, copy(length_vector)))
      index_vector::Vector{T} = copy(start_vector)
      current_index::Int64 = leftward ? 1 : length(length_vector)
      trig_idx::Int64 = leftward ? length(length_vector) : 1
      rldirection::Int64 = leftward ? -1 : 1
      fbdirction::Int64 = forward ? 1 : -1
      total_iteration::Int64 = forward ? 1 : ∏(length_vector)

      @label reset
      for i in (forward ? (index_vector[current_index]:length_vector[current_index]) : 1:index_vector[current_index])
            if current_index != trig_idx
                  current_index -= rldirection
                  @goto reset
            else
                  iteration_lambda(index_vector, total_iteration) && return true
                  index_vector[current_index] += fbdirction
                  total_iteration += fbdirction
            end
      end
      if (leftward && current_index != 1) || (!leftward && (current_index < length(length_vector)))
            index_vector[current_index] = start_vector[current_index]
            current_index += rldirection
            index_vector[current_index] += fbdirction
            @goto reset
      end

      return false
end

""" Uniform High Dimensional Iteration
    Performs an iteration over param_1 dimensions, param_2 times
    Return true when it should end
"""
uhditer(num_dimensions, num_times, iteration_lambda::HIterator{T}; start_vec=nothing, left_dir=true, forward_dir=true) where T = hditer(fill(num_times, num_dimensions), iteration_lambda, start_vec = start_vec, leftward = left_dir, forward = forward_dir)

""" Inverse function of hyper index flattener
"""
function unflatten_uhditer(base, flattenedID::T) where T <: Integer
      indexes = T[]
      id = floor(flattenedID / base)
      new_id = 0
      while (new_id = id / base) > 0
            push!(indexes, id % base + 1)
            id = floor(new_id)
      end
      indexes
end

""" Inverse function of hyper index flattener
    Does not work atm
"""
function unflatten_hditer(dim_count, flattenedID::T; prime_list=HSMath.gen_n_primes(dim_count + 1)) where T <: Integer
      indexes = zeros(T, dim_count)
      id = floor(flattenedID / prime_list[end])
      new_id = 0
      index = 0
      while index < dim_count && (new_id = id / prime_list[end - index]) > 0
            indexes[index + 1] = id % prime_list[end - index] + 1
            id = floor(new_id)
            index += 1
      end
      indexes
end

function flatten_hditer(index_vector::AbstractArray{T}; prime_list=HSMath.gen_n_primes(length(index_vector) + 1)) where T
      sum = 0
      for i in 1:length(index_vector)
            sum += (index_vector[i] - 1) * prime_list[i + 1] ^ i
      end
      sum
end

function flatten_uhditer(base, index_vector::AbstractArray{T}) where T
      sum = 0
      for i in 1:length(index_vector)
            sum += (index_vector[i] - 1) * base ^ i
      end
      sum
end

total_count(base, count) = (base ^ count)

"""   Function that iterates a n length vector such that the combination sum is always equal to k and all elements are the natural numbers
    . Invokes lambda function on every iteration
      Returns if it was stopped or not
"""
function partition(n, iteration_lambda::HIterator{T})::Bool where T
      max_vector = zeros(T, n)
      sum_vector = zeros(T, n)
      index_vector = zeros(T, n)
      for i in 1:n
            partition(n, i, iteration_lambda, max_vector=max_vector, sum_vector=sum_vector, index_vector=index_vector) && return true
            fill(0, max_vector)
            fill(0, sum_vector)
            fill(0, index_vector)
      end
      return false
end

"""   Function that iterates a n length vector such that the combination sum is always equal to k and all elements are the natural numbers
    . Invokes lambda function on every iteration
      Returns if it was stopped or not
"""
function partition(k, n, iteration_lambda::HIterator{T}; max_vector=nothing, sum_vector=nothing, index_vector=nothing)::Bool where T
      if n > 0
            n == 1 && return iteration_lambda([k], 1)
            max_vector = max_vector === nothing ? zeros(T, n) : max_vector
            sum_vector = sum_vector === nothing ? zeros(T, n) : sum_vector
            index_vector = index_vector === nothing ? zeros(T, n) : index_vector
            current_index_index::T = 1
            total_iteration::T = 1
            max_vector[1] = k
            index_vector[1] = max(0, -(k * (n - 1)))

            @label reset
            if index_vector[current_index_index] <= max_vector[current_index_index]
                  if current_index_index != n
                        current_index_index += 1
                        sum_vector[current_index_index] = sum_vector[current_index_index - 1] + index_vector[current_index_index - 1]
                        index_vector[current_index_index] = max(0, -(k * (n - current_index_index - 1) + sum_vector[current_index_index]))
                        max_vector[current_index_index] = k - sum_vector[current_index_index]
                  else
                        iteration_lambda(index_vector, total_iteration) && return true
                        total_iteration += 1
                        index_vector[end] += 1
                  end
                  @goto reset
            end
            if current_index_index != 1
                  current_index_index -= 1
                  index_vector[current_index_index] += 1
                  @goto reset
            end
      end
      return false
end

function partition_test(sum)
      t = time()
      out_array = [1]
      partition(sum, function(x, idx)
          out_array[1] = x[1]
          return false
      end)
      println("", "HyperSphere Completed in " * string(time() - t))
      t = time()
      x = collect(partitions(sum))
      println("", "Paritition Completed in " * string(time() - t))
      return (out_array, x)
end


text:
#
# Dimensions Testing
#
@testset "Dimension" begin

@testset "Basics" begin
    c = Dimension(:C, ["as", "ab", "cd", "sd"])
    o = Dimension(:O, ['2', 'J', 'Q', 'K', 'A'])
    s = Dimension(:S, 'a':2:'z')
    u = Dimension(:U, 0:15)
    l = Dimension(:L, 31)

    @test_throws ArgumentError Dimension(:X, ['a', 'a', 'b'])

    @test spttype(c) == Vector{String}
    @test spttype(o) == Vector{Char}
    @test spttype(s) == StepRange{Char, Int64}
    @test spttype(u) == UnitRange{Int64}
    @test spttype(l) == Base.OneTo{Int64}

    @test eltype(c) == String
    @test eltype(o) == Char
    @test eltype(s) == Char
    @test eltype(u) == Int
    @test eltype(l) == Int
    @test all(values(c) .== ["as", "ab", "cd", "sd"])
    @test all(values(o) .== ['2', 'J', 'Q', 'K', 'A'])
    @test all(values(s) .== 'a':2:'z')
    @test all(values(u) .== 0:15)
    @test all(values(l) .== 1:31)

    @test name(c) == :C
    @test name(o) == :O
    @test name(s) == :S
    @test name(u) == :U
    @test name(l) == :L

    @test length(c) == 4
    @test length(o) == 5
    @test length(s) == 13
    @test length(u) == 16
    @test length(l) == 31

    @test indexin(["ab"], c) == [2]
    @test indexin(["waldo"], c) == [0]
    @test indexin("waldo", c) == 0

    @test indexin(['2'], o) == [1]
    @test indexin('K', o) == 4
    @test indexin(['3'], o) == [0]

    @test indexin(['c'], s) == [2]
    @test indexin(['b'], s) == [0]
    @test indexin('a':'c', s) == [1, 0, 2]

    @test indexin([0, 15, 16], u) == [1, 16, 0]

    @test indexin(0, l) == 0
    @test indexin([16], l)== [16]
    @test indexin([1, 31, 58], l) == [1, 31, 0]

    @test first(c) == "as"
    @test first(o) == '2'
    @test first(s) == 'a'
    @test first(u) == 0
    @test first(l) == 1

    @test last(c) == "sd"
    @test last(o) == 'A'
    @test last(s) == 'y'
    @test last(u) == 15
    @test last(l) == 31

    @test "ab" in c
    @test !("asc" in c)

    @test 'J' in o
    @test !('3' in o)

    @test 'c' in s
    @test !('b' in s)

    @test 13 in u
    @test !(16 in u)

    @test 16 in l
    @test !(0 in l)

    @test all((c .== "cd") .== BitArray([false, false, true, false]))
    @test all((c .== "waldo") .== falses(4))

    @test all((o .== 'Q') .== BitArray([false, false, true, false, false]))
    @test all((o .== '3') .== falses(5))

    @test all((s .== 'b') .== falses(13))

    @test all((u .== 16) .== falses(16))

    @test all((l .== 33) .== falses(31))

    @test all((o .<= 'Q') .== BitArray([true, true, true, false, false]))
    @test all((o .>= 'J') .== BitArray([false, true, true, true, true]))
    @test all((o .> 'A') .== falses(5))
    @test all((o .> 'N') .== falses(5))
    @test all((o .> '0') .== falses(5))
end

@testset "Empty" begin
    d = Dimension(:X, [])
    @test length(d) == 0
    @test name(d) == :X
    @test values(d) == []
    @test isempty(d)
    @test size(d) == (0,)
end

@testset "Equality (& Tuple Constructor)" begin
    X1 = Dimension(:X, 3:5)
    X2 = Dimension(:X, 3:5)
    X3 = Dimension(:X, 3:6)
    s = Dimension(:S, 'a':'y')
    fakeX = Dimension(:X, (3, 4, 5))

    @test X1 == X2
    @test X1 != X3
    @test X2 != X3
    @test X2 != s
    # I don't know how I feel about this
    @test X1 == fakeX
end

@testset "indexin" begin
    @testset "int ranges" begin
        rs = [-10:-1:-15, -15:1:-10, 15:-1:10, 15:-1:10,
            20:-2:10, -10:-2:-20, -21:3:-9]
        ds = map(d -> Dimension(:S, d), rs)

        for d in ds
            @test indexin([-12], values(d)) == indexin([-12], d)
        end
    end

    @testset "char ranges" begin
        rs = ['a':1:'d', 'e':-1:'a', 'a':3:'x']
        ds = map(d -> Dimension(:S, d), rs)

        for d in ds
            @test indexin(['c'], values(d)) == indexin(['c'], d)
            @test indexin(['^'], values(d)) == indexin(['^'], d)
        end
    end

    @testset "more tests" begin
        d = Dimension(:x, 21:3:40)
        xs = [1, 2, 5, 7, 201, 3, 0, -3, 4, 19, 24]
        @test indexin(xs, d) == indexin(xs, values(d))
    end
end

#                                   promotion to OneTo
@testset "OneTo promotion" begin
    d = Dimension(:X, 1:3)
    @test typeof(values(d)) == Base.OneTo{Int64}
end

@testset "update" begin
    x = Dimension(:X, 0:3:15)

    @test_throws ArgumentError update(x, [1, 3])

    i, d = update(x, 0:3:6)
    @test d == Dimension(:X, [0, 3, 6])
    @test i == [1, 2, 3]

    i, d = update(x, [true, true, false, false, true, true])
    @test d == Dimension(:X, [0, 3, 12, 15])
    @test i == [1, 2, 5, 6]

    i, d = update(x, [6])
    @test d == Dimension(:X, [6])
    @test i == [3]

    i, d = update(x, 6)
    @test d == Dimension(:X, [])
    @test i == 3
end

@testset "iteration" begin
    x = Dimension(:X, 0:3:15)
    @test collect(x) == [0, 3, 6, 9, 12, 15]
end

end


text:
using EquivalenceClassesConstructor
using Test

include("../src/Mappings.jl")
include("../src/GraphTools.jl")
include("../src/EquivalenceClasses.jl")

@testset "EquivalenceClassesConstructor.jl" begin
    # Write your tests here.
    @testset "predicates" begin
        include("./Mappings.jl")
    end
    @testset "graph" begin
        include("./GraphTools.jl")
        include("./EquivalenceClasses.jl")
    end
end


text:
function energy(planets)
    eng = 0.0
    # kinetic energy
    for p in planets
        eng += 1/2*p.m*norm2(p.v)
    end
    # potential energy
    for j in 1:nplanets
        pj = planets[j]
        for k in j+1:nplanets
            pk = planets[k]
            eng -= G*pj.m*pk.m/sqdist(pj.r, pk.r)
        end
    end
    eng
end

function barycenter(m,mTot,x,y,z) # Find Barycenter
    #    m : mass of planet
    # mTot : total mass of system
    # x,y,z: position coordinates of planet
    #
    baryX = (m*x)/mTot
    baryY = (m*y)/mTot
    baryZ = (m*z)/mTot
    return baryX,baryY,baryZ;
end

function momentum(m,x,y,z,vx,vy,vz)
    r = [x,y,z];
    vel = [vx,vy,vz];
    H = m*cross(r,vel);
    return H;
end


text:
using Distributions, Random, Test
using EconLH.ProductionFunctionsLH
using EconLH

function pv_test()
	@testset "present value" begin
		R = 1.05;
		T = 3;

		function f(x)
			return 0.7 .* x
		end
		pv = EconLH.present_value(f, R, T);
		@test length(pv) == 1
		@test isa(pv, Float64)
		@test pv ≈ f(1) + f(2) ./ R + f(3) ./ R^2
	end
end

@testset "econLH" begin
	pv_test();
	include("extreme_value_decision_test.jl");
	include("prod_fct_test.jl");
	include("ces_test.jl")
end

# ----------


text:
using MacroTools

"Representable is a functor from a category into a computation "
abstract type StatefulFunctor{C}; end
export StatefulFunctor, @interpretation, @interpret, @functor

macro interpretation(name, variance, category, state_def=Expr(:dummy))
    esc(quote
        mutable struct $name <: StatefulFunctor{$category.Arrow}
        $(state_def.args...)
        end

        (s::$name)(m::$category.Composed, inp...) = s(m.g, s(m.f, inp...)...)
        (s::$name)(m::$category.Product, inp...) = tuple([s(x, inp...) for x in m.factors]...)
        (s::$name)(m::$category.Proj, inp...) = inp[m.m]
        (s::$name)(m::$category.Constant, inp...) = m.val
        (s::$name)(m::$category.Identity, inp) = inp
        (s::$name)(m::$category.Terminal, inp...) = ()
        interp_state_hook(s::$name, m::$category.Arrow, value_expr) = value_expr()
        end)
end

macro interpret(fn)
    @capture(fn, function (s_::functor_)(m_::arrow_, args__) expr__ end) || error("Can't destructure input expression")
    esc(quote
          function ($s::$functor)($m::$arrow, $(args...))
            value_expr = () -> begin $(expr...) end
            interp_state_hook($s, $m, value_expr)
          end
        end)
end

macro functor(sig, obj_map)
    @capture(obj_map, var_ -> body_)
    comp_body =
        if @capture(sig, name_ :: src_ => tgt_)
            :($name(m::$src.Composed) = compose($name(m.g), $name(m.f)))
        elseif @capture(sig, name_ :: src_ <= tgt_)
            :($name(m::$src.Composed) = compose($name(m.f), $name(m.g)))
        else
        end
    esc(quote
        struct $name
        end

        function $name($var::Type)
           $(body.args...)
        end
        $comp_body
        # Functions are (src morphism × input) -> (tgt morphism)
        $name(m::$src.Identity{T}) where {T} = $tgt.Identity{$name(T)}()
        # Cartesian
        $name(m::$src.Product) = $tgt.Product([$name(x) for x in m.factors]...)
        $name(m::$src.Proj{A, B}) where {A,B} = $tgt.Proj{$name(A), $name(B)}(m.m)
        # need to define this
        # $name(m::$src.Constant) = $tgt.Constant(m.val)
        $name(m::$src.Terminal{T}) where {T} = $tgt.Terminal{$name(T)}()
        end)
end


text:
@testset "976.largest-perimeter-triangle.jl" begin
    @test largest_perimeter([2, 1, 2]) == 5
    @test largest_perimeter([1, 2, 1]) == 0
    @test largest_perimeter([3, 2, 3, 4]) == 10
    @test largest_perimeter([3, 6, 2, 3]) == 8
end