text:
using CARPData
using Test

@testset "kshs" begin
    data = load(:kshs3)
    @test data.lb == 9320
    @test data.ub == 9320
    @test data.name == "kshs3"
    @test length(data.vertices) == 6
    @test length(data.edges) == 15
    @test length(data.requireds) == 15
end

@testset "gdb" begin
    data = load(:gdb8)
    @test data.lb == 348
    @test data.ub == 348
    @test data.name == "gdb8"
    @test length(data.vertices) == 27
    @test length(data.edges) == 46
    @test length(data.requireds) == 46
end

@testset "val" begin
    data = load(:val9B)
    @test data.lb == 326
    @test data.lb == 326
    @test data.name == "val9B"
    @test length(data.vertices) == 50
    @test length(data.edges) == 92
    @test length(data.requireds) == 92
end

@testset "egl" begin
    data = load(:egl_s2_A)
    @test data.lb == 9838
    @test data.ub == 9875
    @test data.name == "egl-s2-A"
    @test length(data.vertices) == 140
    @test length(data.edges) == 190
    @test length(data.requireds) == 147
end

@testset "beullens" begin
    data = load(:F16)
    @test data.lb == 2725
    @test data.ub == 2725
    @test data.name == "F16"
    @test length(data.vertices) == 60
    @test length(data.edges) == 80
    @test length(data.requireds) == 54
end

@testset "egl-large" begin
    data = load(:egl_g1_E)
    @test data.lb == 1461469
    @test data.ub == 1506493
    @test data.name == "egl-g1-E"
    @test length(data.vertices) == 255
    @test length(data.edges) == 375
    @test length(data.requireds) == 347
end

@testset "other" begin
    data = load(:A10A)
    @test data.lb == 80
    @test data.ub == 93
    @test data.name == "A10A"
    @test length(data.vertices) == 10
    @test length(data.edges) == 15
    @test length(data.requireds) == 11
end

@testset "Load string" begin
    data = load(joinpath(@__DIR__, "../data/kshs1.dat"))
    @test data !== nothing
    data = load("not a file")
    @test data === nothing
end

@testset "Base.show methods" begin
    data = load(:kshs3)
    show(stderr, data)
    print("\n")
    show(stderr, data.edges[1])
    print("\n")
    show(stderr, data.vertices[1])
    print("\n")
end


text:
function try_fetch(future :: Future)
  ret = fetch(future)
  if ret isa RemoteException
    throw(ret)
  end
  return ret
end

ext: jl


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using DecisionTrees
using RLESUtils, Observers, Loggers
import DecisionTreeVis: get_tree, get_metric

function DecisionTrees.get_truth{T}(members::Vector{Int64},
                                 Dl::DFSetLabeled{T}, otherargs...) #userargs...
  return labels(Dl, members)
end

function classify(problem::ACASXClustering, result::SAESResult, Ds::Vector{DataFrame})
  f = to_function(problem, result.expr)
  return map(f, Ds)
end

function DecisionTrees.get_splitter{T}(members::Vector{Int64},
                                       Dl::DFSetLabeled{T}, problem::ACASXClustering,
                                       psa_params::PSAESParams, logs::TaggedDFLogger, loginterval::Int64) #userargs...
  set_observers!(psa_params.observer, psa_params.sa_params.observer, logs, loginterval)

  problem.Dl = Dl_sub = Dl[members] #fitness function uses problem.Dl
  result = exprsearch(psa_params, problem)

  push_members!(logs, problem, result.expr)
  @notify_observer(psa_params.observer, "expression",
                   [string(result.expr),
                    pretty_string(result.tree, FMT_PRETTY),
                    pretty_string(result.tree, FMT_NATURAL, true)])

  predicts = classify(problem, result, getrecords(Dl_sub))
  info_gain, _, _ = gini_metrics(predicts, labels(Dl_sub))

  return info_gain > 0 ? result : nothing #split object
end

function push_members!{T}(logs::TaggedDFLogger, problem::ACASXClustering{T}, expr)
  decision_id = nrow(logs["members"]) > 0 ?
    maximum(logs["members"][:decision_id]) + 1 : 1
  members_true, members_false = get_members(problem, expr)
  push!(logs, "members", [join(members_true, ","), join(members_false, ","), decision_id])
end

function DecisionTrees.get_labels{T}(result::SearchResult, members::Vector{Int64},
                                     Dl::DFSetLabeled{T}, problem::ACASXClustering, otherargs...) #userargs...
  return classify(problem, result, getrecords(Dl, members))
end

DecisionTreeVis.get_tree(result::SAESResult) = result.tree
DecisionTreeVis.get_metric(result::SAESResult) = result.fitness


text:
# Model
include_model("hopper")

# Dimensions
nq = 4 # configuration dimension
nu = 2 # control dimension
nc = 1 # number of contact points
nf = 2 # number of faces for friction cone
nb = nc * nf
ns = nq

# Parameters
g = 9.81 # gravity
μ = 1.0  # coefficient of friction
mb = 10.0 # body mass
ml = 1.0  # leg mass
Jb = 2.5 # body inertia
Jl = 0.25 # leg inertia

n = 2 * nq
m = nu + nc + nb + nc + nb + ns
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

model = Hopper(n, m, d,
			   mb, ml, Jb, Jl,
			   μ, g,
			   qL, qU,
			   nq,
		       nu,
		       nc,
		       nf,
		       nb,
		   	   ns,
		       idx_u,
		       idx_λ,
		       idx_b,
		       idx_ψ,
		       idx_η,
		       idx_s)

function fd(model::Hopper, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)
	s = view(u, model.idx_s)

	[q2⁺ - q2⁻;
	((1.0 / h) * (M_func(model, q1) * (SVector{4}(q2⁺) - SVector{4}(q1))
	- M_func(model, q2⁺) * (SVector{4}(q3) - SVector{4}(q2⁺)))
	+ transpose(B_func(model, q3)) * SVector{2}(u_ctrl)
	+ transpose(N_func(model, q3)) * SVector{1}(λ)
	+ transpose(P_func(model, q3)) * SVector{2}(b)
	- h * G_func(model, q2⁺)) + s]
end

# Horizon
T = 11

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_s] .= 0.0
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_ul[model.idx_s] .= 0.0
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
z_h = 0.5
q1 = [0.0, 0.5, 0.0, 0.5]
v1 = [0.0, 0.0, 0.0, 0.0]
v2 = v1 - G_func(model,q1) * h
q2 = q1 + 0.5 * h * (v1 + v2)

x1 = [q1; q1]
qT = [1.0, 0.5, 0.0, 0.5]

xl, xu = state_bounds(model, T, [model.qL; model.qL], [model.qU; model.qU],
    x1 = x1, xT = [Inf*ones(model.nq); qT])

# Objective
include_objective(["velocity", "nonlinear_stage"])
obj_velocity = velocity_objective([Diagonal(ones(model.nq)) for t = 1:T],
	model.nq, idx_angle = (3:3), h = h)
obj_tracking = quadratic_tracking_objective(
    [Diagonal(zeros(model.n)) for t = 1:T],
    [Diagonal([1.0, 1.0, zeros(model.m - model.nu)...]) for t = 1:T-1],
    [zeros(model.n) for t = 1:T],
    [zeros(model.m) for t = 1:T-1])
function l_stage_fh1(x, u, t)
	if t > 1
		return 10.0 * (kinematics(model, view(x, 5:8))[2] - 0.25)^2.0
	else
		return 0.0
	end
end
l_terminal_fh1(x) = 0.0 #* (kinematics(model, view(x, 5:8)) - 0.025)^2.0
obj_fh1 = nonlinear_stage_objective(l_stage_fh1, l_terminal_fh1)

function l_stage_b(x, u, t)
	if t > 1
		return 10.0 * (x[6] - 0.5)^2.0
	else
		return 0.0
	end
end
l_terminal_b(x) = 0.0
obj_b = nonlinear_stage_objective(l_stage_b, l_terminal_b)

# obj_penalty = PenaltyObjective(1.0e5, model.m)
obj = MultiObjective([obj_tracking, obj_velocity, obj_fh1, obj_b])#, obj_penalty])

# Constraints
include_constraints("contact_al")
con_contact = contact_al_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact)

# Trajectory initialization
x0 = configuration_to_state(linear_interpolation(q1, qT, T)) # linear interpolation on state
u0 = [1.0e-3 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
# include_snopt()
# @time z̄ , info = solve(prob, copy(z0),
#  	nlp = :SNOPT7,
# 	tol = 1.0e-3, c_tol = 1.0e-3)
#
# check_slack(z̄, prob)
# x̄, ū = unpack(z̄, prob)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h)


text:
abstract type AbstractOS end
abstract type Unix <: AbstractOS end
abstract type BSD <: Unix end

abstract type Windows <: AbstractOS end
abstract type MacOS <: BSD end
abstract type Linux <: BSD end

if Sys.iswindows()
    const os = Windows
elseif Sys.isapple()
    const os = MacOS
else
    const os = Linux
end

abstract type TestData end

function unzip(::Type{<:BSD}, filename, directory)
    @assert success(`tar -xvf $filename -C $directory`) "Unable to extract $filename to $directory"
end

function unzip(::Type{Windows}, filename, directory)
    path_7z = if Base.VERSION < v"0.7-"
        "$JULIA_HOME/7z"
    else
        sep = Sys.iswindows() ? ";" : ":"
        withenv(
            "PATH" => string(
                joinpath(Sys.BINDIR, "..", "libexec"),
                sep,
                Sys.BINDIR,
                sep,
                ENV["PATH"],
            ),
        ) do
            Sys.which("7z")
        end
    end
    @assert success(`$path_7z x $filename -y -o$directory`) "Unable to extract $filename to $directory"
end

"""
Download SLiDE Data from `branch="master"` name into a "data" folder in given argument path.
Skip the actual download if the folder already exists and force=false.
Defaults to the root of the PowerSystems package.

Returns the downloaded folder name.
"""
function Base.download(
    ::Type{TestData};
    folder::AbstractString = abspath(joinpath(@__DIR__, "..")),
    branch::String = "master",
    force::Bool = false,
)

    if Sys.iswindows()
        SLIDEDATA_URL = "https://github.com/NREL/SLiDEData/archive/$branch.zip"
    else
        SLIDEDATA_URL = "https://github.com/NREL/SLiDEData/archive/$branch.tar.gz"
    end
    directory = abspath(normpath(folder))
    data = joinpath(directory, "data")
    if !isdir(data) || force
        tempfilename = Base.download(SLIDEDATA_URL)
        mkpath(directory)
        unzip(os, tempfilename, directory)
        mv(joinpath(directory, "SLiDEData-$branch"), data, force = force)
    end

    return data
end

download(TestData; branch = "master")

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using TensorOperators
using Base.Test

input_size = 40
output_size = 50
sample_size = 300

m = LinearLayer(cpu, Float64, input_size, output_size)
parameters, gradient = linearize_parameters!(cpu, Float64, m)
init!(m)

# Computes with module and directly

x = randn(sample_size, input_size)

# --- Test forward

y = forward!(m, x)
yp = broadcast(+, m.bias.values, x * m.weight.values)

@test_approx_eq yp m.output

# --- Test of gradient (input)

dx = randn(sample_size, input_size)

y = copy(y)
∂x_∂y = copy(compute_inputgradient!(m, x, eye(output_size)))
@test_approx_eq_eps ∂x_∂y transpose(m.weight.values) (1e-16 * length(m.weight.values))

# Approximate gradient
epsilon = 1e-5
dx = epsilon * randn(size(x))
xdx = x + dx
delta = forward!(m, xdx) - y
@test_approx_eq_eps(dx * transpose(∂x_∂y) - delta, zeros(delta), 1e-16 * length(delta))

# --- Test of gradient (parameters)

init_gradient!(m)
∂x_∂y = copy(compute_inputgradient!(m, x, eye(output_size)))


text:
function MeshCatMechanisms.MechanismVisualizer(mech_scen::MechanismScenario, vis)
	return MechanismVisualizer(mech_scen.float.state.mechanism, vis)
end

function MeshCatMechanisms.set_configuration!(mech_scen::MechanismScenario, mvis::MechanismVisualizer)
	set_configuration!(mvis, mech_scen.float.state.q)
end

function set_body_mesh_visual!(mvis::MechanismVisualizer, mech_scen::MechanismScenario, m_id::MeshID, color)
    mesh = mech_scen.MeshCache[m_id]
    body = bodies(mech_scen.float.state.mechanism)[mesh.BodyID]
    return set_body_mesh_visual!(mvis, mech_scen, body, m_id, color)
end

function set_body_mesh_visual!(mvis::MechanismVisualizer, mech_scen::MechanismScenario, body::RigidBody{Float64},
        m_id::MeshID, color)

    color = as_rgba(color)
    mc = mech_scen.MeshCache[m_id]
	mesh = mc.mesh
	mesh_32 = HomogenousMesh_32(mesh)
    my_vis_ele = VisualElement(default_frame(body), mesh_32, color, Translation(0,0,0))
    setelement!(mvis, my_vis_ele, mc.name)
    return nothing
end

function set_mesh_visual!(mvis::MechanismVisualizer, mech_scen::MechanismScenario, m_id::MeshID, color)
    mc = mech_scen.MeshCache[m_id]
    body = root_body(mech_scen.float.state.mechanism)
    color = as_rgba(color)
	mesh = mc.mesh
	mesh_32 = HomogenousMesh_32(mesh)
    my_vis_ele = VisualElement(default_frame(body), mesh_32, color, Translation(0,0,0))
    setelement!(mvis, my_vis_ele, mc.name)
    return nothing
end

as_rgba(color::RGBA{Float32}) = color
as_rgba(color) = RGBA{Float32}(color...)

get_vertices_32(e_mesh::eMesh{T1,T2})  where {T1,T2} = [Point{3,Float32}(k) for k = e_mesh.point]
get_faces_32(   e_mesh::eMesh{Tri,T2}) where {T2}    = [Face{3,Int32}(k)    for k = e_mesh.tri]

function HomogenousMesh_32(eM::eMesh)
	eM = as_tri_eMesh(eM, false)
	vertices = get_vertices_32(eM)
    faces = get_faces_32(eM)
    return HomogenousMesh(vertices=vertices, faces=faces)
end

# function HomogenousMesh_32(e_mesh::eMesh{Tri,T2}; color=RGBA{Float32}(0.5, 0.5, 0.5, 1.0)) where {T2}
#     vertices = get_vertices_32(e_mesh)
#     faces = get_faces_32(e_mesh)
#     return HomogenousMesh(vertices=vertices, faces=faces, color=color)
# end
#
# function HomogenousMesh_32(eM::eMesh{Nothing,Tet})
# 	i3 = Vector{SVector{3,Int32}}()
# 	for k = 1:n_tet(eM)
# 		i_tet = eM.tet[k]
# 		ϵ_tet = eM.ϵ[i_tet]
# 		i_tet_new = sort_so_big_ϵ_last(ϵ_tet, i_tet)
# 		push!(i3, i_tet_new[1:3])
# 	end
# 	vertices=get_vertices_32(eM)
# 	faces = [Face{3,Int32}(k) for k = i3]
# 	return HomogenousMesh(vertices=vertices, faces=faces)
# end

function play_recorded_data(mvis::MechanismVisualizer, mech_scen::MechanismScenario, data_time::Vector{Float64},
        data_state::Matrix{Float64}; dt::Float64=1/60, slowdown::Float64=1.0, t0::Float64=-Inf, t1::Float64=Inf)

    ### ERROR checking
    (size(data_state,1) == length(data_time)) || error("the length of the time vector ($(length(data_time))) needs to be the same as the rows of the state matrix ($(size(data_state,1)))")
    issorted(data_time) || error("data_time (vector of times) is not sorted")
    @assert(t0 < t1)
    i0 = findfirst(t0 .<= data_time)
    i1 = findlast(data_time .<= t1)
    delta_t_sim = data_time[i1] - data_time[i0]
    (60.0 < (slowdown * delta_t_sim) ) && error("total video time is greater than 1 minute you probably didn't mean to do this.")

    t_last_frame = data_time[i0] - dt
    for k = i0:i1
        while (t_last_frame + dt) < data_time[k]  # wait until enough real time has passed
            t_last_frame += dt
            sleep(dt * slowdown)
			copyto!(mech_scen.float, data_state[k, :])
			set_configuration!(mech_scen, mvis)
        end
    end
end


text:
##################################
##
## Prices
##
##################################
"""
Type that contains information of prices

## Prices

- `rᴷ`      : rental rate of capital
- `rᴬ`      : return on the illiquid asset
- `rᴮ`      : return on the liquid asset
- `wedge`   :
- `w`       : wage rate

"""
type Prices
    rᴷ ::Float64
    rᴬ ::Float64
    rᴮ ::Float64
    wedge::Float64
    w  ::Float64
end

function update_prices!(pc::Prices, P::Vector{Float64})
    pc.rᴷ, pc.rᴬ, pc.rᴮ, pc.w = P
end

price_partial() = Prices(0.0, 0.04, 0.03, 0.09,4.0)

function Base.show(io::IO, pr::Prices)
    @printf io "\n"
    @printf io "    Prices \n"
    @printf io "      rᴷ    :    %.3f \n" pr.rᴷ
    @printf io "      rᴬ    :    %.3f \n" pr.rᴬ
    @printf io "      rᴮ    :    %.3f \n" pr.rᴮ
    @printf io "      wedge :    %.3f \n" pr.wedge
    @printf io "      wage  :    %.3f \n" pr.w
end

##################################
##
## SteadyState and Transition
##
##################################
"""
Description of the SteadyState
# Endogenous

### Prices
- `rcapital`
- `ra`
- `rb`
- `rborr`
- `wage `

### Interest rate
- `rnom`

### Aggregates to clear
- `Eb`      : aggregate liq saving from households
- `bond`    : total amount of bonds in the market
- `Ea`      : aggregate illiq saving from households
- `capital` : capital demanded by firms
- `labor`   :

### Firm Stats
- `KYratio`
- `KNratio`

- `π`
- `mc`
- `priceadjust`

- `output`
- `profit`
- `dividend`

### Fund
- `divrate`
- `investment`
- `deprec`
- `caputil`

### Government
- `G`
- `lumptransfer`
- `τ`
- `taxrev`
- `govbond`

"""

function initializeSS(KYratio = 8., rb = 0.015)

    ssvar = [:rcapital, :ra, :rb, :wedge, :wage,            # prices
             :rnom,                                         #
             :Eb, :bond, :Ea, :capital, :labor,             # aggregates
             :KYratio, :KNratio, :π, :mc, :priceadjust,     #
             :output, :profit, :dividend,                   #
             :divrate, :investment, :deprec, :caputil,      #
             :G, :lumptransfer, :τ, :taxrev, :govbond,      #
             :tfp]                                          #

    equmSS = (Symbol => Float64)[var => 0.0 for var in ssvar]

    equmSS[:tfp]      = 1.
    equmSS[:KYratio]  = KYratio
    equmSS[:τ]        = 0.25
    equmSS[:rb]       = rb    #having problem with 0.02/4

    equmSS[:wedge]    = 0.10 #0.0161134

    equmSS[:KNratio] = ( equmSS[:tfp] * equmSS[:KYratio] )^(1./(1.-α))
    equmSS[:mc] = (ϵ-1) / ϵ

    #== Prices ==#
    equmSS[:rcapital] = equmSS[:mc] * α / equmSS[:KYratio]
    equmSS[:wage]     = equmSS[:mc] * equmSS[:tfp] * (1. - α) * equmSS[:KNratio]^α

    #== Firm ==#
    equmSS[:priceadjust] = 0
    equmSS[:labor]       =  ( (1 - equmSS[:τ])*equmSS[:wage] /ψ )^σ

    equmSS[:capital] = equmSS[:KNratio] * equmSS[:labor]
    equmSS[:investment] = δbar * equmSS[:capital]

    equmSS[:profit]  = (1. - equmSS[:mc])* equmSS[:capital] / equmSS[:KYratio] - equmSS[:priceadjust]
    equmSS[:divrate] = equmSS[:profit] / equmSS[:capital]

    #== Illiquid return ==#
    equmSS[:ra] = equmSS[:rcapital] - δbar + equmSS[:divrate]

    #== Government ==#
    equmSS[:lumptransfer] = 0.20 * equmSS[:wage] * equmSS[:labor]

    ### Solution ###

    return equmSS

end

_unpackprices(equmSS::Dict{Symbol, Float64}) =
    Prices(equmSS[:rcapital], equmSS[:ra], equmSS[:rb], equmSS[:wedge] , equmSS[:wage])

##################################
## OLD AGGREGATE vars
##################################

"""
Holds the aggregate variables of the econmy

### Prices

- `prices`

### Government
- `Rev::Float64`
- `τ::Float64`
- `T::Float64`

### Firm
- `tfp::Float64`
- `KY::Float64`
- `KN::Float64`

- `output::Float64`
- `profit::Float64`
- `capital::Float64`
- `investment::Float64`

"""
type AggVar

    prices::Prices

    ##Gov
    Rev::Float64
    τ  ::Float64
    T  ::Float64
    G  ::Float64

    ## Firm
    tfp::Float64
    KY ::Float64
    KN ::Float64

    output::Float64
    profit::Float64
    capital::Float64
    investment::Float64

end

function AggVar(;tfp::Float64 = 1., KY::Float64 = 5.,
                τ::Float64 = 0.25, rᴮ::Float64 = 0.03, wedge::Float64 = 0.09)

    KN = (tfp * KY)^(1./(1.-α))

    #== Prices ==#
    rᴷ = mc * α / KY
    w  = mc * tfp * (1. - α) * KN^α
    priceadjust = 0
    ℓ = 1.0
    capital = KN * ℓ
    investment = δbar * capital
    profit  = (1. - mc)* capital / KY - priceadjust
    divr = profit / capital

    rᴬ = rᴷ - δbar + divr

    # Define this so that a fixed fraction receives more transfer than pays
    T = 0.8 * w

    output = 0.
    Rev = 0.0
    G = 0.0

    AggVar(Prices(rᴷ,rᴬ,rᴮ,wedge,w), Rev, τ, T, G, tfp, KY, KN,
    output, profit, capital, investment)
end


text:
using CSV

function filename_to_hilbert_dims(fn)
   digits = parse.(Int, split(splitext(basename(fn))[1], "_"))
   popfirst!(digits), digits
end

function hamilton_example(fn)
    n, ms = filename_to_hilbert_dims(fn)
    cnt = prod(1 << x for x in ms)
    h = zeros(UInt32, cnt)
    X = zeros(UInt8, n, cnt)
    for (i, row) in enumerate(CSV.File(fn))
        h[i] = convert(UInt32, row[1])
        for j in 1:length(ms)
            X[j, i] = convert(UInt8, row[1 + j])
        end
    end
    hsort = sortperm(h)
    h[hsort], X[:, hsort], ms
end

function steps_are_one_away(h, X)
    n = size(X, 1)
    Y = copy(X)
    Y .= zero(UInt8)
    for i in eachindex(h)
        tdiff = 0
        for j in 1:n
            if X[j, i] > Y[j, i]
                tdiff += X[j, i] - Y[j, i]
            else Y[j] > X[j]
                tdiff += Y[j, i] - X[j, i]
            end
        end
        if tdiff > 1
            @show tdiff, i, X[:, i], Y[:, i]
        end
        Y .= X
    end
end


text:
# define FQ
const FQ = Integer

# GMP lib dependency
addmod(q::Integer, a::Integer, b::Integer) = mod(a + b, q)
submod(q::Integer, a::Integer, b::Integer) = mod(a - b, q)
mulmod(q::Integer, _::Function, a::Integer, b::Integer) = mod(a * b, q)
invmod(a::Integer, q::Integer) = Base.invmod(a, q)

sub(a::Integer, b::Integer) = a - b
sub(a::FQ2, b::FQ2) = (sub(a[1], b[1]), sub(a[2], b[2]))
sub(a::FQ3, b::FQ3) = (sub(a[1], b[1]), sub(a[2], b[2]), sub(a[3], b[3]))

add(a::Integer, b::Integer) = a + b
add(a::FQ2, b::FQ2) = (add(a[1], b[1]), add(a[2], b[2]))
add(a::FQ3, b::FQ3) = (add(a[1], b[1]), add(a[2], b[2]), add(a[3], b[3]))

# modular addition function
addmod(q::Integer, a::FQ2, b::FQ2) = (
    addmod(q, a[1], b[1]),
    addmod(q, a[2], b[2]),
)
addmod(q::Integer, a::FQ3, b::FQ3) = (
    addmod(q, a[1], b[1]),
    addmod(q, a[2], b[2]),
    addmod(q, a[3], b[3]),
)

# modular subtraction
submod(q::Integer, a::FQ2, b::FQ2) = (
    submod(q, a[1], b[1]),
    submod(q, a[2], b[2]),
)
submod(q::Integer, a::FQ3, b::FQ3) = (
    submod(q, a[1], b[1]),
    submod(q, a[2], b[2]),
    submod(q, a[3], b[3]),
)

# modular multiplication
sqmod(q::Integer, mulroot::Function, a::Integer) = mulmod(q, mulroot, a, a)

function sqmod(q::Integer, mulroot::Function, a::FQ2)
    # Squaring: Guide to PBC Algo. 5.17
    a0, a1 = a
    v0 = sub(a0, a1)
    v3 = mulroot(a1)
    v3 = sub(a0, v3)
    v2 = mulmod(q, mulroot, a0, a1)
    v0 = add(mulmod(q, mulroot, v0, v3), v2)
    c1 = addmod(q, v2, v2)
    c0 = addmod(q, v0, mulroot(v2))
    return (c0, c1)
end

mulmod(q::Integer, mulroot::Function, a::FQ2, b::Integer) = (
    mulmod(q, mulroot, b, a[1]),
    mulmod(q, mulroot, b, a[2]),
)
mulmod(q::Integer, mulroot::Function, a::Integer, b::FQ2) = (
    mulmod(q, mulroot, a, b[1]),
    mulmod(q, mulroot, a, b[2]),
)
function mulmod(q::Integer, mulroot::Function, a::FQ2, b::FQ2)
    # Multiplication: Guide to PBC Algo. 5.16
    a0, a1 = a
    b0, b1 = b
    v0 = mulmod(q, mulroot, a0, b0)
    v1 = mulmod(q, mulroot, a1, b1)
    c0 = addmod(q, v0, mulroot(v1))
    c1 = mulmod(q, mulroot, add(a0, a1), add(b0, b1))
    c1 = submod(q, c1, add(v0, v1))
    return (c0, c1)
end

mulmod(q::Integer, mulroot::Function, a::FQ3, b::Integer) = (
    mulmod(q, mulroot, b, a[1]),
    mulmod(q, mulroot, b, a[2]),
    mulmod(q, mulroot, b, a[3]),
)
mulmod(q::Integer, mulroot::Function, a::Integer, b::FQ3) = (
    mulmod(q, mulroot, a, b[1]),
    mulmod(q, mulroot, a, b[2]),
    mulmod(q, mulroot, a, b[3]),
)
function sqmod(q::Integer, mulroot::Function, a::FQ3)
    # Squaring: Guide to PBC Algo. 5.22
    a0, a1, a2 = a
    v4 = mulmod(q, mulroot, a0, a1)
    v4 = add(v4, v4)
    v5 = sqmod(q, mulroot, a2)
    c1 = addmod(q, mulroot(v5), v4)
    v2 = sub(v4, v5)
    v3 = sqmod(q, mulroot, a0)
    v4 = add(sub(a0, a1), a2)
    v5 = mulmod(q, mulroot, a1, a2)
    v5 = add(v5, v5)
    v4 = sqmod(q, mulroot, v4)
    c0 = addmod(q, mulroot(v5), v3)
    c2 = addmod(q, v2, add(v4, sub(v5, v3)))
    return (c0, c1, c2)
end

function mulmod(q::Integer, mulroot::Function, a::FQ3, b::FQ3)
    # Multiplication: Guide to PBC Algo. 5.21
    a0, a1, a2 = a
    b0, b1, b2 = b
    v0 = mulmod(q, mulroot, a0, b0)
    v1 = mulmod(q, mulroot, a1, b1)
    v2 = mulmod(q, mulroot, a2, b2)
    c0 = mulmod(q, mulroot, add(a1, a2), add(b1, b2))
    c0 = sub(c0, add(v1, v2))
    c0 = addmod(q, mulroot(c0), v0)
    c1 = mulmod(q, mulroot, add(a0, a1), add(b0, b1))
    c1 = sub(c1, add(v0, v1))
    c1 = addmod(q, c1, mulroot(v2))
    c2 = mulmod(q, mulroot, add(a0, a2), add(b0, b2))
    c2 = addmod(q, sub(c2, add(v0, v2)), v1)
    return (c0, c1, c2)
end

invmod(q::Integer, _::Function, a::Integer) = invmod(a, q)
function invmod(q::Integer, mulroot::Function, a::FQ2)
    # Inversion: Guide to PBC Algo. 5.19
    a0, a1 = a
    v0 = sqmod(q, mulroot, a0)
    v1 = sqmod(q, mulroot, a1)
    v0 = sub(v0, mulroot(v1))
    v1 = invmod(q, mulroot, v0)
    c0 = mulmod(q, mulroot, a0, v1)
    c1 = mulmod(q, mulroot, sub(zero(a1), a1), v1)
    return (c0, c1)
end

function invmod(q::Integer, mulroot::Function, a::FQ3)
    # Inversion: Guide to PBC Algo. 5.23
    a0, a1, a2 = a
    v0 = sqmod(q, mulroot, a0)
    v1 = sqmod(q, mulroot, a1)
    v2 = sqmod(q, mulroot, a2)
    v3 = mulmod(q, mulroot, a0, a1)
    v4 = mulmod(q, mulroot, a0, a2)
    v5 = mulmod(q, mulroot, a1, a2)
    A = sub(v0, mulroot(v5))
    B = sub(mulroot(v2), v3)
    C = sub(v1, v4)
    v6a = mulmod(q, mulroot, a0, A)
    v6b = mulroot(mulmod(q, mulroot, a2, B))
    v6c = mulroot(mulmod(q, mulroot, a1, C))
    v6 = add(v6a, add(v6b, v6c))
    F = invmod(q, mulroot, v6)
    c0 = mulmod(q, mulroot, A, F)
    c1 = mulmod(q, mulroot, B, F)
    c2 = mulmod(q, mulroot, C, F)
    return (c0, c1, c2)
end

# modular exponentiation
# TODO: Optimise
function powmod(q::Integer, mulroot::Function, base, exp::Integer)
    if exp < 2
        if exp < 0
            return invmod(q, mulroot, powmod(q, mulroot, base, -exp))
        end
        return exp == 0 ? one(base) : base
    end
    c = powmod(q, mulroot, sqmod(q, mulroot, base), fld(exp, 2))
    return isodd(exp) ? mulmod(q, mulroot, c, base) : c
end

frobmod(q::Integer, frobmap::Function, mulroot::Function, a::Integer) = a
function frobmod(q::Integer, frobmap::Function, mulroot::Function, a::FQ2)
    g = frobmap(a)[1]
    a1 = frobmod(q, frobmap, mulroot, a[1])
    a2 = frobmod(q, frobmap, mulroot, a[2])
    a = (a1, mulmod(q, mulroot, a2, g))
    return a
end
function frobmod(q::Integer, frobmap::Function, mulroot::Function, a::FQ3)
    g, g2 = frobmap(a)
    a1 = frobmod(q, frobmap, mulroot, a[1])
    a2 = frobmod(q, frobmap, mulroot, a[2])
    a3 = frobmod(q, frobmap, mulroot, a[3])
    a = (a1, mulmod(q, mulroot, a2, g), mulmod(q, mulroot, a3, g2))
    return a
end


text:
#################################################################
# The policy is a mapping from a belief to an action in a POMDP,
# and it maps states to actions in an MDP.
# The policy is extracted through calls to the action() function.
#################################################################

"""
Base type for a policy (a map from every possible belief, or more abstract policy state, to an optimal or suboptimal action)
"""
abstract type Policy end

"""
    action{B}(policy::Policy, x::B)

Fills and returns action based on the current state or belief, given the policy.
`B` is a generalized information state - can be a state in an MDP, a distribution in POMDP,
or any other representation needed to make a decision using the given policy.
"""
function action end

"""
    updater(policy::Policy)

Returns a default Updater appropriate for a belief type that policy `p` can use
"""
function updater end

"""
    value(p::Policy, s)
    value(p::Policy, s, a)

Returns the utility value from policy `p` given the state (or belief), or state-action (or belief-action) pair.

The state-action version is commonly referred to as the Q-value.
"""
function value end


text:
# Test cropping OSM data
module TestCropMap

using OpenStreetMap
using Base.Test
using Compat

import OpenStreetMap: Bounds, LLA, aspectRatio, getX, getY, inBounds

MAP_FILENAME = "tech_square.osm"

bounds = Bounds(42.3637, 42.3655, -71.0919, -71.0893)

nodes, hwys, builds, feats = getOSMData(MAP_FILENAME)
cropMap!(nodes, bounds, highways=hwys, buildings=builds, features=feats, delete_nodes=true)
@test length(nodes) == 198
@test length(hwys) == 16
@test length(builds) == 1
@test length(feats) == 0

# cropHighways removes out-of-bounds points and interpolates points on boundary
let bounds = Bounds(0, 1, 0, 1)
    coords = [
        LLA(-0.1, -0.02),
        LLA(0.1, 0.1),
        LLA(0.2, 0.0),
        LLA(0.3, -0.1),
        LLA(0.4, 0.1),
        LLA(1.0, 0.3),
        LLA(2.0, 0.4)
    ]

    hwy = first(values(hwys))
    hwy.nodes = collect(1:length(coords))
    highways = Compat.@Dict(1 => hwy)

    nodes = Dict{Int,LLA}([(hwy.nodes[i], coords[i]) for i in 1:length(coords)])
    cropMap!(nodes, bounds, highways=highways)

    hwy_nodes = highways[1].nodes

    # same changes made to nodes and highways
    @test isempty(symdiff(keys(nodes), hwy_nodes))

    # everything ends up in bounds
    if VERSION.minor < 4
        @test all(map(x -> inBounds(x, bounds), values(nodes)))
    else
        @test all(x -> inBounds(x, bounds), values(nodes))
    end

    # proper interpolation
    @test getY(nodes[hwy_nodes[1]]) == 0.0
    @test_approx_eq getX(nodes[hwy_nodes[1]]) 0.04

    @test_approx_eq getY(nodes[hwy_nodes[4]]) 0.35
    @test getX(nodes[hwy_nodes[4]]) == 0.0

    @test getY(nodes[hwy_nodes[end]]) == 1.0
    @test_approx_eq getX(nodes[hwy_nodes[end]]) 0.3
end

end # module TestCropMap


text:
module SFrontiers

export sfmodel_spec, sfmodel_init, sfmodel_opt,
       sfmodel_fit, sfmodel_predict, drawTrun,
       genTrun, sf_demean, sfmodel_boot_marginal,
       # likelihood functions
        LL_T,
       # macros for _spec;
        depvar, timevar, idvar,
        sfdist, sftype,
        @sfdist, @sftype,
        @model, @depvar,
        @frontier,
        @μ, @mu,  @σᵤ², @sigma_u_2,  @σᵥ², @sigma_v_2,
        @σₐ², @sigma_a_2,
        @hscale,
        @sfpanel,
        sfpanel,
        @timevar, @idvar,
        @gamma,
        @eq,
       # functions for sfmodel_init
        frontier,
        μ, mu, σᵤ², sigma_u_2, σᵥ², sigma_v_2,
        σₐ², sigma_a_2,
        misc,
        hscale,  gamma,
        all_init,
       # functions for sfmodel_opt
        warmstart_solver, warmstart_maxIT,
        main_solver, main_maxIT, tolerance, verbose, banner,
        ineff_index, marginal, table_format,
       # functions for sfmodel_fit
        useData,
        sfmodel_CI,
       # MoM Test
        sfmodel_MoMTest,
       # functions for JLMS and BC index
        jlmsbc, jlmsbc_marg,
       # the table for regular and mixed Chi-square test
        sfmodel_MixTable, sfmodel_ChiSquareTable,
       # struct
        Sfmodeltype,
        Trun, trun, truncated, t,
        Half, half, h,
        Expo, expo, e,
        Trun_Scale, trun_scale, trun_scaling, s,
        Half_Scale, half_scale, half_scaling,
        MoM,
        production, cost, #* not prod, b/c conflict with Base
        text, html, latex,
        PFEWHH, PFEWHT, PTREH,  PTRET,  PFECSWH, PanDecay, get_marg, PanKumb90,  #* export for testing purpose
        TFE_WH2010, TFE_CSW2014, TFE_G2005, TRE, TimeDecay, Kumbhakar1990,
      # Optim's algorithms
        NelderMead, SimulatedAnnealing, SAMIN, ParticleSwarm,
        ConjugateGradient, GradientDescent, BFGS, LBFGS,
        Newton, NewtonTrustRegion, IPNewton,
      # sfmodel_mtest()
        pickConsNameFromDF

using DataFrames
using DataStructures             # for OrderedDict
using Distributions              # for TDist, Normal
using FLoops                     # multithreading
using ForwardDiff                # for marginal effect
using HypothesisTests            # for pvalue()
using KahanSummation             # for time decay model, true random effect model
using LinearAlgebra              # extract diagnol and Matrix(I,...)
using NLSolversBase              # for hessian!
using Optim
using PrettyTables               # making tables
using QuadGK                     # for TFE_CSW2014 model
using Random                     # for sfmodel_boot_marginal
using RowEchelon                 # for checkCollinear, check multi-collinearity
using SpecialFunctions           # for erfi used in sfmodel_MoMTest
using StatsFuns                  # for normlogpdf(), normlogcdf()
using Statistics                 #

#############################
##   Define Model Types    ##
#############################

abstract type Sfmodeltype end
  struct Trun      <: Sfmodeltype end
  struct truncated <: Sfmodeltype end
  struct trun      <: Sfmodeltype end
  struct t         <: Sfmodeltype end
  struct Half      <: Sfmodeltype end
  struct half      <: Sfmodeltype end
  struct h         <: Sfmodeltype end
  struct Expo      <: Sfmodeltype end
  struct expo      <: Sfmodeltype end
  struct e <: Sfmodeltype end
  struct Trun_Scale   <: Sfmodeltype end
  struct trun_scale   <: Sfmodeltype end
  struct trun_scaling <: Sfmodeltype end
  struct s <: Sfmodeltype end
  struct Half_Scale   <: Sfmodeltype end
  struct half_scale   <: Sfmodeltype end
  struct half_scaling <: Sfmodeltype end
  struct MoM          <: Sfmodeltype end

  struct PFEWHT   <: Sfmodeltype end # panel fixed-effet of Wang and Ho 2010, truncated normal
  struct PFEWHH   <: Sfmodeltype end # panel fixed-effet of Wang and Ho 2010, half normal
  struct PTREH    <: Sfmodeltype end # panel true random effect model, half normal
  struct PTRET    <: Sfmodeltype end # panel true random effect model, truncated normal
  struct PFECSWH <: Sfmodeltype end
  struct PanDecay <: Sfmodeltype end # Panel Time Decay Model of Battese and Coelli (1992)
  struct PanKumb90   <: Sfmodeltype end # Kumbhakar (1990)

abstract type PorC end
  struct production <: PorC end
  struct cost <: PorC end

abstract type PanelModel end
  struct TFE_WH2010  <: PanelModel end
  struct TFE_CSW2014 <: PanelModel end
  struct TFE_G2005   <: PanelModel end
  struct TRE         <: PanelModel end
  struct TimeDecay   <: PanelModel end
  struct Kumbhakar1990 <: PanelModel end

abstract type TableFormat end
  struct text   <: TableFormat end
  struct html   <: TableFormat end
  struct latex  <: TableFormat end

################################################
##    include other files; order important    ##
################################################

include("SFmacfun.jl")
include("SFloglikefun.jl")
include("SFutil.jl")
include("SFmtest.jl")
include("SFgetvars.jl")
include("SFindex.jl")
include("SFpredict.jl")
include("SFmarginal.jl")
include("SFmainfun.jl")

end # SFrontiers module


text:
using Decimals
using Test

@testset "Arithmetic" begin

@testset "Addition" begin
    @test Decimal(0.1) + 0.2 == 0.1 + Decimal(0.2) == Decimal(0.1) + Decimal(0.2) == Decimal(0.3)
    @test Decimal.([0.1 0.2]) .+ [0.3 0.1] == Decimal.([0.4 0.3])
    @test Decimal(2147483646) + Decimal(1) == Decimal(2147483647)
    @test Decimal(1,3,-2) + parse(Decimal, "0.2523410412138103") == Decimal(0,2223410412138103,-16)
end

@testset "Subtraction" begin
    @test Decimal(0.3) - 0.1 == 0.3 - Decimal(0.1)
    @test 0.3 - Decimal(0.1) == Decimal(0.3) - Decimal(0.1)
    @test Decimal(0.3) - Decimal(0.1) == Decimal(0.2)
    @test Decimal.([0.3 0.1]) .- [0.1 0.5] == Decimal.([0.2 -0.4])
end

@testset "Negation" begin
    @test -Decimal.([0.3 0.2]) == [-Decimal(0.3) -Decimal(0.2)]
    @test -Decimal(0.3) == zero(Decimal) - Decimal(0.3)
    @test iszero(decimal(12.1) - decimal(12.1))
end

@testset "Multiplication" begin
    @test Decimal(12.21) * Decimal(2.12) == Decimal(0,258852,-4)
    @test Decimal(12.2112543) * Decimal(2.121352) == Decimal(0,259043687318136,-13)
    @test Decimal(0.2) * 0.1 == 0.2 * Decimal(0.1)
    @test 0.2 * Decimal(0.1) == Decimal(0.02)
    @test Decimal(12.34) * 0.1234 == 12.34 * Decimal(0.1234)
    @test 12.34 * Decimal(0.1234) == Decimal(1.522756)
    @test Decimal(0.21084210) * -2 == -2 * Decimal(0.21084210)
    @test -2 * Decimal(0.21084210) == Decimal(-0.4216842)
    @test Decimal(0, 2, -1) * 0.0 == zero(Decimal)
    @test Decimal.([0.3, 0.6]) .* 5 == [Decimal(0.3)*5, Decimal(0.6)*5]
    @test one(Decimal) * 1 == Decimal(0, 1, 0)
end

@testset "Inversion" begin
    @test inv(Decimal(0, 1, -1)) == Decimal(0, 1, 1)
    @test inv(Decimal(0, 1, 1)) == Decimal(0, 1, -1)
    @test inv(Decimal(1, 2, -1)) == Decimal(1, 5, 0)
    @test inv(Decimal(1, 5, 0)) == Decimal(1, 2, -1)
    @test inv(Decimal(0, 2, -2)) == Decimal(0, 5, 1)
    @test inv(Decimal(0, 5, 1)) == Decimal(0, 2, -2)
    @test inv(Decimal(1, 4, -1)) == Decimal(1, 25, -1)
    @test inv(Decimal(1, 25, -1)) == Decimal(1, 4, -1)
end

@testset "Division" begin
    @test Decimal(0.2) / Decimal(0.1) == Decimal(2)
    @test Decimal(0.3) / Decimal(0.1) == Decimal(0,3,0)
    @test [Decimal(0.3) / Decimal(0.1), Decimal(0.6) / Decimal(0.1)] == [Decimal(0.3), Decimal(0.6)] ./ Decimal(0.1)
    @test [Decimal(0.3) / 0.1, Decimal(0.6) / 0.1] == [Decimal(0.3), Decimal(0.6)] ./ 0.1
end

end


text:
"""
```julia
vmean(A; dims, multithreaded=:auto)
```
Compute the mean of all elements in `A`, optionally over dimensions specified by `dims`.
As `Statistics.mean`, but vectorized and (optionally) multithreaded.

## Examples
```julia
julia> using VectorizedStatistics

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> vmean(A, dims=1)
1×2 Matrix{Float64}:
 2.0  3.0

julia> vmean(A, dims=2)
2×1 Matrix{Float64}:
 1.5
 3.5
```
"""
function vmean(A; dims=:, multithreaded=:auto)
    if (multithreaded===:auto && length(A) > 4095) || multithreaded===true
        _vtmean(A, dims)
    else
        _vmean(A, dims)
    end
end
vtmean(A; dims=:) = _vtmean(A, dims)
export vmean

## Singlethreaded implementation
# Reduce one dim
_vmean(A, dims::Int) = _vmean(A, (dims,))

# Reduce some dims
function _vmean(A::AbstractArray{T,N}, dims::Tuple) where {T,N}
    sᵢ = size(A)
    sₒ = ntuple(Val(N)) do d
        ifelse(d ∈ dims, 1, sᵢ[d])
    end
    Tₒ = Base.promote_op(/, T, Int)
    B = similar(A, Tₒ, sₒ)
    _vmean!(B, A, dims)
end

# Reduce all the dims!
function _vmean(A, ::Colon)
    # Promote type of accumulator to avoid overflow
    Tₒ = Base.promote_op(/, eltype(A), Int)
    Σ = zero(Tₒ)
    @turbo for i ∈ eachindex(A)
        Σ += A[i]
    end
    return Σ / length(A)
end

# Chris Elrod metaprogramming magic:
# Generate customized set of loops for a given ndims and a vector
# `static_dims` of dimensions to reduce over
function staticdim_mean_quote(static_dims::Vector{Int}, N::Int)
  M = length(static_dims)
  # `static_dims` now contains every dim we're taking the mean over.
  Bᵥ = Expr(:call, :view, :B)
  reduct_inds = Int[]
  nonreduct_inds = Int[]
  # Firstly, build our expressions for indexing each array
  Aind = :(A[])
  Bind = :(Bᵥ[])
  inds = Vector{Symbol}(undef, N)
  len = Expr(:call, :*)
  for n ∈ 1:N
    ind = Symbol(:i_,n)
    inds[n] = ind
    push!(Aind.args, ind)
    if n ∈ static_dims
      push!(reduct_inds, n)
      push!(Bᵥ.args, :(firstindex(B,$n)))
      push!(len.args, :(size(A, $n)))
    else
      push!(nonreduct_inds, n)
      push!(Bᵥ.args, :)
      push!(Bind.args, ind)
    end
  end
  # Secondly, build up our set of loops
  if !isempty(nonreduct_inds)
    firstn = first(nonreduct_inds)
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = indices((A,B),$firstn)), block)
    if length(nonreduct_inds) > 1
      for n ∈ @view(nonreduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = indices((A,B),$n)), newblock))
        block = newblock
      end
    end
    rblock = block
    # Push more things here if you want them at the beginning of the reduction loop
    push!(rblock.args, :(Σ = zero(eltype(Bᵥ))))
    # Build the reduction loop
    for n ∈ reduct_inds
      newblock = Expr(:block)
      push!(block.args, Expr(:for, :($(inds[n]) = axes(A,$n)), newblock))
      block = newblock
    end
    # Push more things here if you want them in the innermost loop
    push!(block.args, :(Σ += $Aind))
    # Push more things here if you want them at the end of the reduction loop
    push!(rblock.args, :($Bind = Σ * invdenom))
    # Put it all together
    return quote
      invdenom = inv($len)
      Bᵥ = $Bᵥ
      @turbo $loops
      return B
    end
  else
    firstn = first(reduct_inds)
    # Build the reduction loop
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = axes(A,$firstn)), block)
    if length(reduct_inds) > 1
      for n ∈ @view(reduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = axes(A,$n)), newblock))
        block = newblock
      end
    end
    # Push more things here if you want them in the innermost loop
    push!(block.args, :(Σ += $Aind))
    # Put it all together
    return quote
      invdenom = inv($len)
      Bᵥ = $Bᵥ
      Σ = zero(eltype(Bᵥ))
      @turbo $loops
      Bᵥ[] = Σ * invdenom
      return B
    end
  end
end

# Chris Elrod metaprogramming magic:
# Turn non-static integers in `dims` tuple into `StaticInt`s
# so we can construct `static_dims` vector within @generated code
function branches_mean_quote(N::Int, M::Int, D)
  static_dims = Int[]
  for m ∈ 1:M
    param = D.parameters[m]
    if param <: StaticInt
      new_dim = _dim(param)::Int
      @assert new_dim ∉ static_dims
      push!(static_dims, new_dim)
    else
      t = Expr(:tuple)
      for n ∈ static_dims
        push!(t.args, :(StaticInt{$n}()))
      end
      q = Expr(:block, :(dimm = dims[$m]))
      qold = q
      ifsym = :if
      for n ∈ 1:N
        n ∈ static_dims && continue
        tc = copy(t)
        push!(tc.args, :(StaticInt{$n}()))
        qnew = Expr(ifsym, :(dimm == $n), :(return _vmean!(B, A, $tc)))
        for r ∈ m+1:M
          push!(tc.args, :(dims[$r]))
        end
        push!(qold.args, qnew)
        qold = qnew
        ifsym = :elseif
      end
      # Else, if dimm ∉ 1:N, drop it from list and continue
      tc = copy(t)
      for r ∈ m+1:M
        push!(tc.args, :(dims[$r]))
      end
      push!(qold.args, Expr(:block, :(return _vmean!(B, A, $tc))))
      return q
    end
  end
  return staticdim_mean_quote(static_dims, N)
end

# Efficient @generated in-place mean
@generated function _vmean!(B::AbstractArray{Tₒ,N}, A::AbstractArray{T,N}, dims::D) where {Tₒ,T,N,M,D<:Tuple{Vararg{Integer,M}}}
  branches_mean_quote(N, M, D)
end
@generated function _vmean!(B::AbstractArray{Tₒ,N}, A::AbstractArray{T,N}, dims::Tuple{}) where {Tₒ,T,N}
  :(copyto!(B, A); return B)
end

## As above, but multithreaded
# Reduce one dim
_vtmean(A, dims::Int) = _vtmean(A, (dims,))

# Reduce some dims
function _vtmean(A::AbstractArray{T,N}, dims::Tuple) where {T,N}
    sᵢ = size(A)
    sₒ = ntuple(Val(N)) do d
        ifelse(d ∈ dims, 1, sᵢ[d])
    end
    Tₒ = Base.promote_op(/, T, Int)
    B = similar(A, Tₒ, sₒ)
    _vtmean!(B, A, dims)
end

# Reduce all the dims!
function _vtmean(A, ::Colon)
    # Promote type of accumulator to avoid overflow
    Tₒ = Base.promote_op(/, eltype(A), Int)
    Σ = zero(Tₒ)
    @tturbo for i ∈ eachindex(A)
        Σ += A[i]
    end
    return Σ / length(A)
end

# Chris Elrod metaprogramming magic:
# Generate customized set of loops for a given ndims and a vector
# `static_dims` of dimensions to reduce over
function staticdim_tmean_quote(static_dims::Vector{Int}, N::Int)
  M = length(static_dims)
  # `static_dims` now contains every dim we're taking the mean over.
  Bᵥ = Expr(:call, :view, :B)
  reduct_inds = Int[]
  nonreduct_inds = Int[]
  # Firstly, build our expressions for indexing each array
  Aind = :(A[])
  Bind = :(Bᵥ[])
  inds = Vector{Symbol}(undef, N)
  len = Expr(:call, :*)
  for n ∈ 1:N
    ind = Symbol(:i_,n)
    inds[n] = ind
    push!(Aind.args, ind)
    if n ∈ static_dims
      push!(reduct_inds, n)
      push!(Bᵥ.args, :(firstindex(B,$n)))
      push!(len.args, :(size(A, $n)))
    else
      push!(nonreduct_inds, n)
      push!(Bᵥ.args, :)
      push!(Bind.args, ind)
    end
  end
  # Secondly, build up our set of loops
  if !isempty(nonreduct_inds)
    firstn = first(nonreduct_inds)
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = indices((A,B),$firstn)), block)
    if length(nonreduct_inds) > 1
      for n ∈ @view(nonreduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = indices((A,B),$n)), newblock))
        block = newblock
      end
    end
    rblock = block
    # Push more things here if you want them at the beginning of the reduction loop
    push!(rblock.args, :(Σ = zero(eltype(Bᵥ))))
    # Build the reduction loop
    for n ∈ reduct_inds
      newblock = Expr(:block)
      push!(block.args, Expr(:for, :($(inds[n]) = axes(A,$n)), newblock))
      block = newblock
    end
    # Push more things here if you want them in the innermost loop
    push!(block.args, :(Σ += $Aind))
    # Push more things here if you want them at the end of the reduction loop
    push!(rblock.args, :($Bind = Σ * invdenom))
    # Put it all together
    return quote
      invdenom = inv($len)
      Bᵥ = $Bᵥ
      @tturbo $loops
      return B
    end
  else
    firstn = first(reduct_inds)
    # Build the reduction loop
    block = Expr(:block)
    loops = Expr(:for, :($(inds[firstn]) = axes(A,$firstn)), block)
    if length(reduct_inds) > 1
      for n ∈ @view(reduct_inds[2:end])
        newblock = Expr(:block)
        push!(block.args, Expr(:for, :($(inds[n]) = axes(A,$n)), newblock))
        block = newblock
      end
    end
    # Push more things here if you want them in the innermost loop
    push!(block.args, :(Σ += $Aind))
    # Put it all together
    return quote
      invdenom = inv($len)
      Bᵥ = $Bᵥ
      Σ = zero(eltype(Bᵥ))
      @tturbo $loops
      Bᵥ[] = Σ * invdenom
      return B
    end
  end
end
# Chris Elrod metaprogramming magic:
# Turn non-static integers in `dims` tuple into `StaticInt`s
# so we can construct `static_dims` vector within @generated code
function branches_tmean_quote(N::Int, M::Int, D)
  static_dims = Int[]
  for m ∈ 1:M
    param = D.parameters[m]
    if param <: StaticInt
      new_dim = _dim(param)::Int
      @assert new_dim ∉ static_dims
      push!(static_dims, new_dim)
    else
      t = Expr(:tuple)
      for n ∈ static_dims
        push!(t.args, :(StaticInt{$n}()))
      end
      q = Expr(:block, :(dimm = dims[$m]))
      qold = q
      ifsym = :if
      for n ∈ 1:N
        n ∈ static_dims && continue
        tc = copy(t)
        push!(tc.args, :(StaticInt{$n}()))
        qnew = Expr(ifsym, :(dimm == $n), :(return _vtmean!(B, A, $tc)))
        for r ∈ m+1:M
          push!(tc.args, :(dims[$r]))
        end
        push!(qold.args, qnew)
        qold = qnew
        ifsym = :elseif
      end
      # Else, if dimm ∉ 1:ndims(A), drop it from list and continue
      tc = copy(t)
      for r ∈ m+1:M
        push!(tc.args, :(dims[$r]))
      end
      push!(qold.args, Expr(:block, :(return _vtmean!(B, A, $tc))))
      return q
    end
  end
  return staticdim_tmean_quote(static_dims, N)
end

# Efficient @generated in-place mean
@generated function _vtmean!(B::AbstractArray{Tₒ,N}, A::AbstractArray{T,N}, dims::D) where {Tₒ,T,N,M,D<:Tuple{Vararg{Integer,M}}}
  branches_tmean_quote(N, M, D)
end
@generated function _vtmean!(B::AbstractArray{Tₒ,N}, A::AbstractArray{T,N}, dims::Tuple{}) where {Tₒ,T,N}
  :(copyto!(B, A); return B)
end
##


text:
#=
Module with functions for matrices in finite fields.

Author: Ivan A. Moreno Soto
Last updated: 2019/July/15
=#

module FiniteFieldMatrix

using Random
using Nemo
using LinearAlgebra

include("./ExtensionField.jl")
using .ExtensionField

export
    multiplyPolyMatrices!,
    gaussianEliminationColumnPivoting!,
    GeneratorInverse!,
    getPermutation,
    getScrambler

"""
    multiplyPolyMatrices!(A, B, g)

Multiplies two matrices with coefficients in GF(2^m).
"""
function multiplyPolyMatrices!(A, B, g)
    C = zeros(UInt32, (size(A)[1], size(B)[2]))

    for i in 1:size(A)[1]
        for j in 1:size(A)[2]
            for k in 1:size(B)[2]
                temp = multiply(A[i, j], B[j, k])
                C[i, k] = add(C[i, k], temp)
            end
        end
    end

    # Applying modulo to every entry.
    C = modulo.(C, g)

    return C
end

"""
    searchPivot!(H, i, n)

Searches for the ith pivot in H in the n-i columns left.
"""
function searchPivot!(H, i, n)
    pivot = i
    for j in (i+1):n
        if H[i, j] == 1
            pivot = j
            break
        end
    end

    return pivot
end

"""
    permuteColumns!(H, i, pivot, k)

Permutes columns i and pivot with k entries each one in H.
"""
function permuteColumns!(H, i, pivot, k)
    for j in 1:k
        temp = H[j, i]
        H[j, i] = H[j, pivot]
        H[j, pivot] = temp
    end
end

"""
    pivotColumns!(H, i, n)

Pivots column i of matrix H. The function can search up to n-i columns
right.
"""
function pivotColumns!(H, i, n)
    k = size(H)[1]
    pivot = searchPivot!(H, i, n)

    if i == pivot
        error("Pivoting failed.")
    end

    permuteColumns!(H, i, pivot, k)

    return pivot
end

"""
    swapCodeSupport!(codeSupport, i, j)

Swaps two elements of a code support.
"""
function swapCodeSupport!(codeSupport, i, j)
    t = codeSupport[i]
    codeSupport[i] = codeSupport[j]
    codeSupport[j] = t
end

"""
    gaussianEliminationColumnPivoting!(H, n, codeSupport)

Applies Gauss-Jordan elimination to matrix H of dimension nxk. It also permutes
a given code support to keep a coherent goppa code.
"""
function gaussianEliminationColumnPivoting!(H, n, codeSupport)
    k = size(H)[2]
    # In GF(2)
    for i in 1:n
        # Pivoting the matrix to ensure that it can be reduced.
        if H[i, i] != 1
            pivot = pivotColumns!(H, i, k)
            swapCodeSupport!(codeSupport, i, pivot)
        end

        if H[i, i] != 1
            error("The matrix is not full rank.")
        end

        for j in (i+1):n
            if H[j, i] != 0
                for l in i:k
                    H[j, l] = (H[j, l] + H[i, l]) % 2
                end
            end
        end
    end

    for i in 2:n
        for j in 1:(i-1)
            if H[j, i] != 0
                for l in i:k
                    H[j, l] = (H[j, l] + H[i, l]) % 2
                end
            end
        end
    end
end

"""
    GeneratorInverse!(G)

Finds the right inverse of a generator matrix.
"""
function GeneratorInverse!(G)
    A = (G * transpose(G)) .% 2
    AInverse = BinaryInverse!(A)
    return (transpose(G) * AInverse) .% 2
end

"""
    BinaryInverse!(A)

Finds the inverse of a binary matrix A. Can throw an exception
if the inverse doesn't exist.
"""
function BinaryInverse!(A)
    B = Array{Int}(A)
    n = size(B)[1]
    R = ResidueRing(ZZ, 2)
    S = MatrixSpace(R, n, n)
    NemoInverse = inv(S(B))
    Inverse = zeros(UInt32, (n, n))
    for i in 1:n
        for j in 1:n
            if NemoInverse[i, j] == 1
                Inverse[i, j] = 1
            end
        end
    end
    return Inverse
end

"""
    getPermutation(n)

Returns a random permutation matrix of dimension nxn.
"""
function getPermutation(n)
    P = Matrix{UInt32}(I, n, n)
    return P[randperm(n), :]
end

"""
    fillBinaryMatrixRandomly!(M)

Fills randomly the entries of a matrix M with 0s and 1s.
"""
function fillBinaryMatrixRandomly!(M, k)
    for i in 1:k
        M[i, i] = 1
    end

    for i in 1:k
        for j in i:k
            if rand() <= 0.5
                M[i, j] = 1
            end
        end
    end

    for i in 2:k
        for j in 1:k
            M[i, j] = (M[i, j] + M[1, j]) % 2
        end
    end
end

"""
    getScrambler(k)

Returns a random binary non-singular matrix of dimension kxk.
"""
function getScrambler(k)
    S = SInverse = zeros(UInt32, (k, k))
    fillBinaryMatrixRandomly!(S, k)
    SInverse = BinaryInverse!(S)

    return S, SInverse
end

end # FiniteFieldMatrix

max_stars_repo_name: rexemin/McEliece

max_issues_repo_name: rexemin/McEliece

max_forks_repo_name: rexemin/McEliece


text:
using DelimitedFiles

function read_numbers(input_filename)
    number_string = readlines(input_filename)[1]
    number_strings = split(number_string, ",")
    numbers = map(x->parse(Int64, x), number_strings)
    return numbers
end

function read_bingo_cards(input_filename)
    bigmatrix = readdlm(input_filename, Int, skipstart=1)
    cards = []
    # Assume cards are square
    card_size = size(bigmatrix)[2]
    for i in 1:Int64((size(bigmatrix)[1]/card_size))
        push!(cards, bigmatrix[(i-1)*card_size+1:i*card_size, :])
    end
    return cards
end

function check_row_or_column(row_or_column, nums)
    return sum(in(nums).(row_or_column)) == length(row_or_column)
end

function check_card_win(card, nums)
    card_size = size(card)[1]
    for i in 1:card_size
        if (check_row_or_column(card[i,:],nums) ||
            check_row_or_column(card[:,i],nums))
            return true
        end
    end
    return false
end

function test_card_win()
    card = [1 2 3
            4 5 6
            7 8 9]
    @assert check_card_win(card, [1,2]) == false
    @assert check_card_win(card, [1,2,3]) == true
    @assert check_card_win(card, [1,4,7]) == true
    @assert check_card_win(card, [9,5,8,1,6,4]) == true
end

function get_card_score(card, nums)
    last_number_called = nums[length(nums)]
    # flatten the card
    card_flat = vec(card)
    # get a bit vector of numbers that are not called
    matches = (!).(in(nums).(card_flat))
    # multiply by the numbers in the card
    total = sum(matches .* card_flat)
    return total*last_number_called
end

function test_card_score()
    card = [1 2 3
            4 5 6
            7 8 9]
    # should be (1+2+3+4)*9 = 90
    @assert get_card_score(card, [5,6,7,8,9]) == 90
    # should be (2+7)*4 = 36
    @assert get_card_score(card, [9,8,1,3,6,23,5,4]) == 36
end

test_card_win()
test_card_score()

numbers = read_numbers("input.txt")
cards = read_bingo_cards("input.txt")

num_cards = length(cards)
card_size = size(cards[1])[1]

have_winning_card = false
for i in card_size:length(numbers)
    local nums_so_far = numbers[1:i]
    for icard in 1:num_cards
        if check_card_win(cards[icard], nums_so_far)
            global have_winning_card = true
            println("Card ",icard, " wins ", nums_so_far)
            println("Winning score ",get_card_score(cards[icard], nums_so_far))
            break
        end
    end
    if have_winning_card
        break
    end
end


text:
using LinearAlgebra
using Statistics: mean, std

Base.signbit(z::Complex{T}) where T<:Number = signbit(real(z))

"""
    integrate is the main entry point

    input:
    ------
    eq: a Symbolics expression to integrate
    x: the independent variable (optional)

    abstol: the desired tolerance
    num_steps: the number of different steps with expanding basis to be tried
    num_trials: the number of trials in each step (no changes to the basis)
    radius: the radius of the disk in the complex plane to generate random test points
    show_basis: if true, the basis (list of candidate terms) is printed
    opt: the sparse regression optimizer
    bypass: if true, do not integrate terms separately but consider all at once
    symbolic: try symbolic integration first
    max_basis: the maximum number of candidate terms to consider
    verbose: print a detailed report
    complex_plane: generate random test points on the complex plane (if false, the points will be on real axis)
    homotomy: use the homotopy algorithm to generat the basis

    output:
    -------
    solved, unsolved, err

    solved is the solved integral and unsolved is the residual unsolved portion of the input
    err is the numerical error in reaching the solution
"""
function integrate(eq, x=nothing; abstol=1e-6, num_steps=2, num_trials=5, radius=1.0,
                   show_basis=false, opt = STLSQ(exp.(-10:1:0)), bypass=false,
                   symbolic=true, max_basis=100, verbose=false, complex_plane=true,
                   homotopy=true)
    eq = expand(eq)
    eq = apply_div_rule(eq)

    if x == nothing
        x = var(eq)
        if x == nothing
            @syms 𝑥
            x = 𝑥
        end
    end

    l = Logger(verbose)

    # eq is a constant
    if !isdependent(eq, x)
        return x * eq, 0, 0
    end

    # homotopy = homotopy && !bypass

    return integrate_sum(eq, x, l; bypass, abstol, num_trials, num_steps,
                              radius, show_basis, opt, symbolic,
                              max_basis, verbose, complex_plane, homotopy)
end

"""
    integrate_sum applies the integral summation rule ∫ Σᵢ fᵢ(x) dx = Σᵢ ∫ fᵢ(x) dx

    eq: the integrand
    x: the indepedent variable
    l: a logger

    output is the same as integrate
"""
function integrate_sum(eq, x, l; bypass=false, kwargs...)
    solved = 0
    unsolved = 0
    ϵ₀ = 0
    ts = bypass ? [eq] : terms(eq)

    if length(ts) > 1
        inform(l, "Integrating sum", ts)
    end

    for p in ts
        s, u, ϵ = integrate_term(p, x, l; kwargs...)
        solved += s
        unsolved += u
        ϵ₀ = max(ϵ₀, ϵ)
    end

    if !isequal(unsolved, 0)
        eq = factor_rational(simplify_trigs(unsolved))

        if !isequal(eq, unsolved)
            eq = expand(eq)
            unsolved = 0
            ϵ₀ = 0
            ts = bypass ? [eq] : terms(eq)

            if length(ts) > 1
                inform(l, "Integrating transformed sum", ts)
            else
                inform(l, "Transforming the expression", ts[1])
            end

            for p in ts
                s, u, ϵ = integrate_term(p, x, l; kwargs...)
                solved += s
                unsolved += u
                ϵ₀ = max(ϵ₀, ϵ)

                if !isequal(u, 0)   # premature termination on the first failure
                    return 0, eq, ϵ₀
                end

            end
        end
    end

    return expand(solved), unsolved, ϵ₀
end

"""
    integrate_term is the central part of the code that tries different
    methods to integrate eq

    eq: the integrand
    x: the indepedent variable
    l: a logger

    output is the same as integrate
"""
function integrate_term(eq, x, l; kwargs...)
    args = Dict(kwargs)
    abstol, num_steps, num_trials, show_basis, symbolic, verbose, max_basis,
        radius, homotopy = args[:abstol], args[:num_steps],
        args[:num_trials], args[:show_basis], args[:symbolic], args[:verbose],
        args[:max_basis], args[:radius], args[:homotopy]

    attempt(l, "Integrating term", eq)

    if is_number(eq)
        y = eq * x
        result(l, "Successful", y)
        return y, 0, 0
    end

    basis = generate_basis(eq, x; homotopy)

    if show_basis
        inform(l, "Generating basis (|β| = $(length(basis)))", basis)
    end

    if length(basis) > max_basis
        result(l, "|β| = $(length(basis)) is too large")
        return 0, eq, Inf
    end

    D = Differential(x)
    ϵ₀ = Inf
    y₀ = 0

    # rescue
    ϵᵣ = Inf
    yᵣ = 0

    for i = 1:num_steps
        if length(basis) > max_basis break end

        Δbasis = [expand_derivatives(D(f)) for f in basis]

        if symbolic
            y, ϵ = try_symbolic(Float64, eq, x, basis, Δbasis; kwargs...)

            if !isequal(y, 0) && accept_solution(eq, x, y, radius) < abstol
                result(l, "Successful symbolic", y)
                return y, 0, 0
            else
                inform(l, "Failed symbolic")
            end
        end

        for j = 1:num_trials
            r = radius #*sqrt(2)^j
            y, ϵ = try_integrate(Float64, eq, x, basis, Δbasis, r; kwargs...)

            ϵ = accept_solution(eq, x, y, r)
            if ϵ < abstol
                result(l, "Successful numeric (attempt $j out of $num_trials)", y)
                return y, 0, ϵ
            elseif ϵ < ϵᵣ
                ϵᵣ = ϵ
                yᵣ = y
            end
        end

        inform(l, "Failed numeric")

        if i < num_steps
            basis = expand_basis(basis, x)
            if show_basis
                inform(l, "Expanding the basis (|β| = $(length(basis)))", basis)
            elseif verbose
                inform(l, "Expanding the basis (|β| = $(length(basis)))")
            end
        end
    end

    if ϵᵣ < abstol * 10
        result(l, "Accepting numeric (rescued)", yᵣ)
        return yᵣ, 0, ϵᵣ
    else
        result(l, "Unsucessful", eq)
        return 0, eq, ϵ₀
    end
end

###############################################################################

"""
    the core of the randomized parameter-fitting algorithm

    `try_integrate` tries to find a linear combination of the basis, whose
    derivative is equal to eq

    output
    -------
    integral, error
"""
function try_integrate(T, eq, x, basis, Δbasis, radius; kwargs...)
    args = Dict(kwargs)
    abstol, opt, complex_plane, verbose =
        args[:abstol], args[:opt], args[:complex_plane], args[:verbose]

    basis = basis[2:end]    # remove 1 from the beginning
    Δbasis = Δbasis[2:end]
    n = length(basis)

    # A is an nxn matrix holding the values of the fragments at n random points
    A = zeros(Complex{T}, (n, n))
    X = zeros(Complex{T}, n)

    init_basis_matrix!(T, A, X, x, eq, Δbasis, radius, complex_plane; abstol)

    y₁, ϵ₁ = sparse_fit(T, A, x, basis, Δbasis, opt; abstol)
    if ϵ₁ < abstol
        return y₁, ϵ₁
    end

    y₂, ϵ₂ = find_singlet(T, A, basis; abstol)
    if ϵ₂ < abstol
        return y₂, ϵ₂
    end

    if n < 8    # here, 8 is arbitrary and signifies a small basis
        y₃, ϵ₃ = find_dense(T, A, basis; abstol)
        if ϵ₃ < abstol
            return y₃, ϵ₃
        end
    end

    # moving toward the poles
    ∂eq = expand_derivatives(Differential(x)(eq))
    modify_basis_matrix!(T, A, X, x, eq, ∂eq, Δbasis, radius; abstol)
    y₄, ϵ₄ = sparse_fit(T, A, x, basis, Δbasis, opt; abstol)

    if ϵ₄ < abstol || ϵ₄ < ϵ₁
        return y₄, ϵ₄
    else
        return y₁, ϵ₁
    end
end

function init_basis_matrix!(T, A, X, x, eq, Δbasis, radius, complex_plane; abstol=1e-6)
    n = size(A, 1)
    k = 1
    i = 1

    while k <= n
        try
            x₀ = test_point(complex_plane, radius)
            X[k] = x₀
            d = Dict(x => x₀)

            b₀ = Complex{T}(substitute(eq, d))
            if is_proper(b₀)
                for j = 1:n
                    A[k,j] = Complex{T}(substitute(Δbasis[j], d)) / b₀
                end
                if all(is_proper, A[k,:])
                    k += 1
                end
            end
        catch e
            println("Error from init_basis_matrix!: ", e)
        end
    end
end

function modify_basis_matrix!(T, A, X, x, eq, ∂eq, Δbasis, radius; abstol=1e-6)
    n = size(A, 1)
    k = 1
    for k = 1:n
        d = Dict(x => X[k])
        # One Newton iteration toward the poles
        # note the + sign instead of the usual - in Newton-Raphson's method. This is
        # because we are moving toward the poles and not zeros.
        x₀ = X[k] + Complex{T}(substitute(eq, d)) / Complex{T}(substitute(∂eq, d))
        X[k] = x₀
        d = Dict(x => x₀)
        b₀ = Complex{T}(substitute(eq, d))
        for j = 1:n
            A[k,j] = Complex{T}(substitute(Δbasis[j], d)) / b₀
        end
    end
end

function sparse_fit(T, A, x, basis, Δbasis, opt; abstol=1e-6)
    n = length(basis)
    # find a linearly independent subset of the basis
    l = find_independent_subset(A; abstol)
    A, basis, Δbasis, n = A[l,l], basis[l], Δbasis[l], sum(l)

    try
        b = ones(n)
        # q₀ = A \ b
        q₀ = Optimize.init(opt, A, b)
        @views Optimize.sparse_regression!(q₀, A, permutedims(b)', opt, maxiter = 1000)
        ϵ = rms(A * q₀ - b)
        q = nice_parameter.(q₀)
        if sum(iscomplex.(q)) > 2 return nothing, Inf end   # eliminating complex coefficients
        return sum(q[i]*basis[i] for i = 1:length(basis) if q[i] != 0; init=zero(x)), abs(ϵ)
    catch e
        println("Error from sparse_fit", e)
        return nothing, Inf
    end
end

function find_singlet(T, A, basis; abstol)
    σ = vec(std(A; dims=1))
    μ = vec(mean(A; dims=1))
    l = (σ .< abstol) .* (abs.(μ) .> abstol)
    if sum(l) == 1
        k = findfirst(l)
        return nice_parameter(1/μ[k]) * basis[k], σ[k]
    else
        return nothing, Inf
    end
end

function find_dense(T, A, basis; abstol=1e-6)
    n = size(A, 1)
    b = ones(T, n)

    try
        q = A \ b
        if minimum(abs.(q)) > abstol
            ϵ = maximum(abs.(A*q .- b))
            if ϵ < abstol
                y = sum(nice_parameter.(q) .* basis)
                return y, ϵ
            end
        end
    catch e
        #
    end
    return nothing, Inf
end


text:
import Base: INVALID_OS_HANDLE
import Mmap

const PAGESIZE = Mmap.PAGESIZE

# Helpers to infer default parameters given a file handle/fd
import Mmap: Anonymous, gethandle, grow!

fileposition(::Anonymous) = Int64(0)
fileposition(io::IO) = Int64(position(io))

function filedim(io::IO, T::Type, offset::Union{Integer,Nothing} = nothing)
    offset = offset === nothing ? position(io) : offset
    return (div(filesize(io) - offset, sizeof(T)),)
end
function filedim(file::AbstractString, T::Type, offset::Union{Integer,Nothing} = nothing)
    offset = offset === nothing ? 0 : offset
    return (div(filesize(file) - offset, sizeof(T)),)
end

function filemode(fd::RawFD)
    F_GETFL = Cint(3)
    O_ACCMODE = Cint(3)
    O_RDONLY = 0
    O_WRONLY = 1
    O_RDWR   = 2
    fmode = ccall(:fcntl, Cint, (RawFD, Cint, Cint...), fd, F_GETFL)
    Base.systemerror("fcntl F_GETFL", fmode == -1)
    fmode &= O_ACCMODE
    return fmode == O_RDWR ? PROT_READ|PROT_WRITE :
           fmode == O_WRONLY ? PROT_WRITE : PROT_READ
end
filemode(io::IO) = filemode(gethandle(io))
filemode(::Anonymous) = PROT_READ | PROT_WRITE

fileflags(::Anonymous) = MAP_SHARED | MAP_ANONYMOUS
fileflags(::IO) = MAP_SHARED

###

# Raw call to mmap syscall
function _mmap(ptr::Ptr{Cvoid}, len::Int,
                                   prot::MmapProtection, flags::MmapFlags,
                                   fd::RawFD, offset::Int64)

    len < typemax(Int) - PAGESIZE || throw(
            ArgumentError("requested size must be < $(typemax(Int)-PAGESIZE), got $len"))
    offset >= 0 || throw(ArgumentError("requested offset must be ≥ 0, got $offset"))

    ret = ccall(:jl_mmap, Ptr{Cvoid}, (Ptr{Cvoid}, Csize_t, Cint, Cint, RawFD, Int64),
                ptr, len, prot, flags, fd, offset)
    Base.systemerror("mmap", reinterpret(Int, ret) == -1)
    return ret
end

function _unmap!(ptr::Ptr{Cvoid}, len::Int)
    ret = ccall(:munmap, Cint, (Ptr{Cvoid}, Csize_t), ptr, len)
    Base.systemerror("munmap", ret != 0)
    return
end

# Low-level form which mirrors a raw mmap, but constructs a Julia array of given
# dimension(s) at a specific offset within a file (includes accounting for page alignment
# requirement).
function mmap(::Type{Array{T}}, dims::NTuple{N,Integer},
              prot::MmapProtection, flags::MmapFlags,
              fd::RawFD, offset::Integer) where {T, N}
    isbitstype(T) || throw(ArgumentError("unable to mmap type $T; must satisfy `isbitstype(T) == true`"))

    len = prod(dims) * sizeof(T)
    iszero(len) && return Array{T}(undef, ntuple(x->0, Val(N)))
    len > 0 || throw(ArgumentError("requested size must be ≥ 0, got $len"))
    len = Int(len)

    page_pad = rem(Int64(offset), PAGESIZE)
    mmaplen::Int = len + page_pad

    ptr = _mmap(C_NULL, mmaplen, prot, flags, fd, Int64(offset) - page_pad)
    aptr = convert(Ptr{T}, ptr + page_pad)
    array = unsafe_wrap(Array{T,N}, aptr, dims)
    finalizer(_ -> _unmap!(ptr, mmaplen), array)
    return array
end
function mmap(::Type{Array{T}}, len::Int, prot::MmapProtection, flags::MmapFlags,
              fd::RawFD, offset::Int) where {T}
    return mmap(Array{T}, (Int(len),), prot, flags, fd, offset)
end

# Higher-level interface which takes an IO object and sets default flag values.
function mmap(io::IO, ::Type{Array{T}}, dims::NTuple{N,Integer};
              offset::Union{Integer,Nothing} = nothing,
              prot::Union{MmapProtection,Nothing} = nothing,
              flags::Union{MmapFlags,Nothing} = nothing,
              grow::Bool = true
             ) where {T, N}
    isbitstype(T) || throw(ArgumentError("unable to mmap type $T; must satisfy `isbitstype(T) == true`"))
    isopen(io) || throw(ArgumentError("$io must be open to mmap"))

    len = prod(dims) * sizeof(T)
    iszero(len) && return Array{T}(undef, ntuple(x->0, Val(N)))
    len > 0 || throw(ArgumentError("requested size must be ≥ 0, got $len"))
    len = Int(len)

    offset = offset === nothing ? fileposition(io) : Int(offset)
    prot   = prot === nothing ? filemode(io) : prot
    flags  = flags === nothing ? fileflags(io) : flags

    grow && iswritable(io) && grow!(io, offset, len)

    return mmap(Array{T}, dims, prot, flags, gethandle(io), offset)
end
function mmap(io::IO, ::Type{Array{T}}, len::Integer;
              offset::Union{Integer,Nothing} = nothing,
              prot::Union{MmapProtection,Nothing} = nothing,
              flags::Union{MmapFlags,Nothing} = nothing,
              grow::Bool = true
             ) where {T}
    return mmap(io, Array{T}, (len,);
                offset = offset, prot = prot, flags = flags, grow = grow)
end
function mmap(io::IO, ::Type{Array{T}} = Array{UInt8};
              offset::Union{Integer,Nothing} = nothing,
              prot::Union{MmapProtection,Nothing} = nothing,
              flags::Union{MmapFlags,Nothing} = nothing,
              grow::Bool = true
             ) where {T}
    return mmap(io, Array{T}, filedim(io, T, offset);
                offset = offset, prot = prot, flags = flags, grow = grow)
end

# Mapping of files
function mmap(file::AbstractString, ::Type{Array{T}}, dims::NTuple{N,Integer};
              offset::Union{Integer,Nothing} = nothing,
              prot::Union{MmapProtection,Nothing} = nothing,
              flags::Union{MmapFlags,Nothing} = nothing,
              grow::Bool = true
             ) where {T, N}
    if prot === nothing
        # Default to read-only if file exists, otherwise read/write/create
        openmode = isfile(file) ? "r" : "w+"
    else
        # If protection requested write...
        if (prot & PROT_WRITE) == PROT_WRITE
            # read/write if file exists, otherwise read/write/create
            openmode = isfile(file) ? "r+" : "w+"
        else
            # otherwise read-only
            openmode = "r"
        end
    end

    return open(file, openmode) do io
        mmap(io, Array{T}, dims;
             offset = offset, prot = prot, flags = flags, grow = grow)
    end
end
function mmap(file::AbstractString, ::Type{Array{T}}, len::Integer;
              offset::Union{Integer,Nothing} = nothing,
              prot::Union{MmapProtection,Nothing} = nothing,
              flags::Union{MmapFlags,Nothing} = nothing,
              grow::Bool = true
             ) where {T}
    return mmap(file, Array{T}, (len,);
                offset = offset, prot = prot, flags = flags, grow = grow)
end
# Default mapping of the [rest of] given file
function mmap(file::AbstractString, ::Type{Array{T}} = Array{UInt8};
              offset::Union{Integer,Nothing} = nothing,
              prot::Union{MmapProtection,Nothing} = nothing,
              flags::Union{MmapFlags,Nothing} = nothing,
              grow::Bool = true
             ) where {T}
    return mmap(file, Array{T}, filedim(file, T, offset);
                offset = offset, prot = prot, flags = flags, grow = grow)
end

# form to construct anonymous memory maps
function mmap(::Type{Array{T}}, dims::NTuple{N,Integer};
              prot::Union{MmapProtection,Nothing} = nothing,
              flags::Union{MmapFlags,Nothing} = nothing
             ) where {T, N}
    prot = prot === nothing ? PROT_READ | PROT_WRITE : prot
    flags = MAP_ANONYMOUS | (flags === nothing ? MAP_SHARED : flags)
    return mmap(Anonymous(), Array{T}, dims;
                offset = Int64(0), prot = prot, flags = flags, grow = false)
end
function mmap(::Type{Array{T}}, len::Integer;
              prot::Union{MmapProtection,Nothing} = nothing,
              flags::Union{MmapFlags,Nothing} = nothing
             ) where {T}
    return mmap(Array{T}, (len,);
                prot = prot, flags = flags)
end

"""
    mmap(file::Union{IO,AbstractString}[, ::Type{Array{T}}, dims]; kws...)
    mmap(::Type{Array{T}}, dims; kws...)

Memory maps an `Array` from a file `file` of size `dims` (either a scalar `Integer` or
a tuple of dimension lengths), or creates an anonymous mapping not backed by a file if
no file is given.

If not specified, the array type defaults to `Array{UInt8}`, and `dim` defaults to a
vector with length equal to the number of elements remaining in the file (accounting for
a non-zero position in the file stream `io`).

# Extended help

This function provides a relatively simple wrapper around the underlying `mmap` system
call. In particular, it is the user's responsibility to ensure the combination(s) of
stream state (e.g. read/write or read-only, file length), protection flags, and memory map
flags must form a valid `mmap` request, as they are not validated here.

See your system `mmap` man page for details on the behaviors of each flag.

## Keywords

The following keywords are available:

* `flags::MmapFlags` — any of the `MAP_*` system-specific constants. For files, this
  defaults to `MAP_SHARED`, and anonymous mappings default to `MAP_SHARED | MAP_ANONYMOUS`.

* `prot::MmapProtection` — any of the `PROT_*` flags. The default is `PROT_READ |
  PROT_WRITE` for all anonymous maps and memory maps of non-existent files (in which case
  the file will be created). If the file already exists, it is opened read-only and the
  default is `PROT_READ` only.

* `offset::Integer` — offset _in bytes_ from the beginning of the file, defaulting to the current
  stream position. This keyword is not valid for anonymous maps.

* `grow::Bool` — Whether to grow a file to accomodate the memory map, if the file is
  writable. Defaults to `true`. This keyword is not valid for anonymous maps.
"""
function mmap end

function pagepointer(array::Array)
    ptr = pointer(array)
    off = convert(UInt, ptr) % PAGESIZE
    ptr -= off
    len = off + sizeof(array)
    return (ptr, off, len)
end

"""
    mincore(array::Array)

Returns a boolean array (of length `fld1(sizeof(m), PAGESIZE)`) indicating whether the
pages of the memory-mapped `array` are resident in RAM. Pages corresponding to `false`
values will cause a fault if referenced.

!!! note
    Memory maps are always page-aligned, so an offset which is not a multiple of `PAGESIZE`
    will result in a return array where the first element does _not_ correspond to the
    first `PAGESIZE ÷ sizeof(eltype(array))` elements of `array` — i.e. the first page
    may cover memory addresses before the first element of `array`.
"""
function mincore(array::Array)
    ptr, _, len = pagepointer(array)
    sz = fld1(len, PAGESIZE)
    mask = Vector{Bool}(undef, sz)
    GC.@preserve array begin
        Base.systemerror("mincore",
                ccall(:mincore, Cint, (Ptr{Cvoid}, Csize_t, Ptr{UInt8}),
                      ptr, len, mask) != 0)
    end
    return mask
end

"""
    madvise!(array, flag::AdviseFlags = MADV_NORMAL)

Advises the kernel on the intended usage of the memory-mapped `array`, with the intent
`flag` being one of the available `MADV_*` constants.
"""
function madvise!(array::Array, flag::AdviseFlags = MADV_NORMAL)
    GC.@preserve array begin
        ptr, off, len = pagepointer(array)
        Base.systemerror("madvise",
                ccall(:madvise, Cint, (Ptr{Cvoid}, Csize_t, Cint),
                      ptr, len, flag) != 0)
    end
    return array
end


text:
using Distributed
@everywhere using MCPhylo, SpecialFunctions

## Data
magnesium = Dict{Symbol, Any}(
  :rt => [1, 9, 2, 1, 10, 1, 1, 90],
  :nt => [40, 135, 200, 48, 150, 59, 25, 1159],
  :rc => [2, 23, 7, 1, 8, 9, 3, 118],
  :nc => [36, 135, 200, 46, 148, 56, 23, 1157]
)

magnesium[:rtx] = permutedims(hcat([magnesium[:rt] for i in 1:6]...))
magnesium[:rcx] = permutedims(hcat([magnesium[:rc] for i in 1:6]...))
magnesium[:s2] = 1 ./ (magnesium[:rt] .+ 0.5) +
                 1 ./ (magnesium[:nt] .- magnesium[:rt] .+ 0.5) +
                 1 ./ (magnesium[:rc] .+ 0.5) +
                 1 ./ (magnesium[:nc] .- magnesium[:rc] .+ 0.5)
magnesium[:s2_0] = 1 / mean(1 ./ magnesium[:s2])

## Model Specification
model = Model(

  rcx = Stochastic(2,
    (nc, pc) ->
      UnivariateDistribution[Binomial(nc[j], pc[i, j]) for i in 1:6, j in 1:8],
    false
  ),

  pc = Stochastic(2,
    () -> Uniform(0, 1),
    false
  ),

  rtx = Stochastic(2,
    (nt, pc, theta) ->
      UnivariateDistribution[(
        phi = logit(pc[i, j]);
        pt = invlogit(theta[i, j] + phi);
        Binomial(nt[j], pt)) for i in 1:6, j in 1:8
      ],
    false
  ),

  theta = Stochastic(2,
    (mu, tau) ->
      UnivariateDistribution[Normal(mu[i], tau[i]) for i in 1:6, j in 1:8],
    false
  ),

  mu = Stochastic(1,
    () -> Uniform(-10, 10),
    false
  ),

  OR = Logical(1,
    mu -> exp.(mu)
  ),

  tau = Logical(1,
    (priors, s2_0) ->
      Float64[
        sqrt(priors[1]),
        sqrt(priors[2]),
        priors[3],
        sqrt(s2_0 * (1 / priors[4] - 1)),
        sqrt(s2_0) * (1 / priors[5] - 1),
        sqrt(priors[6]) ]
  ),

  priors = Stochastic(1,
    s2_0 ->
      UnivariateDistribution[
        InverseGamma(0.001, 0.001),
        Uniform(0, 50),
        Uniform(0, 50),
        Uniform(0, 1),
        Uniform(0, 1),
        Truncated(Normal(0, sqrt(s2_0 / erf(0.75))), 0, Inf)
      ],
    false
  )

)

## Initial Values
inits = [
  Dict(:rcx => magnesium[:rcx], :rtx => magnesium[:rtx],
       :theta => zeros(6, 8), :mu => fill(-0.5, 6),
       :pc => fill(0.5, 6, 8), :priors => [1, 1, 1, 0.5, 0.5, 1]),
  Dict(:rcx => magnesium[:rcx], :rtx => magnesium[:rtx],
       :theta => zeros(6, 8), :mu => fill(0.5, 6),
       :pc => fill(0.5, 6, 8), :priors => [1, 1, 1, 0.5, 0.5, 1])
]

## Sampling Scheme
scheme = [AMWG(:theta, 0.1),
          AMWG(:mu, 0.1),
          Slice(:pc, 0.25, Univariate),
          Slice(:priors, [1.0, 5.0, 5.0, 0.25, 0.25, 5.0], Univariate)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, magnesium, inits, 12500, burnin=2500, thin=2, chains=2)
describe(sim)


text:
# Command: New Window
"""
    window!(session::Session)

Create a new top-level browsing context.
"""
function window!(session::Session)::NamedTuple{(:type, :handle),NTuple{2, String}}
    @unpack addr, id = session
    response =
        HTTP.post("$addr/session/$id/window/new", [("Content-Type" => "application/json; charset=utf-8")])
    @assert response.status == 200
    result = JSON3.read(response.body, NamedTuple)
    @assert result.status == 0
    (type = result.value["type"], handle = result.value["handle"])
end


text:
include("../../../utils.jl")

# parameters that govern counterfactual


text:
using LinearAlgebra
using BaryRational
using Test
using Random

include("test_aaa.jl")
include("test_fh.jl")

@testset "BaryRational.jl" begin
    @testset "aaa_rational" begin
        @test test_aaa_gamma_poles()
        @test test_aaa_spiral()
        @test test_aaa_exp()
        @test test_aaa_length()
        @test test_aaa_tan()
        @test test_aaa_gamma()
        @test test_aaa_ranlog()
        @test test_aaa_infinp()
        @test test_aaa_naninp()
        @test test_aaa_residuals()
        #@test test_aaa_case2()
        #@test test_aaa_scale_invar()
    end
    @testset "FH_rational_interpolation" begin
        @test test_fh_runge()
        @test test_abs_x()
    end
end


text:
using KitePodSimulator
using Documenter

DocMeta.setdocmeta!(KitePodSimulator, :DocTestSetup, :(using KitePodSimulator); recursive=true)

makedocs(;
    modules=[KitePodSimulator],
    authors="Uwe Fechner <fechner@aenarete.eu> and contributors",
    repo="https://github.com/ufechner7/KitePodSimulator.jl/blob/{commit}{path}#{line}",
    sitename="KitePodSimulator.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ufechner7.github.io/KitePodSimulator.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Functions" => "functions.md",
        "Parameters" => "parameters.md",
    ],
)

deploydocs(;
    repo="github.com/ufechner7/KitePodSimulator.jl",
    devbranch="main",
)


text:
using WriteVTK

function vtk_connectivity_map_highorder(Nq1, Nq2 = 1, Nq3 = 1)
    connectivity = Array{Int, 1}(undef, Nq1 * Nq2 * Nq3)
    L = LinearIndices((1:Nq1, 1:Nq2, 1:Nq3))

    corners = (
        (1, 1, 1),
        (Nq1, 1, 1),
        (Nq1, Nq2, 1),
        (1, Nq2, 1),
        (1, 1, Nq3),
        (Nq1, 1, Nq3),
        (Nq1, Nq2, Nq3),
        (1, Nq2, Nq3),
    )
    edges = (
        (2:(Nq1 - 1), 1:1, 1:1),
        (Nq1:Nq1, 2:(Nq2 - 1), 1:1),
        (2:(Nq1 - 1), Nq2:Nq2, 1:1),
        (1:1, 2:(Nq2 - 1), 1:1, 1:1),
        (2:(Nq1 - 1), 1:1, Nq3:Nq3),
        (Nq1:Nq1, 2:(Nq2 - 1), Nq3:Nq3),
        (2:(Nq1 - 1), Nq2:Nq2, Nq3:Nq3),
        (1:1, 2:(Nq2 - 1), Nq3:Nq3),
        (1:1, 1:1, 2:(Nq3 - 1)),
        (Nq1:Nq1, 1:1, 2:(Nq3 - 1)),
        (1:1, Nq2:Nq2, 2:(Nq3 - 1)),
        (Nq1:Nq1, Nq2:Nq2, 2:(Nq3 - 1)),
    )
    faces = (
        (1:1, 2:(Nq2 - 1), 2:(Nq3 - 1)),
        (Nq1:Nq1, 2:(Nq2 - 1), 2:(Nq3 - 1)),
        (2:(Nq1 - 1), 1:1, 2:(Nq3 - 1)),
        (2:(Nq1 - 1), Nq2:Nq2, 2:(Nq3 - 1)),
        (2:(Nq1 - 1), 2:(Nq2 - 1), 1:1),
        (2:(Nq1 - 1), 2:(Nq2 - 1), Nq3:Nq3),
    )
    if Nq2 == Nq3 == 1
        @assert Nq1 > 1
        corners = (corners[1:2]...,)
        edges = (edges[1],)
        faces = ()
    elseif Nq3 == 1
        @assert Nq1 > 1
        @assert Nq2 > 1
        corners = (corners[1:4]...,)
        edges = (edges[1:4]...,)
        faces = (faces[5],)
    end

    # corners
    iter = 1
    for (i, j, k) in corners
        connectivity[iter] = L[i, j, k]
        iter += 1
    end
    # edges
    for (is, js, ks) in edges
        for k in ks, j in js, i in is
            connectivity[iter] = L[i, j, k]
            iter += 1
        end
    end
    # faces
    for (is, js, ks) in faces
        for k in ks, j in js, i in is
            connectivity[iter] = L[i, j, k]
            iter += 1
        end
    end
    # interior
    for k in 2:(Nq3 - 1), j in 2:(Nq2 - 1), i in 2:(Nq1 - 1)
        connectivity[iter] = L[i, j, k]
        iter += 1
    end
    return connectivity
end

#=
This is the 1D WriteMesh routine
=#
function writemesh_highorder(
    base_name,
    x1;
    x2 = nothing,
    x3 = nothing,
    fields = (),
    realelems = 1:size(x1)[end],
)
    (Nq1, _) = size(x1)

    con = vtk_connectivity_map_highorder(Nq1)

    M = MeshCell{VTKCellTypes.VTKCellType, Array{Int, 1}}
    cells = Array{M, 1}(undef, length(realelems))

    for (i, e) in enumerate(realelems)
        offset = (e - 1) * Nq1
        cells[i] = MeshCell(VTKCellTypes.VTK_LAGRANGE_CURVE, offset .+ con)
    end

    if x2 == nothing
        @assert isnothing(x3)
        vtkfile =
            vtk_grid("$(base_name)", @view(x1[:]), cells; compress = false)
    elseif x3 == nothing
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            cells;
            compress = false,
        )
    else
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            @view(x3[:]),
            cells;
            compress = false,
        )
    end
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

#=
This is the 2D WriteMesh routine
=#
function writemesh_highorder(
    base_name,
    x1,
    x2;
    x3 = nothing,
    fields = (),
    realelems = 1:size(x1)[end],
)
    @assert size(x1) == size(x2)
    (Nq1, Nq2, _) = size(x1)
    @assert Nq1 == Nq2
    con = vtk_connectivity_map_highorder(Nq1, Nq2)

    M = MeshCell{VTKCellTypes.VTKCellType, Array{Int, 1}}
    cells = Array{M, 1}(undef, length(realelems))
    for (i, e) in enumerate(realelems)
        offset = (e - 1) * Nq1 * Nq2
        cells[i] =
            MeshCell(VTKCellTypes.VTK_LAGRANGE_QUADRILATERAL, offset .+ con[:])
    end

    if x3 == nothing
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            cells;
            compress = false,
        )
    else
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            @view(x3[:]),
            cells;
            compress = false,
        )
    end
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

#=
This is the 3D WriteMesh routine
=#
function writemesh_highorder(
    base_name,
    x1,
    x2,
    x3;
    fields = (),
    realelems = 1:size(x1)[end],
)
    (Nq1, Nq2, Nq3, _) = size(x1)
    @assert Nq1 == Nq2 == Nq3
    con = vtk_connectivity_map_highorder(Nq1, Nq2, Nq3)
    M = MeshCell{VTKCellTypes.VTKCellType, Array{Int, 1}}
    cells = Array{M, 1}(undef, length(realelems))
    for (i, e) in enumerate(realelems)
        offset = (e - 1) * Nq1 * Nq2 * Nq3
        cells[i] = MeshCell(VTKCellTypes.VTK_LAGRANGE_HEXAHEDRON, offset .+ con)
    end

    vtkfile = vtk_grid(
        "$(base_name)",
        @view(x1[:]),
        @view(x2[:]),
        @view(x3[:]),
        cells;
        compress = false,
    )
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

#=
This is the 1D WriteMesh routine
=#
function writemesh_raw(
    base_name,
    x1;
    x2 = nothing,
    x3 = nothing,
    fields = (),
    realelems = 1:size(x1)[end],
)
    (Nq1, _) = size(x1)
    Nsubcells = (Nq1 - 1)

    M = MeshCell{VTKCellTypes.VTKCellType, Array{Int, 1}}
    cells = Array{M, 1}(undef, Nsubcells * length(realelems))
    for e in realelems
        offset = (e - 1) * Nq1
        for i in 1:(Nq1 - 1)
            cells[i + (e - 1) * Nsubcells] =
                MeshCell(VTKCellTypes.VTK_LINE, offset .+ [i, i + 1])
        end
    end

    vtkfile = vtk_grid("$(base_name)", @view(x1[:]), cells; compress = false)
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

#=
This is the 2D WriteMesh routine
=#
function writemesh_raw(
    base_name,
    x1,
    x2;
    x3 = nothing,
    fields = (),
    realelems = 1:size(x1)[end],
)
    @assert size(x1) == size(x2)
    (Nq1, Nq2, _) = size(x1)
    Nsubcells = (Nq1 - 1) * (Nq2 - 1)

    M = MeshCell{VTKCellTypes.VTKCellType, Array{Int, 1}}
    cells = Array{M, 1}(undef, Nsubcells * length(realelems))
    ind = LinearIndices((1:Nq1, 1:Nq2))
    for e in realelems
        offset = (e - 1) * Nq1 * Nq2
        for j in 1:(Nq2 - 1)
            for i in 1:(Nq1 - 1)
                cells[i + (j - 1) * (Nq1 - 1) + (e - 1) * Nsubcells] = MeshCell(
                    VTKCellTypes.VTK_PIXEL,
                    offset .+ ind[i:(i + 1), j:(j + 1)][:],
                )
            end
        end
    end

    if x3 == nothing
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            cells;
            compress = false,
        )
    else
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            @view(x3[:]),
            cells;
            compress = false,
        )
    end
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

#=
This is the 3D WriteMesh routine
=#
function writemesh_raw(
    base_name,
    x1,
    x2,
    x3;
    fields = (),
    realelems = 1:size(x1)[end],
)
    (Nq1, Nq2, Nq3, _) = size(x1)
    (N1, N2, N3) = (Nq1 - 1, Nq2 - 1, Nq3 - 1)
    Nsubcells = N1 * N2 * N3

    M = MeshCell{VTKCellTypes.VTKCellType, Array{Int, 1}}
    cells = Array{M, 1}(undef, Nsubcells * length(realelems))
    ind = LinearIndices((1:Nq1, 1:Nq2, 1:Nq3))
    for e in realelems
        offset = (e - 1) * Nq1 * Nq2 * Nq3
        for k in 1:N3
            for j in 1:N2
                for i in 1:N1
                    cells[i + (j - 1) * N1 + (k - 1) * N1 * N2 + (e - 1) * Nsubcells] =
                        MeshCell(
                            VTKCellTypes.VTK_VOXEL,
                            offset .+ ind[i:(i + 1), j:(j + 1), k:(k + 1)][:],
                        )
                end
            end
        end
    end

    vtkfile = vtk_grid(
        "$(base_name)",
        @view(x1[:]),
        @view(x2[:]),
        @view(x3[:]),
        cells;
        compress = false,
    )
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end


text:
module LibSDL2

using CEnum: @cenum

# This is temporary, don't worry.
if isfile("/usr/lib/libSDL2.so")
    @info "Using system libSDL2"
    # Override this with `@eval LibSDL2 const libsdl2 = $(SDL2_jll.libsdl2)`
    const libsdl2 = "libSDL2"
else
    @info "Using SDL2_jll"
    @eval using SDL2_jll
    const libsdl2 = SDL2_jll.libsdl2
end

# Forward declare some methods.
function UpperBlit end
function UpperBlitScaled end
function BlitMap end

# This type is not supported by Clang.jl or Julia yet.
const AudioCVT = Cvoid

include(joinpath(@__DIR__, "..", "gen", "bindings", "libsdl2_types.jl"))
include(joinpath(@__DIR__, "..", "gen", "bindings", "libsdl2_funcs.jl"))
# Low stakes gripe: This adds 2s to precompilation time, but I could fully generate it
# ahead of time with a bit of effort.
include(joinpath(@__DIR__, "..", "gen", "bindings", "libsdl2_docs.jl"))

# Macros ported by NHDaly
MUSTLOCK(S::Ptr{Surface}) = ((unsafe_load(S).flags & RLEACCEL) != 0)
LoadBMP(file) = LoadBMP_RW(RWFromFile(file, "rb"), 1)
SaveBMP(surface, file) = SaveBMP_RW(surface, RWFromFile(file, "wb" ), 1)

WINDOWPOS_CENTERED_DISPLAY(X) = ( WINDOWPOS_CENTERED_MASK | (X) )
WINDOWPOS_CENTERED() = WINDOWPOS_CENTERED_DISPLAY(0)
WINDOWPOS_ISCENTERED(X) = ( ( (X) & 0xFFFF0000 ) == WINDOWPOS_CENTERED_MASK )

SCANCODE_TO_KEYCODE(X) = (X | SDLK_SCANCODE_MASK)

end # module