text:
using BQCESubroutine
using BQCESubroutine.Utilities
using LinearAlgebra
using YaoLocations: YaoLocations, plain, Locations, CtrlLocations, merge_locations
using Test
using BenchmarkTools
using YaoArrayRegister
using LoopVectorization
using PaddedMatrices
using MatchCore
using MLStyle
using BQCESubroutine: BitContext, BitRoutine, multi_broutine2x2!
using BQCESubroutine: diag_kernel, perm_kernel, generic_kernel, threaded_subspace_loop
using BQCESubroutine: codegen_broutine2x2_multi, codegen_broutine, hoist_entries!, codegen_basic_broutine
using BQCESubroutine: threaded_basic_broutine!, basic_broutine!, threaded_multi_broutine2x2!

st = rand(ComplexF64, 1 << 10);
locs = Locations((9, 10))
threaded_basic_broutine!(copy(st), Val(:PSWAP), locs, 2.0) ≈ basic_broutine!(copy(st), Val(:PSWAP), locs, 2.0)
broutine!(copy(st), Val(:Z), locs) ≈ broutine!(copy(st), BQCESubroutine.BZ, locs)
threaded_basic_broutine!(copy(st), Val(:Z), locs) ≈ broutine!(copy(st), BQCESubroutine.BZ, locs)
threaded_basic_broutine!(copy(st), Val(:Z), locs) ≈ basic_broutine!(copy(st), Val(:Z), locs)

threaded_multi_broutine2x2!(copy(st), Val(:Z), locs) ≈ broutine!(copy(st), BQCESubroutine.BZ, locs)
threaded_multi_broutine2x2!(copy(st), Val(:Z), locs) ≈ multi_broutine2x2!(copy(st), Val(:Z), locs)

@code_lowered broutine!(st, Val(:Z), Locations(24))
BQCESubroutine.set_nthreads(12)
@benchmark broutine!(st, Val(:Z), Locations(24)) setup=(st=rand(ComplexF64, 1 << 25))
@benchmark threaded_basic_broutine!(st, Val(:Z), Locations(25)) setup=(st=rand(ComplexF64, 1 << 25))
@benchmark basic_broutine!(st, Val(:Z), Locations(25)) setup=(st=rand(ComplexF64, 1 << 25))

BQCESubroutine.disable_threads()
@benchmark broutine!(st, Val(:Z), Locations(25)) setup=(st=rand(ComplexF64, 1 << 25))
@benchmark instruct!(st, Val(:Z), (25, )) setup=(st=rand(ComplexF64, 1 << 25))

@benchmark threaded_basic_broutine!(st, Val(:Z), Locations(25)) setup=(st=rand(ComplexF64, 1 << 25))

@benchmark threaded_multi_broutine2x2!(st, Val(:Z), Locations((23, 24))) setup=(st=rand(ComplexF64, 1 << 25))
@benchmark multi_broutine2x2!(st, Val(:Z), Locations((23, 24))) setup=(st=rand(ComplexF64, 1 << 25))
@benchmark instruct!(st, Val(:Z), (23, 24)) setup=(st=rand(ComplexF64, 1 << 25))

@broutine X = [0 1;1 0]

threaded_basic_broutine!(copy(st), Val(:X), locs) ≈ broutine!(copy(st), BX, locs)
multi_broutine2x2!(copy(st), Val(:X), locs)≈ broutine!(copy(st), BX, locs)
threaded_multi_broutine2x2!(copy(st), Val(:X), locs) ≈ multi_broutine2x2!(copy(st), Val(:X), locs)

ctx = BitContext(;threading=true, expand_sz = 0)
brt = BitRoutine(Main, :(Z = [1 0;0 -1]))
codegen_basic_broutine(ctx, brt)|>flatten_block

function BQCESubroutine.threaded_basic_broutine!(var"##st#297"::AbstractVector{T} where T, ::Val{:Z}, var"##locs#298"::Locations)
    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:162 =#
    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:258 =#
    nothing
    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:259 =#
    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:296 =#
    var"##plain_locs#299" = (YaoLocations.plain)(var"##locs#298")
    var"##nqubits#300" = (BQCESubroutine.Utilities.log2i)(size(var"##st#297", 1))
    var"##nlocs_needed#301" = (BQCESubroutine.Utilities.log2i)(BQCESubroutine.BQCES_nthreads[] - 1) + 1
    var"##step_1_l#305" = 1 << (var"##plain_locs#299"[1] - 1)
    var"##step_1_h#306" = 1 << var"##plain_locs#299"[1]
    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:297 =#
    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:298 =#
    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:561 =# @fastmath #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:560 =# @inbounds(begin
                if var"##nlocs_needed#301" ≤ var"##nqubits#300" - var"##plain_locs#299"[1]
                    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:455 =#
                    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:414 =#
                    if var"##step_1_l#305" == 1
                        #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:415 =#
                        Threads.@threads for var"##idx#302_1" = 0:var"##step_1_h#306":size(var"##st#297", 1) - var"##step_1_h#306"
                                #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:312 =#
                                nothing
                                #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:313 =#
                                var"##I#307_2" = ((var"##idx#302_1" + 0) + var"##step_1_l#305") + 1
                                var"##T#308_2" = -(var"##st#297"[var"##I#307_2"])
                                var"##st#297"[var"##I#307_2"] = var"##T#308_2"
                            end
                        #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:416 =#
                        return var"##st#297"
                    end
                    Threads.@threads for var"##idx#302_1" = 0:var"##step_1_h#306":size(var"##st#297", 1) - var"##step_1_h#306"
                            for var"##idx#302_2" = var"##idx#302_1":1:(var"##idx#302_1" + var"##step_1_l#305") - 1
                                #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:426 =#
                                #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:312 =#
                                nothing
                                #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:313 =#
                                var"##I#309_2" = ((var"##idx#302_2" + 0) + var"##step_1_l#305") + 1
                                var"##T#310_2" = -(var"##st#297"[var"##I#309_2"])
                                var"##st#297"[var"##I#309_2"] = var"##T#310_2"
                            end
                        end
                    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:456 =#
                    return var"##st#297"
                end
                #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:478 =#
                var"##m#304" = (var"##nqubits#300" - var"##nlocs_needed#301") - 1
                #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:479 =#
                Threads.@threads for var"##base#303" = (BQCESubroutine.Utilities.bsubspace)(var"##nqubits#300", (1:var"##m#304"..., var"##plain_locs#299"...))
                        #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:480 =#
                        for var"##idx#302_1" = 0:1 << var"##m#304" - 1
                            #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:481 =#
                            #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:312 =#
                            nothing
                            #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:313 =#
                            var"##I#311_2" = ((var"##base#303" + var"##idx#302_1") + var"##step_1_l#305") + 1
                            var"##T#312_2" = -(var"##st#297"[var"##I#311_2"])
                            var"##st#297"[var"##I#311_2"] = var"##T#312_2"
                        end
                    end
                return var"##st#297"
            end)
    #= /home/roger/julia_code/BQCESubroutine/src/codegen/broutine.jl:163 =#
    return var"##st#297"
end


text:
using LinearAlgebra

function interpolateCellsToFaces(field, method)
    if method == "harmonicMean"
        for iFace in 1:mesh.nInternalFaces
            face = mesh.faces[iFace]

            ownerValue = field.cellValues[face.iOwner]
            neighbourValue = field.cellValues[face.iNeighbour]

            field.faceValues[iFace] = 2 * ownerValue * neighbourValue / (ownerValue + neighbourValue)
        end
    elseif method == "linear"
        for iFace in 1:mesh.nInternalFaces
            face = mesh.faces[iFace]

            ownerValue = field.cellValues[face.iOwner]
            neighbourValue = field.cellValues[face.iNeighbour]

            field.faceValues[iFace] = (1 - face.gf) * ownerValue + face.gf * neighbourValue
        end
    end
end

ext: jl


text:
#! format: off
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    # pretty
    Base.precompile(Tuple{typeof(pretty),DefaultStyle,CSTParser.EXPR,State})
    Base.precompile(Tuple{typeof(pretty),YASStyle,CSTParser.EXPR,State})
    Base.precompile(Tuple{typeof(pretty),BlueStyle,CSTParser.EXPR,State})

end


text:
# Silverman's rule of thumb for KDE bandwidth selection
function bandwidth(data::Vector{Float64}, alpha::Float64 = 0.9)
    # Determine length of data
    ndata = length(data)

    # Calculate width using variance and IQR
    var_width = std(data)
    q25, q75 = quantile(data, [0.25, 0.75])
    quantile_width = (q75 - q25) / 1.34

    # Deal with edge cases with 0 IQR or variance
    width = min(var_width, quantile_width)
    if width == 0.0
        if var_width == 0.0
            width = 1.0
        else
            width = var_width
        end
    end

    # Set bandwidth using Silverman's rule of thumb
    return alpha * width * ndata^(-0.2)
end


text:
if basename(pwd()) == "aoc"
    cd("2016/22")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
##### problem_stats.jl
# Created by Alexander Dowling

############################
# ProblemStats

type ProblemStats
	nCont::Int64
	nContFixed::Int64
	nContFree::Int64
	nContLowerOnly::Int64
	nContUpperOnly::Int64
	nContTwoBounds::Int64
	nContUnbounded::Int64
	nContActiveBounds::Int64
	nBin::Int64
	nBinFixed::Int64
	nBinFree::Int64

	linear::Int64
	linEql::Int64
	linIneqlAct::Int64
	linIneqlInact::Int64
	quad::Int64
	quadEql::Int64
	quadIneqlAct::Int64
	quadIneqlInact::Int64
	nonlinear::Int64
	nonlinEql::Int64
	nonlinIneqlAct::Int64
	nonlinIneqlInact::Int64

	status::Symbol
	solveTime::Float64
end

function ProblemStats(m::Model, status::Symbol=:Unsolved, solveTime::Float64=0.0)
	return assembleProblemStats(m, status, solveTime)
end

function assembleProblemStats(m::Model, status=:Unsolved, solveTime=0.0)

	epsilon = 1E-6

	fixed = m.colUpper - m.colLower .< epsilon
	cont = m.colCat .== :Cont
	lwr = m.colLower .> -Inf
	upr = m.colUpper .< Inf
	bnry = m.colCat .== :Bin

	n = MathProgBase.numconstr(m)

	nLin = MathProgBase.numlinconstr(m)
	nQuad = MathProgBase.numquadconstr(m)
	nNonLin = n - nLin - nQuad

	(gLB, gUB) = JuMP.constraintbounds(m)

	# This portion of code converts m into an NLP
	d = JuMP.NLPEvaluator(m)
	MathProgBase.initialize(d, [:ExprGraph, :Jac])

	x = get_x(m, status)

	g = zeros(n,1)
	MathProgBase.eval_g(d,g,x)

	eql = (gUB - gLB) .< epsilon
	ineql = .!eql
	active = ((gUB - g) .< epsilon) .| ((g - gLB) .< epsilon)

	N = 1:n

	lin = N .<= nLin
	quad = (N .> nLin) .& (N .<= nLin + nQuad)
	nonlin = N .> nLin + nQuad

	activeBnds = ((x - m.colLower) .< epsilon) .| ((m.colUpper - x) .< epsilon)

	return ProblemStats( sum(cont), # Continuous
						sum(cont .& fixed), # Fixed
						sum(cont .& .!fixed), # Free
						sum(cont .& lwr .& .!upr), # Only lower bounded
						sum(cont .& .!lwr .& upr), # Only upper bounded
						sum(cont .& lwr .& upr .& .!fixed), # Lower and upper bounded
						sum(cont .& .!lwr .& .!upr), # Unbounded
						sum(activeBnds .& .!fixed), # Number of active bounds
						sum(bnry), # Binary
						sum(bnry .& fixed), # Binary fixed
						sum(bnry .& .!fixed), # Binary free
						nLin, # Linear constraints
						sum(lin .& eql),
						sum(lin .& ineql .& active),
						sum(lin .& ineql .& .!active),
						nQuad, # Quadratic constraints
						sum(quad .& eql),
						sum(quad .& ineql .& active),
						sum(quad .& ineql .& .!active),
						nNonLin, # Nonlinear constraints
						sum(nonlin .& eql),
						sum(nonlin .& ineql .& active),
						sum(nonlin .& ineql .& .!active),
						status, # Solver status
						solveTime # Solver time
						)

end

function stringProblemStats(ps::ProblemStats)

	return join([ps.nCont, ps.nContFixed, ps.nContFree, ps.nContLowerOnly, ps.nContUpperOnly,ps.nContTwoBounds,ps.nContUnbounded,ps.nContActiveBounds,
					ps.nBin,ps.nBinFixed,ps.nBinFree,
					ps.linear,ps.linEql,ps.linIneqlAct,ps.linIneqlInact,
					ps.quad,ps.quadEql,ps.quadIneqlAct,ps.quadIneqlInact,
					ps.nonlinear,ps.nonlinEql,ps.nonlinIneqlAct,ps.nonlinIneqlInact,
					string(ps.status),ps.solveTime],",")

end

function stringProblemStats()

	return string("nContinous,nContFixed,nContFree,nContLowrOnly,nContUpprOnly,nContBndBoth,nContUnbnd,nContActiveBnds,nBinary,nBinFixed,nBinFree,nLinear,nLinEql,nLinIneqlAct,nLinIneqlInact,nQuadratic,nQuadEql,nQuadIneqlAct,nQuadIneqlInact,nNonlinear,nNLEql,nNLIneqlAct,nNLIneqlInact,status,solveTime")

end

function printProblemStats(ps::ProblemStats, f=STDOUT)

	println(f," ")
	println(f,"Continous Variables (total): ",ps.nCont)
	println(f,"\t Fixed: ",ps.nContFixed)
	println(f,"\t Free: ",ps.nContFree)
	println(f,"\t\t Lower Bounded Only: ",ps.nContLowerOnly)
	println(f,"\t\t Upper Bounded Only: ",ps.nContUpperOnly)
	println(f,"\t\t Lower \& Upper Bounded: ",ps.nContTwoBounds)
	println(f,"\t\t Unbounded: ",ps.nContUnbounded)
	println(f,"\t\t Active Bounds: ",ps.nContActiveBounds)
	println(f," ")
	println(f,"Binary Variable (total): ",ps.nBin)
	println(f,"\t Fixed: ",ps.nBinFixed)
	println(f,"\t Free: ",ps.nBinFree)
	println(f," ")

	println(f,"Linear Constraints: ",ps.linear)
	println(f,"\t Equality: ",ps.linEql)
	println(f,"\t Ineq. (Active): ",ps.linIneqlAct)
	println(f,"\t Ineq. (Inactive): ",ps.linIneqlInact)
	println(f," ")

	println(f,"Quadratic Constraints: ",ps.quad)
	println(f,"\t Equality: ",ps.quadEql)
	println(f,"\t Ineq. (Active): ",ps.quadIneqlAct)
	println(f,"\t Ineq. (Inactive): ",ps.quadIneqlInact)
	println(f," ")

	println(f,"Nonlinear Constraints: ",ps.nonlinear)
	println(f,"\t Equality: ",ps.nonlinEql)
	println(f,"\t Ineq. (Active): ",ps.nonlinIneqlAct)
	println(f,"\t Ineq. (Inactive): ",ps.nonlinIneqlInact)
	println(f," ")

	println(f,"Solver Exit Status: ",string(ps.status))
	println(f,"Solve Time: ",string(ps.solveTime), "s")
	println(f," ")

	return nothing

end


text:
@testset "Dictionary" begin
    @test Dictionary() isa Dictionary{Any, Any}

    d = Dictionary{Int64, Int64}()

    @test eltype(d) === Int64
    @test keytype(d) === Int64
    @test tokentype(d) === Int
    @test Base.IteratorSize(d) === Base.HasLength()
    @test isinsertable(d)
    @test length(d) == 0
    @test isempty(d)
    @test isempty(keys(d))
    @test d == d
    @test d == copy(d)
    @test d == Dictionary(copy(keys(d)), d)
    @test isequal(d, d)
    @test isequal(copy(d), d)
    @test isequal(Dictionary(copy(keys(d)), d), d)
    @test !isless(d, d)
    @test !isless(copy(d), d)
    @test !isless(Dictionary(copy(keys(d)), d), d)
    @test cmp(d, d) == 0
    @test cmp(copy(d), d) == 0
    @test cmp(Dictionary(copy(keys(d)), d), d) == 0
    @test_throws IndexError d[10]
    @test get(d, 10, 15) == 15
    @test length(unset!(d, 10)) == 0
    io = IOBuffer(); print(io, d); @test String(take!(io)) == "{}"
    if VERSION < v"1.6-"
        io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "0-element Dictionary{Int64,Int64}"
    else
        io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "0-element Dictionary{Int64, Int64}"
    end
    @test_throws IndexError d[10] = 11
    @test_throws IndexError delete!(d, 10)

    insert!(d, 10.0, 11.0)

    @test d[10] == 11
    @test haskey(d, 10)
    @test !haskey(d, 100)
    @test get(d, 10, 15) == 11
    @test get!(d, 10, 15) == 11
    @test length(d) == 1
    @test keys(d) === keys(keys(d))
    @test unique(d)::Indices{Int64} == Indices([11])
    @test !isempty(d)
    @test d == d
    @test d == copy(d)
    @test d == Dictionary(copy(keys(d)), d)
    @test d != empty(d)
    @test empty(d) != d
    @test fill(0, d) != d
    @test d != fill(0, d)
    @test fill(0, copy(keys(d))) != d
    @test d != fill(0, copy(keys(d)))
    @test isequal(d, d)
    @test isequal(copy(d), d)
    @test isequal(Dictionary(copy(keys(d)), d), d)
    @test !isequal(d, empty(d))
    @test !isequal(empty(d), d)
    @test !isequal(fill(0, d), d)
    @test !isequal(d, fill(0, d))
    @test !isequal(fill(0, copy(keys(d))), d)
    @test !isequal(d, fill(0, copy(keys(d))))
    @test !isless(d, d)
    @test !isless(copy(d), d)
    @test !isless(Dictionary(copy(keys(d)), d), d)
    @test !isless(d, empty(d))
    @test isless(empty(d), d)
    @test isless(fill(0, d), d)
    @test !isless(d, fill(0, d))
    @test isless(fill(0, copy(keys(d))), d)
    @test !isless(d, fill(0, copy(keys(d))))
    @test cmp(d, d) == 0
    @test cmp(copy(d), d) == 0
    @test cmp(Dictionary(copy(keys(d)), d), d) == 0
    @test cmp(d, empty(d)) == 1
    @test cmp(empty(d), d) == -1
    @test cmp(fill(0, d), d) == -1
    @test cmp(d, fill(0, d)) == 1
    @test cmp(fill(0, copy(keys(d))), d) == -1
    @test cmp(d, fill(0, copy(keys(d)))) == 1
    @test_throws IndexError insert!(d, 10, 12)
    @test d[10.0] == 11
    set!(d, 10.0, 12.0)
    @test length(d) == 1
    @test d[10] == 12
    d[10.0] = 13.0
    @test length(d) == 1
    @test d[10] == 13
    io = IOBuffer(); print(io, d); @test String(take!(io)) == "{10 = 13}"
    if VERSION < v"1.6-"
        io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "1-element Dictionary{Int64,Int64}\n 10 │ 13"
    else
        io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "1-element Dictionary{Int64, Int64}\n 10 │ 13"
    end
    @test !isequal(d, empty(d))
    @test isequal(d, copy(d))
    @test isempty(empty(d))

    delete!(d, 10.0)
    @test isequal(d, Dictionary{Int64, Int64}())

    @test get!(d, 10, 14.0) == 14
    @test d[10] == 14
    delete!(d, 10)

    for i = 2:2:1000
        insert!(d, i, i+1)
        @test d[i] == i + 1 ? true : (@show i; false)
    end
    @test all(in(i, d) == !iseven(i) for i in 2:2:1000)
    @test all(in(i, keys(d)) == iseven(i) for i in 2:2:1000)
    @test isempty(empty!(d))

    @test get!(() -> 15, d, 10) == 15
    @test get!(() -> 16, d, 10) == 15

    d = Dictionary([:a, :b], [1, 2])
    @test isequal(d, d)
    @test d == d
    @test !isless(d, d)
    d2 = Dictionary((a=1, b=2))
    @test isequal(d, d2)
    @test d == d
    d3 = dictionary([:a=>1, :b=>2])
    @test isequal(d, d3)
    d4 = dictionary(zip([:a, :b], [1, 2]))
    @test isequal(d, d4)
    @test !isless(d, d4)
    @test !isless(d4, d)
    @test hash(d) == hash(d4)

    @test isdictequal(d, copy(d))
    @test isdictequal(Dictionary(['a','b'],[1,2]), Dictionary(['b','a'],[2,1]))
    @test !isdictequal(Dictionary(['a','b'],[1,2]), Dictionary(['b','c'],[2,1]))
    @test !isdictequal(Dictionary(['a','b'],[1,2]), Dictionary(['a','b','c'],[1,2,3]))
    @test !isdictequal(Dictionary(['a','b'],[1,2]), Dictionary(['b','a'],[2,3]))

    d5 = Dictionary(['a','b'],[1,missing])
    @test isdictequal(d5, d5) === missing
    @test (d5 == d5) === missing
    d6 = Dictionary(['a','b'],[1,missing])
    @test isdictequal(d5, d5) === missing
    @test (d5 == d5) === missing

    @test isequal(merge(d, d), d)
    @test isequal(merge(d, d2), d)

    @test isequal(merge(d, Dictionary([:c], [3])), Dictionary([:a, :b, :c], [1, 2, 3]))
    @test isequal(merge(d, Dictionary([:b, :c], [4, 3])), Dictionary([:a, :b, :c], [1, 4, 3]))

    @test isequal(index(first, ["Alice", "Bob", "Charlie"]), Dictionary(['A', 'B', 'C'], ["Alice", "Bob", "Charlie"]))
    @test isequal(index(first, ["Alice", "Bob", "Charlie", "Conner"]), Dictionary(['A', 'B', 'C'], ["Alice", "Bob", "Charlie"]))

    d7 = Dictionary(Int32[1, 2], UInt32[1, 2])
    @test convert(Dictionary{Int64, Int64}, d7)::Dictionary{Int64, Int64} == Dictionary(Int64[1, 2], Int64[1, 2])

    rd = Dictionary([:b, :a], [2, 1])
    @test reverse(d)::Dictionary == rd
    @test Iterators.reverse(d)::ReverseDictionary == rd
    @test first(d) == 1
    @test last(d) == 2
    @test firstindex(d) == :a
    @test lastindex(d) == :b
    #@test d[begin] == 1 # Parsing issues on earlier versions of Julia...
    @test d[end] == 2

    # TODO token interface

    @testset "filter!" begin
        d = Dictionary([1,2,3,4,5], [1,3,2,4,5])
        filter!(iseven, d)
        @test d == Dictionary([3,4], [2,4])
    end

    @testset "Dict tests from Base" begin
        h = Dictionary{Int, Int}()

        for i in 1:10000
            insert!(h, i, i+1)
        end
        for i in 1:10000
            @test h[i] == i+1
        end
        for i in 1:2:10000
            delete!(h, i)
        end
        for i in 1:10000
            if iseven(i)
                @test h[i] == i+1
            else
                @test_throws IndexError h[i]
            end
        end
        for i in 1:2:10000
            insert!(h, i, i+1)
        end
        for i in 1:10000
            @test h[i] == i+1
        end
        for i in 1:10000
            delete!(h, i)
        end
        @test isempty(h)
        insert!(h, 77, 100)
        @test h[77] == 100
        for i in 1:10000
            set!(h, i, i+1)
        end
        for i in 1:10000
            @test h[i] == i+1
        end
        for i in 1:2:10000
            delete!(h, i)
        end
        for i in 1:10000
            if iseven(i)
                @test h[i] == i+1
            else
                @test_throws IndexError h[i]
            end
        end
        for i in 10001:20000
            insert!(h, i, i+1)
        end
        for i in 1:10000
            if iseven(i)
                @test h[i] == i+1
            else
                @test_throws IndexError h[i]
            end
        end
        for i in 10000:20000
            @test h[i] == i+1
        end
    end

    @testset "dictionary" begin
        res = Dictionary(['a','b','c'], [1,2,3])
        @test isequal(dictionary(pairs(res)), res)
        @test isequal(dictionary(['a'=>1, 'b'=>2, 'c'=>3]), res)
        @test isequal(dictionary(['a'=>1, 'b'=>2, 'c'=>3, 'a'=>4]), res)
        @test isequal(dictionary((k,v) for (k,v) in pairs(res)), res)
    end

    @testset "index" begin
        res = Dictionary(['A','B','C'], ["Alice","Bob","Charlie"])
        @test isequal(index(first, ["Alice", "Bob", "Charlie"]), res)
        @test isequal(index(first, ["Alice", "Bob", "Charlie", "Conner"]), res)
    end

    @testset "Factories" begin
        d = Dictionary(['a','b','c'], [1,2,3])
        @test similar(d) isa Dictionary{Char, Int}
        @test similar(d, Float64) isa Dictionary{Char, Float64}
        @test sharetokens(d, similar(d))

        @test isempty(empty(d)::Dictionary{Char, Int})
        @test isempty(empty(d, Float64)::Indices{Float64})
        @test isempty(empty(d, String, Float64)::Dictionary{String, Float64})

        @test isequal(zeros(d)::Dictionary{Char, Float64}, Dictionary(['a','b','c'],[0.0,0.0,0.0]))
        @test isequal(zeros(Int64, d)::Dictionary{Char, Int64}, Dictionary(['a','b','c'],[0,0,0]))

        @test isequal(ones(d)::Dictionary{Char, Float64}, Dictionary(['a','b','c'],[1.0,1.0,1.0]))
        @test isequal(ones(Int64, d)::Dictionary{Char, Int64}, Dictionary(['a','b','c'],[1,1,1]))

        @test isequal(keys(rand(1:10, d)::Dictionary{Char, Int}), Indices(['a','b','c']))
        @test isequal(keys(randn(d)::Dictionary{Char, Float64}), Indices(['a','b','c']))
    end

    @testset "rand" begin
        dict = Dictionary(["a", "b", "c", "d", "e"], 1:5)
        for i = 1:100
            @test rand(dict) in 1:5
        end

        delete!(dict, "b") # Test handling of deleted slots
        for i = 1:100
            @test rand(dict) in [1, 3, 4, 5]
        end
    end

    @testset "sort" begin
        dict = Dictionary([1, 3, 2], ['c', 'a', 'b'])
        @test sort(dict)::Dictionary == Dictionary([3, 2, 1], ['a', 'b', 'c'])
        @test sort(dict; rev=true)::Dictionary == Dictionary([1, 2, 3], ['c', 'b', 'a'])

        @test sortkeys(dict)::Dictionary == Dictionary([1, 2, 3], ['c', 'b', 'a'])
        @test sortkeys(dict; rev=true)::Dictionary == Dictionary([3, 2, 1], ['a', 'b', 'c'])

        @test sortpairs(dict)::Dictionary == Dictionary([1, 2, 3], ['c', 'b', 'a'])
        @test sortpairs(dict; rev=true)::Dictionary == Dictionary([3, 2, 1], ['a', 'b', 'c'])
        @test sortpairs(dict; by=kv->kv.second=>kv.first)::Dictionary == Dictionary([3, 2, 1], ['a', 'b', 'c'])

        dictcopy = deepcopy(dict)
        sort!(dictcopy)
        @test dictcopy == Dictionary([3, 2, 1], ['a', 'b', 'c'])

        dictcopy = deepcopy(dict)
        sort!(dictcopy; rev = true)
        @test dictcopy == Dictionary([1, 2, 3], ['c', 'b', 'a'])

        dictcopy = deepcopy(dict)
        sortkeys!(dictcopy)
        @test dictcopy == Dictionary([1, 2, 3], ['c', 'b', 'a'])

        dictcopy = deepcopy(dict)
        sortkeys!(dictcopy; rev = true)
        @test dictcopy == Dictionary([3, 2, 1], ['a', 'b', 'c'])

        dictcopy = deepcopy(dict)
        sortpairs!(dictcopy)
        @test dictcopy == Dictionary([1, 2, 3], ['c', 'b', 'a'])

        dictcopy = deepcopy(dict)
        sortpairs!(dictcopy; rev = true)
        @test dictcopy == Dictionary([3, 2, 1], ['a', 'b', 'c'])

        dictcopy = deepcopy(dict)
        sortpairs!(dictcopy; by = kv->kv.second=>kv.first)
        @test dictcopy == Dictionary([3, 2, 1], ['a', 'b', 'c'])
    end
end


text:
using MathGL

# Create the graph, using a sufficient resolution for 4x2 pictures here
gr = MathGL.Graph(2400, 1200)

# Standard sine and cosine plots with grid ect.
@mglplot gr [
    # Give position of the next plot: In a grid of 4x2 plots take the first (top left)
    subplot 4 2 1
    # The title for this subplot
    title "Sine and cosine plot"
    # The xticks, given as array of values and array of corresponding strings
    xtick [-2pi, -pi, 0, pi, 2pi] ["-2\\pi", "-\\pi", "0", "\\pi", "2\\pi"]
    # Set the ranges for x and y respectively. Each range is -1 to 1 by default
    xrange -2pi 2pi; yrange -1 1
    # Plot the axis (meaning: plot tics and tic labels) and the grid with style "h="
    # 'h' is a lighter gray, '=' is a type of dashed line
    axis; grid stl="h="
    # Labels for the x and y axis
    xlabel "x"; ylabel "y"
    # The x-values (as array, not as linspace!) that shall be used for plotting
    x = Array(linspace(-2pi,2pi,100))
    # The values for the two curves
    y1 = 0.5sin(x); y2 = 0.5cos(x)
    # Plot with style options "B" (dark Blue) and "G" (dark green)
    plot x y1 stl="B"
    plot x y2 stl="G"
    # Reset the tics
    xtick
]

# Some log-log plot
@mglplot gr [
    # Same as before
    subplot 4 2 2
    title "Log-log plot"
    # This time: set x and y range in one command, ranges, instead of
    # xrange and yrange
    ranges 1e-2 1e2 7e-1 1e1
    # Set logscale axis by defining the coordinate transformations for the
    # x respectively y values. So we have  φ: xrange x yrange --> ℝ²,
    # φ(x,y) = (lg(x), lg(y)) as transformation applied before plotting
    coords X="lg(x)" Y="lg(y)"
    # Plot the ticks and tick labels, plot the grid for minor ticks
    # (subgrid) and the ordinary grid
    axis; subgrid; grid
    # Once again the labels. MathGL supports some light form of support for
    # latex-like commands
    xlabel "x"; ylabel "y = \\sqrt{1+x^2}"
    # Generate the data and plot
    x = 1logspace(-2, 2, 200)
    plot x sqrt(1+x)
]

# Using cylindrical coordinates
@mglplot gr [
    subplot 4 2 3
    title "Cylindrical plot"
    # This range will be fed to the transformation of cylindrical coordinates
    # The xrange has the role of the angle / times 2pi, so it is between 0 and 1
    # The yrange has the role of the radius, it is also between 0 and 1
    # The zrange stays the same, so it stays from -1 to 1
    ranges 0 1 0 1 -1 1
    coords X="y*sin(2*pi*x)" Y="y*cos(2*pi*x)"
    # Rotate the picture. 50 degrees away from top view, -60 degrees around
    # z axis
    rotate 50 -60
    # Plot 6 circles
    x = Array(linspace(0,1,50))
    o = ones(50)
    # The numbers in the style string determine the width of the line to be
    # drawn
    plot x 0.80o 0.7o  stl="b9"
    plot x 0.65o 0.3o  stl="B7"
    plot x 0.50o 0.0o  stl="r5"
    plot x 0.30o -0.3o stl="G3"
    plot x 0.10o -0.7o stl="g1"
    # Set the origin to a margin value so that it does not disturb drawing the grid
    origin -1 1 -1
    # Plot a grid with light grey color ('h') and undashed lines ('-')
    grid stl="h-"
    # reset coordinates
    coords
]

# A plot marking the area between the y values and the x-axis
@mglplot gr [
    subplot 4 2 4
    title "Area plot"
    # Set ticks and ranges like before
    xtick [-2pi, -pi, 0, pi, 2pi] ["-2\\pi", "-\\pi", "0", "\\pi", "2\\pi"]
    xrange -2pi 2pi; yrange -1 1; zrange 0 0
    # As an alternative: set the origin to the top of the drawing. This
    # affects the axis command and the xlabel / ylabel commands so that the
    # ticks ect. for the x-axis are now drawn above the graph
    origin -2pi 1
    box; axis
    xlabel "x"; ylabel "y"
    # Set the origin back to 0. This is important for the area-plot, since
    # the area marked is the area between the x-axis defined by the origin
    origin 0 0
    # Prepare data
    x = Array(linspace(-2pi,2pi,100))
    y1 = 0.5sin(x)+0.4sin(3x); y2 = 0.8cos(x); y3 = -0.3cos(x)
    # Do the plots. g,G: light,dark green; k: black; y: yellow, c: cyan, b:
    # blue, #: special style for the area drawing: draw lines instead of
    # filling the full area
    area x y1 stl="gG#"
    plot x y1 stl="k"
    area x y2 stl="yr#"
    plot x y2 stl="k"
    area x y3 stl="cb#"
    plot x y3 stl="k"
    # reset the xticks
    xtick
]

# A default chart plot (as implemented in MathGL)
@mglplot gr [
    subplot 4 2 5
    title "Chart plot"
    # Use light-effects. One can also set up different light sources of different color and type
    light true
    rotate 50 60
    # Prepare data by constructing a random 7x2 matrix. So we compare
    # 2 datasets with 7 values each
    data = rand(7,2)+0.1
    # Here '#' means that the edges shall be plotted as lines (black by default)
    chart data stl="#"
]

# A ring like chart plot
@mglplot gr [
    subplot 4 2 6
    title "Ring chart plot"
    light true
    # Will again use something like cylinder coordinates, but won't let the radius go to 0,
    # so y ∈ [-1, 1] with r = y+2
    xrange 0 1
    coords X="(y+2)/3*cos(2*pi*x)" Y="(y+2)/3*sin(2*pi*x)"
    rotate 50 100
    # Prepare data by constructing a random 7x2 matrix. So we compare
    # 2 datasets with 7 values each
    data = rand(7,2)+0.1
    # The last values will be made invisible (by the space ' ' in the style string),
    # define some fixed values for it so that one can see the rings internals
    data[7,:] = [0.25,0.5]
    # "bgrcmy ": use colors blue,green,red,cyan,mangenta,yellow,invisible
    # in this order for the 7 values supplied. '#': plot the edges as black lines
    chart data stl="bgrcmy #"
    # Reset the coordinate system
    coords
]

# Prepare some data that will be used in the next two plots
v = Array(linspace(-1,1,9))
x = Array(0:49)/49
y = Array(0:39)/39
a = 0.6sin(2pi*x).*sin(3pi*y') + 0.4cos(3pi*x.*y')
b = 0.6cos(2pi*x).*cos(3pi*y') + 0.4cos(3pi*x.*y')

# A surface plot with density plot at the bottom
@mglplot gr [
    subplot 4 2 7
    title "Surface plot"
    ranges -1 1 -1 1 -1 1
    light true
    rotate 50 60
    # Plots a bounding box around the 3d figure to be plotted
    box
    # Surface plot with "wire"-option '#' and optional parameter meshnum=12
    surf a stl="#" opt="meshnum 12"
    # density plot with "wire"-option '#' and optional parameter meshnum=12
    # and alpha=0.75 (meaning transparency). (The density plot can of course
    # also be given at other positions than at the bottom).
    dens a stl="#" opt="meshnum 12; alpha 0.75"
]

# Create a very nice contour plot using 3 different ways of drawing contours
@mglplot gr [
    subplot 4 2 8
    title "Contour plot"
    light true
    rotate 50 60
    box
    # This is responsible for the vertical areas
    contv a
    # This is responsible for the horizontal areas
    contf a
    # This is responsible for drawing black lines at the meeting edges
    # between vertical + horizontal areas
    cont  a stl="k"
]

MathGL.write(gr, "../graphs/main_example.png")


text:
#!/usr/local/bin/julia
using go

debug_log = joinpath(splitdir(@__FILE__)[1], "./gtp_debug.log")
model_dir = joinpath(splitdir(@__FILE__)[1], "../models")

f = open(debug_log, "w")  # stderr log
redirect_stderr(f)

function usage()
    println(STDOUT, "Usage: ./gtp_bot.jl {random|keras [model name]}")
    exit()
end

nargs = length(ARGS)
if nargs != 0
    if ARGS[1] == "random"
        policy = go.RandomPolicy()
    elseif ARGS[1] == "keras"
        if nargs != 2
            usage()
        end
        policy = go.KerasNetwork(model_dir, ARGS[2])
    else
        usage()
    end
else
    usage()
end

go.gtp(policy)


text:
using Fastback
using Test
using Dates

inst = Instrument("TICKER");

# generate data
data = Dict{Instrument, Vector{BidAsk}}();
dt = DateTime(2018, 1, 2, 9, 30, 0);
data[inst] = [
    BidAsk(dt + Second(1), 100.0, 101.0),
    BidAsk(dt + Second(2), 100.5, 102.0),
    BidAsk(dt + Second(3), 102.5, 103.0),
    BidAsk(dt + Second(4), 101.0, 102.0),
    BidAsk(dt + Second(5), 99.5, 100.0),
    BidAsk(dt + Second(6), 101.5, 102.0),
    BidAsk(dt + Second(7), 103.5, 104.0),
    BidAsk(dt + Second(8), 104.5, 105.0)
];

@testset "Dummy trading" begin

    # create trading account
    acc = Account(10_000.0)
    collect_balance, balance_curve = periodic_collector(Float64, Second(1));
    collect_equity, equity_curve = periodic_collector(Float64, Second(1));

    # dummy backtest
    for (i, ba) in enumerate(data[inst])
        if i == 1
            execute_order!(acc, OpenOrder(inst, 500.0, Long::TradeDir), ba)
        end

        update_account!(acc, inst, ba)

        if i == 7
            execute_order!(acc, CloseAllOrder(), ba)
        end

        # collect data for analysis
        collect_balance(ba.dt, acc.balance)
        collect_equity(ba.dt, acc.equity)
    end

    # display(acc)

end

@testset "Dummy trading #2" begin

    # create trading account
    acc = Account(100_000.0)
    execute_order!(acc, OpenOrder(inst, 100.0, Long::TradeDir), data[inst][1])  # 100.0 101.0
    update_account!(acc, inst, data[inst][2]) # 100.5, 102.0
    execute_order!(acc, CloseOrder(acc.open_positions[1]), data[inst][3]) # 102.5, 103.0
    update_account!(acc, inst, data[inst][4]) # 101.0, 102.0

    # display(acc)

    p::Position = acc.closed_positions[1]
    @test p.open_price == 101.0
    @test p.last_price == 102.5

    @test pnl_gross(p) ≈ p.size * (midprice(p.last_quote) - midprice(p.open_quote))
    @test pnl_net(p) ≈ p.size * (p.last_price - p.open_price)

    @test return_gross(p) ≈ (midprice(p.last_quote) - midprice(p.open_quote)) / midprice(p.open_quote)
    @test return_net(p) ≈ (p.last_price - p.open_price) / p.open_price

    @test length(acc.open_positions) == 0
    @test length(acc.closed_positions) == 1
    @test acc.initial_balance == 100_000.0
    @test acc.equity ≈ acc.initial_balance + pnl_net(p)
    @test acc.balance ≈ acc.equity

end


text:
using JLD

function costdatafilename()
    return "./data/cost_data.jld"
end

function savecost(cost)
    save(costdatafilename(), "cost", cost)
end

function loadcost()
    raw = load(costdatafilename())
    return raw["cost"]
end

max_stars_repo_name: qlete/ANManagement

max_issues_repo_name: qlete/ANManagement

max_forks_repo_name: qlete/ANManagement


text:
# Kafka wire protocol. See details at:
# https://cwiki.apache.org/confluence/display/KAFKA/A+Guide+To+The+Kafka+Protocol#AGuideToTheKafkaProtocol-Messagesets

struct RequestHeader
    api_key::Int16
    api_version::Int16
    correlation_id::Int32
    client_id::String
end

struct ResponseHeader
    correlation_id::Int32
end

struct Broker
    node_id::Int32
    host::String
    port::Int32
end

struct PartitionMetadata
    partition_error_code::Int16
    partition_id::Int32
    leader::Int32
    replicas::Vector{Int32}
    isr::Vector{Int32}
end

struct TopicMetadata
    topic_error_code::Int16
    topic_name::String
    partition_metadata::Vector{PartitionMetadata}
end

struct TopicMetadataRequest
    # RequestHeader omitted, write manually header::RequestHeader
    topics::Vector{String}
end

# to obtain metadata for all topics, we need to pass an empty array
# however, Kafka mixes up empty arrays and nulls and rest of API
# transforms Julia's T[] to null (i.e. array with lenght -1)
# to overcome it in this specific case we use special request that
# emulates empty array (with length 0)
struct AllTopicsMetadataRequest
    # RequestHeader omitted, write manuallyheader::RequestHeader
    pseudo_length::Int32
    AllTopicsMetadataRequest() = new(0)
end

struct TopicMetadataResponse
    # ResponseHeader omitted, read manually
    brokers::Vector{Broker}
    topic_metadata::Vector{TopicMetadata}
end

struct Message_v0
    crc::Int32
    magic_byte::Int8
    attributes::Int8
    key::Vector{UInt8}
    value::Vector{UInt8}
end
struct Message_v1
    crc::Int32
    magic_byte::Int8
    attributes::Int8
    timestamp::Int64
    key::Vector{UInt8}
    value::Vector{UInt8}
end
const Message = Message_v0

struct OffsetMessage
    offset::Int64
    message_size::Int32
    message::Message
end

# NOTE: MessageSet is serialized differently than other types,
# see io.jl for details
struct MessageSet
    elements::Vector{OffsetMessage}
end

# produce

struct ProduceRequestPartitionData
    partition::Int32
    message_set_size::Int32
    message_set::MessageSet
end

struct ProduceRequestTopicData
    topic_name::String
    partition_data::Vector{ProduceRequestPartitionData}
end

struct ProduceRequest
    # RequestHeader omitted, write manually
    required_acks::Int16
    timeout::Int32
    topic_data::Vector{ProduceRequestTopicData}
end

struct ProduceResponse_v0
    # ResponseHeader omitted, read it manually
    # responses format: [TopicName [Partition ErrorCode Offset]]
    responses::Vector{Tuple{String, Vector{Tuple{Int32,Int16,Int64}}}}
end
struct ProduceResponse_v1 # (supported in 0.9.0 or later)
    # ResponseHeader omitted, read it manually
    responses::Vector{Tuple{String, Vector{Tuple{Int32,Int16,Int64}}}}
    throttle_time::Int32
end
struct ProduceResponse_v2 # (supported in 0.10.0 or later)
    # ResponseHeader omitted, read it manually
    responses::Vector{Tuple{String, Vector{Tuple{Int32,Int16,Int64,Int64}}}}
    throttle_time::Int32
end
const ProduceResponse = ProduceResponse_v0

# fetch

struct FetchRequestPartitionData
    partition::Int32
    offset::Int64
    max_bytes::Int32
end

struct FetchRequestTopicData
    topic_name::String
    partition_fetches::Vector{FetchRequestPartitionData}
end

struct FetchRequest
    # RequestHeader omitted, write it manually
    replica_id::Int32 # should always be -1 for clients
    max_wait_time::Int32
    min_bytes::Int32
    topic_fetches::Vector{FetchRequestTopicData}
end

struct FetchResponsePartitionData
    partition::Int32
    error_code::Int16
    highwater_mark_offset::Int64
    message_set_size::Int32
    message_set::MessageSet
end

struct FetchResponseTopicData
    topic_name::String
    partition_results::Vector{FetchResponsePartitionData}
end

struct FetchResponse_v0
    # ResponseHeader omitted, read it manually
    topic_results::Vector{FetchResponseTopicData}
end

struct FetchResponse_v1
    # ResponseHeader omitted, read it manually
    throttle_time::Int32
    topic_results::Vector{FetchResponseTopicData}
end

const FetchResponse = FetchResponse_v0

# offset listing

struct OffsetRequestPartitionData
    partition::Int32
    time::Int64
    max_number_of_offsets::Int32
end

struct OffsetRequestTopicData
    topic_name::String
    partition_data::Vector{OffsetRequestPartitionData}
end

struct OffsetRequest
    replica_id::Int32
    topic_data::Vector{OffsetRequestTopicData}
end

struct OffsetResponsePartitionData
    partition::Int32
    error_code::Int16
    offsets::Vector{Int64}
end

struct OffsetResponseTopicData
    topic_name::String
    partition_data::Vector{OffsetResponsePartitionData}
end

struct OffsetResponse
    topic_data::Vector{OffsetResponseTopicData}
end

# API versions

struct ApiVersion
    api_key::Int16
    min_version::Int16
    max_version::Int16
end

struct ApiVersionsResponse
    error_code::Int16
    api_versions::Vector{ApiVersion}
    throttle_time_ms::Int32
end


text:
using HWfuncapp
using Test

@testset "HWfuncapp.jl" begin
	# test that q1 with 15 nodes has an error smaller than 1e-9

    # test that the integral of function h [-10,0] ≈ 1.46039789878568

    @testset "q1" begin
		out = q1()
		@test out["error"] < 1e-9
end


text:
@testset "Selections" begin

    rng = StableRNG(42)

    @testset "Rank" begin
        fitness = [5, 1, 3, 4, 2.0]

        s = ranklinear(1.0)
        Random.seed!(rng, 2)
        @test s(fitness, 2, rng=rng) == [5,1]

        s = ranklinear(2.0)
        Random.seed!(rng, 2)
        @test s([1,2], 2, rng=rng) == [2,2]
    end

    @testset "Uniform" begin
        s = uniformranking(2)
        Random.seed!(rng, 2);
        @test sort(unique(s([1.0,2.0,3.0], 10, rng=rng))) == [1,2]
        Random.seed!(rng, 2);
        @test sort(unique(s([5,2,3], 5, rng=rng))) == [2,3]
        @test_throws AssertionError s([1.], 4)
    end

    @testset "Roulette" begin
        Random.seed!(rng, 2);
        @test roulette([30.0, -0.1,-0.2, -30.0], 2, rng=rng) == [4,1]
        Random.seed!(rng, 2);
        @test roulette([0.0,0.0,3.0], 2, rng=rng) == [3, 3]
        Random.seed!(rng, 2);
        @test roulette([0,2,0], 2, rng=rng) == [2,2]
        Random.seed!(rng, 2);
        @test rouletteinv([0.0,0.0,3.0], 2, rng=rng) == [1, 1]
        Random.seed!(rng, 2);
        @test rouletteinv([1,2,0], 2, rng=rng) == [3,3]
    end

    @testset "Truncation" begin
        @test truncation([0.0,2.0,3.0], 2) == [1,2]
        @test truncation([0,2,1], 1) == [1]
        @test_throws AssertionError truncation([1.,2.], 5)
    end

    @testset "Tournament" begin
        @test_throws AssertionError tournament(0)
        t = tournament(3, select=argmax)
        Random.seed!(rng, 2);
        @test all(t([0,2,0],100,rng=rng) .== 2)
        t = tournament(3)
        Random.seed!(rng, 2);
        @test all(i->i∈[1,3], t([0,2,0],100,rng=rng))
        t = tournament(2)
        Random.seed!(rng, 2);
        @test all(t([0.0,0.0,1.0],100,rng=rng) .< 3.0)
        fitness = [0 0 1 1; 0 1 0 1]
        t = tournament(2, select=Evolutionary.twowaycomp)
        Random.seed!(rng, 2);
        @test all(t(fitness,100,rng=rng) .< 4)
        @test mean(t(fitness,100,rng=rng)) < 2
    end

    @testset "SUS" begin
        Random.seed!(rng, 2);
        @test sort(unique(sus([1.0,0.0,2.0], 5, rng=rng))) == [1,3]
        Random.seed!(rng, 2);
        @test sort(unique(sus([0,1,2], 5, rng=rng))) == [2,3]
    end

    @testset "DE" begin
        N = 5
        col = collect(1:10)

        Random.seed!(rng, 1)
        idxs = random(col, N, rng=rng)
        @test length(idxs) == N
        @test all(i-> i ∈ col, idxs)
        @test !all(length(unique(random(col, N, rng=rng))) == N for i in 1:100)

        Random.seed!(rng, 1)
        idxs = randomoffset(col, N, rng=rng)
        @test length(idxs) == N
        @test all(i-> i ∈ col, idxs)
        @test all( let idxs = randomoffset(col, N, rng=rng); let l = (idxs[1] - 1 + N)%10; idxs[N] == (l == 0 ? 10 : l) end end for i in 1:100 )

        idxs = best(reverse(col), N)
        @test length(idxs) == N
        @test all(i-> i ∈ col, idxs)
        @test idxs == fill(10, N)

        Random.seed!(rng, 1)
        idxs = permutation(col, N, rng=rng)
        @test length(idxs) == N
        @test all(i-> i ∈ col, idxs)
        @test all(length(unique(permutation(col, N, rng=rng))) == N for i in 1:100)

        Random.seed!(rng, 1)
        @testset "Mutually Exclusive Indices" for i in idxs
            targets = Evolutionary.randexcl(rng, col, [i], 4)
            @test length(unique(push!(targets, i))) == N
        end
    end

end


text:
# # Train Kernel Parameters

# Here we show a few ways to train (optimize) the kernel (hyper)parameters at the example of kernel-based regression using KernelFunctions.jl.
# All options are functionally identical, but differ a little in readability, dependencies, and computational cost.

# We load KernelFunctions and some other packages. Note that while we use `Zygote` for automatic differentiation and `Flux.optimise` for optimization, you should be able to replace them with your favourite autodiff framework or optimizer.

using KernelFunctions
using LinearAlgebra
using Distributions
using Plots
using BenchmarkTools
using Flux
using Flux: Optimise
using Zygote
using Random: seed!
seed!(42);

# ## Data Generation
# We generate a toy dataset in 1 dimension:

xmin, xmax = -3, 3  # Bounds of the data
N = 50 # Number of samples
x_train = rand(Uniform(xmin, xmax), N)  # sample the inputs
σ = 0.1
y_train = sinc.(x_train) + randn(N) * σ  # evaluate a function and add some noise
x_test = range(xmin - 0.1, xmax + 0.1; length=300)
nothing #hide

# Plot the data

scatter(x_train, y_train; label="data")
plot!(x_test, sinc; label="true function")

# ## Manual Approach
# The first option is to rebuild the parametrized kernel from a vector of parameters
# in each evaluation of the cost function. This is similar to the approach taken in
# [Stheno.jl](https://github.com/JuliaGaussianProcesses/Stheno.jl).

# To train the kernel parameters via [Zygote.jl](https://github.com/FluxML/Zygote.jl),
# we need to create a function creating a kernel from an array.
# A simple way to ensure that the kernel parameters are positive
# is to optimize over the logarithm of the parameters.

function kernelcall(θ)
    return (exp(θ[1]) * SqExponentialKernel() + exp(θ[2]) * Matern32Kernel()) ∘
           ScaleTransform(exp(θ[3]))
end
nothing #hide

# From theory we know the prediction for a test set x given
# the kernel parameters and normalization constant:

function f(x, x_train, y_train, θ)
    k = kernelcall(θ[1:3])
    return kernelmatrix(k, x, x_train) *
           ((kernelmatrix(k, x_train) + exp(θ[4]) * I) \ y_train)
end
nothing #hide

# Let's look at our prediction.
# With starting parameters `p0` (picked so we get the right local
# minimum for demonstration) we get:

p0 = [1.1, 0.1, 0.01, 0.001]
θ = log.(p0)
ŷ = f(x_test, x_train, y_train, θ)
scatter(x_train, y_train; label="data")
plot!(x_test, sinc; label="true function")
plot!(x_test, ŷ; label="prediction")

# We define the following loss:

function loss(θ)
    ŷ = f(x_train, x_train, y_train, θ)
    return norm(y_train - ŷ) + exp(θ[4]) * norm(ŷ)
end
nothing #hide

# The loss with our starting point:

loss(θ)

# Computational cost for one step:

@benchmark let
    θ = log.(p0)
    opt = Optimise.ADAGrad(0.5)
    grads = only((Zygote.gradient(loss, θ)))
    Optimise.update!(opt, θ, grads)
end

# ### Training the model

# Setting an initial value and initializing the optimizer:
θ = log.(p0) # Initial vector
opt = Optimise.ADAGrad(0.5)
nothing #hide

# Optimize

anim = Animation()
for i in 1:15
    grads = only((Zygote.gradient(loss, θ)))
    Optimise.update!(opt, θ, grads)
    scatter(
        x_train, y_train; lab="data", title="i = $(i), Loss = $(round(loss(θ), digits = 4))"
    )
    plot!(x_test, sinc; lab="true function")
    plot!(x_test, f(x_test, x_train, y_train, θ); lab="Prediction", lw=3.0)
    frame(anim)
end
gif(anim, "train-kernel-param.gif"; show_msg=false, fps=15);
nothing; #hide

# ![](train-kernel-param.gif)

# Final loss
loss(θ)

# ## Using ParameterHandling.jl
# Alternatively, we can use the [ParameterHandling.jl](https://github.com/invenia/ParameterHandling.jl) package
# to handle the requirement that all kernel parameters should be positive.
# The package also allows arbitrarily nesting named tuples that make the parameters
# more human readable, without having to remember their position in a flat vector.

using ParameterHandling

raw_initial_θ = (
    k1=positive(1.1), k2=positive(0.1), k3=positive(0.01), noise_var=positive(0.001)
)

flat_θ, unflatten = ParameterHandling.value_flatten(raw_initial_θ)
flat_θ #hide

# We define a few relevant functions and note that compared to the previous `kernelcall` function, we do not need explicit `exp`s.

function kernelcall(θ)
    return (θ.k1 * SqExponentialKernel() + θ.k2 * Matern32Kernel()) ∘ ScaleTransform(θ.k3)
end
nothing #hide

function f(x, x_train, y_train, θ)
    k = kernelcall(θ)
    return kernelmatrix(k, x, x_train) *
           ((kernelmatrix(k, x_train) + θ.noise_var * I) \ y_train)
end
nothing #hide

function loss(θ)
    ŷ = f(x_train, x_train, y_train, θ)
    return norm(y_train - ŷ) + θ.noise_var * norm(ŷ)
end
nothing #hide

initial_θ = ParameterHandling.value(raw_initial_θ)
nothing #hide

# The loss at the initial parameter values:

(loss ∘ unflatten)(flat_θ)

# Cost per step

@benchmark let
    θ = flat_θ[:]
    opt = Optimise.ADAGrad(0.5)
    grads = (Zygote.gradient(loss ∘ unflatten, θ))[1]
    Optimise.update!(opt, θ, grads)
end

# ### Training the model

# Optimize

opt = Optimise.ADAGrad(0.5)
for i in 1:15
    grads = (Zygote.gradient(loss ∘ unflatten, flat_θ))[1]
    Optimise.update!(opt, flat_θ, grads)
end
nothing #hide

# Final loss

(loss ∘ unflatten)(flat_θ)

# ## Flux.destructure
# If we don't want to write an explicit function to construct the kernel, we can alternatively use the `Flux.destructure` function.
# Again, we need to ensure that the parameters are positive. Note that the `exp` function is now part of the loss function, instead of part of the kernel construction.

# We could also use ParameterHandling.jl here.
# To do so, one would remove the `exp`s from the loss function below and call `loss ∘ unflatten` as above.

θ = [1.1, 0.1, 0.01, 0.001]

kernel = (θ[1] * SqExponentialKernel() + θ[2] * Matern32Kernel()) ∘ ScaleTransform(θ[3])

params, kernelc = Flux.destructure(kernel);

# This returns the trainable `params` of the kernel and a function to reconstruct the kernel.
kernelc(params)

# From theory we know the prediction for a test set x given
# the kernel parameters and normalization constant

function f(x, x_train, y_train, θ)
    k = kernelc(θ[1:3])
    return kernelmatrix(k, x, x_train) * ((kernelmatrix(k, x_train) + (θ[4]) * I) \ y_train)
end
nothing #hide

function loss(θ)
    ŷ = f(x_train, x_train, y_train, exp.(θ))
    return norm(y_train - ŷ) + exp(θ[4]) * norm(ŷ)
end
nothing #hide

# Cost for one step

@benchmark let θt = θ[:], optt = Optimise.ADAGrad(0.5)
    grads = only((Zygote.gradient(loss, θt)))
    Optimise.update!(optt, θt, grads)
end

# ### Training the model

# The loss at our initial parameter values:
θ = log.([1.1, 0.1, 0.01, 0.001]) # Initial vector
loss(θ)

# Initialize optimizer

opt = Optimise.ADAGrad(0.5)
nothing #hide

# Optimize

for i in 1:15
    grads = only((Zygote.gradient(loss, θ)))
    Optimise.update!(opt, θ, grads)
end
nothing #hide

# Final loss
loss(θ)


text:
abstract type AbstractDevice end

struct CPU <: AbstractDevice
end

struct GPU <: AbstractDevice
    id::Int
end

GPU() = GPU(1)

to_device(device::CPU, x::CuArray) = convert(Array, x)

function to_device(device::GPU, x)
    T = typeof(x)
    flds = fieldnames(T)
    if isa(x, CuArray)
        return x
    elseif isa(x, AbstractFloat)
        return Float32(x)
    elseif isa(x, Tuple)
        return ((to_device(device, el) for el in x)...,)
    elseif isempty(flds)
        # primitive or array
        return cu(x)
    else
        # struct, recursively convert and construct type from fields
        fld_vals = [to_device(device, getfield(x, fld)) for fld in flds]
        return T(fld_vals...)
    end
end

# is_cuarray(x) = x isa CuArray

## currently GPU's ID is just a placeholder
# guess_device(args) = any(is_cuarray, args) ? GPU(1) : CPU()
# device_of(A) = A isa CuArray ? GPU(1) : CPU()

# """
# Retrieve function compatible with specified device
# See also: to_device(device, f)
# """
# device_function(device::CPU, f) = f

"""
Convert object to a compatible with the specified device.
For CPU it's usually no-op. For GPU behavior differs between object types:
 * Arrays are converted to CuArrays
 * structs are converted recursively
 * functions are looked up using `device_function()` or transformed using tracer
 * all other objects are returned as is
"""
to_device(device::CPU, x) = x
to_device(device::CPU, f::Function, args) = f

(device::AbstractDevice)(x) = to_device(device, x)

# to_same_device(A, example) = device_of(example)(A)


text:
push!(LOAD_PATH, joinpath(pwd(), "src"))
using Test
using JuML

train1Mcsv = joinpath("data", "airlinetrain1m.csv")
testcsv = joinpath("data", "airlinetest.csv")

importcsv(train1Mcsv; path2 = testcsv, outname = "airlinetraintest")

traintest_df = DataFrame(joinpath("data", "airlinetraintest"))

@test length(traintest_df.covariates) == 2

@test length(traintest_df.factors) == 7

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
### PLugin
#=
Whenever the plugin is updated, run this function and commit the plugin files.
upgrade_plugin()
=#

const khepri_grasshopper_dlls = ["KhepriGrasshopper.gha", "FastColoredTextBox.dll"]
const julia_khepri = dirname(dirname(abspath(@__FILE__)))

upgrade_plugin() =
  let # 1. The dlls are updated in VisualStudio after compilation of the plugin, and they are stored in the folder
      #    contained inside the Plugins folder, which has a specific location regarding this file itself
      plugin_folder = joinpath(dirname(dirname(julia_khepri)), "Plugins", "KhepriGrasshopper", "KhepriGrasshopper", "bin")
      # 2. The bundle needs to be copied to the current folder
      local_folder = joinpath(julia_khepri, "Plugin")
      # 3. Now we copy the dlls to the local folder
      for dll in khepri_grasshopper_dlls
          src = joinpath(plugin_folder, dll)
          dst = joinpath(local_folder, dll)
          rm(dst, force=true)
          cp(src, dst)
      end
  end

#

update_plugin() =
  let grasshopper_user_plugins = joinpath(ENV["APPDATA"], "Grasshopper", "Libraries"),
      local_khepri_plugin = joinpath(julia_khepri, "Plugin")
    for dll in khepri_grasshopper_dlls
      let local_path = joinpath(local_khepri_plugin, dll),
          grasshopper_path = joinpath(grasshopper_user_plugins, dll)
          cp(local_path, grasshopper_path, force=true)
      end
    end
  end

#=
We want to support these two syntaxes

a < Number
a < Number("Number")
a < Number("Number", "N")
a < Number("Number", "N", "Parameter N")
a < Number("Number", "N", "Parameter N", 0.0)

a < Number(0.0)
a < Number("Number", 0.0)
a < Number("Number", "N", 0.0)
a < Number("Number", "N", "Parameter N", 0.0)

a < String()
a < String("String")
a < String("String")
a < Number(0.0, "a")

=#

in_gh(sym) = Symbol("GH$(sym)")
kgh_io_function_names = Symbol[]
is_kgh_io_function_name(sym) =
  sym in kgh_io_function_names

macro ghdef(name, init)
  let str = string(name),
      ghname = esc(in_gh(name))
    quote
      push!(kgh_io_function_names, $(QuoteNode(name)))
      $(ghname)(description::Base.String, short_description=description[1:1], message=description*" parameter", value=$init) =
        [$(str), description, short_description, message, value]
      $(ghname)(value::Base.Any, description=$(str), short_description=description[1:1], message=description*" parameter") =
        [$(str), description, short_description, message, value]
    end
  end
end

#Strings require conversion to avoid SubStrings!!!!!!!
@ghdef(String, "")
@ghdef(Path, "")
@ghdef(Boolean, false)
@ghdef(Number, 0.0)
@ghdef(Integer, 0)
@ghdef(Point, u0())
@ghdef(Vector, vx(1))
@ghdef(Any, nothing)
@ghdef(Eval, nothing)
@ghdef(JL, nothing)
@ghdef(Strings, [])
@ghdef(Paths, [])
@ghdef(Booleans, [])
@ghdef(Numbers, [])
@ghdef(Integers, [])
@ghdef(Points, [])
@ghdef(Vectors, [])
@ghdef(Many, [])
@ghdef(Evals, [])
@ghdef(JLs, [])
@ghdef(Stringss, [])
@ghdef(Pathss, [])
@ghdef(Booleanss, [])
@ghdef(Numberss, [])
@ghdef(Integerss, [])
@ghdef(Pointss, [])
@ghdef(Vectorss, [])
@ghdef(Manies, [])
@ghdef(Evalss, [])
@ghdef(JLss, [])

export define_kgh_function

kgh_forms(text, idx=1) =
  let (expr, idx) = Meta.parse(text, idx, greedy=true, depwarn=false)
    isnothing(expr) ?
      [] :
      [expr, kgh_forms(text, idx)...]
  end

is_kgh_io_function_call(e) =
  e isa Expr && e.head === :call && is_kgh_io_function_name(e.args[1])

match_expr(e1, e2) =
  e2 === :_ ||
  e1 == e2 ||
  e1 isa Expr && e2 isa Expr &&
  match_expr(e1.head, e2.head) &&
  length(e1.args) == length(e2.args) &&
  all([match_expr(e1, e2) for (e1, e2) in zip(e1.args, e2.args)])

is_kgh_input(form) =
  match_expr(form, :(_ < _)) &&
  is_kgh_io_function_call(form.args[3])
is_kgh_output(form) =
  (match_expr(form, :(_ > _)) &&
   is_kgh_io_function_call(form.args[3])) ||
  is_kgh_io_function_call(form)

kgh_io_param(form) =
  form.args[2] == :_ ? :__result : form.args[2]
kgh_io_call(form) =
  let param = form.args[2],
      form = form.args[3],
      func = Expr(:., :KhepriGrasshopper, QuoteNode(in_gh(form.args[1])))
    match_expr(form, :(_())) ||
    (match_expr(form, :(_(_))) && !(form.args[2] isa String)) ?
        :($(func)($(form.args[2:end]...), $(string(param)))) :
        :($(func)($(form.args[2:end]...)))
  end

kgh_inputs(forms) = filter(is_kgh_input, forms)
kgh_outputs(forms) = filter(is_kgh_output, forms)

create_kgh_function(name::String, body::String) =
  let forms = kgh_forms(body),
      inputs = filter(is_kgh_input, forms),
      outputs = filter(is_kgh_output, forms),
      forms = filter(f -> !(f in inputs || f in outputs), forms),
      inp_params = map(kgh_io_param, inputs),
      out_params = map(kgh_io_param, outputs),
      inp_forms = map(kgh_io_call, inputs),
      out_forms = map(kgh_io_call, outputs),
      inps = Symbol("__inps_$name"),
      inp_inits = [:($(inp) = $(inps)[$(i)]) for (i, inp) in enumerate(inp_params)],
      outs = Symbol("__outs_$name"),
      out_inits = [:($(outs)[$(i)] = $(out)) for (i, out) in enumerate(out_params)],
      inp_docs = Symbol("__doc_inps_$name"),
      out_docs = Symbol("__doc_outs_$name"),
      shapes = Symbol("__shapes_$name"),
      prev_collected_shapes = gensym("prev_collected_shapes")
    if isempty(inp_forms) && isempty(out_forms)
      quote
        $(inp_docs) = Any[$(inp_forms...)]
        $(out_docs) = Any[$(out_forms...)]
        $(inps) = Array{Any,1}(undef, $(length(inp_inits)))
        $(outs) = Array{Any,1}(undef, $(length(out_inits)))
        $(shapes) = Shape[]
        function $(Symbol("__func_$name"))()
          nothing
        end
      end
    else
      quote
        $(inp_docs) = Any[$(inp_forms...)]
        $(out_docs) = Any[$(out_forms...)]
        $(inps) = Array{Any,1}(undef, $(length(inp_inits)))
        $(outs) = Array{Any,1}(undef, $(length(out_inits)))
        $(shapes) = Shape[]
        function $(Symbol("__func_$name"))()
          $(inp_inits...)
          $prev_collected_shapes = collected_shapes()
          #collected_shapes($(shapes))
          __result = #try
            begin
            $(forms...)
            end
          #finally
            #collected_shapes($prev_collected_shapes)
          #end
          $(out_inits...)
          nothing
        end
      end
    end
  end

define_kgh_function(name::String, body::String) =
  Base.eval(Main, create_kgh_function(name, body))

#=
fn = """
  a < Number("Number", "N", "Parameter N", 0.0)
  b < Any()
  c < Number(5)
  d < String("Foo")

  f(a + b)

  c = a + 1

  c > Any()
  d > Strings("Bar")
"""

create_kgh_function("foo", fn)
=#

#=
fn = """
a < Number("Number", "N", "Parameter N", 0.0)
b < Any()
c < Number(5)
d < String("Foo")

f(a + b)

c = a + 1

c > Any()
d > String("Bar")
"""

create_kgh_function("foo", fn)

=#

#=
fn = """
a < Number("Number", "N", "Parameter N", 0.0)
b < Any()
c < Number(5)
d < String("Foo")
e > String("Bar")
_ > Number()

e = f(a + b)
g(b*c)
"""

create_kgh_function("foo", fn)

fn2 = """
e = f(a + b)
g(b*c)
"""

create_kgh_function("foo", fn2)
using InteractiveUtils
InteractiveUtils.@code_lowered KhepriBase.foo()
=#
#=
fn = """
a < Number("Number", "N", "Parameter N", 0.0)
b < Any()
c < Number(5)
d < String("Foo")

f(a + b)

c = a + 1

c > Any()
d > String("Bar")
"""

eval(create_kgh_function("foo", fn))
=#

#=
create_kgh_function("foo", raw"""
path < String()
autos < Numbers()
costs < Numbers()

open(path, "w") do f
  for (a, c) in zip(autos, costs)
    println(f, "$a $c")
  end
end
""")

dump(:("foo $bar"))

=#

#=
str1 = """# v < Type("Input", "I", "Parameter I", default)
a < Number(2)
b < Number(3)
_ > Number()

sqrt(a^2 + b^2)"""

create_kgh_function("zzz", str1)
=#

#=
a < Number()
b < Number()
sphere(x(a), b)

create_kgh_function("zzz", str1)
=#
#=
example = raw"""
generate_outline_vertices(p, area, prop) =
  let l = sqrt(area/prop)
      w = area/l
      [add_xy(p, 0, 0), add_xy(p, w, 0), add_xy(p, w, l), add_xy(p, 0, l)]
  end

areas = 50:50:200
wwrs = 0.2:0.2:0.4
rots = 0:pi/4:pi/2
props = 1:2:5
epocas = [1, 3, 4, 5, 6, 7, 8, 9]
kinds = [1, 2]
stors = 2:3:11
ddy = sqrt(maximum(areas))

Wmaterials = [[[0, "Reboco - 2cm", "Pedra", "Estuque_Claro_1.5"], [15.7, "Reboco - 2cm", "OSB - 2cm", "AirGap", "Isolamento XPS - 4cm", "Pedra", "Estuque_Claro_1.5" ]],
              [[0, "Reboco - 2cm", "Pedra", "Estuque_Claro_1.5"], [15.7, "Reboco - 2cm", "OSB - 2cm", "AirGap", "Isolamento XPS - 4cm", "Pedra", "Estuque_Claro_1.5" ]],
              [[0, "Reboco - 2cm", "Pedra", "AirGap", "Tijolo Furado_11", "Estuque_Claro_1.5"], [7, "Reboco - 2cm", "Pedra", "AirGap", "Isolamento XPS - 4cm", "Tijolo Furado_11", "Estuque_Claro_1.5"]],
              [[0, "Reboco - 2cm", "Tijolo Furado_15", "AirGap", "Tijolo Furado_11", "Estuque_Claro_1.5"], [7, "Reboco - 2cm", "Tijolo Furado_15", "AirGap", "Isolamento XPS - 4cm", "Tijolo Furado_11", "Estuque_Claro_1.5"]],
              [[0, "Reboco - 2cm", "Tijolo Furado_11", "AirGap", "Tijolo Furado_11", "Estuque_Claro_1.5"], [7, "Reboco - 2cm", "Tijolo Furado_11", "AirGap", "Isolamento XPS - 4cm", "Tijolo Furado_11", "Estuque_Claro_1.5"]],
              [[0, "Reboco - 2cm", "Tijolo Furado_15", "AirGap", "Isolamento XPS - 4cm", "Tijolo Furado_11", "Estuque_Claro_1.5"], [7, "Reboco - 2cm", "Tijolo Furado_15", "AirGap", "Isolamento XPS - 4cm", "Isolamento XPS - 4cm", "Tijolo Furado_11", "Estuque_Claro_1.5"]],
              [[0, "Reboco - 2cm", "OSB - 2cm", "AirGap", "Isolamento XPS - 4cm", "Tijolo Furado_22", "Estuque_Claro_1.5"], [7, "Reboco - 2cm", "OSB - 2cm", "AirGap", "Isolamento XPS - 4cm", "Isolamento XPS - 4cm", "Tijolo Furado_22", "Estuque_Claro_1.5"]],
              [[0, "Reboco - 2cm", "Isolamento EPS - 6cm", "Tijolo Furado_22", "Estuque_Claro_1.5"], [7, "Reboco - 2cm", "OSB - 2cm", "Isolamento XPS - 4cm", "Isolamento EPS - 6cm", "Tijolo Furado_22", "Estuque_Claro_1.5"]],
              [[0, "Chapa Metalica", "Isolamento EPS - 4cm", "Isolamento EPS - 6cm", "Isolamento EPS - 6cm", "Estuque_Claro_1.5"], [7, "Chapa Metalica", "Isolamento XPS - 4cm", "Isolamento EPS - 4cm", "Isolamento EPS - 6cm", "Isolamento EPS - 6cm", "Estuque_Claro_1.5"]],
              [[0, "Chapa Metalica", "Isolamento EPS - 4cm", "Isolamento EPS - 6cm", "Isolamento EPS - 6cm", "Estuque_Claro_1.5"], [7, "Chapa Metalica", "Isolamento XPS - 4cm", "Isolamento EPS - 4cm", "Isolamento EPS - 6cm", "Isolamento EPS - 6cm", "Estuque_Claro_1.5"]]]

Fmaterials = [[[0, "Paineis de Madeira_12", "Estuque_Claro_1.5"]],
              [[0, "Paineis de Madeira_12", "Estuque_Claro_1.5"]],
              [[0, "Paineis de Madeira_12", "Estuque_Claro_1.5"]],
              [[0, "Ceramica vidrada - 1cm", "Betonilha de Acentamento_8", "Laje Betao_15", "Estuque_Claro_1.5"]],
              [[0, "Ceramica vidrada - 1cm", "Betonilha de Acentamento_8", "Laje Betao_15", "Estuque_Claro_1.5"]],
              [[0, "Ceramica vidrada - 1cm", "Betonilha de Acentamento_8", "Laje Betao_15", "Estuque_Claro_1.5"]],
              [[0, "Ceramica vidrada - 1cm", "Betonilha de Acentamento_8", "Laje Betao_15", "Estuque_Claro_1.5"]],
              [[0, "Ceramica vidrada - 1cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Ceramica vidrada - 1cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Ceramica vidrada - 1cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]]]

Rmaterials = [[[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]],
              [[0, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"], [7.5, "Betonilha de Acentamento_8", "Tela impermeabilizacao - 2mm", "Isolamento XPS - 4cm", "Betonilha de Acentamento_8", "Laje Aligeirada_0.25", "Estuque_Claro_1.5"]]]

Windowmaterials = [[[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]],
                   [[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]],
                   [[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]],
                   [[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]],
                   [[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]],
                   [[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]],
                   [[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]],
                   [[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]],
                   [[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]],
                   [[[0, "2.69, 0.75, 0.8"], [50, "1.70, 0.38, 0.7"]]]]
vertss = []
grsn = []
grse = []
grss = []
grsw = []
Wmats = []
Rmats = []
Fmats = []
Windowmats = []
Storeys = []
dy = 0

for epoca in epocas
  dx = 0
  for wall in kinds, roof in kinds, window in kinds
    for wwr in wwrs
      for area in areas
        for rot in rots
          for prop in props, stor in stors
            for wn in [0, wwr], we in [0, wwr], ws in [0, wwr], ww in [0, wwr]
              push!(vertss, generate_outline_vertices(loc_from_o_phi(xy(dx, dy), rot), area, prop))
              push!(grsn, wn)
              push!(grse, we)
              push!(grss, ws)
              push!(grsw, ww)
              push!(Wmats, Wmaterials[epoca][wall])
              push!(Rmats, Rmaterials[epoca][roof])
              push!(Windowmats, Windowmaterials[epoca][window])
              push!(Fmats, Fmaterials[epoca][1])
              push!(Storeys, stor)
              dx += sqrt(area) + 10
            end
          end
        end
      end
    end
  end
  dy += ddy + 10
end

vertss > JLs("Verticess JL")
grsn > Numbers("WWR'North")
grse > Numbers("WWR'East")
grss > Numbers("WWR'South")
grsw > Numbers("WWR'West")
Storeys > Integers("Number of floors")
Wmats > JLs("Wall materials JL")
Rmats > JLs("Roof materials JL")
Fmats > JLs("Floor materials JL")
Windowmats > JLs("Window materials JL")
"""

print(create_kgh_function("zzz", example))
=#


text:
using TriangleMesh
using Plots
using Triplot

poly = polygon_unitSquare()
# poly = polygon_regular(5)
# poly = polygon_Lshape()
mesh = create_mesh(poly, set_area_max=true, quality_meshing=true)

VX,VY = mesh.point[1,:],mesh.point[2,:]
EToV  = mesh.cell # element-to-vertex mapping

function uniform_tri_nodes(N)
    r1D = LinRange(-1,1,N+1)
    Np = (N+1)*(N+2)÷2
    r,s = zeros(Np),zeros(Np)
    sk = 1
    for i = 0:N
        for j = 0:N-i
            r[sk] = r1D[i+1]
            s[sk] = r1D[j+1]
            sk += 1
        end
    end
    return r,s
end
r,s = uniform_tri_nodes(20)

λ1(r,s) = -(r+s)/2
λ2(r,s) = (1+s)/2
λ3(r,s) = (1+r)/2
V(r,s) = [λ1.(r,s) λ2.(r,s) λ3.(r,s)]
dVr() = [-.5 .5 0.0]
dVs() = [-.5 0.0 .5]
function map_triangle_pts(r,s,x,y)
    return V(r,s)*x,V(r,s)*y
end

function compute_geometric_terms(x,y)
    dxdr,dydr = dVr()*x, dVr()*y
    dxds,dyds = dVs()*x, dVs()*y
    G = [dxdr dxds; dydr dyds] # G*r = x, r = inv(G)*x
    J = det(G)
    drdx, dsdx, drdy, dsdy = inv(G)
    return J, drdx, dsdx, drdy, dsdy
end

# A = spzeros(length(VX),length(VX))
# for e = 1:mesh.n_cell
#     xv,yv = x[mesh.cell[:,e]],y[mesh.cell[:,e]]
#     J,rx,sx,ry,sy = compute_geometric_terms(xv,yv)
#
#     Vx = rx*dVr() + sx*dVs()
#     Vy = ry*dVr() + sy*dVs()
#     ids = mesh.cell[:,e]
#     @. A[ids,ids] += 2.0*J*(Vx'*Vx + Vy'*Vy)
#     # @show J
# end

scatter(VX,VY)
for e = 1:mesh.n_cell
    xv,yv = VX[mesh.cell[:,e]],VY[mesh.cell[:,e]]
    xe,ye = map_triangle_pts(r,s,xv,yv)
    zz = V(r,s)*(yv)
    scatter!(xe,ye,leg=false,ms=1)
end
display(plot!())


text:
using Base: UUID, SHA1
using TOML
using ..Versions: VersionSpec, VersionRange
using ..LazilyInitializedFields

# The content of a registry is assumed to be constant during the
# lifetime of a `Registry`. Create a new `Registry` if you want to have
# a new view on the current registry.

# See loading.jl
const TOML_CACHE = Base.TOMLCache(TOML.Parser(), Dict{String, Dict{String, Any}}())
const TOML_LOCK = ReentrantLock()
parsefile(toml_file::AbstractString) = Base.parsed_toml(toml_file, TOML_CACHE, TOML_LOCK)

# Info about each version of a package
@lazy mutable struct VersionInfo
    git_tree_sha1::Base.SHA1
    yanked::Bool
    @lazy uncompressed_compat::Union{Dict{UUID, VersionSpec}}
end
VersionInfo(git_tree_sha1::Base.SHA1, yanked::Bool) = VersionInfo(git_tree_sha1, yanked, uninit)

# This is the information that exists in e.g. General/A/ACME
struct PkgInfo
    # Package.toml:
    repo::Union{String, Nothing}
    subdir::Union{String, Nothing}

    # Versions.toml:
    version_info::Dict{VersionNumber, VersionInfo}

    # Compat.toml
    compat::Dict{VersionRange, Dict{String, VersionSpec}}

    # Deps.toml
    deps::Dict{VersionRange, Dict{String, UUID}}
end

isyanked(pkg::PkgInfo, v::VersionNumber) = pkg.version_info[v].yanked
treehash(pkg::PkgInfo, v::VersionNumber) = pkg.version_info[v].git_tree_sha1

function uncompress(compressed::Dict{VersionRange, Dict{String, T}}, vsorted::Vector{VersionNumber}) where {T}
    @assert issorted(vsorted)
    uncompressed = Dict{VersionNumber, Dict{String, T}}()
    for v in vsorted
        uncompressed[v] = Dict{String, T}()
    end
    for (vs, data) in compressed
        first = length(vsorted) + 1
        # We find the first and last version that are in the range
        # and since the versions are sorted, all versions in between are sorted
        for i in eachindex(vsorted)
            v = vsorted[i]
            v in vs && (first = i; break)
        end
        last = 0
        for i in reverse(eachindex(vsorted))
            v = vsorted[i]
            v in vs && (last = i; break)
        end
        for i in first:last
            v = vsorted[i]
            uv = uncompressed[v]
            for (key, value) in data
                if haskey(uv, key)
                    # Change to an error?
                    error("Overlapping ranges for $(key) in $(repr(path)) for version $v.")
                else
                    uv[key] = value
                end
            end
        end
    end
    return uncompressed
end

const JULIA_UUID = UUID("1222c4b2-2114-5bfd-aeef-88e4692bbb3e")
function initialize_uncompressed!(pkg::PkgInfo, versions = keys(pkg.version_info))
    # Only valid to call this with existing versions of the package
    # Remove all versions we have already uncompressed
    versions = filter!(v -> !isinit(pkg.version_info[v], :uncompressed_compat), collect(versions))

    sort!(versions)

    uncompressed_compat = uncompress(pkg.compat, versions)
    uncompressed_deps   = uncompress(pkg.deps,   versions)

    for v in versions
        vinfo = pkg.version_info[v]
        compat = Dict{UUID, VersionSpec}()
        uncompressed_deps_v = uncompressed_deps[v]
        # Everything depends on Julia
        uncompressed_deps_v["julia"] = JULIA_UUID
        uncompressed_compat_v = uncompressed_compat[v]
        for (pkg, uuid) in uncompressed_deps_v
            vspec = get(uncompressed_compat_v, pkg, nothing)
            compat[uuid] = vspec === nothing ? VersionSpec() : vspec
        end
        @init! vinfo.uncompressed_compat = compat
    end
    return pkg
end

function compat_info(pkg::PkgInfo)
    initialize_uncompressed!(pkg)
    return Dict(v => info.uncompressed_compat for (v, info) in pkg.version_info)
end

@lazy struct PkgEntry
    # Registry.toml:
    path::String
    registry_path::String
    name::String
    uuid::UUID

    # Version.toml / (Compat.toml / Deps.toml):
    @lazy info::PkgInfo
end

registry_info(pkg::PkgEntry) = init_package_info!(pkg)

function init_package_info!(pkg::PkgEntry)
    # Already uncompressed the info for this package, return early
    @isinit(pkg.info) && return pkg.info
    path = joinpath(pkg.registry_path, pkg.path)

    path_package = joinpath(path, "Package.toml")
    d_p = parsefile(path_package)
    name = d_p["name"]::String
    name != pkg.name && error("inconsistend name in Registry.toml and Package.toml for pkg at $(path)")
    repo = get(d_p, "repo", nothing)::Union{Nothing, String}
    subdir = get(d_p, "subdir", nothing)::Union{Nothing, String}

    # Versions.toml
    path_vers = joinpath(path, "Versions.toml")
    d_v = isfile(path_vers) ? parsefile(path_vers) : Dict{String, Any}()
    version_info = Dict{VersionNumber, VersionInfo}(VersionNumber(k) =>
        VersionInfo(SHA1(v["git-tree-sha1"]::String), get(v, "yanked", false)::Bool) for (k, v) in d_v)

    # Compat.toml
    compat_file = joinpath(path, "Compat.toml")
    compat_data_toml = isfile(compat_file) ? parsefile(compat_file) : Dict{String, Any}()
    # The Compat.toml file might have string or vector values
    compat_data_toml = convert(Dict{String, Dict{String, Union{String, Vector{String}}}}, compat_data_toml)
    compat = Dict{VersionRange, Dict{String, VersionSpec}}()
    for (v, data) in compat_data_toml
        vr = VersionRange(v)
        d = Dict{String, VersionSpec}(dep => VersionSpec(vr_dep) for (dep, vr_dep) in data)
        compat[vr] = d
    end

    # Deps.toml
    deps_file = joinpath(path, "Deps.toml")
    deps_data_toml = isfile(deps_file) ? parsefile(deps_file) : Dict{String, Any}()
    # But the Deps.toml only have strings as values
    deps_data_toml = convert(Dict{String, Dict{String, String}}, deps_data_toml)
    deps = Dict{VersionRange, Dict{String, UUID}}()
    for (v, data) in deps_data_toml
        vr = VersionRange(v)
        d = Dict{String, UUID}(dep => UUID(uuid) for (dep, uuid) in data)
        deps[vr] = d
    end
    # All packages depend on julia
    deps[VersionRange()] = Dict("julia" => JULIA_UUID)

    @init! pkg.info = PkgInfo(repo, subdir, version_info, compat, deps)

    return pkg.info
end

struct RegistryInstance
    path::String
    name::String
    uuid::UUID
    url::Union{String, Nothing}
    repo::Union{String, Nothing}
    description::Union{String, Nothing}
    pkgs::Dict{UUID, PkgEntry}
    tree_info::Union{Base.SHA1, Nothing}
    # various caches
    name_to_uuids::Dict{String, Vector{UUID}}
end

function RegistryInstance(path::AbstractString; parse_packages::Bool=true)
    d = parsefile(joinpath(path, "Registry.toml"))
    pkgs = Dict{UUID, PkgEntry}()
    if parse_packages
        for (uuid, info) in d["packages"]::Dict{String, Any}
            uuid = UUID(uuid::String)
            info::Dict{String, Any}
            name = info["name"]::String
            pkgpath = info["path"]::String
            pkg = PkgEntry(pkgpath, path, name, uuid, uninit)
            pkgs[uuid] = pkg
        end
    end
    tree_info_file = joinpath(path, ".tree_info.toml")
    tree_info = if isfile(tree_info_file)
        Base.SHA1(parsefile(tree_info_file)["git-tree-sha1"]::String)
    else
        nothing
    end
    return RegistryInstance(
        path,
        d["name"]::String,
        UUID(d["uuid"]::String),
        get(d, "url", nothing)::Union{String, Nothing},
        get(d, "repo", nothing)::Union{String, Nothing},
        get(d, "description", nothing)::Union{String, Nothing},
        pkgs,
        tree_info,
        Dict{String, UUID}(),
    )
end

function Base.show(io::IO, ::MIME"text/plain", r::RegistryInstance)
    println(io, "Registry: $(repr(r.name)) at $(repr(r.path)):")
    println(io, "  uuid: ", r.uuid)
    println(io, "  repo: ", r.repo)
    if r.tree_info !== nothing
        println(io, "  git-tree-sha1: ", r.tree_info)
    end
    println(io, "  packages: ", length(r.pkgs))
end

function uuids_from_name(r::RegistryInstance, name::String)
    create_name_uuid_mapping!(r)
    return get(Vector{UUID}, r.name_to_uuids, name)
end

function create_name_uuid_mapping!(r::RegistryInstance)
    isempty(r.name_to_uuids) || return
    for (uuid, pkg) in r.pkgs
        uuids = get!(Vector{UUID}, r.name_to_uuids, pkg.name)
        push!(uuids, pkg.uuid)
    end
    return
end

function reachable_registries(; depots::Union{String, Vector{String}}=Base.DEPOT_PATH, parse_packages::Bool=true)
    # collect registries
    if depots isa String
        depots = [depots]
    end
    registries = RegistryInstance[]
    for d in depots
        isdir(d) || continue
        reg_dir = joinpath(d, "registries")
        isdir(reg_dir) || continue
        for name in readdir(reg_dir)
            file = joinpath(reg_dir, name, "Registry.toml")
            isfile(file) || continue
            push!(registries, RegistryInstance(joinpath(reg_dir, name); parse_packages))
        end
    end
    return registries
end

# Dict interface

function Base.haskey(r::RegistryInstance, uuid::UUID)
    return haskey(r.pkgs, uuid)
end

function Base.keys(r::RegistryInstance)
    return keys(r.pkgs)
end

Base.getindex(r::RegistryInstance, uuid::UUID) = r.pkgs[uuid]

Base.get(r::RegistryInstance, uuid::UUID, default) = get(r.pkgs, uuid, default)

Base.iterate(r::RegistryInstance) = iterate(r.pkgs)
Base.iterate(r::RegistryInstance, state) = iterate(r.pkgs, state)


text:
using Pkg
Pkg.activate(joinpath(@__DIR__))
Pkg.instantiate()
include("advanced_cli.jl")

max_stars_repo_name: Novirmor/MocosSimLauncher

max_issues_repo_name: Novirmor/MocosSimLauncher

max_forks_repo_name: Novirmor/MocosSimLauncher


text:
kernelM(l::Int,x, y) = x * diffsbesselj(l,x) * sbesselj(l,y) - y * sbesselj(l,x) * diffsbesselj(l,y)

function scattering_field(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},WithoutSymmetry{3}}) where T
    L = wavemode.basis_order
    L1 = wavemode.basis_field_order

    i2s = 1:size(wavemode.eigenvectors,2)
    i3s = 1:size(wavemode.eigenvectors,3)

    nn1_to_n = [
            basisorder_to_basislength(Acoustic{T,3},l-1) + m + l + 1
    for l = 0:L for m = -l:l for l1 = 0:L1 for m1 = -l1:l1];

    nn1_to_n1 = [
            basisorder_to_basislength(Acoustic{T,3},l1-1) + m1 + l1 + 1
    for l = 0:L for m = -l:l for l1 = 0:L1 for m1 = -l1:l1];

    indices_same_n = collect(groupby(nn1 -> nn1_to_n[nn1] , 1:length(nn1_to_n)))

    v = regular_basis_function(wavemode.wavenumber, wavemode.medium)

    function scat_field(x::AbstractVector{T})
        vs =  v(L1, x)

        vecs = [
            sum(wavemode.eigenvectors[ind,i2,i3] .* vs[nn1_to_n1[ind]]) # note here:  vs[nn1_to_n1[ind]]] == vs
        for ind in indices_same_n, i2 in i2s, i3 in i3s];

        # Sum over species and different eigvectors
        vecs = sum(vecs, dims = 3)[:,:];

        # for only 1 species the code below gives the same
        # nn1_to_n1 = [
        #         basisorder_to_basislength(Acoustic{T,3},l1-1) + m1 + l1 + 1
        # for l = 0:L for m = -l:l for l1 = 0:L1 for m1 = -l1:l1];
        #
        # vecs2 = [
        #     wavemode.eigenvectors[i[1],i[2],i[3]] * vs[nn1_to_n1[i[1]]]
        # for i in CartesianIndices(wavemode.eigenvectors)];
        #
        # vecs2 = sum(vecs2, dims = 3);
        # vecs2 = [sum(vecs2[inds]) for inds in indices_same_n]
        # vecs ≈ vecs2
        return vecs
    end

end

function scattering_field(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},Sym}) where {T, Sym<:AbstractAzimuthalSymmetry{3}}
    L = wavemode.basis_order
    L1 = wavemode.basis_field_order

    i2s = 1:size(wavemode.eigenvectors,2)
    i3s = 1:size(wavemode.eigenvectors,3)

    nl1_to_n = [
            basisorder_to_basislength(Acoustic{T,3},l-1) + m + l + 1
    for l = 0:L for m = -l:l for l1 = abs(m):L1];

    nl1_to_n1 = [
            basisorder_to_basislength(Acoustic{T,3},l1-1) - m + l1 + 1
    for l = 0:L for m = -l:l for l1 = abs(m):L1];

    indices_same_n = collect(groupby(nl1 -> nl1_to_n[nl1] , 1:length(nl1_to_n)))

    v = regular_basis_function(wavemode.wavenumber, wavemode.medium)

    function scat_field(x::AbstractVector{T})
        vs =  v(L1, x)

        vecs = [
            sum(wavemode.eigenvectors[ind,i2,i3] .* vs[nl1_to_n1[ind]])
        for ind in indices_same_n, i2 in i2s, i3 in i3s];

        # Sum over different eigvectors
        vecs = sum(vecs, dims = 3)[:,:];

        return vecs
    end

end

function material_scattering_coefficients(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},WithoutSymmetry{3}}) where T

    # Unpacking parameters
    k = wavemode.ω / wavemode.medium.c

    species = wavemode.material.species
    S = length(species)
    rs = outer_radius.(species)

    R = outer_radius(wavemode.material.shape)
    k_eff = wavemode.wavenumber

    L = wavemode.basis_order
    L1 = wavemode.basis_field_order
    dL = L + L1

    # Pre-calculating
    Ml1s = [
        (R - rs[j]) * kernelM(l1, k*(R - rs[j]), k_eff*(R - rs[j])) * number_density(species[j])
    for l1 = 0:L1, j in eachindex(species)] ./ (k_eff^T(2) - k^T(2))

    nn1_to_l1 = [l1 for l = 0:L for m = -l:l for l1 = 0:L1 for m1 = -l1:l1];

    # Sum over species and different eigvectors
    vecs = [
        wavemode.eigenvectors[i] * Ml1s[nn1_to_l1[i[1]]+1,i[2]]
    for i in CartesianIndices(wavemode.eigenvectors)]

    vec = reshape(sum(vecs, dims = (2,3)),:)

    # Sum over wave basis components
    gaunt = [
            gaunt_coefficient(dl,dm,l,m,l1,m1)
        for l = 0:L for m = -l:l
        for l1 = 0:L1 for m1 = -l1:l1
    for dl = 0:dL for dm = -dl:dl]

    # gaunt[(dn),(n,l1)] = gaunt_coefficient(dl,dm,l,m,l1,m1)
    gaunt = reshape(gaunt, (:,length(vec)))

    return gaunt * vec
end

function material_scattering_coefficients(wavemode::EffectiveRegularWaveMode{T,3,Acoustic{T,3},Sym}) where {T,Sym<:AbstractAzimuthalSymmetry{3}}

    # Unpacking parameters
    k = wavemode.ω / wavemode.medium.c

    species = wavemode.material.species
    S = length(species)
    rs = outer_radius.(species)

    R = outer_radius(wavemode.material.shape)
    k_eff = wavemode.wavenumber

    L = wavemode.basis_order
    L1 = wavemode.basis_field_order
    dL = L + L1

    # Pre-calculating
    Ml1s = [
        (R - rs[j]) * kernelM(l1, k*(R - rs[j]), k_eff*(R - rs[j])) * number_density(species[j])
    for l1 = 0:L1, j in eachindex(species)] ./ (k_eff^T(2) - k^T(2))

    nl1_to_l1 = [l1 for l = 0:L for m = -l:l for l1 = abs(m):L1];

    # Sum over species and different eigvectors
    vecs = [
        wavemode.eigenvectors[i] * Ml1s[nl1_to_l1[i[1]]+1,i[2]]
    for i in CartesianIndices(wavemode.eigenvectors)]

    vec = reshape(sum(vecs, dims = (2,3)),:)

    # Sum over wave basis components
    gaunt = [
            gaunt_coefficient(dl,dm,l,m,l1,-m)
        for l = 0:L for m = -l:l
        for l1 = abs(m):L1
    for dl = 0:dL for dm = -dl:dl]

    gaunt = reshape(gaunt, (:,length(vec)))

    return gaunt * vec
end


text:
#include("../get_home_dir.jl")

using DataFrames
using CSV
using PyPlot
using Random
using StatsBase
using Statistics

#include(string(PROJECT_HOME,"/assim_code/time_averaging.jl"))

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams");
rcParams["lines.linewidth"] = 1;
rcParams["font.size"] = 24;

pmean = CSV.read("dlwp3/DL_means2.csv",DataFrame, header=0);
p25 = CSV.read("dlwp3/DL_lower2.csv",DataFrame,header=0);
p75 = CSV.read("dlwp3/DL_upper2.csv",DataFrame,header=0);

ptrue = CSV.read("dlwp3/DL_true2.csv",DataFrame,header=0);

colors_list = ["tab:blue","tab:green","tab:orange","tab:red"];
label_list = ["All obs.", "1 AM/PM obs.", "6 AM/PM obs."];

year_x = collect(1:24);

#=
figure()
plot(year_x, ptrue[:,1],color="k",linewidth=3, label="Truth")
for i in 1:3
	plot(year_x, pmean[:,i],color=colors_list[i],linewidth=3, label=label_list[i])
end
for i in 1:3
	fill_between(year_x, p25[:,i], p75[:,i] , color=colors_list[i],alpha=0.33)
end
=#

figure()
plot(year_x, (ptrue[:,1] .- mean(ptrue[:,1])) / std(ptrue[:,1]),color="k",linewidth=3, label="Truth")
for i in 1:3
	plot(year_x, (pmean[:,i] .- mean(pmean[:,i])) / std(pmean[:,i]),color=colors_list[i],linewidth=3, label=label_list[i])
end

legend()
xlabel("Time of day")
ylabel("Leaf water potential, diurnal mean (MPa)")

#need to start at 0.3 soil moisture in year 3


text:
mutable struct Grid <: AbstractGrid
  player::AbstractPlayer

  rows::Int
  cols::Int

  hidden_rows::Int
  height::Int

  table::Array{AbstractString,2}
end

function Grid(cur_player::AbstractPlayer=Player())
  cur_rows = 20
  cur_cols = 10

  cur_hidden_rows = 2
  cur_height = 0

  cur_grid = Grid(
    cur_player,
    cur_rows,
    cur_cols,
    cur_hidden_rows,
    cur_height,
    fill!(
      Array{AbstractString}(
        cur_rows + cur_hidden_rows,
        cur_cols
      ),
      ""
    )
  )

  cur_player.grid = cur_grid

  cur_grid
end


text:
@testset "solve_burgers.jl" begin
    @test_broken solve()
end