text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -6.019960436071491 <= q <= 11.696817216382863)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.9042935565441375 + 0.31236426273908524*1/(1 + exp(-(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4])))) + -0.4128166127064836*1/(1 + exp(-(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4])))) + -0.2476169775927386*1/(1 + exp(-(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))))))) + 1/(1 + exp(-(0.8600467319105785 + -0.4006275013725942*1/(1 + exp(-(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4])))) + 0.22106887260345687*1/(1 + exp(-(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4])))) + 0.7210213948433983*1/(1 + exp(-(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))))))) + 1/(1 + exp(-(0.982756228863964 + 0.34792310275142935*1/(1 + exp(-(0.5868252904954097 + 0.04712499510197299*$(x[1]) + 0.4005456947892885*$(x[2]) + 0.34544434364772814*$(x[3]) + 0.48957180414623247*$(x[4])))) + 0.7924515917531365*1/(1 + exp(-(0.5797949667252156 + -0.7245574143823843*$(x[1]) + -0.08207866517990237*$(x[2]) + 0.4631961331947827*$(x[3]) + 0.9179781125579716*$(x[4])))) + -0.9399735455074207*1/(1 + exp(-(0.8773151656273224 + 0.003673243009740812*$(x[1]) + 0.698013040735662*$(x[2]) + 0.6787695040522581*$(x[3]) + -0.9128403322093082*$(x[4]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#:D why not....
nth = [:second, :third, :fourth, :fifth, :sixth, :seventh, :eigth, :ninth, :tenth,
    :eleventh, :twelfth, :thirteenth, :fourteenth, :fifteenth, :sixteenth, :seventeenth,
    :eighteenth, :nineteenth, :twentieth]

for (i,n) in enumerate(nth)
    @eval $n(x) = x[$(i+1)]
end

ext: jl


text:
using OrdinalGWAS, Test, CSV, SnpArrays, DataFrames

const datadir = joinpath(dirname(@__FILE__), "..", "data")
const covfile = datadir * "/covariate.txt"
const plkfile = datadir * "/hapmap3"
const snpsetfile = datadir * "/hapmap_snpsetfile.txt"
const vcfcovfile = datadir * "/vcf_example.csv"
const vcffile = datadir * "/vcf_test"
const vcfsnpsetfile = datadir * "/snpsetfile_vcf.txt"
const bgencovfile = datadir * "/bgen_ex.csv"
const bgenfile = datadir * "/bgen_test"
const bgensnpsetfile = datadir * "/bgen_snpsetfile.txt"

@testset "score test" begin
    @time ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:score)
    @test isfile("ordinalgwas.null.txt")
    @test isfile("ordinalgwas.pval.txt")
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, 6][1:5]
    @test isapprox(scorepvals, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-4)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    # VCF
    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile; geneticformat = "VCF",
        vcftype = :DS, geneticrowinds = 1:190, snpinds = [86; 656],
    test = :score, covrowinds = 1:190)
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:2]
    @test isapprox(scorepvals, [0.00762272, 0.000668338], rtol=1e-4)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    # BGEN
    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        test = :score)
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:2]
    @test isapprox(scorepvals, [0.12449778, 0.00055727], rtol=1e-4)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)
end

@testset "LRT test" begin
    @time ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:LRT)
    @test isfile("ordinalgwas.null.txt")
    @test isfile("ordinalgwas.pval.txt")
    lrtpvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(lrtpvals, [1.0, 1.91858366e-3, 1.80505056e-5, 5.87338471e-6, 8.08102258e-3], rtol=1e-4)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    # VCF
    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile; geneticformat = "VCF",
        vcftype = :GT, snpinds = [86; 656], test = :LRT)
    lrtpvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:2]
    @test isapprox(lrtpvals, [0.00955468405473856, 0.0007086063489553798], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    # BGEN
    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        snpinds = [3, 25], test = :LRT)
        lrtpvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:2]
    @test isapprox(lrtpvals, [0.000445619, 1.660631e-6], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)
end

@testset "snpmodel" begin
    # dominant model
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:score, snpmodel=DOMINANT_MODEL)
    @test isfile("ordinalgwas.null.txt")
    @test isfile("ordinalgwas.pval.txt")
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [1.0, 0.14295, 0.000471942, 0.00555348, 0.000652844], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)
    # recessive model
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:score, snpmodel=RECESSIVE_MODEL)
    @test isfile("ordinalgwas.null.txt")
    @test isfile("ordinalgwas.pval.txt")
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [1.0, 0.00673612, 0.000279908, 4.15322e-5, 0.167642], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)
end

@testset "link" begin
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, link=ProbitLink(), pvalfile="opm.pval.txt")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("opm.pval.txt")
    scorepvals = CSV.read("opm.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [1.0, 1.00769167e-2, 2.62725649e-5, 1.08974849e-5, 5.10288399e-3], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("opm.pval.txt", force=true)
end

@testset "snp mask" begin
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, snpinds=1:5, pvalfile="first5snps.pval.txt")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("first5snps.pval.txt")
    @test countlines("first5snps.pval.txt") == 6
    scorepvals = CSV.read("first5snps.pval.txt", DataFrame)[!, end]
    @test isapprox(scorepvals, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("first5snps.pval.txt", force=true)
end

@testset "sub samples" begin
    # only use first 300 samples
    @time ordinalgwas(@formula(trait ~ sex), covfile, plkfile, test=:score, covrowinds=1:300, geneticrowinds=1:300)
    @test isfile("ordinalgwas.null.txt")
    @test isfile("ordinalgwas.pval.txt")
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [1.0, 0.00355969, 0.000123604, 5.2213e-6, 0.00758234], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)
end

@testset "test formula" begin
    # score test
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, pvalfile="GxE.pval.txt",
        testformula=@formula(trait ~ snp + snp & sex))
    @test isfile("ordinalgwas.null.txt")
    @test isfile("GxE.pval.txt")
    scorepvals = CSV.read("GxE.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [1.0, 1.74460104e-2, 1.66707324e-4, 4.76376246e-5, 2.91384712e-2], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("GxE.pval.txt", force=true)
    # LRT, only first 5 SNPs
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, pvalfile="GxE.pval.txt",
        testformula=@formula(trait ~ snp + snp & sex), test=:LRT, snpinds=1:5)
    @test isfile("ordinalgwas.null.txt")
    @test isfile("GxE.pval.txt")
    lrtpvals = CSV.read("GxE.pval.txt", DataFrame)[!, end]
    @test isapprox(lrtpvals, [1.0, 7.22410973e-3, 1.01730983e-4, 1.88174211e-5, 2.88295705e-2], rtol=1e-4)
    rm("ordinalgwas.null.txt", force=true)
    rm("GxE.pval.txt", force=true)

    # BGEN
    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        snpinds = [3, 25], test = :LRT, testformula = @formula(trait ~ snp + snp & sex))
    lrtpvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:2]
    @test isapprox(lrtpvals, [0.002084, 8.135205e-6], rtol = 1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)
end

@testset "snpset" begin
    #window
    #score test
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, pvalfile = "snpset.pval.txt",
        snpset=250, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    scorepvals = CSV.read("snpset.pval.txt", DataFrame)[!, end][1:5]
    #@test isapprox(scorepvals, [1.0, 1.74460104e-2, 1.66707324e-4, 4.76376246e-5, 2.91384712e-2], rtol=1e-4)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)
    #lrt
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, pvalfile = "snpset.pval.txt",
        snpset=25, test=:LRT, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    lrtpvals = CSV.read("snpset.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(lrtpvals, [2.1817554071810948e-13, 0.2865769729670889, 0.32507802233937966,
    0.3344823237332578, 0.42948375949508427], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)

    # VCF
    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile; geneticformat = "VCF",
        vcftype = :DS, pvalfile = "snpset.pval.txt",
    snpset=250, test=:score, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    scorepvals = CSV.read("snpset.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [0.4278366599084349,
    0.42781616067453476, 0.4519573757701432,
    0.4278763804444088, 0.4345883185481474], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)

    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile; geneticformat = "VCF",
        vcftype = :DS, pvalfile = "snpset.pval.txt",
    snpset=25, test=:LRT, analysistype = "snpset")
    lrtpvals = CSV.read("snpset.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(lrtpvals, [1.0
    0.9999996378252629
    1.0
    0.9999999999996334
    0.9999999976994737], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)

    # BGEN
    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        test = :score, snpset = 8, analysistype = "snpset")
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:2]
    @test isapprox(scorepvals, [0.00767237, 0.5512827], rtol = 1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        test = :LRT, snpset = 8, analysistype = "snpset")
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:2]
    @test isapprox(scorepvals, [0.05257402, 0.5706598], rtol = 1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    #snpset file
    #score test
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, pvalfile = "snpset.pval.txt",
        snpset = snpsetfile, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    scorepvals = CSV.read("snpset.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [1.72134e-5, 0.036925, 0.747855,
     0.0276508, 0.611958], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)
    #lrt
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, pvalfile = "snpset.pval.txt",
        snpset = snpsetfile, test = :lrt, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    lrtpvals = CSV.read("snpset.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(lrtpvals, [6.75377e-13, 0.000256566, 0.359382,
     0.000163268, 0.0867508], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)

    # VCF
    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile; geneticformat = "VCF",
        vcftype = :DS, pvalfile = "snpset.pval.txt",
        snpset = vcfsnpsetfile, test = :score, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    scorepvals = CSV.read("snpset.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [0.06814002639277685, 0.5566664123188036,
    0.520381855174413, 0.07557764137466122, 0.5620803022597403], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)

    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile; geneticformat = "VCF",
        vcftype = :DS, pvalfile = "snpset.pval.txt",
        snpset = vcfsnpsetfile, test = :lrt, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    lrtpvals = CSV.read("snpset.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(lrtpvals, [0.09069975735216675, 0.6465153355309161,
    0.6307411986741357, 0.06275888993714969, 0.50252192003468], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)

    # BGEN
    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        test = :score, snpset = bgensnpsetfile, analysistype = "snpset")
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:2]
    @test isapprox(scorepvals, [0.013317595, 0.454769036], rtol = 1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    # BGEN
    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        test = :lrt, snpset = bgensnpsetfile, analysistype = "snpset")
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:2]
    @test isapprox(scorepvals, [0.08986309, 0.4688198], rtol = 1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    #specific snp (one snpset)
    #score test
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, pvalfile = "snpset.pval.txt",
        snpset = 50:55, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    scorepvals = open("snpset.pval.txt")
    scorepval = split(readline(scorepvals))[end]
    close(scorepvals)
    @test isapprox(parse(Float64, scorepval), 0.3647126536663949, rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)
    #lrt
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, pvalfile = "snpset.pval.txt",
        snpset = collect(1:15), test=:LRT, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    lrtpvals = open("snpset.pval.txt")
    lrtpval = split(readline(lrtpvals))[end]
    close(lrtpvals)
    @test isapprox(parse(Float64, lrtpval), 7.525696044086955e-15, rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)

    # VCF
    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile; geneticformat = "VCF",
        vcftype = :DS, pvalfile = "snpset.pval.txt",
        snpset=85:90, test=:score, analysistype = "snpset")
    @test isfile("ordinalgwas.null.txt")
    @test isfile("snpset.pval.txt")
    scorepvals = open("snpset.pval.txt")
    scorepval = split(readline(scorepvals))[end]
    close(scorepvals)
    @test isapprox(parse(Float64, scorepval), 0.0965927460813927, rtol=1e-2)

    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile; geneticformat = "VCF",
        vcftype = :DS, pvalfile = "snpset.pval.txt",
        snpset=85:90, test=:lrt, analysistype = "snpset")
    lrtpvals = open("snpset.pval.txt")
    lrtpval = split(readline(lrtpvals))[end]
    close(lrtpvals)
    @test isapprox(parse(Float64, lrtpval), 0.0732485446883825, rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("snpset.pval.txt", force=true)

    # BGEN
    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        test = :score, snpset = [3, 25, 100], analysistype = "snpset")
    scorepvals = open("ordinalgwas.pval.txt")
    scorepval = split(readline(scorepvals))[end]
    close(scorepvals)
    @test isapprox(parse(Float64, scorepval), 4.422683e-8, rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        test = :LRT, snpset = [3, 25, 100], analysistype = "snpset")
    scorepvals = open("ordinalgwas.pval.txt")
    scorepval = split(readline(scorepvals))[end]
    close(scorepvals)
    @test isapprox(parse(Float64, scorepval), 2.853246e-8, rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)
end

@testset "GxE snp in null" begin
    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, e = :sex, pvalfile = "gxe_snp.pval.txt",
        snpinds=1:5, test=:score, analysistype = "gxe")
    @test isfile("gxe_snp.pval.txt")
    scorepvals = CSV.read("gxe_snp.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [1.0, 0.637742242597749, 0.9667114198051628,
    0.26352674694121003, 0.7811133315582837], rtol=1e-2)
    rm("gxe_snp.pval.txt", force=true)

    ordinalgwas(@formula(trait ~ sex), covfile, plkfile, e = "sex", pvalfile = "gxe_snp.pval.txt",
        snpinds=1:5, test=:LRT, analysistype = "gxe")
    @test isfile("gxe_snp.pval.txt")
    lrtpvals = CSV.read("gxe_snp.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(lrtpvals, [1.0, 0.6279730133445315, 0.9671662821946985,
    0.26693502209463904, 0.7810214899265426], rtol=1e-2)
    rm("gxe_snp.pval.txt", force=true)

    # VCF
    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile, e = :sex; geneticformat = "VCF",
        vcftype = :DS, pvalfile = "gxe_snp.pval.txt",
        snpinds=1:5, test=:score, analysistype = "gxe")
    scorepvals = CSV.read("gxe_snp.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [0.45861769035708144, 1.0, 0.03804677528312195,
     0.18254151103030725, 0.34454453512541156], rtol=1e-2)
    rm("gxe_snp.pval.txt", force=true)

    ordinalgwas(@formula(y ~ sex), vcfcovfile, vcffile, e = :sex; geneticformat = "VCF",
        vcftype = :DS, pvalfile = "gxe_snp.pval.txt",
        snpinds=1:5, test=:lrt, analysistype = "gxe")
    @test isfile("gxe_snp.pval.txt")
    lrtpvals = CSV.read("gxe_snp.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(lrtpvals, [0.526667096902957, 1.0, 0.008073040021982156,
    0.10590569987122991, 0.3557829099471382], rtol=1e-2)
    rm("gxe_snp.pval.txt", force=true)

    # BGEN
    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        test = :score, snpinds = 1:5, analysistype = "gxe", e = :sex)
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [0.415677
        0.92019145
        0.8975205
        0.4947529
        0.4947529], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)

    ordinalgwas(@formula(y ~ sex), bgencovfile, bgenfile; geneticformat = "BGEN",
        test = :lrt, snpinds = 1:5, analysistype = "gxe", e = :sex)
    scorepvals = CSV.read("ordinalgwas.pval.txt", DataFrame)[!, end][1:5]
    @test isapprox(scorepvals, [0.419214
        0.9207400
        0.8981011
        0.4947914
        0.4947914], rtol=1e-2)
    rm("ordinalgwas.null.txt", force=true)
    rm("ordinalgwas.pval.txt", force=true)
end

@testset "split, gz" begin
    # split hapmap3 by chromosome
    SnpArrays.split_plink(plkfile, :chromosome; prefix = datadir * "/hapmap3.chr.")
    # compress to gz
    for chr in 1:23
        plinkfile = plkfile * ".chr." * string(chr)
        SnpArrays.compress_plink(plinkfile)
        @test isfile(plinkfile * ".bed.gz")
        @test isfile(plinkfile * ".fam.gz")
        @test isfile(plinkfile * ".bim.gz")
    end
    # fit null model
    @time nm = ordinalgwas(@formula(trait ~ sex), covfile, nothing)
    @test isfile("ordinalgwas.null.txt")
    # gwas by chromosome, refit null model each time, use uncompressed Plink set
    @time for chr in 1:23
        plinkfile = plkfile * ".chr." * string(chr)
        pvalfile = plkfile * ".chr." * string(chr) * ".pval.txt"
        ordinalgwas(@formula(trait ~ sex), covfile, plinkfile, pvalfile = pvalfile)
        @test isfile(pvalfile)
        if chr == 1
            pvals_chr1 = CSV.read(pvalfile, DataFrame)[!, end][1:5]
            @test isapprox(pvals_chr1, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-2)
        end
        rm(plinkfile * ".pval.txt", force=true)
    end
    # gwas by chromosome, use fitted null model each time, use uncompressed Plink set
    @time for chr in 1:23
        plinkfile = plkfile * ".chr." * string(chr)
        pvalfile = plkfile * ".chr." * string(chr) * ".pval.txt"
        ordinalgwas(nm, plinkfile, pvalfile = pvalfile)
        @test isfile(pvalfile)
        if chr == 1
            pvals_chr1 = CSV.read(pvalfile, DataFrame)[!, end][1:5]
            @test isapprox(pvals_chr1, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-2)
        end
        rm(pvalfile, force=true)
    end
    # gwas by chromosome, use fitted null model each time, use compressed bed and bim files
    @time for chr in 1:23
        bedfile = plkfile * ".chr." * string(chr) * ".bed.gz"
        bimfile = plkfile * ".chr." * string(chr) * ".bim.gz"
        pvalfile = plkfile * ".chr." * string(chr) * ".pval.txt"
        ordinalgwas(nm, bedfile, bimfile, 324; pvalfile = pvalfile)
        @test isfile(pvalfile)
        if chr == 1
            pvals_chr1 = CSV.read(pvalfile, DataFrame)[!, end][1:5]
            @test isapprox(pvals_chr1, [1.0, 4.56531284e-3, 3.10828383e-5, 1.21686724e-5, 8.20686005e-3], rtol=1e-2)
        end
        rm(pvalfile, force=true)
    end
    # clean up
    # delete result files
    isfile("ordinalgwas.null.txt") && rm("ordinalgwas.null.txt")
    for chr in 1:26
        plinkfile = plkfile * ".chr." * string(chr)
        # delete uncompressed chromosome Plink files
        rm(plinkfile * ".bed", force=true)
        rm(plinkfile * ".fam", force=true)
        rm(plinkfile * ".bim", force=true)
        # delete compressed chromosome Plink files
        rm(plinkfile * ".bed.gz", force=true)
        rm(plinkfile * ".fam.gz", force=true)
        rm(plinkfile * ".bim.gz", force=true)
        # delete pval files
        rm(plinkfile * ".pval.txt", force=true)
    end
end

max_stars_repo_name: Hua-Zhou/PolrGWAS

max_issues_repo_name: Hua-Zhou/PolrGWAS


text:
@testset "Test domain distributions" begin
    @test size(ln_real) == size(Ω)
    @test size(mvln_real) == size(Ω)

    err_msg = "Passed MidpointDisplacement as a type, need an instance e.g. MidpointDisplacement()"
    @test_throws ErrorException(err_msg) DomainDistribution(MidpointDisplacement, Ω)

    @test size(rand(dd1)) == size(Ω)
    @test size(rand(cdd1)) == size(Ω)
    @test size(rand(cdd2)) == size(Ω)
    @test length(rand(mdd1)) == 2
    @test size(rand(bdd2)) == size(Ω)
    @test size(rand(bdd3)) == size(Ω)
end


text:
include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/common_functions.jl")
include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/FK_sampling_functions.jl")
include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/SB_sampling_functions.jl")

using JLD, GibbsTypePriors, DataFramesMeta, RCall, Optim, OptimTestProblems, StatsBase

function EV_of_number_of_clusters_NGG(n,H,par_vec)
    pk_ngg = GibbsTypePriors.Pkn_NGG.(1:n, n, par_vec[1], par_vec[2])
    E_ngg =  pk_ngg|> ar -> map(*, ar, 1:n) |> sum
    x = ((1:n).-E_ngg).^2
    V_ngg = pk_ngg |> ar -> map(*, ar, x) |> sum
    pk_nggm = Pkn_NGGM_precomp.(1:n,n,H,par_vec[1], par_vec[2], [pk_ngg])
    E_nggm =  pk_nggm|> ar -> map(*, ar, 1:n) |> sum
    x_nggm = ((1:n).-E_nggm).^2
    V_nggm = pk_nggm |> ar -> map(*, ar, x_nggm) |> sum
    return [E_ngg, V_ngg] , [E_nggm, V_nggm]
end

function EV_of_number_of_clusters_NGG_approx(n,par_vec)
    pk_ngg_approx = Pkn_NGG_pred_approx(n,par_vec[1], par_vec[2])
    E_ngg_approx =  pk_ngg_approx|> ar -> map(*, ar, 1:n) |> sum
    x = ((1:n).-E_ngg_approx).^2
    V_ngg_approx = pk_ngg_approx |> ar -> map(*, ar, x) |> sum
    return E_ngg_approx, V_ngg_approx
end

function EV_of_number_of_clusters_NGG_FK(n,par_vec,M)
    pk_ngg_fk = Pkn_NGG_FK_fast(n,par_vec[1], par_vec[2],M)
    E_ngg_fk =  pk_ngg_fk|> ar -> map(*, ar, 1:n) |> sum
    x = ((1:n).-E_ngg_fk).^2
    V_ngg_fk = pk_ngg_fk |> ar -> map(*, ar, x) |> sum
    return E_ngg_fk, V_ngg_fk
end

function EV_of_number_of_clusters_NGG_FK_slow(n,par_vec,M)
    println([par_vec[2],par_vec[1]])
    pk_ngg_fk = Pkn_NGG_FK(n,par_vec[1], par_vec[2], 250; runs = 2*10^2)
    E_ngg_fk =  pk_ngg_fk|> ar -> map(*, ar, 1:n) |> sum
    x = ((1:n).-E_ngg_fk).^2
    V_ngg_fk = pk_ngg_fk |> ar -> map(*, ar, x) |> sum
    return E_ngg_fk, V_ngg_fk
end

function EV_of_number_of_clusters_NGG_FK_fast(n,par_vec,M)
    println([par_vec[2],par_vec[1]])
    pk_ngg_fk = Pkn_NGG_FK_fast(n,par_vec[1], par_vec[2], 250; runs = 2*10^2)
    E_ngg_fk =  pk_ngg_fk|> ar -> map(*, ar, 1:n) |> sum
    x = ((1:n).-E_ngg_fk).^2
    V_ngg_fk = pk_ngg_fk |> ar -> map(*, ar, x) |> sum
    return E_ngg_fk, V_ngg_fk
end

sigma = collect(range(0.05,0.99, length=10))
alpha = collect(exp.(range(log(1), log(200), length =10)))
grid = collect(Iterators.product(alpha, sigma))
n=100
grid_vec = vec(grid)
H = 250

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_6.jld", "grid_vec" ,grid_vec)

EV = EV_of_number_of_clusters_NGG.(n,H,grid_vec)
EV_NGG= first.(EV)
EV_NGGM = last.(EV)

## Expectation and variance for NGG
save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_1.jld", "NGG", EV_NGG)

## Expectation and variance for NGG multinomial
save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_3.jld", "NGGM", EV_NGGM)

## Expectation and variance for NGG predictive approximation

EV_NGG_approx = EV_of_number_of_clusters_NGG_approx.(n,grid_vec)

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_2.jld", "NGG_approx", EV_NGG_approx)

function EV_of_number_of_clusters_NGG_SB(par_vec,n,M)
    println([par_vec[1], par_vec[2]])
    pk_ngg_sb= Pkn_NGG_SB_1(n,par_vec[1], par_vec[2],M)
    E_ngg_sb =  pk_ngg_sb|> ar -> map(*, ar, 1:n) |> sum
    x = ((1:n).-E_ngg_sb).^2
    V_ngg_sb = pk_ngg_sb |> ar -> map(*, ar, x) |> sum
    return E_ngg_sb, V_ngg_sb
end

sigma = collect(range(0.05,0.99, length=10))
#alpha = collect(range(0.05,20, length=5))
alpha = collect(exp.(range(log(1), log(200), length =10)))
grid = collect(Iterators.product(alpha, sigma))
#grid_borders = vcat(collect(Iterators.product(alpha[1], sigma)),collect(Iterators.product(alpha[nb], sigma)), collect(Iterators.product(alpha, sigma[1])),collect(Iterators.product(alpha, sigma[ns])))
n=100
grid_vec = vec(grid)[11:100]
H = 250

EV_NGG_SB = EV_of_number_of_clusters_NGG_SB.(grid_vec,n,H)

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_5.jld", "NGG_FSB", EV_NGG_SB)

## Expectation and variance for NGG FK approximation
## As FK is computationally expensive algorithm we use the reduced grid for computing Expectation and Variance
reduced_grid = filter(x -> (x[2] > 0.25)&&(x[1] > 3), grid_vec)

EV_NGG_FK_ = EV_of_number_of_clusters_NGG_FK_fast.(n,reduced_grid,250)
EV_NGG_FK__fast = DataFrame(EV_NGG_FK_)
DF_NGG_FK_fast = DataFrame(Exp = EV_NGG_FK__fast[1],Var =EV_NGG_FK__fast[2], beta = first.(reduced_grid), sigma = last.(reduced_grid))
DF_NGG_FK_fast.std = sqrt.(DF_NGG_FK_fast.Var)
DF_NGG_FK_fast.beta_log= log.(DF_NGG_FK_fast.beta)

DF_NGG_FK_reduced = EV_NGG_FK__fast[9:64,:]
reduced_grid_red =  reduced_grid[9:64]

DF_NGG_FK = DataFrame(Exp = DF_NGG_FK_reduced[1],Var =DF_NGG_FK_reduced[2], beta = first.(reduced_grid_red), sigma = last.(reduced_grid_red))
DF_NGG_FK.std = sqrt.(DF_NGG_FK.Var)
DF_NGG_FK.beta_log= log.(DF_NGG_FK.beta)

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_4.jld", "DF_NGG_FK", DF_NGG_FK)

### more tine consuming computatoin for FK (using reduced grid and coombination of fast&slow algorithms)

reduced_grid = filter(x -> (x[2] > 0.25)&&(x[1] > 3), grid_vec)
EV_NGG_FK_slow = EV_of_number_of_clusters_NGG_FK_slow.(n,reduced_grid,250)

EV_NGG_FK__slow = DataFrame(EV_NGG_FK_slow)
EV_NGG_FK__slow = DataFrame(Exp = EV_NGG_FK__slow[1],Var =EV_NGG_FK__slow[2], beta = first.(reduced_grid), sigma = last.(reduced_grid))
EV_NGG_FK__slow.std = sqrt.(EV_NGG_FK__slow.Var)
EV_NGG_FK__slow.beta_log= log.(EV_NGG_FK__slow.beta)

DF_NGG_FK_reduced_slow = EV_NGG_FK__slow[9:64,:]
reduced_grid_red =  reduced_grid[9:64]

DF_NGG_FK_slow = DataFrame(Exp = DF_NGG_FK_reduced_slow[1],Var =DF_NGG_FK_reduced_slow[2], beta = first.(reduced_grid_red), sigma = last.(reduced_grid_red))
DF_NGG_FK_slow.std = sqrt.(DF_NGG_FK_slow.Var)
DF_NGG_FK_slow.beta_log= log.(DF_NGG_FK_slow.beta)

save("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_6_4_2.jld", "DF_NGG_FK_slow", DF_NGG_FK_slow)

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
"""
This module is an example of how to create and import a submodule.

Functions can be exported using `export` and are then usable when the parent
module has a line like `using .ExampleSubmodule`.
"""
module ExampleSubmodule

export anexamplefunction

function anexamplefunction()
  print("some information I suppose")
end

end # end ExampleSubmodule

max_stars_repo_name: BradLyman/learn-julia

max_issues_repo_name: BradLyman/learn-julia

max_forks_repo_name: BradLyman/learn-julia


text:
#
# Provides concrete types under `ParameterizedLabelledArray`
# which contain parameters to describe astrodynamic systems.
#

"""
A supertype for parameter representations in astrodynamics.
"""
abstract type ParameterVector{F, MU, LU, TU, AU, N, T, B} <: ParameterizedLabelledArray{F, 1, T, SLArray{Tuple{N},F,1,N,T}} end

"""
Returns the mass unit of the parameter vector.
"""
Base.@pure massunit(::ParameterVector{F, MU, LU, TU, AU, N}) where {F, MU, LU, TU, AU, N} = MU

"""
Returns the length unit of the parameter vector.
"""
Base.@pure lengthunit(::ParameterVector{F, MU, LU, TU, AU, N}) where {F, MU, LU, TU, AU, N} = LU

"""
Returns the time unit of the parameter vector.
"""
Base.@pure timeunit(::ParameterVector{F, MU, LU, TU, AU, N}) where {F, MU, LU, TU, AU, N} = TU

"""
Returns the angular unit of the parameter vector.
"""
Base.@pure angularunit(::ParameterVector{F, MU, LU, TU, AU, N}) where {F, MU, LU, TU, AU, N} = AU

"""
Returns the velocity unit of the state vector.
"""
velocityunit(::ParameterVector{F, MU, LU, TU, AU}) where {F, MU, LU, TU, AU} = LU/TU

"""
Returns the mass-parameter unit of the state vector.
"""
massparamunit(::ParameterVector{F, MU, LU, TU, AU}) where {F, MU, LU, TU, AU} = LU^3/TU^2

"""
Returns the length of the parameter vector.
"""
Base.@pure Base.length(::ParameterVector{F, MU, LU, TU, AU, N}) where {F, MU, LU, TU, AU, N} = N

"""
Returns the size of the parameter vector.
"""
Base.@pure Base.size(::ParameterVector{F, MU, LU, TU, AU, N}) where {F, MU, LU, TU, AU, N} = (N,)

"""
Returns the name of the parameter vector.
"""
name(::ParameterVector{F, MU, LU, TU, AU, N, T, B}) where {F, MU, LU, TU, AU, N, T, B} = B

"""
All parameters required for the Restricted Two-body Problem.
"""
struct R2BPParameters{F, MU, LU, TU, AU, B} <: ParameterVector{F, MU, LU, TU, AU, 1, (:μ,), B}
    __rawdata::SLArray{Tuple{1}, F, 1, 1, (:μ,)}

    function R2BPParameters(μ; massunit=u"kg", lengthunit=u"km", timeunit=u"s", angularunit=u"°", name=:Primary)
        μμ = (μ isa Unitful.AbstractQuantity) ? ustrip(lengthunit^3/timeunit^2, μ) : μ
        F = eltype(μμ)
        F isa AbstractFloat || (F = Float64)
        B = Symbol(name)
        return new{F, massunit, lengthunit, timeunit, angularunit, B}(SLArray{Tuple{1}, F, 1, 1, (:μ,)}(μμ))
    end
end

"""
Returns the normalized mass parameter, `μ`.
"""
get_μ(sys::R2BPParameters) = sys[1]

"""
Displays `R2BPParameters`.
"""
function Base.show(io::IO, state::R2BPParameters; showfloats=true, space="")
    println(io, space, "R2BP parameters ", name(state) == "Primary" ? "" : "for the $(name(state)) system", showfloats ? " with eltype $(eltype(state))" : "", "\n")
    println(io, space, "  ", "μ = ", get_μ(state), " ", massparamunit(state))
end
Base.show(io::IO, ::MIME"text/plain", state::R2BPParameters; kwargs...) = Base.show(io, state; kwargs...)

"""
Converts types and units for a `R2BPParameters`.
"""
function Base.convert(::Type{R2BPParameters{F, MU, LU, TU, AU}}, system::R2BPParameters) where {F, MU, LU, TU, AU}
    B = name(system)
    μ = ustrip(LU^3/TU^2, get_μ(system) * massparamunit(system)) |> F
    return R2BPParameters(μ; massunit=MU, lengthunit=LU, timeunit=TU, angularunit=AU, name=B)
end

"""
All parameters required for the Circular Restricted Three-body Problem.
"""
struct CR3BPParameters{F, MU, LU, TU, AU, B} <: ParameterVector{F, MU, LU, TU, AU, 1, (:μ,), B}
    __rawdata::SLArray{Tuple{1}, F, 1, 1, (:μ,)}

    function CR3BPParameters(μ::Real; massunit=u"kg", lengthunit=missing, timeunit=missing, angularunit=u"°", primary=:Primary, secondary=:Secondary)
        @assert μ ≤ 1//2 "Nondimensional mass parameter must be less than 1/2, by definition! Received $μ."
        F = typeof(μ)
        F isa AbstractFloat || (F = Float64)
        B = (Symbol(primary), Symbol(secondary))
        return new{F, massunit, lengthunit, timeunit, angularunit, B}(SLArray{Tuple{1}, F, 1, 1, (:μ,)}(μ))
    end

    function CR3BPParameters(μ₁::Unitful.AbstractQuantity, μ₂::Unitful.AbstractQuantity, a::Unitful.Length; massunit=u"kg", angularunit=u"°", primary=:Primary, secondary=:Secondary)
        ∑μᵢ = μ₁ + μ₂
        μ  = min(μ₁, μ₂) / ∑μᵢ
        TU = 2π * upreferred(√(a^3 / ∑μᵢ))
        F = typeof(μ)
        F isa AbstractFloat || (F = Float64)
        B = (Symbol(primary), Symbol(secondary))
        return new{F, massunit, a, TU, angularunit, B}(SLArray{Tuple{1}, F, 1, 1, (:μ)}(F(μ)))
    end
end

"""
Returns the normalized mass parameter, `μ`.
"""
get_μ(sys::CR3BPParameters) = sys[1]

"""
Displays a `CR3BPParameters` instance.
"""
function Base.show(io::IO, sys::CR3BPParameters; showfloats=true, space="")
    sysname   = all(name(sys) .== (:Primary, :Secondary)) ? "" : begin (n1, n2) = name(sys); string(" for the ", n1, "-", n2, " system") end
    println(io, space, "CR3BP parameters", sysname, showfloats ? " with eltype $(eltype(sys))" : "", "\n")
    println(io, space, "  ", "μ = ", get_μ(sys))
end
Base.show(io::IO, ::MIME"text/plain", sys::CR3BPParameters; kwargs...) = Base.show(io::IO, sys::CR3BPParameters; kwargs...)

"""
Converts types and units for a `CR3BPParameters`.
"""
function Base.convert(::Type{CR3BPParameters{F, MU, LU, TU, AU}}, system::CR3BPParameters) where {F, MU, LU, TU, AU}
    B = name(system)
    μ  = get_μ(system)
    μ₁, μ₂ = let
        ∑μᵢ = TU^3 / (TU / 2π)^2
        μ₂ = μ * ∑μᵢ
        μ₁ = ∑μᵢ - μ₂

        μ₁, μ₂
    end

    μₙ = min(μ₁, μ₂) / (μ₁ + μ₂)
    return CR3BPParameters(μₙ; massunit=MU, lengthunit=LU, timeunit=TU, angularunit=AU, primary=name(system)[1], secondary=name(system)[2])
end


text:
using HomotopyContinuation2.ModelKit

@testset "ModelKit" begin
    @testset "SymEngine" begin
        @test Expression(MathConstants.catalan) isa Expression
        @test Expression(MathConstants.e) isa Expression
        @test Expression(MathConstants.pi) isa Expression
        @test Expression(MathConstants.γ) isa Expression
        @test Expression(Float16(2.1)) isa Expression
        @test Expression(Float32(2.1)) isa Expression
        @test Expression(big(2.1)) isa Expression
        @test Expression(Int32(2)) isa Expression
        @test Expression(Int16(2)) isa Expression
        @test Expression(Int128(2)) isa Expression
        @test Expression(2 // 3) isa Expression
        @test convert(BigFloat, Expression(big(2.1))) == big(2.1)
        @test ModelKit.is_number(Expression(2))
        @test ModelKit.to_number(Expression(big(2.1))) isa BigFloat
        @test ModelKit.to_number(Expression(2 // 3)) isa Rational{Int}
        @test ModelKit.to_number(Expression(2) // 3) isa Rational{Int}
        @test ModelKit.to_number(Expression(2) // Expression(3)) isa Rational{Int}
        @test ModelKit.to_number(Expression(big(2)^129)) isa BigInt
        @test complex(Expression(1), Expression(2)) == Expression(1) + im * 2
        @test complex(Expression(1), 2) == Expression(1) + im * 2
        @test complex(1, Expression(2)) == Expression(1) + im * 2
        @test +Expression(2) == 2
        @test Expression(HC2.DoubleDouble.DoubleF64(2)) == big(2.0)
    end

    @testset "Variables" begin
        @var a b x[1:2] y[1:2, 1:3]

        @test a isa Variable
        @test b isa Variable
        @test x isa Vector{Variable}
        @test length(x) == 2
        @test y isa Matrix{Variable}
        @test size(y) == (2, 3)
        @test join(x, " ") == "x₁ x₂"

        @test reshape(sort(vec(y)), 2, 3) == y
        @var x[9:11]
        @test sort(x) == x

        @var c, d
        @test c isa Variable
        @test d isa Variable
        @test variables(c) == [c]
        @test c === c
        @test c !== copy(c)

        let
            c2 = @unique_var c
            @test c != c2
        end

        f = a + b
        @test variables(f) == [a, b]
        @test nvariables(f) == 2
        g = x[1]
        @test variables([f, g]) == [a, b, x[1]]
        @test nvariables([f, g]) == 3
    end

    @testset "Subs" begin
        @var x y w z u
        f = x^2 * (x + y * w)

        @test subs(f, x => z) == z^2 * (z + w * y)
        @test subs([f], x => z) == [z^2 * (z + w * y)]
        @test subs(f, [x, y] => [z^2, z + 2]) == z^4 * (w * (2 + z) + z^2)
        @test subs(f, [x, y] => [z^2, z + 2], w => u) == z^4 * (u * (2 + z) + z^2)
        @test subs(f, x => z^2, y => 3, w => u) == z^4 * (3 * u + z^2)
        @test subs(f, Dict(x => z^2, y => 3, w => u)) == z^4 * (3 * u + z^2)
    end

    @testset "Evaluation" begin
        @var x y w
        f = x^2 * (x + y * w)
        @test evaluate([f, f], [x, y, w] => [2, 3, -5]) == [-52, -52]
        @test [f, 2f]([x, y, w] => [2, 3, -5]) == [-52, -104]
    end

    @testset "Linear Algebra" begin
        @var x[1:2, 1:2]
        @test det(x) == -x[2, 1] * x[1, 2] + x[2, 2] * x[1, 1]
        @test det(x') == -x[2, 1] * x[1, 2] + x[2, 2] * x[1, 1]
        @test det(transpose(x)) == -x[2, 1] * x[1, 2] + x[2, 2] * x[1, 1]
    end

    @testset "Differentiation" begin
        @var x y

        f = x^2 + y^2
        g = x^3 + y^3

        @test differentiate(f, x) == 2x
        @test differentiate(f, [x, y]) == [2x, 2y]

        @test differentiate([f, g], x) == [2x, 3 * x^2]
        @test differentiate([f, g], [x, y]) == [2x 2y; 3 * x^2 3 * y^2]
    end

    @testset "Expand" begin
        @var x y
        @test expand((x + y)^2) == 2 * x * y + x^2 + y^2
    end

    @testset "to_dict" begin
        @var x y a
        @test ModelKit.to_dict(x^2 + y * x * a + y * x, [x, y]) ==
              Dict([2, 0] => Expression(1), [1, 1] => a + 1)
        @test ModelKit.degree(x * y^5 + y^2, [x, y]) == 6
    end

    @testset "Modeling" begin
        @testset "Bottleneck" begin
            @var x y z
            f = [
                (0.3 * x^2 + 0.5z + 0.3x + 1.2 * y^2 - 1.1)^2 +
                (0.7 * (y - 0.5x)^2 + y + 1.2 * z^2 - 1)^2 - 0.3,
            ]

            I = let x = [x, y, z]
                n, m = length(x), length(f)
                @unique_var y[1:n] v[1:m] w[1:m]
                J = differentiate(f, x)
                f′ = f(x => y)
                J′ = J(x => y)
                Nx = (x - y) - J' * v
                Ny = (x - y) - J′' * w
                System([f; f′; Nx; Ny], [x; y; v; w])
            end
            @test I isa System
            @test size(I) == (8, 8)
        end

        @testset "Steiner" begin
            @var x[1:2] a[1:5] c[1:6] y[1:2, 1:5]

            #tangential conics
            f = sum([a; 1] .* monomials(x, 2))
            ∇ = differentiate(f, x)
            #5 conics
            g = sum(c .* monomials(x, 2))
            ∇_2 = differentiate(g, x)
            #the general system
            #f_a_0 is tangent to g_b₀ at x₀
            function Incidence(f, a₀, g, b₀, x₀)
                fᵢ = f(x => x₀, a => a₀)
                ∇ᵢ = ∇(x => x₀, a => a₀)
                Cᵢ = g(x => x₀, c => b₀)
                ∇_Cᵢ = ∇_2(x => x₀, c => b₀)

                [fᵢ; Cᵢ; det([∇ᵢ ∇_Cᵢ])]
            end
            @var v[1:6, 1:5]
            I = vcat(map(i -> Incidence(f, a, g, v[:, i], y[:, i]), 1:5)...)
            F = System(I, [a; vec(y)], vec(v))
            @test size(F) == (15, 15)
        end

        @testset "Reach plane curve" begin
            @var x y
            f = (x^3 - x * y^2 + y + 1)^2 * (x^2 + y^2 - 1) + y^2 - 5
            ∇ = differentiate(f, [x; y]) # the gradient
            H = differentiate(∇, [x; y]) # the Hessian

            g = ∇ ⋅ ∇
            v = [-∇[2]; ∇[1]]
            h = v' * H * v
            dg = differentiate(g, [x; y])
            dh = differentiate(h, [x; y])

            ∇σ = g .* dh - ((3 / 2) * h) .* dg

            F = System([v ⋅ ∇σ; f], [x, y])
            @test size(F) == (2, 2)
        end
    end

    @testset "Horner" begin
        # use fano scheme as test
        n = 4
        s = 1
        @var x[1:n]
        ν = monomials(x, 5)
        N = length(ν)
        q₀ = randn(ComplexF64, N - 1)
        @var q[1:N-1]
        F = sum(q[i] * ν[i] for i = 1:N-1) + 1

        @var a[1:n-1] b[1:n-1] t
        L = [a; 1] .* t + [b; 0]

        G = subs(F, x => L)

        FcapL = collect(values(ModelKit.to_dict(expand(G), [t])))

        @test all(FcapL) do f
            g = horner(f)
            expand(g) == f
        end
        @test all(FcapL) do f
            g = horner(f, [a; b])
            expand(g) == f
        end
    end

    @testset "Rand / dense poly" begin
        @var x y
        f, c = dense_poly([x, y], 3)
        @test length(c) == 10
        g = rand_poly(Float64, [x, y], 3)
        @test subs(f, c => coefficients(g, [x, y])) == g
        _, coeffs = exponents_coefficients(g, [x, y])
        @test subs(f, c => coeffs) == g
    end

    @testset "System" begin
        @var x y a b
        f = [(x + y)^3 + x^2 + x + 5y + 3a, 2 * x^2 + b]
        F = System(f, [x, y], [b, a])
        @test nvariables(F) == 2
        @test nparameters(F) == 2
        @test variables(F) == [x, y]
        @test parameters(F) == [b, a]
        show_F = """
        System of length 2
         2 variables: x, y
         2 parameters: b, a

         3*a + x + 5*y + x^2 + (x + y)^3
         b + 2*x^2"""
        @test sprint(show, F) == show_F
        @test degrees(F) == [3, 2]

        T = ModelKit.compile(F; optimizations = false)
        F2 = ModelKit.interpret(T)
        @test F == F2
        @test sprint(show, T) == "Compiled: $show_F"
        @test size(T) == size(F) == (2, 2)

        T2 = ModelKit.compile(F; optimizations = true)
        F3 = ModelKit.interpret(T2)
        @test expand.(F.expressions) == expand.(F3.expressions)
    end

    @testset "System variables groups + homogeneous" begin
        @var x y z
        f = System([
            (x^2 + y^2 + x * y - 3 * z^2) * (x + 3z),
            (x^2 + y^2 + x * y - 3 * z^2) * (y - x + 2z),
            2x + 5y - 3z,
        ])
        @test is_homogeneous(f)
        multi_degrees(f) == [1 1; 2 0]

        @var x y z v w
        g = System([x * y - 2v * w, x^2 - 4 * v^2], variable_groups = [[x, v], [y, w]])
        @test is_homogeneous(g)
        multi_degrees(g) == [1 1; 2 0]
    end

    @testset "Homotopy" begin
        @var x y z t

        h = [x^2 + y + z + 2t, 4 * x^2 * z^2 * y + 4z - 6x * y * z^2]
        H = Homotopy(h, [x, y, z], t)

        @test nvariables(H) == 3
        @test nparameters(H) == 0
        @test variables(H) == [x, y, z]
        @test parameters(H) == []
        show_H = """
        Homotopy in t of length 2
         3 variables: x, y, z

         2*t + y + z + x^2
         4*z - 6*x*y*z^2 + 4*x^2*y*z^2"""

        @test sprint(show, H) == show_H

        T = ModelKit.compile(H)
        H2 = ModelKit.interpret(T)
        @test H == H2
        @test sprint(show, T) == "Compiled: $show_H"
        @test size(T) == size(H) == (2, 3)
    end

    @testset "Codegen helpers" begin
        @test ModelKit.sqr(3 + 2im) == (3 + 2im)^2
        @test ModelKit.sqr(3) == 3^2
    end

    @testset "Instructions" begin

        @testset "Simple dt" begin
            @var x a b t
            a = 0.5
            b = 0.1
            h = (x - (t - (a + im * b))^2) * (x + (t - (a + im * b))^2)
            list, _ = ModelKit.instruction_list([h])
            D = ModelKit.DiffMap()
            D[:t, 1] = 1
            slp = @eval ModelKit begin
                function dt_test1(x, t)
                    $(ModelKit.to_expr(ModelKit.univariate_diff!(list, 1, D)))
                end
            end
            expr1 = ModelKit.expand(ModelKit.dt_test1(x, t))
            expr2 = ModelKit.expand(ModelKit.differentiate(h, t))
            @test expr1 == expr2
        end
        @testset "Higher order pow diff" begin
            for d = 2:10
                @var x
                f = x^d
                list, _ = ModelKit.instruction_list([f])

                D = ModelKit.DiffMap()
                D[:x, 1] = :x1
                D[:x, 2] = :x2
                D[:x, 3] = :x3

                @eval ModelKit begin
                    function __diff_4_pow(x, x1, x2, x3)
                        $(ModelKit.to_expr(ModelKit.univariate_diff!(list, 4, D)))
                    end
                end

                @var x1 x2 x3 λ

                @test expand(
                    subs(
                        differentiate(
                            subs(f, x => x .+ λ .* x1 .+ λ^2 .* x2 + λ^3 .* x3),
                            λ,
                            4,
                        ),
                        λ => 0,
                    ) / 24,
                ) == expand(ModelKit.__diff_4_pow(x, x1, x2, x3))
            end
        end

        @testset "Higher order mul" begin
            @var x y
            f = x * y
            list, _ = ModelKit.instruction_list([f])
            D = ModelKit.DiffMap()
            D[:x, 1] = :x1
            D[:x, 2] = :x2
            D[:x, 3] = :x3
            D[:y, 1] = :y1
            D[:y, 2] = :y2
            D[:y, 3] = :y3

            @eval ModelKit begin
                function __diff_4_mul__(x, y, (x1, x2, x3), (y1, y2, y3))
                    $(ModelKit.to_expr(ModelKit.univariate_diff!(list, 4, D)))
                end
            end
            @var x1 x2 x3 y1 y2 y3 λ

            @test expand(
                subs(
                    differentiate(
                        subs(
                            f,
                            x => x .+ λ .* x1 .+ λ^2 .* x2 + λ^3 .* x3,
                            y => y .+ λ .* y1 .+ λ^2 .* y2 + λ^3 .* y3,
                        ),
                        λ,
                        4,
                    ),
                    λ => 0,
                ) / 24,
            ) == expand(ModelKit.__diff_4_mul__(x, y, (x1, x2, x3), (y1, y2, y3)))
        end

        @testset "Higher order plus" begin
            @var x y
            f = x + y
            list, _ = ModelKit.instruction_list([f])

            D = ModelKit.DiffMap()
            D[:x, 1] = :x1
            D[:x, 2] = :x2
            D[:x, 3] = :x3
            D[:y, 1] = :y1
            D[:y, 2] = :y2
            D[:y, 3] = :y3

            @eval ModelKit begin
                function __diff_3_plus__(x, y, (x1, x2, x3), (y1, y2, y3))
                    $(ModelKit.to_expr(ModelKit.univariate_diff!(list, 3, D)))
                end
            end

            @var x1 x2 x3 y1 y2 y3 λ

            @test expand(
                subs(
                    differentiate(
                        subs(
                            f,
                            x => x .+ λ .* x1 .+ λ^2 .* x2 + λ^3 .* x3,
                            y => y .+ λ .* y1 .+ λ^2 .* y2 + λ^3 .* y3,
                        ),
                        λ,
                        3,
                    ),
                    λ => 0,
                ) / 6,
            ) == expand(ModelKit.__diff_3_plus__(x, y, (x1, x2, x3), (y1, y2, y3)))
        end
    end

    @testset "Homotopy codegen (Katsura(3))" begin
        n = 3
        @var x[0:n] ẋ[0:n] ẍ[0:n] x3[0:n] t γ
        K = [
            (
                sum(x[abs(l)+1] * x[abs(m - l)+1] for l = -n:n if abs(m - l) <= n) - x[m+1] for m = 0:n-1
            )...,
            x[1] + 2 * sum(x[i+1] for i = 1:n) - 1,
        ]

        h = γ .* t .* [x[1:n] .^ 2 .- 1; x[n+1] - 1] + (1 - t) .* K
        H = ModelKit.Homotopy(h, x, t, [γ])
        TH = ModelKit.compile(H)

        tx3 = TaylorVector{4}(zeros(Expression, 4, 4))
        y, y1, y2, y3 = vectors(tx3)
        y .= x
        y1 .= ẋ
        y2 .= ẍ
        y3 .= x3

        dt1 =
            ModelKit.taylor!(zeros(Expression, 4), Val(1), TH, TaylorVector{1}(tx3), t, [γ])
        dt2 =
            ModelKit.taylor!(zeros(Expression, 4), Val(2), TH, TaylorVector{2}(tx3), t, [γ])
        dt3 =
            ModelKit.taylor!(zeros(Expression, 4), Val(3), TH, TaylorVector{3}(tx3), t, [γ])
        dt4 =
            ModelKit.taylor!(zeros(Expression, 4), Val(4), TH, TaylorVector{4}(tx3), t, [γ])

        @var λ
        Hd1 = subs(H.expressions, t => t + λ)
        true_dt1 = subs(differentiate(Hd1, λ, 1), λ => 0)

        Hd2 = subs(H.expressions, x => x .+ λ .* ẋ, t => t + λ)
        true_dt2 = subs(differentiate(Hd2, λ, 2), λ => 0) / 2

        Hd3 = subs(H.expressions, x => x .+ λ .* ẋ .+ λ^2 .* ẍ, t => t + λ)
        true_dt3 = subs(differentiate(Hd3, λ, 3), λ => 0) / 6

        Hd4 = subs(H.expressions, x => x .+ λ .* ẋ .+ λ^2 .* ẍ .+ λ^3 .* x3, t => t + λ)
        true_dt4 = subs(differentiate(Hd4, λ, 4), λ => 0) / 24

        @test expand.(ModelKit.evaluate(TH, x, t, [γ])) == expand.(h)
        @test expand.(ModelKit.jacobian(TH, x, t, [γ])) == expand.(differentiate(h, x))
        @test expand.(dt1) == expand.(true_dt1)
        @test expand.(dt2) == expand.(true_dt2)
        @test expand.(dt3) == expand.(true_dt3)
        @test expand.(dt4) == expand.(true_dt4)
    end

    @testset "taylor! - system" begin
        @var x[1:2] ẋ[1:2] ẍ[1:2] x3[1:2] p[1:2] ṗ[1:2]
        f = [(x[1] + x[2])^3 + x[1]^2 + x[1] + 5 * x[2] + 3 * p[1], 2 * x[1]^2 + p[2]]
        F = System(f, x, p)
        TF = ModelKit.compile(F)

        tx = TaylorVector{4}([x'; ẋ'; ẍ'; x3'])
        tv = TaylorVector{5}(Expression, 2)
        v, v1, v2, v3, v4 = vectors(tv)
        @var λ

        ModelKit.taylor!(TaylorVector{3}(tv), TF, TaylorVector{2}(tx), p)
        Fd2 = subs(F.expressions, x => x .+ λ .* ẋ)
        true_v1 = subs(differentiate(Fd2, λ, 1), λ => 0)
        true_v2 = subs(differentiate(Fd2, λ, 2), λ => 0) / 2
        @test expand.(v) == expand.(f)
        @test expand.(v1) == expand.(true_v1)
        @test expand.(v2) == expand.(true_v2)

        ModelKit.taylor!(TaylorVector{4}(tv), TF, TaylorVector{3}(tx), p)
        Fd3 = subs(F.expressions, x => x .+ λ .* ẋ .+ λ^2 .* ẍ)
        true_v1 = subs(differentiate(Fd3, λ, 1), λ => 0)
        true_v2 = subs(differentiate(Fd3, λ, 2), λ => 0) / 2
        true_v3 = subs(differentiate(Fd3, λ, 3), λ => 0) / 6
        @test expand.(v) == expand.(f)
        @test expand.(v1) == expand.(true_v1)
        sub = variables(v2) => rand(1:10_000, 6)
        @test v2(sub) == true_v2(sub)
        @test v3(sub) == true_v3(sub)
    end
end


text:
using ElectronDisplay
using Electron
using VegaLite
using DataFrames
using Test

@testset "ElectronDisplay" begin

p1 =  DataFrame(x=rand(10),y=rand(10)) |> @vlplot(:point, x=:x, y=:y)

f = display(p1)

@test f isa Electron.Window

end


text:
# Jonas August, 2/2022

#=
This problem was asked by Amazon.

Given a N by M matrix of numbers, print out the matrix in a clockwise spiral.

For example, given the following matrix:

[[1,  2,  3,  4,  5],
 [6,  7,  8,  9,  10],
 [11, 12, 13, 14, 15],
 [16, 17, 18, 19, 20]]

You should print out the following:

1
2
3
4
5
10
15
20
19
18
17
16
11
6
7
8
9
14
13
12

=#

function line(m)
    for x in m
        println(x)
    end
end

function spiral(A::AbstractMatrix)
    m, n = size(A)
    if m == 0 || n == 0 return end
    if m > 0 line(A[begin,begin:end-1]) end
    if n > 1 line(A[begin:end-1,end]) end
    if m > 1 line(A[end,end:-1:begin+1]) end
    if n > 0 line(A[end:-1:begin+1,begin]) end
    spiral(A[begin+1:end-1,begin+1:end-1])
end

A = reshape(1:20, 5, 4)'
spiral(A)

max_stars_repo_path: Spiral printout of matrix - Problem 65.jl


text:
# ==============================================================================
# Max-Plus eigenvalues eigenvectors (Howard algorithm).
#
# The Max-Plus version of Howard's algorithm is described in the paper:
#    "Numerical computation of spectral elements in max-plus algebra"
# by Jean Cochet-Terrasson, Guy Cohen, Stephane Gaubert, Michael Mc Gettrick,
# Jean-Pierre Quadrat
#
# 1997: Stephane.Gaubert@inria.fr: Original C version for ScicosLab
#       http://www.cmap.polytechnique.fr/~gaubert/HOWARD2.html
# 2021: quentin.quadrat@gmail.com: Julia portage of the C version.

"""
    sparse_ij(S::SpaMP)

Concat indices of a sparse matrix line by line.

# Examples
```julia-repl
julia> A = [1 2; 3 4; 5 6]
3×2 Matrix{Int64}:
 1  2
 3  4
 5  6

julia> vcat(A'...)
6-element Vector{Int64}:
 1
 2
 3
 4
 5
 6

julia> sparse_ij(sparse(MP([1.0 2; 3 4])))
```
"""
function sparse_ij(S::SparseMatrixCSC)
    (I,J,) = findnz(S)
    # FIXME: can we Vector(reshape([J I]', (1, 2 * size(I,1)))) ?
    return Matrix(reshape([J I]', (1, 2 * size(I,1))))
end

# -----------------------------------------------------------------------------
# Howard context
# -----------------------------------------------------------------------------
struct Context
    # INPUTS: Sparse matrix given as input in howard() function
    ij::Matrix{Int64}  # Index of the sparse matrix to egein FIXME shall be Vector
    A::Vector{Float64} # Values in classic algebra of the sparse matrix
    nnodes::Int      # sparse matrix dimension
    narcs::Int       # sparse matrix dimension

    # INTERNALS
    newpi::Vector{Int64}
    piinv_idx::Vector{Int64}
    piinv_succ::Vector{Int64}
    piinv_elem::Vector{Int64}
    piinv_last::Vector{Int64}
    visited::Vector{Bool}
    component::Vector{Int64}
    c::Vector{Float64}
    newc::Vector{Float64}
    vaux::Vector{Float64}
    newchi::Vector{Float64}

    # OUTPUTS
    v::Vector{Float64}     # Bias vectors
    chi::Vector{Float64}   # Cycle time vector
    pi::Vector{Int64}      # Optimal policy

    function Context(S::SparseMatrixCSC{T,U}) where {T,U}
        A = convert(SparseMatrixCSC{Float64,U},transpose(S)).nzval
        nnodes = size(S,1)
        narcs = size(A,1)
        new(## INPUTS: Scilab: [ij,a,s]=spget(sparse(M))
            sparse_ij(sparse(S')), # FIXME sparse(S') can be simplified but needed for findnz() ?
            A,
            nnodes,
            narcs,
            # INTERNALS
            zeros(Int, nnodes),
            zeros(Int, nnodes),
            zeros(Int, nnodes),
            zeros(Int, nnodes),
            zeros(Int, nnodes),
            zeros(Bool, nnodes),
            zeros(Int, nnodes),
            zeros(Float64, nnodes),
            zeros(Float64, nnodes),
            fill(-Inf, (nnodes,)), # epsilon value
            zeros(Float64, nnodes),
            # OUTPUTS
            zeros(Float64, nnodes),
            zeros(Float64, nnodes),
            fill(673720360, nnodes))
    end
end

# -----------------------------------------------------------------------------
function show_info(con::Context, iteration::Int)
    (@printf "*** ITERATION %d of Max Plus Howard Algorithm *** \n" iteration)
    (@printf "PI=")
    show(stdout, con.pi)
    (@printf "\nc=")
    show(stdout, con.c)
    (@printf "\nchi=")
    show(stdout, con.chi)
    (@printf "\nv=")
    show(stdout, con.v)
    (@printf "\n")
end

# -----------------------------------------------------------------------------
function show_info_improve_bias(con::Context, i::Int)
    #(@printf "type 2 improvement\n")
    I = con.ij[2i-1]
    J = con.ij[2i]
    #(@printf "Improvement of the BIAS at node %d\n" I)
    #(@printf "A[%d] + v[%d] - chi[%d] - v[%d]= %f + %f - %f - %f = %f > 0\n" i J I I con.A[i] con.v[J] con.chi[I] con.v[I] con.A[i]+con.v[J]-con.chi[I]-con.v[I])
end

# -----------------------------------------------------------------------------
function sanity_checks(con::Context)
    # Check matrix dimension
    if ((con.nnodes < 1) || (con.narcs < 1))
        throw(DimensionMismatch("Sparse matrix shall not be empty"))
    end

    # Check if nodes have successor
    u = zeros(Bool, con.nnodes)
    #(@printf "sanity_checks %d:\n" con.narcs)
    for i = 1:con.narcs
        #(@printf "  %d\n" i)
        #(@printf " => ij[%d]: %d\n"  2i-1  con.ij[2i-1])
        u[con.ij[2i-1]] = true
    end

    for i = 1:con.nnodes
        if (u[i] == false)
            throw(InitError("Node X has no successor"))
        end
    end
    #(@printf "U=")
    #show(stdout, u)
    #(@printf "\n")
end

# -----------------------------------------------------------------------------
# This epsilon value is somehow arbitrary.
function epsilon(con::Context)
    m, M = extrema(con.A)
    (M - m + 1.0) * 0.000000001
end

# -----------------------------------------------------------------------------
function initial_policy(con::Context)
    for i = 1:con.narcs
        #(@printf "\ncon.vaux[%d]=%f A[%d]=%f\n" con.ij[2i-1] con.vaux[con.ij[2i-1]] i con.A[i])
        if (con.vaux[con.ij[2i-1]] <= con.A[i])
            #(@printf "=> pi[ij[2i-1]] = ij[2i] => pi[%d] = %d\n" con.ij[2i-1] con.ij[2i])
            con.pi[con.ij[2i-1]] = con.ij[2i]
            con.c[con.ij[2i-1]] = con.A[i]
            con.vaux[con.ij[2i-1]] = con.A[i]
            #(@printf "=> c[ij[2i-1]] = A[i] => c[%d] = %f\n" con.ij[2i-1] con.A[i])
            #(@printf "=> vaux[ij[2i-1]] = A[i] => vaux[%d] = %f\n" con.ij[2i-1] con.A[i])
        end
    end

    #(@printf "PI=")
    #show(stdout, con.pi)
    #(@printf "\nC=")
    #show(stdout, con.c)
    #(@printf "\nvaux=")
    #show(stdout, con.vaux)
end

# -----------------------------------------------------------------------------
function update_policy(con::Context)
    con.pi .= con.newpi
    con.c .= con.newc
    con.vaux .= con.v
end

# -----------------------------------------------------------------------------
function build_inverse(con::Context)

    con.piinv_idx .= -1
    con.piinv_last .= -1

    ptr = 1
    for i = 1:con.nnodes
        j = con.pi[i]
        if (con.piinv_idx[j] == -1)
            con.piinv_succ[ptr] = -1
            con.piinv_elem[ptr] = i
            con.piinv_last[j] = ptr
            con.piinv_idx[j] = ptr
        else
            con.piinv_succ[ptr] = -1
            con.piinv_elem[ptr] = i
            locus = con.piinv_last[j]
            con.piinv_succ[locus] = ptr
            con.piinv_last[j] = ptr
        end
        ptr += 1
    end

    #for i = 1:con.nnodes
    #    (@printf "%d:\n" i)
    #    (@printf "  piinv_succ=%d\n" con.piinv_succ[i])
    #    (@printf "  piinv_elem=%d\n" con.piinv_elem[i])
    #    (@printf "  piinv_last=%d\n" con.piinv_last[i])
    #    (@printf "  piinv_idx=%d\n" con.piinv_idx[i])
    #end
end

# -----------------------------------------------------------------------------
function depth_first_label(con::Context, lambda::Float64, color::Int, i::Int)
    a = con.piinv_idx[i]
    while ((a != -1) && (con.visited[con.piinv_elem[a]] == false))
        nexti = con.piinv_elem[a]
        con.visited[nexti] = true
        con.v[nexti] = -lambda + con.c[nexti] + con.v[i]
        con.component[nexti] = color
        con.chi[nexti]= lambda
        depth_first_label(con, lambda, color, nexti)
        a = con.piinv_succ[a]
    end
end

# -----------------------------------------------------------------------------
function visit_from(con::Context, initialpoint::Int, color::Int)
    #(@printf "Visiting from node %d color=%d\n" initialpoint color)
    index = initialpoint
    con.component[index] = color
    newindex = con.pi[index]
    while true
        con.component[newindex] = color
        index = newindex
        newindex = con.pi[index]
        (con.component[newindex] == 0) || break;
    end

    # A cycle has been detected, since newindex is already visited
    weight::Float64 = 0
    length::Float64 = 0
    i = index

    while true
        weight += con.c[i]
        length += 1
        i = con.pi[i]
        (i != index) || break
    end
    lambda = weight / length
    con.v[i] = con.vaux[i]
    con.chi[i] = lambda
    depth_first_label(con, lambda, color, index)
end

# -----------------------------------------------------------------------------
function value(con::Context, components::Int)
    initialpoint = 1
    color = 1

    # init_depth_first()
    con.visited .= false
    con.component .= 0

    while true
        visit_from(con, initialpoint, color);
        while ((initialpoint <= con.nnodes) && (con.component[initialpoint] != 0))
            initialpoint += 1
        end
        color += 1
        (initialpoint <= con.nnodes) || break
    end

    components = color - 1
end

# -----------------------------------------------------------------------------
function first_order_improvement(con::Context, improved::Bool)
    for i = 1:con.narcs
        if (con.chi[con.ij[2i]] > con.newchi[con.ij[2i-1]])
            show_info_improve_chi(i)
            improved = true;
            con.newpi[con.ij[2i-1]] = con.ij[2i]
            con.newchi[con.ij[2i-1]] = con.chi[con.ij[2i]]
            con.newc[con.ij[2i-1]] = con.A[i];
        end
    end
end

# -----------------------------------------------------------------------------
function second_order_improvement(con::Context, components::Int, improved::Bool, ϵ::Float64)
    if (components > 1)
        for i = 1:con.narcs
            if (con.chi[con.ij[2i]] == con.newchi[con.ij[2i-1]])
                w = a[i] + v[con.ij[2i]] - con.chi[con.ij[2i-1]]
                if (w > con.vaux[con.ij[2i-1]] + ϵ)
                    show_info_improve_bias(con, i)
                    improved = true
                    con.vaux[ij[2i-1]] = w
                    con.newpi[ij[2 - 1]] = con.ij[2i]
                    con.newc[ij[2i-1]] = con.A[i]
                end
            end
        end
    else
        for i = 1:con.narcs
            w = con.A[i]+ con.v[con.ij[2i]] - con.chi[con.ij[2i-1]]
            if (w > con.vaux[con.ij[2i-1]] + ϵ)
                show_info_improve_bias(con, i)
                improved = true
                con.vaux[con.ij[2i-1]] = w
                con.newpi[con.ij[2i-1]] = con.ij[2i]
                con.newc[con.ij[2i-1]] = con.A[i]
            end
        end
    end
end

# -----------------------------------------------------------------------------
function improve(con::Context, components::Int, improved::Bool, ϵ::Float64)
    improved = false

    con.newchi .= con.chi
    con.vaux .= con.v
    con.newpi .= con.pi
    con.newc .= con.c

    if (components > 1)
        first_order_improvement(con, improved, ϵ)
    end

    if (improved == false)
        second_order_improvement(con, components, improved, ϵ)
    end
end

# -----------------------------------------------------------------------------
"""
    howard(S::SpaMP)

TODO

# Example 1: Irreducible matrix, only 1 eigenvalue.
```julia-repl
julia> using SparseArrays

julia> S = sparse(MP([1 2; 3 4]))
2×2 Max-Plus Sparse Matrix with 4 stored entries:
  1   2
  3   4

julia> λ,v = howard(S)
(MP[4, 4], MP[2, 4])

# λ is constant
julia> (S * v) == (λ[1] * v)
true
```

# Example 2: Two blocks diagonal matrix.
```julia-repl
julia> using SparseArrays, LinearAlgebra

julia> S = sparse([mp0 2 mp0; mp1 mp0 mp0; mp0 mp0 2])
3×3 Max-Plus Sparse Matrix with 3 stored entries:
  .   2   .
  0   .   .
  .   .   2

julia> λ,v = howard(S)
(MP[1, 1, 2], MP[1, 0, 2])

# The entries of λ take two values
julia> (S / Matrix(Diagonal(λ))) * v == v
true
```

# Example 3: Block triangular matrix with 2 eigenvalues.
```julia-repl
julia> S = sparse([1 1; mp0 2])
2×2 Max-Plus sparse matrix with 3 stored entries:
  1   1
  .   2

julia> λ,v = howard(S)
(MP[2, 2], MP[1, 2])

julia> (S * v) == (λ[1] * v)
true

# But MP(1) is also eigen value
S * [0; %0] == MP(1) * [0; mp0]
```

# Example 4: Block triangular matrix with 1 eigenvalue
```julia-repl
julia> using SparseArrays, LinearAlgebra

julia> S = sparse([2 1; mp0 mp1])
2×2 Max-Plus sparse matrix with 3 stored entries:
  2   1
  .   0

julia> λ,v = howard(S)
(MP[2, 0], MP[2, 0])

# λ is not constant λ[1] is eigen value
# with eigen vector [v(1);0]
julia> (S / Matrix(Diagonal(λ))) * v == v
true
```
"""
function howard(S::SpaMP)
    improved::Bool = false
    iteration::Int = 0
    components::Int = 0
    MAX_ITERATIONS::Int = 1000

    con = Context(plustimes(S))
    #(@printf "IJ=")
    #show(stdout, con.ij)
    #(@printf "\nA=")
    #show(stdout, con.A)
    #(@printf "\n")
    sanity_checks(con)
    initial_policy(con)
    ϵ = epsilon(con)
    build_inverse(con)

    while true
        value(con, components)
        #show_info(con, iteration)
        improve(con, components, improved, ϵ)
        update_policy(con)
        build_inverse(con)
        iteration += 1
        ((improved != false) && (iteration < MAX_ITERATIONS)) || break
    end

    if (iteration >= MAX_ITERATIONS)
        throw(AssertiobError("Max iteration reached"))
    end

    MP(con.chi), MP(con.v)
end


text:
@testset "basic functionality" begin
    @testset "null outputs" begin
        # `Nothing` output is treated specially: MATLAB `varargout` output is empty, requesting output will error
        @test is_eqq(mx_wrap_jlcall(0, "() -> nothing"), nothing)
        @test_throws MEngineError mx_wrap_jlcall(1, "() -> nothing")

        # `Missing` output corresponds to empty Matrix{Float64}, i.e. with size 0x0
        @test is_eq(mx_wrap_jlcall(1, "() -> missing"), zeros(Float64, 0, 0))
    end

    @testset "nested args/kwargs" begin
        ps = deeply_nested_pairs(; roundtrip = true)
        jl_vec = Any[jl for (jl, mx) in ps]
        mx_vec = Any[mx for (jl, mx) in ps]
        jl_dict = Dict{Symbol, Any}(Symbol("k", i) => v for (i, v) in enumerate(jl_vec))
        mx_dict = Dict{String, Any}("k$i" => v for (i, v) in enumerate(mx_vec))

        @test is_eq(mx_wrap_jlcall(length(jl_vec), "(args...; kwargs...) -> args", (jl_vec...,)), (mx_vec...,))
        @test is_eq(mx_wrap_jlcall(1, "(args...; kwargs...) -> kwargs", (), (; jl_dict...,)), mx_dict)
    end
end

@testset "local project" begin
    @test is_eq(mx_wrap_jlcall(1, "TestProject.dot", ([1.0, 2.0, 3.0],); project = jlcall_test_project(), modules = ["TestProject"], restart = true), 14.0)
    @test is_eq(mx_wrap_jlcall(1, "() -> dirname(Base.active_project())"; restart = false), jlcall_test_project())
end

@testset "setting threads" begin
    @test is_eq(mx_wrap_jlcall(1, "() -> Threads.nthreads()"; threads = 3, restart = true), 3) # Restart Julia with --threads=3
    @test is_eq(mx_wrap_jlcall(1, "() -> Threads.nthreads()"; threads = 4), 3) # Setting threads shouldn't change active session
end

@testset "shared server environments" begin
    # Initialize shared server environment
    @test is_eqq(mx_wrap_jlcall(0; setup = joinpath(@__DIR__, "setup.jl"), shared = true, restart = true), nothing)

    # Call custom library code in persistent stateful environment
    @test is_eq(mx_wrap_jlcall(1, "Setup.mul2", ([1, 2],); shared = true), [2, 4])
    @test is_eq(mx_wrap_jlcall(1, "Setup.mul2", ([3.0, 4.0],); shared = true), [6.0, 8.0])
    @test is_eq(mx_wrap_jlcall(1, "LinearAlgebra.norm", ([3.0, 4.0],); modules = ["LinearAlgebra", "Statistics"], shared = true), 5.0)
    @test is_eq(mx_wrap_jlcall(1, "LinearAlgebra.det", ([1.0 2.0; 3.0 4.0],); shared = true), -2.0)
    @test is_eq(mx_wrap_jlcall(1, "x -> Statistics.mean(Setup.mul2(x))", ([1.0, 2.0, 3.0],); shared = true), 4.0)
end

@testset "unique server environments" begin
    # Initialize unique server environments
    @test is_eqq(mx_wrap_jlcall(0; shared = false, restart = true), nothing)

    # Run custom code in each environment, requiring re-initialization each time
    @test is_eq(mx_wrap_jlcall(1, "Setup.mul2", ([1, 2],); setup = joinpath(@__DIR__, "setup.jl"), shared = false), [2, 4])
    @test_throws MEngineError mx_wrap_jlcall(1, "Setup.mul2", ([1, 2],); shared = false)

    @test is_eq(mx_wrap_jlcall(1, "LinearAlgebra.norm", ([3.0, 4.0],); modules = ["LinearAlgebra", "Statistics"], shared = false), 5.0)
    @test_throws MEngineError mx_wrap_jlcall(1, "LinearAlgebra.det", ([1.0 2.0; 3.0 4.0],); shared = false)
    @test_throws MEngineError mx_wrap_jlcall(1, "x -> Statistics.mean(Setup.mul2(x))", ([1.0, 2.0, 3.0],); shared = false)
end

@testset "server-free" begin
    # Run unique local Julia process for each `jlcall`
    @test is_eq(mx_wrap_jlcall(1, "LinearAlgebra.norm", ([3.0, 4.0],); modules = ["LinearAlgebra"], server = false), 5.0)
    @test_throws MEngineError mx_wrap_jlcall(1, "LinearAlgebra.norm", ([3.0, 4.0],); server = false)
end

@testset "port number" begin
    for port in [2345, 3456]
        @test is_eqq(mx_wrap_jlcall(0, "() -> nothing"; port = port, restart = true), nothing)
    end
end

@testset "extending matlabify" begin
    setup_script = tempname() * ".jl"
    open(setup_script; write = true) do io
        println(io, "MATDaemon.matlabify(v::Base.VersionNumber) = string(v)")
    end
    # Note: restart = true is important, as it tests the ability to call dynamically defined `matlabify` methods
    @test is_eq(mx_wrap_jlcall(1, "() -> Base.VERSION"; restart = true, setup = setup_script), string(Base.VERSION))
end


text:
abstract type QuerySynthesisStrategy <: QueryStrategy end

abstract type DataBasedQss <: QuerySynthesisStrategy end
abstract type ModelBasedQss <: QuerySynthesisStrategy end
abstract type HybridQss <: QuerySynthesisStrategy end

function get_query_objects(qs::QuerySynthesisStrategy,
                           query_data::Array{T, 2},
                           pools::Vector{Symbol},
                           history::Vector{Array{T, 2}})::Array{T, 2} where T <: Real
    return query_synthesis_optimize(qs_score_function(qs, query_data, MLLabelUtils.labelmap(pools)),
                                    qs.optimizer,
                                    query_data,
                                    pools)
end


text:
@test get_scope("""
  function foo()
    |
  end
  """) == "function"

@test get_scope("""
  let x = 1, y = 2
    # comment
    Pkg.add(x, |)
  end
  """) == "Pkg.add"

@test get_scope("""
  let x = 1, y = 2
    # |comment
    Pkg.add(x, )
  end
  """).kind == :comment

@test get_scope("""
  for i in 1:10
    foo(|, y)
  end
  """) == "foo"

@test get_scope("""
  try
    fo|o(, y)
  end
  """) == "try"

# @test get_scope("""
#   function foo()
#     le|t x = 1, y = 2
#       foo()
#     end
#   end
#   """) == "function"

@test get_scope("""
  try
    foo()
  catch e
    |bar()
  end
  """) == "catch"

@test get_scope("""
  if a
    foo()
  else
    |bar()
  end
  """) == "else"

@test Jewel.codemodule(to_cursor("""
  module Foo
  module Bar|
  end
  end
  """)...) == "Foo.Bar"

# Treating end correctly in arrays
@test get_scope("""
  [end|]
  """) == Jewel.Scope(:array, "[")

# Failing gracefully with too many ends
@test get_scope("""
  end|
  """).kind == :toplevel

# Ignore the :end keyword
@test get_scope("""
  if foo
    :end|
  """) == "if"


text:
# a Widget is expected to have a field "widget::PyObject" which is a QWidget
abstract Widget

roundint(x::AbstractFloat) = round(Int, x)

# methods for the python widget
showwidget(w::PyObject) = (w[:showNormal](); w[:raise_](); w[:activateWindow](); nothing)
hidewidget(w::PyObject) = (w[:hide](); nothing)
widgetpos(w::PyObject) = (point = w[:pos](); (point[:x](), point[:y]))
widgetsize(w::PyObject) = (sz = w[:size](); (sz[:width](), sz[:height]))
movewidget(w::PyObject, x::Int, y::Int) = w[:move](x,y)
movewidget(w::PyObject, pos::P2) = movewidget(w, map(roundint,pos)...)
resizewidget(w::PyObject, width::Int, height::Int) = w[:resize](width, height)
resizewidget(w::PyObject, sz::P2) = resizewidget(w, map(roundint,sz)...)
move_resizewidget(w::PyObject, x::Int, y::Int, width::Int, height::Int) = (movewidget(w, x, y); resizewidget(w, width, height))
move_resizewidget(w::PyObject, pos::P2, sz::P2) = (movewidget(w, pos); resizewidget(w, sz))
savepng(w::PyObject, filename::AbstractString) = QT.QPixmap()[:grabWidget](w)[:save](filename, "PNG")
windowtitle(w::PyObject, title::AbstractString) = (w[:setWindowTitle](title); nothing)
Base.close(w::PyObject) = hidewidget(w)  # TODO: clean up python objects properly
moveToScreen(w::PyObject, screenNum::Int = 1) = movewidget(w, screenPosition(screenNum))
moveToLastScreen(w::PyObject) = moveToScreen(w, screenCount())

# methods for the Julia widget (generally just pass to the python version)
showwidget(w::Widget) = showwidget(w.widget)
hidewidget(w::Widget) = hidewidget(w.widget)
widgetpos(w::Widget) = widgetpos(w.widget)
widgetsize(w::Widget) = widgetsize(w.widget)
movewidget(w::Widget, args...) = movewidget(w.widget, args...)
resizewidget(w::Widget, args...) = resizewidget(w.widget, args...)
move_resizewidget(w::Widget, args...) = move_resizewidget(w.widget, args...)
savepng(w::Widget, args...) = savepng(w.widget, args...)
windowtitle(w::Widget, args...) = windowtitle(w.widget, args...)
Base.close(w::Widget) = close(w.widget)
moveToScreen(w::Widget, screenNum::Int = 1) = moveToScreen(w.widget, screenNum)
moveToLastScreen(w::Widget) = moveToLastScreen(w.widget)

# methods for desktop info
desktop() = QAPP[:desktop]()
screenCount() = desktop()[:screenCount]()
screenGeometry(screenNum::Int = 1) = desktop()[:screenGeometry](screenNum)
screenPosition(screenNum::Int = 1) = position(screenGeometry(screenNum))
screenSize(screenNum::Int = 1) = size(screenGeometry(screenNum))


text:
"""
    HoD <: Transform

Get the hour of day corresponding to the data.
"""
struct HoD <: Transform end
cardinality(::HoD) = OneToOne()

_apply(x, ::HoD; kwargs...) = hour.(x)


text:
module HDiscontinuousGalerkin

using Tensors
using FastGaussQuadrature
using SparseArrays
using LinearAlgebra

import Base:@propagate_inbounds

abstract type AbstractRefShape end
abstract type AbstractQuadratureRule end

abstract type Interpolation{dim,shape,order} end
abstract type DiscreteFunctionSpace{dim,T,FE} end

include("utils.jl")
include("shapes.jl")
include("mesh.jl")
include("generate_mesh.jl")
include("triangle_mesh.jl")
include("basis.jl")
include("FiniteElement.jl")
include("LagrangeFE.jl")
include("quadrature.jl")
include("GrundmannMoellerQuad.jl")
include("StrangQuad.jl")
include("ScalarFunctionSpaces.jl")
include("VectorFunctionSpaces.jl")
include("TraceFunctionSpaces.jl")
#include("ParametricFunctionSpaces.jl")
include("assembler.jl")
include("DiscreteFunctions.jl")
include("dofhandler.jl")
include("boundary.jl")
include("iterator.jl")

# Function exports
# mesh
export parse_mesh_triangle, cell_diameter, get_coordinates
export nodes, faces, cells, getcells, get_faces, getnodes, node
export Cell, Node, PolygonalMesh
export getncells, getnfaces, n_faces_per_cell, get_normal
export face_orientation, getfaceset
export getnodeset, getnodesets
export rectangle_mesh
export TriangleCell, RectangleCell
export reference_edge_nodes, getnnodes, get_cell_name
export get_vertices_matrix, getcells_matrix
export n_nodes_per_cell, get_cell_coordinates!
export cell_centroid, cell_volume

# Boundaries
export Dirichlet
export apply!

# Shapes
export get_nodal_points, volume
export get_num_faces
export get_num_vertices

# Quadratures
export QuadratureRule
export RefTetrahedron
export DefaultQuad, GrundmannMoeller, Strang, GaussLegendre
export getpoints, getweights, integrate

# basis
export Dubiner, Lagrange, jacobi, Legendre
export dubiner_basis, ∇dubiner_basis
export value, derivative, gradient_value
export getnbasefunctions
export get_default_geom_interpolator
export getorder, getlowerdiminterpol
export gettopology, get_interpolation

# finite Elements
export ContinuousLagrange, GenericFiniteElement

#utils
export get_affine_map

#Assembler
export start_assemble, assemble!, end_assemble

#FunctionSpaces
export VectorFunctionSpace, ScalarFunctionSpace, ScalarTraceFunctionSpace
export getnquadpoints, getdetJdV, shape_value
export shape_gradient, shape_divergence
export getnfacequadpoints, getfacedetJdS
export face_shape_value
export spatial_coordinate, reference_coordinate
export getnlocaldofs
export reinit!

#Discrete Functions
export function_value
export TrialFunction
export errornorm
export function_value, nodal_avg

#Handlers
export DofHandler, ndofs, ndofs_per_cell, celldofs!
export create_sparsity_pattern, reconstruct!
export dof_range

#ITerator
export CellIterator

end # module


text:
# Unit tests for judiComposite
# Mathias Louboutin (mlouboutin3@gatech.edu)
# July 2021

# number of sources/receivers
nsrc = 2
nrec = 120
ns = 251

ftol = 1f-6

@testset "judiVStack Unit Tests with $(nsrc) sources" for nsrc=[1, 2]
    @timeit TIMEROUTPUT "judiVStack (nsrc=$(nsrc))" begin
        # set up judiVector from array
        n = (10, 10)   # (x,y,z) or (x,z)
        dsize = (nsrc*nrec*ns, 1)
        rec_geometry = example_rec_geometry(nsrc=nsrc, nrec=nrec)
        data = randn(Float32, ns, nrec)
        d_obs = judiVector(rec_geometry, data)
        w0 = judiWeights([randn(Float32,n) for i = 1:nsrc])

        # Composite objs
        c1 = [d_obs; w0]
        c1_b = deepcopy(c1)
        c1_z = similar(c1)
        c2 = [w0; d_obs]
        c2_b = deepcopy(c2)

        @test firstindex(c1) == 1
        @test lastindex(c1) == 2
        @test isapprox(length(c1), length(d_obs) + length(w0))
        @test eltype(c1) == Float32
        @test isfinite(c1)

        @test isapprox(c1[1], c1.components[1])
        @test isapprox(c1.components[1], d_obs)
        @test isapprox(c2.components[2], d_obs)
        @test isapprox(c2.components[1], w0)
        @test isapprox(c1.components[2], w0)

        @test isapprox(c1.components[1], c1_b.components[1])
        @test isapprox(c1.components[2], c1_b.components[2])
        @test isapprox(c2.components[1], c2_b.components[1])
        @test isapprox(c2.components[2], c2_b.components[2])

        @test size(c1, 1) == length(d_obs) + length(w0)
        @test size(c2, 1) == length(d_obs) + length(w0)
        @test size(c1, 2) == 1
        @test size(c2, 2) == 1

        # Transpose
        c1_a = transpose(c1)
        @test size(c1_a, 1) == 1
        @test size(c1_a, 2) == length(d_obs) + length(w0)
        @test isapprox(c1_a.components[1],d_obs)
        @test isapprox(c1_a.components[2], w0)

        # Adjoint
        c1_a = adjoint(c1)
        @test size(c1_a, 1) == 1
        @test size(c1_a, 2) == length(d_obs) + length(w0)
        @test isapprox(c1_a.components[1], d_obs)
        @test isapprox(c1_a.components[2], w0)

        # Conj
        c1_a = conj(c1)
        @test size(c1_a, 2) == 1
        @test size(c1_a, 1) == length(d_obs) + length(w0)
        @test isapprox(c1_a.components[1], d_obs)
        @test isapprox(c1_a.components[2], w0)

        # Test arithithmetic
        u =  [d_obs; w0]
        v =  [2f0 * d_obs; w0 + 1f0]
        w =  [d_obs + 1f0; 2f0 * w0]
        a = .5f0 + rand(Float32)
        b = .5f0 + rand(Float32)

        @test isapprox(u + (v + w), (u + v) + w; rtol=ftol)
        @test isapprox(2f0*u, 2f0.*u; rtol=ftol)
        @test isapprox(u + v, v + u; rtol=ftol)
        @test isapprox(u, u + 0; rtol=ftol)
        @test iszero(norm(u + u*(-1)))
        @test isapprox(-u, -1f0 * u; rtol=ftol)
        @test isapprox(a .* (b .* u), (a * b) .* u; rtol=ftol)
        @test isapprox(u, u .* 1; rtol=ftol)
        @test isapprox(a .* (u + v), a .* u + a .* v; rtol=1f-5)
        @test isapprox((a + b) .* v, a .* v + b .* v; rtol=1f-5)

        # Test the norm
        u_scale = deepcopy(u)
        @test isapprox(norm(u_scale, 2), sqrt(norm(d_obs, 2)^2 + norm(w0, 2)^2))
        @test isapprox(norm(u_scale, 2), sqrt(dot(u_scale, u_scale)); rtol=1f-6)
        @test isapprox(norm(u_scale, 1), norm(d_obs, 1) + norm(w0, 1))
        @test isapprox(norm(u_scale, Inf), max(norm(d_obs, Inf), norm(w0, Inf)))
        @test isapprox(norm(u_scale - 1f0, 1), norm(u_scale .- 1f0, 1))
        @test isapprox(norm(1f0 - u_scale, 1), norm(1f0 .- u_scale, 1))
        @test isapprox(norm(u_scale/2f0, 1), norm(u_scale, 1)/2f0)
        # Test broadcasting
        u_scale = deepcopy(u)
        v_scale = deepcopy(v)

        u_scale .*= 2f0
        @test isapprox(u_scale, 2f0 * u)
        v_scale .+= 2f0
        @test isapprox(v_scale, 2f0 + v)
        u_scale ./= 2f0
        @test isapprox(u_scale, u)
        u_scale .= 2f0 .* u_scale .+ v_scale
        @test isapprox(u_scale, 2f0 * u + (2f0 + v))
        u_scale .= u .+ v
        @test isapprox(u_scale, u + v)
        u_scale .= u .- v
        @test isapprox(u_scale, u - v)
        u_scale .= u .* v
        @test isapprox(u_scale[1], u[1].*v[1])
        u_scale .= u ./ v
        @test isapprox(u_scale[1], u[1]./v[1])

        # Test multi-vstack
        u1 = [u; 2f0*d_obs]
        u2 = [3f0*d_obs; u]

        @test size(u2) == size(u1)
        @test size(u2, 1) == u.m + length(d_obs)

        @test isapprox(u1.components[1], d_obs)
        @test isapprox(u1.components[2], w0)
        @test isapprox(u1.components[3], 2f0*d_obs)

        @test isapprox(u2.components[1], 3f0*d_obs)
        @test isapprox(u2.components[2], d_obs)
        @test isapprox(u2.components[3], w0)

        v1 = [u; 2f0 * w0]
        v2 = [3f0 * w0; u]

        @test size(v2) == size(v1)
        @test size(v2, 1) == u.m + length(w0)

        @test isapprox(v1.components[1], d_obs)
        @test isapprox(v1.components[2], w0)
        @test isapprox(v1.components[3], 2f0*w0)

        @test isapprox(v2.components[1], 3f0*w0)
        @test isapprox(v2.components[2], d_obs)
        @test isapprox(v2.components[3], w0)

        w1 = [c1; 2f0.*c2]
        w2 = [c2./2f0; c1]

        @test size(w1) == size(w2)
        @test size(w2, 1) == c1.m + c2.m

        @test isapprox(w1.components[1], d_obs)
        @test isapprox(w1.components[2], w0)
        @test isapprox(w1.components[3], 2f0*w0)
        @test isapprox(w1.components[4], 2f0*d_obs)

        @test isapprox(w2.components[1], w0 / 2f0)
        @test isapprox(w2.components[2], d_obs / 2f0)
        @test isapprox(w2.components[3], d_obs)
        @test isapprox(w2.components[4], w0)

        @test isapprox(w2.components[1], w0 / 2f0)
        @test isapprox(w2.components[2], d_obs / 2f0)
        @test isapprox(w2.components[3], d_obs)
        @test isapprox(w2.components[4], w0)

        # Test joDirac pertains judiWeights structure

        I = joDirac(nsrc*prod(n),DDT=Float32,RDT=Float32)
        @test isapprox(I*w0, w0)
        lambda = randn(Float32)
        @test isapprox(lambda*I*w0, lambda*w0)
        @test isapprox(I'*w0, w0)
        @test isapprox((lambda*I)'*w0, lambda * w0)

        # Test Forward and Adjoint joCoreBlock * judiVStack
        J = joOnes(nsrc*prod(n),DDT=Float32,RDT=Float32)

        a = [I;J]*w0
        b = [w0; J*w0]
        @test isapprox(a[1], b[1])
        @test isapprox(a[2:end], b[2:end])
    end
end


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ f38350b0-b0e9-11eb-14c1-0b087fed8be4
md"""
# [Day 22: Wizard Simulator 20XX](https://adventofcode.com/2015/day/22)

Little Henry Case decides that defeating bosses with [swords and stuff](https://adventofcode.com/2015/day/21) is boring. Now he's playing the game with a wizard. Of course, he gets stuck on another boss and needs your help again.

In this version, combat still proceeds with the player and the boss taking alternating turns. The player still goes first. Now, however, you don't get any equipment; instead, you must choose one of your spells to cast. The first character at or below `0` hit points loses.

Since you're a wizard, you don't get to wear armor, and you can't attack normally. However, since you do magic damage, your opponent's armor is ignored, and so the boss effectively has zero armor as well. As before, if armor (from a spell, in this case) would reduce damage below `1`, it becomes `1` instead - that is, the boss' attacks always deal at least `1` damage.

On each of your turns, you must select one of your spells to cast. If you cannot afford to cast any spell, you lose. Spells cost mana; you start with 500 mana, but have no maximum limit. You must have enough mana to cast a spell, and its cost is immediately deducted when you cast it. Your spells are Magic Missile, Drain, Shield, Poison, and Recharge.

-   Magic Missile costs `53` mana. It instantly does `4` damage.
-   Drain costs `73` mana. It instantly does `2` damage and heals you for `2` hit points.
-   Shield costs `113` mana. It starts an effect that lasts for `6` turns. While it is active, your armor is increased by `7`.
-   Poison costs `173` mana. It starts an effect that lasts for `6` turns. At the start of each turn while it is active, it deals the boss `3` damage.
-   Recharge costs `229` mana. It starts an effect that lasts for `5` turns. At the start of each turn while it is active, it gives you `101` new mana.

Effects all work the same way. Effects apply at the start of both the player's turns and the boss' turns. Effects are created with a timer (the number of turns they last); at the start of each turn, after they apply any effect they have, their timer is decreased by one. If this decreases the timer to zero, the effect ends. You cannot cast a spell that would start an effect which is already active. However, effects can be started on the same turn they end.

For example, suppose the player has `10` hit points and `250` mana, and that the boss has `13` hit points and `8` damage:

```
-- Player turn --
- Player has 10 hit points, 0 armor, 250 mana
- Boss has 13 hit points
Player casts Poison.

-- Boss turn --
- Player has 10 hit points, 0 armor, 77 mana
- Boss has 13 hit points
Poison deals 3 damage; its timer is now 5.
Boss attacks for 8 damage.

-- Player turn --
- Player has 2 hit points, 0 armor, 77 mana
- Boss has 10 hit points
Poison deals 3 damage; its timer is now 4.
Player casts Magic Missile, dealing 4 damage.

-- Boss turn --
- Player has 2 hit points, 0 armor, 24 mana
- Boss has 3 hit points
Poison deals 3 damage. This kills the boss, and the player wins.

```

Now, suppose the same initial conditions, except that the boss has `14` hit points instead:

```
-- Player turn --
- Player has 10 hit points, 0 armor, 250 mana
- Boss has 14 hit points
Player casts Recharge.

-- Boss turn --
- Player has 10 hit points, 0 armor, 21 mana
- Boss has 14 hit points
Recharge provides 101 mana; its timer is now 4.
Boss attacks for 8 damage!

-- Player turn --
- Player has 2 hit points, 0 armor, 122 mana
- Boss has 14 hit points
Recharge provides 101 mana; its timer is now 3.
Player casts Shield, increasing armor by 7.

-- Boss turn --
- Player has 2 hit points, 7 armor, 110 mana
- Boss has 14 hit points
Shield's timer is now 5.
Recharge provides 101 mana; its timer is now 2.
Boss attacks for 8 - 7 = 1 damage!

-- Player turn --
- Player has 1 hit point, 7 armor, 211 mana
- Boss has 14 hit points
Shield's timer is now 4.
Recharge provides 101 mana; its timer is now 1.
Player casts Drain, dealing 2 damage, and healing 2 hit points.

-- Boss turn --
- Player has 3 hit points, 7 armor, 239 mana
- Boss has 12 hit points
Shield's timer is now 3.
Recharge provides 101 mana; its timer is now 0.
Recharge wears off.
Boss attacks for 8 - 7 = 1 damage!

-- Player turn --
- Player has 2 hit points, 7 armor, 340 mana
- Boss has 12 hit points
Shield's timer is now 2.
Player casts Poison.

-- Boss turn --
- Player has 2 hit points, 7 armor, 167 mana
- Boss has 12 hit points
Shield's timer is now 1.
Poison deals 3 damage; its timer is now 5.
Boss attacks for 8 - 7 = 1 damage!

-- Player turn --
- Player has 1 hit point, 7 armor, 167 mana
- Boss has 9 hit points
Shield's timer is now 0.
Shield wears off, decreasing armor by 7.
Poison deals 3 damage; its timer is now 4.
Player casts Magic Missile, dealing 4 damage.

-- Boss turn --
- Player has 1 hit point, 0 armor, 114 mana
- Boss has 2 hit points
Poison deals 3 damage. This kills the boss, and the player wins.

```

You start with 50 hit points and 500 mana points. The boss's actual stats are in your puzzle input. **What is the least amount of mana you can spend and still win the fight?** (Do not include mana recharge effects as "spending" negative mana.)
"""

# ╔═╡ 77d2e090-c3da-46cd-b6e5-a3cb403d73e0
puzzle_input = "Hit Points: 58
Damage: 9"

# ╔═╡ Cell order:
# ╟─f38350b0-b0e9-11eb-14c1-0b087fed8be4
# ╟─77d2e090-c3da-46cd-b6e5-a3cb403d73e0


text:
"""
    FreeBody{R,T}

A rigid body moving freely in 3-dimensional space, with direct control of both force and torque
at the center of mass, for a total of 6 controls.

# Constructors
    FreeBody{R,T}(mass, J)
    FreeBody(R=UnitQuaternion{Float64}; kwargs...)

where `R <: Rotation{3}` with keyword arguments
* `mass` - mass of the body
* `J` - a vector of the three diagonal elements of the inertia matrix
"""
struct FreeBody{R,T} <: RigidBody{R}
    mass::T
    J::Diagonal{T,SVector{3,T}}
    Jinv::Diagonal{T,SVector{3,T}}
    function FreeBody{R,T}(mass::Real, Jdiag::AbstractVector) where {R,T}
        @assert length(Jdiag) == 3 "inertia diagonal must have only 3 entries"
        @assert mass > 0 "mass must be positive"
        J = SVector{3}(Jdiag)
        new{R,T}(mass, J, inv(J))
    end
end
RobotDynamics.control_dim(::FreeBody) = 6

FreeBody(R=UnitQuaternion{Float64}; mass=1.0, J=@SVector ones(3)) =
    FreeBody{R,promote_type(typeof(mass), eltype(J))}(mass, J)

function forces(model::FreeBody, x::StaticVector, u::StaticVector)
    q = orientation(model, x)
    F = @SVector [u[1], u[2], u[3]]
    q*F  # world frame
end

function moments(model::FreeBody, x::StaticVector, u::StaticVector)
    return @SVector [u[4], u[5], u[6]]  # body frame
end

function RobotDynamics.wrench_jacobian!(F, model::FreeBody, z)
    x = state(z)
    u = control(z)
    q = orientation(model, x)
    ir, iq, iv, iω, iu = RobotDynamics.gen_inds(model)
    iF = SA[1,2,3]
    iM = SA[4,5,6]
    F[iF, iq] .= Rotations.∇rotate(q, u[iF])
    F[iF, iu[iF]] .= RotMatrix(q)
    for i = 1:3
        F[iM[i], iu[i+3]] = 1
    end
    return F
end

function RobotDynamics.wrench_sparsity(::FreeBody)
    SA[false true  false false true;
       false false false false true]
end

inertia(model::FreeBody, x, u) = model.J
inertia_inv(model::FreeBody, x, u) = model.Jinv
mass(model::FreeBody, x, u) = model.mass


text:
function generate_decomposition!(db::SQLite.DB, id::Int, name::AbstractString, scenario::Union{Int, AbstractString},
                                 cliques_path::AbstractString, cliquetrees_path::AbstractString, graphs_path::AbstractString,
                                 extension_alg::AbstractString, option::AbstractDict,
                                 preprocess_path::AbstractString, graph_path::AbstractString;
                                 seed, rng, kwargs...)
    println("Generating decomposition: ($name $scenario) ($extension_alg)")
    g = loadgraph(graph_path)

    # Preprocessing
    option = deepcopy(option)
    option[:nb_edges_to_add] != 0 && add_edges!(g, pop!(option, :nb_edges_to_add), pop!(option, :how); option...)

    # Chordal extension
    chordal_g, data = chordal_extension(g, extension_alg; kwargs...)

    # Features extraction
    features = get_features_graph(chordal_g)
    features["nb_added_edge_dec"] = ne(chordal_g) - ne(g)
    clique = maximal_cliques(chordal_g)
    cliquetree, nb_lc = get_cliquetree_array(clique)
    merge!(features, get_clique_features(clique))

    # Save cliques and cliquetree
    uuid = uuid1(rng)
    clique_path = joinpath(cliques_path, "$(name)_$(scenario)_$(uuid)_clique.csv")
    cliquetree_path = joinpath(cliquetrees_path, "$(name)_$(scenario)_$(uuid)_cliquetree.csv")
    graph_path_dec = joinpath(graphs_path, "$(name)_$(scenario)_$(uuid)_graphs.lgz")
    save_clique(clique, clique_path)
    # Case when the graph is complete (1 clique)
    if isempty(cliquetree)
        cliquetree = [[1, 1]]
    end
    save_cliquetree(cliquetree, cliquetree_path)
    savegraph(graph_path_dec, chordal_g)

    # Other columns
    date = Dates.now()

    features = Dict(Symbol(k) => v for (k, v) in features)
    insert_decomposition!(db, id, uuid, name, scenario, extension_alg, preprocess_path, date, clique_path, cliquetree_path, graph_path_dec; features...)
end

function generate_decomposition_dfrow!(db::SQLite.DB, row, cliques_path::AbstractString, cliquetrees_path::AbstractString,
                                       graphs_path::AbstractString,
                                       extension_alg::AbstractString, option::AbstractDict, preprocess_path::AbstractString;
                                       seed, rng, kwargs...)
    generate_decomposition!(db, row[:id], row[:name], row[:scenario], cliques_path, cliquetrees_path, graphs_path,
                            extension_alg, option,
                            preprocess_path, row[:graph_path];
                            seed=seed, rng=rng, kwargs...)
end

function generate_decompositions!(db::SQLite.DB;
                                  cliques_path, cliquetrees_path, graphs_path,
                                  extension_alg::AbstractString, preprocess_path::AbstractString,
                                  min_nv=typemin(Int), max_nv=typemax(Int), subset=nothing,
                                  seed=MersenneTwister(42), rng=MersenneTwister(42),
                                  kwargs...)
    println("Generate decompositions config: extension_alg=$extension_alg, preprocess_path=$preprocess_path, min_nv=$min_nv, max_nv=$max_nv, seed=$seed, rng=$rng")
    println("subset\n$subset\nend subset")
    !isdir(cliques_path) && mkpath(cliques_path)
    !isdir(cliquetrees_path) && mkpath(cliquetrees_path)
    query = "SELECT id, name, scenario, graph_path FROM instances WHERE nbus >= $min_nv AND nbus <= $max_nv"
    if !isnothing(subset)
        query *= " AND id IN ($(join(subset, ',')))"
    end
    option = JSON.parse(read(open(preprocess_path, "r"), String))
    option = Dict(Symbol(k) => v for (k, v) in option)
    results = DBInterface.execute(db, query) |> DataFrame
    generate_function!(row) = generate_decomposition_dfrow!(db, row, cliques_path, cliquetrees_path, graphs_path,
                                                            extension_alg, option, preprocess_path;
                                                            seed=seed, rng=rng, kwargs...)
    generate_function!.(eachrow(results[!, [:id, :name, :scenario, :graph_path]]))
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using Unums
using Base.Test

@unumbers

a = Unum{3,5}(0x0000000000000040, 0xFFFFFFFE00000000, 0x0002, 0x0007, 0x001E)
b = zero(Unum{3,5})
Unums.sum_inexact(a, b, -63, -63)

a = Unum{3,5}(0x0000000000000000, 0x0000000000000000, 0x0003, 0x0006, 0x001E)
b = Unum{3,5}(0x0000000000000040, 0xFFFFFFFE00000000, 0x0001, 0x0007, 0x001E)
a + b

@test 1 == 0

include("./test-infrastructure.jl")
include("./test-operations.jl")
include("./test-slice.jl")
include("./test-warlpiri.jl")

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl