text:
module DependencyTrees

export
    DependencyTree, isprojective,
    Treebank,
    TypedDependency, UntypedDependency, CoNLLU,

    typed, untyped,

    ArcEager, ArcEagerReduce, ArcEagerShift,
    ArcStandard, ArcHybrid, ArcSwift, ListBasedNonProjective,

    StaticOracle, static_oracle, static_oracle_prefer_shift,

    DynamicOracle, dynamic_oracle,

    AlwaysExplore, NeverExplore, ExplorationPolicy,

    initconfig, isfinal, possible_transitions,
    stacktoken, buffertoken,

    labeled_accuracy, unlabeled_accuracy

using Random
using LightGraphs

import Base.==
import Base.getindex
import Base.reduce
import Base: iterate, IteratorSize, length

include("errors.jl")
include("dependencies.jl")
include("conllu.jl")
include("trees.jl")

include("treebanks/treebanks.jl")

include("transition_parsing/systems.jl")
include("transition_parsing/oracles.jl")

include("evaluation/accuracy.jl")

end # module


text:
using FixedEffectModels, DataFrames, Statistics, CSV, Test

df = CSV.read(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv"))
df.pState = categorical(df.State)
df.pYear = categorical(df.Year)

@test convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ NDI))[1])[1:5, :] ≈ [ -37.2108  9.72654; -35.5599  9.87628; -32.309   8.82602; -34.2826  9.64653; -35.0526  8.84143] atol = 1e-3
@test convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ NDI + fe(pState)))[1])[1:5, :] ≈ [ -21.3715  -0.642783; -19.6571  -0.540335; -16.3427  -1.63789; -18.2631  -0.856975; -18.9784  -1.70283] atol = 1e-3
@test convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ 1 + fe(pState)))[1])[1:5, :] ≈  [-13.5767   -40.5467; -12.0767   -39.3467; -8.97667  -39.3467; -11.0767   -37.6467; -11.9767   -37.5467] atol = 1e-3
@test convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ 1))[1])[1:5, :] ≈ [ -30.0509  -40.0999; -28.5509  -38.8999; -25.4509  -38.8999; -27.5509  -37.1999; -28.4509  -37.0999] atol = 1e-3
@test mean(convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ NDI), add_mean = true)[1]), dims = 1[1]) ≈ [123.951  68.6999] atol = 1e-3
@test mean(convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ NDI + fe(pState)), add_mean = true)[1]), dims = 1[1]) ≈ [123.951  68.6999] atol = 1e-3
@test mean(convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ 1 + fe(pState)), add_mean = true)[1]), dims = 1[1]) ≈ [123.951  68.6999] atol = 1e-3
@test mean(convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ 1), add_mean = true)[1]), dims = 1[1]) ≈ [123.951  68.6999] atol = 1e-3
@test convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ NDI), weights = :Pop)[1])[1:5, :] ≈[ -37.5296  11.8467; -35.8224  11.9922; -32.5151  10.9377; -34.4416  11.7546; -35.163   10.9459] atol = 1e-3
@test convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ NDI + fe(pState)), weights = :Pop)[1])[1:5, :] ≈ [ -22.2429  -1.2635 ; -20.5296  -1.1515 ; -17.2164  -2.23949; -19.1378  -1.45057; -19.854   -2.28819] atol = 1e-3
@test convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ 1 + fe(pState)), weights = :Pop)[1])[1:5, :] ≈ [ -14.0383   -43.1224; -12.5383   -41.9224; -9.43825  -41.9224; -11.5383   -40.2224; -12.4383   -40.1224] atol = 1e-3
@test convert(Matrix{Float64}, partial_out(df, @formula(Sales + Price ~ 1), weights = :Pop)[1])[1:5, :] ≈ [ -26.3745  -44.9103; -24.8745  -43.7103; -21.7745  -43.7103; -23.8745  -42.0103; -24.7745  -41.9103] atol = 1e-3

ext: jl


text:
module Runner

export @runit

"""

"""
macro runit(command::String)
    # find the filename to include
    # TODO: Remove multiple and trailing whitespaces
    args_list = split(command)
    inc_file = popfirst!(args_list)

    # check if args need to be modified
    _ARGS = [] # to backup existing ARGS
    if length(args_list) > 0
        _ARGS = [ popfirst!(ARGS) for i ∈ 1:length(ARGS) ]
        for a in args_list
            push!(ARGS, string(a))
        end
    end

    # include/run the script
    include(inc_file)

    # replace orginal contents of ARGS
    # delete current contents of ARGS
    # TODO: Is there a better way to do this?
    [ pop!(ARGS) for i ∈ 1:length(ARGS) ]
    for a in _ARGS
         push!(ARGS, a)
    end
end

end # module


text:
######################################################################
# mutual_information.jl
# has things like mutual information and stuff
######################################################################

export mutual_information3
export p_obs
# sum over all possible jammer locations
function p_obs(df::DF, xp::Pose, o)
    prob = 0.0
    for txi = 1:df.n, tyi = 1:df.n
        if df.b[txi,tyi] > 0.0
            tx = (txi-0.5) * df.cell_size
            ty = (tyi-0.5) * df.cell_size
            prob += df.b[txi, tyi] * O(df.sensor, (tx,ty,0.,0.), xp, o)
        end
    end
    return prob
end

function p_obs(pf::PF, xp::Pose, o)
    prob = 0.0
    for i = 1:pf.n
        prob += O(pf.model.sensor, particle(pf,i), xp, o) * weight(pf, i)
    end
    return prob / weight_sum(pf)
end

#function mutual_information(pf::PF, xp::Pose)
#    H_o = 0.0
#    H_o_t = 0.0
#
#    # computing H_o
#    for o in pf.obs_list
#        po = p_obs(pf, xp, o)
#        if po != 0
#            H_o -= po * log(po)
#        end
#    end
#
#    # computing H_o_t
#    ws = weight_sum(pf.b)
#    for i in pf.n
#        weight(pf.b, i) / ws
#        for o in pf.obs_list
#            po =
#            O(pf.model.sensor, pf.b.particles[i], xp, o)
#        end
#    end
#
#    return H_o - H_o_t
#end

# This is just a copy (and modification) of df's version
# The commented out version above is an attempt from scratch
function mutual_information(pf::PF, xp::Pose)
	H_o = 0.0
	H_o_t = 0.0

    for o in pf.obs_list
		po = p_obs(pf, xp, o)
		if po > 0.0
			H_o -= po * log(po)
		end

		# sum over possible jammer locations
        ws = weight_sum(pf)
        for i = 1:pf.n
            pot = O(pf.model.sensor, particle(pf,i), xp, o)
            if pot > 0.0
                H_o_t -= pot * weight(pf,i) * log(pot) / ws
            end
        end
	end
	return H_o - H_o_t
end

# computes mutual information for a specific vehicle location
# xp is a proposed pose
# really need to loop over all possible observations
export mutual_information
# this is the one that currently gets called with fov and greedy policy
function mutual_information(df::DF, xp::Pose)
	H_o = 0.0
	H_o_t = 0.0

    for o in df.obs_list
		po = p_obs(df, xp, o)
		if po > 0.0
			H_o -= po * log(po)
		end

		# sum over possible jammer locations
		for txi = 1:df.n, tyi = 1:df.n
            tx = (txi - 0.5) * df.cell_size
            ty = (tyi - 0.5) * df.cell_size

            if df.b[txi, tyi] > 0.0
                pot = O(df.sensor, (tx,ty,0.,0.), xp, o)
                if pot > 0.0
                    H_o_t -= pot * df.b[txi, tyi] * log(pot)
                end
            end
		end
	end
	return H_o - H_o_t
end

# computes mutual information for all locations
# the provided discrete filter gives the domain size and belief
# I need a Vehicle or Sensor to pass into the observation function
function mutual_information(df::DF)
    mut_info = zeros(df.n, df.n)
    for xi = 1:df.n, yi = 1:df.n
        x = (xi - 0.5) * df.cell_size
        y = (yi - 0.5) * df.cell_size
        # TODO: really, should have some other heading
        xp = (x, y, 0.0)
        mut_info[xi,yi] = mutual_information(df, xp)
    end
    return mut_info
end

export mutual_information_se2
function mutual_information_se2(df::DF)
    mut_info = zeros(df.n, df.n)
    for xi = 1:df.n, yi = 1:df.n, hi = 1:36
        x = (xi - 0.5) * df.cell_size
        y = (yi - 0.5) * df.cell_size
        h = (hi - 0.5) * 10.0
        # TODO: really, should have some other heading
        xp = (x, y, h)
        mut_info[xi,yi] = mutual_information(df, xp)
    end
    return mut_info
end

# TODO: I have no idea if anything below is valid

# this version is valid if the sensor noise is the same regardless
#  of the vehicle and sensor locations
# TODO: is this remotely valid
export mutual_information2
function mutual_information2(df::DF, xp::Pose)
	H_o = 0.0

	for o in df.obs_list
		po = p_obs(df, xp, o)
		if po > 0.0
			H_o -= po * log(po)
		end

	end
	return H_o
end

# for FOV stuff... this is ugly
function mutual_information3(m::SearchDomain, uav::Vehicle, df::DF)
	mut_info = zeros(df.n, df.n, df.n)
	for xv = 1:df.n
		x = (xv - 0.5) * df.cell_size
		for yv = 1:df.n
			y = (yv - 0.5) * df.cell_size
			for zv = 1:df.n
				z = (zv - 0.5) * (360./21.0)
				xp = (x, y, z)
				mut_info[xv,yv,zv] = mutual_information(m, uav, df, xp)
			end
		end
	end
	return mut_info
end

function mutual_information2(m::SearchDomain, uav::Vehicle, df::DF)
	mut_info = zeros(df.n, df.n)
	for xv = 1:df.n
		x = (xv - 0.5) * df.cell_size
		for yv = 1:df.n
			y = (yv - 0.5) * df.cell_size
			# TODO: really, should have some other heading
			xp = (x, y, 0.0)
			mut_info[xv,yv] = mutual_information2(m, uav, df, xp)
		end
	end
	return mut_info
end


text:
@testset "EarthEnv LandCover" begin
    using SimpleSDMDataSources: rastername, rasterurl, rasterpath

    @test rastername(EarthEnv{LandCover}, 2, true) == "landcover_complete_2.tif"
    landcover_path = joinpath(ENV["ECODATASOURCES_PATH"], "EarthEnv/consensus_landcover")
    @test rasterpath(EarthEnv{LandCover}) == landcover_path
    @test rasterpath(EarthEnv{LandCover}, 2, true) == joinpath(landcover_path, "landcover_complete_2.tif")
    @test rasterurl(EarthEnv{LandCover}, 2, true) ==
        "https://data.earthenv.org/consensus_landcover/with_DISCover/consensus_full_class_2.tif"
    download_raster(EarthEnv{LandCover}; layer=2)
    @test isfile(joinpath(landcover_path, "landcover_partial_2.tif"))
end


text: abstract AbstractFixedPointSolver{T} <: NonlinearSolver{T}


text:
"""
naive implementation
"""
function naiveform!(C, A, B)
    row, col  = size(C)
    colA = size(A)[2]
    for i in 1:row
        for j in 1:col
            for k in 1:colA
                C[i, j] = C[i, j] + A[i, k] * B[k, j]
            end
        end
    end
    C
end

"""
inner-product-form
"""
function ipform!(C, A, B)
    row, col  = size(C)
    colA = size(A)[2]
    for i in 1:row
        for j in 1:col
            dataC = 0.0
            for k in 1:colA
                dataC = dataC + A[i, k] * B[k, j]
            end
            C[i, j] = dataC
        end
    end
end

"""
outer-product-form
"""
function opform!(C, A, B)
    row, col  = size(C)
    colA = size(A)[2]
    for k in 1:colA
        for j in 1:col
            dataB = B[k, j]
            for i in 1:row
                C[i, j] = C[i, j] + A[i, k] * dataB
            end
        end
    end
end

"""
middle-product-form
"""
function mpform!(C, A, B)
    row, col  = size(C)
    colA = size(A)[2]
    for j in 1:col
        for k in 1:colA
            dataB = B[k, j]
            for i in 1:row
                C[i, j] = C[i, j] + A[i, k] * dataB
            end
        end
    end
end


text:
"""
Depricated file
"""

using Statistics

###Input Layer
function layerForProp!(
    cLayer::Input,
    X::AbstractArray = Array{Any,1}(undef,0);
    kwargs...,
)
    if length(X) != 0
        cLayer.A = X
        cLayer.inputS = cLayer.outputS = size(X)
    end
    cLayer.forwCount += 1
    # Base.GC.gc()
    return nothing
end

###FCLayer forprop

function layerForProp!(
    cLayer::FCLayer,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    kwargs...,
)
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    cLayer.inputS = cLayer.prevLayer.outputS
    cLayer.Z = cLayer.W * Ai .+ cLayer.B
    actFun = cLayer.actFun
    Z = cLayer.Z
    cLayer.outputS = size(Z)
    cLayer.A = eval(:($actFun($Z)))
    cLayer.forwCount += 1
    # Base.GC.gc()
    return nothing
end #function layerForProp!(cLayer::FCLayer)

###AddLayer forprop

function layerForProp!(cLayer::AddLayer; kwargs...,)
    cLayer.A = similar(cLayer.prevLayer[1].A)
    cLayer.A .= 0
    for prevLayer in cLayer.prevLayer
        cLayer.A .+= prevLayer.A
    end
    cLayer.forwCount += 1
    # Base.GC.gc()
    return nothing
end #function layerForProp!(cLayer::AddLayer)

###Activation forprop

function layerForProp!(
    cLayer::Activation,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    kwargs...,
)
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    actFun = cLayer.actFun
    cLayer.A = eval(:($actFun($Ai)))
    cLayer.inputS = cLayer.outputS = size(Ai)
    cLayer.forwCount += 1
    # Ai = nothing
    # Base.GC.gc()
    return nothing
end #function layerForProp!(cLayer::Activation)

###BatchNorm

function layerForProp!(
    cLayer::BatchNorm,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    kwargs...,
)
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end

    initWB!(cLayer)

    cLayer.μ = mean(Ai, dims = 1:cLayer.dim)
    Ai_μ = Ai .- cLayer.μ
    N = prod(size(Ai)[1:cLayer.dim])
    cLayer.Ai_μ_s = Ai_μ .^ 2
    cLayer.var = sum(cLayer.Ai_μ_s, dims = 1:cLayer.dim) ./ N
    cLayer.Z = Ai_μ ./ sqrt.(cLayer.var .+ cLayer.ϵ)
    cLayer.A = cLayer.W .* cLayer.Z .+ cLayer.B
    cLayer.forwCount += 1
    Ai_μ = nothing
    cLayer.inputS = cLayer.outputS = size(Ai)
    # Base.GC.gc()
    return nothing
end #function layerForProp!(cLayer::BatchNorm)

export layerForProp!


text:
module SolversWrapper

using JavaCall

const AllowedSolution = @jimport org.hipparchus.analysis.solvers.AllowedSolution
const BaseAbstractUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BaseAbstractUnivariateSolver
const BaseSecantSolver = @jimport org.hipparchus.analysis.solvers.BaseSecantSolver
const BaseUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BaseUnivariateSolver
const BisectionSolver = @jimport org.hipparchus.analysis.solvers.BisectionSolver
const BracketedRealFieldUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver
const BracketedRealFieldUnivariateSolver_Interval = @jimport org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver$Interval
const BracketedUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BracketedUnivariateSolver
const BracketedUnivariateSolver_Interval = @jimport org.hipparchus.analysis.solvers.BracketedUnivariateSolver$Interval
const BracketingNthOrderBrentSolver = @jimport org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver
const BrentSolver = @jimport org.hipparchus.analysis.solvers.BrentSolver
const Class = @jimport java.lang.Class
const Enum = @jimport java.lang.Enum
const FieldBracketingNthOrderBrentSolver = @jimport org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver
const IllinoisSolver = @jimport org.hipparchus.analysis.solvers.IllinoisSolver
const JComplex = @jimport org.hipparchus.complex.Complex
const JString = @jimport java.lang.String
const LaguerreSolver = @jimport org.hipparchus.analysis.solvers.LaguerreSolver
const MullerSolver = @jimport org.hipparchus.analysis.solvers.MullerSolver
const MullerSolver2 = @jimport org.hipparchus.analysis.solvers.MullerSolver2
const NewtonRaphsonSolver = @jimport org.hipparchus.analysis.solvers.NewtonRaphsonSolver
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const PegasusSolver = @jimport org.hipparchus.analysis.solvers.PegasusSolver
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const RealFieldUnivariateFunction = @jimport org.hipparchus.analysis.RealFieldUnivariateFunction
const RegulaFalsiSolver = @jimport org.hipparchus.analysis.solvers.RegulaFalsiSolver
const RiddersSolver = @jimport org.hipparchus.analysis.solvers.RiddersSolver
const SecantSolver = @jimport org.hipparchus.analysis.solvers.SecantSolver
const UnivariateDifferentiableFunction = @jimport org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction
const UnivariateFunction = @jimport org.hipparchus.analysis.UnivariateFunction
const UnivariateSolverUtils = @jimport org.hipparchus.analysis.solvers.UnivariateSolverUtils

export AbstractPolynomialSolver
export AbstractUnivariateDifferentiableSolver
export AbstractUnivariateSolver
export AllowedSolution
export BaseAbstractUnivariateSolver
export BaseSecantSolver
export BaseUnivariateSolver
export BisectionSolver
export BracketedRealFieldUnivariateSolver
export BracketedUnivariateSolver
export BracketingNthOrderBrentSolver
export BrentSolver
export FieldBracketingNthOrderBrentSolver
export IllinoisSolver
export LaguerreSolver
export MullerSolver
export MullerSolver2
export NewtonRaphsonSolver
export PegasusSolver
export PolynomialSolver
export RegulaFalsiSolver
export RiddersSolver
export SecantSolver
export UnivariateDifferentiableSolver
export UnivariateSolver
export UnivariateSolverUtils
export bracket
export compare_to
export describe_constable
export do_solve
export equals
export force_side
export get_absolute_accuracy
export get_class
export get_declaring_class
export get_evaluations
export get_function_value_accuracy
export get_max
export get_max_evaluations
export get_maximal_order
export get_min
export get_relative_accuracy
export get_start_value
export hash_code
export is_bracketing
export is_sequence
export midpoint
export name
export notify
export notify_all
export ordinal
export solve
export solve_all_complex
export solve_complex
export solve_interval
export to_string
export value_of
export values
export verify_bracketing
export verify_interval
export verify_sequence
export wait

include("abstract_polynomial_solver.jl")
include("abstract_univariate_differentiable_solver.jl")
include("abstract_univariate_solver.jl")
include("allowed_solution.jl")
include("base_abstract_univariate_solver.jl")
include("base_secant_solver.jl")
include("base_univariate_solver.jl")
include("bisection_solver.jl")
include("bracketed_real_field_univariate_solver.jl")
include("bracketed_univariate_solver.jl")
include("bracketing_nth_order_brent_solver.jl")
include("brent_solver.jl")
include("field_bracketing_nth_order_brent_solver.jl")
include("illinois_solver.jl")
include("laguerre_solver.jl")
include("muller_solver.jl")
include("muller_solver2.jl")
include("newton_raphson_solver.jl")
include("pegasus_solver.jl")
include("polynomial_solver.jl")
include("regula_falsi_solver.jl")
include("ridders_solver.jl")
include("secant_solver.jl")
include("univariate_differentiable_solver.jl")
include("univariate_solver.jl")
include("univariate_solver_utils.jl")

end


text:
const VoxelId = NTuple{3, Int}

"""
    SparseVoxelGrid(points, voxel_size)
    SparseVoxelGrid(points, voxel_size::NTuple{3, AbstractFloat})

Creates a sparse spatial grid by organising 3D points into voxels. `points` can either
be a 3xN matrix or a `PointCloud`. `voxel_size` will either create uniformly sized voxels in each
axis or if it is a 3D tuple the size of each axis can be specified.

### Example

To create a spatial grid with voxel side length of 10 metres for arbitrary points:
```julia
using PointClouds
points = rand(3, 100000) * 20.0
grid = SparseVoxelGrid(points, 10.0)
```

The created grid is an iteratable object which returns a `Voxel` in each iteration.
Each voxel can be accessed directly with a `for` loop or all voxels can be `collect`ed into an array.
Likewise, the returned `Voxel` is an iterable object that returns the point indices:
```julia
# Iterate through each voxel in grid
for voxel in grid
    # Get each point index in voxel
    for idx in voxel
        # Do stuff with points[:,idx]
    end
    # Or, you may want all point indices in a voxel
    all_point_indices = collect(voxel)
end
```
"""
struct SparseVoxelGrid{T <: Real}
    voxel_size::SVector{3, T}
    voxel_info::Dict{VoxelId, UnitRange{Int}}
    point_indices::Vector{Int}
end

function SparseVoxelGrid(points::Vector{T1}, voxel_size::SVector{3, T2}) where {T1 <: AbstractVector, T2 <: Real}
    npoints = length(points)

    # In order to avoid allocating a vector for each voxel, we construct the data structure in a backward-looking order.

    # Assign each point to a voxel id
    voxel_ids = Vector{VoxelId}(undef, npoints)
    for j = 1:npoints
        @inbounds voxel_ids[j] = make_voxel_id(points[j], voxel_size)
    end

    # Count the number of points in each voxel
    group_counts = Dict{VoxelId, Int}()
    for id in voxel_ids
        group_counts[id] = get(group_counts, id, 0) + 1
    end

    # Allocate ranges for the indices of points in each voxel based on the counts
    voxel_info = Dict{VoxelId, UnitRange{Int}}()
    current_index = 1
    # TODO: Using keys(group_counts) is a workaround for inference in julia-0.5 - let's just iterate over the dictionary in the future
    for group_id in keys(group_counts)
        group_size = group_counts[group_id]
        voxel_info[group_id] = UnitRange(current_index, current_index+group_size-1)
        current_index += group_size
    end

    # Place indices for points into the appropriate index range for the associated voxel
    point_indices = Vector{Int}(undef, npoints)
    for j = 1:npoints
        id = voxel_ids[j]
        index_in_group = group_counts[id]
        group_counts[id] = index_in_group - 1
        point_indices[first(voxel_info[id]) + index_in_group-1] = j
    end
    return SparseVoxelGrid(voxel_size, voxel_info, point_indices)
end

function SparseVoxelGrid(points::Matrix{T1}, voxel_size) where T1 <: Real
    ndims, npoints = size(points)
    @assert ndims == 3
    new_data = reshape(reinterpret(SVector{3, T1}, vec(points)), (length(points) ÷ 3, ))
    SparseVoxelGrid(new_data[1:end], get_voxel_size(voxel_size))
end

function SparseVoxelGrid(points::Vector{T1}, voxel_size) where T1 <: AbstractVector
    SparseVoxelGrid(points, get_voxel_size(voxel_size))
end

get_voxel_size(voxel_size::NTuple{3, T}) where T <: Real  = SVector{3, T}(voxel_size)
get_voxel_size(voxel_size::T) where T <: Real             = SVector{3, T}(voxel_size, voxel_size, voxel_size)
get_voxel_size(voxel_size::SVector{3, T}) where T <: Real = voxel_size

Base.length(grid::SparseVoxelGrid) = length(grid.voxel_info)
Base.isempty(grid::SparseVoxelGrid) = isempty(grid.voxel_info)
Base.haskey(grid::SparseVoxelGrid, k) = haskey(grid.voxel_info, k)

function Base.show(io::IO, grid::SparseVoxelGrid)
    println(io, typeof(grid))
    println(io, "  Number of voxels: ", length(grid))
    println(io, "  Number of points in grid: ", length(grid.point_indices))
    print(io, "  Side length per dimension: ", collect(grid.voxel_size))
end

"""
    make_voxel_id(point::AbstractVector, voxel_size::SVector{3,AbstractFloat})

Create the voxel id for a given point and voxel size.
"""
@inline function make_voxel_id(point::AbstractVector, voxel_size::SVector{3, T}) where T <: Real
    (floor(Int, point[1] / voxel_size[1]), floor(Int, point[2] / voxel_size[2]),
     floor(Int, point[3] / voxel_size[3]))
end

"An iterator type to return point indices in a voxel. See SparseVoxelGrid() for usage."
struct Voxel
    id::VoxelId
    point_index_range::UnitRange{Int}
    all_point_indices::Vector{Int}
end

function Base.iterate(v::SparseVoxelGrid, state=(v.voxel_info, 1))
    state = iterate(v.voxel_info, state[2])
    state == nothing && return nothing
    id = state[1][1]
    point_index_range = state[1][2]
    return Voxel(id, point_index_range, v.point_indices), state
end

Base.eltype(::SparseVoxelGrid) = Voxel

function Base.getindex(grid::SparseVoxelGrid, id::VoxelId)
    Voxel(id, grid.voxel_info[id], grid.point_indices)
end

function Base.iterate(v::Voxel, state=(v,1))
    if state > length(v.point_index_range)
        return nothing
    end
    return v.all_point_indices[v.point_index_range[state[2]]], state[2] + 1
end

Base.eltype(::Voxel) = Int
Base.length(v::Voxel) = length(v.point_index_range)

function Base.show(io::IO, v::Voxel)
    print(io, typeof(v), " ", v.id, " with ", length(v.point_index_range), " points")
end

"Voxel iterator that returns the `Voxel`s. See `in_cuboid()` for usage."
struct VoxelCuboid
    grid::SparseVoxelGrid
    voxel_id::VoxelId
    range::CartesianIndices{3,Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}}}

end

# TODO the `do` syntax for in_cuboid is faster than the iterator - can the iterator be improved?

"""
    in_cuboid(grid::SparseVoxelGrid, voxel::Voxel, radius::Int)
    in_cuboid(grid::SparseVoxelGrid, voxel_id::NTuple{3,Int}, radius::Int)

Search for neighbouring voxels within a `radius` around the reference `voxel` or `voxel_id`.
Returns a `Voxel` in each iteration.

### Example
The `in_cuboid` function can be implemented using the `do` block syntax:

```julia
radius = 1
query_voxel = (1,1,1)
in_cuboid(grid, query_voxel, radius) do voxel
    for index in voxel
        # Do stuff with point[:, index]
    end
    # Or, collect all indices into an array
    indices = collect(voxel)
end
```

Alternatively, you may use a `for` loop which returns a voxel in each iteratation:
```julia
for voxel in in_cuboid(grid, query_voxel, radius)
    # do stuff with the `Voxel` (i.e. collect(voxel) or for index in voxel etc.)
end
```
"""
in_cuboid(grid::SparseVoxelGrid, voxel::Voxel, radius::Int) = in_cuboid(grid, voxel.id, radius)

function in_cuboid(grid::SparseVoxelGrid, voxel::VoxelId, radius::Int)
    index_start = (-radius+voxel[1], -radius+voxel[2], -radius+voxel[3])
    index_stop = (radius+voxel[1], radius+voxel[2], radius+voxel[3])
    indices = CartesianIndices((index_start[1]:index_stop[1], index_start[2]:index_stop[2], index_start[3]:index_stop[3]))
    VoxelCuboid(grid, voxel, indices)
end

in_cuboid(f::Function, grid::SparseVoxelGrid, voxel::Voxel, radius::Int) = in_cuboid(f, grid, voxel.id, radius)

function in_cuboid(f::Function, grid::SparseVoxelGrid, voxel::VoxelId, radius::Int)
    for i=-radius+voxel[1]:radius+voxel[1], j=-radius+voxel[2]:radius+voxel[2], k=-radius+voxel[3]:radius+voxel[3]
        id = (i, j, k)
        if haskey(grid, id) && id != voxel
           f(grid[id])
        end
    end
end

function Base.getindex(c::VoxelCuboid, id::CartesianIndex{3})
    Voxel(id.I, c.grid.voxel_info[id.I], c.grid.point_indices)
end

#find the next viable Voxel
function next_voxel(c::VoxelCuboid, state::Tuple{CartesianIndex{3}, Int})
    next_state = state
    next_state[1] == nothing && return nothing , 0
    voxel = c[next_state[1]]
    while next_state != nothing
        next_state = iterate(c.range, next_state[1])
        next_state == nothing && return voxel, (next_state, 0)
        id = next_state[1]

        if haskey(c.grid.voxel_info, next_state[1].I) && c.voxel_id != next_state[1].I
            # return current voxel and the state for the next voxel
            return voxel, (next_state[1], 1)
        end
    end
    # Next voxel does not exist exists
    return voxel, (next_state, 0)
end

function Base.iterate(c::VoxelCuboid, state=(c.range[1], 1))

    if state[1] == c.range[1]
        if !haskey(c.grid, state[1].I) # first voxel id is not in grids
            # find the next voxel in grid
            while state != nothing
                state = iterate(c.range, state[1])
                state == nothing && return nothing
                id = state[1]
                if haskey(c.grid, id.I) && c.voxel_id != id.I
                    state =  id, 1
                    @show state
                    return next_voxel(c, state)
                end
            end
            # no voxel id was found set value to quit iterations
            return nothing
        end
        # return the starting voxel
        return next_voxel(c, state)

    else
        state[1] == nothing && return nothing
        return next_voxel(c, state)
    end

end

Base.eltype(::VoxelCuboid) = Voxel

if VERSION >= v"0.5.0-dev+3305"
    # See https://github.com/JuliaLang/julia/issues/15977
    # Possibly could implement length() instead, but it's nontrivial work to compute.
    Base.IteratorSize(::Type{VoxelCuboid}) = Base.SizeUnknown()
end

function Base.show(io::IO, c::VoxelCuboid)
    len = length(c.range)
    print(io, typeof(c), " ID iteration range: ", c.range[1].I, " -> ", c.range[len].I)
end

"""
    voxel_center(grid::SparseVoxelGrid, voxel_id::NTuple{3,Int})

Calculate the centre point for the `voxel_id` in the spatial grid.
"""
@inline function voxel_center(grid::SparseVoxelGrid, voxel_id::VoxelId)
    center = SVector{3, Float64}(voxel_id) .* grid.voxel_size - grid.voxel_size * 0.5
end


text:
using Distributed, StatsBase, Distributions, PyPlot
using TransitionalMCMC

lb  = -15
ub  = 15

fT(x) = logpdf(Uniform(lb,ub), x[1]) .+ logpdf(Uniform(lb,ub), x[2])
sample_fT(Nsamples) = rand(Uniform(lb,ub), Nsamples, 2)

log_fD_T(x) = log.(pdf(MvNormal([0,0],[1 -0.5; -0.5 1]), x) + pdf(MvNormal([5,5],[1 0.5; 0.5 1]), x) + pdf(MvNormal([-5,5],[1 0.9; 0.9 1]), x))

samps, acc =tmcmc(log_fD_T, fT, sample_fT, 2000)

plt.scatter(samps[:,1], samps[:,2])


text:
################################################################################
#
#  Line enumeration
#
################################################################################

# Iterate over the lines in K^n, that is, over the points of projective
# space P^(n-1)(K).
#
# Important: In the prime case, this must always be lexicographically ordered

mutable struct LineEnumCtx{T, S}
  K::T
  a::S # primitive element
  dim::Int
  depth::Int
  v::Vector{S}
  length::BigInt
end

# We iterate through K by adding 1 in the prime field case and by multiplying
# with a primitive element in the general case.

function LineEnumCtx(K::T, n) where {T}
  a = primitive_element(K)
  v = Vector{elem_type(K)}(undef, n)
  for i in 1:n
    v[i] = zero(K)
  end
  depth = n + 1
  dim = n
  q = order(K)
  length = divexact(q^n - 1, q - 1)
  return LineEnumCtx{T, elem_type(T)}(K, a, dim, depth, v, length)
end

function LineEnumCtx(K::T, n::Int) where {T <: Union{GaloisField, GaloisFmpzField}}
  a = zero(K)
  v = Vector{elem_type(T)}(undef, n)
  for i in 1:n
    v[i] = zero(K)
  end
  depth = n + 1
  dim = n
  q = order(K)
  length = divexact(q^n - 1, q - 1)
  return LineEnumCtx{T, elem_type(T)}(K, a, dim, depth, v, length)
end

function enumerate_lines(K, n)
  return LineEnumCtx(K, n)
end

function Base.show(io::IO, P::LineEnumCtx)
  print(io, "Iterator for affine lines in K^n with n = ", dim(P), "\n")
  print(io, "over ", P.K)
end

Base.length(P::LineEnumCtx) = P.length

Base.eltype(::Type{LineEnumCtx{T, S}}) where {T, S} = Vector{S}

depth(P::LineEnumCtx) = P.depth

dim(P::LineEnumCtx) = P.dim

primitive_element(P::LineEnumCtx) = P.a

################################################################################
#
#  Iteration
#
################################################################################

# custom iteration
function next(P::LineEnumCtx)
  if depth(P) > 0
    i = dim(P)
    while true
      @show i, P.v, depth(P)
      if i == depth(P)
        P.v[i] = zero!(P.v[i])
        i = i - 1
      elseif i < depth(P)
        P.depth = i
        if i >= 1
          P.v[i] = one(P.K)
        end
        break
      elseif iszero(P.v[i])
        P.v[i] = one(P.K)
        break
      else
        P.v[i] = P.v[i] * primitive_element(P)
        if isone(P.v[i])
          P.v[i] = zero!(P.v[i])
          i = i - 1
        else
          break
        end
      end
    end
  end
  return P.v
end

function next(P::LineEnumCtx{T, S}) where {T <: Union{GaloisField, GaloisFmpzField}, S}
  if depth(P) > 0
    i = dim(P)
    while true
      @show i, P.v, depth(P)
      if i == depth(P)
        P.v[i] = zero!(P.v[i])
        i = i - 1
      elseif i < depth(P)
        P.depth = i
        if i >= 1
          P.v[i] = one(P.K)
        end
        break
      elseif iszero(P.v[i])
        P.v[i] = one(P.K)
        break
      else
        P.v[i] = P.v[i] + 1
        if iszero(P.v[i])
          i = i - 1
        else
          break
        end
      end
    end
  end
  return P.v
end

################################################################################
#
#  Iterator interface
#
################################################################################

function Base.iterate(P::LineEnumCtx)
  P.v[dim(P)] = one(P.K)
  P.depth = dim(P)
  return P.v, one(BigInt)
end

# For prime fields
function Base.iterate(P::LineEnumCtx{T, S}, state) where
                                  {T <: Union{GaloisField, GaloisFmpzField}, S}
  if state >= P.length
    return nothing
  end

  onee = one(P.K)

  if depth(P) > 0
    i = dim(P)
    while true
      if i == depth(P)
        P.v[i] = zero(P.K)
        i = i - 1
      elseif i < depth(P)
        P.depth = i
        if i >= 1
          P.v[i] = onee
        end
        break
      elseif iszero(P.v[i])
        P.v[i] = onee
        break
      else
        P.v[i] = P.v[i] + 1
        if iszero(P.v[i])
          i = i - 1
        else
          break
        end
      end
    end
  end
  return P.v, inc!(state)
end

# For non-prime fields
function Base.iterate(P::LineEnumCtx, state)
  if state >= P.length
    return nothing
  end

  if depth(P) > 0
    i = dim(P)
    while true
      if i == depth(P)
        P.v[i] = zero!(P.v[i])
        i = i - 1
      elseif i < depth(P)
        P.depth = i
        if i >= 1
          P.v[i] = one(P.K)
        end
        break
      elseif iszero(P.v[i])
        P.v[i] = one(P.K)
        break
      else
        P.v[i] = P.v[i] * primitive_element(P)
        if isone(P.v[i])
          P.v[i] = zero!(P.v[i])
          i = i - 1
        else
          break
        end
      end
    end
  end
  return P.v, inc!(state)
end

################################################################################
#
#  Line orbits
#
################################################################################

line_orbits(G::Vector{gfp_mat}) = _line_orbits(G)

function line_orbits(G::Vector{T}) where {T <: MatElem{gfp_fmpz_elem}}
  F = base_ring(G[1])
  p = modulus(F)
  k = nrows(G[1])
  if fits(UInt, p)
    Fsmall = GF(UInt(p), cached = false)
    GG = Vector{dense_matrix_type(Fsmall)}(undef, length(G))
    let Fsmall = Fsmall
      for i in 1:length(G)
        GG[i] = map_entries(x -> Fsmall(x.data), G[i])::dense_matrix_type(Fsmall)
      end
    end
    _res = line_orbits(GG)
    res = Vector{Tuple{Vector{elem_type(F)}, Int}}(undef, length(_res))
    for i in 1:length(_res)
      t = _res[i][1]
      y = Vector{elem_type(F)}(undef, k)
      for j in 1:k
        y[j] = F(t[j].data)
      end
      res[i] = y, _res[i][2]
    end
    return res
  else
    return _line_orbits(G)
  end
end

function line_orbits(G::Vector{fq_nmod_mat})
  F = base_ring(G[1])
  d = degree(F)
  k = nrows(G[1])
  if d == 1
    p = characteristic(F)
    @assert fits(UInt, p)
    GFp = GF(UInt(p), cached = false)
    GG = Vector{gfp_mat}(undef, length(G))
    let GFp = GFp
      for i in 1:length(G)
        GG[i] = map_entries(x -> GFp(coeff(x, 0)), G[i])
      end
    end
    _res = line_orbits(GG)
    res = Vector{Tuple{Vector{elem_type(F)}, Int}}(undef, length(_res))
    for i in 1:length(_res)
      t = _res[i][1]
      y = Vector{elem_type(F)}(undef, k)
      for j in 1:k
        y[j] = F(t[j].data)
      end
      res[i] = y, _res[i][2]
    end
    return res
  else
    return _line_orbits(G)
  end
end

function line_orbits(G::Vector{fq_mat})
  F = base_ring(G[1])
  a = gen(F)
  d = degree(F)
  p = characteristic(F)
  k = nrows(G[1])
  if fits(UInt, p)
    f = defining_polynomial(F)
    GFp = GF(UInt(p), cached = false)
    GFpx, x = PolynomialRing(GFp, "x", cached = false)
    local fp
    let GFp = GFp
      fp = map_coeffs(a -> GFp(a.data), f, parent = GFpx)
    end
    FF, aa = FlintFiniteField(fp, "aa", cached = false)
    GG = Vector{dense_matrix_type(FF)}(undef, length(G))
    for i in 1:length(G)
      GG[i] = map_entries(b -> sum(coeff(b, i) * aa^i for i in 0:(d-1)), G[i])::dense_matrix_type(FF)
    end
    _res = line_orbits(GG)
    res = Vector{Tuple{Vector{fq}, Int}}(undef, length(_res))
    for i in 1:length(res)
      t = _res[i][1]
      y = Vector{elem_type(F)}(undef, k)
      for j in 1:k
        y[j] = sum(coeff(t[j], i) * a^i for i in 0:(d - 1))
      end
      res[i] = y, _res[i][2]
    end
    return res
  else
    return _line_orbits(G)
  end
end

function _line_orbits(G::Vector)
  K = base_ring(G[1])
  n = nrows(G[1])
  P = enumerate_lines(K, n)
  l = length(P)
  @vprint :GenRep 1 "Computing orbits of lines of total length $l over $K\n"
  lines = Vector{eltype(P)}(undef, l)
  i = 1
  for v in P
    lines[i] = deepcopy(v)
    i += 1
  end

  if !(K isa GaloisField && K isa GaloisFmpzField)
    sort!(lines, lt = _isless)
  end

  res = Tuple{eltype(P), Int}[]

  visited = trues(l)
  sofar = zero(BigInt)
  newline = zero_matrix(K, 1, n)
  newline2 = zero_matrix(K, 1, n)
  ___isless = __isless(elem_type(K)) # the comparison needs some
                                     # temporary variables in some cases
  while sofar < l
    pt = findfirst(visited)
    @assert pt !== nothing
    visited[pt] = false
    norb = 1
    cnd = 1
    orb = Int[pt]
    while cnd <= norb
      set!(newline, lines[orb[cnd]])
      for i in 1:length(G)
        newline2 = mul!(newline2, newline, G[i])
        _normalize!(newline2)
        m = searchsortedfirst(lines, newline2, lt = ___isless)
        @assert m !== nothing
        if visited[m]
          visited[m] = false
          norb += 1
          push!(orb, m)
        end
      end
      cnd += 1
    end
    push!(res, (lines[pt], norb))
    sofar = sofar + norb
  end
  return res
end

################################################################################
#
#  Cheap increment of a BigInt
#
################################################################################

function inc!(a::BigInt)
  ccall((:__gmpz_add_ui, :libgmp), Nothing,
        (Ref{BigInt}, Ref{BigInt}, UInt), a, a, UInt(1))
  return a
end

################################################################################
#
#  Helper functions for sorting and comparing finite field matrices and arrays
#
################################################################################

# gfp
function set!(a::gfp_mat, z::Vector{gfp_elem})
  @assert ncols(a) == length(z)
  for i in 1:length(z)
    Nemo.setindex_raw!(a, z[i].data, 1, i) # no reduction necessary
  end
end

function _normalize!(x::gfp_mat)
  R = base_ring(x)
  @GC.preserve x begin
    piv = 0
    local ell
    for j in 1:ncols(x)
      el = x[1, j]
      if !iszero(el)
        piv = j
        ell = inv(el)
        break
      end
    end

    @assert piv != 0

    for j in (piv+1):ncols(x)
      el = x[j] * ell
      Nemo.setindex_raw!(x, el.data, 1, j)
    end
    Nemo.setindex_raw!(x, UInt(1), 1, piv)
  end
  return x
end

function _isless(x::Vector{gfp_elem}, y::gfp_mat)
  d = length(x)
  for i in 1:d
    xi = x[i]
    yi = y[1, i]
    if xi.data != yi.data
      return xi.data < yi.data
    end
  end
  return false
end

function _isless(x::Vector{gfp_elem}, y::Vector{gfp_elem})
  d = length(x)
  for i in 1:d
    xi = x[i]
    yi = y[i]
    if xi.data != yi.data
      return xi.data < yi.data
    end
  end
  return false
end

# gfp_fmpz

function set!(a::MatElem, z::Vector{gfp_fmpz_elem})
  @assert ncols(a) == length(z)
  for i in 1:length(z)
    a[1, i] = z[i]
  end
end

function _normalize!(x::MatElem{gfp_fmpz_elem})
  R = base_ring(x)
  @GC.preserve x begin
    piv = 0
    local ell
    for j in 1:ncols(x)
      el = x[1, j]
      if !iszero(el)
        piv = j
        ell = inv(el)
        break
      end
    end

    @assert piv != 0

    for j in (piv+1):ncols(x)
      t = x[1, j]
      mul!(t, t, ell)
    end
    one!(x[1, piv].data)
  end
  return x
end

function _isless(x::Vector{gfp_fmpz_elem}, y::MatElem)
  d = length(x)
  for i in 1:d
    xi = x[i]
    yi = y[1, i]
    if xi.data != yi.data
      return xi.data < yi.data
    end
  end
  return false
end

function _isless(x::Vector{gfp_fmpz_elem}, y::Vector{gfp_fmpz_elem})
  d = length(x)
  for i in 1:d
    xi = x[i]
    yi = y[i]
    if xi.data != yi.data
      return xi.data < yi.data
    end
  end
  return false
end

# fq_nmod_mat
function set!(a::fq_nmod_mat, z::Vector{fq_nmod})
  @assert ncols(a) == length(z)
  for i in 1:length(z)
    ccall((:fq_nmod_mat_entry_set, libflint), Nothing,
          (Ref{fq_nmod_mat}, Int, Int, Ref{fq_nmod}, Ref{FqNmodFiniteField}),
          a, 0, i - 1, z[i], base_ring(a))
  end
end

function _isequal(x::fq_nmod_mat, y::Vector{fq_nmod})
  R = base_ring(x)
  @GC.preserve x begin
    for i in 1:length(y)
      el = ccall((:fq_nmod_mat_entry, libflint), Ptr{fq_nmod}, (Ref{fq_nmod_mat}, Int, Int), x, 0, i - 1)
      b = ccall((:fq_nmod_equal, libflint), Cint, (Ref{fq_nmod}, Ptr{fq_nmod}, Ref{FqNmodFiniteField}), y[i], el, R)
      bb = Bool(b)
      if !bb
        return false
      end
    end
  end
  return true
end

function _muleq!(x::fq_nmod_mat, y::fq_nmod)
  R = base_ring(x)
  @GC.preserve x begin
    for i in 1:nrows(x)
      for j in 1:ncols(x)
        el = ccall((:fq_nmod_mat_entry, libflint), Ptr{fq_nmod}, (Ref{fq_nmod_mat}, Int, Int), x, i - 1, j - 1)
        ccall((:fq_nmod_mul, libflint), Cvoid, (Ptr{fq_nmod}, Ptr{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), el, el, y, R)
      end
    end
  end
  return x
end

function _normalize!(x::fq_nmod_mat)
  R = base_ring(x)
  @GC.preserve x begin
    piv = 0
    local ell
    for j in 1:ncols(x)
      el = ccall((:fq_nmod_mat_entry, libflint), Ptr{fq_nmod}, (Ref{fq_nmod_mat}, Int, Int), x, 0, j - 1)
      b = ccall((:fq_nmod_is_zero, libflint), Cint, (Ptr{fq_nmod}, Ref{FqNmodFiniteField}), el, R)
      if !Bool(b)
        piv = j
        ccall((:fq_nmod_inv, libflint), Cvoid, (Ptr{fq_nmod}, Ptr{fq_nmod}, Ref{FqNmodFiniteField}), el, el, R)
        ell = el
        break
      end
    end

    @assert piv != 0

    for j in (piv+1):ncols(x)
      el = ccall((:fq_nmod_mat_entry, libflint), Ptr{fq_nmod},
                 (Ref{fq_nmod_mat}, Int, Int), x, 0, j - 1)
      ccall((:fq_nmod_mul, libflint), Cvoid,
            (Ptr{fq_nmod}, Ptr{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}),
            el, el, ell, R)
    end
    ccall((:fq_nmod_one, libflint), Cvoid,
          (Ptr{fq_nmod}, Ref{FqNmodFiniteField}), ell, R)
  end
  return x
end

function _isless(x::fq_nmod, y::fq_nmod)
  d = degree(parent(x)) - 1
  for i in 0:d
    xi = coeff(x, i)
    yi = coeff(y, i)
    if xi != yi
      return xi < yi
    end
  end
  return false
end

@inline function _isless(x::fq_nmod, y::Ptr{fq_nmod})
  d = degree(parent(x)) - 1
  for i in 0:d
    xi = coeff(x, i)
    yi = ccall((:nmod_poly_get_coeff_ui, libflint), UInt,
               (Ptr{fq_nmod}, Int), y, i)
    if xi != yi
      return xi < yi
    end
  end
  return false
end

function _isless(x::Vector{fq_nmod}, y::Vector{fq_nmod})
  d = length(x)
  for i in 1:d
    xi = x[i]
    yi = y[i]
    if xi != yi
      return _isless(xi, yi)
    end
  end
  return false
end

function _isless(x::Vector{fq_nmod}, y::fq_nmod_mat)
  d = length(x)
  R = base_ring(y)
  @GC.preserve y begin
    for i in 1:d
      xi = x[i]
      el = ccall((:fq_nmod_mat_entry, libflint), Ptr{fq_nmod},
                 (Ref{fq_nmod_mat}, Int, Int), y, 0, i - 1)
      b = ccall((:fq_nmod_equal, libflint), Cint,
                (Ref{fq_nmod}, Ptr{fq_nmod}, Ref{FqNmodFiniteField}), xi, el, R)
      if !Bool(b)
        return _isless(xi, el)
      end
    end
    return false
  end
end

# fq
function set!(a::fq_mat, z::Vector{fq})
  @assert ncols(a) == length(z)
  for i in 1:length(z)
    ccall((:fq_mat_entry_set, libflint), Nothing,
          (Ref{fq_mat}, Int, Int, Ref{fq}, Ref{FqFiniteField}),
          a, 0, i - 1, z[i], base_ring(a))
  end
end

function _isequal(x::fq_mat, y::Vector{fq})
  R = base_ring(x)
  @GC.preserve x begin
    for i in 1:length(y)
      el = ccall((:fq_mat_entry, libflint), Ptr{fq},
                 (Ref{fq_mat}, Int, Int), x, 0, i - 1)
      b = ccall((:fq_equal, libflint), Cint,
                (Ref{fq}, Ptr{fq}, Ref{FqFiniteField}), y[i], el, R)
      bb = Bool(b)
      if !bb
        return false
      end
    end
  end
  return true
end

function _muleq!(x::fq_mat, y::fq)
  R = base_ring(x)
  @GC.preserve x begin
    for i in 1:nrows(x)
      for j in 1:ncols(x)
        el = ccall((:fq_mat_entry, libflint), Ptr{fq},
                   (Ref{fq_mat}, Int, Int), x, i - 1, j - 1)
        ccall((:fq_mul, libflint), Cvoid,
              (Ptr{fq}, Ptr{fq}, Ref{fq}, Ref{FqFiniteField}), el, el, y, R)
      end
    end
  end
  return x
end

function _normalize!(x::fq_mat)
  R = base_ring(x)
  @GC.preserve x begin
    piv = 0
    local ell
    for j in 1:ncols(x)
      el = ccall((:fq_mat_entry, libflint), Ptr{fq},
                 (Ref{fq_mat}, Int, Int), x, 0, j - 1)
      b = ccall((:fq_is_zero, libflint), Cint,
                (Ptr{fq}, Ref{FqFiniteField}), el, R)
      if !Bool(b)
        piv = j
        ccall((:fq_inv, libflint), Cvoid,
              (Ptr{fq}, Ptr{fq}, Ref{FqFiniteField}), el, el, R)
        ell = el
        break
      end
    end

    @assert piv != 0

    for j in (piv+1):ncols(x)
      el = ccall((:fq_mat_entry, libflint), Ptr{fq},
                 (Ref{fq_mat}, Int, Int), x, 0, j - 1)
      ccall((:fq_mul, libflint), Cvoid,
            (Ptr{fq}, Ptr{fq}, Ref{fq}, Ref{FqFiniteField}), el, el, ell, R)
    end
    ccall((:fq_one, libflint), Cvoid, (Ptr{fq}, Ref{FqFiniteField}), ell, R)
  end
  return x
end

function _isless(x::fq, y::fq, xi::fmpz = fmpz(), yi::fmpz = fmpz())
  d = degree(parent(x)) - 1
  for i in 0:d
    ccall((:fmpz_poly_get_coeff_fmpz, libflint), Cvoid,
          (Ref{fmpz}, Ref{fq}, Int), xi, x, i)
    ccall((:fmpz_poly_get_coeff_fmpz, libflint), Cvoid,
          (Ref{fmpz}, Ref{fq}, Int), yi, y, i)
    if xi != yi
      return xi < yi
    end
  end
  return false
end

@inline function _isless(x::fq, y::Ptr{fq}, xi::fmpz = fmpz(),
                                            yi::fmpz = fmpz())
  d = degree(parent(x)) - 1
  for i in 0:d
    ccall((:fmpz_poly_get_coeff_fmpz, libflint), Cvoid,
          (Ref{fmpz}, Ref{fq}, Int), xi, x, i)
    ccall((:fmpz_poly_get_coeff_fmpz, libflint), Cvoid,
          (Ref{fmpz}, Ptr{fq}, Int), yi, y, i)
    if xi != yi
      return xi < yi
    end
  end
  return false
end

function _isless(x::Vector{fq}, y::Vector{fq}, tx::fmpz = fmpz(),
                                               ty::fmpz = fmpz())
  d = length(x)
  for i in 1:d
    xi = x[i]
    yi = y[i]
    if xi != yi
      return _isless(xi, yi, tx, ty)
    end
  end
  return false
end

function _isless(x::Vector{fq}, y::fq_mat, tx::fmpz = fmpz(), ty::fmpz = fmpz())
  d = length(x)
  R = base_ring(y)
  @GC.preserve y begin
    for i in 1:d
      xi = x[i]
      el = ccall((:fq_mat_entry, libflint), Ptr{fq},
                 (Ref{fq_mat}, Int, Int), y, 0, i - 1)
      b = ccall((:fq_equal, libflint), Cint,
                (Ref{fq}, Ptr{fq}, Ref{FqFiniteField}), xi, el, R)
      if !Bool(b)
        return _isless(xi, el, tx, ty)
      end
    end
    return false
  end
end

# In some cases the _isless needs some temporary variables

function __isless(::Type{fq})
  t1, t2 = fmpz(), fmpz()
  return (x, y) -> _isless(x, y, t1, t2)
end

__isless(::Type{T}) where {T} = _isless

################################################################################
#
#  Defining polynomial for finite fields
#
################################################################################

defining_polynomial(F::FqFiniteField) = minpoly(gen(F))


text:
using Documenter, ArpackMKL

makedocs(
    format = :html,
    sitename = "ArpackMKL.jl",
    modules = [ArpackMKL],
    pages = [
        "index.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaLinearAlgebra/ArpackMKL.jl.git",
    target = "build",
    julia  = "nightly",
    deps = nothing,
    make = nothing
)


text:
using Cbc
using Hungarian
using JuMP
using Kombinator
using LightGraphs
using Munkres
using Test

is_travis = "GITHUB_ACTIONS" in keys(ENV)
# TODO: need for optional dependencies for this to work, I suppose.
# https://github.com/JuliaLang/Pkg.jl/issues/1285
if !is_travis
    using JuMP
    # Why Gurobi?
    # - Need support for lazy constraints (elementary paths):
    #   - Mosek does not support them
    # - Need support for MISOCP
    #   - Pajarito is not yet ported to MOI
    using Gurobi
end

@testset "Kombinator.jl" begin
    include("um.jl")
    include("st.jl")
    include("ep.jl")
    # include("matching.jl")

    include("nlcop.jl")
end


text:
struct CuDeviceArray{T,N} <: AbstractArray{T,N}
    shape::NTuple{N,Int}
    ptr::Ptr{T}
end

@inline function Base.setindex!(A::CuDeviceArray{T}, x, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    Base.unsafe_store!(A.ptr, x, index)
end

@inline function Base.getindex(A::CuDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    Base.unsafe_load(A.ptr, index)::T
end

Base.size(g::CuDeviceArray) = g.shape

Base.checkbounds(::CuDeviceArray, I...) = nothing

function kernel_vadd(a, b, c, i)
    c[i] = a[i] + b[i]

    return nothing
end

code_llvm(stdout, kernel_vadd, Tuple{CuDeviceArray{Float32,1},
                                     CuDeviceArray{Float32,1},
                                     CuDeviceArray{Float32,1},
                                     Int})

using TinyCassette
struct Ctx end

code_llvm(stdout, TinyCassette.execute,
          Tuple{Ctx, typeof(kernel_vadd),
                CuDeviceArray{Float32,1},
                CuDeviceArray{Float32,1},
                CuDeviceArray{Float32,1},
                Int})