text:
"""
    FirstPassageAbstract

Types ihneriting from this struct define the type of information stored about
first-passage times
"""
abstract type FirstPassageAbstract end

"""
    NoFirstPassageTimes <: FirstPassageAbstract

Compile-time indicator that the observation stores no first-passage time data
"""
struct NoFirstPassageTimes <: FirstPassageAbstract end

"""
    FirstPassageTimeInfo{C,L,U,A,R} <: FirstPassageAbstract

Compile-time indicator for the first-passage time information conveyed by the
data-point. `C` lists the affected coordinates (of the observations, not the
process), `L` indicates the level, `U` are indicators for whether the respective
first-passage times are up-crossings and `A` is and indicator for whether
an additional reset needs to be reached before first-passage time (note that
reset time is defined as a first-passage time to a corresponding reset level
stored in `R` that happens anytime before the actual first-passage time,
the reset-time happens in the direction opposite to the direction of
first-passage time crossing; additionally, the coordinate can reach the
first-passage time level prior to being reset).

    FirstPassageTimeInfo(
        coords,
        levels,
        upcrossings,
        additional_reset_required,
        reset_levels=tuple()
    )

Base constructor.
"""
struct FirstPassageTimeInfo{C,L,U,A,R} <: FirstPassageAbstract

    FirstPassageTimeInfo{C,L,U,A,R}() where {C,L,U,A,R} = new{C,L,U,A,R}()

    function FirstPassageTimeInfo(
            coords,
            levels,
            upcrossings,
            additional_reset_required,
            reset_levels=tuple()
        )
        N = length(coords)
        @assert length(levels) == length(upcrossings) == N
        @assert length(additional_reset_required) == N
        @assert all( map(c->(typeof(c)<:Integer), coords) )
        @assert all( map(l->(typeof(l)<:AbstractFloat), levels) )
        @assert all( map(u->(typeof(u)==Bool), upcrossings) )
        @assert all( map(a->(typeof(a)==Bool), additional_reset_required) )
        C, L, U = Tuple(coords), Tuple(levels), Tuple(upcrossings)
        A = Tuple(additional_reset_required)

        M = length(reset_levels)
        any(A) && @assert (M == N || M == sum(A))
        R_vec = Union{Float64,Nothing}[]
        if M == N
            # just to keep with the convention
            for i in 1:N
                A[i] && push!(R_vec, reset_levels[i])
                !A[i] && push!(R_vec, nothing)
            end
        else
            counter = 1
            for i in 1:N
                A[i] && (push!(R_vec, reset_levels[counter]); counter += 1)
                !A[i] && push!(R_vec, nothing) # this is never used
            end
        end
        R = Tuple(R_vec)

        new{C,L,U,A,R}()
    end
end

strip_fpt_info(::Type{FirstPassageTimeInfo{C,L,U,A,R}}) where {C,L,U,A,R} = (C,L,U,A,R)

cpad(s, n) = cpad(string(s), n)

function cpad(s::String, n::Integer)
    rpad_n = div(max(0, n-length(s)), 2)+length(s)
    lpad(reverse(lpad(reverse(s), rpad_n)),n)
end

Base.summary(io::IO, fpt::Type{<:FirstPassageTimeInfo}) = _summary_fpt(io, fpt)
Base.summary(fpt::Type{<:FirstPassageTimeInfo}) = summary(Base.stdout, fpt)

function _summary_fpt(io::IO, ::Type{NoFirstPassageTimes}; prepend="")
    println(io, prepend, "No first passage times recorded.")
end

function _summary_fpt(io::IO, fpt; prepend="")
    (
        coords,
        levels,
        upcrossings,
        additional_reset_required,
        reset_levels
    ) = strip_fpt_info(fpt)
    pad_len = 14
    println(io, prepend, "First passage times of the observation `v`: ")
    println(io, prepend, repeat("-", pad_len*5+6))
    println(
        io,
        prepend,
        "|", cpad("coordinate", pad_len), "|", cpad("level", pad_len),
        "|", cpad("up-crossing", pad_len), "|", cpad("extra reset", pad_len),
        "|", cpad("reset level", pad_len), "|"
    )
    println(io, prepend, repeat("-", pad_len*5+6))
    for (coord, lev, upcross, addres, reslvl) in zip(strip_fpt_info(fpt)...)
        println(
            io,
            prepend,
            "|", cpad(coord, pad_len), "|", cpad(lev, pad_len),
            "|", cpad((upcross ? "up-crossing" : "down-crossing"), pad_len),
            "|", cpad((addres ? "✔" : "✘"), pad_len),
            "|", cpad((addres ? reslvl : "✘"), pad_len), "|"
        )
    end
end


text:
##############################
#   BasicSociety Type   #
##############################

"""
    type BasicSociety{N, K, T}

Type representing a generic Bayesian Agent Society

### Fields
* links  (Matrix{Bool}): N×N-matrix that describes the interactions of the agents of the society
* agents (Vector{MoralVector}): N-vector of the MoralVector{K, T}s of the BasicSociety
* ρ (Float64): Cognitive Style of the agents
* ε (Float64): Distrust of the agents
"""
struct BasicSociety{N, K, T <: Real} <: Society{N, K, T}
    links::Matrix{Bool}
    agents::Vector{MoralVector{K, T}}
    ρ::Vector{Float64}
    ε::Matrix{Float64}
end

"""
    BasicSociety()

Construct a random BasicSociety with default size NSOC and parameters (ρDEF, εDEF).
The agents have the default number of components (KMORAL)
"""
BasicSociety() =
    BasicSociety{NSOC, KMORAL, Float64}(Matrix{Bool}(1 - eye(NSOC)), MoralVector{KMORAL, Float64}[MoralVector() for i in 1:NSOC], ρDEF*ones(NSOC), εDEF*(ones(NSOC, NSOC) - eye(NSOC, NSOC)))

"""
    BasicSociety(n::Int)

Construct a random BasicSociety with size n and parameters (ρDEF, εDEF).
The agents have the default number of components (KMORAL)
"""
BasicSociety(n::Int) =
    BasicSociety{n, KMORAL, Float64}(Matrix{Bool}(1 - eye(n)), MoralVector{KMORAL, Float64}[MoralVector() for i in 1:n], ρDEF*ones(n), εDEF*(ones(n, n) - eye(n, n)))

"""
    BasicSociety(n::Int, k::Int)

Construct a random BasicSociety with size n and parameters (ρDEF, εDEF).
The agents have k number of components
"""
BasicSociety(n::Int, k::Int) =
    BasicSociety{n, k, Float64}(Matrix{Bool}(1 - eye(n)), MoralVector{k, Float64}[MoralVector(k) for i in 1:n], ρDEF*ones(n), εDEF*(ones(n, n) - eye(n, n)))

"""
    BasicSociety(n::Int, ρ::Vector{Float64}, ε::Matrix{Float64})

Construct a random BasicSociety with size n and parameters (ρ, ε).
The agents have the default number of components (KMORAL)
"""
BasicSociety(n::Int, ρ::Vector{Float64}, ε::Matrix{Float64}) =
    BasicSociety{n, KMORAL, Float64}(Matrix{Bool}(1 - eye(n)), MoralVector{KMORAL, Float64}[MoralVector() for i in 1:n], ρ, ε)

"""
    BasicSociety(agentArray::Vector{MoralVector{K, T}}, ρ::Vector{Float64}, ε::Matrix{Float64})

Construct a BasicSociety with agents agentArray and parameters (ρ, ε).
"""
BasicSociety{K, T}(agentArray::Vector{MoralVector{K, T}}, ρ::Vector{Float64}, ε::Matrix{Float64}) =
    BasicSociety{length(agentArray), K, T}(Matrix{Bool}(1 - eye(length(agentArray))), deepcopy(agentArray), ρ, ε)

######################################
#   BasicSociety Redefinitions  #
######################################

function show(io::IO, soc::BasicSociety)
    N, K = size(soc)
    println(io, N, "-sized BasicSociety on a ", K, "-dimensional Moral space")
end

######################################
#    BasicSociety Definitions   #
######################################

"""
    epssoc(soc::BasicSociety, i::Int, j::Int)

Computes `ε` agent `i` in a Society has about agent `j`
"""
epssoc(soc::BasicSociety, i::Int, j::Int) = soc.ε[i, j]

"""
    epssoc(soc::BasicSociety)

Computes the matrix of distrust `ε` in a Society
"""
epssoc{N, K, T}(soc::BasicSociety{N, K, T}) = [epssoc(soc, i, j) for i in 1:N, j in 1:N]

"Computes `ρ` of agent `i` in a Society"
rhosoc(soc::BasicSociety, i::Int) = soc.ρ[i]

"Computes `γ` of agent `i` in a Society"
gammasoc(soc::BasicSociety, i::Int) = gammasoc(rhosoc(soc, i))

"Computes `γ` of agent `i` in a Society given MoralVector `x`"
gammasoc{N, K, T}(soc::BasicSociety{N, K, T}, i::Int, x::MoralVector{K, T}) = sqrt(x[:]' * soc.C[i] * x[:])

"Computes `ρ` of agent `i` in a Society given MoralVector `x`"
rhosoc{N, K, T}(soc::BasicSociety{N, K, T}, i::Int, x::MoralVector{K, T}) = rhosoc(gammasoc(soc, i, x))

"""
    cogcost{K,T}(soc::BasicSociety{N, K, T}, i::Int, j::Int, x::MoralVector{K,T})

Cognitive cost MoralVector `soc[i]` suffers when learning MoralVector `soc[j]`'s opinion about MoralVector `x`
"""
function cogcost{N, K,T}(soc::BasicSociety{N, K, T}, i::Int, j::Int, x::MoralVector{K,T})
    γi  = gammasoc(soc, i)
    εij = epssoc(soc, i, j)
    agi = soc[i]
    agj = soc[j]

    return - log( εij + (1 - 2εij) * phi(sign(agj ⋅ x) * (agi ⋅ x) / γi) )
end

"""
    cogcost{K,T}(soc::BasicSociety{N, K, T}, i::Int, j::Int, x::MoralVector{K,T})

Cognitive cost MoralVector `agi` suffers when learning MoralVector `soc[i]`'s opinion about MoralVector `x`
"""
function cogcost{N, K,T}(agi::MoralVector{K,T}, soc::BasicSociety{N, K, T}, i::Int, j::Int, x::MoralVector{K,T})
    γi  = gammasoc(soc, i)
    εij = epssoc(soc, i, j)
    agj = soc[j]

    return - log( εij + (1 - 2εij) * phi(sign(agj ⋅ x) * (agi ⋅ x) / γi) )
end

ext: jl


text:
import Pkg

Pkg.instantiate()
srcpath = joinpath(@__DIR__(), "Courses")
buildpath = joinpath(@__DIR__(), "Notebooks")

if isempty(ARGS)
    courses = filter(x->isdir(joinpath(srcpath), x) && x !== "Template", readdir(srcpath))
else
    courses = ARGS
end

rm(buildpath, recursive=true, force=true)
mkdir(buildpath)

notebooks_failed = String[]

Base.julia_cmd()
for dir in courses
    @info "building $dir..."
    src_course = joinpath(srcpath, dir)
    build_course = mkdir(joinpath(buildpath, dir))
    cp(src_course, build_course; force=true)

    Pkg.activate(build_course)
    Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaAcademy/JuliaAcademyData.jl.git"))
    Pkg.resolve()
    test = false
    Pkg.instantiate()
    for file in readdir(build_course)
        if endswith(file, ".jl")
            course = joinpath(build_course, file)
            academy_environment = @__DIR__
            # Create a stacked environment with the JuliaAcademy environment in
            # it so we can load Literate for each notebook
            script = """
            pushfirst!(LOAD_PATH, $(repr(academy_environment)));
            import Literate;
            Literate.notebook($(repr(course)), $(repr(build_course)); credit=false, execute=false)
            """
            try
                run(`$(Base.julia_cmd()) -e $script`)
            catch
                err_src_course = relpath(joinpath(src_course, file), @__DIR__)
                @error "Script failed to build: $err_src_course)"
                push!(notebooks_failed, err_src_course)
            end
            rm(course)
        end
    end
end

if !isempty(notebooks_failed)
    error("The following notebooks failed: \n",
          join(notebooks_failed, '\n'))
end

max_stars_repo_name: pnavaro/JuliaAcademyMaterials

max_issues_repo_name: pnavaro/JuliaAcademyMaterials

max_forks_repo_name: pnavaro/JuliaAcademyMaterials


text:
using YAML
using ACE
using ACE: PIBasis, PIBasisFcn, PIPotential
using ACE.OrthPolys: TransformedPolys
using ACE: rand_radial, cutoff, numz
using JuLIP: energy, bulk, i2z, z2i, chemical_symbol

function _basis_groups(inner, coeffs)
    ## grouping the basis functions
    NL = []
    M = []
    C = []
    for b in keys(inner.b2iAA)
       if coeffs[ inner.b2iAA[b] ] != 0
          push!(NL, ( [b1.n for b1 in b.oneps], [b1.l for b1 in b.oneps] ))
          push!(M, [b1.m for b1 in b.oneps])
          push!(C, coeffs[ inner.b2iAA[b] ])
       end
    end
    ords = length.(M)
    perm = sortperm(ords)
    NL = NL[perm]
    M = M[perm]
    C = C[perm]
    @assert issorted(length.(M))
    bgrps = []
    alldone = fill(false, length(NL))
    for i = 1:length(NL)
       if alldone[i]; continue; end
       nl = NL[i]
       Inl = findall(NL .== Ref(nl))
       alldone[Inl] .= true
       Mnl = M[Inl]
       Cnl = C[Inl]
       pnl = sortperm(Mnl)
       Mnl = Mnl[pnl]
       Cnl = Cnl[pnl]
       push!(bgrps, Dict("n" => nl[1], "l" => nl[2],
                         "M" => Mnl, "C" => Cnl))
    end
    return bgrps
end

function _write_group(g; ndensitymax = 1)
    ## write functions to function dict in YACE format
    order = length(g["l"])
    mu0 = 0
    #mu = order
    num_ms = length(g["M"])
    ns = g["n"]
    ls = g["l"]
    c_tildes = (g["C"] ./ (4*π)^(order/2))
    return "{mu0: $(mu0), rank: $(order), ndensitymax: $(ndensitymax), num_ms_combs: $(num_ms), mus!!, ns: $(ns), ls: $(ls), ms_combs!!!, c_tildes: $(c_tildes))}"
end

function export_ACE(fname, IP)
    #decomposing into V1, V2, V3 (One body, two body and ACE bases)
    V1 = IP.components[1]
    V2 = IP.components[2]
    V3 = IP.components[3]

    #grabbing the species and making a dict "0" => "Al", "1" => "Ti"
    #needs to be consistent everywhere!
    species = string.(collect(keys(IP.components[1].E0)))
    nspecies = length(species)

    #0 index instead of 1 index like Julia
    species_dict = Dict(zip(collect(0:length(species)-1), species))

    #creating "data" dict where we'll store everything
    data = Dict()

    # not sure about both keys below... ask Yury?
    # both are not in the previous verison of PACE
    data["deltaSplineBins"] = "0.001"
    data["embeddings"] = "none?" ##

    # grabbing the elements key and E0 key from the onebody (V1)
    elements, E0 = export_ACE_V1(V1, species_dict)
    data["elements"] = elements
    data["E0"] = E0

    #storing all info with PACE
    bonds = export_ACE_V2(V2, species_dict)
    data["bonds"] = bonds

    functions = export_ACE_V3(V3, species_dict)
    data["functions"] = functions
    ## 1, 2 indices for multispecies support

    YAML.write_file(fname, data)
end

function YACE_format(l)
    return "[" * join(l, ",") * "]"
end

function export_ACE_V1(V1, species_dict)
    E0 = []
    elements = []
    for species_ind in keys(species_dict)
        push!(E0, V1(Symbol(species_dict[species_ind])))
        push!(elements, species_dict[species_ind])
    end
    return YACE_format(elements), YACE_format(E0)
end

function export_ACE_V2(V2, species_dict)
    #these are ALL coeffs, but these are species dependent [0,0], [1,1]
    #we need to decompose coeffs and assign them in the right bond key
    #I'm not sure how..
    coeffs = V2.coeffs

    #grabbing the transform and basis
    transbasis = V2.basis.J
    Pr = V2.basis

    #grabbing all the required params
    p = transbasis.trans.p
    r0 = transbasis.trans.r0
    xr = transbasis.J.tr
    xl = transbasis.J.tl
    pr = transbasis.J.pr
    pl = transbasis.J.pl
    rcut = cutoff(Pr)
    maxn = length(Pr)

    #guessing "radbasname" is that just "polypairpots"
    radbasename = "polypairpots"

    bonds = Dict()
    #this does not respect the coefficient decompositions required per pair
    #need to figure out how to get the right coeffs per pair
    for species_ind1 in sort(collect(keys(species_dict)))
        for species_ind2 in sort(collect(keys(species_dict)))
            bonds["[$(species_ind1), $(species_ind2)]"] = "{p: $(p), r0: $(r0), xr: $(xr), xl: $(xl), pr: $(pr), pl: $(pl), rcut: $(rcut), radbasename: $(radbasename), maxn: $(maxn), coefficients: $(coeffs)}"
        end
    end

    return bonds
end

function export_ACE_V3(V3, species_dict)
    #we'll need to species dict here to also sort the coeffs per interation
    inner = V3.pibasis.inner[1]
    coeffs = V3.coeffs[1]

    #grabbing the basis functions, without interation specified
    groups = _basis_groups(inner, coeffs)

    #dumping all coeffs under "0", this is wrong
    #functions needs to respect the species_dict
    functions = Dict("0" => [])
    for group in groups
        push!(functions["0"],_write_group(group))
    end
    #HOW DO WE GET NSPECIES?
    return functions
end

## TO DO
#-  quotations might be an issue, C library allows to export with/without quotes,
#   Julia does not, maybe we can ask Yury to include them everywhere
#-  grouping of coefficients per interaction
#-  2B, correct radbasename and required params, as well as maybe repulsive core?
#-  embeddings, if it's linear what do we do?

export_ACE("test.yace", IP)


text:
"""
`SolversArgs` module aggregates the methods necessary for making solvers
available to use with `GuillotineModels.CommandLine` methods (especially
`run`).

The package `Requires.jl` is used to implement
`GuillotineModels.SolversArgs.empty_configured_model(Val(:SOME_SOLVER))`
only if `SOME_SOLVER` package is loaded.

The solver-specific methods this module implements for each supported solver
are: `GuillotineModels.Utilities.Args.{accepted_arg_list,
throw_if_incompatible_options}`, and `GuillotineModels.CommandLine.SolversArgs.
empty_configured_model`. The `GuillotineModels.Utilities` methods are always
available, not just when the respective solver package was imported.

The supported solvers are `CPLEX`, `Gurobi`, `Cbc`, and `GLPK`.
Supporting new solvers is just a question of implementing a specialization
of the three methods mentioned above for the specific solver (see the
code for examples).
"""
module SolversArgs

export empty_configured_model

import TimerOutputs.@timeit
import ...TIMER
import Requires.@require
import JuMP

using ArgParse

import ...Utilities
import ...Utilities.Args: Arg
import ...Utilities.Args: accepted_arg_list, throw_if_incompatible_options

# CPLEX is not define in this scope, so we return a symbol with the
# name of the constant and use `getfield` inside the scope CPLEX
# is defined.
const CPLEX_NAME2CODE_LP_ALG = Dict{String, Symbol}(
	"CPX_ALG_AUTOMATIC" => :CPX_ALG_AUTOMATIC,
	"auto" => :CPX_ALG_AUTOMATIC,
	"automatic" => :CPX_ALG_AUTOMATIC,
	"0" => :CPX_ALG_AUTOMATIC,
	"CPX_ALG_PRIMAL" => :CPX_ALG_PRIMAL,
	"primal" => :CPX_ALG_PRIMAL,
	"1" => :CPX_ALG_PRIMAL,
	"CPX_ALG_DUAL" => :CPX_ALG_DUAL,
	"dual" => :CPX_ALG_DUAL,
	"2" => :CPX_ALG_DUAL,
	"CPX_ALG_NET" => :CPX_ALG_NET,
	"net" => :CPX_ALG_NET,
	"network" => :CPX_ALG_NET,
	"3" => :CPX_ALG_NET,
	"CPX_ALG_BARRIER" => :CPX_ALG_BARRIER,
	"barrier" => :CPX_ALG_BARRIER,
	"4" => :CPX_ALG_BARRIER,
	"CPX_ALG_SIFTING" => :CPX_ALG_SIFTING,
	"sifting" => :CPX_ALG_SIFTING,
	"5" => :CPX_ALG_SIFTING,
	"CPX_ALG_CONCURRENT" => :CPX_ALG_CONCURRENT,
	"concurrent" => :CPX_ALG_CONCURRENT,
	"6" => :CPX_ALG_CONCURRENT,
	"CPX_ALG_CONCURRENT" => :CPX_ALG_CONCURRENT,
	"parallel" => :CPX_ALG_CONCURRENT,
	"7" => :CPX_ALG_CONCURRENT
)
function __init__()
	#function cplex_empty_configured_model(p_args)
	@require CPLEX="a076750e-1247-5638-91d2-ce28b192dca0" begin
	@timeit TIMER function empty_configured_model(::Val{:CPLEX}, p_args)
		scrind_value = p_args["no-output"] ? CPLEX.CPX_OFF : CPLEX.CPX_ON
		root_relax_method = CPLEX_NAME2CODE_LP_ALG[p_args["root-relax-method"]]
		lp_method = CPLEX_NAME2CODE_LP_ALG[p_args["LP-method"]]
		# https://www.ibm.com/support/pages/cplex-performance-tuning-linear-programs
		configuration = Pair{String, Any}[
			"CPX_PARAM_EPGAP" => 1e-6,
			"CPX_PARAM_TILIM" => p_args["time-limit"],
			"CPX_PARAM_STARTALG" => getfield(CPLEX, root_relax_method),
			"CPX_PARAM_LPMETHOD" => getfield(CPLEX, lp_method),
			# "Sifting is a simple form of column generation well suited for models
			# where the number of variables dramatically exceeds the number of
			# constraints."
			#"CPX_PARAM_STARTALG" => CPLEX.CPX_ALG_SIFTING,
			#"CPX_PARAM_LPMETHOD" => CPLEX.CPX_ALG_SIFTING,
			# "The barrier method tends to work well on problems where the product
			# of the constraint matrix multiplied by its transpose is sparse. "
			#"CPX_PARAM_STARTALG" => CPLEX.CPX_ALG_BARRIER,
			#"CPX_PARAM_LPMETHOD" => CPLEX.CPX_ALG_BARRIER,
			#"CPX_PARAM_STARTALG" => CPLEX.CPX_ALG_NET,
			#"CPX_PARAM_LPMETHOD" => CPLEX.CPX_ALG_NET,
			#"CPX_PARAM_BARDISPLAY" => 2, # 2 == diagnostic information level
			#"CPX_PARAM_SIMDISPLAY" => 2, # 2 == diagnostic information level
			# For the LPs of the iterative pricing of PPG2KP we need to avoid
			# numerical instability problems.
			#"CPX_PARAM_BARALG" => 1,
			# "[...] the computation time for the simplex method depends more on
			# the number of constraints than the number of variables."
			# And the dual is the opposite (CPLEX.CPX_ALG_DUAL).
			#"CPX_PARAM_STARTALG" => CPLEX.CPX_ALG_PRIMAL,
			#"CPX_PARAM_LPMETHOD" => CPLEX.CPX_ALG_PRIMAL,
			#"CPX_PARAM_PPRIIND" => CPLEX.CPX_PPRIIND_FULL, # pricing inside simplex
			#"CPX_PARAM_PERIND" => CPLEX.CPX_ON, # start using perturbations
			#"CPX_PARAM_PERLIM" => 1000000, # num degenerate iters until perturbation
			# Group parameter to help with numerical instability without the need
			# of fine-tuning.
			#"CPX_PARAM_NUMERICALEMPHASIS" => CPLEX.CPX_ON,
			#"CPX_PARAM_BAREPCOMP" => 1e-10,
			#"CPX_PARAM_EPMRK" => 0.9, # last measure against numerical instability
			"CPX_PARAM_SCRIND" => scrind_value,
			"CPX_PARAM_THREADS" => p_args["threads"],
			"CPX_PARAM_RANDOMSEED" => p_args["seed"]
		]
		if p_args["no-tricks"]
			extra_conf = [
				"CPX_PARAM_PROBE" => -1
				, "CPX_PARAM_HEURFREQ" => -1
				, "CPX_PARAM_REPEATPRESOLVE" => -1
				, "CPX_PARAM_DIVETYPE" => 1
			]
			configuration = [configuration; extra_conf]
		end
		raw_parameters = eval(
			Meta.parse(p_args["raw-parameters"])
		) :: Vector{Pair{String, Any}}
		configuration = [configuration; raw_parameters]
		model = JuMP.direct_model(CPLEX.Optimizer())
		JuMP.set_optimizer_attributes(model, configuration...)
		model
	end
	end # @require CPLEX

	@require Gurobi="2e9cd046-0924-5485-92f1-d5272153d98b" begin
	@timeit TIMER function empty_configured_model(::Val{:Gurobi}, p_args)
		# Note: without the type prefix below some values are casted and Gurobi.jl
		# fails to find the correct method to call (it differs from Int and Float).
		configuration = Pair{String, Any}[
			"Method" => p_args["LP-method"],
			# "NodeMethod" => 2, # use barrier for MIP non-root nodes
			"Threads" => p_args["threads"],
			"OutputFlag" => p_args["no-output"] ? 0 : 1,
			"TimeLimit" => p_args["time-limit"],
			"Seed" => p_args["seed"],
			"MIPGap" => 1e-6
		]
		# TODO: implement a no-trick flag for Gurobi?
		raw_parameters = eval(
			Meta.parse(p_args["raw-parameters"])
		) :: Vector{Pair{String, Any}}
		configuration = [configuration; raw_parameters]
		model = JuMP.direct_model(Gurobi.Optimizer(Gurobi.Env()))
		JuMP.set_optimizer_attributes(model, configuration...)
		model
	end
	end # @require Gurobi

	@require Cbc="9961bab8-2fa3-5c5a-9d89-47fab24efd76" begin
	@timeit TIMER function empty_configured_model(::Val{:Cbc}, p_args)
		configuration = Pair{String, Any}[
			# TODO: discover if the barrier algorithm is being correctly used.
			# The `barrier` command-line "parameter" is in fact an action/command
			# that solves the model relaxation with the barrier algorithm, and
			# not an argument for which linear algorithm the branch-and-bound
			# should use in the root node. The `logLevel` parameter, for example,
			# does not work with barrier, which will output even if `logLevel` is
			# set to zero.
			"logLevel" => p_args["no-output"] ? 0 : 1
			, "threads" => p_args["threads"]
			, "randomCbcSeed" => p_args["seed"]
			, "seconds" => p_args["time-limit"]
			, "ratioGap" => 1e-6
			, "barrier" => ""
		]
		# TODO: implement a no-trick flag for Cbc?
		raw_parameters = eval(
			Meta.parse(p_args["raw-parameters"])
		) :: Vector{Pair{String, Any}}
		configuration = [configuration; raw_parameters]
		model = JuMP.Model(JuMP.with_optimizer(Cbc.Optimizer))
		JuMP.set_optimizer_attributes(model, configuration...)
		model
	end
	end # @require Cbc

	@require GLPK="60bf3e95-4087-53dc-ae20-288a0d20c6a6" begin
	@timeit TIMER function empty_configured_model(::Val{:GLPK}, p_args)
		configuration = Pair{String, Any}[
			# TODO: find how configure: # of threads, tolerance gap, barrier, etc...
			# TODO: specially: does GLPK have a RANDOMSEED option?!
			"msg_lev" => p_args["no-output"] ? GLPK.MSG_OFF : GLPK.MSG_ON
			, "tm_lim" => p_args["time-limit"] * 1000 # is milliseconds not seconds
		]
		# TODO: implement a no-trick flag for GLPK? does GLPK even have tricks?
		raw_parameters = eval(
			Meta.parse(p_args["raw-parameters"])
		) :: Vector{Pair{String, Any}}
		configuration = [configuration; raw_parameters]
		model = JuMP.direct_model(GLPK.Optimizer())
		JuMP.set_optimizer_attributes(model, configuration...)
		model
	end
	end # @require GLPK
end

function Utilities.Args.accepted_arg_list(::Val{:CPLEX}) :: Vector{Arg}
	return [
		Arg(
			"root-relax-method", "automatic",
			"Changes CPX_PARAM_LPMETHOD. Accept simplified single-word names (primal, dual, etc...), CPLEX names (CPX_ALG_PRIMAL, CPX_ALG_DUAL, ...), and codes (1, 2, ...)."
		),
		Arg(
			"LP-method", "automatic",
			"Changes CPX_PARAM_STARTALG/CPXPARAM_MIP_Strategy_StartAlgorithm. Accept simplified single-word names (primal, dual, etc...), CPLEX names (CPX_ALG_PRIMAL, CPX_ALG_DUAL, ...), and codes (1, 2, ...)."
		),
		Arg(
			"threads", 1,
			"The value of CPXPARAM_Threads. If zero and no callbacks is the number os cores, if zero and callbacks is one. If a positive number, is that number of cores."
		),
		Arg(
			"time-limit", 31536000.0,
			"BROKEN, DO NOT USE, ALWAYS OVERWRITTEN BY `--generic-time-limit`. The value of CPX_PARAM_TILIM. Depends on CPX_PARAM_CLOCKTYPE which is by default wall clock time (not CPU time). Time limit in seconds for solver B&B (not root solving). Our default is one year."
		),
		Arg(
			"seed", 1,
			"The value of CPX_PARAM_RANDOMSEED. The random seed used by CPLEX."
		),
		Arg(
			"no-tricks", false,
			"Set CPX_PARAM_PROBE, CPX_PARAM_HEURFREQ, and CPX_PARAM_REPEATPRESOLVE to -1. Also, set CPX_PARAM_DIVETYPE to 1. Basically, disable many tricks used by CPLEX to speedup the solving (but that sometimes have the opposite effect)."
		),
		Arg(
			"no-output", false,
			"Set CPX_PARAM_SCRIND to false. Disables the solver output."
		),
		Arg(
			"raw-parameters", "Pair{String, Any}[]",
			"A string of Julia code to be evaluated to CPLEX parameters. For example: 'Pair{String, Any}[\"CPX_PARAM_SCRIND\" => CPLEX.CPX_OFF]' will have the same effect as --CPLEX-no-output. Obviously, this is a security breach."
		)
	]
end

function Utilities.Args.throw_if_incompatible_options(::Val{:CPLEX}, p_args)
	warn_about_time_limit(Val(:CPLEX), p_args)
	valid_lp_alg_names = sort!(collect(keys(CPLEX_NAME2CODE_LP_ALG)))
	lpm = "LP-method"
	rrm = "root-relax-method"
	Utilities.throw_if_unrecognized(lpm, p_args[lpm], valid_lp_alg_names)
	Utilities.throw_if_unrecognized(rrm, p_args[rrm], valid_lp_alg_names)
end

function Utilities.Args.accepted_arg_list(::Val{:Gurobi}) :: Vector{Arg}
	return [
		Arg(
			"LP-method", -1,
			"Gurobi has a parameter called 'Method' that defines the algorithm used to solve continuous models (including MIP root node continuous relaxations). The options (described in Gurobi 9.0 reference) are: automatic (-1), primal simplex (0), dual simplex (1), barrier (2), concurrent (3), deterministic concurrent (4), deterministic concurrent simplex (5)."
		),
		Arg(
			"threads", 1,
			"Number of threads for all Gurobi parallelizable algorithms. Zero is automatic, probably the number of cores but may be fewer. If a positive number, is that number of cores."
		),
		Arg(
			"time-limit", 31536000.0,
			"BROKEN, DO NOT USE, ALWAYS OVERWRITTEN BY `--generic-time-limit`. Set Gurobi parameter: TimeLimit. Total time limit in seconds. Our default is one year."
		),
		Arg(
			"seed", 1,
			"The random seed used by Gurobi. Our default (one) is different from Gurobi (which is zero)."
		),
		Arg(
			"no-output", false,
			"Set OutputFlag to zero. Disables the solver output."
		),
		Arg(
			"raw-parameters", "Pair{String, Any}[]",
			"A string of Julia code to be evaluated to Gurobi parameters. For example: 'Pair{String, Any}[\"OutputFlag\" => 0]' will have the same effect as --Gurobi-no-output. Obviously, this is a security breach."
		)
	]
end

function Utilities.Args.throw_if_incompatible_options(::Val{:Gurobi}, p_args)
	warn_about_time_limit(Val(:Gurobi), p_args)
	# TODO: should we throw if the option does not exist?
end

function Utilities.Args.accepted_arg_list(::Val{:Cbc}) :: Vector{Arg}
	return [
		Arg(
			"threads", 1,
			"Number of threads for \"parallel branch-and-bound\"."
		),
		Arg(
			"time-limit", 31536000.0,
			"BROKEN, DO NOT USE, ALWAYS OVERWRITTEN BY `--generic-time-limit`. Set Cbc parameter: seconds. Total time limit in seconds (? not very well documented)."
		),
		Arg(
			"seed", 1,
			"Set Cbc parameter: randomCbcSeed. Our default is 1 (different from Cbc default that is -1)."
		),
		Arg(
			"no-output", false,
			"Set logLevel to zero. Disables the solver output."
		),
		Arg(
			"raw-parameters", "Pair{String, Any}[]",
			"A string of Julia code to be evaluated to Cbc parameters. For example: 'Pair{String, Any}[\"logLevel\" => 0]' will have the same effect as --Cbc-no-output. Obviously, this is a security breach. The complete list of parameters can be found by running the cbc executable and typing ? at the prompt."
		)
	]
end

function Utilities.Args.throw_if_incompatible_options(::Val{:Cbc}, p_args)
	warn_about_time_limit(Val(:Cbc), p_args)
	# TODO: should we throw if the option does not exist?
end

function Utilities.Args.accepted_arg_list(::Val{:GLPK}) :: Vector{Arg}
	return [
		Arg(
			"time-limit", 2097152,
			"BROKEN, DO NOT USE, ALWAYS OVERWRITTEN BY `--generic-time-limit`. Set GLPK parameter: tm_lim. The original parameter is in milliseconds, but to keep it similar to the other solvers this option takes seconds. To set this parameter with milliseconds precision use the --raw-parameter option. The default is a little over 24 days because GLPK uses a Int32 for milliseconds and 24 days is close to the maximum time-limit supported."
		),
		Arg(
			"no-output", false,
			"Set msg_lev to GLPK.MSG_OFF. Disables the solver output."
		),
		Arg(
			"raw-parameters", "Pair{String, Any}[]",
			"A string of Julia code to be evaluated to GLPK parameters. For example: 'Pair{String, Any}[\"msg_lev\" => GLPK.MSG_OFF]' will have the same effect as --GLPK-no-output. Obviously, this is a security breach. If you find the complete list of GLPK parameters please send it to me (henriquebecker91@gmail.com)."
		)

	]
end

function warn_about_time_limit(::Val{SOLVER_SYM}, p_args) where {SOLVER_SYM}
	if haskey(p_args, "time-limit")
		tl = p_args["time-limit"]
		args = Utilities.Args.accepted_arg_list(Val(SOLVER_SYM))
		tl_arg_idx = findfirst(arg -> arg.name == "time-limit", args)
		if !isnothing(tl_arg_idx) && args[tl_arg_idx].default != tl
			@warn "--$SOLVER_SYM-time-limit is BROKEN, it is always overwritten" *
				" by --generic-time-limit even when this one is not passed (i.e., it" *
				" is overwritten by its default value in such case)."
		end
	end

	return
end

function Utilities.Args.throw_if_incompatible_options(::Val{:GLPK}, p_args)
	warn_about_time_limit(Val(:GLPK), p_args)
	# TODO: should we throw if the option does not exist?
end

"""
    empty_configured_model(::Val{T}, p_args) :: JuMP-like Model

Creates an empty (no variables or constraints) and configured (`p_args` is used
to set the parameters of an attached solver) model of the solver `::Val{T}`
(for example, passing `::Val{:CPLEX}` will call the specialized method for
CPLEX solver).

The options available for `p_args` are available in `GuillotineModels.
Utilities.Args.accepted_arg_list(::Val{SOLVER_PACKAGE_SYMBOL})`.
"""
function empty_configured_model(
	::Val{T}, p_args
) where {T}
	@error(
		"Solver " * string(T) * " is not supported (i.e., there is not" *
		" an implementation of Julia method `GuillotineModels.CommandLine" *
		".SolversArgs.empty_configured_model(::Val{Symbol(\"$(string(T))\")}," *
		"p_args)` for it) or the solver package was not imported before this" *
		"method was called."
	)
end

function empty_configured_model(
	::Val{:NoSolver}, p_args
) where {T}
	return JuMP.Model()
end

function Utilities.Args.accepted_arg_list(::Val{:NoSolver}) :: Vector{Arg}
	return Arg[
		Arg(
			"no-output", false,
			"Provided to avoid breaking the assumption every solver provides this" *
			" flag. Does notthing, as NoSolver should not print anything anyway."
		),
	]
end

function Utilities.Args.throw_if_incompatible_options(::Val{:NoSolver}, p_args)
	warn_about_time_limit(Val(:NoSolver), p_args)
	return nothing
end

end # module


text:
using Test
using GatherShot

@testset "reset mutated file" begin
    relpath = "sub/file.txt"
    mktempdir() do parent
        mkdir(joinpath(parent, dirname(relpath)))
        write(joinpath(parent, relpath), "original")
        mktempdir() do child
            mkdir(joinpath(child, dirname(relpath)))
            other = joinpath(child, relpath)
            write(other, "mutated")

            @test readline(other) == "mutated"
            GatherShot.reset(parent, child, other)
            @test readline(other) == "original"
        end
    end
end

@testset "new xml is one more than old xml" begin
    @test GatherShot.next_log(["f017.xml", "g001.xml"]) == 18
end


text:
struct Fda
    T::Array{Float64}
    P::Array{Float64}
    Tcenters::Array{Float64}
    eig::Vector{Float64}
    sstot::Number
    W::Matrix{Float64}
    xmeans::Vector{Float64}
    lev::AbstractVector
    ni::AbstractVector
end

"""
    fda(X, y; nlv, pseudo = false)
Factorial discriminant analysis (FDA).
* `X` : X-data.
* `y` : y-data (class membership).
* `nlv` : Nb. discriminant components.
* `pseudo` : If true, a MP pseudo-inverse is used (instead
    of a usual inverse) for inverting W.

Eigen factorization of Inverse(W) * B.

The functions maximize the compromise p'Bp / p'Wp, i.e. max p'Bp with
constraint p'Wp = 1. Vectors p (columns of P) are the linear discrimant
coefficients "LD".

`X` is internally centered.
"""
function fda(X, y; nlv, pseudo = false)
    fda!(copy(X), y; nlv = nlv, pseudo = pseudo)
end

function fda!(X, y; nlv, pseudo = false)
    X = ensure_mat(X)
    n, p = size(X)
    xmeans = colmeans(X)
    center!(X, xmeans)
    res = matW(X, y)
    lev = res.lev
    nlev = length(lev)
    ni = res.ni
    res.W .= res.W * n / (n - nlev)
    zres = matB(X, y)
    !pseudo ? Winv = LinearAlgebra.inv!(cholesky!(Hermitian(res.W))) : Winv = pinv(res.W)
    # Winv * B is not symmetric
    fm = eigen!(Winv * zres.B; sortby = x -> -abs(x))
    nlv = min(nlv, n, p, nlev - 1)
    P = fm.vectors[:, 1:nlv]
    eig = fm.values
    P = real.(P)
    eig = real.(eig)
    sstot = sum(eig)
    norm_P = sqrt.(diag(P' * res.W * P))
    scale!(P, norm_P)
    T = X * P
    Tcenters = zres.ct * P
    Fda(T, P, Tcenters, eig, sstot, res.W, xmeans, lev, ni)
end

"""
    fdasvd(X, y; nlv, pseudo = false)
Factorial discriminant analysis (FDA).
* `X` : X-data.
* `y` : Univariate class membership.
* `nlv` : Nb. discriminant components.
* `pseudo` : If true, a MP pseudo-inverse is used (instead
    of a usual inverse) for inverting W.

Weighted SVD factorization of the matrix of the class centers.

`X` is internally centered.
"""
function fdasvd(X, y; nlv, pseudo = false)
    fdasvd!(copy(X), y; nlv = nlv, pseudo = pseudo)
end

function fdasvd!(X, y; nlv, pseudo = false)
    X = ensure_mat(X)
    n, p = size(X)
    xmeans = colmeans(X)
    center!(X, xmeans)
    res = matW(X, y)
    lev = res.lev
    nlev = length(lev)
    ni = res.ni
    res.W .= res.W * n / (n - nlev)
    !pseudo ? Winv = inv(res.W) : Winv = pinv(res.W)
    ct = aggstat(X, y; fun = mean).res
    Ut = cholesky!(Hermitian(Winv)).U'
    Zct = ct * Ut
    nlv = min(nlv, n, p, nlev - 1)
    fm = pcasvd(Zct, ni; nlv = nlv)
    Pz = fm.P
    Tcenters = Zct * Pz
    eig = (fm.sv).^2
    sstot = sum(eig)
    P = Ut * Pz[:, 1:nlv]
    T = X * P
    Tcenters = ct * P
    Fda(T, P, Tcenters, eig, sstot, res.W, xmeans, lev, ni)
end

"""
    summary(object::Fda, X)
Summarize the fitted model.
* `object` : The fitted model.
* `X` : The X-data that was used to fit the model.
"""
function Base.summary(object::Fda)
    nlv = size(object.T, 2)
    eig = object.eig[1:nlv]
    pvar =  eig ./ sum(object.eig)
    cumpvar = cumsum(pvar)
    explvar = DataFrame(lv = 1:nlv, var = eig, pvar = pvar,
        cumpvar = cumpvar)
    (explvar = explvar,)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Readline"
version = v"8.0.4"

# Collection of sources required to build Readline
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/readline/readline-$(version.major).$(version.minor).tar.gz",
                  "e339f51971478d369f8a053a330a190781acb9864cf4c541060f12078948e461"),
    FileSource("https://ftp.gnu.org/gnu/readline/readline-$(version.major).$(version.minor)-patches/readline80-001",
               "d8e5e98933cf5756f862243c0601cb69d3667bb33f2c7b751fe4e40b2c3fd069"),
    FileSource("https://ftp.gnu.org/gnu/readline/readline-$(version.major).$(version.minor)-patches/readline80-002",
               "36b0febff1e560091ae7476026921f31b6d1dd4c918dcb7b741aa2dad1aec8f7"),
    FileSource("https://ftp.gnu.org/gnu/readline/readline-$(version.major).$(version.minor)-patches/readline80-003",
               "94ddb2210b71eb5389c7756865d60e343666dfb722c85892f8226b26bb3eeaef"),
    FileSource("https://ftp.gnu.org/gnu/readline/readline-$(version.major).$(version.minor)-patches/readline80-004",
               "b1aa3d2a40eee2dea9708229740742e649c32bb8db13535ea78f8ac15377394c"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/readline-*/

atomic_patch -p0 ${WORKSPACE}/srcdir/readline80-001
atomic_patch -p0 ${WORKSPACE}/srcdir/readline80-002
atomic_patch -p0 ${WORKSPACE}/srcdir/readline80-003
atomic_patch -p0 ${WORKSPACE}/srcdir/readline80-004

export CPPFLAGS="-I${prefix}/include"
if [[ "${target}" == *-mingw* ]]; then
    # Only on Windows the library name embeds the ABI version
    export NCURSES_ABI_VER=6
    export LDFLAGS="-L${libdir}"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-curses
make -j${nproc} SHLIB_LIBS="-lncurses${NCURSES_ABI_VER}"
make install
install_license COPYING
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = Product[
    LibraryProduct(["libhistory", "libhistory8"], :libhistory),
    LibraryProduct(["libreadline", "libreadline8"], :libreadline),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "Ncurses_jll",
]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
function citation_path(pkg)
    bib_path = joinpath(pkg.source, "CITATION.bib")
    if isfile(bib_path)
        bib_path
    end
end

function get_citation(pkg)
    bib_path = citation_path(pkg)
    if !isnothing(bib_path)
        @debug "Reading CITATION.bib for $(pkg.name)"
        try
            import_bibtex(bib_path)
        catch e
            @warn("There was an error reading the CITATION.bib file for $(pkg.name)")
            @debug e
        end
    end
end

"""
    collect_citations(only_direct::Bool)

Collect the citations from all the dependencies in the current environment.
"""
function collect_citations(only_direct::Bool)
    @debug "Generating citation report for the current environment"
    deps = Pkg.dependencies()
    pkg_citations = Dict{String, DataStructures.OrderedDict{String,Entry}}()
    for pkg in values(deps)
        if only_direct && !pkg.is_direct_dep
            continue
        end
        c = get_citation(pkg)
        if !isnothing(c)
            push!(pkg_citations, pkg.name=>c)
        end
    end

    return pkg_citations
end

function bibliography(pkg_citations)
    merge(values(pkg_citations)...)
end

function cited_packages(pkg_citations)
    keys(pkg_citations)
end

function export_citations(filename, pkg_citations)
    citations = bibliography(pkg_citations)
    if isfile(filename)
        @warn "Overwriting $filename"
    end
    if isempty(citations)
        return nothing
    end
    export_bibtex(filename, citations)
    pkgs = join(cited_packages(pkg_citations), ", ", " and ")
    @info "A $filename file with the citations for $pkgs" *
    " was generated in the current working directory ($(pwd()))."
end

"""
    get_citations(; only_direct=false, filename="julia_citations.bib")

This function will create a .bib file with all the citations collected form
the CITATION.bib files corresponding to the dependecies of
the current active environment. Use `filename` to change the name of the
file. To include just the direct dependencies use `only_direct=true`
"""
function get_citations(;only_direct=false, filename="julia_citations.bib")
    pkg_citations = collect_citations(only_direct)

    if isempty(pkg_citations)
        @warn "No citations found in current environment"
    else
        export_citations(filename, pkg_citations)
    end

    return nothing
end


text:
using Statistics

Statistics._var(A::CuArray, corrected::Bool, mean, dims) =
    sum((A .- something(mean, Statistics.mean(A, dims=dims))).^2, dims=dims)/(prod(size(A)[dims])-corrected)

Statistics._var(A::CuArray, corrected::Bool, mean, dim::Integer) =
    sum((A .- something(mean, Statistics.mean(A, dims=dim))).^2, dims=dim)/(size(A)[dim]-corrected)

Statistics._var(A::CuArray, corrected::Bool, mean, ::Colon) =
    sum((A .- something(mean, Statistics.mean(A))).^2)/(length(A)-corrected)

Statistics._std(A::CuArray, corrected::Bool, mean, dims) =
    Base.sqrt.(Statistics.var(A; corrected=corrected, mean=mean, dims=dims))

Statistics._std(A::CuArray, corrected::Bool, mean, ::Colon) =
    Base.sqrt.(Statistics.var(A; corrected=corrected, mean=mean, dims=:))

# Revert https://github.com/JuliaLang/Statistics.jl/pull/25
Statistics._mean(A::CuArray, ::Colon)    = sum(A) / length(A)
Statistics._mean(f, A::CuArray, ::Colon) = sum(f, A) / length(A)
Statistics._mean(A::CuArray, dims)    = mean!(Base.reducedim_init(t -> t/2, +, A, dims), A)
Statistics._mean(f, A::CuArray, dims) = sum(f, A, dims=dims) / mapreduce(i -> size(A, i), *, unique(dims); init=1)


text:
"""
Given the number of rows and columns, return an NTuple{4,Float64} containing
`(width, height, vspace, hspace)`, where `width` and `height` are the
width and height of each subplot and `vspace` and `hspace` are the vertical
and horizonal spacing between subplots, respectively.
"""
function sizes(nr::Int, nc::Int, subplot_titles::Bool=false)
    # NOTE: the logic of this function was mostly borrowed from plotly.py
    dx = 0.2 / nc
    dy = subplot_titles ? 0.55 / nr : 0.3 / nr
    width = (1. - dx * (nc - 1)) / nc
    height = (1. - dy * (nr - 1)) / nr
vspace = nr == 1 ? 0.0 : (1 - height * nr) / (nr - 1)
    hspace = nc == 1 ? 0.0 : (1 - width * nc) / (nc - 1)
    width, height, vspace, hspace
end

function gen_layout(rows::Tuple{Vararg{Int}}, subplot_titles::Bool=false)

    x = 0.0  # start from left
    y = 1.0  # start from top

    out = Layout()

    subplot = 1
    nr = length(rows)
    for nc in rows
        w, h, dy, dx = sizes(nr, nc, subplot_titles)
        x = 0.0 # reset x as we start a new row
        for col in 1:nc

            out["xaxis$subplot"] = attr(domain=[x, x + w], anchor="y$subplot")
            out["yaxis$subplot"] = attr(domain=[y - h, y], anchor="x$subplot")

            x += nc == 1 ? 0.0 : w + dx
            subplot += 1
         end
         y -= nr == 1 ? 0.0 : h + dy
    end

    out

end

gen_layout(nr::Int, nc::Int, subplot_titles::Bool=false) = gen_layout(tuple(fill(nc, nr)), subplot_titles)

function handle_titles!(big_layout, sub_layout, ix::Int)
    # don't worry about it if the sub_layout doesn't have a title
    if _isempty(sub_layout.title)
        return big_layout
    end

    # check for title or title text
    text = sub_layout.title
    if text isa Dict
        text = get(text, :text, missing)
        if ismissing(text)
            return big_layout
        end
    end

    ann = attr(
        font_size=16,
        showarrow=false,
        text=text,
        x=mean(big_layout["xaxis$(ix).domain"]),
        xanchor="center",
        xref="paper",
        y=big_layout["yaxis$(ix).domain"][2],
        yanchor="bottom",
        yref="paper"
    )
    anns = get(big_layout.fields, :annotations, [])
    push!(anns, ann)
    big_layout[:annotations] = anns
    big_layout
end

# plots are 3d if any of their traces have a 3d type. This should flow down
# the methods as ordered here
_is3d(p::Plot) = any(_is3d, p.data)
_is3d(t::GenericTrace) = _is3d(t[:type])
_is3d(t::Symbol) = _is3d(string(t))
_is3d(s::AbstractString) = s in ["surface", "mesh3d", "scatter3d"]

# else (maybe if trace didn't have a :type field set)
_is3d(x::Any) = false

function _cat(rows::Tuple{Vararg{Int}}, ps::Plot...)
    copied_plots = Plot[copy(p) for p in ps]
    subplot_titles = any(map(x -> haskey(x.layout.fields, :title) ||
                                  haskey(x.layout.fields, "title"), ps))
    layout = gen_layout(rows, subplot_titles)

    for ix in 1:sum(rows)
        handle_titles!(layout, copied_plots[ix].layout, ix)
        layout["xaxis$ix"] = merge(copied_plots[ix].layout["xaxis"], layout["xaxis$ix"])
        layout["yaxis$ix"] = merge(copied_plots[ix].layout["yaxis"], layout["yaxis$ix"])

        if _is3d(copied_plots[ix])
            # need to remove (x|y)axis$ix and move their domains into scene$ix here
            xaxis = pop!(layout, "xaxis$(ix)")
            yaxis = pop!(layout, "yaxis$(ix)")
            layout["scene$ix"] = attr(
                domain = attr(
                    x = xaxis[:domain],
                    y = yaxis[:domain]
                )
            )
            for trace in copied_plots[ix].data
                trace["scene"] = "scene$ix"
            end
        else
            for trace in copied_plots[ix].data
                trace["xaxis"] = "x$ix"
                trace["yaxis"] = "y$ix"
            end
        end

    end

    Plot(vcat([p.data for p in copied_plots]...), layout)
end

_cat(nr::Int, nc::Int, ps::Plot...) = _cat(Tuple(fill(nc, nr)), ps...)

Base.hcat(ps::Plot...) = _cat(1, length(ps), ps...)
Base.vcat(ps::Plot...) = _cat(length(ps), 1,  ps...)

Base.hvcat(rows::Tuple{Vararg{Int}}, ps::Plot...) = _cat(rows, ps...)

# implementation of algorithms from plotly.py/packages/python/plotly/plotly/subplots.py

function _init_subplot_xy!(
        layout::Layout, secondary_y::Bool, domain::NamedTuple{(:x, :y)}, max_subplot_ids::Dict
    )
    # Get axis label and anchor
    x_count = max_subplot_ids["xaxis"] + 1
    y_count = max_subplot_ids["xaxis"] + 1

    # Compute x/y labels (the values of trace.xaxis/trace.yaxis
    x_label = "x$(x_count > 1 ? x_count : "")"
    y_label = "y$(y_count > 1 ? y_count : "")"

    # Anchor x and y axes to each other
    x_anchor, y_anchor = y_label, x_label

    # Build layout.xaxis/layout.yaxis containers
    xaxis_name = "xaxis$(x_count > 1 ? x_count : "")"
    yaxis_name = "yaxis$(y_count > 1 ? y_count : "")"

    #  copy over existing xaxis/yaxis objects, but delete titles
    x_axis = attr(;layout.xaxis...)
    y_axis = attr(;layout.yaxis...)
    x_axis.domain = domain.x
    x_axis.anchor = x_anchor
    y_axis.domain = domain.y
    y_axis.anchor = y_anchor

    pop!(x_axis, :title_text)
    pop!(y_axis, :title_text)
    layout[xaxis_name] = x_axis
    layout[yaxis_name] = y_axis

    subplot_refs = [
        SubplotRef(
            subplot_kind="xy",
            layout_keys=Symbol.([xaxis_name, yaxis_name]),
            trace_kwargs=attr(xaxis=x_label, yaxis=y_label)
        )
    ]

    if secondary_y
        y_count += 1
        secondary_yaxis_name = "yaxis$(y_count > 1 ? y_count : "")"
        secondary_y_label = "y$(y_count > 1 ? y_count : "")"
        push!(
            subplot_refs,
            SubplotRef(
                subplot_kind="xy",
                layout_keys=Symbol.([xaxis_name, secondary_yaxis_name]),
                trace_kwargs=attr(xaxis=x_label, yaxis=secondary_y_label),
            )
        )
        secondary_y_axis = attr(anchor=y_anchor, overlaying=y_label, side="right")
        layout[secondary_yaxis_name] = secondary_y_axis
    end

    max_subplot_ids["xaxis"] = x_count
    max_subplot_ids["yaxis"] = y_count

    subplot_refs
end

function _init_subplot_single!(
        layout::Layout, subplot_kind::String, domain::NamedTuple{(:x, :y)}, max_subplot_ids::Dict
    )
    count = max_subplot_ids[subplot_kind] + 1
    label = "$(subplot_kind)$(count > 1 ? count : "")"
scene = attr(domain=attr(;domain...))
    layout[label] = scene

    trace_key = subplot_kind in _subplot_prop_named_subplot ? "subplot" : subplot_kind
    subplot_ref = SubplotRef(
        subplot_kind=subplot_kind,
        layout_keys=[Symbol(label)],
        trace_kwargs=attr(;Dict(Symbol(trace_key) => label)...)
    )
    max_subplot_ids[subplot_kind] = count
    [subplot_ref]
end

function _init_subplot_domain!(domain::NamedTuple{(:x, :y)})
    [SubplotRef(
        subplot_kind="domain",
        layout_keys=[],
        trace_kwargs=attr(domain=attr(;domain...))
    )]
end

function get_subplotkind_from_trace_type(k::Symbol)
    schema = get_plotschema()

    trace_attrs = schema.traces[k][:attributes]
    trace_attr_names = keys(trace_attrs)
    if :xaxis in trace_attr_names && :yaxis in trace_attr_names
        return "xy"
    elseif :domain in trace_attr_names
        return "domain"
    elseif :geo in trace_attr_names
        return "geo"
    elseif :scene in trace_attr_names
        return "scene"
    elseif k === :splom
        return "xy"
    elseif :subplot in trace_attr_names
        return trace_attrs[:subplot][:dflt]
    else
        @error "Unknown subplot type for trace $k. Please open issue"
    end
end

function _init_subplot!(
        layout::Layout, subplot_kind::String, secondary_y::Bool,
        domain::NamedTuple{(:x, :y)}, max_subplot_ids::Dict
    )
    sksymbol = Symbol(subplot_kind)
    if sksymbol in _TRACE_TYPES
        subplot_kind = get_subplotkind_from_trace_type(sksymbol)
    end
    if subplot_kind == "xy"
        return _init_subplot_xy!(layout, secondary_y, domain, max_subplot_ids)
    elseif subplot_kind in _single_subplot_types
        return _init_subplot_single!(layout, subplot_kind, domain, max_subplot_ids)
    elseif subplot_kind == "domain"
        return _init_subplot_domain!(domain)
    end
    error("Unknown subplot_kind $subplot_kind")
end

function _configure_shared_axes!(layout::Layout, grid_ref::Matrix, specs::Matrix{<:_Maybe{Spec}}, x_or_y::String, shared::Union{Bool,String}, row_dir::Int)
    rows, cols = size(grid_ref)

    layout_key_ind = x_or_y == "x" ? 1 : 2

    rows_iter = collect(1:rows)
    row_dir < 0 && reverse!(rows_iter)

    function update_axis_matches!(
            first_axis_id, subplot_ref::SubplotRef, spec::Spec, remove_label::Bool
        )
        span = x_or_y == "x" ? spec.colspan : spec.rowspan

        if subplot_ref.subplot_kind == "xy" && span == 1
            if ismissing(first_axis_id)
                first_axis_name = subplot_ref.layout_keys[layout_key_ind]
                first_axis_id = replace(String(first_axis_name), "axis" => "")
            else
                axis_name = subplot_ref.layout_keys[layout_key_ind]
                axis_to_match = layout[axis_name]
                axis_to_match[:matches] = first_axis_id
                if remove_label
                    axis_to_match[:showticklabels] = false
                end
            end
        end
        return first_axis_id
    end

    if shared == "columns" || (x_or_y == "x" && shared == true)
        for c in 1:cols
            first_axis_id = missing
            ok_to_remove_label = x_or_y == "x"
            for r in rows_iter
                subplot_ref = grid_ref[r, c][1]
                spec = specs[r, c]
                first_axis_id = update_axis_matches!(
                    first_axis_id, subplot_ref, spec, ok_to_remove_label
                )
            end
        end
    elseif shared == "rows" || (x_or_y == "y" && shared == true)
        for r in rows_iter
            first_axis_id = missing
            ok_to_remove_label = x_or_y == "y"
            for c in 1:cols
                subplot_ref = grid_ref[r, c][1]
                spec = specs[r, c]
                first_axis_id = update_axis_matches!(
                    first_axis_id, subplot_ref, spec, ok_to_remove_label
                )
            end
        end
    elseif shared == "all"
        first_axis_id = missing
        for c in 1:cols, (ri, r) in enumerate(rows_iter)
            subplot_ref = grid_ref[r, c][1]
            spec = specs[r, c]
            ok_to_remove_label = let
                if x_or_y == "y"
                    c > 0
                else
                    row_dir > 0 ? ri > 0 : r < rows
                end
            end
            first_axis_id = update_axis_matches!(
                first_axis_id, subplot_ref, spec, ok_to_remove_label
            )
        end
    end
end

_build_subplot_title_annotations(::Missing, list_of_domains; kw...) = []

function _build_subplot_title_annotations(
        subplot_titles::Array{<:_Maybe{String}},
        list_of_domains::Vector{<:NamedTuple{(:x, :y)}};
        title_edge="top", offset=0
    )
    # If shared_axes is false (default) use list_of_domains
    # This is used for insets and irregular layouts
    # if not shared_xaxes and not shared_yaxes:
    x_dom = [dom.x for dom in list_of_domains]
    y_dom = [dom.y for dom in list_of_domains]
    subtitle_pos_x = []
    subtitle_pos_y = []

    if title_edge == "top"
        text_angle = 0
        xanchor = "center"
        yanchor = "bottom"

        for x_domains in x_dom
            push!(subtitle_pos_x, sum(x_domains) / 2)
        end
        for y_domains in y_dom
            push!(subtitle_pos_y, y_domains[2])
        end
        yshift = offset
        xshift = 0
    elseif title_edge == "bottom"
        text_angle = 0
        xanchor = "center"
        yanchor = "top"

        for x_domains in x_dom
            push!(subtitle_pos_x, sum(x_domains) / 2)
        end
        for y_domains in y_dom
            push!(subtitle_pos_y, y_domains[1])
        end
        yshift = -offset
        xshift = 0
    elseif title_edge == "right"
        text_angle = 90
        xanchor = "left"
        yanchor = "middle"

        for x_domains in x_dom
            push!(subtitle_pos_x, x_domains[2])
        end
        for y_domains in y_dom
            push!(subtitle_pos_y, sum(y_domains) / 2.0)
        end
        yshift = 0
        xshift = offset
    elseif title_edge == "left"
        text_angle = -90
        xanchor = "right"
        yanchor = "middle"

        for x_domains in x_dom
            push!(subtitle_pos_x, x_domains[1])
        end
        for y_domains in y_dom
            push!(subtitle_pos_y, sum(y_domains) / 2.0)
        end
        yshift = 0
        xshift = -offset
    else
        error("title_edge must be one of [top, bottom, left, right]")
    end

    plot_titles = []
    for index in 1:length(subplot_titles)
        title = subplot_titles[index]
        should_continue = (
            ismissing(title) || title == false || length(title) == 0 ||
            index > length(subtitle_pos_y)
        )
        if should_continue
            continue
        end
        annotation = attr(
            y=subtitle_pos_y[index], yref="paper", yanchor=yanchor,
            x=subtitle_pos_x[index], xref="paper", xanchor=xanchor,
            text=title, font_size=16,
            showarrow=false,
        )
        xshift != 0 && setindex!(annotation, xshift, :xshift)
        yshift != 0 && setindex!(annotation, yshift, :yshift)
        text_angle != 0 && setindex!(annotation, text_angle, :textangle)
        push!(plot_titles, annotation)
    end
    plot_titles
end

function Layout(sp::Subplots; kw...)
    layout = Layout(;kw...)

    @unpack_Subplots sp

    row_dir = start_cell == "top-left" ? -1 : 1

    col_seq = collect(1:cols) .- 1
    row_seq = collect(1:rows) .- 1
    row_dir < 0 && reverse!(row_seq)

    grid = [
        (
            sum(_widths[1:c]) + c * horizontal_spacing,
            sum(_heights[1:r]) + r * vertical_spacing
        )
        for r in row_seq, c in col_seq
    ]
    domains_grid = _Maybe{NamedTuple{(:x, :y)}}[missing for _ in 1:rows, __ in 1:cols]
    list_of_domains = NamedTuple{(:x, :y)}[]
    max_subplot_ids = _get_initial_max_subplot_ids()

    for r in 1:rows, c in 1:cols
        spec = specs[r, c]
        if ismissing(spec)
            continue
        end
        c_spanned = c + spec.colspan - 1
        r_spanned = r + spec.rowspan - 1

        c_spanned > cols && error("Some colspan value is too large for this subplot grid")
        r_spanned > rows && error("Some rowspan value is too large for this subplot grid")

        # grid x c_spanned -> x_domain
        x_s = grid[r, c][1]
        x_e = grid[r, c_spanned][1] + _widths[c_spanned] - spec.r
        x_domain = (max(0.0, x_s), min(1.0, x_e))

        # grid x r_spanned x row_dir -> yaxis_domain
        y_domain = let
            if row_dir > 1
                (grid[r, c][2] + spec.b, grid[r_spanned, c][2] + _heights[r_spanned] - spec.t)
            else
                (grid[r_spanned, c][2], grid[r, c][2] + _heights[end - r + 1] - spec.t)
            end
        end
        domain = (x = x_domain, y = y_domain)
        domains_grid[r, c] = domain

        push!(list_of_domains, domain)

        grid_ref[r, c] = _init_subplot!(layout, spec.kind, spec.secondary_y, domain, max_subplot_ids)

    end

    _configure_shared_axes!(layout, grid_ref, specs, "x", shared_xaxes, row_dir)
    _configure_shared_axes!(layout, grid_ref, specs, "y", shared_yaxes, row_dir)

    # handle insets

    insets_ref = ismissing(insets) ? missing : Any[missing for _ in 1:length(insets)]
    if !ismissing(insets)
        for (i_inset, inset) in enumerate(insets)
            r = inset.cell[1]
            c = inset.cell[2]

            0 <= r <= rows || error("Some `cell` out of range")
            0 <= c <= cols || error("Some `cell` out of range")

            x_s = grid[r, c][1] + inset.l * _widths[c]
            x_e = inset.w == "to_end" ? grid[r, c][1] + _widths[c] : (x_s + inset.w * _widths[c])
            x_domain = (x_s, x_e)

            y_s = grid[r, c][2] + inset.b * _heights[end - r + 1]
            y_e = inset.h == "to_end" ? (grid[r, c][2] + _heights[end - r + 1]) : (y_s + inset.h * _heights[end - r + 1])
            y_domain = (y_s, y_e)
            domain = (x = x_domain, y = y_domain)

            push!(list_of_domains, domain)

            insets_ref[i_inset] = _init_subplot!(
                layout, inset.kind, false, domain, max_subplot_ids
            )
        end
    end

    plot_title_annotations = _build_subplot_title_annotations(
        subplot_titles, list_of_domains
    )
    layout[:annotations] = plot_title_annotations

    if !ismissing(column_titles)
        domains_list = let
            if row_dir > 0
                [domains_grid[end, c] for c in 1:cols]
            else
                [domains_grid[1, c] for c in 1:cols]
            end
        end
        append!(
            layout[:annotations],
            _build_subplot_title_annotations(column_titles, domains_list)
        )
    end

    if !ismissing(row_titles)
        domains_list = [domains_grid[r, end] for r in 1:rows]
        append!(
            layout[:annotations],
            _build_subplot_title_annotations(row_titles, domains_list, title_edge="right")
        )
    end

    if !ismissing(x_title)
        domains_list = [(x = (0, max_width), y = (0, 1))]
        append!(
            layout[:annotations],
            _build_subplot_title_annotations(
                [x_title], domains_list, title_edge="bottom", offset=30
            )
        )
    end
    if !ismissing(y_title)
        domains_list = [(x = (0, 1), y = (0, 1))]
        append!(
            layout[:annotations],
            _build_subplot_title_annotations(
                [y_title], domains_list, title_edge="left", offset=40
            )
        )
    end
    layout.subplots = sp
    layout
end


text:
# This file contains code that was formerly part of Julia. License is MIT: http://julialang.org/license

using QuadGK, Test

@testset "quadgk" begin
    @test quadgk(cos, 0,0.7,1)[1] ≈ sin(1)
    @test quadgk(x -> exp(im*x), 0,0.7,1)[1] ≈ (exp(1im)-1)/im
    @test quadgk(x -> exp(im*x), 0,1im)[1] ≈ -1im*expm1(-1)
    @test isapprox(quadgk(cos, 0,BigFloat(1),order=40)[1], sin(BigFloat(1)),
                   atol=1000*eps(BigFloat))
    @test quadgk(x -> exp(-x), 0,0.7,Inf)[1] ≈ 1.0
    @test quadgk(x -> exp(x), -Inf,0)[1] ≈ 1.0
    @test quadgk(x -> exp(-x^2), -Inf,Inf)[1] ≈ sqrt(pi)
    @test quadgk(x -> [exp(-x), exp(-2x)], 0, Inf)[1] ≈ [1,0.5]
    @test quadgk(cos, 0,0.7,1, norm=abs)[1] ≈ sin(1)

    # Test a function that is only implemented for Float32 values
    cos32(x::Float32) = cos(20x)
    @test quadgk(cos32, 0f0, 1f0)[1]::Float32 ≈ sin(20f0)/20

    # test integration of a type-unstable function where the instability is only detected
    # during refinement of the integration interval:
    @test quadgk(x -> x > 0.01 ? sin(10(x-0.01)) : 1im, 0,1.01, rtol=1e-4, order=3)[1] ≈ (1 - cos(10))/10+0.01im rtol=1e-4
end

module Test19626
    using QuadGK
    using Test

    # Define a mock physical quantity type
    struct MockQuantity <: Number
        val::Float64
    end

    # Following definitions needed for quadgk to work with MockQuantity
    import Base: +, -, *, abs, isnan, isinf, isless, float
    +(a::MockQuantity, b::MockQuantity) = MockQuantity(a.val+b.val)
    -(a::MockQuantity, b::MockQuantity) = MockQuantity(a.val-b.val)
    *(a::MockQuantity, b::Number) = MockQuantity(a.val*b)
    abs(a::MockQuantity) = MockQuantity(abs(a.val))
    float(a::MockQuantity) = a
    isnan(a::MockQuantity) = isnan(a.val)
    isinf(a::MockQuantity) = isinf(a.val)
    isless(a::MockQuantity, b::MockQuantity) = isless(a.val, b.val)

    # isapprox only needed for test purposes
    Base.isapprox(a::MockQuantity, b::MockQuantity) = isapprox(a.val, b.val)

    # Test physical quantity-valued functions
    @test QuadGK.quadgk(x->MockQuantity(x), 0.0, 1.0, atol=MockQuantity(0.0))[1] ≈
        MockQuantity(0.5)
end

@testset "inference" begin
    @test @inferred(QuadGK.cachedrule(Float16, 3)) == (Float16[-0.96, -0.7744, -0.434, 0.0], Float16[0.1047, 0.269, 0.4014, 0.4504], Float16[0.555, 0.8896])
    @test @inferred(QuadGK.cachedrule(Complex{BigFloat}, 3)) isa NTuple{3,Vector{BigFloat}}
    @test @inferred(quadgk(x -> exp(-x^2), 0, Inf, rtol=1e-8)) isa Tuple{Float64,Float64}
    @test @inferred(quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)) isa Tuple{Float64,Float64}
    @test @inferred(quadgk(x -> 1, 0, 1im)) === (1.0im, 0.0)
    @test @inferred(quadgk(x -> sin(10x), 0,1))[1] ≈ (1 - cos(10))/10
end

@testset "gauss" begin
    x,w = gauss(10, -1, 1)
    x′,w′ = @inferred gauss(x->1, 10, -1, 1)
    @test x ≈ x′ ≈ [-0.9739065285171717, -0.8650633666889845, -0.6794095682990244, -0.4333953941292472, -0.14887433898163124, 0.14887433898163124, 0.4333953941292472, 0.6794095682990244, 0.8650633666889845, 0.9739065285171717]
    @test w ≈ w′ ≈ [0.06667134430868811, 0.14945134915058064, 0.2190863625159821, 0.26926671930999635, 0.29552422471475276, 0.29552422471475276, 0.26926671930999635, 0.2190863625159821, 0.14945134915058064, 0.06667134430868811]

    # Gauss–Hermite quadrature:
    # we don't support infinite intervals, but any sufficiently broad interval
    # for a decaying weight function should be equivalent to infinite (to any desired precision).
    xH,wH = @inferred gauss(x->exp(-x^2), 5, -10, 10, rtol=1e-10)
    @test sum(xH.^4 .* wH) ≈ quadgk(x -> x^4 * exp(-x^2), -Inf, Inf)[1] ≈ 3sqrt(π)/4
    @test xH ≈ [-2.020182870456085632929,  -0.9585724646138185071128, 0, 0.9585724646138185071128, 2.020182870456085632929]
    @test wH ≈ [0.01995324205904591320774, 0.3936193231522411598285, 0.9453087204829418812257, 0.393619323152241159829, 0.01995324205904591320774]

    x,w = gauss(200, 2, 3)
    x′,w′ = gauss(x->1, 200, 2, 3, rtol=1e-11)
    @test x ≈ x′
    @test w ≈ w′

    # check for underflow bug: https://discourse.julialang.org/t/nan-returned-by-gauss/68260
    x,w = gauss(1093)
    @test all(isfinite, x) && all(isfinite, w)
end

≅(x::Tuple, y::Tuple) = all(a -> isapprox(a[1],a[2]), zip(x,y))

@testset "inplace" begin
    I = [0., 0.]
    I′,E′ = quadgk!(I, 0, 1) do r,x
        r[1] = cos(100x)
        r[2] = sin(30x)
    end
    @test quadgk(x -> [cos(100x), sin(30x)], 0, 1) ≅ (I′,E′) ≅ ([-0.005063656411097513, 0.028191618337080532], 4.2100180879009775e-10)
    @test I === I′ # result is written in-place to I
end

# This is enough for allocation currently caused by the do-lambda in quadgk(...)
const smallallocbytes = 500

@testset "segbuf" begin
    # Should not need subdivision
    function id(x::Float64)::Float64
        1.0
    end
    # Should need subdivision
    function osc(x::Float64)::Float64
        (x - 0.3)^2 * sin(87(x + 0.07))
    end
    no_subdiv() = @timed quadgk(id, -1.0, 1.0)
    subdiv_alloc() = @timed quadgk(osc, -1.0, 1.0)
    segbuf = alloc_segbuf(size=1)
    subdiv_alloc_segbuf() = @timed quadgk(osc, -1.0, 1.0, segbuf=segbuf)
    no_subdiv() # warmup
    @test no_subdiv()[3] < smallallocbytes # [3] == .bytes starting in Julia 1.5
    subdiv_alloc() # warmup
    @test subdiv_alloc()[3] > smallallocbytes
    subdiv_alloc_segbuf() # warmup
    @test subdiv_alloc_segbuf()[3] < smallallocbytes
end


text:
# These are all the systems used for the general running of Dioramas
struct Timer <: System end

Overseer.requested_components(::Timer) = (TimingData,)

function Overseer.update(::Timer, m::AbstractLedger)
    for t in m[TimingData]
        nt = time()
        t.dtime = t.reversed ? -nt + t.time : nt - t.time
        t.time = nt
        t.frames += 1
    end
end

struct Sleeper <: System end
Overseer.requested_components(::Sleeper) = (TimingData, Canvas)

function Overseer.update(::Sleeper, m::AbstractLedger)
    sd = m[TimingData]
    @timeit sd[1].timer "swapping" swapbuffers(m[Canvas][1])
    curtime    = time()
    dt         = (curtime - sd[1].time)
    sleep_time = 1 / sd[1].preferred_fps - dt
    st         = sleep_time - 0.002
    if st > 0.0
        sleep(st)
    else
        yield()
    end
end

struct Resizer <: System end
Overseer.requested_components(::Resizer) = (Canvas, IOTarget)

function Overseer.update(::Resizer, m::AbstractLedger)
    c = singleton(m, Canvas)
    iofbo = singleton(m, IOTarget)

    fwh = c.framebuffer_size
    resize!(c, fwh)
    for c in components(m, RenderTarget)
        for rt in c
            resize!(rt, fwh)
        end
    end
    bind(iofbo)
    draw(iofbo)
    return clear!(iofbo)
end


text:
using Distributions

struct PIDController
  roll
  pitch
  yaw
  x
  y
  z
  trimThrust
end

struct LQRController
  gains
end

struct PIDDataDistributions
  rollP
  rollI
  rollD
  pitchP
  pitchI
  pitchD
  yawP
  yawI
  yawD
  xP
  xI
  xD
  yP
  yI
  yD
  zP
  zI
  zD
end

struct LQRDataDistributions
  x
  y
  z
  u
  v
  w
  p
  q
  r
  phi
  theta
  psi
end

const NUM_DATAPOINTS = 100000
const MOTOR_MAX_PWM = 2000

function bound(lowerLim,upperLim,value)
  return max(min(value,upperLim),lowerLim)
end

function setController(rollPid,pitchPid,yawPid,xPid,yPid,zPid)
  return PIDController(rollPid,pitchPid,yawPid,xPid,yPid,zPid,900,590)
end

function setController(gains)
  return LQRController(gains)
end

function setDistributions(bw,distType)
  if distType == "PID"
    return PIDDataDistributions(Distributions.Uniform(bw["roll"][1]...),
                                Distributions.Uniform(bw["roll"][2]...),
                                Distributions.Uniform(bw["roll"][3]...),
                                Distributions.Uniform(bw["pitch"][1]...),
                                Distributions.Uniform(bw["pitch"][2]...),
                                Distributions.Uniform(bw["pitch"][3]...),
                                Distributions.Uniform(bw["yaw"][1]...),
                                Distributions.Uniform(bw["yaw"][2]...),
                                Distributions.Uniform(bw["yaw"][3]...),
                                Distributions.Uniform(bw["x"][1]...),
                                Distributions.Uniform(bw["x"][2]...),
                                Distributions.Uniform(bw["x"][3]...),
                                Distributions.Uniform(bw["y"][1]...),
                                Distributions.Uniform(bw["y"][2]...),
                                Distributions.Uniform(bw["y"][3]...),
                                Distributions.Uniform(bw["z"][1]...),
                                Distributions.Uniform(bw["z"][2]...),
                                Distributions.Uniform(bw["z"][3]...))
  elseif distType == "LQR"
     return LQRDataDistributions(Distributions.Uniform(bw["x"]...),
                                 Distributions.Uniform(bw["y"]...),
                                 Distributions.Uniform(bw["z"]...),
                                 Distributions.Uniform(bw["p"]...),
                                 Distributions.Uniform(bw["q"]...),
                                 Distributions.Uniform(bw["r"]...),
                                 Distributions.Uniform(bw["u"]...),
                                 Distributions.Uniform(bw["v"]...),
                                 Distributions.Uniform(bw["w"]...),
                                 Distributions.Uniform(bw["phi"]...),
                                 Distributions.Uniform(bw["theta"]...),
                                 Distributions.Uniform(bw["psi"]...))
  end
end

function generatePoints(controller::PIDController,
                        distributions::PIDDataDistributions,
                        numPoints)
  # Simulate measured data
  data = [transpose([rand(distributions.rollP),
                     rand(distributions.rollI),
                     rand(distributions.rollD),
                     rand(distributions.pitchP),
                     rand(distributions.pitchI),
                     rand(distributions.pitchD),
                     rand(distributions.yawP),
                     rand(distributions.yawI),
                     rand(distributions.yawD),
                     rand(distributions.xP),
                     rand(distributions.xI),
                     rand(distributions.xD),
                     rand(distributions.yP),
                     rand(distributions.yI),
                     rand(distributions.yD),
                     rand(distributions.zP),
                     rand(distributions.zI),
                     rand(distributions.zD)])
          for i=1:numPoints]
  # Determine control for measured data
  label = [pidControl(p,controller) for p in data]
  return Dict("data" => data,
              "label" => label)
end

function generatePoints(controller::LQRController,
                        distributions::LQRDataDistributions,
                        numPoints)
  data = [[[rand(distributions.x),
            rand(distributions.y),
            rand(distributions.z),
            rand(distributions.u),
            rand(distributions.v),
            rand(distributions.w),
            rand(distributions.p),
            rand(distributions.q),
            rand(distributions.r),
            rand(distributions.phi),
            rand(distributions.theta),
            rand(distributions.psi)];
           zeros(12,1)]
          for i=1:numPoints]
  data = map(p->transpose(p),data)
  label = [lqrControl(p,controller) for p in data]
  return Dict("data" => data,
              "label" => label)
end

#FIXME ensure correct
# Assumes motors as such: forward left, forward right, back left, back right
# Coordinate system:
#   positive x right, y back, z down
#   positive pitch back, roll left, yaw right
function pidControl(pt,controller::PIDController)
  rollDesired =  pt[13]*controller.x[1] +
                 pt[14]*controller.x[2] +
                 pt[15]*controller.x[3]
  pitchDesired = pt[10]*controller.y[1] +
                 pt[11]*controller.y[2] +
                 pt[12]*controller.y[3]

  pitchControl = (pt[4] - pitchDesired)*controller.pitch[1] +
                 pt[5]*controller.pitch[2] +
                 pt[6]*controller.pitch[3]
  rollControl = (pt[1] - rollDesired)*controller.roll[1] +
                 pt[2]*controller.roll[2] +
                 pt[3]*controller.roll[3]
  thrustControl = pt[16]*controller.z[1] +
                  pt[17]*controller.z[2] +
                  pt[18]*controller.z[3] +
                  controller.trimThrust +
                  controller.baseThrust
  yawControl = pt[7]*controller.yaw[1] +
               pt[8]*controller.yaw[2] +
               pt[9]*controller.yaw[3]

  ff = open("./data/tmp.txt","a+")
  write(ff,"$thrustControl,$pitchControl,$rollControl,$yawControl\n")
  close(ff)

  lf = thrustControl+rollControl-pitchControl-yawControl
  rf = thrustControl-rollControl-pitchControl+yawControl
  lb = thrustControl+rollControl+pitchControl+yawControl
  rb = thrustControl-rollControl+pitchControl-yawControl
  return map(u->bound(-MOTOR_MAX_PWM,MOTOR_MAX_PWM,u),[u1 u2 u3 u4])
end

function lqrControl(pt,controller::LQRController)
  return controller.gains*pt[1:12]
end

function writePoints(f,p,controlType)
  data = p["data"]
  label = p["label"]
  if controlType == "PID"
    header = "[rollP rollI rollD pitchP pitchI pitchD yawP yawI yawD xP xI xD yP yI yD zP zI zD], [lf rf lb rb]\n"
  elseif controlType == "LQR"
    header = "[xd yd zd ud vd wd pd qd rd phid thetad psid x y z u v w p q r phi theta psi], [del_lat del_lon del_yaw del_thr]\n"
  end
  write(f,header)
  #for i in length(data)
  #  writedlm(f, [data label], ", ")
  #end
  writedlm(f, [data label], ", ")
end

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

token = "LQR"
if token == "test"
  f = open("./data/test.txt","w")

  rollPid = [-2,0,40*.017]
  pitchPid = [2,0,40*.017]
  yawPid = [23.22,.192,80*.017]
  xPid = [.0122,.0041,.026]
  yPid = [.0122,.0041,.026]
  zPid = [.2330,.069,.0188]

  controller = setController(rollPid,pitchPid,yawPid,xPid,yPid,zPid)

  rollBw = [[-pi pi],[-1000 1000],[-5 5]]
  pitchBw = [[-pi pi],[-1000 1000],[-5 5]]
  yawBw = [[-pi pi],[-1000 1000],[-5 5]]
  xBw = [[-2000 2000],[-4000 4000],[-3000 3000]]
  yBw = [[-2000 2000],[-4000 4000],[-3000 3000]]
  zBw = [[-2000 2000],[-4000 4000],[-3000 3000]]
  bw = Dict("roll" => rollBw,
            "pitch" => pitchBw,
            "yaw" => yawBw,
            "x" => xBw,
            "y" => yBw,
            "z" => zBw)

  distributions = setDistributions(bw,"PID")

  numPoints = NUM_DATAPOINTS

  points = generatePoints(controller,distributions,numPoints)
  writePoints(f,points,"PID")

  close(f)
elseif token == "LQR"
  f = open("./data/lqrData.txt","w")
  gains = [0.0020 1.0000 -0.0000 -0.0418 1.3856 -0.0000 1.1039 0.0001 -0.0000 4.6346 0.0395 -0.0000;
           -1.0000 0.0020 -0.0043 -1.3583 -0.0458 -0.0015 0.0004 1.1348 -0.0028 -0.0247 4.6926 -0.0028;
           -0.0035 0.0000 0.1773 -0.0070 -0.0001 -0.0188 -0.0000 -0.0010 1.0131 -0.0002 -0.0142 0.9842;
           0.0037 -0.0000 -0.9842 0.0147 0.0001 -0.9920 0.0000 0.0044 0.2049 0.0001 0.2483 0.1773]
  controller = setController(gains)

  xBw = [-5 5]
  yBw = [-5 5]
  zBw = [-5 5]
  uBw = [-5 5]
  vBw = [-5 5]
  wBw = [-5 5]
  pBw = [-7 7]
  qBw = [-7 7]
  rBw = [-7 7]
  phiBw = [-pi pi]
  thetaBw = [-pi pi]
  psiBw = [-pi pi]
  bw = Dict("x" => xBw,
            "y" => yBw,
            "z" => zBw,
            "u" => uBw,
            "v" => vBw,
            "w" => wBw,
            "p" => pBw,
            "q" => qBw,
            "r" => rBw,
            "phi" => phiBw,
            "theta" => thetaBw,
            "psi" => psiBw)

  distributions = setDistributions(bw,"LQR")

  numPoints = NUM_DATAPOINTS
  points = generatePoints(controller,distributions,numPoints)
  writePoints(f,points,"LQR")
  close(f)
elseif token == "PID"
end

max_stars_repo_name: umd-agrc/SimpleControlSim

max_issues_repo_name: umd-agrc/SimpleControlSim

max_forks_repo_name: umd-agrc/SimpleControlSim


text:
module OverEngineering

include("fields.jl")
include("field_operators.jl")

end # module Domain

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
# TODO: Make this fast
function expl1(a::AbstractSet, b::AbstractSet)
	d = length(union(a, b)) - length(intersect(a, b))
	return exp(-d)
end

expl1(a::AbstractVector, b::AbstractVector) = exp(-cityblock(a, b))


text:
using FFTW, Jets, JetPackTransforms, Test

@testset "DCT, 1D" begin
    A = JopDct(Float64, 128)
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs

    m = rand(domain(A))
    d = A*m
    d_check = dct(m)
    d_check[1] *= sqrt(2)
    d_check[:] *= sqrt(1/length(m))
    @test d ≈ d_check
end

@testset "DCT, 2D" begin
    A = JopDct(Float64, 128, 64)
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs

    m = rand(domain(A))
    d = A*m
    d_check = dct(m)
    d_check[1] *= 2*sqrt(2)
    d_check[:] *= sqrt(1/length(m))
    @test d_check ≈ d

    A = JopDct(Float64, 128, 64; dims=(1,))
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs

    m = rand(domain(A))
    d = A*m
    d_check = dct(m,1)
    d_check[1] *= sqrt(2)
    d_check[:] *= sqrt(1/size(m,1))
    @test d ≈ d_check

    A = JopDct(Float64, 128, 64; dims=(2,))
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs

    m = rand(domain(A))
    d = A*m
    d_check = dct(m,2)
    d_check[1] *= sqrt(2)
    d_check .*= sqrt(1/size(m,2))
    @test d ≈ d_check
end

@testset "DCT, 3D" begin
    A = JopDct(Float64, 128, 64, 32)
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs

    m = rand(domain(A))
    d = A*m
    d_check = dct(m)
    d_check[1] *= 3*sqrt(2)
    d_check[:] *= sqrt(1/length(m))
    @test d ≈ d_check

    for dims in ((1,2), (1,3), (2,1), (2,3), (1,), (2,), (3,))
        A = JopDct(Float64, 128, 64, 32; dims=dims)
        lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
        @test lhs ≈ rhs

        m = rand(domain(A))
        d = A*m
        d_check = dct(m,dims)
        d_check[1] *= length(dims)*sqrt(2)
        d_check[:] *= sqrt(1/mapreduce(dim->size(m,dim), *, dims))
        @test d ≈ d_check
    end
end


text:
# Different methods of generating PSF

using FFTW

field(amplitude, phase) = amplitude .* exp.( 1im * phase)

"""
psf(amplitude, phase) -> psfimage
psf(pupilfield) -> psfimage

Calculate psf for given amplitude and phase.

Examples
====
```jldoctest

julia> ap,_ = PhaseRetrieval.aperture(-1:.2:1,-1:.2:1,.8)
([0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], [NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN 1.0 NaN NaN NaN NaN NaN; NaN NaN NaN NaN 1.0 1.0 1.0 NaN NaN NaN NaN; NaN NaN NaN 1.0 1.0 1.0 1.0 1.0 NaN NaN NaN; NaN NaN NaN NaN 1.0 1.0 1.0 NaN NaN NaN NaN; NaN NaN NaN NaN NaN 1.0 NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN])

julia> psfimage = psf(ap)
[17.720305330423212 17.720305330423212 2.830830026003772 0.7990467619253477 1.7153703234534299 0.002318499057199975 0.8566413660014238 0.0023184990571999695 1.7153703234534299 0.7990467619253475 2.830830026003772; 17.720305330423212 17.720305330423212 2.830830026003772 0.7990467619253475 1.7153703234534299 0.0023184990571999695 0.8566413660014238 0.002318499057199975 1.7153703234534299 0.7990467619253477 2.830830026003772; 2.830830026003773 2.830830026003773 0.2240431494891658 5.881503709461158 4.671643508435733 0.0810140527710051 0.6181197482998197 0.0810140527710051 4.671643508435733 5.881503709461155 0.2240431494891658; 0.7990467619253475 0.7990467619253477 5.881503709461158 11.063720826850961 3.6825070656623606 0.6902785321094302 4.96008586865757 0.6902785321094302 3.6825070656623606 11.063720826850961 5.881503709461158; 1.7153703234534294 1.7153703234534299 4.671643508435735 3.6825070656623615 0.5365498748309357 19.64548752112056 38.22662768629444 19.64548752112056 0.5365498748309356 3.682507065662362 4.671643508435733; 0.00231849905719993 0.002318499057199919 0.08101405277100532 0.6902785321094301 19.64548752112056 78.45538081840574 118.33852533074673 78.45538081840574 19.64548752112056 0.6902785321094301 0.08101405277100517; 0.8566413660014239 0.8566413660014239 0.6181197482998197 4.96008586865757 38.22662768629444 118.33852533074673 169.0 118.33852533074673 38.22662768629444 4.96008586865757 0.6181197482998197; 0.002318499057199919 0.00231849905719993 0.08101405277100517 0.6902785321094301 19.64548752112056 78.45538081840574 118.33852533074673 78.45538081840574 19.64548752112056 0.6902785321094301 0.08101405277100532; 1.7153703234534299 1.7153703234534294 4.671643508435733 3.682507065662362 0.5365498748309356 19.64548752112056 38.22662768629444 19.64548752112056 0.5365498748309357 3.6825070656623615 4.671643508435735; 0.7990467619253477 0.7990467619253475 5.881503709461158 11.063720826850961 3.6825070656623606 0.6902785321094302 4.96008586865757 0.6902785321094302 3.6825070656623606 11.063720826850961 5.881503709461158; 2.830830026003773 2.830830026003773 0.2240431494891658 5.881503709461155 4.671643508435733 0.0810140527710051 0.6181197482998197 0.0810140527710051 4.671643508435733 5.881503709461158 0.2240431494891658]

```
"""
psf(field::AbstractArray) = abs.(toimageplane(field)) .^2

psf(amplitude, phase) = psf(field(amplitude, phase))

subpsf(field::AbstractArray,Q::Integer) = psf(subdivide_sum(field,Q))
subpsf(amplitude, phase,Q::Integer) = psf(subdivide_sum(field(amplitude, phase),Q))

toimageplane(field) = ifftshift(fft(fftshift(field)))


text:
_basename(v::AbstractArray) = basename.(v)
_basename(::Nothing) = nothing
_basename(v) = basename(v)

"""
`filepicker(label="Choose a file..."; multiple=false, accept="*")`

Create a widget to select files.
If `multiple=true` the observable will hold an array containing the paths of all
selected files. Use `accept` to only accept some formats, e.g. `accept=".csv"`
"""
function filepicker(::WidgetTheme, lbl="Choose a file..."; attributes=PropDict(),
    label=lbl, className="", multiple=false, value=multiple ? String[] : "",  kwargs...)

    (value isa AbstractObservable) || (value = Observable{Any}(value))
    filename = Observable{Any}(_basename(value[]))

    if multiple
        onFileUpload = js"""function (data, e){
            var files = e.target.files;
            var fileArray = Array.from(files);
            this.filename(fileArray.map(function (el) {return el.name;}));
            return this.path(fileArray.map(function (el) {return el.path;}));
        }
        """
    else
        onFileUpload = js"""function(data, e) {
            var files = e.target.files;
            this.filename(files[0].name);
            return this.path(files[0].path);
        }
        """
    end
    multiple && (attributes=merge(attributes, PropDict(:multiple => true)))
    attributes = merge(attributes, PropDict(:type => "file", :style => "display: none;",
        Symbol("data-bind") => "event: {change: onFileUpload}"))
    className = mergeclasses(getclass(:input, "file"), className)
    template = dom"div[style=display:flex; align-items:center;]"(
        node(:label, className=getclass(:input, "file", "label"))(
            node(:input; className=className, attributes=attributes, kwargs...),
            node(:span,
                node(:span, (node(:i, className = getclass(:input, "file", "icon"))), className=getclass(:input, "file", "span", "icon")),
                node(:span, label, className=getclass(:input, "file", "span", "label")),
                className=getclass(:input, "file", "span"))
        ),
        node(:span, attributes = Dict("data-bind" => " text: filename() == '' ? 'No file chosen' : filename()"),
            className = getclass(:input, "file", "name"))
    )

    observs = ["path" => value, "filename" => filename]
    ui = knockout(template, observs, methods = ["onFileUpload" => onFileUpload])
    slap_design!(ui)
    Widget{:filepicker}(observs, scope = ui, output = ui["path"], layout = dom"div.field"∘Widgets.scope)
end

_parse(::Type{S}, x) where{S} = parse(S, x)
function _parse(::Type{Dates.Time}, x)
    h, m = parse.(Int, split(x, ':'))
    Dates.Time(h, m)
end

"""
`datepicker(value::Union{Dates.Date, Observable, Nothing}=nothing)`

Create a widget to select dates.
"""
function datepicker end

"""
`timepicker(value::Union{Dates.Time, Observable, Nothing}=nothing)`

Create a widget to select times.
"""
function timepicker end

for (func, typ, str, unit) in [(:timepicker, :(Dates.Time), "time", Dates.Second), (:datepicker, :(Dates.Date), "date", Dates.Day) ]
    @eval begin
        function $func(::WidgetTheme, val=nothing; value=val, kwargs...)
            (value isa AbstractObservable) || (value = Observable{Union{$typ, Nothing}}(value))
            f = x -> x === nothing ? "" : split(string(x), '.')[1]
            g = t -> _parse($typ, t)
            pair = ObservablePair(value, f=f, g=g)
            ui = input(pair.second; typ=$str, kwargs...)
            Widget{$(Expr(:quote, func))}(ui, output = value)
        end

        function $func(T::WidgetTheme, vals::AbstractRange, val=medianelement(vals); value=val, kwargs...)
            f = x -> x === nothing ? "" : split(string(x), '.')[1]
            fs = x -> x === nothing ? "" : split(string(convert($unit, x)), ' ')[1]
            $func(T; value=value, min=f(minimum(vals)), max=f(maximum(vals)), step=fs(step(vals)), kwargs...)
        end
    end
end

"""
`colorpicker(value::Union{Color, Observable}=colorant"#000000")`

Create a widget to select colors.
"""
function colorpicker(::WidgetTheme, val=colorant"#000000"; value=val, kwargs...)
    (value isa AbstractObservable) || (value = Observable{Color}(value))
    f = t -> "#"*hex(t)
    g = t -> parse(Colorant,t)
    pair = ObservablePair(value, f=f, g=g)
    ui = input(pair.second; typ="color", kwargs...)
    Widget{:colorpicker}(ui, output = value)
end

"""
`spinbox([range,] label=""; value=nothing)`

Create a widget to select numbers with placeholder `label`. An optional `range` first argument
specifies maximum and minimum value accepted as well as the step.
"""
function spinbox(::WidgetTheme, label=""; value=nothing, placeholder=label, isinteger=nothing, kwargs...)
    isinteger = something(isinteger, isa(_val(value), Integer))
    T = isinteger ? Int : Float64
    (value isa AbstractObservable) || (value = Observable{Union{T, Nothing}}(value))
    ui = input(value; isnumeric=true, placeholder=placeholder, typ="number", kwargs...)
    Widget{:spinbox}(ui, output = value)
end

spinbox(T::WidgetTheme, vals::AbstractRange, args...; value=first(vals), isinteger=(eltype(vals) <: Integer), kwargs...) =
    spinbox(T, args...; value=value, isinteger=isinteger, min=minimum(vals), max=maximum(vals), step=step(vals), kwargs...)

"""
`autocomplete(options, label=""; value="")`

Create a textbox input with autocomplete options specified by `options`, with `value`
as initial value and `label` as label.
"""
function autocomplete(::WidgetTheme, options, args...; attributes=PropDict(), kwargs...)
    (options isa AbstractObservable) || (options = Observable{Any}(options))
    option_array = _js_array(options)
    s = gensym()
    attributes = merge(attributes, PropDict(:list => s))
    t = textbox(args...; extra_obs=["options_js" => option_array], attributes=attributes, kwargs...)
    Widgets.scope(t).dom = node(:div,
        Widgets.scope(t).dom,
        node(:datalist, node(:option, attributes=Dict("data-bind"=>"value : key"));
            attributes = Dict("data-bind" => "foreach : options_js", "id" => s))
    )
    w = Widget{:autocomplete}(t)
    w[:options] = options
    w
end

"""
`input(o; typ="text")`

Create an HTML5 input element of type `type` (e.g. "text", "color", "number", "date") with `o`
as initial value.
"""
function input(::WidgetTheme, o; extra_js=js"", extra_obs=[], label=nothing, typ="text", wdgtyp=typ,
    className="", style=Dict(), isnumeric=Knockout.isnumeric(o),
    computed=[], attributes=Dict(), bind="value", bindto="value", valueUpdate="input", changes=0, kwargs...)

    (o isa AbstractObservable) || (o = Observable(o))
    (changes isa AbstractObservable) || (changes = Observable(changes))
    isnumeric && (bind == "value") && (bind = "numericValue")
    data = Pair{String, Any}["changes" => changes, bindto => o]
    append!(data, (string(key) => val for (key, val) in extra_obs))
    attrDict = merge(
        attributes,
        Dict(:type => typ, Symbol("data-bind") => "$bind: $bindto, valueUpdate: '$valueUpdate', event: {change : function () {this.changes(this.changes()+1)}}")
    )
    className = mergeclasses(getclass(:input, wdgtyp), className)
    template = node(:input; className=className, attributes=attrDict, style=style, kwargs...)()
    ui = knockout(template, data, extra_js; computed=computed)
    (label != nothing) && (ui.dom = flex_row(wdglabel(label), ui.dom))
    slap_design!(ui)
    Widget{:input}(data, scope = ui, output = ui[bindto], layout = dom"div.field"∘Widgets.scope)
end

function input(::WidgetTheme; typ="text", kwargs...)
    if typ in ["checkbox", "radio"]
        o = false
    elseif typ in ["number", "range"]
        o = 0.0
    else
        o = ""
    end
    input(o; typ=typ, kwargs...)
end

function input(T::WidgetTheme, ::Type{S}, args...; isinteger=nothing, kwargs...) where {S<:Number}
    (isinteger === nothing) && (isinteger = S<:Integer ? true : S<:AbstractFloat ? false : nothing)
    spinbox(T, args...; isinteger=isinteger, kwargs...)
end

input(T::WidgetTheme, ::Type{<:Bool}, args...; kwargs...) = toggle(T, args...; kwargs...)
input(T::WidgetTheme, ::Type{<:AbstractString}, args...; kwargs...) = textbox(T, args...; kwargs...)
input(T::WidgetTheme, ::Type{<:Dates.Date}, args...; kwargs...) = datepicker(T, args...; kwargs...)
input(T::WidgetTheme, ::Type{<:Dates.Time}, args...; kwargs...) = timepicker(T, args...; kwargs...)
input(T::WidgetTheme, ::Type{<:Color}, args...; kwargs...) = colorpicker(T, args...; kwargs...)

"""
`button(content... = "Press me!"; value=0)`

A button. `content` goes inside the button.
Note the button `content` supports a special `clicks` variable, that gets incremented by `1`
with each click e.g.: `button("clicked {{clicks}} times")`.
The `clicks` variable is initialized at `value=0`
"""
function button(::WidgetTheme, content...; label = "Press me!", value = 0, style = Dict{String, Any}(),
    className = getclass(:button, "primary"), attributes=Dict(), kwargs...)

    isempty(content) && (content = (label,))
    (value isa AbstractObservable) || (value = Observable(value))
    className = "delete" in split(className, ' ') ? className : mergeclasses(getclass(:button), className)
    attrdict = merge(
        Dict("data-bind"=>"click : function () {this.clicks(this.clicks()+1)}"),
        attributes
    )
    template = node(:button, content...; className=className, attributes=attrdict, style=style, kwargs...)
    button = knockout(template, ["clicks" => value])
    slap_design!(button)
    Widget{:button}(scope = button, output = value, layout = dom"div.field"∘Widgets.scope)
end

for wdg in [:toggle, :checkbox]
    @eval begin
        $wdg(::WidgetTheme, value, lbl::AbstractString=""; label=lbl, kwargs...) =
            $wdg(gettheme(); value=value, label=label, kwargs...)

        $wdg(::WidgetTheme, label::AbstractString, val=false; value=val, kwargs...) =
            $wdg(gettheme(); value=value, label=label, kwargs...)

        $wdg(::WidgetTheme, value::AbstractString, label::AbstractString; kwargs...) =
            error("value cannot be a string")

        function $wdg(::WidgetTheme; bind="checked", valueUpdate="change", value=false, label="", labelclass="", kwargs...)
            s = gensym() |> string
            (label isa Tuple) || (label = (label,))
            widgettype = $(Expr(:quote, wdg))
            wdgtyp = string(widgettype)
            labelclass = mergeclasses(getclass(:input, wdgtyp, "label"), labelclass)
            ui = input(value; bind=bind, typ="checkbox", valueUpdate="change", wdgtyp=wdgtyp, id=s, kwargs...)
            Widgets.scope(ui).dom = dom"div.field"(Widgets.scope(ui).dom, dom"label[className=$labelclass, for=$s]"(label...))
            Widget{widgettype}(ui)
        end
    end
end

"""
`checkbox(value::Union{Bool, AbstractObservable}=false; label)`

A checkbox.
e.g. `checkbox(label="be my friend?")`
"""
function checkbox end

"""
`toggle(value::Union{Bool, AbstractObservable}=false; label)`

A toggle switch.
e.g. `toggle(label="be my friend?")`
"""
function toggle end

"""
`textbox(hint=""; value="")`

Create a text input area with an optional placeholder `hint`
e.g. `textbox("enter number:")`. Use `typ=...` to specify the type of text. For example
`typ="email"` or `typ=password`. Use `multiline=true` to display a `textarea` spanning
several lines.
"""
function textbox(::WidgetTheme, hint=""; multiline=false, placeholder=hint, value="", typ="text", kwargs...)
    multiline && return textarea(gettheme(); placeholder=placeholder, value=value, kwargs...)
    Widget{:textbox}(input(value; typ=typ, placeholder=placeholder, kwargs...))
end

"""
`textarea(hint=""; value="")`

Create a textarea with an optional placeholder `hint`
e.g. `textarea("enter number:")`. Use `rows=...` to specify how many rows to display
"""
function textarea(::WidgetTheme, hint=""; label=nothing, className="",
    placeholder=hint, value="", attributes=Dict(), style=Dict(), bind="value", valueUpdate = "input", kwargs...)

    (value isa AbstractObservable) || (value = Observable(value))
    attrdict = convert(PropDict, attributes)
    attrdict[:placeholder] = placeholder
    attrdict["data-bind"] = "$bind: value, valueUpdate: '$valueUpdate'"
    className = mergeclasses(getclass(:textarea), className)
    template = node(:textarea; className=className, attributes=attrdict, style=style, kwargs...)
    ui = knockout(template, ["value" => value])
    (label != nothing) && (ui.dom = flex_row(wdglabel(label), ui.dom))
    slap_design!(ui)
    Widget{:textarea}(scope = ui, output = ui["value"], layout = dom"div.field"∘Widgets.scope)
end

function wdglabel(T::WidgetTheme, text; padt=5, padr=10, padb=0, padl=10,
    className="", style = Dict(), kwargs...)

    className = mergeclasses(getclass(:wdglabel), className)
    padding = Dict(:padding=>"$(padt)px $(padr)px $(padb)px $(padl)px")
    node(:label, text; className=className, style = merge(padding, style), kwargs...)
end

function flex_row(a,b,c=dom"div"())
    dom"div.[class=interact-flex-row]"(
        dom"div[class=interact-flex-row-left]"(a),
        dom"div[class=interact-flex-row-center]"(b),
        dom"div[class=interact-flex-row-right]"(c)
    )
end

flex_row(a) = dom"div.[class=interact-flex-row]"(a)


text:
using DimReduc
using LazySets
import DimReduc: forward_network, forward_affine_map

nnet1 = read_nnet("nnet/mnist_1000.nnet")
nnet2 = read_nnet("nnet/mnist_large.nnet")
nnet3 = read_nnet("nnet/mnist3.nnet")
nnet4 = read_nnet("nnet/mnist4.nnet")

solver = FastGrid(1)

center = zeros(784)
radius = 0.004
inputSet = Hyperrectangle(low=center .- radius, high=center .+ radius)

file = open("results/affinemap.txt", "a")

print(file, "mnist1\n")
W, b = nnet1.layers[1].weights, nnet1.layers[1].bias
Q = forward_affine_map(solver, W, b, inputSet)
print(file, "width = $(sort(Q.radius) .* 2)\n")

print(file, "mnist2\n")
W, b = nnet2.layers[1].weights, nnet2.layers[1].bias
Q = forward_affine_map(solver, W, b, inputSet)
print(file, "width = $(sort(Q.radius) .* 2)\n")

print(file, "mnist3\n")
W, b = nnet3.layers[1].weights, nnet3.layers[1].bias
Q = forward_affine_map(solver, W, b, inputSet)
print(file, "width = $(sort(Q.radius) .* 2)\n")

print(file, "mnist4\n")
W, b = nnet4.layers[1].weights, nnet4.layers[1].bias
Q = forward_affine_map(solver, W, b, inputSet)
print(file, "width = $(sort(Q.radius) .* 2)\n")

close(file)


text:
function single(xi, qi, ti, mu, eta, kappa)

    ei = ones(length(ti))
    fi = zeros(size(xi))

    return ei, fi
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
#!/usr/bin/env julia
println("--> Julia $VERSION")

try
  using BenchmarkTools
catch
  import Pkg
  Pkg.add("BenchmarkTools")
  Pkg.update("BenchmarkTools")
  using BenchmarkTools
end

if length(ARGS) < 1
  N = 1000
else
  N = parse(Int, ARGS[1]);
end

A = randn(N,N);
B = randn(N,N);

o=@benchmark A*B;
println("N: $N")
println(minimum(o).time/1e9, " seconds")

max_stars_repo_name: scienceopen/numba-examples

max_issues_repo_name: scienceopen/numba-examples

max_forks_repo_name: scienceopen/numba-examples


text:
module ScatteringTransform
using Distributed, SharedArrays
using LinearAlgebra, Interpolations, Wavelets, FFTW
using SpecialFunctions, LinearAlgebra
using HDF5, JLD
using Plots
# todo remove nScales controls

# TODO: integrate path methods
struct pathType{T}
  m::Int64
  Idxs::Array{T,1}
  function pathType(m::Int64, Idxs)
    @assert m == length(Idxs)
    @assert minimum([true, [tp <: Integer for tp in typeof.(Idxs[1:end-1])]...])
    new{eltype(Idxs)}(m, Idxs)
  end
end

"""
    path = pathType(m::Int64, indexInShear::Array{Int64,1}, layers::stParallel)

path constructor that uses the fixed index in each layer (i.e. the index of the shearlet filter in shearletIdxs)
"""
function pathType(indexInLayers::Array)
  pathType(length(indexInLayers), indexInLayers)
end

import Base: show
function Base.show(io::IO, p::pathType{T}) where {T}
  pl = p.Idxs[end]
  if typeof(pl) <: Colon
    print(io, "pathType[$(p.Idxs[1:end-1]...), :]")
  else
    print(io, "pathType$(p.Idxs)")
  end
end

export pathType
include("subsampling.jl")
export resample, sizes, bsplineType, bilinearType, autocorrType
include("modifiedTransforms.jl")
export WT, wavelet, cwt, getScales, computeWavelets
include("basicTypes.jl")
export stParallel, Scattered
include("Utils.jl")
export calculateThinStSizes, getPadBy, pad, outputSize, createFFTPlans, remoteMultiply, createRemoteFFTPlan, computeAllWavelets, plotAllWavelets
include("nonlinearities.jl")
export nonlinearity, absType, ReLUType, tanhType, softplusType, piecewiseType, spInverse, aTanh, Tanh, ReLU, piecewiseLinear, plInverse
include("transform.jl")
export st, transformMidLayer!, transformFinalLayer!
include("pathMethods.jl")
export pathToThinIndex
export flatten, MatrixAggrigator, plotCoordinate, reshapeFlattened,
  numberSkipped, logabs, maxPooling, numScales, incrementKeeper,
  numInLayer
include("postProcessing.jl")
export logabs, ReLU, MatrixAggrigator, reshapeFlattened,
  loadSyntheticMatFile, transformFolder, flatten

# TODO make a way to access a Scattered2D by the index of (depth, scale,shearingFactor)
end