text:
# === PRIVATE FUNCTIONS THAT ARE NOT EXPORTED =============================================================== #
# =========================================================================================================== #

# === PUBLIC FUNCTIONS THAT ARE EXPORTED ==================================================================== #
"""
    z_score_transform_data_set(dataFrame::DataFrame,colNameArray::Array{String,1})::VLResult
"""
function z_score_transform_data_set(dataFrame::DataFrame,colNameArray::Array{String,1})::VLResult

    # initialize -
    (number_of_rows, number_of_cols) = size(dataFrame)
    z_score_array = Array{Float64,2}(undef,number_of_rows,length(colNameArray))

    try

        # ok:  for each col, lets z-center the data -
        for (col_index,col_name) in enumerate(colNameArray)

            # grab the data -
            data_col = dataFrame[!,col_name]

            # transform -
            tmp = parse.(Float64,data_col)

            # compute the mean, std -
            µ = mean(tmp)
            σ = std(tmp)

            # compute the z-score -
            z_score_col = (tmp .- μ)./σ

            # add -
            for (row_index,value) in enumerate(z_score_col)
                z_score_array[row_index,col_index] = value
            end
        end

        # return -
        return VLResult(z_score_array)
    catch error
        return VLResult(error)
    end
end

"""
    z_score_transform_array(dataArray::Array{Float64,2})::VLResult
"""
function z_score_transform_array(dataArray::Array{Float64,2})::VLResult

    # initialize -
    (number_of_rows, number_of_cols) = size(dataArray)
    z_score_array = Array{Float64,2}(undef, number_of_rows, number_of_cols)

    try

        for col_index = 1:number_of_cols

            # grab the data col -
            data_col = dataArray[:,col_index]

            # zscale this col -
            zscale_vector_result = z_score_transform_vector(data_col)
            if (isa(zscale_vector_result.value, Exception) == true)
                throw(zscale_vector_result.value)
            end
            scaled_col_data = zscale_vector_result.value

            # package, and go around again ...
            for (row_index, scaled_value) in enumerate(scaled_col_data)
                z_score_array[row_index,col_index] = scaled_value
            end
        end

        # return -
        return VLResult(z_score_array)
    catch error
        return VLResult(error)
    end
end

"""
    z_score_transform_vector(dataVector::Array{Float64,1})::VLResult
"""
function z_score_transform_vector(dataVector::Array{Float64,1})::VLResult

    # initalize -
    number_of_elements = length(dataVector)
    z_score_array = Array{Float64,1}(undef, number_of_elements)

    try

        # compute the mean, std -
        µ = mean(dataVector)
        σ = std(dataVector)

        # scale -
        for index = 1:number_of_elements

            # get the raw value -
            raw_value = dataVector[index]

            # scale -
            z_score_array[index] =  (raw_value - µ)/σ
        end

        # rerurn -
        return VLResult(z_score_array)

    catch error
        return VLResult(error)
    end
end
# =========================================================================================================== #

max_stars_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_issues_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_forks_repo_name: varnerlab/PregnancyCoagDataAnalysisKit


text:
using Dierckx, QuadGK

function int_tabulated(λ, f; int_k=3)
    @assert issorted(λ)
    itp = Spline1D(λ, f, k=int_k, bc="error")
    return quadgk(itp, λ[1], λ[end])
end

ext: jl


text:
function mgc2sp(mgc::AbstractVector, α=0.35, γ=0.0, fftlen=256)
    # transform to cesptrum
    c = mgc2mgc(mgc, α, γ, fftlen>>1, 0.0, 0.0)

    # zero padding
    buf = zeros(eltype(c), fftlen)
    copy!(buf, 1, c, 1, length(c))

    # FFT
    rfft(buf)
end

function mgc2sp{T<:AbstractFloat}(mgc::AbstractMatrix{T}, α=0.0, γ=0.0, fftlen=256)
    sp = Array{Complex{T}}(fftlen>>1 + 1, size(mgc, 2))
    for i in 1:size(mgc, 2)
        @inbounds sp[:,i] = mgc2sp(view(mgc, :, i), α, γ, fftlen)
    end
    sp
end

function mgc2sp{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T}, fftlen=256)
    assert_not_ready_to_filt(state)

    def = paramdef(state)
    α = allpass_alpha(def)
    γ = glog_gamma(def)
    data = rawdata(state)
    mgc2sp(data, α, γ, fftlen)
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(load, (File{format"2DM"},))
    @assert precompile(load, (File{format"MSH"},))
    @assert precompile(load, (File{format"OBJ"},))
    @assert precompile(load, (File{format"OFF"},))
    @assert precompile(load, (File{format"PLY_ASCII"},))
    @assert precompile(load, (File{format"STL_ASCII"},))
    @assert precompile(load, (File{format"STL_BINARY"},))
end


text:
module MayOptimizeLinearAlgebraBenchmarks

export
    run_benchmarks,
    load_benchmarks,
    save_benchmarks

const USE_AVX = false

using HDF5
using LinearAlgebra
using Statistics
using StaticArrays

using BenchmarkTools
using BenchmarkTools: Trial

using MayOptimize
using MayOptimize:
    AVX,
    Standard,
    CholeskyBanachiewiczLowerI,
    CholeskyBanachiewiczLowerII,
    CholeskyBanachiewiczUpper,
    CholeskyCroutLower,
    CholeskyCroutUpperI,
    CholeskyCroutUpperII,
    exec!

using MayOptimize.LinearAlgebraMethods:
    AbstractAlgorithm,
    CholeskyFactorization,
    Floats

import MayOptimize.LinearAlgebraMethods:
    exec!

struct LDiv{opt} <: AbstractAlgorithm
    operand::String
end

exec!(alg::LDiv{opt}, args...) where {opt} = ldiv!(opt, args...)

struct CholeskyBLAS <: AbstractAlgorithm end
struct CholeskyStatic <: AbstractAlgorithm end

group(T::Type{<:Floats}, alg::Type{<:AbstractAlgorithm}) = group(T, alg())

@noinline group(T::Type{<:Floats}, alg::AbstractAlgorithm) =
    string(T, "/", group(alg))

@noinline group(alg::LDiv{<:Standard}) = "ldiv!-$(alg.operand)/BLAS"
@noinline group(alg::LDiv{opt}) where {opt} =
    string("ldiv!-", alg.operand, "/", nameof(opt))

group(::CholeskyBLAS  ) = "Cholesky/BLAS"
group(::CholeskyStatic) = "Cholesky/StaticArrays"

const CholeskyMethod = Union{CholeskyFactorization,CholeskyBLAS,CholeskyStatic}
nops(::CholeskyMethod, n::Integer) = (n^3)/3

# Assume triangular matrix. Number of perations is: n divisions plus 1 - 2n +
# n² additions and multiplications, which is approximately n² (because
# divisions cost more than additions and multiplications).
nops(::LDiv, n::Integer) = n^2

for A in (:CholeskyBanachiewiczLowerII,
          :CholeskyBanachiewiczLowerI,
          :CholeskyBanachiewiczUpper,
          :CholeskyCroutLower,
          :CholeskyCroutUpperII,
          :CholeskyCroutUpperI)
    @eval @noinline group(::$A{opt}) where {opt} =
        string($(string(A)), "/", nameof(opt))
end

@static if USE_AVX
    include("../src/linalg-avx.jl")
end

# struct BenchmarkTools.Trial:
#  params  - BenchmarkTools.Parameters
#  times   - measured times in nanoseconds
#  gctimes - garbage collector times in nanoseconds
#  memory  - memory allocation in bytes
#  allocs  - number of allocations

struct BenchmarkEntry
    group::String
    len::Vector{Int}
    allocs::Vector{Int}
    memory::Vector{Int}
    nops::Vector{Float64}
    min::Vector{Float64}
    max::Vector{Float64}
    avg::Vector{Float64}
    std::Vector{Float64}
    median::Vector{Float64}
    BenchmarkEntry(group::String) = new(group, Int[], Int[], Int[],
                                        Float64[], Float64[], Float64[],
                                        Float64[], Float64[], Float64[])
end

function push_result!(e::BenchmarkEntry,
                      len::Integer,
                      nops::Union{Real,Function},
                      trial::Trial)
    push!(e.len, len)
    push!(e.allocs, trial.allocs)
    push!(e.memory, trial.memory)
    if isa(nops, Real)
        push!(e.nops, nops)
    else
        push!(e.nops, nops(len))
    end
    push!(e.min, minimum(trial.times))
    push!(e.max, maximum(trial.times))
    push!(e.avg, mean(trial.times))
    push!(e.std, std(trial.times; corrected = true))
    push!(e.median, median(trial.times))
    return e
end

get_L(A::Cholesky) = LowerTriangular(convert(Matrix, A.L))
get_U(A::Cholesky) = UpperTriangular(convert(Matrix, A.U))

function save_benchmarks(filename::String,
                         database::AbstractDict{String,BenchmarkEntry})
    # Auxiliary function.
    @inline function write_field!(dst::HDF5.Group,
                                  src::BenchmarkEntry,
                                  ::Val{sym}) where {sym}
        dst[string(sym)] = getfield(src, sym)
        nothing
    end
    temp = filename*".tmp"
    h5open(temp, "w") do f
        for group in keys(database)
            e = database[group]
            g = create_group(f, group)
            write_field!(g, e, Val(:len))
            write_field!(g, e, Val(:allocs))
            write_field!(g, e, Val(:memory))
            write_field!(g, e, Val(:nops))
            write_field!(g, e, Val(:min))
            write_field!(g, e, Val(:max))
            write_field!(g, e, Val(:avg))
            write_field!(g, e, Val(:std))
            write_field!(g, e, Val(:median))
            close(g)
        end
    end
    mv(temp, filename; force=true)
    nothing
end

load_benchmarks(filename::String) =
    load_benchmarks(Dict{String,BenchmarkEntry}(), filename)

function load_benchmarks(database::AbstractDict{String,BenchmarkEntry},
                         filename::String)
    # Auxiliary function.
    @inline function read_field!(dst::BenchmarkEntry,
                                 src::HDF5.Group,
                                 ::Val{sym},
                                 I::Vector{Int}) where {sym}
        A = read(src, string(sym))
        B = resize!(getfield(dst, sym), length(A))
        for i in eachindex(B, I)
            B[i] = A[I[i]]
        end
        nothing
    end
    h5open(filename, "r") do f
        for type in keys(f)
            type_grp = f[type]
            for alg in keys(type_grp)
                alg_grp = type_grp[alg]
                for opt in keys(alg_grp)
                    g = alg_grp[opt]
                    len = read(g, "len")
                    i = sortperm(len)
                    e = BenchmarkEntry(type*"/"*alg*"/"*opt)
                    read_field!(e, g, Val(:len),    i)
                    read_field!(e, g, Val(:allocs), i)
                    read_field!(e, g, Val(:memory), i)
                    read_field!(e, g, Val(:nops),   i)
                    read_field!(e, g, Val(:min),    i)
                    read_field!(e, g, Val(:max),    i)
                    read_field!(e, g, Val(:avg),    i)
                    read_field!(e, g, Val(:std),    i)
                    read_field!(e, g, Val(:median), i)
                    database[e.group] = e
                end
            end
        end
    end
    return database
end

function run_benchmarks(filename::AbstractString,
                        alg::Union{AbstractString,Symbol}; kwds...)
    return run_benchmarks(filename, Symbol(alg); kwds...)
end

function run_benchmarks(filename::AbstractString,
                        alg::Symbol; kwds...)
    if alg === :ldiv
        return run_benchmarks(
            filename,
            map(x->LDiv{x}, (Standard,Debug,InBounds,Vectorize))...;
            kwds...)
    elseif alg === :Cholesky
        algs = DataType[]
        push!(algs, CholeskyBLAS)
        #push!(algs, CholeskyStatic)
        for A in (CholeskyBanachiewiczLowerII,
                  CholeskyBanachiewiczLowerI,
                  CholeskyBanachiewiczUpper,
                  CholeskyCroutLower,
                  CholeskyCroutUpperII,
                  CholeskyCroutUpperI)
            for opt in (Debug, InBounds, Vectorize)
                push!(algs, A{opt})
            end
        end
        return run_benchmarks(filename, algs...; kwds...)
    else
        error("unknown algorithm ", alg)
    end
end

function run_benchmarks(filename::AbstractString,
                        args::Type{<:AbstractAlgorithm}...;
                        sizes = 10:10:300,
                        types = (Float32,))
    db = Dict{String,BenchmarkEntry}()
    if isfile(filename)
        load_benchmarks(db, filename)
    end
    for T in types, n in sizes
        run_benchmarks(filename, db, T, Int(n), args...)
    end
    return db
end

function run_benchmarks(filename::AbstractString,
                        db::AbstractDict{String,BenchmarkEntry},
                        T::Type{<:Floats},
                        n::Int,
                        args::Type{<:AbstractAlgorithm}...)
    A = generate_positive_definite_matrix(T, n)
    C = cholesky(A)
    x = generate_array(T, n)
    for X in args
        run_benchmarks(filename, db, X, A, C, x)
    end
    return db
end

instances(X::Type{<:AbstractAlgorithm}) = (X(),)
instances(X::Type{<:LDiv}) = (X("L"), X("L'"), X("R"), X("R'"))

function run_benchmarks(filename::AbstractString,
                        db::AbstractDict{String,BenchmarkEntry},
                        X::Type{<:AbstractAlgorithm},
                        A::Matrix{T},
                        C::Cholesky{T},
                        x::Vector{T}) where {T<:Floats}
    n = length(x)
    for alg in instances(X)
        grp = group(T, alg)
        if haskey(db, grp)
            if n ∈ db[grp].len
                io = stderr
                printstyled(io, "warning:"; color=:yellow)
                print(io, " skipping test for ")
                printstyled(io, grp; color=:cyan)
                print(io, " with ")
                printstyled(io, "n=", n; color=:cyan)
                print(io, " element(s)\n")
                continue
            end
            e = db[grp]
        else
            e = BenchmarkEntry(grp)
            db[grp] = e
        end
        run_benchmark(e, alg, A, C, x)
        save_benchmarks(filename, db)
    end
    return db
end

function run_benchmark(e::BenchmarkEntry,
                       alg::LDiv{opt},
                       A::Matrix{T},
                       C::Cholesky{T},
                       x::Vector{T}) where {T<:Floats,opt}
    n = length(x)
    w = similar(x)
    if alg.operand == "L"
        L = get_L(C)
        push_result!(e, n, nops(alg, n),
                     @benchmark ldiv!($opt, $w, $(L),  $x))
    elseif alg.operand == "L'"
        L = get_L(C)
        push_result!(e, n, nops(alg, n),
                     @benchmark ldiv!($opt, $w, $(L'), $x))
    elseif alg.operand == "R"
        R = get_U(C)
        push_result!(e, n, nops(alg, n),
                     @benchmark ldiv!($opt, $w, $(R),  $x))
    elseif alg.operand == "R'"
        R = get_U(C)
        push_result!(e, n, nops(alg, n),
                     @benchmark ldiv!($opt, $w, $(R'), $x))
    end
    return e
end

function run_benchmark(e::BenchmarkEntry,
                       alg::CholeskyBLAS,
                       A::Matrix{T},
                       C::Cholesky{T},
                       x::Vector{T}) where {T<:Floats}
    n = length(x)
    B = similar(A)
    opt = Standard
    push_result!(e, n, nops(alg, n), @benchmark cholesky!($opt, $B, $A))
end

function run_benchmark(e::BenchmarkEntry,
                       alg::CholeskyStatic,
                       A::Matrix{T},
                       C::Cholesky{T},
                       x::Vector{T}) where {T<:Floats}
    run_benchmark(e, alg, A, C, x, Val(length(x)))
end

function run_benchmark(e::BenchmarkEntry,
                       alg::CholeskyStatic,
                       A::Matrix{T},
                       C::Cholesky{T},
                       x::Vector{T},
                       ::Val{n}) where {T<:Floats,n}
    Asm = SMatrix{n,n}(A)
    push_result!(e, n, nops(alg, n), @benchmark cholesky($Asm))
end

function run_benchmark(e::BenchmarkEntry,
                       alg::CholeskyFactorization,
                       A::Matrix{T},
                       C::Cholesky{T},
                       x::Vector{T}) where {T<:Floats}
    n = length(x)
    B = similar(A)
    push_result!(e, n, nops(alg, n), @benchmark cholesky!($alg,  $B, $A))
end

# Generate random values uniformly distributed in [-1,1].
generate_array(T::Type{<:Real}, dims...) = 2*rand(T, dims...) .- 1
generate_array(T::Type{<:Complex{<:Real}}, dims...) =
    2*rand(T, dims...) .- Complex(1,1)

# Build a random Hermitian positive definite matrix.
function generate_positive_definite_matrix(T::Type, n::Integer)
    H = generate_array(T, 10n, n)
    A = H'*H
    return (A + A')/2
end

#plt.figure(2); plt.clf();for lab in ("Standard","Debug","InBounds","Vectorize"); dat = q["ldiv!,$lab,L'"];plt.plot(dat.len, dat.nops ./ dat.median, label=lab); end; plt.legend();plt.xlabel("size");plt.ylabel("Gflops (median)");plt.title("Left division by transpose of lower triangular matrix");savefig("ldiv-Lt-median.png",bbox_inches="tight")

#figure(2); clf(2); db = BenchmarkingMayOptimizeLinearAlgebraMethods.load_benchmarks("benchmarks-1.h5");clf(); for o in ("BLAS","Debug","InBounds","Vectorize"); g="Float32/ldiv!-R/$o"; plot(db[g].len, db[g].nops./db[g].median,label=g); end; legend();

# Cholesky:
# BLAS is:
# 11.3 Gflops for Float64 for n=200
# 10.7 Gflops for Float32 for n=200

# CholeskyCroutUpperII
# 8.77 Gflops for Float64 for n=200
# 11.3 Gflops for Float32 for n=200

# CholeskyBanachiewiczUpper (this is the fastest):
# 9.34 Gflops for Float64 for n=200
# 13.0 Gflops for Float32 for n=200

make_plot(filename::AbstractString, args...; kwds...) =
     make_plot(load_benchmarks(filename), args...; kwds...)

function make_plot(db::AbstractDict{String,BenchmarkEntry}, plt;
                   type = "Float32", op = :ldiv, field = :median,
                   first_figure = 1,
                   save::Bool = false, dir = ".", ext = ".png",
                   dpi = 120, pad_inches = 5/dpi)

    fig = first_figure - 1
    if op === :ldiv
        for A in ("L", "L'", "R", "R'")
            fig += 1
            plt.figure(fig)
            plt.clf()
            plt.tight_layout()
            for o in ("BLAS", "Debug", "InBounds", "Vectorize")
                g = "$type/ldiv!-$A/$o"
                plt.plot(db[g].len, db[g].nops./getfield(db[g], field), label=o);
            end
            plt.legend()
            plt.title(string("Left division by ",
                             (endswith(A, "'") ? "transpose of " : ""),
                             (endswith(A, r"L'?") ? "lower" : "upper"),
                             " triangular matrix"))
            if save
                filename = joinpath(
                    dir, string("ldiv-", replace(A, r"'$" => "t"), "-",
                                field, ext))
                plt.savefig(filename; dpi=dpi, pad_inches=pad_inches)
            end
        end
    elseif op === :cholesky
        for opt in (:Debug, :InBounds, :Vectorize)
            fig += 1
            plt.figure(fig)
            plt.clf()
            plt.tight_layout()
            for f in (:BLAS,
                      :CholeskyBanachiewiczLowerI,
                      :CholeskyBanachiewiczLowerII,
                      :CholeskyBanachiewiczUpper,
                      :CholeskyCroutLower,
                      :CholeskyCroutUpperI,
                      :CholeskyCroutUpperII)
                g = string(type, (f === :BLAS ? "/Cholesky/$f" : "/$f/$opt"))
                plt.plot(db[g].len, db[g].nops./getfield(db[g], field), label="$f");
            end
            plt.legend()
            plt.title("Cholesky decomposition (optimization level: `$opt`)")
            if save
                filename = joinpath(
                    dir, string("cholesky-", lowercase(string(opt)), "-",
                                field, ext))
                plt.savefig(filename; dpi=dpi, pad_inches=pad_inches)
            end
        end
    else
        error("unknown operation: ", op)
    end
end

end # module

max_stars_repo_name: emmt/ConditionallyOptimize


text:
using Distributions
using Turing
using Test

obs12 = [1,2,1,2,2,2,2,2,2,2]

@model constrained_simplex_test(obs12) = begin
  ps ~ Dirichlet(2, 3)
  for i = 1:length(obs12)
    obs12[i] ~ Categorical(ps)
  end
  return ps
end

chain = sample(constrained_simplex_test(obs12), HMC(1000, 0.75, 2))

check_numerical(chain, [:ps], [[5/16; 11/16]], eps=0.015)


text:
mutable struct NodeManager
    masterID::Union{ Missing, Int }
    localmeta::Union{ Missing, Dict{Int64, MrPhelps.WorkerMetaData} }
    computemeta::Union{ Missing, Dict{Int64, MrPhelps.WorkerMetaData} }
    machinenodemap::Union{ Missing, Dict{String, Vector{Int64}} }
end

"""
    NodeManager()

Constructs a NodeManager objects. This object contains metadata about the processes
available to Distributed.jl.

"""
function NodeManager()
    if Distributed.nworkers() == 1
        @info   "NodeManager instance created with only a local worker and no Distributed.jl workers. \n " *
                "You can attach workers later, and update the NodeManager with" *
                " `update!(nm::NodeManager)`."
    end
    local_metadata, compute_metadata = MrPhelps.worker_meta()
    return NodeManager( first( [ k for (k,v) in local_metadata ] ),
                        local_metadata,
                        compute_metadata,
                        availablemachines( compute_metadata ))
end

"""
    update!(nm::NodeManager)

Updates a NodeManager instance `nm` via the processes present in Distributed.jl .

"""
function update!(nm::NodeManager)
    if Distributed.nworkers() == 1
        @info   "NodeManager instance created with only a local worker and no Distributed.jl workers. \n " *
                "You can attach workers later, and update the NodeManager with" *
                " `update!(nm::NodeManager)`."
    end
    local_metadata, compute_metadata = MrPhelps.worker_meta()
    nm.masterID         = first( [ k for (k,v) in local_metadata ] )
    nm.localmeta        = local_metadata
    nm.computemeta      = compute_metadata
    nm.machinenodemap   = availablemachines(compute_metadata)
    return nm
end

"""
    availablemachines( computemetadict::Dict )

Takes a dictionary of process ids mapped to WorkerMetaData(`computemetadict`) as input and returns a
dictionary of whose keys are unique hosts and whose values are vectors of their Process ID numbers.

"""
function availablemachines( computemetadict::Dict )
    addresses = unique( [ worker_meta.address for ( id, worker_meta ) in computemetadict ] )
    groupedby = Dict()
    for addr in addresses
        groupedby[ addr ] = [ k for (k, meta) in computemetadict if ( meta.address == addr ) ]
    end
    return groupedby
end

availablemachines( computemetadict::Missing ) = missing

"""
    machinenames( nm::NodeManager )

Takes a NodeManager(`nm`) as input and returns a vector of available machine names.

"""
machinenames( nm::NodeManager ) = return [k for k in keys(nm.machinenodemap)]

"""
    worker_count(nm::NodeManager, machine::String)

Retrieves the number of workers allocated to a specific machine.

"""
worker_count(nm::NodeManager, machine::String) = length( nm.machinenodemap[ machine ] )

"""
    total_worker_counts(nm::NodeManager, machines::Vector{String} = machinenames( nm ) )

Retrieves the total number of workers in either a NodeManager instance, or in a vector of machines.

"""
function total_worker_counts(nm::NodeManager, machines::Vector{String} = machinenames( nm ) )
    return sum( [ length( nm.machinenodemap[ machine ] ) for machine in machines ] )
end

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps


text:
module Database

using Revise
using YAML, DataFrames, Logging
using SearchLight, SearchLight.Configuration

const ADAPTERS_FOLDER = joinpath(@__DIR__, "database_adapters")

function setup_adapter(adapter = SearchLight.config.db_config_settings["adapter"] * "DatabaseAdapter") :: Union{Bool,Nothing}
  Core.eval(@__MODULE__, Meta.parse("""include(joinpath(ADAPTERS_FOLDER, "$adapter.jl"))"""))
  Core.eval(@__MODULE__, Meta.parse("using .$adapter"))

  Core.eval(@__MODULE__, :(db_adapter = Symbol($adapter)))
  ! Core.isdefined(@__MODULE__, :DatabaseAdapter) && Core.eval(@__MODULE__, :(const DatabaseAdapter = $db_adapter))
  ! Core.isdefined(@__MODULE__, :DatabaseHandle) && Core.eval(@__MODULE__, :(const DatabaseHandle = $db_adapter.DatabaseHandle))
  ! Core.isdefined(@__MODULE__, :ResultHandle) && Core.eval(@__MODULE__, :(const ResultHandle = $db_adapter.ResultHandle))

  Core.eval(@__MODULE__, :(export DatabaseAdapter))

  true
end

"""
    connect()::DatabaseHandle
    connect(conn_settings::Dict{String,Any})::DatabaseHandle

Connects to the DB and returns a database handler. If used without arguments, it defaults to using `config.db_config_settings`

# Examples
```julia
julia> Database.connect()
PostgreSQL.PostgresDatabaseHandle(Ptr{Nothing} @0x00007fbf3839f360,0x00000000,false)

julia> dict = config.db_config_settings
Dict{String,Any} with 6 entries:
  "host"     => "localhost"
  "password" => "adrian"
  "username" => "adrian"
  "port"     => 5432
  "database" => "blogjl_dev"
  "adapter"  => "PostgreSQL"

julia> Database.connect(dict)
PostgreSQL.PostgresDatabaseHandle(Ptr{Nothing} @0x00007fbf3839f360,0x00000000,false)
```
"""
@inline function connect() #::DatabaseHandle
  connect(SearchLight.config.db_config_settings)
end
function connect(conn_settings::Dict)
  isdefined(@__MODULE__, :DatabaseAdapter) || connect!(conn_settings)

  c = Base.invokelatest(DatabaseAdapter.connect, conn_settings) #::DatabaseHandle
  Core.eval(@__MODULE__, :(CONNECTION = $c))

  c
end

@inline function connect!(conn_settings::Dict)
  SearchLight.config.db_config_settings["adapter"] = conn_settings["adapter"]
  setup_adapter() && Database.connect(conn_settings) #::DatabaseHandle
end

"""

"""
@inline function disconnect(conn)
  DatabaseAdapter.disconnect(conn)
end
@inline function disconnect()
  DatabaseAdapter.disconnect(CONNECTION)
end

"""
    create_database()::Bool
    create_database(db_name::String)::Bool

Invokes the database adapter's create database method. If invoked without param, it defaults to the
database name defined in `config.db_config_settings`
"""
@inline function create_database() :: Bool
  create_database(SearchLight.config.db_config_settings["database"])
end
@inline function create_database(db_name::String)::Bool
  DatabaseAdapter.create_database(db_name)
end

"""
    create_migrations_table(table_name::String)::Bool

Invokes the database adapter's create migrations table method. If invoked without param, it defaults to the
database name defined in `config.db_migrations_table_name`
"""
@inline function create_migrations_table(table_name::String) :: Bool
  DatabaseAdapter.create_migrations_table(table_name)
end

"""
    init()::Bool

Sets up the DB tables used by SearchLight.
"""
@inline function init() :: Bool
  DatabaseAdapter.create_migrations_table(SearchLight.config.db_migrations_table_name)
end

@inline function escape_column_name(c::String)
  result =  try
              DatabaseAdapter.escape_column_name(c, CONNECTION)::String
            catch ex
              @error ex
            end

  result
end

@inline function escape_value(v::T)::T where {T}
  result =  try
              DatabaseAdapter.escape_value(v, CONNECTION)
            catch ex
              @error ex
            end

  result
end

@inline function table_columns(table_name::String) :: DataFrames.DataFrame
  query(DatabaseAdapter.table_columns_sql(table_name), suppress_output = true)
end

"""
    query(sql::String; suppress_output::Bool = false)::DataFrames.DataFrame

Executes the `sql` query against the database adapter and returns a DataFrame result.
Optionally logs the result DataFrame.

# Examples:
```julia
julia> SearchLight.to_fetch_sql(Article, SQLQuery(limit = 5)) |> Database.query;

2017-01-16T21:33:40.079 - info: SQL QUERY: SELECT \"articles\".\"id\" AS \"articles_id\", \"articles\".\"title\" AS \"articles_title\", \"articles\".\"summary\" AS \"articles_summary\", \"articles\".\"content\" AS \"articles_content\", \"articles\".\"updated_at\" AS \"articles_updated_at\", \"articles\".\"published_at\" AS \"articles_published_at\", \"articles\".\"slug\" AS \"articles_slug\" FROM \"articles\" LIMIT 5

  0.001172 seconds (16 allocations: 576 bytes)

2017-01-16T21:33:40.089 - info: 5×7 DataFrames.DataFrame
...
```
"""
@inline function query(sql::String; suppress_output::Bool = false, system_query::Bool = false) :: DataFrames.DataFrame
  df::DataFrames.DataFrame =  DatabaseAdapter.query(sql, (suppress_output || system_query), CONNECTION)
  (! suppress_output && ! system_query && SearchLight.config.log_db) && @info(df)

  df
end

"""

"""
@inline function relation_to_sql(m::T, rel::Tuple{SQLRelation,Symbol})::String where {T<:AbstractModel}
  DatabaseAdapter.relation_to_sql(m, rel)
end

"""

"""
@inline function to_find_sql(m::Type{T}, q::SQLQuery, joins::Vector{SQLJoin{N}})::String where {T<:AbstractModel, N<:AbstractModel}
  DatabaseAdapter.to_find_sql(m, q, joins)
end
@inline function to_find_sql(m::Type{T}, q::SQLQuery)::String where {T<:AbstractModel}
  DatabaseAdapter.to_find_sql(m, q)
end

const to_fetch_sql = to_find_sql

"""

"""
@inline function to_store_sql(m::T; conflict_strategy = :error)::String where {T<:AbstractModel} # upsert strateygy = :none | :error | :ignore | :update
  DatabaseAdapter.to_store_sql(m, conflict_strategy = conflict_strategy)
end

"""

"""
@inline function delete_all(m::Type{T}; truncate::Bool = true, reset_sequence::Bool = true, cascade::Bool = false)::Nothing where {T<:AbstractModel}
  DatabaseAdapter.delete_all(m, truncate = truncate, reset_sequence = reset_sequence, cascade = cascade)
end

"""

"""
@inline function delete(m::T)::T where {T<:AbstractModel}
  DatabaseAdapter.delete(m)
end

"""

"""
@inline function count(m::Type{T}, q::SQLQuery = SQLQuery())::Int where {T<:AbstractModel}
  DatabaseAdapter.count(m, q)
end

"""

"""
@inline function update_query_part(m::T)::String where {T<:AbstractModel}
  DatabaseAdapter.update_query_part(m)
end

"""

"""
@inline function to_select_part(m::Type{T}, cols::Vector{SQLColumn}, joins = SQLJoin[])::String where {T<:AbstractModel}
  DatabaseAdapter.to_select_part(m, cols, joins)
end
"""

"""
function _to_select_part(m::Type{T}, cols::Vector{SQLColumn}, joins = SQLJoin[])::String where {T<:AbstractModel}
  _m::T = m()

  joined_tables = []

  if has_relation(_m, RELATION_HAS_ONE)
    rels = _m.has_one
    joined_tables = vcat(joined_tables, map(x -> is_lazy(x) ? nothing : (x.model_name)(), rels))
  end

  if has_relation(_m, RELATION_HAS_MANY)
    rels = _m.has_many
    joined_tables = vcat(joined_tables, map(x -> is_lazy(x) ? nothing : (x.model_name)(), rels))
  end

  if has_relation(_m, RELATION_BELONGS_TO)
    rels = _m.belongs_to
    joined_tables = vcat(joined_tables, map(x -> is_lazy(x) ? nothing : (x.model_name)(), rels))
  end

  filter!(x -> x != nothing, joined_tables)

  if ! isempty(cols)
    table_columns = []
    cols = vcat(cols, columns_from_joins(joins))

    for column in cols
      push!(table_columns, prepare_column_name(column, _m))
    end

    return join(table_columns, ", ")
  else
    table_columns = join(to_fully_qualified_sql_column_names(_m, persistable_fields(_m), escape_columns = true), ", ")
    table_columns = isempty(table_columns) ? String[] : vcat(table_columns, map(x -> prepare_column_name(x, _m), columns_from_joins(joins)))

    related_table_columns = String[]
    for rels in map(x -> to_fully_qualified_sql_column_names(x, persistable_fields(x), escape_columns = true), joined_tables)
      for col in rels
        push!(related_table_columns, col)
      end
    end

    return join([table_columns ; related_table_columns], ", ")
  end
end

"""

"""
@inline function to_from_part(m::Type{T})::String where {T<:AbstractModel}
  DatabaseAdapter.to_from_part(m)
end

"""

"""
@inline function to_where_part(m::Type{T}, w::Vector{SQLWhereEntity}, scopes::Vector{Symbol})::String where {T<:AbstractModel}
  DatabaseAdapter.to_where_part(m, w, scopes)
end
@inline function to_where_part(w::Vector{SQLWhereEntity})::String
  DatabaseAdapter.to_where_part(w)
end

"""

"""
@inline function required_scopes(m::Type{T})::Vector{SQLWhereEntity} where {T<:AbstractModel}
  s = scopes(m)
  haskey(s, :required) ? s[:required] : SQLWhereEntity[]
end

"""

"""
@inline function scopes(m::Type{T})::Dict{Symbol,Vector{SQLWhereEntity}} where {T<:AbstractModel}
  in(:scopes, fieldnames(m)) ? getfield(m()::T, :scopes) :  Dict{Symbol,Vector{SQLWhereEntity}}()
end

"""

"""
@inline function to_order_part(m::Type{T}, o::Vector{SQLOrder})::String where {T<:AbstractModel}
  DatabaseAdapter.to_order_part(m, o)
end

"""

"""
@inline function to_group_part(g::Vector{SQLColumn}) :: String
  DatabaseAdapter.to_group_part(g)
end

"""

"""
@inline function to_limit_part(l::SQLLimit) :: String
  DatabaseAdapter.to_limit_part(l)
end

"""

"""
@inline function to_offset_part(o::Int) :: String
  DatabaseAdapter.to_offset_part(o)
end

"""

"""
@inline function to_having_part(h::Vector{SQLHaving}) :: String
  DatabaseAdapter.to_having_part(h)
end

"""

"""
@inline function to_join_part(m::Type{T}, joins = SQLJoin[])::String where {T<:AbstractModel}
  DatabaseAdapter.to_join_part(m, joins)
end

"""
  columns_from_joins(joins::Vector{SQLJoin})::Vector{SQLColumn}

Extracts columns from joins param and adds to be used for the SELECT part
"""
function columns_from_joins(joins::Vector{SQLJoin})::Vector{SQLColumn}
  jcols = SQLColumn[]
  for j in joins
    jcols = vcat(jcols, j.columns)
  end

  jcols
end

"""

"""
function prepare_column_name(column::SQLColumn, _m::T)::String where {T<:AbstractModel}
  if column.raw
    column.value |> string
  else
    column_data::Dict{Symbol,Any} = SearchLight.from_literal_column_name(column.value)
    if ! haskey(column_data, :table_name)
      column_data[:table_name] = table_name(_m)
    end
    if ! haskey(column_data, :alias)
      column_data[:alias] = ""
    end

    DatabaseAdapter.column_data_to_column_name(column, column_data)
  end
end

"""

"""
@inline function rand(m::Type{T}; limit = 1)::Vector{T} where {T<:AbstractModel}
  DatabaseAdapter.rand(m, limit = limit)
end
@inline function rand(m::Type{T}, scopes::Vector{Symbol}; limit = 1)::Vector{T} where {T<:AbstractModel}
  DatabaseAdapter.rand(m, scopes, limit = limit)
end

"""

"""
@inline function index_name(table_name::Union{String,Symbol}, column_name::Union{String,Symbol}) :: String
  string(table_name) * "__" * "idx_" * string(column_name)
end

end


text:
# Check that interfaces of Evaluators are well implemented
@testset "API of evaluator $Evaluator" for Evaluator in [
    ExaPF.ReducedSpaceEvaluator,
    ExaPF.AugLagEvaluator,
    ExaPF.ProxALEvaluator,
    ExaPF.SlackEvaluator,
    ExaPF.FeasibilityEvaluator,
]
    datafile = joinpath(INSTANCES_DIR, "case9.m")
    # Default constructor: should take as input path to instance
    nlp = Evaluator(datafile)

    # Test printing
    println(devnull, nlp)

    # Test consistence
    n = ExaPF.n_variables(nlp)
    m = ExaPF.n_constraints(nlp)
    u = ExaPF.initial(nlp)

    u_min, u_max = ExaPF.bounds(nlp, ExaPF.Variables())
    g_min, g_max = ExaPF.bounds(nlp, ExaPF.Constraints())
    buffer = get(nlp, ExaPF.PhysicalState())

    @testset "Evaluator's API" begin
        @test n == length(u)
        @test length(u_min) == length(u_max) == n
        @test u_min <= u_max
        @test length(g_min) == length(g_max) == m
        if m > 0
            @test g_min <= g_max
        end

        @test isa(get(nlp, ExaPF.Constraints()), Array{Function})
        @test isa(get(nlp, State()), AbstractVector)
        @test isa(buffer, ExaPF.AbstractBuffer)
        @test ExaPF.constraints_type(nlp) in [:bound, :equality, :inequality]

        @test isa(ExaPF.has_hessian(nlp), Bool)
        @test isa(ExaPF.has_hessian_lagrangian(nlp), Bool)

        # setters
        nbus = get(nlp, PS.NumberOfBuses())
        loads = similar(u, nbus) ; fill!(loads, 1)
        ExaPF.setvalues!(nlp, PS.ActiveLoad(), loads)
        ExaPF.setvalues!(nlp, PS.ReactiveLoad(), loads)

        ExaPF.reset!(nlp)
    end

    @testset "Evaluator's callbacks" begin
        # 1/ update! function
        conv = ExaPF.update!(nlp, u)
        @test isa(conv, ExaPF.ConvergenceStatus)
        @test conv.has_converged

        # 2/ objective function
        c = ExaPF.objective(nlp, u)
        @test isa(c, Real)

        # 3/ gradient! function
        function reduced_cost(u_)
            ExaPF.update!(nlp, u_)
            return ExaPF.objective(nlp, u_)
        end
        g = similar(u) ; fill!(g, 0)
        ExaPF.gradient!(nlp, g, u)
        grad_fd = FiniteDiff.finite_difference_gradient(reduced_cost, u)
        @test isapprox(g, grad_fd, rtol=1e-5)

        # 4/ Constraint
        ## Evaluation of the constraints
        if m > 0
            cons = similar(g_min) ; fill!(cons, 0.0)
            ExaPF.constraint!(nlp, cons, u)

            # Vector product
            jv = similar(u_min) ; fill!(jv, 0.0)
            v = similar(g_min) ; fill!(v, 1.0)
            ExaPF.jtprod!(nlp, jv, u, v)
            function reduced_cons(u_)
                ExaPF.update!(nlp, u_)
                ExaPF.constraint!(nlp, cons, u_)
                return dot(v, cons)
            end
            jv_fd = FiniteDiff.finite_difference_gradient(reduced_cons, u)

            @test isapprox(jv, jv_fd, rtol=1e-5)

            # Jacobian
            J =  ExaPF.jacobian(nlp, u)
            @test J' * v ≈ jv
        end
    end
end


text:
using Genie.Router, Genie.WebChannels
using Lib

route("/") do
  serve_static_file("welcome.html")
end

route("/error500") do
  error_500("Something went wrong")
end

route("/error404") do
  error_404("the page you want")
end

route("/mychannel/subscribe") do
  serve_static_file("test.html")
end

channel("/mychannel/subscribe") do
  WebChannels.subscribe(@params(:WS_CLIENT), "mychannel")
  @show "Subscription OK"
end

channel("/mychannel/ping") do
  WebChannels.message(@params(:WS_CLIENT), "pong")

  msg = Lib.test()
  WebChannels.message(@params(:WS_CLIENT), msg)

  @show "pong OK"
end

function do_your_stuff_here()
  # push updates
  #WebChannels.message("mychannel", json_payload)

  WebChannels.message("mychannel", "test")
end


text:
using Latexify
using DiffEqBiological
using ParameterizedFunctions
using Markdown
using Test

#inline
@test latexify(:(x * y); env=:inline, cdot=false) == raw"$x y$"

@test latexify(:(x * y); env=:inline, cdot=true) == raw"$x \cdot y$"

@test latexify(:(x*(y+z)*y*(z+a)*(z+b)); env=:inline, cdot=false) ==
raw"$x \left( y + z \right) y \left( z + a \right) \left( z + b \right)$"

@test latexify(:(x*(y+z)*y*(z+a)*(z+b)); env=:inline, cdot=true) ==
raw"$x \cdot \left( y + z \right) \cdot y \cdot \left( z + a \right) \cdot \left( z + b \right)$"

# raw
@test latexify(:(x * y); env=:raw, cdot=false) == raw"x y"

@test latexify(:(x * y); env=:raw, cdot=true) == raw"x \cdot y"

@test latexify(:(x * (y + z) * y * (z + a) * (z + b)); env=:raw, cdot=false) ==
raw"x \left( y + z \right) y \left( z + a \right) \left( z + b \right)"

@test latexify(:(x * (y + z) * y * (z + a) * (z + b)); env=:raw, cdot=true) ==
raw"x \cdot \left( y + z \right) \cdot y \cdot \left( z + a \right) \cdot \left( z + b \right)"

# array
@test latexify( [:(x*y), :(x*(y+z)*y*(z+a)*(z+b))]; env=:array, transpose=true, cdot=false) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
x y & x \left( y + z \right) y \left( z + a \right) \left( z + b \right) \\
\end{array}
\right]
\end{equation}
"

@test latexify( [:(x*y), :(x*(y+z)*y*(z+a)*(z+b))]; env=:array, transpose=true, cdot=true) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
x \cdot y & x \cdot \left( y + z \right) \cdot y \cdot \left( z + a \right) \cdot \left( z + b \right) \\
\end{array}
\right]
\end{equation}
"

# align
f = @ode_def TestAlignFeedback2 begin
    dx = y*c_1 - x
    dy = x^c_2 - y
end c_1 c_2

@test latexify(f; env=:align, cdot=false) ==
raw"\begin{align}
\frac{dx}{dt} =& y c_{1} - x \\
\frac{dy}{dt} =& x^{c_{2}} - y
\end{align}
"

@test latexify(f; env=:align, cdot=true) ==
raw"\begin{align}
\frac{dx}{dt} =& y \cdot c_{1} - x \\
\frac{dy}{dt} =& x^{c_{2}} - y
\end{align}
"

# reactions
if VERSION >= v"1.0"

rn = @reaction_network begin
    hillr(D₂,α,K,n), ∅ --> m₁
    hillr(D₁,α,K,n), ∅ --> m₂
    (δ,γ), m₁ ↔ ∅
    (δ,γ), m₂ ↔ ∅
    β, m₁ --> m₁ + P₁
    β, m₂ --> m₂ + P₂
    μ, P₁ --> ∅
    μ, P₂ --> ∅
    (k₊,k₋), 2P₁ ↔ D₁
    (k₊,k₋), 2P₂ ↔ D₂
    (k₊,k₋), P₁+P₂ ↔ T
end α K n δ γ β μ k₊ k₋;

@test_broken latexify(rn; cdot=false) ==
raw"\begin{align}
\frac{dm_1}{dt} =& \frac{\alpha K^{n}}{K^{n} + D_2^{n}} - \delta m_1 + \gamma \\
\frac{dm_2}{dt} =& \frac{\alpha K^{n}}{K^{n} + D_1^{n}} - \delta m_2 + \gamma \\
\frac{dP_1}{dt} =& \beta m_1 - \mu P_1 -2 k_+ \frac{P_1^{2}}{2} + 2 k_- D_1 - k_+ P_1 P_2 + k_- T \\
\frac{dP_2}{dt} =& \beta m_2 - \mu P_2 -2 k_+ \frac{P_2^{2}}{2} + 2 k_- D_2 - k_+ P_1 P_2 + k_- T \\
\frac{dD_1}{dt} =& k_+ \frac{P_1^{2}}{2} - k_- D_1 \\
\frac{dD_2}{dt} =& k_+ \frac{P_2^{2}}{2} - k_- D_2 \\
\frac{dT}{dt} =& k_+ P_1 P_2 - k_- T
\end{align}
"

@test_broken latexify(rn; cdot=true) ==
raw"\begin{align}
\frac{dm_1}{dt} =& \frac{\alpha \cdot K^{n}}{K^{n} + D_2^{n}} - \delta \cdot m_1 + \gamma \\
\frac{dm_2}{dt} =& \frac{\alpha \cdot K^{n}}{K^{n} + D_1^{n}} - \delta \cdot m_2 + \gamma \\
\frac{dP_1}{dt} =& \beta \cdot m_1 - \mu \cdot P_1 -2 \cdot k_+ \cdot \frac{P_1^{2}}{2} + 2 \cdot k_- \cdot D_1 - k_+ \cdot P_1 \cdot P_2 + k_- \cdot T \\
\frac{dP_2}{dt} =& \beta \cdot m_2 - \mu \cdot P_2 -2 \cdot k_+ \cdot \frac{P_2^{2}}{2} + 2 \cdot k_- \cdot D_2 - k_+ \cdot P_1 \cdot P_2 + k_- \cdot T \\
\frac{dD_1}{dt} =& k_+ \cdot \frac{P_1^{2}}{2} - k_- \cdot D_1 \\
\frac{dD_2}{dt} =& k_+ \cdot \frac{P_2^{2}}{2} - k_- \cdot D_2 \\
\frac{dT}{dt} =& k_+ \cdot P_1 \cdot P_2 - k_- \cdot T
\end{align}
"

@test latexify(rn; env=:chemical, cdot=false) ==
raw"\begin{align}
\require{mhchem}
\ce{ \varnothing &->[\frac{\alpha K^{n}}{K^{n} + D_2^{n}}] m_{1}}\\
\ce{ \varnothing &->[\frac{\alpha K^{n}}{K^{n} + D_1^{n}}] m_{2}}\\
\ce{ m_{1} &<=>[\delta][\gamma] \varnothing}\\
\ce{ m_{2} &<=>[\delta][\gamma] \varnothing}\\
\ce{ m_{1} &->[\beta] m_{1} + P_{1}}\\
\ce{ m_{2} &->[\beta] m_{2} + P_{2}}\\
\ce{ P_{1} &->[\mu] \varnothing}\\
\ce{ P_{2} &->[\mu] \varnothing}\\
\ce{ 2 P_1 &<=>[k_{+}][k_{-}] D_{1}}\\
\ce{ 2 P_2 &<=>[k_{+}][k_{-}] D_{2}}\\
\ce{ P_{1} + P_{2} &<=>[k_{+}][k_{-}] T}
\end{align}
"

@test latexify(rn; env=:chemical, cdot=true) ==
raw"\begin{align}
\require{mhchem}
\ce{ \varnothing &->[\frac{\alpha \cdot K^{n}}{K^{n} + D_2^{n}}] m_{1}}\\
\ce{ \varnothing &->[\frac{\alpha \cdot K^{n}}{K^{n} + D_1^{n}}] m_{2}}\\
\ce{ m_{1} &<=>[\delta][\gamma] \varnothing}\\
\ce{ m_{2} &<=>[\delta][\gamma] \varnothing}\\
\ce{ m_{1} &->[\beta] m_{1} + P_{1}}\\
\ce{ m_{2} &->[\beta] m_{2} + P_{2}}\\
\ce{ P_{1} &->[\mu] \varnothing}\\
\ce{ P_{2} &->[\mu] \varnothing}\\
\ce{ 2 \cdot P_1 &<=>[k_{+}][k_{-}] D_{1}}\\
\ce{ 2 \cdot P_2 &<=>[k_{+}][k_{-}] D_{2}}\\
\ce{ P_{1} + P_{2} &<=>[k_{+}][k_{-}] T}
\end{align}
"

# tabular
@test_broken latexify(rn.symjac; env=:tabular, cdot=false) ==
raw"\begin{tabular}{ccccccc}
$ - \delta$ & $0$ & $0$ & $0$ & $0$ & $\frac{ - K^{n} n \alpha D_2^{-1 + n}}{\left( K^{n} + D_2^{n} \right)^{2}}$ & $0$\\
$0$ & $ - \delta$ & $0$ & $0$ & $\frac{ - K^{n} n \alpha D_1^{-1 + n}}{\left( K^{n} + D_1^{n} \right)^{2}}$ & $0$ & $0$\\
$\beta$ & $0$ & $ - \mu - 2 k_+ P_1 - k_+ P_2$ & $ - k_+ P_1$ & $2 k_-$ & $0$ & $k_{-}$\\
$0$ & $\beta$ & $ - k_+ P_2$ & $ - \mu - k_+ P_1 - 2 k_+ P_2$ & $0$ & $2 k_-$ & $k_{-}$\\
$0$ & $0$ & $k_+ P_1$ & $0$ & $ - k_-$ & $0$ & $0$\\
$0$ & $0$ & $0$ & $k_+ P_2$ & $0$ & $ - k_-$ & $0$\\
$0$ & $0$ & $k_+ P_2$ & $k_+ P_1$ & $0$ & $0$ & $ - k_-$\\
\end{tabular}
"

@test_broken latexify(rn.symjac; env=:tabular, cdot=true) ==
raw"\begin{tabular}{ccccccc}
$ - \delta$ & $0$ & $0$ & $0$ & $0$ & $\frac{ - K^{n} \cdot n \cdot \alpha \cdot D_2^{-1 + n}}{\left( K^{n} + D_2^{n} \right)^{2}}$ & $0$\\
$0$ & $ - \delta$ & $0$ & $0$ & $\frac{ - K^{n} \cdot n \cdot \alpha \cdot D_1^{-1 + n}}{\left( K^{n} + D_1^{n} \right)^{2}}$ & $0$ & $0$\\
$\beta$ & $0$ & $ - \mu - 2 \cdot k_+ \cdot P_1 - k_+ \cdot P_2$ & $ - k_+ \cdot P_1$ & $2 \cdot k_-$ & $0$ & $k_{-}$\\
$0$ & $\beta$ & $ - k_+ \cdot P_2$ & $ - \mu - k_+ \cdot P_1 - 2 \cdot k_+ \cdot P_2$ & $0$ & $2 \cdot k_-$ & $k_{-}$\\
$0$ & $0$ & $k_+ \cdot P_1$ & $0$ & $ - k_-$ & $0$ & $0$\\
$0$ & $0$ & $0$ & $k_+ \cdot P_2$ & $0$ & $ - k_-$ & $0$\\
$0$ & $0$ & $k_+ \cdot P_2$ & $k_+ \cdot P_1$ & $0$ & $0$ & $ - k_-$\\
\end{tabular}
"
end

# mdtable
arr = ["x*(y-1)", 1.0, 3*2, :(x-2y), :symb]

@test latexify(arr; env=:mdtable, cdot=false) ==
Markdown.md"| $x \left( y - 1 \right)$ |
| ------------------------:|
|                    $1.0$ |
|                      $6$ |
|                $x - 2 y$ |
|                   $symb$ |
"

@test latexify(arr; env=:mdtable, cdot=true) ==
Markdown.md"| $x \cdot \left( y - 1 \right)$ |
| ------------------------------:|
|                          $1.0$ |
|                            $6$ |
|                $x - 2 \cdot y$ |
|                         $symb$ |
"


text:
@inline _sortperm(X::AbstractVecOrMat) = size(X,2)==1 ? sortperm(ndims(X)==1 ? X : vec(X), alg=RadixSort) : sortperm(collect(eachrow(X)))  # sort a data matrix

function Init!(o::StrBootTest{T}) where T  # for efficiency when varying r repeatedly to make CI, do stuff once that doesn't depend on r
	o.kX = o.kX₁ + o.kX₂
	o.kX₂==0 && (o.X₂ = zeros(T,o.Nobs,0))
  iszero(o.kY₂) && (o.Y₂ = zeros(T,o.Nobs,0))
  o.kZ = o.kX₁ + o.kY₂
  if o.liml && o.kX₂==o.kY₂  # exactly identified liml = 2SLS
  	o.κ = one(T)
  	o.liml = false
  end
  if !(o.REst = length(o.R₁)>0)  # base model contains no restrictions?
    o.R₁ = zeros(T,0,o.kZ)
    o.r₁ = zeros(T,0)
  end
  isnan(o.κ) && (o.κ = o.kX₂>0 ? one(T) : zero(T))  # if κ in κ-class estimation not specified, it's 0 or 1 for OLS or 2SLS
  o.WRE = !(iszero(o.κ) || o.scorebs) || o.arubin

  iszero(o.B) && (o.scorebs = true)

  o.haswt = !iszero(length(o.wt))
	if o.haswt
		o.sumwt = sum(o.wt)
		o._Nobs = o.fweights ? o.sumwt : T(o.Nobs)
	else
		o.sumwt = one(T)
		o._Nobs = T(o.Nobs)
		o.sqrtwt = T[]
	end

	o.subcluster = o.NClustVar - o.NErrClustVar

	if !(iszero(o.NClustVar) || o.issorted)
		p = _sortperm(view(o.ID, :, [collect(o.subcluster+1:o.NClustVar); collect(1:o.subcluster)]))  # sort by err cluster vars, then remaining boot cluster vars
		o.ID = ndims(o.ID)==1 ? o.ID[p] : o.ID[p,:]
		o.X₁ = ndims(o.X₁)==1 ? o.X₁[p] : o.X₁[p,:]
		o.X₂ = ndims(o.X₂)==1 ? o.X₂[p] : o.X₂[p,:]
		o.y₁ = o.y₁[p]
		o.Y₂ = ndims(o.Y₂)==1 ? o.Y₂[p] : o.Y₂[p,:]
		isdefined(o, :FEID) && nrows(o.FEID)>0 && (o.FEID = o.FEID[p])

		if o.haswt
			o.wt = o.wt[p]
			o.sqrtwt = sqrt.(o.wt)
			o.y₁ .*= o.sqrtwt  # can overwrite sorted copy of user's data
			length(o.Y₂)>0 && (o.Y₂ .*= o.sqrtwt)
			length(o.X₁)>0 && (o.X₁ .*= o.sqrtwt)
			length(o.X₂)>0 && (o.X₂ .*= o.sqrtwt)
		end
	elseif o.haswt
		o.sqrtwt = sqrt.(o.wt)
		o.y₁ = o.y₁ .* o.sqrtwt  # don't overwrite user's data
		length(o.Y₂)>0 && (o.Y₂ = o.Y₂ .* o.sqrtwt)
		length(o.X₁)>0 && (o.X₁ = o.X₁ .* o.sqrtwt)
		length(o.X₂)>0 && (o.X₂ = o.X₂ .* o.sqrtwt)
	end

  if o.WREnonARubin
	  if iszero(o.NClustVar)
	    o.info⋂ = o.info✻ = [i:i for i in 1:o.Nobs]  # no clustering, so no collapsing by cluster
	  else
	    o.info✻, o.ID✻ = panelsetupID(o.ID, 1:o.NBootClustVar)

	  end
  elseif iszero(o.NClustVar)
	  o.info✻ = [i:i for i in 1:o.Nobs]  # causes no collapsing of data in panelsum() calls
  else
	  o.info✻ = panelsetup(o.ID, 1:min(o.NClustVar,o.NBootClustVar))  # bootstrap cluster grouping defs rel to original data
  end
	o.N✻ = nrows(o.info✻)

	if o.NClustVar > o.NBootClustVar  # info for grouping by intersections of all bootstrap && clustering vars wrt data; used to speed crosstab UXAR wrt bootstrapping cluster && intersection of all error clusters
		o.info✻⋂ = panelsetup(o.ID, 1:o.NClustVar)
		if o.subcluster>0 && nrows(o.info✻) ≠ nrows(o.info✻⋂)
			throw(ErrorException("\nThis program can only perform the subcluster bootstrap when the bootstrap clusters are nested within the (intersections of the) error clusters.\n"))
		end
	else
		o.info✻⋂ = o.info✻  # info for grouping by intersections of all bootstrap && clustering vars wrt data; used to speed crosstab UXAR wrt bootstrapping cluster && intersection of all error clusters
		o.WREnonARubin && (o._ID✻⋂ = o.ID✻)
	end
	o.N✻⋂ = nrows(o.info✻⋂)

	if o.NClustVar > o.NErrClustVar  # info for intersections of error clustering wrt data
		o.info⋂ = panelsetup(o.ID, o.subcluster+1:o.NClustVar)
		ID⋂_✻⋂ = length(o.info⋂)==o.Nobs ? o.ID : o.ID[first.(o.info⋂ ),:]  # version of ID matrix with one row for each all-error-cluster-var intersection instead of 1 row for each obs; gets resorted
		o.ID✻⋂ = o.N✻⋂==o.Nobs ? o.ID : o.ID[first.(o.info✻⋂),:]  # version of ID matrix with one row for each all-bootstrap && error cluster-var intersection instead of 1 row for each obs
	else
		o.info⋂ = o.info✻⋂  # info for intersections of error clustering wrt data
		o.ID✻⋂ = ID⋂_✻⋂ = nrows(o.info⋂)==o.Nobs ? o.ID : o.ID[first.(o.info⋂),:]  # version of ID matrix with one row for each all-error-cluster-var intersection instead of 1 row for each obs; gets resorted
	end
	o.WREnonARubin && (o.info✻_✻⋂ = iszero(o.NBootClustVar) ? o.info✻⋂ : panelsetup(o.ID✻⋂, 1:o.NBootClustVar))
	o.N⋂ = nrows(o.info⋂)

	if o.bootstrapt
    if o.NClustVar>0
	    minN = T(Inf)

	    combs = [x & 2^y > 0 for x in 2^o.NErrClustVar-1:-1:1, y in o.NErrClustVar-1:-1:0]  # represent all error clustering combinations. First is intersection of all error clustering vars
	    o.clust = Vector{StrClust{T}}(undef, nrows(combs))  # leave out no-cluster combination
	    o.NErrClustCombs = length(o.clust)

	    o.BootClust = 2^(o.NClustVar - o.NBootClustVar)  # location of bootstrap clustering within list of cluster combinations

	    @inbounds for c ∈ 1:o.NErrClustCombs  # for each error clustering combination
		    ClustCols = o.subcluster .+ findall(@view combs[c,:])
		    even = isodd(length(ClustCols))  # not a typo

		    if isone(c)
		      if iszero(o.subcluster)
		    	  order = Vector{Int64}(undef,0)
		    	  info  = Vector{UnitRange{Int64}}(undef, 0)  # causes no collapsing of data in panelsum() calls
						N = o.N✻⋂
		      else
		    	  order = _sortperm(@view ID⋂_✻⋂[:,ClustCols])
		    	  info  = panelsetup(view(ID⋂_✻⋂,order,:), ClustCols)
						N = nrows(info)
		      end
		    else
		      if any(combs[c, minimum(findall(combs[c,:] .≠ combs[c-1,:])):end])  # if this sort ordering same as last to some point and missing thereafter, no need to re-sort
		    	  order = _sortperm(@view ID⋂_✻⋂[:,ClustCols])
		    	  info = panelsetup(view(ID⋂_✻⋂,order,:), ClustCols)
		      else
		    	  order = Vector{Int64}(undef,0)
						info = panelsetup(ID⋂_✻⋂, ClustCols)
		      end
					N = nrows(info)
		    end

		    minN = min(minN,N)
				multiplier = o.clusteradj && !o.clustermin ? T(N / (N-1)) : one(T)
				o.clust[c] = StrClust{T}(N, multiplier, even, order, info)
	    end

	    (o.scorebs || !o.WREnonARubin) &&
		  	(o.ClustShare = o.haswt ? @panelsum(o, o.wt, o.info⋂)/o.sumwt : T.(length.(o.info⋂)./o.Nobs)) # share of observations by group

    else  # if no clustering, cast "robust" as clustering by observation
      o.clust = [StrClust{T}(o.Nobs, o.small ? o._Nobs / (o._Nobs - one(T)) : one(T), true, Vector{Int64}(undef,0), Vector{UnitRange{Int64}}(undef,0))]
      o.NErrClustCombs = one(Int16)
      (o.scorebs || !o.WREnonARubin) &&
    		(o.ClustShare = o.haswt ? o.wt/o.sumwt : [one(T)/o.Nobs])
    end

		o.purerobust = o.robust && !o.scorebs && iszero(o.subcluster) && o.N✻==o.Nobs  # do we ever error-cluster *and* bootstrap-cluster by individual?
		o.granular   = o.WREnonARubin ? 2*o.Nobs*o.B*(2*o.N✻+1) < o.N✻*(o.N✻*o.Nobs+o.N⋂*o.B*(o.N✻+1)) :
											o.robust && !o.scorebs && (o.purerobust || (o.N⋂+o.N✻)*o.kZ*o.B + (o.N⋂-o.N✻)*o.B + o.kZ*o.B < o.N⋂*o.kZ^2 + o.Nobs*o.kZ + o.N⋂ * o.N✻ * o.kZ + o.N⋂ * o.N✻)

		if o.robust && !o.purerobust
			(o.subcluster>0 || o.granular) && !o.WREnonARubin &&
				(o.info⋂_✻⋂ = panelsetup(o.ID✻⋂, o.subcluster+1:o.NClustVar))  # info for error clusters wrt data collapsed to intersections of all bootstrapping && error clusters; used to speed crosstab UXAR wrt bootstrapping cluster && intersection of all error clusterings
			((o.scorebs && o.B>0) || (o.WREnonARubin && !o.granular && o.bootstrapt)) &&
				(o.JN⋂N✻ = zeros(T, o.N⋂, o.N✻))
		end

		if o.WREnonARubin && o.bootstrapt && !o.granular && o.NClustVar > o.NBootClustVar
			_, o._ID✻⋂ = panelsetupID(o.ID, 1:o.NClustVar)
		end
  else
	  minN = T(nrows(o.info✻))
	end

	InitFEs(o)
	if o.B>0 && o.robust && o.granular && !o.purerobust && o.bootstrapt && !o.WREnonARubin
		if o.NFE>0 && !o.FEboot
			_, o.ID✻    = panelsetupID(o.ID  , 1:o.NBootClustVar)
		else
			_, o.ID✻_✻⋂ = panelsetupID(o.ID✻⋂, 1:o.NBootClustVar)
		end
	end

	o.enumerate = o.B>0 && o.auxtwtype == :rademacher && o.N✻*log(2) < log(o.B)+1e-6  # generate full Rademacher set?
	o.enumerate && (o.maxmatsize = 0)

	o.Nw = iszero(o.maxmatsize) ? one(Int64) : ceil(Int64, (o.B+1) * Float64(max(nrows(o.ID✻), length(o.ID✻_✻⋂), o.N✻) * sizeof(T)) / o.maxmatsize / 1073741824) # 1073741824 = giga(byte)
	if isone(o.Nw)
		MakeWildWeights!(o, o.B, first=true)  # make all wild weights, once
		o.enumerate && (o.B = ncols(o.v) - 1)  # replications reduced to 2^G
		o.WeightGrp = [1:ncols(o.v)]
	else
    o.seed = rand(o.rng, UInt64)
		_B = ceil(Int64, (o.B+1) / o.Nw)
		o.Nw = ceil(Int64, (o.B+1) / _B)
		o.WeightGrp = [(i-1)*_B+1:i*_B for i ∈ 1:o.Nw]
		o.WeightGrp[end] = first(o.WeightGrp[end]):o.B+1
	end

	if o.ml
		o.dof = nrows(o.R)
	else
		if o.arubin
			o.R = hcat(zeros(o.kX₂,o.kX₁), Matrix(I(o.kX₂)))  # attack surface is all endog vars
			o.R₁ = o.kX₁>0 && nrows(o.R₁)>0 ? hcat(o.R₁[:,1:o.kX₁], zeros(nrows(o.R₁),o.kX₂)) : zeros(0, o.kX)  # and convert model constraints from referring to X₁, Y₂ to X₁, X₂
		end
		o.dof = nrows(o.R)

		if !o.WRE && iszero(o.κ)  # regular OLS
			o.DGP = StrEstimator{T}(true, o.liml, o.fuller, o.κ)
			o.Repl = StrEstimator{T}(true, false, zero(T), zero(T))  # XXX isDGP=1 for Repl? doesn't matter?
			setR!(o.DGP, o, o.null ? [o.R₁ ; o.R] : o.R₁)  # DGP constraints: model constraints + null if imposed
			setR!(o.Repl, o, o.R₁)  # model constraints only
			InitVarsOLS!(o.DGP, o, o.Repl.R₁perp)
			InitTestDenoms!(o.DGP, o)
			o.M = o.DGP  # StrEstimator object from which to get A, AR, XAR
		elseif o.arubin
			if o.willplot  # for plotting/CI purposes get original point estimate since not normally generated
				o.DGP = StrEstimator{T}(true, o.liml, o.fuller, o.κ)
				setR!(o.DGP, o, o.R₁, zeros(T,0,o.kZ))  # no-null model
				InitVarsIV!(o.DGP, o)
				EstimateIV!(o.DGP, o, o.r₁)
				o.confpeak = o.DGP.β̈  # estimated coordinate of confidence peak
			end

			o.DGP = StrEstimator{T}(true, false, zero(T), zero(T))
			setR!(o.DGP, o, o.R₁)
			InitVarsARubin!(o.DGP, o)
			InitTestDenoms!(o.DGP, o)
			o.M = o.DGP  # StrEstimator object from which to get A, AR, XAR
			o.kZ = o.kX

		elseif o.WREnonARubin
			o.Repl = StrEstimator{T}(false, o.liml, o.fuller, o.κ)
			setR!(o.Repl, o, o.R₁, o.R)
			InitVarsIV!(o.Repl, o)
			EstimateIV!(o.Repl, o, o.r₁)

			o.DGP = StrEstimator{T}(true, T(o.kX₂ ≠ o.kY₂), zero(T), one(T))
			setR!(o.DGP, o, o.null ? [o.R₁ ; o.R] : o.R₁, zeros(T,0,o.kZ))  # DGP constraints: model constraints + null if imposed
			InitVarsIV!(o.DGP, o)
			if !o.null  # if not imposing null, then DGP constraints, κ, Hessian, etc. do not vary with r and can be set now
				EstimateIV!(o.DGP, o, o.r₁)
				MakeResidualsIV!(o.DGP, o)
			end

			InitWRE!(o)

		else  # the score bootstrap for IV/GMM uses a IV/GMM DGP but then masquerades as an OLS test because most factors are fixed during the bootstrap. To conform, need DGP and Repl objects with different R, R₁, one with FWL, one not

			o.DGP = StrEstimator{T}(true, o.liml, o.fuller, o.κ)
			setR!(o.DGP, o, o.null ? [o.R₁ ; o.R] : o.R₁, zeros(T,0,o.kZ))  # DGP constraints: model constraints + null if imposed
			InitVarsIV!(o.DGP, o)
			o.Repl = StrEstimator{T}(true, o.liml, o.fuller, o.κ)
			setR!(o.Repl, o, o.R₁, I)  # process replication restraints = model constraints only
			InitVarsIV!(o.Repl, o, o.Repl.R₁perp)
			EstimateIV!(o.Repl, o, o.r₁)  # bit inefficient to estimate in both objects, but maintains the conformity
			InitTestDenoms!(o.Repl, o)
			o.M = o.Repl  # StrEstimator object from which to get A, AR, XAR; DGP follows WRE convention of using FWL, Repl follows OLS convention of not; scorebs for IV/GMM mixes the two
			if !o.null  # if not imposing null, then DGP constraints, κ, Hessian, etc. do not vary with r and can be set now
				EstimateIV!(o.DGP, o, o.r₁)
				MakeResidualsIV!(o.DGP, o)
			end
		end
  end

  if !o.WREnonARubin && o.bootstrapt
		o.denom = [Matrix{T}(undef,0,0) for _ in 1:o.dof, _ in 1:o.dof]
		if o.robust
			o.Kcd =                       Matrix{Matrix{T}}(undef, o.NErrClustCombs, o.dof)
			o.Jcd = iszero(o.B) ? o.Kcd : Matrix{Matrix{T}}(undef, o.NErrClustCombs, o.dof)  # if B = 0, Kcd will be multiplied by v, which is all 1's, and will constitute Jcd
		end

		if o.robust && o.granular<o.NErrClustCombs && o.B>0
			inds = o.subcluster>0 ?
				        [CartesianIndex(j, i) for (j,v) ∈ enumerate(o.info⋂_✻⋂) for i ∈ v] :  # crosstab ∩,* is wide
								o.NClustVar == o.NBootClustVar ?
										[CartesianIndex(i, i) for i ∈ 1:o.N✻⋂] :  # crosstab ∩,* is square
										[CartesianIndex(i, j) for (j,v) ∈ enumerate(o.clust[o.BootClust].info) for i ∈ v]  # crosstab ∩,* is tall
			o.crosstab⋂✻ind = LinearIndices(FakeArray(Tuple(max(inds...))...))[inds]
		end
	end

  o.sqrt = isone(o.dof)  # work with t/z stats instead of F/chi2?

  if o.small
		_dof = o._Nobs - o.kZ + (o.ml ? 0 : ncols(o.Repl.R₁invR₁R₁)) - o.FEdfadj
		o.dof_r = o.NClustVar>0 ? minN - one(T) : _dof  # floating-point dof_r allows for fractional fweights, FWIW...
		o.smallsample = _dof / (o._Nobs - o.robust)
	else
		o.smallsample = one(T)
		o.dof_r = zero(T)
	end
	o.clustermin && (o.smallsample *= (minN - 1) / minN)  # ivreg2-style adjustment when multiway clustering
	o.multiplier = o.small ? o.smallsample / o.dof : o.smallsample  # divide by # of constraints because F stat is so defined
  !(o.robust || o.ml) && (o.multiplier *= o._Nobs)  # will turn sum of squared errors in denom of t/z into mean
  o.sqrt && (o.multiplier = √o.multiplier)

	o.dist = fill(T(NaN), 1, o.B+1)
  (o.Nw>1 || o.WREnonARubin || (!o.null && o.dof≤2)) && (o.numer = fill(T(NaN), o.dof, o.B+1))

  if !o.WREnonARubin
		o.poles = o.anchor = zeros(T,0)
		o.interpolable = o.bootstrapt && o.null && o.Nw==1 && (iszero(o.κ) || o.arubin)
		o.interpolate_u = !(o.robust || o.ml)
		if o.interpolable
			o.∂numer∂r = Vector{Matrix{T}}(undef, o.q)
			o.interpolate_u && (o.∂u∂r = Vector{Matrix{T}}(undef, o.q))
			if o.robust
				o.∂denom∂r   = fill(Matrix{T}(undef,0,0), o.q, o.dof, o.dof)
				o.∂²denom∂r² = fill(Matrix{T}(undef,0,0), o.q, o.q, o.dof, o.dof)
				o.∂Jcd∂r     = fill(Matrix{T}(undef,0,0), o.q, o.NErrClustCombs, o.dof)
			end
		end
  end
	o.initialized = true
	nothing
end

function InitFEs(o::StrBootTest{T}) where T
	if isdefined(o, :FEID) && length(o.FEID)>0
		p = _sortperm(o.FEID)
		sortID = o.FEID[p]
		i_FE = 1; o.FEboot = o.B>0 && !o.WREnonARubin && o.NClustVar>0; j = o.Nobs; o._FEID = ones(Int64, o.Nobs)
		o.invFEwt = zeros(T, o.NFE>0 ? o.NFE : o.Nobs)
		o.FEs = Vector{StrFE{T}}(undef, o.NFE>0 ? o.NFE : o.Nobs)
		_sqrtwt = T[]
		@inbounds for i ∈ o.Nobs-1:-1:1
			if sortID[i] ≠ sortID[i+1]
				is = @view p[i+1:j]
				if o.haswt
					_sqrtwt  = @view o.sqrtwt[is]
					wt = _sqrtwt / (sumFEwt = sum(@view o.wt[is]))
				else
					sumFEwt = T(j - i)
					wt = [one(T)/sumFEwt]
				end
				o.FEs[i_FE] = StrFE{T}(is, wt, _sqrtwt)
				((o.B>0 && o.robust && o.granular < o.NErrClustVar) || (o.WREnonARubin && o.robust && o.granular && o.bootstrapt)) &&
					(o.invFEwt[i_FE] = one(T) / sumFEwt)

				j = i

				if o.FEboot  # are all of this FE's obs in same bootstrapping cluster? (But no need to check if B=0 for then CT(W.*E) in 2nd term of (62) orthogonal to v = col of 1's)
					tmp = o.ID[is, 1:o.NBootClustVar]
					o.FEboot = all(tmp .== view(tmp, 1,:)')
				end
				i_FE += 1
			end
			o._FEID[p[i]] = i_FE
		end
		is = @view p[1:j]
		if o.haswt
			_sqrtwt  = @view o.sqrtwt[is]
			wt = _sqrtwt / (sumFEwt = sum(@view o.wt[is]))
		else
			sumFEwt = T(j)
			wt = T[1/sumFEwt]
		end
		o.FEs[i_FE] = StrFE{T}(is, wt, _sqrtwt)
		o.robust && ((o.B>0 && o.granular < o.NErrClustVar) || (o.WREnonARubin && o.granular && o.bootstrapt)) &&
			(o.invFEwt[i_FE] = 1 / sumFEwt)

		if iszero(o.NFE)
			o.NFE = i_FE
			resize!(o.invFEwt, o.NFE)
			resize!(o.FEs    , o.NFE)
		end
		o.FEdfadj==-1 && (o.FEdfadj = o.NFE)

		if o.FEboot  # are all of this FE's obs in same bootstrapping cluster?
			tmp = o.ID[is, 1:o.NBootClustVar]
			o.FEboot = all(tmp .== @view tmp[1,:])
		end

		if o.robust && o.B>0 && o.bootstrapt && !o.FEboot && o.granular < o.NErrClustVar
			o.infoBootAll = panelsetup(o.ID✻⋂, 1:o.NBootClustVar)  # info for bootstrapping clusters wrt data collapsed to intersections of all bootstrapping && error clusters
		end

		o.X₁ = partialFE(o, o.X₁)  # don't overwrite caller's data
		o.X₂ = partialFE(o, o.X₂)
		o.y₁ = partialFE(o, o.y₁)
		o.Y₂ = partialFE(o, o.Y₂)
	else
		o.FEdfadj = 0
	end
	nothing
end

# draw wild weight matrix of width _B. If first=true, insert column of 1s at front. For non-Anderson-Rubin WRE, subtract 1 from all
const ϕ = (1 + √5)/2

function MakeWildWeights!(o::StrBootTest{T}, _B::Integer; first::Bool=true) where T
  if _B>0  # in scoretest or waldtest WRE, still make v a col of 1's
    o.v = o.enumerate ? o.WREnonARubin ? [zeros(o.N✻) count_binary(o.N✻, -2, 0)] :  # complete Rademacher set
								                         [ones( o.N✻) count_binary(o.N✻, -1, 1)] :
					o.auxtwtype == :normal ? randn(o.rng, T, o.N✻, _B+first) .- o.WREnonARubin :
					o.auxtwtype == :gamma  ? quantile.(Gamma{T}(4,.5),  rand(o.rng, T, o.N✻, _B+first)) .- (2 - o.WREnonARubin) :
					o.auxtwtype == :webb   ? rand(o.rng, T.([-√1.5, -1, -√.5, √.5, 1, √1.5] .- o.WREnonARubin), o.N✻, _B+first) :
				  o.auxtwtype == :mammen ? getindex.(Ref(T.([1-ϕ; ϕ] .- o.WREnonARubin)), ceil.(Int16, rand(o.rng, o.N✻, _B+first) ./ (ϕ/√5))) :
		      o.WREnonARubin         ? -2rand(o.rng, Bool, o.N✻, _B+first) : # Rademacher
			                             (o.v_sd = .5; rand(o.rng, Bool, o.N✻, _B+first) .- T(.5))
		first && !(o.enumerate && isone(o.v_sd)) && (o.v[:,1] .= o.WREnonARubin ? zero(T) : o.v_sd)  # keep original residuals in first entry to compute base model stat
  else
		o.v = Matrix{T}(undef,0,1)  # in places, ncols(v) indicates B -- 1 for classical tests
  end
	nothing
end


text:
import LinearAlgebra: UniformScaling

# ==========================
# Concrete Addition Operator
# ==========================
+ᶜ(X::LazySet, Y::LazySet) = minkowski_sum(X,Y)
+ᶜ(v::AbstractVector, X::LazySet) = translate(X,v)
+ᶜ(X::LazySet, v::AbstractVector) = translate(X,v)
⊕ᶜ = +ᶜ
+ᶜ(Z1::AbstractZonotope{N}, Z2::AbstractZonotope{N}) where {N<:Real} =
    minkowski_sum(Z1, Z2) #; remove_zero_generators=false)

# =============================
# Concrete Subtraction Operator
# =============================
-ᶜ(X::LazySet, Y::LazySet) = minkowski_difference(X,Y)
-ᶜ(v::AbstractVector, X::LazySet) = translate(reflect(X),v)
-ᶜ(X::LazySet, v::AbstractVector) = translate(X,-v)
⊖ᶜ = -ᶜ

# ================================
# Concrete Multiplication Operator
# ================================
*ᶜ(α::Real, P::LazySet) = LazySets.scale(α, P)
*ᶜ(P::LazySet, α::Real) = LazySets.scale(α, P)
*ᶜ(U::UniformScaling, P::LazySet) = linear_map(U, P)

"""
    *ᶜ(M::AbstractMatrix{N}, P::AbstractPolyhedron{N}) where {N<:Real}

Concrete linear map of a matrix `M` and polyhedral set `P`.

### Input
- `M` -- matrix
- `P` -- polyhedral set

### Output
Linear map of `P` with `M`.

### Note
If performance is of outmost importance. The computation can be sped up by using the method [`linear_map`](@ref) directly and adjusting the optinal parameter.
"""
*ᶜ(M::AbstractMatrix{N}, P::AbstractPolyhedron{N}) where {N<:Real} =
    linear_map(M, P)
*ᶜ(M::AbstractMatrix{N}, S::AbstractSingleton{N}) where {N<:Real} =
    linear_map(M, S)
*ᶜ(M::AbstractMatrix{N}, Z::AbstractZonotope{N}) where {N<:Real} =
    linear_map(M, Z)
*ᶜ(M::Array{T,2} where T, pz::PolynomialZonotope)  =
    linear_map(M, pz)

*ᶜ(α::Real, Z::Zonotope) = LazySets.scale(α, Z)
*ᶜ(Z::Zonotope, α::Real) = LazySets.scale(α, Z)
*ᶜ(α::Number, Z::PolynomialZonotope) = LazySets.scale(α, pz)
*ᶜ(Z::PolynomialZonotope, α::Number) = LazySets.scale(α, pz)

# ==============================
# Concrete intersection Operator
# ==============================
∩ᶜ(X::LazySet, Y::LazySet) = intersection(X,Y)

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
using BatchedRoutines, Test

@testset "Testing batched_scal! with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
    A = rand(elty, 100)
    B = rand(elty, 10, 10, 100)
    expect_B = copy(B)
    batched_scal!(A, B)

    for k in 1:100
        expect_B[:, :, k] .= A[k] * expect_B[:, :, k]
    end

    @test expect_B ≈ B
end

@testset "Testing batched_scal with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
    A = rand(elty, 100)
    B = rand(elty, 10, 10, 100)
    expect_B = copy(B)
    test_B = batched_scal(A, B)

    for k in 1:100
        expect_B[:, :, k] .= A[k] * expect_B[:, :, k]
    end

    @test expect_B ≈ test_B
end

@testset "Testing batched_gemm! with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]

    alpha, beta = rand(elty), rand(elty)
    A = rand(elty, 10, 10, 100)
    B = rand(elty, 10, 10, 100)
    C = rand(elty, 10, 10, 100)
    test_C = copy(C)

    batched_gemm!('N', 'N', alpha, A, B, beta, C)

    for k in 1:100
        test_C[:, :, k] = alpha * (A[:, :, k]) * B[:, :, k] + beta * test_C[:, :, k]
    end

    @test C ≈ test_C
end

@testset "Testing batched_gemm with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
    alpha = rand(elty)
    A = rand(elty, 10, 10, 100)
    B = rand(elty, 10, 10, 100)

    C = batched_gemm('N', 'N', alpha, A, B)

    test_C = zeros(elty, 10, 10, 100)
    for k in 1:100
        test_C[:, :, k] = alpha * (A[:, :, k]) * B[:, :, k]
    end

    @test test_C ≈ C
end


text:
struct LocalRuleGPUCache{R,T}
    rng::R
    σ_cpu::T
end

function RuleSamplerCache(r::LocalRule, s::MetropolisSampler, v::Union{gpuAStateBatch,gpuADoubleStateBatch},
                          net, part)
    σ_cpu = state_collect(v)
    rng = build_rng_generator_T(σ_cpu, s.seed)
    return LocalRuleGPUCache(rng, σ_cpu)
end

function init_sampler_rule_cache!(rc::LocalRuleGPUCache, s, net, σ, c)

end

# Efficient step for Local Metropolis sampling on the gpu
# avoids scalar operation on gpu by copying the state to the cpu,
# modyfing it there, then copying it back to gpu
#TODO Improve by creating a gpu kernel that does not require copies
# @timbesard suggests to pre-generate the random data with
# CuArrays.rand
function propose_step!(σp::Union{AStateBatch,ADoubleStateBatch}, s::MetropolisSampler{LocalRule},
                       net::Union{MatrixNet,KetNet}, c, rc::LocalRuleGPUCache)
    σ_cpu = rc.σ_cpu
    copy!(σ_cpu, σp)

    #switch the cpu buffer
    sites_to_flip = 1:nsites(c.hilb)
    for i=1:num_batches(σp)
        flipat = rand(rc.rng, sites_to_flip)
        flipat!(rc.rng, unsafe_get_batch(σ_cpu, i), c.hilb, flipat)
    end

    copy!(σp, σ_cpu)
end


text:
#!/usr/bin/env -S julia --compile=min -O0
using JuliaFormatter

help = """
JuliaFormatter formats Julia (.jl) programs. The formatter is width-sensitive.

Without an explicit file or path, this help message is written to stdout.
Given a file, it operates on that file; given a directory, it operates on
all .jl files in that directory, recursively.  By default, JuliaFormatter overwrites
files with the reformatted source.

Usage:

    format.jl [flags] [path ...]

Flags:

    -i, --indent
        The number of spaces used for an indentation.

    -m, --margin
        The maximum number of characters of code on a single line.  Lines over the
        limit will be wrapped if possible. There are cases where lines cannot be wrapped
        and they will still end up wider than the requested margin.

    -v, --verbose
        Print the name of the files being formatted with relevant details.

    -h, --help
        Print this message.

    -c, --check
        Check if the files are correctly formatted.  Return code 0 means nothing
        would change.  Return code 1 means some files would be reformatted.

    --always_for_in
        Always replaces `=` with `in` for `for` loops.
        Example: `for i = 1:10` will be transformed to `for i in 1:10`.

    --whitespace_typedefs
        Add whitespace in type definitions.
        Example: `Union{A <: B, C}` to `Union{A<:B,C}`.

    --whitespace_ops_in_indices`
        Add whitespace to binary ops in indices.
        Example: `arr[a + b]` to `arr[a+b]`.
        Additionally, if there's a colon `:` involved, parenthesis will be added to the LHS and RHS.
        Example: `arr[(i1 + i2):(i3 + i4)]` instead of `arr[i1+i2:i3+i4]`.
"""

function parse_opts!(args::Vector{String})
    i = 1
    opts = Dict{Symbol,Union{Int,Bool}}()
    while i ≤ length(args)
        arg = args[i]
        if arg[1] != '-'
            i += 1
            continue
        end
        if arg == "-i" || arg == "--indent"
            opt = :indent
        elseif arg == "-m" || arg == "--margin"
            opt = :margin
        elseif arg == "-v" || arg == "--verbose"
            opt = :verbose
        elseif arg == "-h" || arg == "--help"
            opt = :help
        elseif arg == "-c" || arg == "--check"
            opt = :overwrite
        elseif arg == "--always_for_in"
            opt = :always_for_in
        elseif arg == "--whitespace_typedefs"
            opt = :whitespace_typedefs
        elseif arg == "--whitespace_ops_in_indices"
            opt = :whitespace_ops_in_indices
        else
            error("invalid option $arg")
        end
        if opt in (
            :verbose,
            :help,
            :always_for_in,
            :whitespace_typedefs,
            :whitespace_ops_in_indices,
        )
            opts[opt] = true
            deleteat!(args, i)
        elseif opt == :overwrite
            opts[opt] = false
            deleteat!(args, i)
        else
            i < length(args) || error("option $arg requires and argument")
            val = tryparse(Int, args[i+1])
            val != nothing || error("invalid value for option $arg: $(args[i+1])")
            opts[opt] = val
            deleteat!(args, i:i+1)
        end
    end
    return opts
end

opts = parse_opts!(ARGS)
if isempty(ARGS) || haskey(opts, :help)
    write(stdout, help)
    exit(0)
end
exit(format(ARGS; opts...) ? 0 : 1)


text:
# To run it locally, BenchmarkCI should be added to root project
using BenchmarkCI
on_CI = haskey(ENV, "GITHUB_ACTIONS")

judgement = BenchmarkCI.judge()
if on_CI
    BenchmarkCI.postjudge()
else
    ciresult = BenchmarkCI.CIResult(; judgement=judgement)
    output_file = "benchmark_result.md"
    @info "benchmark finished" output_file
    open(output_file, "w") do io
        BenchmarkCI.printcommentmd(io, ciresult)
    end
end


text:
module EGraphs

include("../docstrings.jl")

using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

import ..Rule
using ..TermInterface

using ..Util

include("enode.jl")
include("eclass.jl")

include("egg.jl")
export find
export EClass
export ENode
export ariety
export EGraph
export AbstractAnalysis
export merge!
export addexpr!
export addanalysis!
export rebuild!

include("analysis.jl")

include("ematch.jl")
include("Schedulers/Schedulers.jl")

include("saturation_report.jl")
include("saturation.jl")
export saturate!
include("equality.jl")
export areequal
export @areequal
export @areequalg

include("extraction.jl")
export extract!
export ExtractionAnalysis
export astsize
export astsize_inv
export @extract

end


text:
abstract type lossInfo end

struct studentTeacherInfo{S, T, F<:Function} <: lossInfo
    student::S
    teacher::T
    loss_fn::F
end

get_model(l::lossInfo) = @warn "specify the function get_model for your concrete lossInfo type"
get_model(s::studentTeacherInfo) = s.student
get_teacher(s::studentTeacherInfo) = s.teacher
function get_num_weights(s::studentTeacherInfo)
    ## don't add activation layers to the number of weights
    return sum([prod(size(aa)) for aa in params(s.student) if length(size(aa)) > 1 && sum(aa) != 0])
end
# get_num_weights(s::studentTeacherInfo) = sum(prod([size(a)...]) for a in params(s.student))

function get_loss(s::studentTeacherInfo, data)
    return s.loss_fn(s.student(data), s.teacher(data))
end

"""
    get_loss_and_grad(s::studentTeacherInfo, data, params)

params are the parameters with which to get gradient with respect to
"""
function get_loss_and_grad(s::studentTeacherInfo, data, pms)
    # specific for student-teacher type losses. pms is the tracked params with which to take the gradient with respect to
    # specific for training regimes that require gradients
    loss_holder = s.loss_fn(s.student(data), s.teacher(data)) #this assigns a
    grad_holder = Tracker.gradient(() -> loss_holder, pms) # this should not use allocations
    return loss_holder, grad_holder
end

function get_loss_and_grad!(store, data; nest=false)
    # specific for student-teacher type losses. pms is the tracked params with which to take the gradient with respect to
    # specific for training regimes that require gradients
    store[:loss] = store[:lossInfo].loss_fn(store[:lossInfo].student(data), store[:lossInfo].teacher(data))
    store[:grad] = Tracker.gradient(() -> store[:loss], store[:params]; nest = nest )
    store[:loss] = Tracker.data(store[:loss])
    return store
end


text:
mutable struct Model{T, O<:MOI.AbstractOptimizer}
    params::Vector{Parameter}
    backend::ParametronMOIModel{T}
    optimizer::O
    initialized::Bool
    objective::Objective # abstract type to support different types of objectives. RTTI used in update!.
    constraints::Constraints{T}
    model_var_to_optimizer::Vector{MOI.VariableIndex}

    function Model{T}(optimizer::O) where {T, O}
        VI = MOI.VariableIndex
        params = Parameter[]
        backend = ParametronMOIModel{T}()
        initialized = false
        objective = Objective(T, @expression zero(AffineFunction{T}))
        MOI.set(backend, MOI.ObjectiveSense(), MOI.OptimizationSense(Minimize))
        MOI.set(backend, MOI.ObjectiveFunction{typeof(objective.f)}(), objective.f)
        constraints = Constraints{T}()
        model_var_to_optimizer = Vector{MOI.VariableIndex}()
        new{T, O}(params, backend, optimizer, initialized, objective, constraints, model_var_to_optimizer)
    end
end

"""
$(SIGNATURES)

Create a new `Model`, representing an optimization problem to be solved
by the optimizer `optimizer` (a `MathOptInterface.AbstractOptimizer`).
"""
Model(optimizer::MOI.AbstractOptimizer) = Model{Float64}(optimizer)

Base.show(io::IO, m::Model{T, O}) where {T, O} = print(io, "Model{$T, $O}(…)")

"""
$(SIGNATURES)

Mark all parameters associated with the model as 'dirty' (out of date),
meaning they must be updated upon their next evaluation.
"""
setdirty!(model::Model) = foreach(setdirty!, model.params)

addparameter!(model::Model, param::Parameter) = (push!(model.params, param); param)

"""
$(SIGNATURES)

Create a new decision variable (`Variable`) associated with the model.
"""
function Variable(m::Model)
    m.initialized && error("Model has already been initialized.")
    index = MOI.add_variable(m.backend)
    Variable(index)
end

"""
$(SIGNATURES)

Set the objective function and optimization sense (`Minimize` or `Maximize`).
"""
function setobjective!(m::Model{T}, sense::Sense, expr) where T
    m.initialized && error("Model was already initialized. setobjective! can only be called before initialization.")
    m.objective = Objective(T, expr)
    MOI.set(m.backend, MOI.ObjectiveSense(), MOI.OptimizationSense(sense))
    MOI.set(m.backend, MOI.ObjectiveFunction{typeof(m.objective.f)}(), m.objective.f)
    nothing
end

function add_constraint(m::Model, constraint::Constraint)
    m.initialized && error("Model was already initialized. add_constraint can only be called before initialization.")
    constraint.modelindex = MOI.add_constraint(m.backend, constraint.f, constraint.set)
    push!(m.constraints, constraint)
    nothing
end

for F in (:AffineFunction, QuadraticFunction)
    @eval begin
        _add_nonnegative_constraint!(m::Model{T}, expr, val::$F) where {T} =
            add_constraint(m, Constraint(T, expr, MOI.GreaterThan(zero(T))))
        _add_nonpositive_constraint!(m::Model{T}, expr, val::$F) where {T} =
            add_constraint(m, Constraint(T, expr, MOI.LessThan(zero(T))))
        _add_zero_constraint!(m::Model{T}, expr, val::$F) where {T} =
            add_constraint(m, Constraint(T, expr, MOI.EqualTo(zero(T))))

        _add_nonnegative_constraint!(m::Model{T}, expr, val::AbstractVector{<:$F}) where {T} =
            add_constraint(m, Constraint(T, expr, MOI.Nonnegatives(length(val))))
        _add_nonpositive_constraint!(m::Model{T}, expr, val::AbstractVector{<:$F}) where {T} =
            add_constraint(m, Constraint(T, expr, MOI.Nonpositives(length(val))))
        _add_zero_constraint!(m::Model{T}, expr, val::AbstractVector{<:$F}) where {T} =
            add_constraint(m, Constraint(T, expr, MOI.Zeros(length(val))))
    end
end

add_nonnegative_constraint!(m::Model, expr) = _add_nonnegative_constraint!(m, expr, evalarg(expr))
add_nonpositive_constraint!(m::Model, expr) = _add_nonpositive_constraint!(m, expr, evalarg(expr))
add_zero_constraint!(m::Model, expr) = _add_zero_constraint!(m, expr, evalarg(expr))

add_integer_constraint!(m::Model, x::Variable) = add_constraint(m, Constraint(MOI.SingleVariable(MOI.VariableIndex(x)), MOI.Integer()))
add_binary_constraint!(m::Model, x::Variable) = add_constraint(m, Constraint(MOI.SingleVariable(MOI.VariableIndex(x)), MOI.ZeroOne()))

function mapindices!(m::Model, idxmap)
    mapindices!(m.constraints, idxmap)
    backend_var_indices = MOI.get(m.backend, MOI.ListOfVariableIndices())
    resize!(m.model_var_to_optimizer, length(backend_var_indices))
    for index in backend_var_indices
        m.model_var_to_optimizer[index.value] = idxmap[index]
    end
end

"""
$(SIGNATURES)

Copy the problem to be solved to the optimizer.

Users should generally not need to call this function directly, as it is automatically
called the first time [`solve!`](@ref) is called on a `Model`.
"""
@noinline function initialize!(m::Model)
    indexmap = MOI.copy_to(m.optimizer, m.backend)
    mapindices!(m, indexmap)
    m.initialized = true
    nothing
end

"""
$(SIGNATURES)

Re-evaluate the expressions used to build the constraints and objective function of `Model` `m`.

Users should generally not need to call this function directly, as it is automatically
called in [`solve!`](@ref).
"""
function update!(m::Model{T}) where T
    setdirty!(m)
    if m.objective isa Objective{AffineFunction{T}, MOI.ScalarAffineFunction{T}}
        update!(m.objective::Objective{AffineFunction{T}, MOI.ScalarAffineFunction{T}}, m.optimizer, m.model_var_to_optimizer)
    elseif m.objective isa Objective{QuadraticFunction{T}, MOI.ScalarQuadraticFunction{T}}
        update!(m.objective::Objective{QuadraticFunction{T}, MOI.ScalarQuadraticFunction{T}}, m.optimizer, m.model_var_to_optimizer)
    else
        error("Unhandled objective type")
    end
    update!(m.constraints, m.optimizer, m.model_var_to_optimizer)
    nothing
end

"""
$(SIGNATURES)

Solve the model `m`. (Re-)evaluate constraint and objective expressions, update the
optimizer's internal representation of the problem, and start the optimization procedure.
"""
function solve!(m::Model)
    if !m.initialized
        initialize!(m)
        m.initialized = true
    end
    update!(m)
    MOI.optimize!(m.optimizer)
    nothing
end

"""
$(SIGNATURES)

Return the value of variable `x` as determined by the optimizer.
"""
value(m::Model, x::Variable) = MOI.get(m.optimizer, MOI.VariablePrimal(), m.model_var_to_optimizer[x.index])

"""
$(SIGNATURES)

Return the value of the objective function at the solution found by the optimizer.
"""
objectivevalue(m::Model) = MOI.get(m.optimizer, MOI.ObjectiveValue())

"""
$(SIGNATURES)

Return the termination status of the solver.
"""
terminationstatus(m::Model) = MOI.get(m.optimizer, MOI.TerminationStatus())

"""
$(SIGNATURES)

Return information regarding the primal of the problem.
"""
primalstatus(m::Model) = MOI.get(m.optimizer, MOI.PrimalStatus())

"""
$(SIGNATURES)

Return information regarding the dual of the problem.
"""
dualstatus(m::Model) = MOI.get(m.optimizer, MOI.DualStatus())

"""
$(SIGNATURES)

Add a constraint to the model using operators `==`, `<=`, `>=`, or `in`/`∈`.

`in`/`∈` may only be used for single variables with a right hand side that is
one of:

* `ℤ` or `Integers`
* `{0, 1}` or `ZeroOne`

# Examples

The constraint `x >= zeros(2)` can be added to a model as follows:

```jldoctest constraint-macro; setup = :(using Parametron; model = Parametron.mock_model())
julia> x = [Variable(model) for i = 1 : 2];

julia> @constraint(model, x >= zeros(2))
```

The constraint that variable `x[1]` should be an integer can be
expressed using:

```jldoctest constraint-macro
julia> @constraint(model, x[1] ∈ ℤ)
```
"""
macro constraint(model, expr)
    if !(expr isa Expr) || expr.head != :call || length(expr.args) != 3
        return :(throw(ArgumentError("Expected expression of the form `a relation b`")))
    end
    relation = expr.args[1]
    lhs = expr.args[2]
    rhs = expr.args[3]
    if relation == :(>=)
        ret = :(Parametron.add_nonnegative_constraint!($model, Parametron.@expression $lhs - $rhs))
    elseif relation == :(<=)
        ret = :(Parametron.add_nonpositive_constraint!($model, Parametron.@expression $lhs - $rhs))
    elseif relation == :(==)
        ret = :(Parametron.add_zero_constraint!($model, Parametron.@expression $lhs - $rhs))
    elseif relation ∈ [:∈, :in]
        if rhs ∈ [:ℤ, :Integers]
            ret = :(Parametron.add_integer_constraint!($model, $lhs))
        elseif rhs ∈ [:({0, 1}), :ZeroOne]
            ret = :(Parametron.add_binary_constraint!($model, $lhs))
        else
            return :(throw(ArgumentError("'in' only supports ℤ/Integers and {0, 1}/ZeroOne")))
        end
    else
        return :(throw(ArgumentError("Relation not recognized")))
    end
    esc(ret)
end

"""
$(SIGNATURES)

Set the objective function of the model.

# Examples

Let `model` be a `Model` instance. The objective 'minimize x ⋅ x' can
be added as follows:

```jldoctest; setup = :(using Parametron; using LinearAlgebra; model = Parametron.mock_model())
julia> x = [Variable(model) for i = 1 : 2];

julia> @objective model Minimize x ⋅ x;
```
"""
macro objective(model, sense, expr)
    quote
        Parametron.setobjective!($model, $sense, Parametron.@expression $expr)
    end |> esc
end

# Support constructs like value.(model, x)
Broadcast.broadcastable(m::Model) = Ref(m)


text:
"""
``
eqcond(m::OneAssetHANK)
```

Expresses the equilibrium conditions in canonical form using Γ0, Γ1, C, Ψ, and Π matrices.

### Outputs

* `Γ0` (`n_states` x `n_states`) holds coefficients of current time states.
* `Γ1` (`n_states` x `n_states`) holds coefficients of lagged states.
* `C`  (`n_states` x `1`) is a vector of constants
* `Ψ`  (`n_states` x `n_shocks_exogenous`) holds coefficients of iid shocks.
* `Π`  (`n_states` x `n_states_expectational`) holds coefficients of expectational states.
"""
function eqcond(m::OneAssetHANK)
    nstates = n_states(m)
    n_s_exp = n_shocks_expectational(m)
    n_s_exo = n_shocks_exogenous(m)

    x       = zeros(Float64, 2 * nstates + n_s_exp + n_s_exo)

    # Read in parameters, settings
    niter_hours       = get_setting(m, :niter_hours)
    n_v               = get_setting(m, :n_jump_vars)
    n_g               = get_setting(m, :n_state_vars)
    ymarkov_combined  = get_setting(m, :ymarkov_combined)
    zz                = get_setting(m, :zz)
    a                 = get_setting(m, :a)
    I, J              = size(zz)
    amax              = get_setting(m, :amax)
    amin              = get_setting(m, :amin)
    aa                = repeat(a, 1, J)
    A_switch          = kron(ymarkov_combined, speye(Float64, I))
    daf, dab, azdelta = initialize_diff_grids(a, I, J)

    # Set up steady state parameters
    V_ss = m[:V_ss].value
    g_ss = m[:g_ss].value
    G_ss = m[:G_ss].value
    w_ss = m[:w_ss].value
    N_ss = m[:N_ss].value
    C_ss = m[:C_ss].value
    Y_ss = m[:Y_ss].value
    B_ss = m[:B_ss].value
    r_ss = m[:r_ss].value
    ρ_ss = m[:ρ_ss].value
    θ_MP = m[:θ_MP].value
    σ_MP = m[:σ_MP].value

    h_ss              = reshape(m[:h_ss].value, I, J)
    ceselast          = m[:ceselast].value
    inflation_ss      = m[:inflation_ss].value
    maxhours          = m[:maxhours].value
    govbcrule_fixnomB = m[:govbcrule_fixnomB].value
    priceadjust       = m[:priceadjust].value
    taylor_inflation  = m[:taylor_inflation].value
    taylor_outputgap  = m[:taylor_outputgap].value
    meanlabeff        = m[:meanlabeff].value

    # Necessary for construction of household problem functions
    labtax            = m[:labtax].value
    coefrra           = m[:coefrra].value
    frisch            = m[:frisch].value
    labdisutil        = m[:labdisutil].value

    # Note: make TFP a setting?
    TFP               = 1.0

    @inline function get_residuals(x::Vector{T}) where {T<:Real}
        # Prepare steady state deviations
        V           = reshape(x[1:n_v - 1]       + V_ss, I, J)
        inflation   = x[n_v]                     + inflation_ss
        gg          = x[n_v + 1 : n_v + n_g - 1] + g_ss[1:end-1]
        MP          = x[n_v + n_g]
        w           = x[n_v + n_g + 1]           + w_ss
        hours       = x[n_v + n_g + 2]           + N_ss
        consumption = x[n_v + n_g + 3]           + C_ss
        output      = x[n_v + n_g + 4]           + Y_ss
        assets      = x[n_v + n_g + 5]           + B_ss

        # Set up system of time derivatives and expectational errors, etc.
        V_dot           = x[nstates + 1 : nstates + n_v - 1]
        inflation_dot   = x[nstates + n_v]
        g_dot           = x[nstates + n_v + 1 : nstates + n_v + n_g - 1]
        mp_dot          = x[nstates + n_g + n_v]
        VEErrors        = x[2*nstates + 1 : 2*nstates + n_v - 1]
        inflation_error = x[2*nstates + n_v]
        mp_shock        = x[2*nstates + n_v + 1]

        g_end = (1 - dot(gg, azdelta[1:end-1])) ./ azdelta[end]
        g     = vcat(gg, g_end)

        #-----------------------------------------------------------------
        # Get equilibrium values, given steady state values
        #-----------------------------------------------------------------
        normalized_wage = w / TFP
        profshare       = (zz ./ meanlabeff) .* ((1. - normalized_wage) .* output)
        r_nominal       = r_ss + taylor_inflation * inflation + taylor_outputgap * (log(output)-log(Y_ss)) + MP
        r               = r_nominal - inflation
        lumptransfer    = labtax * w * hours - G_ss -
            (r_nominal - (1-govbcrule_fixnomB) * inflation) * assets

        #-----------------------------------------------------------------
        # Compute one iteration of the HJB
        #-----------------------------------------------------------------

        # Get flow utility, income, and labor hour functions
        util, income, labor = construct_household_problem_functions(V, w, coefrra, frisch, labtax, labdisutil)

        # Initialize other variables, using V to ensure everything is a dual number
        Vaf = copy(V)
        Vab = copy(V)
        cf  = similar(V)
        hf  = similar(V)
        sf  = similar(V)
        cb  = similar(V)
        hb  = similar(V)
        sb  = similar(V)
        c0  = similar(V)
        h0  = similar(V)
        for j=1:J, i=1:I
            h0[i,j] = h_ss[i,j]
        end

        #-----------------------------------------------------------------
        # Construct Initial Difference Matrices
        #-----------------------------------------------------------------
        for j=1:J, i=1:I
            if i==I
                Vaf[end,j] = income(h_ss[end,j], zz[end,j], profshare[end,j], lumptransfer, r, amax)^(-coefrra)
                Vab[i, j]  = (V[i,j] - V[i-1,j]) / dab[i]
            elseif i==1
                Vaf[i, j] = (V[i+1,j] - V[i,j]) / daf[i]
                Vab[1, j] = income(h_ss[1, j], zz[1, j], profshare[1, j], lumptransfer, r, amin)^(-coefrra)
            else
                Vaf[i, j] = (V[i+1,j] - V[i,j]) / daf[i]
                Vab[i, j] = (V[i,j] - V[i-1,j]) / dab[i]
            end

            cf[i,j] = Vaf[i,j] ^ (-1 / coefrra)
            cb[i,j] = Vab[i,j] ^ (-1 / coefrra)

            hf[i,j] = min(labor(zz[i,j], Vaf[i,j]), maxhours)
            hb[i,j] = min(labor(zz[i,j], Vab[i,j]), maxhours)
        end

        #-----------------------------------------------------------------
        # Hours Iteration
        #-----------------------------------------------------------------
        for ih=1:niter_hours
            for j=1:J, i=1:I
                if i == I
                    cf[end, j]  = income(hf[end, j], zz[end, j], profshare[end, j], lumptransfer, r, aa[end, j])
                    hf[end, j]  = min(labor(zz[end, j], cf[end, j]^(-coefrra)), maxhours)

                    if (ih == niter_hours) Vaf[end, j] = cf[I, j] ^ (-coefrra) end

                elseif i == 1
                    cb[1, j]  = income(hb[1, j], zz[1, j], profshare[1, j], lumptransfer, r, aa[1, j])
                    hb[1, j]  = min(labor(zz[1, j], cb[1, j]^(-coefrra)), maxhours)

                    if (ih == niter_hours) Vab[1, j] = cb[1, j] ^ (-coefrra) end
                end

                c0[i, j] = income(h0[i, j], zz[i, j], profshare[i, j], lumptransfer, r, aa[i, j])
                h0[i, j] = min(labor(zz[i, j], c0[i, j]^(-coefrra)), maxhours)
            end
        end

        for i in eachindex(h0)
            c0[i] = income(h0[i], zz[i], profshare[i], lumptransfer, r, aa[i])
            sf[i] = income(hf[i], zz[i], profshare[i], lumptransfer, r, aa[i]) - cf[i]
            sb[i] = income(hb[i], zz[i], profshare[i], lumptransfer, r, aa[i]) - cb[i]
        end

        A, u, h, c, s  = upwind(util, A_switch, cf, cb, c0, hf, hb, h0, sf, sb, Vaf, Vab, daf, dab)

        #-----------------------------------------------------------------
        # Collect/calculate Residuals
        #-----------------------------------------------------------------
        hjb_residual = vec(u) + A * vec(V) + V_dot + VEErrors - ρ_ss * vec(V)
        pc_residual  = -((r - 0) * inflation - (ceselast / priceadjust *
                                                (w / TFP - (ceselast-1) / ceselast) +
                                                inflation_dot - inflation_error))

        g_azdelta            = vec(g) .* vec(azdelta)
        g_intermediate       = spdiagm(0 => 1 ./ azdelta) * A' * g_azdelta
        g_residual           = g_dot - g_intermediate[1:end-1]

        mp_residual          = mp_dot - (-θ_MP * MP + σ_MP * mp_shock)

        realsav              = sum(vec(aa) .* vec(g) .* vec(azdelta))
        realsav_dot          = sum(vec(s)  .* vec(g) .* vec(azdelta))
        bondmarket_residual  = realsav_dot/realsav + govbcrule_fixnomB * inflation

        labmarket_residual   = sum(vec(zz) .* vec(h) .* vec(g) .* vec(azdelta)) - hours
        consumption_residual = sum(vec(c) .* vec(g) .* vec(azdelta)) - consumption
        output_residual      = TFP * hours - output
        assets_residual      = assets - realsav

        # Return equilibrium conditions
        return [hjb_residual; pc_residual; g_residual; mp_residual; bondmarket_residual; labmarket_residual;
                consumption_residual; output_residual; assets_residual]
    end

    derivs = ForwardDiff.jacobian(get_residuals, x)

    Γ1 = -derivs[:, 1:nstates]
    Γ0 =  derivs[:,   nstates +           1:2 * nstates]
    Π  = -derivs[:, 2*nstates +           1:2 * nstates + n_s_exp]
    Ψ  = -derivs[:, 2*nstates + n_s_exp + 1:2 * nstates + n_s_exp + n_s_exo]
    C  = zeros(nstates)

    if typeof(Ψ) == Vector{Float64}
        Ψ = reshape(Ψ, length(Ψ), n_s_exo)
    end

    return Γ0, Γ1, Ψ, Π, C
end


text:
using LibFoo
using Test

@test 1 == 1

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
using Test, Random, Statistics
Random.seed!(0)
using Hyperopt, Plots
using Optim

@testset "BOHB" begin
    f(a;c=10) = sum(@. 100 + (a-3)^2 + (c-100)^2)
    @test_nowarn Hyperband(50)
    let bohb = @hyperopt for i=18, sampler=Hyperband(R=50, η=3, inner=BOHB(dims=[Continuous(), Continuous()])), a = LinRange(1,5,800), c = exp10.(LinRange(-1,3,1800))
            # println(i, "\t", a, "\t", b, "\t", c)
            if !(state === nothing)
                a,c = state
            end
            res = Optim.optimize(x->f(x[1],c=x[2]), [a,c], NelderMead(), Optim.Options(f_calls_limit=floor(Int, i)))
            Optim.minimum(res), Optim.minimizer(res)
        end
        @test length(bohb.history) == 69
        @test length(bohb.results) == 69
        @test minimum(bohb) < 300
        @hyperopt for i=1, ho=bohb, sampler=Hyperband(R=50, η=3, inner=BOHB(dims=[Continuous(), Continuous()])), a = LinRange(1,5,800), c = exp10.(LinRange(-1,3,1800))
            # println(i, "\t", a, "\t", b, "\t", c)
            if !(state === nothing)
                a,c = state
            end
            res = Optim.optimize(x->f(x[1],c=x[2]), [a,c], NelderMead(), Optim.Options(f_calls_limit=floor(Int, i)))
            Optim.minimum(res), Optim.minimizer(res)
        end
        @test length(bohb.history) == 138
        @test length(bohb.results) == 138
    end

    # extra robust option

    ho = @hyperopt for i=18, sampler=Hyperband(R=50, η=3, inner=BOHB(dims=[UnorderedCategorical(5), Continuous(), Continuous()])),
        algorithm = [SimulatedAnnealing(), ParticleSwarm(), NelderMead(), BFGS(), NewtonTrustRegion()],
        a = LinRange(1,5,1800),
        c = exp10.(LinRange(-1,3,1800))
        if !(state === nothing)
            algorithm, a, c = state # should come in same order as they appear in the argument list above.
        end
        println(i, " algorithm: ", typeof(algorithm))
        res = Optim.optimize(x->f(x[1],c=x[2]), [a,c], algorithm, Optim.Options(time_limit=2i+2, show_trace=false, show_every=5))
        Optim.minimum(res), (algorithm, Optim.minimizer(res)...)
    end

end


text:
function get_current_state(obj::OrekitStepInterpolator)
    return jcall(obj, "getCurrentState", SpacecraftState, ())
end

function get_interpolated_state(obj::OrekitStepInterpolator, arg0::AbsoluteDate)
    return jcall(obj, "getInterpolatedState", SpacecraftState, (AbsoluteDate,), arg0)
end

function get_previous_state(obj::OrekitStepInterpolator)
    return jcall(obj, "getPreviousState", SpacecraftState, ())
end

function is_current_state_interpolated(obj::OrekitStepInterpolator)
    return jcall(obj, "isCurrentStateInterpolated", jboolean, ())
end

function is_forward(obj::OrekitStepInterpolator)
    return jcall(obj, "isForward", jboolean, ())
end

function is_previous_state_interpolated(obj::OrekitStepInterpolator)
    return jcall(obj, "isPreviousStateInterpolated", jboolean, ())
end

function restrict_step(obj::OrekitStepInterpolator, arg0::SpacecraftState, arg1::SpacecraftState)
    return jcall(obj, "restrictStep", OrekitStepInterpolator, (SpacecraftState, SpacecraftState), arg0, arg1)
end


text:
## Model to evaluate QMDP Policy

mutable struct QMDPEval <: Policy
    env::CrosswalkEnv
    pomdp::SingleOCPOMDP
    policy::Policy
    a::Float64
    tick::Int64
    Δt::Float64
end

function QMDPEval(env::CrosswalkEnv, pomdp::SingleOCPOMDP, policy::Policy)
    return QMDPEval(env, pomdp, policy, 0., 0, pomdp.ΔT)
end

function reset_policy!(policy::QMDPEval)
    policy.tick = 0
    policy.a = 0.
end

function action(policy::QMDPEval, b::Dict{Int64, SingleOCDistribution},  verbose::Bool = false)
    policy.a = POMDPs.action(policy.policy, b).acc
    policy.tick += 1
    return policy.a
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using Sockets

export AbstractConnection

abstract type AbstractConnection end

"""
    ConnectionPool(outbox::Channel, connections::Set{AbstractConnection}=Set())

Manages the distribution of messages from the `outbox` channel to a set of
connections. The ConnectionPool asynchronously takes messages from its outbox
and sends each message to every connection in the pool. Any closed connections
are automatically removed from the pool.
"""
struct ConnectionPool
    outbox::Channel
    connections::Set{AbstractConnection}
    new_connections::Channel{AbstractConnection}

    function ConnectionPool(outbox::Channel, connections=Set{AbstractConnection}())
        new_connections = Channel{AbstractConnection}(32)
        pool = new(outbox, connections, new_connections)
        @async process_messages(pool)
        pool
    end
end

addconnection!(pool::ConnectionPool, c::AbstractConnection) = put!(pool.new_connections, c)
Sockets.send(pool::ConnectionPool, msg) = put!(pool.outbox, msg)

"""
    ensure_connection(pool::ConnectionPool)

Ensure that the pool has at least one connection, potentially blocking
the current task until that is the case.
"""
function ensure_connection(pool::ConnectionPool)
    if isempty(pool.connections)
        wait(pool.new_connections)
    end
    while isready(pool.new_connections)
        push!(pool.connections, take!(pool.new_connections))
    end
end

Base.wait(pool::ConnectionPool) = ensure_connection(pool)

function process_messages(pool::ConnectionPool)
    while true
        msg = fetch(pool.outbox)  # don't take! yet, since we're not sure msg can be sent
        ensure_connection(pool)
        msg_sent = false
        @sync begin
            for connection in pool.connections
                @async begin
                    if isopen(connection)
                        send(connection, msg)
                        msg_sent = true
                    else
                        delete!(pool.connections, connection)
                    end
                end
            end
        end
        if msg_sent
            # Message was sent, so remove it from the channel
            take!(pool.outbox)
        end
    end
end