text:
#################################
###### Polynomial programs ######
#################################

#  A system of polynomial equations.
#  We convert inequalities into equations with slack variables (squared).
type Program
    slackvars :: Int  # counter
    vars :: Set{Symbol}  # all variables, including slack variables
    constraints :: Array{SoSPoly}  # polynomials that are constrained to 0
    namedconstraints :: Dict{Any,SoSPoly}
    objective :: SoSPoly
    maximize :: Bool
    symmetries :: Array{Dict{Symbol,Symbol}}
end
Program(; maximize=true, minimize=false) = begin
    mx = (maximize && !minimize) # default to maximize; switch to minimize if there's any sign to do so
    Program(0, Set{Symbol}(), [], Dict{Any,SoSPoly}(), SoSPoly(), mx, [])
end

#  Generate a new slack variable (a Symbol)
slackvar(prog :: Program) = symbol( @sprintf("slack%d", prog.slackvars += 1) )

# convenience
monoms(prog :: Program, d :: Int64) = monoms([prog.vars...], d)

typealias Maybe{T} Union{T,Void}

#  Parse a polynomial (in Expr form) into a SoSPoly, by traversing
#    a Julia syntax tree.

parsepoly{T<:Number}(prog :: Maybe{Program}, ex :: T) = Dict( SoSMonom() => convert(Float64, ex) ) :: SoSPoly

function parsepoly(prog :: Maybe{Program}, ex::Symbol)
    prog == nothing || push!(prog.vars, ex)
    Dict(Dict(ex => 1) => 1.0) :: SoSPoly
end

function parsepoly(prog :: Maybe{Program}, ex :: Expr)
    if ex.head == :call
        if(ex.args[1] == :(+))
            ret = parsepoly(prog, ex.args[2])
            for i in 3:length(ex.args)
                ret += parsepoly(prog, ex.args[i])
            end
            return ret
        elseif(ex.args[1] == :(*))
            ret = parsepoly(prog, ex.args[2])
            for i in 3:length(ex.args)
                ret *= parsepoly(prog, ex.args[i])
            end
            return ret
        elseif(ex.args[1] == :(-))
            if(length(ex.args) == 2) # unary minus
                return parsepoly(prog, :( (-1) * $(ex.args[2]) ))
            elseif(length(ex.args) == 3) # binary minus
                return parsepoly(prog, :( $(ex.args[2]) + (-1) * $(ex.args[3]) ))
            end
        elseif(ex.args[1] == :(^))
            if(! isa(ex.args[3],Int64) || ex.args[3] < 0)
                throw(ArgumentError("Bad exponent: $(args[3]). Exponents must be non-negative integers"))
            elseif(ex.args[3] == 0)
                return parsepoly(prog, 1)
            end
            operand = parsepoly(prog, ex.args[2])
            return prod([operand for i in 1:ex.args[3]])
        end
    elseif ex.head == :block # this happens on the RHS of single-equals
        if length(ex.args) > 2
            throw(ArgumentError(@sprintf("Multi-line block in expression: %s", ex)))
        end
        return parsepoly(prog, ex.args[2]) # assume it's a 1-line block
    end
    throw(ArgumentError("Expression not understood as polynomial: $ex"))
end

#  Massage any constraint (given in Expr form) into a polynomial equaling zero.
function relationToPolynomial(prog :: Program, ex :: Expr)
    if ex.head == :comparison
        if(ex.args[2] == :(>=) || ex.args[2] == :(≥))
            polyex = :( $(ex.args[1]) - $(ex.args[3]) - $(slackvar(prog))^2 )
        elseif(ex.args[2] == :(==))
            polyex = :( $(ex.args[1]) - $(ex.args[3]) )
        elseif(ex.args[2] == :(<=) || ex.args[2] == :(≤))
            polyex = :( $(ex.args[3]) - $(ex.args[1]) - $(slackvar(prog))^2 )
        else
            throw(ArgumentError(
                "Comparison operator must be >=, ≥, ==, =, <=, or ≤"))
        end
    elseif ex.head == :(=) # isn't Julia comparison, but we should support
        polyex = :( $(ex.args[1]) - $(ex.args[2]) )
    else
        throw(ArgumentError("Constraint must be an equation or inequality"))
    end

    parsepoly(prog, polyex)
end

#  Perturb the objective by Gaussian noise with standard deviation `tol`.
#  Each moment up to degree `deg`, in the variables known to `prog`, receives
#  noise.
function perturb_objective(prog,tol,deg)
    for md in 1:deg
        for mon in monoms(prog,md)
            poly = Dict( mon => tol * randn() ) :: SoSPoly
            prog.objective += poly
        end
    end
end

#  Perturb the objective by Gaussian noise with standard deviation `tol`.
#  A random choice of `s` terms (with replacement), from the moments up to
#  degree `deg`, in the variables known to `prog`, receives noise.
function perturb_objective_sparse(prog,tol,s,deg)
    mon = [ monoms(prog,i) for i in 0:deg ]
    mon = vcat(mon...)
    for i in 1:s
        m = rand(mon)
        prog.objective += Dict( m => tol*randn() )
    end
end

#  Add a constraint to the program.
function constraint(prog,str)
    poly = relationToPolynomial(prog, parse(str))
    push!(prog.constraints, poly)
end

#  Add a polynomial to a constraint indexed by `key`.
function partconstraint(prog,key,str)
    poly = parsepoly(prog, parse(str))

    if haskey(prog.namedconstraints, key)
        prev = prog.namedconstraints[key]
        addpoly!(prev,poly)
    else
        push!(prog.constraints, poly)
        prog.namedconstraints[key] = poly
    end
end

#  Set a polynomial as the objective.
objective(prog,str) = prog.objective = parsepoly(prog, parse(str))

#  Add a polynomial to the objective.
partobjective(prog,str) = prog.objective += parsepoly(prog, parse(str))

#  `printf`-style macros for constraints and objectives.
#  This section is `esc` hell and was a pain to figure out.
#  Since we're calling a macro from a macro, we have three different
#  scopes in play at once, and we need to make sure everything gets evaluated
#  in the right scope.

#  Add a constraint (given in printf form) to the program.
macro constraint(prog,args...)
    escargs = [esc(x) for x in args[2:end]]
    quote
        str = @sprintf($(args[1]),$(escargs...))
        constraint( $(esc(prog)), str )
    end
end

#  Add a polynomial (given in printf form) to a constraint indexed by `key`.
macro partconstraint(prog,key,args...)
    escargs = [esc(x) for x in args[2:end]]
    quote
        str = @sprintf($(args[1]),$(escargs...))
        partconstraint( $(esc(prog)), $(esc(key)), str )
    end
end

#  Set a polynomial (given in printf form) as the objective.
macro objective(prog,args...)
    escargs = [esc(x) for x in args[2:end]]
    quote
        str = @sprintf($(args[1]),$(escargs...))
        objective( $(esc(prog)), str )
    end
end

#  Add a polynomial (given in printf form) to the objective.
macro partobjective(prog,args...)
    escargs = [esc(x) for x in args[2:end]]
    quote
        str = @sprintf($(args[1]),$(escargs...))
        partobjective( $(esc(prog)), str )
    end
end

max_stars_repo_name: willperry/SumOfSquaresOptimization


text:
using QuizQuestions
using Test

@testset "QuizQuestions.jl" begin
    @test numericq(1).atol == 0
    @test numericq(1, 1e-3; hint="hint").hint == "hint"

    @test stringq(r"abc").re == r"abc"
    @test_throws MethodError stringq("abc") # need regular expression

    @test radioq(["one","two", "three"], 1, keep_order=true).answer == 1
    @test radioq(1:3, 1, keep_order=true).answer == 1

    @test buttonq(["one","two", "three"], 1).answer == ["correct", "incorrect", "incorrect"]

    @test multiq(("one","two","three"), (2, 3), label="Primes?").label == "Primes?"

    @test fillblankq("question ____", r"1", label="one").label == "one"
    @test fillblankq("question ____", ("1","2","3"), 1, label="one").label == "one"
    @test fillblankq("question ____", 1, label="one").label == "one"

    d = Dict("Select a Volvo" => "XC90", "Select a Mercedes" => "GLE 350",
             "Select an Audi" => "A4")
    r = matchq(d)
    i = rand(1:3)
    @test d[r.questions[i]] == r.choices[r.answer[i]]

    r = hotspotq("empty", (0,0), (1,1); explanation="XXX")
    @test r.explanation == "XXX"

    r = plotlylightq("empty"; label="XXX")
    @test r.label == "XXX"

end

ext: jl


text:
# Cost functions
export quaderror, ceerror

#============================#
#      QUADRATIC ERROR       #
#============================#

"Quadratic error between the output of the network and a target"
quaderror(output, target) = 0.5*norm(target - output)^2

"Derivative of the quadratic error with respect to the outputs"
quaderrorprime(output, target) = (output - target)

#============================#
#    CROSS-ENTROPY ERROR     #
#============================#

"Cross-entropy error between the output of the network and a target"
ceerror(output, target) = -output'log(target)

"Derivative of the cross-entropy error with respect to the outputs"
ceerrorprime(output, target) = -target ./ output


text:
using MicroCoverage
using Test

module MicroCoverageTestSuite

using MicroCoverage
using Test

function get_test_directory()::String
    return dirname(@__FILE__)
end

function get_filename(parts...)::String
    return joinpath(get_test_directory(), parts...)
end

function readstring_filename(parts...)::String
    _filename = get_filename(parts...)
    return read(_filename, String)
end

@testset "MicroCoverage.jl" begin
    @testset "assert.jl" begin
        @test MicroCoverage.always_assert(true, "") == nothing
        @test_throws MicroCoverage.AlwaysAssertionError MicroCoverage.always_assert(false, "")
    end
    @testset "instrument.jl" begin
        MicroCoverage.instrument("", Expr(:block), Val(:block)) == Expr(:block)
    end
    @testset "public_interface.jl" begin
        MicroCoverage.with_temp_dir() do tmp_depot
            MicroCoverage.with_temp_dir() do tmp_src_directory
                nonexistent_filename = joinpath(tmp_src_directory, "nonexistent.jl")
                foo_jl_filename = joinpath(tmp_src_directory, "foo.jl")
                test_foo_jl_filename = joinpath(tmp_src_directory, "test_foo.jl")
                foo_jl_microcov_filename = joinpath(tmp_src_directory, "foo.jl.microcov")
                bar_jl_filename = joinpath(tmp_src_directory, "bar.jl")
                test_bar_jl_filename = joinpath(tmp_src_directory, "test_bar.jl")
                bar_jl_microcov_filename = joinpath(tmp_src_directory, "bar.jl.microcov")

                open(foo_jl_filename, "w") do io
                    print(io, readstring_filename("inputs", "foo.jl"))
                end
                open(test_foo_jl_filename, "w") do io
                    print(io, readstring_filename("inputs", "test_foo.jl"))
                end
                open(bar_jl_filename, "w") do io
                    print(io, readstring_filename("inputs", "bar.jl"))
                end
                open(test_bar_jl_filename, "w") do io
                    print(io, readstring_filename("inputs", "test_bar.jl"))
                end

                @test_throws ArgumentError MicroCoverage.start(nonexistent_filename)
                @test_throws ArgumentError MicroCoverage.clean(nonexistent_filename)
                MicroCoverage.start(strip(foo_jl_filename))
                MicroCoverage.start(bar_jl_filename)

                MicroCoverage.preview_coverage(; dump_coverage_io = devnull)
                include(foo_jl_filename)
                MicroCoverage.preview_coverage(; dump_coverage_io = devnull)
                include(bar_jl_filename)
                MicroCoverage.preview_coverage(; dump_coverage_io = devnull)
                include(test_foo_jl_filename)
                MicroCoverage.preview_coverage(; dump_coverage_io = devnull)
                include(test_bar_jl_filename)
                MicroCoverage.preview_coverage(; dump_coverage_io = devnull)

                @test !isfile(foo_jl_microcov_filename)
                @test !ispath(foo_jl_microcov_filename)
                @test !isfile(bar_jl_microcov_filename)
                @test !ispath(bar_jl_microcov_filename)

                MicroCoverage.stop(; dump_coverage = true, dump_coverage_io = devnull)

                @test isfile(foo_jl_microcov_filename)
                @test ispath(foo_jl_microcov_filename)
                @test isfile(bar_jl_microcov_filename)
                @test ispath(bar_jl_microcov_filename)

                if Sys.iswindows()
                    @test chomp(read(foo_jl_microcov_filename, String)) == chomp(readstring_filename("expected_outputs", "foo.jl.microcov.expected"))
                    @test chomp(read(bar_jl_microcov_filename, String)) == chomp(readstring_filename("expected_outputs", "bar.jl.microcov.expected"))
                else
                    @test read(foo_jl_microcov_filename, String) == readstring_filename("expected_outputs", "foo.jl.microcov.expected")
                    @test read(bar_jl_microcov_filename, String) == readstring_filename("expected_outputs", "bar.jl.microcov.expected")
                end
                touch(foo_jl_microcov_filename)
                MicroCoverage.clean(foo_jl_filename)
                touch(foo_jl_microcov_filename)
                MicroCoverage.clean(tmp_src_directory)
            end
        end
    end
end

end # end module MicroCoverageTestSuite


text:
using SharedArrays
using Distributed
export Fitness, GRN, evolve

struct Fitness
    nin::Int64
    nout::Int64
    func::Function
    cacheable::Bool
end

function Fitness(nin::Int64, nout::Int64, func::Function)
    Fitness(nin, nout, func, true)
end

function tournament(fits::Array{Float64}; tsize=3)
    # return the index of the winner of a n-way tournament
    if length(fits) == 1
        return 1
    else
        n = min(tsize, length(fits))
        fshuffle = randperm(length(fits))[1:n]
        winner = indmax(fits[fshuffle])
        return fshuffle[winner]
    end
end

function speciation(population::Array, reprs::Array, config::Config;
                    dfunc::Function=AGRN.distance)
    # return a vector of ints corresponding to the species of each individual in population
    ptype = typeof(population[1])
    species = Array{Int64}(length(population))
    for p in eachindex(population)
        distances = Array{Float64}(length(reprs))
        for r in eachindex(reprs)
            distances[r] = dfunc(population[p], reprs[r], config)
        end
        # distances = [distance(population[p], reprs[r]) for r in eachindex(reprs)]
        if minimum(distances) < config.speciation_thresh
            species[p] = indmin(distances)
        else
            species[p] = length(reprs)+1
            append!(reprs, [ptype(population[p])])
        end
    end
    species_set = sort(unique(species))
    for s in eachindex(species_set)
        species[species.==species_set[s]] = s
    end
    species
end

function species_sizes(fits::Array{Float64}, species::Array{Int64}, config::Config)
    nspecies = length(unique(species))
    spec_fits = map(x->mean(fits[species.==x])-minimum(fits), 1:nspecies)
    @debug("spec fits: $spec_fits")
    spec_sizes = map(x->spec_fits[x]/sum(spec_fits), 1:nspecies)
    spec_sizes = spec_sizes./sum(spec_sizes).*config.ga_population
    spec_sizes[isnan.(spec_sizes)] = 0
    spec_sizes = Int64.(round.(spec_sizes))

    while sum(spec_sizes) > config.ga_population
        spec_sizes[indmax(spec_sizes)] -= 1
    end
    while sum(spec_sizes) < config.ga_population
        spec_sizes[indmin(spec_sizes)] += 1
    end
    @debug("spec sizes: $spec_sizes")
    spec_sizes
end

function evolve(fitness::Fitness, config::Config)
    population = Array{GRN}(config.ga_population)
    fits = Array{Float64}(config.ga_population)
    fits .= -Inf
    species = Array{Int64}(config.ga_population)
    for p in eachindex(population)
        if p <= config.init_species
            population[p] = GRN(fitness.nin, fitness.nout, 1, config)
            species[p] = p
        else
            repr = rand(1:config.init_species)
            population[p] = mutate(population[repr], config)
            species[p] = repr
        end
    end
    best = population[1]
    max_fit = -Inf
    eval_count = 0

    for generation in 1:config.ga_num_generations
        # evaluation
        @debug("evaluation $generation")
        dfits = SharedArray{Float64}(config.ga_population)
        dfits .= fits
        @sync @distributed for p=1:config.ga_population
            dfits[p] = fitness.func(population[p])
        end
        fits .= dfits

        new_best = false
        for p in eachindex(population)
            fit = fits[p]
            eval_count += 1
            if fit > max_fit
                max_fit = fit
                best = population[p]
                new_best = true
            end
        end

        @info(@sprintf("E: %d %d %0.5f %0.5f",
                              generation, eval_count, mean(fits), max_fit))
        if new_best
            @info(@sprintf("B: %d %0.5f|%s", generation, max_fit,
                                  JSON.json(best)))
        end

        # representatives
        @debug("representatives $generation")
        nspecies = length(unique(species))
        reprs = Array{GRN}(nspecies)
        for s in 1:nspecies
            reprs[s] = GRN(population[rand(find(species.==s))])
        end

        # species sizes
        @debug("species sizes $generation $nspecies")
        spec_sizes = species_sizes(fits, species, config)
        new_pop = Array{GRN}(config.ga_population)
        new_fits = -Inf*ones(Float64, config.ga_population)
        popi = 1

        # create new population
        @debug("new population $generation $spec_sizes")
        for s in 1:nspecies
            @debug("species $s")
            sfits = fits[species.==s]
            spec = population[species.==s]
            ncross = round(spec_sizes[s] * config.ga_crossover_rate)
            nmut = round(spec_sizes[s] * config.ga_mutation_rate)
            ncopy = spec_sizes[s] - (ncross+nmut)
            while ncopy < 0
                if nmut > ncross
                    nmut -= 1
                else
                    ncross -= 1
                end
                ncopy = spec_sizes[s] - (ncross+nmut)
            end

            if popi > length(new_pop)
                break
            end

            # crossover
            @debug("Crossover $s popi: $popi, ncross: $ncross")
            for i in 1:ncross
                p1 = spec[tournament(sfits)]
                p2 = spec[tournament(sfits)]
                new_pop[popi] = crossover(p1, p2, config)
                popi += 1
            end

            if popi > length(new_pop)
                break
            end

            # mutation
            @debug("Mutation $s popi: $popi, nmut: $nmut")
            for i in 1:nmut
                parent = spec[tournament(sfits)]
                child = mutate(parent, config)
                if fitness.cacheable
                    if distance(parent, child, config) == 0
                        new_fits[popi] = fits[popi]
                    end
                end
                new_pop[popi] = child
                popi += 1
            end

            if popi > length(new_pop)
                break
            end

            # copy
            @debug("Copy $s popi: $popi, ncopy: $ncopy")
            for i in 1:ncopy
                new_pop[popi] = GRN(spec[tournament(sfits)])
                if fitness.cacheable
                    new_fits[popi] = fits[popi]
                end
                popi += 1
            end

            if popi > length(new_pop)
                break
            end
        end

        @debug("variable set $generation")
        species = speciation(new_pop, reprs, config)
        population = new_pop
        fits = new_fits

        @debug("done $generation")
    end

    max_fit, best
end


text:
module HiveLoader

# https://github.com/JuliaDatabases/Hive.jl v0.3.0
using Hive # HiveSession HiveAuth
using Octo.Repo: ExecuteResult

const current = Dict{Symbol, Any}(
    :sess => nothing,
)

current_sess() = current[:sess]

# db_connect
function db_connect(; host::String="localhost", port::Integer=10000, auth::HiveAuth=HiveAuth(), tprotocol::Symbol=:binary)
    sess = HiveSession(host, port, auth; tprotocol=tprotocol)
    current[:sess] = sess
end

# db_disconnect
function db_disconnect()
    sess = current_sess()
    if sess isa HiveSession
        Hive.close(sess)
        current[:sess] = nothing
    end
end

# query
function query(sql::String)
    sess = current_sess()
    pending = Hive.execute(sess, sql)
    rs = Hive.result(pending)
    sch = Hive.schema(rs)
    column_names = tuple(Symbol.(getproperty.(sch.columns, :columnName))...)
    df = reduce(vcat, Hive.records(rs))
    nts = NamedTuple{column_names}.(df)
    Hive.close(rs)
    nts
end

function query(prepared::String, vals::Vector) # throw UnsupportedError
    throw(UnsupportedError("needs to be implemented"))
end

# execute
function execute(sql::String)::ExecuteResult
    sess = current_sess()
    result = Hive.execute(sess, sql)
    ExecuteResult()
end

function execute(prepared::String, vals::Vector)::ExecuteResult # throw UnsupportedError
    throw(UnsupportedError("needs to be implemented"))
end

function execute(prepared::String, nts::Vector{<:NamedTuple})::ExecuteResult # throw UnsupportedError
    throw(UnsupportedError("needs to be implemented"))
end

end # module Octo.Backends.HiveLoader


text:
using Makie
using MakieLayout

begin
    scene = Scene(camera=campixel!)
    screen = display(scene)

    outergrid = GridLayout(scene; alignmode=Outside(30))

    innergrid = outergrid[1, 1] = GridLayout()

    innergrid[1, 1] = LText(scene, text="Hello", textsize=40)
    innergrid[1, 2] = LText(scene, text="Hello", textsize=40, width=80)
    innergrid[1, 3] = LText(scene, text="Hello", textsize=40, width=120)
    innergrid[1, 4] = LText(scene, text="Hello", textsize=40, width=160)
    innergrid[1, 5] = LText(scene, text="Hello", textsize=40, width=180)
    innergrid[1, :] = [LRect(scene) for _ in 1:5]
end


text:
using Mea.SigmaDC
using Mea.Periodize

using Base.Test

@testset "Testing SigmaDC.jl" begin

    modelvec=Periodize.buildmodelvec("./data/self_ctow0.dat", "./data/statsparams0.json")
    @test isapprox(SigmaDC.calc_sigmadc(modelvec, 12.0), 0.7637, atol=1e-3)

end  # testset


text:
@testset "DictChain" begin
    @testset "Matching keys and value types for get()/set() methods" begin
        dc1 = DictChain{Int,String}()

        @test_throws KeyError dc1[:NotThere]
        @test_throws MethodError (dc1[:NotThere] = 3)
        @test_throws MethodError (dc1[:NotThere] = "abc")
        @test_throws MethodError (dc1[3] = 3)

        dc1[3] = "abc"
        @test dc1[3] == "abc"
    end

    @testset "merging and chaining" begin
        dc1 = DictChain{Symbol,String}()

        # incompatible dictionary types
        @testset "incompatible dictionary types" begin
            dc2 = DictChain{Int,String}()
            @test_throws MethodError chain(dc1, dc2)
            @test_throws MethodError merge(dc1, dc2)
            @test_throws MethodError merge!(dc1, dc2)

            dc3 = DictChain{Symbol,Int}()
            @test_throws MethodError chain(dc1, dc3)
            @test_throws MethodError merge(dc1, dc3)
            @test_throws MethodError merge!(dc1, dc3)
        end

        d1 = Dict{Symbol,Int}(:a => 1)
        d2 = Dict{Symbol,Int}(:a => 2, :b => 4)
        d3 = Dict{Symbol,Int}(:b => 3, :c => 5)

        @testset "using constructor" begin
            dc = DictChain(d1, d2, d3)
            @test typeof(dc) == DictChain{Symbol,Int}
            @test (dc[:a], dc[:b], dc[:c]) == (1, 4, 5)

            dc = DictChain(d2, d1, d3)
            @test (dc[:a], dc[:b], dc[:c]) == (2, 4, 5)

            dc = DictChain(DictChain(d3, d1), d2)
            @test (dc[:a], dc[:b], dc[:c]) == (1, 3, 5)
        end

        @testset "using merge()" begin
            dc = merge(DictChain(d2, d1), d3)
            @test (dc[:a], dc[:b], dc[:c]) == (2, 3, 5)

            dc = merge(d2, DictChain(d3, d1))
            @test (dc[:a], dc[:b], dc[:c]) == (1, 3, 5)

            dc = merge(DictChain(d1, d3), d2)
            @test (dc[:a], dc[:b], dc[:c]) == (2, 4, 5)
        end

        @testset "using merge!()" begin
            dc = DictChain{Int,String}()
            @test_throws MethodError merge!(dc, d1) # incompatible

            dc = DictChain{Symbol,Int}()
            merge!(dc, d1)
            @test dc[:a] == 1
            @test_throws KeyError dc[:b]
            merge!(dc, d2)
            @test (dc[:a], dc[:b]) == (2, 4)
            merge!(dc, d3)
            @test (dc[:a], dc[:b], dc[:c]) == (2, 3, 5)
        end

        @testset "using chain()" begin
            dc = chain(chain(d1, d2), d3)
            @test typeof(dc) == DictChain{Symbol,Int}
            @test (dc[:a], dc[:b], dc[:c]) == (2, 3, 5)

            dc = chain(d2, chain(d1, d3))
            @test (dc[:a], dc[:b], dc[:c]) == (1, 3, 5)

            dc = chain(chain(d3, d1), d2)
            @test (dc[:a], dc[:b], dc[:c]) == (2, 4, 5)

            dc = chain(d1, d2, d3)
            @test (dc[:a], dc[:b], dc[:c]) == (2, 3, 5)

            dc = chain(d3, d1, d2)
            @test (dc[:a], dc[:b], dc[:c]) == (2, 4, 5)
        end

        @testset "Check chaining from real use case" begin
            # This sequence is used in compare_optimizers.jl:repeated_bboptimize
            parameters = BlackBoxOptim.EMPTY_PARAMS
            ftol = 1e-5
            params = chain(parameters, ParamsDict(:FitnessTolerance => ftol))
            @test params[:FitnessTolerance] == ftol
        end
    end

    @testset "converting to Dict" begin
        d1 = Dict{Symbol,Int}(:a => 1)
        d2 = Dict{Symbol,Int}(:a => 2, :b => 4)
        d3 = Dict{Symbol,Int}(:a => 3, :b => 5)

        dc = DictChain(d1, d2, d3)
        d123 = convert(Dict{Symbol,Int}, dc)
        @test typeof(d123) == Dict{Symbol,Int}
        @test length(d123) == 2
        @test d123[:a] == 1
        @test d123[:b] == 4
    end

    @testset "show()" begin
        d1 = Dict{Symbol,Int}(:a => 1)
        d2 = Dict{Symbol,Int}(:a => 2, :b => 4)
        d3 = Dict{Symbol,Int}(:a => 3, :b => 5)

        dc = DictChain(d1, d2, d3)
        iob = IOBuffer()
        show(iob, dc)
        @test replace(String(take!(iob)), ' '=>"") == "BlackBoxOptim.DictChain{Symbol,$Int}[Dict(:a=>1),Dict(:a=>2,:b=>4),Dict(:a=>3,:b=>5)]"
    end

    @testset "flatten" begin
        d1 = Dict{Symbol,Int}(:a => 1)
        d2 = Dict{Symbol,Int}(:a => 2, :b => 4)
        d3 = Dict{Symbol,Int}(:a => 3, :b => 5)

        dc = DictChain(d1, d2, d3)

        fd = flatten(dc)
        @test fd[:a] == 1
        @test fd[:b] == 4
        @test sort(collect(keys(fd))) == [:a, :b]
    end
end

@testset "Parameters" begin

    @testset "When no parameters or key type doesn't match" begin
        ps = ParamsDictChain()
        @test isa(ps, Parameters)
        @test_throws KeyError ps[:NotThere]
        @test_throws KeyError ps["Neither there"]

        ps[:a] = 1
        @test ps[:a] == 1
    end

    @testset "With one parameter in one set" begin
        ps = ParamsDictChain(ParamsDict(:a => 1))
        @test isa(ps, Parameters)

        @test ps[:a] == 1
        @test_throws KeyError ps["a"] # incorrect key

        @test_throws KeyError ps[:A]
        @test_throws KeyError ps[:B]
    end

    @testset "With parameters in multiple sets" begin
        ps = ParamsDictChain(ParamsDict(:a => 1, :c => 4),
                             ParamsDict(:a => 2, :b => 3),
                             ParamsDict(:c => 5))
        @test isa(ps, Parameters)

        @test ps[:a] == 1
        @test ps[:c] == 4
        @test ps[:b] == 3

        @test_throws KeyError ps[:A]
        @test_throws KeyError ps[:B]
    end

    @testset "Updating parameters after construction" begin
        ps = ParamsDictChain(ParamsDict(:a => 1, :c => 4),
                             ParamsDict(:a => 2, :b => 3),
                             ParamsDict(:c => 5))

        ps[:c] = 6
        ps[:b] = 7

        @test ps[:a] == 1
        @test ps[:c] == 6
        @test ps[:b] == 7
    end

    @testset "Constructing from another parameters object" begin
        ps1 = ParamsDictChain(ParamsDict(:a => 1, :c => 4),
                              ParamsDict(:a => 2, :b => 3))
        ps2 = ParamsDictChain(ParamsDict(:a => 5), ps1,
                              ParamsDict(:c => 6))

        @test ps1[:a] == 1
        @test ps2[:a] == 5
        @test ps2[:c] == 4
    end

    @testset "Get key without default" begin
        ps = ParamsDictChain(ParamsDict(:a => 1, :c => 4),
                             ParamsDict(:a => 2, :b => 3))
        @test get(ps, :a) == 1
        @test get(ps, :b) == 3
        @test get(ps, :d) == nothing
    end

    @testset "Get key with default" begin
        ps = ParamsDictChain(ParamsDict(:a => 1, :c => 4),
                             ParamsDict(:a => 2, :b => 3))
        @test get(ps, :d, 10) == 10
    end

    @testset "Merge with Parameters or Dict" begin
        ps = ParamsDictChain(ParamsDict(:a => 1, :c => 4),
                             ParamsDict(:a => 2, :b => 3))
        ps2 = chain(ps, ParamsDict(:d => 5, :a => 20))
        @test ps2[:d] == 5
        @test ps2[:b] == 3
        @test ps2[:a] == 20
    end
end


text:
#------------------------------------------------------------------------------#
#                                   Package                                    #
#------------------------------------------------------------------------------#

"""
# 1. Description

    module IsoOrthoTensor

Isotropic and orthogonal tensors [1 apud 2] (as Julia multi-dimensional  arrays)
in one to three Euclidean space dimensions.

Symbolic Hermite Polynomial Tensors — or N-dimensional Hermite polynomials [1] —
in one to three Euclidean space dimensions.

# 2. Usage

```julia-repl
julia> import HermitePolyTensor

```

# 3. References

[1]: H. Grad, “Note on N-dimensional  Hermite  polynomials,”  Communications  on
Pure and Applied Mathematics, vol. 2, no. 4, pp. 325–330, 1949.

[2]: K. K.  Mattila,  L.  A.  Hegele  Júnior,  P.  C.  Philippi,  “High-Accuracy
Approximation of High-Rank Derivatives: Isotropic Finite  Differences  Based  of
Lattice-Boltzmann Stencils,” The Scientific World Journal, vol. 2014, article ID
142907, 16 pages, 2014.
"""
module IsoOrthoTensor

#------------------------------------------------------------------------------#
#                                   Imports                                    #
#------------------------------------------------------------------------------#

using Reexport

#------------------------------------------------------------------------------#
#                                   Includes                                   #
#------------------------------------------------------------------------------#

include("version.jl")

include("IOT.jl")
@reexport using .IOT

end # module


text:
mutable struct FPS
    quiet::Bool
    damper::Float64
    print_time::Float64
    name::String
    lasttime::Float64
    lastprinttime::Float64
    fps::Float64
    FPS(;quiet=false, damper=5, print_time  = 0.5 , name = "unnamed") =
        new(quiet,damper,print_time,name,time()-0.1,time()-0.1,0);
end
import Base.tic
function tic(t::FPS)
    printed = false
    curtime = time()
    elapsed = max(curtime - t.lasttime,1e-6)
    t.lasttime = curtime
    elapsedprint = curtime - t.lastprinttime
    curfps = (1.0/elapsed);
    factor = exp(-(abs(curfps-t.fps)/t.damper)^2)
    t.fps = t.fps*factor + curfps*(1-factor)
    if !t.quiet && elapsedprint>=t.print_time
        t.lastprinttime = curtime
        print("[$(t.name)]$(t.fps) fps\n");
        printed = true
    end
    return (t.fps,printed)
end
export FPS

(fps_t::FPS)() = tic(fps_t)

function match_images(vargin...)
    n = length(vargin)
    m = Integer(ceil(n/2))
    n = Integer(ceil(n/m))

    ax_base = PyPlot.subplot(m,n,1)
    ax_base[:imshow](vargin[1])
    for i=2:length(vargin)
        ax2 = PyPlot.subplot(m,n,i, sharex=ax_base,sharey=ax_base)
        ax2[:imshow](vargin[i])
    end
end
export match_images


text:
for (f,acb_f) in [(:dft, :acb_dft), (:inverse_dft, :acb_dft_inverse)]
    @eval begin
        function $f(x::Vector{ArbComplex{P}}) where P
            n = length(x)
            x_c = ArblibVector(x)
            transf_c = ArblibVector(ArbComplex{P}, n)

            # Signature: void acb_dft(acb_ptr w, acb_srcptr v, slong n, slong prec)
            ccall(@libarb($acb_f), Cvoid, (Ref{ArbComplex{P}}, Ref{ArbComplex{P}}, Cint, Cint ), transf_c.ptr, x_c.ptr, n, P)

            free!(x_c)

            # Copy data and release memory
            take!(transf_c)
        end
    end
end

dft(x::Array{ArbComplex,1}) = dft(Array{ArbComplex{workingprecision(ArbComplex)}, 1}(x))
dft(x::Array{ArbReal{P},1}) where {P} = dft(Array{ArbComplex{P}, 1}(x))
dft(x::Array{ArbFloat{P},1}) where {P} = dft(Array{ArbComplex{P}, 1}(x))
dft(x::Array{ArbReal,1}) = dft(Array{ArbComplex{workingprecision(ArbReal)}, 1}(x))
dft(x::Array{ArbFloat,1}) = dft(Array{ArbComplex{workingprecision(ArbFloat)}, 1}(x))

inverse_dft(x::Array{ArbComplex,1}) = inverse_dft(Array{ArbComplex{workingprecision(ArbComplex)}, 1}(x))
inverse_dft(x::Array{ArbReal{P},1}) where {P} = dft(Array{ArbComplex{P}, 1}(x))
inverse_dft(x::Array{ArbFloat{P},1}) where {P} = dft(Array{ArbComplex{P}, 1}(x))
inverse_dft(x::Array{ArbReal,1}) = dft(Array{ArbComplex{workingprecision(ArbReal)}, 1}(x))
inverse_dft(x::Array{ArbFloat,1}) = dft(Array{ArbComplex{workingprecision(ArbFloat)}, 1}(x))


text:
"""
    linearsearch(a::Union{AbstractVector{T}, T}, item::T)

Linear search in array, find first occurance. Time complexity O(n).
"""
function linearsearch(a::Union{AbstractVector{T}, T}, item::T)::Int where T
    if length(a) == 0 return 0 end

    for i in 1:length(a)
        if a[i] == first(item)
            return i
        end
    end

    return 0
end

"""
    binarysearch(a::Union{AbstractVector{T}, T}, item::T)

Binary search in array, find any occurance. Time complexity O(logn).
"""
function binarysearch(a::Union{AbstractVector{T}, T}, item::T) where T
    lo::Int = 1
    hi::Int = length(a)
    if hi == 0 return nothing end

    mid::Int = 0
    item = first(item)
    while lo <= hi
        mid = lo + (hi - lo) ÷ 2

        if item < a[mid]
            hi = mid - 1
        elseif item > a[mid]
            lo = mid + 1
        else
            break
        end
    end

    return item == a[mid] ? mid : nothing
end


text:
"""
    Count{S <: MOI.AbstractScalarSet}(dimension::Int, set::MOI.AbstractScalarSet)

``\\{(y, x) \\in \\mathbb{N} \\times \\mathbb{T}^\\mathtt{dimension} : y = |\\{i | x_i \\in S \\}|\\}``

`dimension` is the number of variables that are checked against the `set`.

Also called `among`.

## Example

    [w, x, y, z] in Count(3, MOI.EqualTo(2.0))
    # w == sum([x, y, z] .== 2.0)
"""
struct Count{S <: MOI.AbstractScalarSet} <: MOI.AbstractVectorSet
    dimension::Int
    set::S
end

function Count(dimension::Int, value::Real)
    return Count(dimension, MOI.EqualTo(value))
end

MOI.dimension(set::Count{S}) where {S} = set.dimension + 1
copy(set::Count{S}) where {S} = Count(set.dimension, copy(set.set))
function Base.:(==)(x::Count{S}, y::Count{S}) where {S}
    return x.dimension == y.dimension && x.set == y.set
end

"""
    CountedValuesType

Kind of values to be counted for a `GlobalCardinality` constraint:

* either the values to count are fixed when creating the set:
  `FIXED_COUNTED_VALUES`
* or the values are themselves variables (typically constrained elsewhere):
  `VARIABLE_COUNTED_VALUES`
"""
@enum CountedValuesType begin
    FIXED_COUNTED_VALUES
    VARIABLE_COUNTED_VALUES
end

"""
    CountedValuesClosureType

Whether values that are not counted in `GlobalCardinality` constraint are
allowed in the array whose values are counted:

* either uncounted values are allowed: `OPEN_COUNTED_VALUES`
* or they are not allowed: `CLOSED_COUNTED_VALUES`
"""
@enum CountedValuesClosureType begin
    OPEN_COUNTED_VALUES
    CLOSED_COUNTED_VALUES
end

"""
    GlobalCardinality{CVT, CVCT, T}(dimension::Int, values::Vector{T})

This set represents the large majority of the variants of the
global-cardinality constraint, with the parameters set in `CountedValuesType`
(`CVT` parameter) and `CountedValuesClosureType` (`CVCT` parameter).

## Fixed and open

``\\{(x, y) \\in \\mathbb{T}^\\mathtt{dimension} \\times \\mathbb{N}^d : y_i = |\\{ j | x_j = \\mathtt{values}_i, \\forall j \\}| \\}``

The first `dimension` variables are an array, the last variables are the
number of times that each item of `values` is present in the first array.
Values that are not in `values` are ignored.

Also called [`gcc`](https://sofdem.github.io/gccat/gccat/Cglobal_cardinality.html)
or `count`.

### Example

    [x, y, z, v, w] in GlobalCardinality{FIXED_COUNTED_VALUES, OPEN_COUNTED_VALUES}(3, [2.0, 4.0])
    [x, y, z, v, w] in GlobalCardinality{OPEN_COUNTED_VALUES}(3, [2.0, 4.0])
    [x, y, z, v, w] in GlobalCardinality(3, [2.0, 4.0])
    # v == sum([x, y, z] .== 2.0)
    # w == sum([x, y, z] .== 4.0)

## Variable and open

``\\{(x, y, z) \\in \\mathbb{T}^\\mathtt{dimension} \\times \\mathbb{N}^\\mathtt{n\\_values} \\times \\mathbb{T}^\\mathtt{n\\_values} : y_i = |\\{ j | x_j = z_i, \\forall j \\}| \\}``

The first `dimension` variables are an array, the next `n_values` variables
are the number of times that each item of the last `n_values` variables is
present in the first array. Values of the first array that are not in the
`n_values` are ignored.

Also called `distribute`.

### Example

    [x, y, z, t, u, v, w] in GlobalCardinality{VARIABLE_COUNTED_VALUES, OPEN_COUNTED_VALUES, T}(3, 2)
    [x, y, z, t, u, v, w] in GlobalCardinality{OPEN_COUNTED_VALUES, T}(3, 2)
    [x, y, z, t, u, v, w] in GlobalCardinality{T}(3, 2)
    # t == sum([x, y, z] .== v)
    # u == sum([x, y, z] .== w)

## Fixed and closed

``\\{(x, y) \\in \\mathbb{T}^\\mathtt{dimension} \\times \\mathbb{N}^d : y_i = |\\{ j | x_j = \\mathtt{values}_i, \\forall j \\}| \\}``

The first `dimension` variables are an array, the last variables are the
number of times that each item of `values` is present in the first array.
Each value of the first array must be within `values`.

### Example

    [x, y, z, v, w] in GlobalCardinality{FIXED_COUNTED_VALUES, CLOSED_COUNTED_VALUES, T}(3, [2.0, 4.0])
    # v == sum([x, y, z] .== 2.0)
    # w == sum([x, y, z] .== 4.0)
    # x ∈ [2.0, 4.0], y ∈ [2.0, 4.0], z ∈ [2.0, 4.0]

## Variable and closed

``\\{(x, y, z) \\in \\mathbb{T}^\\mathtt{dimension} \\times \\mathbb{N}^\\mathtt{n\\_values} \\times \\mathbb{T}^\\mathtt{n\\_values} : y_i = |\\{ j | x_j = z_i, \\forall j \\}| \\}``

The first `dimension` variables are an array, the next `n_values` variables
are the number of times that each item of the last `n_values` variables is
present in the first array. Each value of the first array must be within the
next given `n_values`.

Also called `distribute`.

### Example

    [x, y, z, t, u, v, w] in GlobalCardinality{VARIABLE_COUNTED_VALUES, CLOSED_COUNTED_VALUES, T}(3, 2)
    # t == sum([x, y, z] .== v)
    # u == sum([x, y, z] .== w)
    # x ∈ [v, w], y ∈ [v, w], z ∈ [v, w]
"""
struct GlobalCardinality{CVT, CVCT, T <: Real} <: MOI.AbstractVectorSet
    dimension::Int
    values::Vector{T}
    n_values::Int

    function GlobalCardinality{CVT, CVCT, T}(dimension::Int, values::Vector{T}, n_values::Int) where {CVT, CVCT, T <: Real}
        if CVT == FIXED_COUNTED_VALUES && length(values) == 0
            error("Inconsistent GlobalCardinality set: the counted values should be fixed, but array of values to count is $(values).")
        end
        if CVT == VARIABLE_COUNTED_VALUES && n_values <= 0
            error("Inconsistent GlobalCardinality set: the counted values should be variables, but their number is $(n_values).")
        end

        return new{CVT, CVCT, T}(dimension, values, n_values)
    end
end

# Helper functions: always FIXED_COUNTED_VALUES.
function GlobalCardinality(dimension::Int, values::Vector{T}) where {T <: Real}
    return GlobalCardinality{FIXED_COUNTED_VALUES, OPEN_COUNTED_VALUES, T}(dimension, values, -1)
end

function GlobalCardinality{CVCT}(dimension::Int, values::Vector{T}) where {CVCT, T <: Real}
    return GlobalCardinality{FIXED_COUNTED_VALUES, CVCT, T}(dimension, values, -1)
end

function GlobalCardinality{CVCT, T}(dimension::Int, values::Vector{T}) where {CVCT, T <: Real}
    return GlobalCardinality{FIXED_COUNTED_VALUES, CVCT, T}(dimension, values, -1)
end

function GlobalCardinality{CVT, CVCT}(dimension::Int, values::Vector{T}) where {CVT, CVCT, T <: Real}
    return GlobalCardinality{CVT, CVCT, T}(dimension, values, -1)
end

function GlobalCardinality{CVT, CVCT, T}(dimension::Int, values::Vector{T}) where {CVT, CVCT, T <: Real}
    return GlobalCardinality{CVT, CVCT, T}(dimension, values, -1)
end

# Helper functions: always VARIABLE_COUNTED_VALUES.
function GlobalCardinality{T}(dimension::Int, n_values::Int) where {T <: Real}
    return GlobalCardinality{VARIABLE_COUNTED_VALUES, OPEN_COUNTED_VALUES, T}(dimension, T[], n_values)
end

function GlobalCardinality{CVCT, T}(dimension::Int, n_values::Int) where {CVCT, T <: Real}
    return GlobalCardinality{VARIABLE_COUNTED_VALUES, CVCT, T}(dimension, T[], n_values)
end

function GlobalCardinality{CVT, CVCT, T}(dimension::Int, n_values::Int) where {CVT, CVCT, T <: Real}
    return GlobalCardinality{CVT, CVCT, T}(dimension, T[], n_values)
end

function n_values(set::GlobalCardinality)
    if length(set.values) > 0
        return length(set.values)
    else
        return set.n_values
    end
end

function MOI.dimension(set::GlobalCardinality{CVT, CVCT, T}) where {CVT, CVCT, T <: Real}
    dim = set.dimension + n_values(set) # Array to count in, counts.
    if CVT == VARIABLE_COUNTED_VALUES
        dim += n_values(set) # Values to count.
    end
    return dim
end

function copy(set::GlobalCardinality{CVT, CVCT, T})  where {CVT, CVCT, T <: Real}
    return GlobalCardinality{CVT, CVCT, T}(set.dimension, copy(set.values), set.n_values)
end

function Base.:(==)(x::GlobalCardinality{CVT, CVCT, T}, y::GlobalCardinality{CVT, CVCT, T}) where {CVT, CVCT, T <: Real}
    return x.dimension == y.dimension && x.values == y.values && x.n_values == y.n_values
end

# Shortcuts for types.
# const GlobalCardinality{T} = GlobalCardinality{FIXED_COUNTED_VALUES, OPEN_COUNTED_VALUES, T}
const GlobalCardinalityOpen{T} = GlobalCardinality{FIXED_COUNTED_VALUES, OPEN_COUNTED_VALUES, T}
const GlobalCardinalityClosed{T} = GlobalCardinality{FIXED_COUNTED_VALUES, CLOSED_COUNTED_VALUES, T}
const GlobalCardinalityFixed{T} = GlobalCardinality{FIXED_COUNTED_VALUES, OPEN_COUNTED_VALUES, T}
const GlobalCardinalityFixedOpen{T} = GlobalCardinality{FIXED_COUNTED_VALUES, OPEN_COUNTED_VALUES, T}
const GlobalCardinalityFixedClosed{T} = GlobalCardinality{FIXED_COUNTED_VALUES, CLOSED_COUNTED_VALUES, T}
const GlobalCardinalityVariable{T} = GlobalCardinality{VARIABLE_COUNTED_VALUES, OPEN_COUNTED_VALUES, T}
const GlobalCardinalityVariableOpen{T} = GlobalCardinality{VARIABLE_COUNTED_VALUES, OPEN_COUNTED_VALUES, T}
const GlobalCardinalityVariableClosed{T} = GlobalCardinality{VARIABLE_COUNTED_VALUES, CLOSED_COUNTED_VALUES, T}

# Ease dispatch for implementing solvers.
function MOI.supports_constraint(o::MOI.AbstractOptimizer, f::MOI.AbstractVectorFunction, s::Type{GlobalCardinality{CVT, CVCT, T}}) where {CVT, CVCT, T <: Real}
    return MOI.supports_constraint(o, f, s, Val(CVT), Val(CVCT))
end

function MOI.supports_constraint(o::MOI.AbstractOptimizer, f::MOI.AbstractVectorFunction, s::Type{GlobalCardinality{CVT, OPEN_COUNTED_VALUES, T}}, ::Val{CVT}, ::Val{OPEN_COUNTED_VALUES}) where {CVT, T <: Real}
    return MOI.supports_constraint(o, f, s, Val(CVT))
end

function MOI.supports_constraint(o::MOI.AbstractOptimizer, f::MOI.AbstractVectorFunction, s::Type{GlobalCardinality{FIXED_COUNTED_VALUES, CVCT, T}}, ::Val{FIXED_COUNTED_VALUES}, ::Val{CVCT}) where {CVCT, T <: Real}
    return MOI.supports_constraint(o, f, s, Val(CVCT))
end

"""
    CountCompare(dimension::Int)

``\\{(z, x, y) \\in \\mathbb{N} \\times \\mathbb{R}^\\mathtt{dimension} \\times \\mathbb{R}^\\mathtt{dimension} : Z = |\\{i | x_i = y_i\\}|\\}``

The first `dimension` variables are the first array that is compared to the
second one, indicated by the next `dimension` variables. The last variable
is the number of values that are identical in both arrays.

## Example

    [v, w, x, y, z] in Count(2)
    # w == sum([w, x] .== [y, z])
"""
struct CountCompare <: MOI.AbstractVectorSet
    dimension::Int
end

MOI.dimension(set::CountCompare) = 2 * set.dimension + 1

"""
    CountDistinct(dimension::Int)

The first variable in the set is forced to be the number of distinct values in
the rest of the expressions.

This is a relaxed version of `AllDifferent`; it encodes an `AllDifferent`
constraint when the first variable is the number of variables in the set.

Also called `nvalues`.

## Example

    [x, y, z] in CountDistinct(3)
    # x = 1 if y == z, x = 2 if y != z
"""
struct CountDistinct <: MOI.AbstractVectorSet
    dimension::Int
end

MOI.dimension(set::CountDistinct) = set.dimension + 1

# isbits types, nothing to copy
function copy(
    set::Union{CountCompare, CountDistinct},
)
    return set
end


text:
#============================
horizontal maps / (x,y) plots
============================#
@userplot HorizontalPlane
@recipe function f(p::HorizontalPlane)
    x, y, z = p.args
    @series begin
        seriestype --> :heatmap
        xguide --> "Longitude"
        yguide --> "Latitude"
        x, y, z
    end
end
"""
    plothorizontalslice(x, grd; depth)

Plots a heatmap of the horizontal slice of tracer `x` at depth `depth`.
"""
plothorizontalslice(x, grd; depth=nothing, kwargs...) = horizontalplane(grd.lon, grd.lat, horizontalslice(x, grd; depth=depth); kwargs...)
plothorizontalslice!(x, grd; depth=nothing, kwargs...) = horizontalplane!(grd.lon, grd.lat, horizontalslice(x, grd; depth=depth); kwargs...)
plothorizontalslice!(plt, x, grd; depth=nothing, kwargs...) = horizontalplane!(plt, grd.lon, grd.lat, horizontalslice(x, grd; depth=depth); kwargs...)
"""
    surfacemap(x, grd)

Plots a surface heatmap of tracer `x`.
"""
surfacemap(args...; kwargs...) = plothorizontalslice(args...; kwargs..., depth=0)
surfacemap!(args...; kwargs...) = plothorizontalslice!(args...; kwargs..., depth=0)
"""
    plotverticalintegral(x, grd, mask=1)

Plots the vertical integral of tracer `x`.
"""
plot∫dz(x, grd; mask=1, kwargs...) = horizontalplane(grd.lon, grd.lat, ∫dz(x, grd, mask); kwargs...)
plot∫dz!(plt, x, grd; mask=1, kwargs...) = horizontalplane!(plt, grd.lon, grd.lat, ∫dz(x, grd, mask); kwargs...)
plotverticalintegral = plot∫dz
plotverticalintegral! = plot∫dz!
"""
    plotverticalaverage(x, grd, mask=1)

Plots the vertical average of tracer `x`.
"""
plotverticalmean(x, grd; mask=1, kwargs...) = horizontalplane(grd.lon, grd.lat, verticalmean(x, grd, mask); kwargs...)
plotverticalaverage = plotverticalmean
export plotverticalaverage, plotverticalmean, plotverticalintegral, plotverticalintegral!, plot∫dz, plot∫dz!, surfacemap, surfacemap!, plothorizontalslice

"""
    minimap(grd; central_longitude=200°)

Plots a surface map of grd with no ticks, labels, or colorbar.

The goal of this function is to provide a way to plot both a minimap and a
cruise track when plotting a transect:
```
plottransect(dummy, grd; ct=sort(ct))
plot!(inset_subplots=bbox(0.73,0.73,0.15,0.15), subplot=2)
minimap!(grd; subplot=2)
cruisetrackplots!(sort(ct), subplot=2)
```
"""
@userplot MiniMap
@recipe function f(p::MiniMap; central_longitude=200°)
    grd, = p.args
    wlon = central_longitude - 180°
    xgrd = @. mod(grd.lon - wlon, 360°) + wlon
    ix = sortperm(xgrd)
    @series begin
        seriestype --> :heatmap
        seriescolor --> :GnBu_3
        rev --> true
        xguide --> ""
        yguide --> ""
        framestyle  --> false
        showaxis --> false
        legend --> :none
        colorbar --> false
        grid --> false
        xgrd[ix], grd.lat, view(iswet(grd),:,ix,1)
    end
end

#================================
Vertical–meridional / (y,z) plots
================================#
@userplot MeridionalPlane
@recipe function f(p::MeridionalPlane)
    y, z, v = p.args
    @series begin
        seriestype --> :heatmap
        xguide --> "Latitude"
        yguide --> "Depth"
        yflip --> true
        y, z, v
    end
end
"""
    plotmeridionalslice(x, grd; lon)

Plots a meridional slice of tracer `x` at longitude `lon`.
"""
plotmeridionalslice(x, grd; lon=nothing, kwargs...) = meridionalplane(grd.lat, grd.depth, meridionalslice(x, grd; lon=lon); kwargs...)
"""
    plotzonalmean(x, grd; mask=1)

Plots a zonal average of tracer `x`.
"""
plotzonalmean(x, grd; mask=1, kwargs...) = meridionalplane(grd.lat, grd.depth, zonalmean(x, grd, mask)'; kwargs...)
plotzonalaverage = plotzonalmean
"""
    plotzonalintegral(x, grd; mask=1)

Plots a zonal integral of tracer `x`.
"""
plot∫dx(x, grd; mask=1, kwargs...) = meridionalplane(grd.lat, grd.depth, ∫dx(x, grd, mask)'; kwargs...)
plotzonalintegral = plot∫dx
export plotzonalintegral, plot∫dx, plotzonalmean, plotzonalaverage, plotmeridionalslice

#===========================
Vertical–zonal / (x,z) plots
===========================#
@userplot ZonalPlane
@recipe function f(p::ZonalPlane)
    x, z, v = p.args
    @series begin
        seriestype --> :heatmap
        xguide --> "Longitude"
        yguide --> "Depth"
        yflip --> true
        x, z, v
    end
end
"""
    plotzonalslice(x, grd; lat)

Plots a zonal slice of tracer `x` at latitude `lat`.
"""
plotzonalslice(x, grd; lat=nothing, kwargs...) = zonalplane(grd.lon, grd.depth, zonalslice(x, grd; lat=lat); kwargs...)
"""
    plotmeridionalmean(x, grd; mask=1)

Plots a meridional average of tracer `x`.
"""
plotmeridionalmean(x, grd; mask=1, kwargs...) = zonalplane(grd.lon, grd.depth, meridionalmean(x, grd, mask)'; kwargs...)
plotmeridionalaverage = plotmeridionalmean
"""
    plotmeridionalintegral(x, grd; mask=1)

Plots a meridional integral of tracer `x`.
"""
plot∫dy(x, grd; mask=1, kwargs...) = zonalplane(grd.lon, grd.depth, ∫dy(x, grd, mask)'; kwargs...)
plotmeridionalintegral = plot∫dy
export plotmeridionalintegral, plot∫dy, plotmeridionalmean, plotmeridionalaverage, plotzonalslice

#=================
Vertical / z plots
=================#
@userplot ZPlot
@recipe function f(p::ZPlot)
    x, z = p.args
    @series begin
        yguide --> "Depth"
        yflip --> true
        x, z
    end
end

"""
    horizontalintegral(x, grd; mask=1)

Plots a horizontal integral of tracer `x`.
"""
plot∫dxdy(x, grd; mask=1, kwargs...) = zplot(∫dxdy(x, grd, mask), grd.depth; kwargs...)
plot∫dxdy!(x, grd; mask=1, kwargs...) = zplot!(∫dxdy(x, grd, mask), grd.depth; kwargs...)
plot∫dxdy!(plt, x, grd; mask=1, kwargs...) = zplot!(plt, ∫dxdy(x, grd, mask), grd.depth; kwargs...)
plothorizontalintegral = plot∫dxdy
plothorizontalintegral! = plot∫dxdy!
"""
    plothorizontalaverage(x, grd; mask=1)

Plots a horizontal average of tracer `x`.
"""
plothorizontalmean(x, grd; mask=1, kwargs...) = zplot(horizontalmean(x, grd, mask), grd.depth; kwargs...)
plothorizontalmean!(plt, x, grd; mask=1, kwargs...) = zplot!(plt, horizontalmean(x, grd, mask), grd.depth; kwargs...)
plothorizontalaverage = plothorizontalmean
plothorizontalaverage! = plothorizontalmean!
"""
    plotdepthprofile(x, grd; lonlat)

Plots the profile of tracer `x` interpolated at `lonlat=(x,y)` coordinates.
"""
plotdepthprofile(x, grd; lonlat=nothing, kwargs...) = zplot(depthprofile(x, grd; lonlat=lonlat), grd.depth; kwargs...)
export plot∫dxdy, plot∫dxdy!, plothorizontalintegral, plothorizontalintegral!, plothorizontalmean, plothorizontalmean!, plothorizontalaverage, plothorizontalaverage!, plotdepthprofile

#======================================================
Vertical / (dist,z) plots (for transect/section slices)
======================================================#
@userplot VerticalPlane
@recipe function f(p::VerticalPlane)
    dist, z, v = p.args
    @series begin
        seriestype --> :heatmap
        xguide --> "Distance"
        yguide --> "Depth"
        yflip --> true
        grid --> false
        dist, z, v
    end
end
"""
    plottransect(x, grd; ct=ct)

Plots the transect of tracer `x` along transect `ct`.
"""
plottransect(x, grd; ct=nothing, start=:south, show_stations=false, insetmap=false, kwargs...) = verticalplane(verticalsection(x, grd; ct=ct)...; kwargs...)
export plottransect

"""
    RatioAtStation(x, y, grd, station, depthlims=(0,Inf))

Plots a meridional transect of tracer `x` along cruise track `ct`.

The keyword argument `zlims=(ztop, zbottom)` can be provided
if you only want to only plot for depths `z ∈ (ztop, zbottom)`.
(`z` is positive downwards in this case)
"""
@userplot RatioAtStation
@recipe function f(p::RatioAtStation; depthlims=(0,Inf))
    depthlims = OceanGrids.convertdepth.(depthlims)
    x, y, grd, st = p.args
    x3D = rearrange_into_3Darray(x, grd)
    y3D = rearrange_into_3Darray(y, grd)
    ndepth = length(grd.depth)
    knots = (grd.lat, grd.lon, 1:ndepth)
    itpx = interpolate(knots, x3D, (Gridded(Linear()), Gridded(Linear()), NoInterp()))
    itpy = interpolate(knots, y3D, (Gridded(Linear()), Gridded(Linear()), NoInterp()))
    iz = findall(depthlims[1] .≤ grd.depth .≤ depthlims[2])
    ibot = findfirst(grd.depth .> depthlims[2])
    x1D = itpx(convertlat(st.lat), mod(convertlon(st.lon), 360°), iz)
    y1D = itpy(convertlat(st.lat), mod(convertlon(st.lon), 360°), iz)
    @series begin
        seriestype := :scatter
        label --> st.name
        markershape --> :circle
        seriescolor --> :deep
        marker_z --> ustrip.(-grd.depth)
        colorbar_title --> "Depth (m)" # TODO fix in UnitfulRecipes
        markersize --> 4
        legend --> :topleft
        x1D, y1D
    end
end

"""
    PlotParameter(p::AbstractParameters, s)

Plots the PDF of parameter `p` with symbol `s`
"""
@userplot PlotParameter
@recipe function f(plt::PlotParameter)
    d, v, initv, s, u = extract_dvisu(plt.args...)
    xs = default_range(d)
    ys = [pdf(d, x) for x in xs] # the PDF
    xu = xs * upreferred(u) .|> u .|> ustrip
    vu = v * upreferred(u) |> u |> ustrip
    @series begin
        label --> "prior"
        xu, ys
    end
    @series begin
        label --> "initial value"
        initv * [1,1], [0, pdf(d, v)]
    end
    @series begin
        xguide --> "$s ($(string(u)))"
        yguide --> "PDF"
        label --> "value"
        markershape --> :o
        [vu], [pdf(d, v)]
    end
end

"""
    PlotParameters(p::AbstractParameters)

Plots the PDF of all the flattenable parameters in `p`.
"""
@userplot PlotParameters
@recipe function f(plt::PlotParameters)
    p, = plt.args
    layout --> (length(p),1)
    for (i,s) in enumerate(flattenable_symbols(p))
        d, v, initvu, s, u = extract_dvisu(p,s)
        xs = default_range(d)
        ys = [pdf(d, x) for x in xs] # the PDF
        xu = xs * upreferred(u) .|> u .|> ustrip
        vu = v * upreferred(u) |> u |> ustrip
        initv = ustrip(upreferred(initvu * units(p,s)))
        @series begin
            label --> "prior"
            subplot := i
            xu, ys
        end
        @series begin
            label --> "initial value"
            subplot := i
            initvu * [1,1], [0, pdf(d, initv)]
        end
        @series begin
            xguide --> "$s ($(string(u)))"
            yguide --> "PDF"
            label --> "value"
            markershape --> :o
            subplot := i
            [vu], [pdf(d, v)]
        end
    end
end

extract_dvisu(d, v, initv, s, u) = d, v, initv, s, u
extract_dvisu(p::AbstractParameters, s) = prior(p,s), UnPack.unpack(p, Val(s)), initial_value(p,s), s, units(p,s)

function default_range(d::Distribution, n = 4)
    μ, σ = mean(d), std(d)
    xmin = max(μ - n*σ, minimum(support(d)))
    xmax = min(μ + n*σ, maximum(support(d)))
    range(xmin, xmax, length=1001)
end

#======================================================
Diagnostic recipes
======================================================#
"""
    PlotStencil(st)

Plots the stencil `st`.
"""
@userplot PlotStencil
@recipe function f(p::PlotStencil)
    st, = p.args
    @series begin
        linecolor --> :black
        label --> ""
        convert.(Tuple, vcat([[0*x, x, 0*x] for x in st]...))
    end
    @series begin
        marker_z --> [x.I[3] for x in st]
        label --> ""
        xguide --> "lat index"
        yguide --> "lon index"
        zguide --> "depth index"
        seriescolor --> :darkrainbow
        linewidth --> 0
        markershape --> :circle
        convert.(Tuple, st)
    end
end


text:
using BasisFunctions

export FourierDomain
"""
A FourierDomain is a domain that is implicitly defined by a Fourier series.
Example: if relop is x -> x-2 > 0, then the Fourier domain is the domain where the given
Fourier series is greater than 2.
"""
struct FourierDomain{T} <: Domain{T}
    f       ::  Expansion
    relop	::	Function
end

function indomain(x, d::FourierDomain)
	z = d.f(x)
	in(x,domain(d.f)) && d.relop(z)
end

# TODO: add postprocessing so that x outside domain -> NaN
function in(g::AbstractGrid, d::FourierDomain)
	z = d.f(g)[1]
	map(d.relop, z)
end

(<)(f::Expansion{T}, a::Number) where {T} = FourierDomain{convert(SVector,T)}(f, x-> x < a)
(<=)(f::Expansion{T}, a::Number) where {T} = FourierDomain{convert(SVector,T)}(f, x-> x <= a)
(>)(f::Expansion{T}, a::Number) where {T} = FourierDomain{convert(SVector,T)}(f, x-> x > a)
(>=)(f::Expansion{T}, a::Number) where {T} = FourierDomain{convert(SVector,T)}(f, x-> x >= a)

boundingbox(d::FourierDomain) = boundingbox(domain(d.f))

export ComparisonDomain
"""
A ComparisonDomain is a domain that is implicitly defined by two Fourier series.
Example: if binop is (x,y) -> x > y, then the Fourier domain is the domain where the first given
Fourier series is greater than the second.
"""
struct ComparisonDomain{T} <: Domain{T}
    f       ::  Expansion
    g		::	Expansion
    binop	::	Function
end

function indomain(x, d::ComparisonDomain)
	z1 = d.f(x)
	z2 = d.g(x)
	in(x,domain(d.f)) && in(x,domain(d.g)) && d.binop(z1, z2)
end

# TODO: add postprocessing so that x outside domain -> NaN
function in(g::AbstractGrid, d::ComparisonDomain)
	z1 = d.f(g)
	z2 = d.g(g)
	in(g,domain(d.f)) & in(g,domain(d.g)) & map(d.binop, z1, z2)
end

(<)(f::Expansion{T}, g::Expansion{T}) where {T} = ComparisonDomain{convert(SVector,T)}(f, g, (x,y)-> x < y)
(<=)(f::Expansion{T}, g::Expansion{T}) where {T} = ComparisonDomain{convert(SVector,T)}(f, g, (x,y)-> x <= y)
(>)(f::Expansion{T}, g::Expansion{T}) where {T} = ComparisonDomain{convert(SVector),T}(f, g, (x,y)-> x > y)
(>=)(f::Expansion{T}, g::Expansion{T}) where {T} = ComparisonDomain{convert(SVector,T)}(f, g, (x,y)-> x >= y)

boundingbox(d::ComparisonDomain) = boundingbox(domain(d.f))+boundingbox(domain(d.g))


text:
function compress{T <: Float}(
    x :: DenseMatrix{T},
)
    # get size of matrix
    (n,p) = size(x)

    # get designated blocksizes
    blocksize = ((n-1) >> 2) + 1
    tblocksize = ((p-1) >> 2) + 1

    # output vector
    outlen = blocksize*p
    y = zeros(Int8, outlen+3)

    # first three bytes are fixed
    y[1] = MNUM1
    y[2] = MNUM2
    y[3] = ONE8

    # counters for bookkeeping
    ngeno = 0
    ybyte = 4   # start after magic numbers!

    # loop over all columns in x
    for j = 1:p

        # each column starts with fresh counter
        # ybyte is offset by first 3 bookkeeping bytes plus previous columns
        ngeno = 0
        yidx = 4 + (j-1)*blocksize
        ybyte = ZERO8

        # loop over cases
        for i = 1:n
            ngeno += 1
#            println("ngeno = $ngeno, ybyte = $ybyte, i = $i, j = $j, x[i,j] = $(x[i,j])")
            if T == Float64
#                y[yidx] = y[yidx] << 2 + bin64[x[i,j]]
#                ybyte = ybyte << 2 + bin64[x[i,j]]
                ybyte = bin64[x[i,j]] << 2*(ngeno-1) | ybyte
#                println("ybyte = $ybyte")
            else
#                y[yidx] = y[yidx] << 2 + bin32[x[i,j]]
                ybyte = bin32[x[i,j]] << 2*(ngeno-1) | ybyte
            end

            # reset whenver we reach 4 genotypes
            # OR when we reach end of a column
            if ngeno >= 4 || i == n
                y[yidx] = ybyte
                yidx += 1
                ngeno = 0
                ybyte = ZERO8
            end
        end
    end
    return y
end

max_stars_repo_name: klkeys/PLINK

max_issues_repo_name: klkeys/PLINK


text:
module AudioToaster

#=
n == "connectivity.CopperTrace" ||
n == "connectivity.GroundFill" ||

n == "semiconductors.DIL__IC" ||

n == "semiconductors.ICSymbol" ||

elseif n == "passive.PotentiometerPanel" ||
    n == "passive.TrimmerPotentiometer" ||
    n == "passive.PotentiometerSymbol"
=#

export process, loadfile

using LightXML
using ACME
using WAV

struct Point
    x::Real
    y::Real
end
Point(el::XMLElement) = Point(
    parse(Float64, attribute(el, "x")),
    parse(Float64, attribute(el, "y")))
Point(xy::Union{Array,Tuple}) = Point(xy[1], xy[2])

struct Rectangle
    topleft::Point
    bottomright::Point
end

"""
`Transistor(description, pinout)` where description is an Expr returning an
ACME.Element when evaluated and pinout is a permutation that will be
applied to `[:collector, :base, :emitter]`.

Example:
```julia
t = Transistor(:(bjt(:npn, βf=200)), [1, 2, 3])
```
"""
struct Transistor
    description::Expr
    pinout::Array{Int}
    Transistor(d, p) = isperm(p) && length(p) == 3 ? new Transistor(d, p) : error("Invalid permutation")
end

struct VeroBoard
    area::Rectangle
    horizontal::Bool
end
VeroBoard(c::XMLElement) = VeroBoard(
    Rectangle(Point(c["firstPoint"]),
    Point(c["secondPoint"])),
    content(c["orientation"]) == "HORIZONTAL")

struct TriPadBoard
    area::Rectangle
    horizontal::Bool
end
TriPadBoard(c::XMLElement) = TriPadBoard(
    Rectangle(Point(c["firstPoint"]),
    Point(c["secondPoint"])),
    content(c["orientation"]) == "HORIZONTAL")

struct Breadboard
    location::Point
    fullsize::Bool
    offset::Bool
    orientation::String
end
Breadboard(c::XMLElement) = Breadboard(
    Point(c["point"]),
    content(c["breadboardSize"]) != "Half",
    content(c["powerStripPosition"]) != "Inline"
    content(c["orientation"]))

struct TraceCut
    location::Point
    betweentraces::Bool
end

struct Wire
    startpoint::Point
    endpoint::Point
end
Wire(ps::Union{Array{Point},Tuple{Point,Point}}) = Wire(ps[1], ps[2])

struct CopperTrace
    area::Rectangle
    pins::Array{Point}
end
CopperTrace(a::Rectangle) = CopperTrace(a, [])

"""
Given an ACME.DiscreteModel and the filepath of a WAV
applies the circuit to the WAV and returns a new WAV.
"""
function process(model::DiscreteModel, inputfile::String, outputfile::String)
    y, fs, nbits = wavread(inputfile)
    processed = run!(model, reshape(y, 1, :))
    wavwrite(reshape(processed, length(processed), 1), outputfile, Fs=fs, nbits=nbits)
end

function process(circuit::Circuit, inputfile::String, outputfile::String)
    y, fs, nbits = wavread(inputfile)
    model = DiscreteModel(circuit, 1/fs)
    processed = run!(model, reshape(y, 1, :))
    wavwrite(reshape(processed, length(processed), 1), outputfile, Fs=fs, nbits=nbits)
    return model
end

function iswithin(rect::Rectangle, point::Point)::Bool
    return rect.topleft.x < point.x &&
        rect.bottomright.x > point.x &&
        rect.topleft.y < point.y &&
        rect.bottomright.y > point.y
end

function findcomponentbyname(haystack::XMLElement, needle::Union{Regex,String})
    for component in child_elements(haystack)
        if occursin(needle, content(component["name"]))
            return component
        end
    end
end

"""
Connects a single pin on a component to any other pin occupying the same space.
"""
function addconnection!(circ, col, p::Point, pin)
    haskey(col, p) ? connect!(circ, pin, col[p]) : col[p] = pin
end

"""
Connects both pins of a two legged components, positive then negative on polar components.
"""
function addtwoleggedcomponent!(circ, col, comp, ps::Union{Array{Point},Tuple{Point,Point}})
    addconnection!(circ, col, ps[1], (comp, (1, 1)))
    addconnection!(circ, col, ps[2], (comp, (1, -1)))
end

function addtransistor!(circuit, col, transistor, pins, orientation)
    for i in 1:3
        addconnection!(circuit, col, transistor, pins[i], orientation[i])
    end
end

"""
Given an element loaded from a diy file, returns the value of that element
i.e. Returns the resistance of resistors, the capacitance of capacitors, etc
"""
function getvalue(el::XMLElement)::Float64
    return convertunitmagnitude(
        parse(Float64, content(el["value"]["value"])),
        content(el["value"]["unit"])[1])
end

function convertunitmagnitude(value::Real, unit::Char)::Real
    unitdict = Dict{Char, Float64}(
        'K'=>1e3,
        'k'=>1e3,
        'M'=>1e6,
        'p'=>1e-12,
        'n'=>1e-9,
        'u'=>1e-6,
        'm'=>1e-3
    )
    return haskey(unitdict, unit) ? value * unitdict[unit] : value
end

function convertunitmagnitude(value::Real, unit::String)::Real
    return length(unit) == 0 ? value : convertunitmagnitude(value, unit[1])
end

function connectcopperstrips(circ::Circuit, connections, board::VeroBoard; tracecuts)

end

function connectcopperstrips(circ::Circuit, connections, board::TriPadBoard; tracecuts)

end

function connectcopperstrips(circ::Circuit, connections, board::Breadboard; tracecuts)
    strips = []

end

"""
Loads a circuit from a .diy file
"""
function loadfile(filepath::String)
    circ = Circuit()
    boards = []
    wires = []
    tracecuts = []
    connections = Dict{Point, Union{Tuple,Symbol}}()
    transistors = Dict{String, Transistor}(
        "2N3904" => Transistor(:(bjt(:npn, βf=200)), [1, 2, 3]),
        "2N3906" => Transistor(:(bjt(:pnp, βf=200)), [1, 2, 3]),
        "BC546" => Transistor(:(bjt(:npn, βf=400)), [1, 2, 3]),
        "BC546A" => Transistor(:(bjt(:npn, βf=200)), [1, 2, 3]),
        "BC546B" => Transistor(:(bjt(:npn, βf=450)), [1, 2, 3]),
        "BC546C" => Transistor(:(bjt(:npn, βf=600)), [1, 2, 3]),
        "BC547" => Transistor(:(bjt(:npn, βf=400)), [1, 2, 3]),
        "BC547A" => Transistor(:(bjt(:npn, βf=200)), [1, 2, 3]),
        "BC547B" => Transistor(:(bjt(:npn, βf=450)), [1, 2, 3]),
        "BC547C" => Transistor(:(bjt(:npn, βf=600)), [1, 2, 3]),
        "BC548" => Transistor(:(bjt(:npn, βf=400)), [1, 2, 3]),
        "BC548A" => Transistor(:(bjt(:npn, βf=200)), [1, 2, 3]),
        "BC548B" => Transistor(:(bjt(:npn, βf=450)), [1, 2, 3]),
        "BC548C" => Transistor(:(bjt(:npn, βf=600)), [1, 2, 3]),
        "BC549" => Transistor(:(bjt(:npn, βf=400)), [1, 2, 3]),
        "BC549A" => Transistor(:(bjt(:npn, βf=200)), [1, 2, 3]),
        "BC549B" => Transistor(:(bjt(:npn, βf=450)), [1, 2, 3]),
        "BC549C" => Transistor(:(bjt(:npn, βf=600)), [1, 2, 3]),
        "BC550" => Transistor(:(bjt(:npn, βf=400)), [1, 2, 3]),
        "BC550A" => Transistor(:(bjt(:npn, βf=200)), [1, 2, 3]),
        "BC550B" => Transistor(:(bjt(:npn, βf=450)), [1, 2, 3]),
        "BC550C" => Transistor(:(bjt(:npn, βf=600)), [1, 2, 3],
        "BC556" => Transistor(:(bjt(:pnp, βf=400)), [1, 2, 3]),
        "BC556A" => Transistor(:(bjt(:pnp, βf=200)), [1, 2, 3]),
        "BC556B" => Transistor(:(bjt(:pnp, βf=450)), [1, 2, 3]),
        "BC556C" => Transistor(:(bjt(:pnp, βf=600)), [1, 2, 3]),
        "BC557" => Transistor(:(bjt(:pnp, βf=400)), [1, 2, 3]),
        "BC557A" => Transistor(:(bjt(:pnp, βf=200)), [1, 2, 3]),
        "BC557B" => Transistor(:(bjt(:pnp, βf=450)), [1, 2, 3]),
        "BC557C" => Transistor(:(bjt(:pnp, βf=600)), [1, 2, 3]),
        "BC558" => Transistor(:(bjt(:pnp, βf=400)), [1, 2, 3]),
        "BC558A" => Transistor(:(bjt(:pnp, βf=200)), [1, 2, 3]),
        "BC558B" => Transistor(:(bjt(:pnp, βf=450)), [1, 2, 3]),
        "BC558C" => Transistor(:(bjt(:pnp, βf=600)), [1, 2, 3]),
        "BC559" => Transistor(:(bjt(:pnp, βf=400)), [1, 2, 3]),
        "BC559A" => Transistor(:(bjt(:pnp, βf=200)), [1, 2, 3]),
        "BC559B" => Transistor(:(bjt(:pnp, βf=450)), [1, 2, 3]),
        "BC559C" => Transistor(:(bjt(:pnp, βf=600)), [1, 2, 3]),
        "BC560" => Transistor(:(bjt(:pnp, βf=400)), [1, 2, 3]),
        "BC560A" => Transistor(:(bjt(:pnp, βf=200)), [1, 2, 3]),
        "BC560B" => Transistor(:(bjt(:pnp, βf=450)), [1, 2, 3]),
        "BC560C" => Transistor(:(bjt(:pnp, βf=600)), [1, 2, 3])
    )
    file = parse_file(filepath)
    fr = root(file)
    gridspacing = parse(Float64, content(fr["gridSpacing"]["value"]))
    for c in child_elements(fr["components"])
        n = SubString(name(c), 22)

        # store boards until all components are collected
        if n == "boards.VeroBoard"
            push!(boards, VeroBoard(c))

        elseif n == "boards.TriPadBoard"
            push!(boards, TriPadBoard(c))

        elseif n == "boards.Breadboard"
            push!(boards, Breadboard(c))

        elseif startswith("passive.Resistor", n)
            addtwoleggedcomponent!(circ, connections,
                add!(circ, resistor(getvalue(c))),
                map(Point, child_elements(c["points"])))

        elseif n == "passive.RadialCeramicDiskCapacitor" ||
        n == "passive.AxialElectrolyticCapacitor" ||
        n == "passive.RadialElectrolytic" ||
        n == "passive.AxialFilmCapacitor" ||
        n == "passive.RadialFilmCapacitor" ||
        n == "passive.CapacitorSymbol"
            addtwoleggedcomponent!(circ, connections,
                add!(circ, capacitor(getvalue(c))),
                map(Point, child_elements(c["points"])))

        elseif n == "connectivity.HookupWire"
            pointelems = collect(child_elements(c["controlPoints"]))
            points = (Point(pointelems[1]), Point(pointelems[4]))
            if occursin(r"g(rou)?nd"i, content(c["name"]))
                for p in points
                    addconnection!(circ, connections, p, :gnd)
                end
            else
                push!(wires, Wire(points))
            end

        elseif n == "connectivity.Jumper" || n == "connectivity.Line"
            push!(wires, Wire(map(Point, child_elements(c["points"]))))

        # treat terminal strips like a bunch of wires
        elseif n == "boards.TerminalStrip"
            points = map(Point, collect(child_elements(c["controlPoints"])))
            terminalcount = parse(Int, content(c["terminalCount"]))
            for t in 1:terminalcount
                push!(wires, Wire(points[t], points[t + terminalcount]))
            end

        elseif n == "electromechanical.OpenJack1__4" ||
        n == "electromechanical.CliffJack1__4" ||
        n == "electromechanical.ClosedJack1__4"
            jackname = content(c["name"])
            ps = map(Point, collect(child_elements(c["controlPoints"])))
            if occursin(r"output"i, jackname)
                addtwoleggedcomponent!(circ, connections, add!(circ, voltageprobe()), ps)
            elseif occursin(r"input"i, jackname)
                addtwoleggedcomponent!(circ, connections, add!(circ, voltagesource()), ps)
            end

        elseif n == "electromechanical.PlasticDCJack"
            dcmatch = match(r"(-?[\d\.]+)([A-Za-z]*)", content(c["value"]))
            dcval = convertunitmagnitude(parse(Float64, dcmatch[1]), dcmatch[2])
            points = map(Point, collect(child_elements(c["controlPoints"])))
            ps = content(c["polarity"]) == "CENTER_NEGATIVE" ?
                (points[1], points[3]) : (points[3], points[1])
            addtwoleggedcomponent!(circ, connections, add!(circ, voltagesource(dcval)), ps)

        elseif n == "misc.BatterySymbol"
            value = parse(Float64, content(el["voltageNew"]["value"]))
            unit = content(el["voltageNew"]["unit"])[1]
            addtwoleggedcomponent(circ, connections,
                add!(circ, voltagesource(convertunitmagnitude(value, unit))),
                map(Point, child_elements(c["points"])))

        elseif n == "misc.GroundSymbol"
            addconnection!(circ, connections, Point(c["point"]), :gnd)

        elseif n == "semiconductors.DiodeSymbol" ||
        n == "semiconductors.SchottkyDiodeSymbol" ||
        n == "semiconductors.ZenerDiodeSymbol" ||
        n == "semiconductors.DiodeGlass" ||
        n == "semiconductors.DiodePlastic" ||
        n == "semiconductors.LEDSymbol" ||
        n == "semiconductors.LED"
            # TODO: load diode specifics from diode name
            # for now, default diode will do
            addtwoleggedcomponent!(circ, connections,
                add!(circ, diode()),
                map(Point, child_elements(c["points"])))

        elseif n == "connectivity.TraceCut"
            push!(tracecuts, TraceCut(
                Point(c["point"]),
                content(c["cutBetweenHoles"]) == "true"))

        elseif n == "passive.InductorSymbol"
            addtwoleggedcomponent!(circ, connections,
                add!(circ, inductor(getvalue(c))),
                map(Point, child_elements(c["points"])))

        elseif n == "semiconductors.BJTSymbol"
            val = uppercase(content(c["value"]))
            t = if haskey(transistors, val)
                    transistors[val].description
                else
                    bjt(content(c["polarity"]) == "NPN" ? :npn : :pnp)
                end
            addtransistor!(circ, connections,
                add!(circ, t),
                map(Point, collect(child_elements(c["controlPoints"]))),
                content(c["flip"]) == "Y"
                    ? [:base, :collector, :emitter]
                    : [:base, :emitter, :collector])

        elseif n == "semiconductors.TransistorTO92" ||
        n == "semiconductors.TransistorTO1" ||
        n == "semiconductors.TransistorTO220"
            points = map(Point, collect(child_elements(c["points"])))
            pins = [:collector, :base, :emitter]
            value = uppercase(content(c["value"])))
            if haskey(transistors, value)
                tr = add!(circ, transistors[value].description)
                permute!(pins, transistors[value].pinout)
            else
                tr = add!(circ,
                    bjt(occursin(r"pnp"i, value) ? :pnp : :npn,
                        βf = occursin(r"gain ?= ?\d*\.?\d+"i, value) ?
                        parse(Float64, match(r"gain ?= ?(\d*\.?\d+)"i)) : 200))
                if occursin(r"ebc"i, value)
                    permute!(pins, [3, 2, 1])
                elseif occursin(r"bce"i, value)
                    permute!(pins, [2, 1, 3])
                elseif occursin(r"bec"i, value)
                    permute!(pins, [2, 3, 1])
                elseif occursin(r"ceb"i, value)
                    permute!(pins, [1, 3, 2])
                end
            end
            addtransistor!(circ, connections, tr, points, pins)
        end
    end
    # connect the wires
    for w in wires
        if haskey(connections, w.startpoint) && haskey(connections, w.endpoint)
            connect!(circ, connections[w.startpoint], connections[w.endpoint])
        elseif haskey(connections, w.startpoint) && !haskey(connections, w.endpoint)
            addconnection!(circ, connections, w.endpoint, connections[w.startpoint])
        elseif haskey(connections, w.endpoint) && !haskey(connections, w.startpoint)
            addconnection!(circ, connections, w.startpoint, connections[w.endpoint])
        else
            # handle wires connected to boards or other wires here
        end
    end

    for board in boards
        connectcopperstrips(circ, connections, board, tracecuts)
    end

    free(file)
    return circ
end

end # module AudioToaster

max_stars_repo_name: stellartux/diyloader

max_issues_repo_name: stellartux/diyloader

max_forks_repo_name: stellartux/diyloader


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#maximum-speed.jl

#let's test and see what is the best way of finding the maximum of two numbers?

const count = 1000000

#enclose in a function body because in julia there are no constant type globals.
function testy(T::Type)
println("running for $T")

#generate a series of arrays
test1 = [rand(T) for i=1:count]
test2 = [rand(T) for i=1:count]
test3 = [zero(T) for i=1:count]

res = 0

println("using max()")
tic()
for (i=1:count)
  res = max(test1[i], test2[i])
end
toc()
println(res)

println("using max() against zero, uint")
tic()
for (i=1:count)
  res = max(test1[i], test3[i])
end
toc()
println(res)

println("using min() against zero, uint")
tic()
for (i=1:count)
  res = min(test1[i], test3[i])
end
toc()
println(res)

println("using ternary op")
tic()
for (i=1:count)
  res = (test1[i] > test2[i]) ? test1[i] : test2[i]
end
toc()
println(res)

println("using if...then")
tic()
for (i=1:count)
  if (test1[i] > test2[i])
    res = test1[i]
  else
    res = test2[i]
  end
end
toc()
println(res)

println("using comparison against zero")
tic()
for (i=1:count)
  res = (test2[i] == 0) ? test1[i] : test2[i]
end
toc()
println(res)

println("using comparison not zero")
tic()
for (i=1:count)
  res = (test1[i] != 0) ? test1[i] : test2[i]
end
toc()
println(res)

println("----")
end

testy(UInt64)
testy(Int64)
testy(UInt64)
testy(Int64)
testy(UInt64)
testy(Int64)

#20 Aug 2015
#results, last section only, println(s) edited out.

#running for UInt64
#using max()
#elapsed time: 0.015386148 seconds
#using max() against zero, uint
#elapsed time: 0.036181472 seconds
#using min() against zero, uint
#elapsed time: 0.004971114 seconds
#using ternary op
#elapsed time: 0.021031701 seconds
#using if...then
#elapsed time: 0.021056531 seconds
#using comparison against zero
#elapsed time: 0.040860501 seconds
#using comparison not zero
#elapsed time: 0.014127978 seconds

#running for Int64
#elapsed time: 0.00130995 seconds
#using max() against zero, uint
#elapsed time: 0.001307646 seconds
#using min() against zero, uint
#elapsed time: 0.001316235 seconds
#using ternary op
#elapsed time: 0.00957752 seconds
#using if...then
#elapsed time: 0.008764473 seconds
#using comparison against zero
#elapsed time: 0.001603041 seconds
#using comparison not zero
#elapsed time: 0.001535636 seconds

#results.  Int64 max() and min() functions are MUCH faster.  if/then vs. ternary
#has no speed advantage.  Almost certainly the JNZ() directive is making a
#difference, but it seems alwasy better to just use max() or min().  Strangely
#terary op is worse for unsigned ints than it is for signed ints.

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
module Transit
  import JSON
  import Base.getindex
  import Base

  using DataStructures
  using URIParser
  using Decimals

  export Encoder, encode, parse, decode

  include("tagged_value.jl")
  include("tsymbol.jl")
  include("tset.jl")
  include("turi.jl")
  include("link.jl")
  include("constants.jl")
  include("cache.jl")
  include("decoder.jl")
  include("utilities.jl")
  include("emitter.jl")
  include("encoder.jl")
  include("writer.jl")
  include("reader.jl")
end


text:
function [Grid] = HexagonGridNeu(nx,ny,r,OrientFace)
dphi=2*pi/6;
dx=sqrt(3)*r;
dy=2*r;
x0=sqrt(3)/2*r;
y0=r;

NumNodes=2*(2*nx+1);
Nodes(1:NumNodes)=Node([0 0 0],0);
NodeNumber=1;
yShift=y0;
for j=1:ny
  phi=0;
  xShift=x0;
  for i=1:nx
    phi=phi-dphi;
    x=sin(phi)+xShift;
    y=-cos(phi)+yShift;
    z=0;
    Nodes(NodeNumber)=Node([x;y;z],NodeNumber);
    NodeNumber=NodeNumber+1;
    phi=phi+dphi;
    x=sin(phi)+xShift;
    y=-cos(phi)+yShift;
    z=0;
    Nodes(NodeNumber)=Node([x;y;z],NodeNumber);
    NodeNumber=NodeNumber+1;
    xShift=xShift+dx;
  end
  phi=phi-dphi;
  x=sin(phi)+xShift;
  y=-cos(phi)+yShift;
  z=0;
  Nodes(NodeNumber)=Node([x;y;z],NodeNumber);
  NodeNumber=NodeNumber+1;

  xShift=x0;
  phi=pi;
  yShift=yShift;
  for i=1:nx
    phi=phi+dphi;
    x=sin(phi)+xShift;
    y=-cos(phi)+yShift;
    z=0;
    Nodes(NodeNumber)=Node([x;y;z],NodeNumber);
    NodeNumber=NodeNumber+1;
    phi=phi-dphi;
    x=sin(phi)+xShift;
    y=-cos(phi)+yShift;
    z=0;
    Nodes(NodeNumber)=Node([x;y;z],NodeNumber);
    NodeNumber=NodeNumber+1;
    xShift=xShift+dx;
  end
  phi=phi+dphi;
    x=sin(phi)+xShift;
    y=-cos(phi)+yShift;
    z=0;
    Nodes(NodeNumber)=Node([x;y;z],NodeNumber);
    NodeNumber=NodeNumber+1;
end

Grid.Nodes=Nodes;
Grid.NumNodes=NumNodes;

NumEdges=2*2*nx+nx+1;
Edges(1:NumEdges)=Edge([1 2],Grid,0,0);
EdgeNumber=1;

N1=1;
N2=N1+1;
for i=1:2*nx
  Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber);
  EdgeNumber=EdgeNumber+1;
  N1=N2;
  N2=N2+1;
end
N1=2*nx+2;
N2=N1+1;
for j=1:ny
  for i=1:2*nx
    Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber);
    EdgeNumber=EdgeNumber+1;
    N1=N2;
    N2=N2+1;
  end
end
for j=1:ny
  N1=1;
  N2=N1+2*nx+1;
  for i=1:nx+1
    Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber);
    EdgeNumber=EdgeNumber+1;
    N1=N1+2;
    N2=N2+2;
  end
end
Grid.Edges=Edges;
Grid.NumEdges=NumEdges;
Grid.NumEdgesI=NumEdges;
Grid.NumEdgesB=0;

NumFaces=nx;
Faces(1:NumFaces)=Face([0 0 0 0],Grid,0,'x');
FaceNumber=1;
Type='o';
E1=1;
E2=E1+1;
E3=2*2*nx+2;
E4=E2+2*nx;
E5=E4-1;
E6=E3-1;
for i=1:nx
  Faces(FaceNumber)=Face([E1 E2 E3 E4 E5 E6],Grid,FaceNumber,Type,OrientFace);
  E1=E1+2;
  E2=E2+2;
  E3=E3+1;
  E4=E4+2;
  E5=E5+2;
  E6=E6+1;
  FaceNumber=FaceNumber+1;
end
Grid.Faces=Faces;
Grid.NumFaces=NumFaces;

Grid=FacesInEdges(Grid);
Grid=FacesInNodes(Grid);

end


text:
module Lasso

    using DataFrames
    using Lasso
    using LinearAlgebra
    using StatsModels
    using ..DataMod, ..ManualModelMod, ..Utils

    export lasso, sparse_encoded_lasso
    export MinAIC, MinAICc, MinBIC, MinCVmse, MinCV1se, AllSeg

    function getonemodel(y::Symbol, xs::Vector{Symbol}, df::DataFrame, λs::Vector{Float64}, segselect::Type{<:SegSelect})::RegressionModel
        fit(LassoModel, compose(y, xs, intercept = false), df, λ = λs, standardize = false, cd_maxiter = typemax(Int), select = segselect())
    end
    function getfullpath(y::Symbol, xs::Vector{Symbol}, df::DataFrame, λs::Vector{Float64})
        fit(LassoPath, compose(y, xs, intercept = false), df, λ = λs, standardize = false, cd_maxiter = typemax(Int))
    end

    function lasso(data::Data, λs::Vector{Float64} = .5 .^ collect(-5:14);
                   fullpath::Bool = false, xs::Vector{Symbol} = data.xs,
                   segselect::Type{<:SegSelect} = MinBIC)
        if fullpath
            return RegressionModel[getfullpath(y, xs, data.df, λs) for y in data.ys]
        else
            return RegressionModel[getonemodel(y, xs, data.df, λs, segselect) for y in data.ys]
        end
    end

    function sparse_encoded_lasso(data::Data, λs::Vector{Float64} = .5 .^ collect(-5:14))::Vector{ManualModel}
        models = lasso(data, λs)
        deselected = [coef(model) .== 0 for model in models]
        len = eachindex(data.ys)
        nobs = nrow(data.df)
        dat = [hcat(ones(Float64, nobs), data.xmat) for _ in len]
        empt = zeros(nobs)
        @inbounds for i in len
            dat[i][:, deselected[i]] .= empt
        end
        mp = [pinv(mat) for mat in dat]
        yvec = [data.df[!, y] for y in data.ys]
        wts = mp .* yvec
        return [ManualModel(wts[i], true, data.ys[i], data) for i in len]::Vector{ManualModel}
    end

end

using .Lasso
export lasso, sparse_encoded_lasso
export MinAIC, MinAICc, MinBIC, MinCVmse, MinCV1se, AllSeg

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IPRule <: SwaggerModel
    value::Any # spec type: Union{ Nothing, String } # spec name: value
    action::Any # spec type: Union{ Nothing, String } # spec name: action

    function IPRule(;value=nothing, action="Allow")
        o = new()
        validate_property(IPRule, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(IPRule, Symbol("action"), action)
        setfield!(o, Symbol("action"), action)
        o
    end
end # type IPRule

const _property_map_IPRule = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("action")=>Symbol("action"))
const _property_types_IPRule = Dict{Symbol,String}(Symbol("value")=>"String", Symbol("action")=>"String")
Base.propertynames(::Type{ IPRule }) = collect(keys(_property_map_IPRule))
Swagger.property_type(::Type{ IPRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IPRule[name]))}
Swagger.field_name(::Type{ IPRule }, property_name::Symbol) =  _property_map_IPRule[property_name]

const _allowed_IPRule_action = ["Allow"]

function check_required(o::IPRule)
    (getproperty(o, Symbol("value")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IPRule }, name::Symbol, val)
    if name === Symbol("action")
        Swagger.validate_param(name, "IPRule", :enum, val, _allowed_IPRule_action)
    end
end


text:
@testset "Base" begin

@testset "ones/zeros" begin
    @test AMDGPU.ones(4, 3) isa ROCArray
    @test AMDGPU.zeros(3, 4) isa ROCArray
end

@testset "view" begin
    xs = copyto!(ROCVector{Int}(undef, 4), 1, collect(1:4), 1, 4)
    a = view(xs, 1:2)
    b = view(xs, 3:4)
    @test a isa ROCVector{Int}
    @test b isa ROCVector{Int}
    @test collect(a)::Vector{Int} == 1:2
    @test collect(b)::Vector{Int} == 3:4
    @allowscalar begin
        @test a[[1, 2]] == 1:2
        @test b[[1, 2]] == 3:4
    end

    b_cpu = Vector{Int}(undef, 2)
    @test copyto!(b_cpu, 1, b, 1, 2) === b_cpu
    @test b_cpu == 3:4

    c = ROCVector{Int}(undef, 2)
    @test copyto!(c, 1, b, 1, 2) === c
    @test collect(c) == 3:4
end

@testset "unsafe_wrap" begin
    A = rand(4, 3)
    A_orig = copy(A)
    RA = Base.unsafe_wrap(ROCArray, pointer(A), size(A))
    @test RA.buf.agent == get_default_agent()
    @test RA isa ROCArray{Float64,2}

    # GPU pointer works
    RA .+= 1.0

    # Host pointer is updated
    @test A ≈ A_orig .+ 1.0

    # Base.show
    @test (println(RA); true)

    # Mem.download!
    B = zeros(4, 3)
    copyto!(B, RA)
    @test B ≈ Array(RA)

    # Mem.upload!
    C = rand(4, 3)
    copyto!(RA, C)
    @test Array(RA) ≈ C

    # Mem.transfer!
    D = rand(4, 3)
    D_orig = copy(D)
    RD = Base.unsafe_wrap(ROCArray, pointer(D), size(D))
    copyto!(RD, RA)
    @test Array(RD) ≈ Array(RA) ≈ C
end

end