text:
using Base.Test, StatsFuns, DataFrames, GLM

function test_show(x)
    io = IOBuffer()
    show(io, x)
end

const glm_datadir = joinpath(dirname(@__FILE__), "..", "data")

## Formaldehyde data from the R Datasets package
form = DataFrame(Any[[0.1,0.3,0.5,0.6,0.7,0.9],[0.086,0.269,0.446,0.538,0.626,0.782]],
    [:Carb, :OptDen])

@testset "lm" begin
    lm1 = fit(LinearModel, @formula(OptDen ~ Carb), form)
    test_show(lm1)
    @test isapprox(coef(lm1), collect(linreg(form[:Carb], form[:OptDen])))
    Σ = [6.136653061224592e-05 -9.464489795918525e-05
        -9.464489795918525e-05 1.831836734693908e-04]
    @test isapprox(vcov(lm1), Σ)
    @test isapprox(cor(lm1.model), diagm(diag(Σ))^(-1/2)*Σ*diagm(diag(Σ))^(-1/2))
    @test dof(lm1) == 3
    @test isapprox(deviance(lm1), 0.0002992000000000012)
    @test isapprox(loglikelihood(lm1), 21.204842144047973)
    @test isapprox(nulldeviance(lm1), 0.3138488333333334)
    @test isapprox(nullloglikelihood(lm1), 0.33817870295676444)
    @test r²(lm1) == r2(lm1)
    @test isapprox(r²(lm1), 0.9990466748057584)
    @test adjr²(lm1) == adjr2(lm1)
    @test isapprox(adjr²(lm1), 0.998808343507198)
    @test isapprox(aic(lm1), -36.409684288095946)
    @test isapprox(aicc(lm1), -24.409684288095946)
    @test isapprox(bic(lm1), -37.03440588041178)
end

dobson = DataFrame(Counts = [18.,17,15,20,10,20,25,13,12],
    Outcome = pool(repeat(["A", "B", "C"], outer = 3)),
    Treatment = pool(repeat(["a","b", "c"], inner = 3)))

@testset "Poisson GLM" begin
    gm1 = fit(GeneralizedLinearModel, @formula(Counts ~ Outcome + Treatment), dobson, Poisson())
    test_show(gm1)
    @test dof(gm1) == 5
    @test isapprox(deviance(gm1), 5.12914107700115, rtol = 1e-7)
    @test isapprox(loglikelihood(gm1), -23.380659200978837, rtol = 1e-7)
    @test isapprox(aic(gm1), 56.76131840195767)
    @test isapprox(aicc(gm1), 76.76131840195768)
    @test isapprox(bic(gm1), 57.74744128863877)
    @test isapprox(coef(gm1)[1:3],
        [3.044522437723423,-0.45425527227759555,-0.29298712468147375])
end

## Example from http://www.ats.ucla.edu/stat/r/dae/logit.htm
admit = readtable(joinpath(glm_datadir, "admit.csv.gz"))
admit[:rank] = pool(admit[:rank])

@testset "Binomial, Bernoulli, LogitLink" begin
    for distr in (Binomial, Bernoulli)
        gm2 = fit(GeneralizedLinearModel, @formula(admit ~ gre + gpa + rank), admit, distr())
        test_show(gm2)
        @test dof(gm2) == 6
        @test isapprox(deviance(gm2), 458.5174924758994)
        @test isapprox(loglikelihood(gm2), -229.25874623794968)
        @test isapprox(aic(gm2), 470.51749247589936)
        @test isapprox(aicc(gm2), 470.7312329339146)
        @test isapprox(bic(gm2), 494.4662797585473)
        @test isapprox(coef(gm2),
            [-3.9899786606380734, 0.0022644256521549043, 0.8040374535155766,
            -0.6754428594116577, -1.3402038117481079,-1.5514636444657492])
    end
end

@testset "Bernoulli ProbitLink" begin
    gm3 = fit(GeneralizedLinearModel, @formula(admit ~ gre + gpa + rank), admit,
        Binomial(), ProbitLink())
    test_show(gm3)
    @test dof(gm3) == 6
    @test isapprox(deviance(gm3), 458.4131713833386)
    @test isapprox(loglikelihood(gm3), -229.20658569166932)
    @test isapprox(aic(gm3), 470.41317138333864)
    @test isapprox(aicc(gm3), 470.6269118413539)
    @test isapprox(bic(gm3), 494.36195866598655)
    @test isapprox(coef(gm3),
        [-2.3867922998680786, 0.0013755394922972369, 0.47772908362647015,
        -0.4154125854823675, -0.8121458010130356, -0.9359047862425298])
end

@testset "Bernoulli CauchitLink" begin
    gm4 = fit(GeneralizedLinearModel, @formula(admit ~ gre + gpa + rank), admit,
        Binomial(), CauchitLink())
    test_show(gm4)
    @test dof(gm4) == 6
    @test isapprox(deviance(gm4), 459.3401112751141)
    @test isapprox(loglikelihood(gm4), -229.6700556375571)
    @test isapprox(aic(gm4), 471.3401112751142)
    @test isapprox(aicc(gm4), 471.5538517331295)
    @test isapprox(bic(gm4), 495.28889855776214)
end

@testset "Bernoulli CloglogLink" begin
    gm5 = fit(GeneralizedLinearModel, @formula(admit ~ gre + gpa + rank), admit,
        Binomial(), CloglogLink())
    test_show(gm5)
    @test dof(gm5) == 6
    @test isapprox(deviance(gm5), 458.89439629612616)
    @test isapprox(loglikelihood(gm5), -229.44719814806314)
    @test isapprox(aic(gm5), 470.8943962961263)
    @test isapprox(aicc(gm5), 471.1081367541415)
    @test isapprox(bic(gm5), 494.8431835787742)
end

## Example with offsets from Venables & Ripley (2002, p.189)
anorexia = readtable(joinpath(glm_datadir, "anorexia.csv.gz"))
anorexia[:Treat] = pool(anorexia[:Treat])

@testset "Offset" begin
    gm6 = fit(GeneralizedLinearModel, @formula(Postwt ~ Prewt + Treat), anorexia,
        Normal(), IdentityLink(), offset=Array(anorexia[:Prewt]))
    test_show(gm6)
    @test dof(gm6) == 5
    @test isapprox(deviance(gm6), 3311.262619919613)
    @test isapprox(loglikelihood(gm6), -239.9866487711122)
    @test isapprox(aic(gm6), 489.9732975422244)
    @test isapprox(aicc(gm6), 490.8823884513153)
    @test isapprox(bic(gm6), 501.35662813730465)
    @test isapprox(coef(gm6),
        [49.7711090149846,-0.5655388496391,-4.0970655280729,4.5630626529188])
    @test isapprox(GLM.dispersion(gm6.model, true), 48.6950385282296)
    @test isapprox(stderr(gm6),
        [13.3909581420259, 0.1611823618518, 1.8934926069669, 2.1333359226431])
end

@testset "Normal LogLink offset" begin
    gm7 = fit(GeneralizedLinearModel, @formula(Postwt ~ Prewt + Treat), anorexia,
        Normal(), LogLink(), offset=Array(anorexia[:Prewt]), convTol=1e-8)
    test_show(gm7)
    @test isapprox(deviance(gm7), 3265.207242977156)
    @test isapprox(coef(gm7),
        [3.992326787835955, -0.994452693131178, -0.050698258703974, 0.051494029957641])
    @test isapprox(GLM.dispersion(gm7.model, true), 48.017753573192266)
    @test isapprox(stderr(gm7), [0.15716774, 0.0018862835, 0.02258404, 0.023882795], atol=1e-6)
end

## Gamma example from McCullagh & Nelder (1989, pp. 300-2)
clotting = DataFrame(u = log.([5,10,15,20,30,40,60,80,100]),
                     lot1 = [118,58,42,35,27,25,21,19,18])

@testset "Gamma InverseLink" begin
    gm8 = fit(GeneralizedLinearModel, @formula(lot1 ~ 1 + u), clotting, Gamma())
    test_show(gm8)
    @test dof(gm8) == 3
    @test isapprox(deviance(gm8), 0.016729715178484157)
    @test isapprox(loglikelihood(gm8), -15.994961974777247)
    @test isapprox(aic(gm8), 37.989923949554495)
    @test isapprox(aicc(gm8), 42.78992394955449)
    @test isapprox(bic(gm8), 38.58159768156315)
    @test isapprox(coef(gm8), [-0.01655438172784895,0.01534311491072141])
    @test isapprox(GLM.dispersion(gm8.model, true), 0.002446059333495581, atol=1e-6)
    @test isapprox(stderr(gm8), [0.00092754223, 0.000414957683], atol=1e-6)
end

@testset "Gamma LogLink" begin
    gm9 = fit(GeneralizedLinearModel, @formula(lot1 ~ 1 + u), clotting, Gamma(), LogLink(),
        convTol=1e-8)
    test_show(gm9)
    @test dof(gm9) == 3
    @test isapprox(deviance(gm9), 0.16260829451739)
    @test isapprox(loglikelihood(gm9), -26.24082810384911)
    @test isapprox(aic(gm9), 58.48165620769822)
    @test isapprox(aicc(gm9), 63.28165620769822)
    @test isapprox(bic(gm9), 59.07332993970688)
    @test isapprox(coef(gm9), [5.50322528458221, -0.60191617825971])
    @test isapprox(GLM.dispersion(gm9.model, true), 0.02435442293561081)
    @test isapprox(stderr(gm9), [0.19030107482720, 0.05530784660144])
end

@testset "Gamma IdentityLink" begin
    gm10 = fit(GeneralizedLinearModel, @formula(lot1 ~ 1 + u), clotting, Gamma(), IdentityLink(),
        convTol=1e-8)
    test_show(gm10)
    @test dof(gm10) == 3
    @test isapprox(deviance(gm10), 0.60845414895344)
    @test isapprox(loglikelihood(gm10), -32.216072437284176)
    @test isapprox(aic(gm10), 70.43214487456835)
    @test isapprox(aicc(gm10), 75.23214487456835)
    @test isapprox(bic(gm10), 71.02381860657701)
    @test isapprox(coef(gm10), [99.250446880986, -18.374324929002])
    @test isapprox(GLM.dispersion(gm10.model, true), 0.10417373, atol=1e-6)
    @test isapprox(stderr(gm10), [17.864084, 4.297895], atol=1e-4)
end

# Logistic regression using aggregated data and weights
admit_agr = DataFrame(count = [28., 97, 93, 55, 33, 54, 28, 12],
                      admit = repeat([false, true], inner=[4]),
                      rank = pool(repeat([1, 2, 3, 4], outer=[2])))

@testset "Aggregated Binomial LogitLink" begin
    for distr in (Binomial, Bernoulli)
        gm14 = fit(GeneralizedLinearModel, @formula(admit ~ rank), admit_agr, distr(), wts=Array(admit_agr[:count]))
        @test dof(gm14) == 4
        @test nobs(gm14) == 400
        @test isapprox(deviance(gm14), 474.9667184280627)
        @test isapprox(loglikelihood(gm14), -237.48335921403134)
        @test isapprox(aic(gm14), 482.96671842822883)
        @test isapprox(aicc(gm14), 483.0679842510136)
        @test isapprox(bic(gm14), 498.9325766164946)
        @test isapprox(coef(gm14),
            [0.164303051291, -0.7500299832, -1.36469792994, -1.68672866457], atol=1e-5)
    end
end

# Logistic regression using aggregated data with proportions of successes and weights
admit_agr2 = DataFrame(Any[[61., 151, 121, 67], [33., 54, 28, 12], pool([1, 2, 3, 4])],
    [:count, :admit, :rank])
admit_agr2[:p] = admit_agr2[:admit] ./ admit_agr2[:count]

## The model matrix here is singular so tests like the deviance are just round off error
@testset "Binomial LogitLink aggregated" begin
    gm15 = fit(GeneralizedLinearModel, @formula(p ~ rank), admit_agr2, Binomial(), wts=admit_agr2[:count])
    test_show(gm15)
    @test dof(gm15) == 4
    @test nobs(gm15) == 400
# The model matrix is singular so the deviance is essentially round-off error
#    @test isapprox(deviance(gm15), -2.4424906541753456e-15, rtol = 1e-7)
    @test isapprox(loglikelihood(gm15), -9.50254433604239)
    @test isapprox(aic(gm15), 27.00508867208478)
    @test isapprox(aicc(gm15), 27.106354494869592)
    @test isapprox(bic(gm15), 42.970946860516705)
    @test isapprox(coef(gm15),
        [0.1643030512912767, -0.7500299832303851, -1.3646980342693287, -1.6867295867357475])
end

# Weighted Gamma example (weights are totally made up)
@testset "Gamma InverseLink Weights" begin
    gm16 = fit(GeneralizedLinearModel, @formula(lot1 ~ 1 + u), clotting, Gamma(),
        wts=[1.5,2.0,1.1,4.5,2.4,3.5,5.6,5.4,6.7])
    test_show(gm16)
    @test dof(gm16) == 3
    @test nobs(gm16) == 32.7
    @test isapprox(deviance(gm16), 0.03933389380881689)
    @test isapprox(loglikelihood(gm16), -43.35907878769152)
    @test isapprox(aic(gm16), 92.71815757538305)
    @test isapprox(aicc(gm16), 93.55439450918095)
    @test isapprox(bic(gm16), 97.18028280909267)
    @test isapprox(coef(gm16), [-0.017217012615523237, 0.015649040411276433])
end

# Weighted Poisson example (weights are totally made up)
@testset "Poisson LogLink Weights" begin
    gm17 = fit(GeneralizedLinearModel, @formula(Counts ~ Outcome + Treatment), dobson, Poisson(),
        wts = [1.5,2.0,1.1,4.5,2.4,3.5,5.6,5.4,6.7])
    test_show(gm17)
    @test dof(gm17) == 5
    @test isapprox(deviance(gm17), 17.699857821414266)
    @test isapprox(loglikelihood(gm17), -84.57429468506352)
    @test isapprox(aic(gm17), 179.14858937012704)
    @test isapprox(aicc(gm17), 181.39578038136298)
    @test isapprox(bic(gm17), 186.5854647596431)
    @test isapprox(coef(gm17), [3.1218557035404793, -0.5270435906931427,-0.40300384148562746,
                           -0.017850203824417415,-0.03507851122782909])
end

@testset "Sparse GLM" begin
    srand(1)
    X = sprand(1000, 10, 0.01)
    β = randn(10)
    y = Bool[rand() < logistic(x) for x in X * β]
    gmsparse = fit(GeneralizedLinearModel, X, y, Binomial())
    gmdense = fit(GeneralizedLinearModel, full(X), y, Binomial())

    @test isapprox(deviance(gmsparse), deviance(gmdense))
    @test isapprox(coef(gmsparse), coef(gmdense))
    @test isapprox(vcov(gmsparse), vcov(gmdense))
end

@testset "Predict" begin
    srand(1)
    X = rand(10, 2)
    Y = logistic.(X * [3; -3])

    gm11 = fit(GeneralizedLinearModel, X, Y, Binomial())
    @test isapprox(predict(gm11), Y)

    newX = rand(5, 2)
    newY = logistic.(newX * coef(gm11))
    @test isapprox(predict(gm11, newX), newY)

    off = rand(10)
    newoff = rand(5)

    @test_throws ArgumentError predict(gm11, newX, offset=newoff)

    gm12 = fit(GeneralizedLinearModel, X, Y, Binomial(), offset=off)
    @test_throws ArgumentError predict(gm12, newX)
    @test isapprox(predict(gm12, newX, offset=newoff),
        logistic.(newX * coef(gm12) .+ newoff))

        # Prediction from DataFrames
    d = convert(DataFrame, X)
    d[:y] = Y

    gm13 = fit(GeneralizedLinearModel, @formula(y ~ 0 + x1 + x2), d, Binomial())
    @test predict(gm13) == predict(gm13, d[[:x1, :x2]])
    @test predict(gm13) == predict(gm13, d)

    newd = convert(DataFrame, newX)
    predict(gm13, newd)

    Ylm = X * [0.8, 1.6] + 0.8randn(10)
    mm = fit(LinearModel, X, Ylm)
    pred = predict(mm, newX, :confint)

    @test isapprox(pred[1,2], 0.6122189104014528)
    @test isapprox(pred[2,2], -0.33530477814532056)
    @test isapprox(pred[3,2], 1.340413688904295)
    @test isapprox(pred[4,2], 0.02118806218116165)
    @test isapprox(pred[5,2], 0.8543142404183606)
    @test isapprox(pred[1,3], 2.6853964084909836)
    @test isapprox(pred[2,3], 1.2766396685055916)
    @test isapprox(pred[3,3], 3.6617479283005894)
    @test isapprox(pred[4,3], 0.6477623101170038)
    @test isapprox(pred[5,3], 2.564532433982956)
end

@testset "Issue 118" begin
    @inferred nobs(lm(randn(10, 2), randn(10)))
end

@testset "Issue 84" begin
    X = [1 1; 2 4; 3 9]
    Xf = [1 1; 2 4; 3 9.]
    y = [2, 6, 12]
    yf = [2, 6, 12.]
    @test isapprox(lm(X, y).pp.beta0, ones(2))
    @test isapprox(lm(Xf, y).pp.beta0, ones(2))
    @test isapprox(lm(X, yf).pp.beta0, ones(2))
end

@testset "Issue 153" begin
    X = [ones(10) randn(10)]
    Test.@inferred cholfact(DensePredQR{Float64}(X))
end

@testset "Issue 117" begin
    data = DataFrame(x = [1,2,3,4], y = [24,34,44,54])
    @test coef(glm(@formula(y ~ x), data, Normal(), IdentityLink())) == [14;10]
end


text:
using Symbolics, Test
import Symbolics: latex

@testset "LaTeX rendering" begin
    @sym x y z

    Symbolics.abbreviate_power_of_negative_unity(true)
    @test latex((-1)^x) == ("(-)^{x}", 1)

    Symbolics.abbreviate_power_of_negative_unity(false)
    @test latex((-1)^x) == ("\\left(-1\\right)^{x}", 1)

    Symbolics.abbreviate_power_of_negative_unity(true)
    f=sin(x^3)
    @test latex(f) == ("\\sin\\left(x^{3}\\right)", 1)
    g′=x/(-z)
    @test latex(g′) == ("\\frac{x}{-z}", 0)
    g=x/-sin(z)
    @test latex(g) == ("\\frac{x}{-\\left(\\sin{z}\\right)}", 1)
    @test latex(x/((x+y+z)^4)) == ("\\frac{x}{\\left(x+y+z\\right)^{4}}", 1)
    @test latex(x/(sin(z)^4)) == ("\\frac{x}{\\sin^{4}z}", 0)
    @test latex(x/(sin(y+z)^4)) == ("\\frac{x}{\\sin^{4}\\left(y+z\\right)}", 1)
    @test latex(x/(sin(y+z))) == ("\\frac{x}{\\sin\\left(y+z\\right)}", 1)
    @test latex(x/(sin(z)^(y+x))) == ("\\frac{x}{\\sin^{y+x}z}", 0)
    @test latex(x + im) == ("x+\\mathrm{i}", 0)
    h=cos(f+g)
    @test latex(h) == ("\\cos\\left[\\sin\\left(x^{3}\\right)+\\frac{x}{-\\left(\\sin{z}\\right)}\\right]", 2)

    Symbolics.advance_delimiters(false)
    @test latex(cos((-1)^x+sin(h))) == ("\\cos\\left((-)^{x}+\\sin\\left(\\cos\\left(\\sin\\left(x^{3}\\right)+\\frac{x}{-\\left(\\sin{z}\\right)}\\right)\\right)\\right)", 0)

    Symbolics.advance_delimiters(true)
    @test latex(cos((-1)^x+sin(h))) == ("\\cos\\left((-)^{x}+\\sin\\left\\{\\cos\\left[\\sin\\left(x^{3}\\right)+\\frac{x}{-\\left(\\sin{z}\\right)}\\right]\\right\\}\\right)", 4)
    @test latex(cos((-1)^x+f)) == ("\\cos\\left[(-)^{x}+\\sin\\left(x^{3}\\right)\\right]", 2)

    @test latex(cos((-1)^(2x+y*(z+cos(x)))+sin(h))) == ("\\cos\\left((-)^{2x+y\\left(z+\\cos{x}\\right)}+\\sin\\left\\{\\cos\\left[\\sin\\left(x^{3}\\right)+\\frac{x}{-\\left(\\sin{z}\\right)}\\right]\\right\\}\\right)", 4)
    @test latex(cos((-1)^(2x+y*(z+cos(x)))+f)) == ("\\cos\\left[(-)^{2x+y\\left(z+\\cos{x}\\right)}+\\sin\\left(x^{3}\\right)\\right]", 2)
    @test latex(x-y-z+x^(2y+4z*(y+x))) == ("-\\left(y+z\\right)+x+x^{2y+4z\\left(y+x\\right)}", 1)
    @test latex(asinh(x+sin(y))) == ("\\operatorname{arcsinh}\\left(x+\\sin{y}\\right)", 1)
    @test latex(log10(x)+log(x)+log2(x)) == ("\\log_{10}{x}+\\ln{x}+\\log_2{x}", 0)
    @test latex(√x) == ("\\sqrt{x}", 0)
    @test latex(√(x+y)) == ("\\sqrt{x+y}", 0)
    @test latex(sin(x/y)) == ("\\sin\\left(\\frac{x}{y}\\right)", 1)
    @test latex(cos(sin(x/y))) == ("\\cos\\left[\\sin\\left(\\frac{x}{y}\\right)\\right]", 2)
    @test latex(sin(-x)) == ("\\sin\\left(-x\\right)", 1)
    @test latex(exp(-x/y)) == ("\\exp\\left(-\\frac{x}{y}\\right)", 1)
    @test latex(exp(x/-y)) == ("\\exp\\left(\\frac{x}{-y}\\right)", 1)

    @test latex(1e-17x) == ("1.0\\times10^{-17}x",0)
end

ext: jl


text:
using GitLab
using JLD, HttpCommon
using Base.Test

event_request = JLD.load(joinpath(dirname(@__FILE__), "event_request.jld"), "event_request")
event_json = Requests.json(event_request)
event = GitLab.event_from_payload!("Issue", event_json)
## @show event

################
# WebhookEvent #
################

@test get(event.repository.name) == "Calculus"
@test get(event.sender.username) == "mdpradeep"

#################
# EventListener #
#################
## @test !(GitLab.has_valid_secret(event_request, "wrong"))
## @test GitLab.has_valid_secret(event_request, "secret")
@test !(GitLab.is_valid_event(event_request, ["wrong"]))
@test GitLab.is_valid_event(event_request, ["Note Hook"])
@test !(GitLab.from_valid_repo(event, ["Random"]))
@test GitLab.from_valid_repo(event, ["Calculus"])

@test (GitLab.handle_event_request(event_request, x -> true, secret = "secret",
            events = ["Issue"], repos = ["Calculus"])).status == 400

@test begin
    listener = EventListener(x -> true;
                             secret = "secret",
                             repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitLab.jl"],
                             events = ["commit_comment"],
                             forwards = ["http://bob.com", HttpCommon.URI("http://jim.org")])
    r = listener.server.http.handle(HttpCommon.Request(),HttpCommon.Response())
    r.status == 400
end

###################
# CommentListener #
###################

result = GitLab.handle_comment((e, m) -> m, event, GitLab.AnonymousAuth(), r"`runbenchmarks\(.*?\)`", false)

@test result.match == "`runbenchmarks()`"

@test begin
    listener = CommentListener((x, y) -> true, r"trigger";
                               secret = "secret",
                               repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitLab.jl"],
                               forwards = ["http://bob.com", HttpCommon.URI("http://jim.org")],
                               check_collab = false)
    r = listener.listener.server.http.handle(HttpCommon.Request(), HttpCommon.Response())
    r.status == 400
end


text:
using BenchmarkTools
using GeoPhyInv
using Test
#md using Plots
#md gr()

# # Intro

#md # ```@docs
#md # GeoPhyInv.Medium
#md # Medium(::Symbol, ::Real)
#md # ```

# # Examples

# Load a predefined model.
mod=Medium(:marmousi2);

# Get the medium parameters that are stored.
names(mod)

# Look at the reference values.
mod.ref

# Inspect the bounds.
mod.bounds

# Plotting is easy
#md p1=heatmap(mod, :vp)
#md plot(p1,size=(800,400))

# To construct an instance of `Medium`, we need a 2-D grid.
mgrid = [range(0.0, stop=10.,step=0.1), range(0.0, stop=30.,step=0.2)];

# Allocate basic medium parameters on the grid.
mod = Medium(mgrid, [:vp,:rho,:vs]);

# Bounds for these parameters should be input for modeling or inversion. Use `update!`
vpb=[2100.,2200.]; vsb=[1500, 1700]; rhob=[2100., 2300.];
update!(mod, [:vp,:vs,:rho], [vpb, vsb, rhob]);

# Just fill `mod` with average (reference) values.
fill!(mod);

# Once the basic medium parameters are input, we can access some other derived parameters.
mod[:Zp];

# Otherwise, we can manually update parameters of `mod`.
mod[:vp].=3000.;
mod[:vs].=2000.;

# A model can be also be updated by adding random noise.
update!(mod, [:vp,:rho], randn_perc=1.);

# # Methods
#md # ```@docs
#md # Base.getindex(::Medium, ::Symbol)
#md # GeoPhyInv.update!(::GeoPhyInv.Medium, ::Vector{Symbol},)
#md # Base.copyto!(x::AbstractArray, mod::Medium, fields::Vector{Symbol})
#md # Base.vec(mod::Medium, ::Symbol)
#md # ```


text:
using AbstractTrees
using Flux
using Zygote
using Zygote: Grads

"""
    paramsdict(model::Model) -> Dict{Symbol, Array}

Return a Dict with mapping parameter name => parameter value.

Useful for inspecting the distribution of parameters.
"""
function paramsdict(model::Model; filterfn = _ -> true)
    ps = Dict()
    for child::Model in IterModels(model)
        if filterfn(child.layer)
            pslayer = Dict(field => getfield(child.layer, field) for field in child.paramattrs)
            for (name, p) in pslayer
                if !haskey(ps, name)
                    ps[name] = []
                end
                push!(ps[name], p)
            end
        end
    end
    return ps
end

"""
    gradientsdict(model::Model, gs::Grads) -> Dict{Symbol, Array}

Return a Dict with mapping parameter name => gradient value

Useful for inspecting the distribution of gradients.
"""
function gradientsdict(model::Model, gs::Grads; filterfn = _ -> true)
    return Dict(
        name => [gs[p] for p in ps_]
        for (name, ps_) in paramsdict(model, filterfn = filterfn)
    )
end

"""
    layerstats(model, lossfn, batch, [filter_fn])

Return the names, activations and gradients of every layer
where `filter_fn(layer) == true`.

Gradients and activations are calculated on `batch` and with `lossfn`
"""
function layerstats(model, lossfn, batch; filter_fn = l -> true)
    x, y = batch
    hmodel = addhook(model, forward = hooksaveactivations, backward = hooksavegrads, filter_fn = filter_fn)

    gradient(params(hmodel)) do
        return lossfn(hmodel(x), y)
    end

    names = []
    acts = []
    gs = []

    for layer in IterLayers(hmodel)
        if layer isa Hook
            push!(names, string(layer.layer))
            push!(acts, layer.state["activation"])
            push!(gs, layer.state["gradient"])
        end
    end

    return names, acts, gs
end


text:
using UltraDark
using UltraDark: TimeStepOptions
using Test
using MPI

@testset "Actual time step" begin
    time_step_options = TimeStepOptions(update_period=2)
    @test UltraDark.actual_time_step(0.5, 1, time_step_options)[1] ≈ 0.5
    @test UltraDark.actual_time_step(0.5, 1, time_step_options)[2] == 2

    time_step_options = TimeStepOptions(update_period=3)
    @test UltraDark.actual_time_step(0.5, 1, time_step_options)[1] ≈ 0.5
    @test UltraDark.actual_time_step(0.5, 1, time_step_options)[2] == 2

    time_step_options = TimeStepOptions()
    @test UltraDark.actual_time_step(0.01, 1, time_step_options)[1] ≈ 0.01
    @test UltraDark.actual_time_step(0.01, 1, time_step_options)[2] == 10

    time_step_options = TimeStepOptions(multiplier=10.0)
    @test UltraDark.actual_time_step(0.01, 1, time_step_options)[1] ≈ 0.1
    @test UltraDark.actual_time_step(0.01, 1, time_step_options)[2] == 10

    time_step_options = TimeStepOptions(update_period=2)
    @test UltraDark.actual_time_step(0.7, 1, time_step_options)[1] ≈ 0.5
    @test UltraDark.actual_time_step(0.7, 1, time_step_options)[2] == 2

    time_step_options = TimeStepOptions(update_period=3)
    @test UltraDark.actual_time_step(0.7, 2, time_step_options)[1] ≈ 2/3
    @test UltraDark.actual_time_step(0.7, 2, time_step_options)[2] == 3
end

for grid_type in [Grids, PencilGrids]
    @testset "Take n steps, $grid_type" begin
        grids = grid_type(1.0, 16)
        output_config = OutputConfig(mktempdir(), [])
        @test UltraDark.take_steps!(grids, 0, 1, 10, output_config, Config.constant_scale_factor, nothing) == 10.0
    end
end

# Test that evolve to arrives at right time
for grid_type in [Grids, PencilGrids]
    @testset "Evolve from $t_begin to $t_end, $grid_type" for t_begin in [0, 0.123, 1], t_end in [1.23, 2]
        grids = grid_type(1.0, 16)

        a = 1
        grids.ψx .= (grids.x.^2 .+ grids.y.^2 .+ grids.z.^2).^0.5 ./ 1e9 # Set ψx to something non-zero
        grids.ψk .= (grids.fft_plan * grids.ψx)
        grids.ρx .= abs2.(grids.ψx)
        grids.Φk .= -4 * π * (grids.rfft_plan * grids.ρx) ./ (a * grids.rk.^2)
        grids.Φk[1, 1, 1] = 0
        grids.Φx .= grids.rfft_plan \ grids.Φk

        output_config = OutputConfig(mktempdir(), [])
        options = UltraDark.Config.SimulationConfig()
        @test UltraDark.evolve_to!(t_begin, t_end, grids, output_config, options) ≈ t_end
    end
end


text:
using HTTP,JSON

const server_name = "http://titus.phchem.uni-due.de:3000"
const client_id = ""
const client_secret = ""
const redirect_uri = ""
const team_id = "AK Hasselbrink"
const project_id_femto_lab = "9"
const experiment_id_SEC = "53"
const task_id_EC = "297"

function auth_code()
    r = HTTP.request("GET",server_name*"/oauth/authorize?client_id=tim.laemmerzahl@uni-due.de&redirect_uri=$(server_name)&response_type=code")
    String(r.body)

end

authorization_code = auth_code()

"""
token()

Get token with client_id, client_secret, authorization_code, and redirect_uri.
"""
function token()
    resp = HTTP.request("POST", SCINOTE_URL*"/oauth/token",
        ["grant_type" => "authorization_code",
        "client_id" => client_id,
        "client_secret" => client_secret,
        "code" => authorization_code,
        "redirect_uri" => redirect_uri
        ])
    body = String(resp.body)
    JSON.parse(body)
end

"""
refreshtoken()

Refresht the  token with client_id, client_secret, authorization_code, and redirect_uri.
"""
function refreshtoken()
    resp = HTTP.request("POST", SCINOTE_URL*"/oauth/token",
        ["grant_type" => "refresh_token",
        "client_id" => client_id,
        "client_secret" => client_secret,
        "refresh_token" => token,
        "redirect_uri" => redirect_uri
        ])
    body = String(resp.body)
    JSON.parse(body)
end

function api_status()
    r = HTTP.request("GET",server_name*"/api/status")
    JSON.parse(String(r.body))
end

max_stars_repo_name: lambcount/EC

max_issues_repo_name: lambcount/EC

max_forks_repo_name: lambcount/EC


text:
__precompile__()
module ErdosExtras
using Erdos
using JuMP
using GLPK
import BlossomV # matching

LP_OPTIMIZER = GLPK.Optimizer
MIP_OPTIMIZER = GLPK.Optimizer

# export set_lp_solver, set_mip_solver, get_lp_solver, get_mip_solver

# set_lp_solver(solv) = (LP_SOLVER=solv; solv)
# set_mip_solver(solv) = (MIP_SOLVER=solv; solv)
#
# get_lp_solver() = LP_SOLVER
# get_mip_solver() = MIP_SOLVER

include("matching/Matching.jl")
export minimum_weight_perfect_bmatching, minimum_weight_perfect_matching

include("tsp/TSP.jl")
export solve_tsp

end # module


text:
function set_up_operator(p::Param)

    #
    # Variables
    #

    xmin        = p.xmin
    xmax        = p.xmax
    Nx          = p.xnodes

    ymin        = p.ymin
    ymax        = p.ymax
    Ny          = p.ynodes

    #
    # Positional and differential operators
    #

    #Create position vectors
    dx      = (xmax - xmin)/(Nx-1)
    x       = [xmin + (i-1)*dx for i in 1:Nx]

    dy      = (ymax - ymin)/(Ny-1)
    y       = [ymin + (i-1)*dy for i in 1:Ny]

    if(p.Boundaries == :periodic )
        dx      = (xmax - xmin)/(Nx)
        x       = [xmin + (i-1)*dx for i in 1:Nx]

        dy      = (ymax - ymin)/(Ny)
        y       = [ymin + (i-1)*dy for i in 1:Ny]
    end

    #Create x and y matrices
    xmat = x_coord_matrix(p)
    ymat = y_coord_matrix(p)

    #Create Differential Operators
    Dx  =   Diff_Operator_2D(1 ,1, [dx,dy] , [Nx,Ny])
    Dxx =   Diff_Operator_2D(2 ,1, [dx,dy] , [Nx,Ny])

    Dy  =   Diff_Operator_2D(1 ,2, [dx,dy] , [Nx,Ny])
    Dyy =   Diff_Operator_2D(2 ,2, [dx,dy] , [Nx,Ny])

    if(p.Boundaries == :periodic )
        Dx  =   Diff_Operator_2D_Periodic(1 ,1, [dx,dy] , [Nx,Ny])
        Dxx =   Diff_Operator_2D_Periodic(2 ,1, [dx,dy] , [Nx,Ny])

        Dy  =  Diff_Operator_2D_Periodic(1 ,2, [dx,dy] , [Nx,Ny])
        Dyy =  Diff_Operator_2D_Periodic(2 ,2, [dx,dy] , [Nx,Ny])
    end

    #
    # RHS Differential operators (to apply on the field at each iteration)
    #

    # Case with potential operator
    #=

    #Create potential matrix
    # V = potential_matrix(p)

    #RHS of differential equation in non-rotating frame (absoprion done in KG solver)
    RHS_Opp =   (Dxx .+ Dyy) .- (p.μ^2).*sparse(I , Nx*Ny , Nx*Ny) .- V

    =#

    #Case without potential
    RHS_Opp =   (Dxx .+ Dyy) .- (p.μ^2).*sparse(I , Nx*Ny , Nx*Ny)

    #Angular momentum operator
    Lz = p.Ω .* (  (- ymat * Dx) .+ (xmat * Dy )  )

    #
    # Boundary condition operators
    #

    #=
        For the radial case, the border is always set at a small distance
        from the edge. Hence the 0.9 factor bellow.

    =#

    #Boundary Conditions
    BC_Mat  =  sparse(I, Nx*Ny , Ny*Ny)    #Periodic BC (Identity)

    if (p.Boundaries == :radial )
        BC_Mat  =   radial_zero_BC_matrix(p , x , y , 0.95*min(xmax,ymax) )
    elseif (p.Boundaries == :square )
        BC_Mat  =   square_zero_BC_matrix(p)
    end

    #
    # Auxiliary vector (for temporary storage) and absorption matrix
    #

    aux_vector = zeros(Nx*Ny)
    absorption_matrix = spzeros(Nx*Ny , Nx*Ny)

    #Operator structure
    Operators = Opps(x , y , aux_vector , Dx , Dxx , Dy , Dyy , RHS_Opp , BC_Mat , absorption_matrix , Lz)

    #Change absorption matrix according to parameters (at t = 0)
    alpha!(absorption_matrix , p , Operators , 0.0)

    #Return final structure of operators
    return Opps(x , y , aux_vector , Dx , Dxx , Dy , Dyy , RHS_Opp , BC_Mat , absorption_matrix , Lz)

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## Xoshiro RNG
# Lots of implementation is shared with TaskLocalRNG

"""
    Xoshiro

Xoshiro256++ is a fast pseudorandom number generator originally developed by Sebastian Vigna.
Reference implementation is available at http://prng.di.unimi.it

Apart from the high speed, Xoshiro has a small memory footprint, making it suitable for
applications where many different random states need to be held for long time.

Julia's Xoshiro implementation has a bulk-generation mode; this seeds new virtual PRNGs
from the parent, and uses SIMD to generate in parallel (i.e. the bulk stream consists of
multiple interleaved xoshiro instances).
The virtual PRNGs are discarded once the bulk request has been serviced (and should cause
no heap allocations).
"""
mutable struct Xoshiro <: AbstractRNG
    s0::UInt64
    s1::UInt64
    s2::UInt64
    s3::UInt64

    Xoshiro(s0::Integer, s1::Integer, s2::Integer, s3::Integer) = new(s0, s1, s2, s3)
    Xoshiro(seed) = seed!(new(), seed)
end

Xoshiro(::Nothing) = Xoshiro()

function Xoshiro()
    parent = RandomDevice()
    # Constants have nothing up their sleeve, see task.c
    # 0x02011ce34bce797f == hash(UInt(1))|0x01
    # 0x5a94851fb48a6e05 == hash(UInt(2))|0x01
    # 0x3688cf5d48899fa7 == hash(UInt(3))|0x01
    # 0x867b4bb4c42e5661 == hash(UInt(4))|0x01

    Xoshiro(0x02011ce34bce797f * rand(parent, UInt64),
            0x5a94851fb48a6e05 * rand(parent, UInt64),
            0x3688cf5d48899fa7 * rand(parent, UInt64),
            0x867b4bb4c42e5661 * rand(parent, UInt64))
end

copy(rng::Xoshiro) = Xoshiro(rng.s0, rng.s1, rng.s2, rng.s3)

function copy!(dst::Xoshiro, src::Xoshiro)
    dst.s0, dst.s1, dst.s2, dst.s3 = src.s0, src.s1, src.s2, src.s3
    dst
end

function ==(a::Xoshiro, b::Xoshiro)
    a.s0 == b.s0 && a.s1 == b.s1 && a.s2 == b.s2 && a.s3 == b.s3
end

rng_native_52(::Xoshiro) = UInt64

function seed!(rng::Xoshiro, seed::NTuple{4,UInt64})
    s = Base.hash_64_64(seed[1])
    rng.s0 = s
    s += Base.hash_64_64(seed[2])
    rng.s1 = s
    s += Base.hash_64_64(seed[3])
    rng.s2 = s
    s += Base.hash_64_64(seed[4])
    rng.s3 = s
    rng
end

@inline function rand(rng::Xoshiro, ::SamplerType{UInt64})
    s0, s1, s2, s3 = rng.s0, rng.s1, rng.s2, rng.s3
    tmp = s0 + s3
    res = tmp << 23 | tmp >> 41
    t = s1 << 17
    s2 = xor(s2, s0)
    s3 = xor(s3, s1)
    s1 = xor(s1, s2)
    s0 = xor(s0, s3)
    s2 = xor(s2, t)
    s3 = s3 << 45 | s3 >> 19
    rng.s0, rng.s1, rng.s2, rng.s3 = s0, s1, s2, s3
    res
end

## Task local RNG

"""
    TaskLocalRNG

The `TaskLocalRNG` has state that is local to its task, not its thread.
It is seeded upon task creation, from the state of its parent task.
Therefore, task creation is an event that changes the parent's RNG state.

As an upside, the `TaskLocalRNG` is pretty fast, and permits reproducible
multithreaded simulations (barring race conditions), independent of scheduler
decisions. As long as the number of threads is not used to make decisions on
task creation, simulation results are also independent of the number of available
threads / CPUs. The random stream should not depend on hardware specifics, up to
endianness and possibly word size.

Using or seeding the RNG of any other task than the one returned by `current_task()`
is undefined behavior: it will work most of the time, and may sometimes fail silently.
"""
struct TaskLocalRNG <: AbstractRNG end
TaskLocalRNG(::Nothing) = TaskLocalRNG()
rng_native_52(::TaskLocalRNG) = UInt64

function seed!(rng::TaskLocalRNG, seed::NTuple{4,UInt64})
    # TODO: Consider a less ad-hoc construction
    # We can afford burning a handful of cycles here, and we don't want any
    # surprises with respect to bad seeds / bad interactions.
    t = current_task()
    s = Base.hash_64_64(seed[1])
    t.rngState0 = s
    s += Base.hash_64_64(seed[2])
    t.rngState1 = s
    s += Base.hash_64_64(seed[3])
    t.rngState2 = s
    s += Base.hash_64_64(seed[4])
    t.rngState3 = s
    rand(rng, UInt64)
    rand(rng, UInt64)
    rand(rng, UInt64)
    rand(rng, UInt64)
    rng
end

@inline function rand(::TaskLocalRNG, ::SamplerType{UInt64})
    task = current_task()
    s0, s1, s2, s3 = task.rngState0, task.rngState1, task.rngState2, task.rngState3
    tmp = s0 + s3
    res = tmp << 23 | tmp >> 41
    t = s1 << 17
    s2 = xor(s2, s0)
    s3 = xor(s3, s1)
    s1 = xor(s1, s2)
    s0 = xor(s0, s3)
    s2 = xor(s2, t)
    s3 = s3 << 45 | s3 >> 19
    task.rngState0, task.rngState1, task.rngState2, task.rngState3 = s0, s1, s2, s3
    res
end

# Shared implementation between Xoshiro and TaskLocalRNG -- seeding
function seed!(rng::Union{TaskLocalRNG, Xoshiro}, seed::UInt128)
    seed0 = seed % UInt64
    seed1 = (seed>>>64) % UInt64
    seed!(rng, (seed0, seed1, zero(UInt64), zero(UInt64)))
end
seed!(rng::Union{TaskLocalRNG, Xoshiro}, seed::Integer) = seed!(rng, UInt128(seed))

seed!(rng::Union{TaskLocalRNG, Xoshiro}) =
    seed!(rng, (rand(RandomDevice(), UInt64), rand(RandomDevice(), UInt64),
                rand(RandomDevice(), UInt64), rand(RandomDevice(), UInt64)))

function seed!(rng::Union{TaskLocalRNG, Xoshiro}, seed::AbstractVector{UInt64})
    if length(seed) > 4
        throw(ArgumentError("seed should have no more than 256 bits"))
    end
    seed0 = length(seed)>0 ? seed[1] : UInt64(0)
    seed1 = length(seed)>1 ? seed[2] : UInt64(0)
    seed2 = length(seed)>2 ? seed[3] : UInt64(0)
    seed3 = length(seed)>3 ? seed[4] : UInt64(0)
    seed!(rng, (seed0, seed1, seed2, seed3))
end

function seed!(rng::Union{TaskLocalRNG, Xoshiro}, seed::AbstractVector{UInt32})
    if iseven(length(seed))
        seed!(rng, reinterpret(UInt64, seed))
    else
        seed!(rng, UInt64[reinterpret(UInt64, @view(seed[begin:end-1])); seed[end] % UInt64])
    end
end

@inline function rand(rng::Union{TaskLocalRNG, Xoshiro}, ::SamplerType{UInt128})
    first = rand(rng, UInt64)
    second = rand(rng,UInt64)
    second + UInt128(first)<<64
end

@inline rand(rng::Union{TaskLocalRNG, Xoshiro}, ::SamplerType{Int128}) = rand(rng, UInt128) % Int128

@inline function rand(rng::Union{TaskLocalRNG, Xoshiro},
                      T::SamplerUnion(Bool, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64))
    S = T[]
    # use upper bits
    (rand(rng, UInt64) >>> (64 - 8*sizeof(S))) % S
end

function copy(rng::TaskLocalRNG)
    t = current_task()
    Xoshiro(t.rngState0, t.rngState1, t.rngState2, t.rngState3)
end

function copy!(dst::TaskLocalRNG, src::Xoshiro)
    t = current_task()
    t.rngState0, t.rngState1, t.rngState2, t.rngState3 = src.s0, src.s1, src.s2, src.s3
    dst
end

function copy!(dst::Xoshiro, src::TaskLocalRNG)
    t = current_task()
    dst.s0, dst.s1, dst.s2, dst.s3 = t.rngState0, t.rngState1, t.rngState2, t.rngState3
    dst
end

function ==(a::Xoshiro, b::TaskLocalRNG)
    t = current_task()
    a.s0 == t.rngState0 && a.s1 == t.rngState1 && a.s2 == t.rngState2 && a.s3 == t.rngState3
end

==(a::TaskLocalRNG, b::Xoshiro) = b == a

# for partial words, use upper bits from Xoshiro

rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{UInt52Raw{UInt64}}) = rand(r, UInt64) >>> 12
rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{UInt52{UInt64}})    = rand(r, UInt64) >>> 12
rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{UInt104{UInt128}})  = rand(r, UInt104Raw())

rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{CloseOpen01{Float16}}) =
    Float16(Float32(rand(r, UInt16) >>> 5) * Float32(0x1.0p-11))

rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{CloseOpen01{Float32}}) =
    Float32(rand(r, UInt32) >>> 8) * Float32(0x1.0p-24)

rand(r::Union{TaskLocalRNG, Xoshiro}, ::SamplerTrivial{CloseOpen01_64}) =
    Float64(rand(r, UInt64) >>> 11) * 0x1.0p-53

max_stars_repo_name: mateuszbaran/julia

max_issues_repo_name: mateuszbaran/julia

max_forks_repo_name: mateuszbaran/julia


text:
module PyPSAData

using AxisArrays
using NCDatasets
using DataFrames
using CSV
using Base.Iterators: flatten

using ..EnergyModels
const EM = EnergyModels

using ..EnergyModels:
    AbstractNcData, pushaxes!, resolve, @consense, Component, Device,
    DeviceFormulation, attributes, components, naming, astype, typenames

export PypsaNcData

abstract type AbstractAttrDesc end

struct AttrDesc{T} <: AbstractAttrDesc
    pypsaname::String
    indices::Union{Colon,UnitRange{Int64},StepRange{Int64,Int64},Vector{Int64}}
    default::Union{T,Missing}
end

struct SingleAttrDesc{T} <: AbstractAttrDesc
    value::T
end

struct MissingAttrDesc{T} <: AbstractAttrDesc
    default::T
end

struct ComponentDesc
    name::Symbol
    axis::Axis
    componenttype
    attributes::Dict{Symbol,AbstractAttrDesc}
    typeparams::Union{Nothing,Dict{Symbol,Any}}
end

struct PypsaNcData <: AbstractNcData
    dataset::Dataset
    components::Dict{Symbol}{ComponentDesc}
    axes::Dict{Symbol}{Axis}

    function PypsaNcData(dataset, components)
        data = new(dataset, components, Dict{Symbol}{Axis}())
        for cd in values(data.components)
            pushaxes!(data, cd)
        end
        data
    end
end

"""
    _dimnames(::PypsaNcData, ::AbstractAttrDesc)

returns tuples for the dimensions EXCEPT for the main component axis.
"""
_dimnames(attrdesc::AbstractAttrDesc; ds=nothing) = ()
_dimnames(attrdesc::AttrDesc; ds=nothing) = dimnames(ds[attrdesc.pypsaname])[2:end]

function EM.pushaxes!(data::PypsaNcData, cd::ComponentDesc)
    pushaxes!(data, cd.axis, component=cd.name, attribute="primary")

    dims = unique(flatten((_dimnames(attrdesc; ds=data.dataset) for attrdesc in values(cd.attributes))))
    for dim in dims
        pushaxes!(data, Axis{Symbol(dim)}(nomissing(data.dataset[dim][:])), component=cd.name, attribute=dim)
    end
end

function _getattr(ds, pypsaname, attr)
    fullname = string(pypsaname, "_", attr)
    haskey(ds, fullname) ? nomissing(ds[fullname][:]) : Array{Union{Float64,Missing}}(missing, length(ds[string(pypsaname, "_i")]))
end

splitgroup(df, group) = [group]
function splitgroup(df, group, attr, attrs...)
    gdf = groupby(DataFrame(attr=>df[group.idx,attr], :idx=>group.idx, copycols=false), attr);

    if length(gdf) == 1
        ng = attr == :carrier ? (idx=group.idx, name=(group.name..., first(parent(gdf)[!,attr]))) : group
        groups = splitgroup(df, ng, attrs...)
    else
        groups = Any[]
        for sdf in gdf
            ng = (idx=collect(sdf.idx),
                  name=(group.name..., first(sdf[!,attr])))
            append!(groups, splitgroup(df, ng, attrs...))
        end
    end

    groups
end

AttrDesc(::Type{Val{:single}}, attr, name, ds, attrname, attrname_t, idx, axis) =
    SingleAttrDesc(
        haskey(ds, attrname) ?
        astype(attr[:dtype], @consense(ds[attrname][:][idx], "$name must be a single value")) :
        attr[:default]
    )

function AttrDesc(::Type{Val{:static}}, attr, name, ds, attrname, attrname_t, idx, axis)
    haskey(ds, attrname) ?
        AttrDesc{typenames[attr[:dtype]]}(attrname, idx, attr[:default]) :
        MissingAttrDesc(attr[:default])
end

function AttrDesc(::Type{Val{:series}}, attr, name, ds, attrname, attrname_t, idx, axis)
    attrname_t_i = attrname_t * "_i"
    indices_t = haskey(ds, attrname_t_i) ? findall(in(Symbol.(ds[attrname_t_i][:])), axis) : []
    if length(indices_t) == 0
        AttrDesc(Val{:static}, attr, name, ds, attrname, attrname_t, idx, axis)
    else
        @assert(length(indices_t) == length(axis),
                "$name has static and time-dependent $(attr[:PyPSA])")
        AttrDesc(attrname_t, indices_t, attr[:default])
    end
end

AttrDesc(attr, name, ds, pypsaname, idx, axis) =
    AttrDesc(Val{attr[:dimensions]}, attr, name, ds,
             string(pypsaname, "_", attr[:PyPSA]),
             string(pypsaname, "_t_", attr[:PyPSA]),
             idx, axis)

attributedescriptions(args...; attributes=nothing) =
    Dict(attr[:attribute]=>AttrDesc(attr, args...) for attr = eachrow(attributes))

function typeparamdescriptions(name, ds, listname, indices; types=nothing)
    typefield = string(listname, "_type")
    if types === nothing || !haskey(ds, typefield) return nothing end
    typename = @consense(ds[typefield][:][indices], "$name may not have more than one type")
    typ = types[types.name .== typename, :]
    Dict(c => typ[1, c] for c = names(typ) if c != :name)
end

"""
    PypsaNcData

"""
function PypsaNcData(ds)
    # Determine defined components and split them into classes
    components = vcat(
        getcomponents(ds, :buses),
        getcomponents(ds, :generators, withname=false, exp=:p_nom_extendable, uc=:committable),
        getcomponents(ds, :storageunits, withname=false, pypsaname="storage_units", exp=:p_nom_extendable),
        getcomponents(ds, :stores, carrier=:bus, exp=:e_nom_extendable),
        getcomponents(ds, :loads, carrier=:bus),
        getcomponents(ds, :links, carrier=(:bus0, :bus1), exp=:p_nom_extendable),
        getcomponents(ds, :lines, carrier=:bus0, exp=:s_nom_extendable),
        getcomponents(ds, :transformers, carrier=:bus0, exp=:s_nom_extendable)
    )

    compdict = Dict(cd.name => cd for cd in components)
    PypsaNcData(ds, compdict)
end

function _maybe_as_range(x::AbstractArray{T,1}) where T<:Number
    if length(x) < 3 return Vector(x) end
    d = diff(x)
    s = d[1]
    if all(d .== s)
        s == 1 ? (x[1]:x[end]) : (x[1]:s:x[end])
    else
        Vector(x)
    end
end

_join(a...; delim="_") = join(a, delim)
getcarrier(ds, buses_carrier, pypsaname, busattr) =
    getindex.(buses_carrier, nomissing(ds[string(pypsaname, "_", busattr)]))
getcarrier(ds, buses_carrier, pypsaname, busattrs::Tuple) =
    _join.((getcarrier(ds, buses_carrier, pypsaname, busattr) for busattr in busattrs)...)
function getcarrier(ds, pypsaname, busattrs)
    buses_carrier = Dict(nomissing(ds["buses_i"][:]) .=> nomissing(ds["buses_carrier"][:]))
    getcarrier(ds, buses_carrier, pypsaname, busattrs)
end

_length(carrier) = 1
_length(carrier::Tuple) = length(carrier)

function getcomponents(ds, name; pypsaname=false, withname=true, carrier=true, exp=false, uc=false)
    if pypsaname === false
        pypsaname = string(name)
    end

    components = ComponentDesc[]

    if !haskey(ds.dim, pypsaname * "_i") || ds.dim[pypsaname * "_i"] == 0
        return components
    end

    initialname = withname ? (name,) : ()

    df = DataFrame()
    if carrier === true
        df[!, :carrier] = _getattr(ds, pypsaname, :carrier)
    elseif carrier !== false
        if "buses_carrier" in ds
            df[!, :carrier] = getcarrier(ds, pypsaname, carrier)
        else
            initialname = (initialname..., join(repeat(["AC"], _length(carrier)), "_"))
        end
    end

    if exp !== false
        df[!, exp] = recode(_getattr(ds, pypsaname, exp),
                            0=>:fix, false=>:fix, 1=>:ext, true=>:ext)
    end
    if uc !== false
        df[!, uc] = recode(_getattr(ds, pypsaname, uc),
                           0=>:lin, false=>:lin, 1=>:uc, true=>:uc)
    end

    typesfn = joinpath(@__DIR__, string("pypsa.", pypsaname, ".types.csv"))
    types = isfile(typesfn) ? CSV.read(typesfn) : nothing

    componenttypeunion = resolve(Component, name)
    componentattributes = copy(attributes(componenttypeunion))
    componentattributes = componentattributes[componentattributes.quantitytype .!= "Variable", :]
    componentattributes[!, :PyPSA] = string.(componentattributes.attribute)
    componentattributes[!, :dimensions] = recode(length.(componentattributes.dimensions),
                                                 0=>:single, 1=>:static, 2=>:series,
                                                 3:20=>missing)

    idx = nomissing(ds[pypsaname * "_i"][:])
    timeattrs = componentattributes[componentattributes.dimensions .== :series, :attribute]
    for timeattr in timeattrs
        timeidx = string(pypsaname, "_t_", timeattr, "_i")
        if haskey(ds, timeidx)
            label = Symbol(:t_, timeattr)
            df[!, label] = recode(in.(idx, Ref(nomissing(ds[timeidx][:]))),
                                  false=>missing, true=>label)
        end
    end

    groups =
        if isempty(df)
            [(idx=1:size(idx, 1),
              name=Symbol(join(skipmissing(initialname), "_")),
              attrs=Dict{Symbol}{Any}())]
        else
            [(idx=_maybe_as_range(group.idx),
              name=Symbol(join(skipmissing(group.name), "_")),
              attrs=Dict(pairs(df[first(group.idx),:])...))
            for group in splitgroup(df, (idx=1:size(df, 1), name=initialname), names(df)...)]
        end

    for g in groups
        axis = Axis{g.name}(Symbol.(nomissing(ds[pypsaname * "_i"][:][g.idx])))

        formulation = join(skipmissing([
            exp !== false ? ((coalesce(g.attrs[exp], :fix) == :ext) ? "linexp" : missing) : missing,
            uc !== false ?  ((coalesce(g.attrs[uc], :lin) == :uc)   ? "ucdisp" : "lindisp") : "lindisp"
        ]), "_")

        componenttype =
            if componenttypeunion <: Device && !isempty(formulation)
                componenttypeunion{resolve(DeviceFormulation, Symbol(formulation))}
            else
                componenttypeunion
            end

        attributes = attributedescriptions(g.name, ds, pypsaname, g.idx, axis; attributes=componentattributes)
        typeparams = typeparamdescriptions(g.name, ds, pypsaname, g.idx; types=types)

        push!(components, ComponentDesc(g.name, axis, componenttype, attributes, typeparams))
    end

    components
end

Base.get(::PypsaNcData, attrdesc::SingleAttrDesc, _) = attrdesc.value
Base.get(::PypsaNcData, attrdesc::MissingAttrDesc, ax) = AxisArray(fill(attrdesc.default, length(ax)), ax)

function Base.get(data::PypsaNcData, attrdesc::AttrDesc{T}, ax) where T
    da = data.dataset[attrdesc.pypsaname]
    dims = Symbol.(dimnames(da)[2:end])

    # Benchmarking shows that getting the full array first and then
    # subsetting is faster even if we need to get only 1 or 2 lines
    # We would need to improve NCDatasets first!
    AxisArray(astype.(T, da[:][attrdesc.indices,ntuple(i->:,length(dims))...]),
              ax, (axis(data, n) for n = dims)...)
end

function Base.get(data::PypsaNcData, component::Component, param::Symbol)
    componentdesc = data.components[naming(component)]
    get(data, componentdesc.attributes[param], componentdesc.axis)
end

EM.gettypeparams(data::PypsaNcData, device::Device) =
    data.components[naming(device)].typeparams

EM.axis(data::PypsaNcData, e::Component) = data.axes[naming(e)]
EM.axis(data::PypsaNcData, name::Symbol) = data.axes[name]

EM.components(data::PypsaNcData) = EM._components(data)

end

using .PyPSAData

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
abstract type AbstractSampler
end

struct GibbsSampler <: AbstractSampler
   nsteps::Int
end

struct TemperatureDrivenSampler <: AbstractSampler
   nsteps::Int
   autocorcoeff::Float64
   betasd::Float64
end

function TemperatureDrivenSampler(;
      nsteps::Int = 5, autocorcoeff::Float64 = 0.9,
      betasd::Float64 = 0.9)

   TemperatureDrivenSampler(nsteps, autocorcoeff, betasd)
end

struct NoSampler <: AbstractSampler
end

function sample!(vmodel::M1, hmodel::M2, sampler::NoSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0
      ) where {M1 <: AbstractArray{Float64}, M2 <: AbstractArray{Float64}}
   # do nothing
end

function sample!(vmodel::M1, hmodel::M2, sampler::GibbsSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0
      ) where {M1 <: AbstractArray{Float64}, M2 <: AbstractArray{Float64}}

   for i = 1:sampler.nsteps
      samplevisible!(vmodel, rbm, hmodel, downfactor)
      samplehidden!(hmodel, rbm, vmodel, upfactor)
   end
   nothing
end

function sample!(particles::Particles, sampler::GibbsSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0)

   gibbssample!(particles, rbm, sampler.nsteps)
end

function sample!(vv::M1, hh::M2, sampler::TemperatureDrivenSampler, rbm::AbstractRBM,
      upfactor::Float64 = 1.0, downfactor::Float64 = 1.0
      ) where {M1 <: AbstractArray{Float64}, M2 <: AbstractArray{Float64}}
   # TODO respect upfactor and downfactor
   gibbssample_gamma!(vv, hh, rbm, sampler.nsteps;
         autocorcoeff = sampler.autocorcoeff, betasd = sampler.betasd)
end


text:
using ITensors

"""
    trg(T::ITensor; χmax::Int, nsteps::Int) -> κ, T

Perform the TRG algorithm on the partition function composed of the ITensor T.

The indices of T must obey come in pairs `(sₕ => sₕ')` and  `(sᵥ => sᵥ').

χmax is the maximum renormalized bond dimension.

nsteps are the number of renormalization steps performed.

The outputs are κ, the partition function per site, and the final renormalized
ITensor T.
"""
function trg(T::ITensor; χmax::Int, nsteps::Int, cutoff=0.0, svd_alg="divide_and_conquer")
  sₕ, sᵥ = filterinds(T; plev=0)
  @assert hassameinds((sₕ, sₕ', sᵥ, sᵥ'), T)

  # Keep track of the partition function per site
  κ = 1.0
  for n in 1:nsteps
    Fₕ, Fₕ′ = factorize(
      T,
      (sₕ', sᵥ');
      ortho="none",
      maxdim=χmax,
      cutoff=cutoff,
      tags=tags(sₕ),
      svd_alg=svd_alg,
    )

    s̃ₕ = commonind(Fₕ, Fₕ′)
    Fₕ′ *= δ(dag(s̃ₕ), s̃ₕ')

    Fᵥ, Fᵥ′ = factorize(
      T, (sₕ, sᵥ'); ortho="none", maxdim=χmax, cutoff=cutoff, tags=tags(sᵥ), svd_alg=svd_alg
    )

    s̃ᵥ = commonind(Fᵥ, Fᵥ′)
    Fᵥ′ *= δ(dag(s̃ᵥ), s̃ᵥ')

    T =
      (Fₕ * δ(dag(sₕ'), sₕ)) *
      (Fᵥ * δ(dag(sᵥ'), sᵥ)) *
      (Fₕ′ * δ(dag(sₕ), sₕ')) *
      (Fᵥ′ * δ(dag(sᵥ), sᵥ'))

    sₕ, sᵥ = s̃ₕ, s̃ᵥ

    trT = abs((T * δ(sₕ, sₕ') * δ(sᵥ, sᵥ'))[])
    T = T / trT
    κ *= trT^(1 / 2^n)
  end
  return κ, T
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module PkgTests

using Test
import Pkg
import Pkg.PkgError
using LibGit2
using Random: randstring

function capture_stdout(f::Function)
    let fname = tempname()
        try
            open(fname, "w") do fout
                redirect_stdout(fout) do
                    f()
                end
            end
            return read(fname, String)
        finally
            rm(fname, force=true)
        end
    end
end

function temp_pkg_dir(fn::Function, tmp_dir=joinpath(tempdir(), randstring()),
        remove_tmp_dir::Bool=true; initialize::Bool=true)

    # Used in tests below to set up and tear down a sandboxed package directory
    withenv("JULIA_PKGDIR" => tmp_dir) do
        @test !isdir(Pkg.dir())
        try
            if initialize
                Pkg.init()
                @test isdir(Pkg.dir())
                Pkg.resolve()
            else
                mkpath(Pkg.dir())
            end
            fn()
        finally
            remove_tmp_dir && rm(tmp_dir, recursive=true)
        end
    end
end

function write_build(pkg, content)
    build_filename = Pkg.dir(pkg, "deps", "build.jl")
    mkpath(dirname(build_filename))
    write(build_filename, content)
end

# Test basic operations: adding or removing a package, status, free
# Also test for the existence of REQUIRE and META_BRANCH
temp_pkg_dir() do
    @test isfile(joinpath(Pkg.dir(),"REQUIRE"))
    @test isfile(joinpath(Pkg.dir(),"META_BRANCH"))
    @test isempty(Pkg.installed())
    @test sprint(Pkg.status) == "No packages installed\n"
    @test !isempty(Pkg.available())

    # 18325
    cd(Pkg.dir()) do
        avail = Pkg.Read.available()
        avail_copy = Pkg.Query.availcopy(avail);
        delete!(avail_copy["Example"][v"0.0.1"].requires, "julia")
        @test haskey(avail["Example"][v"0.0.1"].requires, "julia")
    end

    @test_throws PkgError Pkg.installed("MyFakePackage")
    @test_throws PkgError Pkg.status("MyFakePackage")
    @test Pkg.installed("Example") === nothing

    # Check that setprotocol! works.
    begin
        try
            Pkg.setprotocol!("notarealprotocol")
            Pkg.add("Example")
            error("unexpected")
        catch ex
            if isa(ex, CompositeException)
                ex = ex.exceptions[1]

                if isa(ex, CapturedException)
                    ex = ex.ex
                end
            end
            @test isa(ex,Pkg.PkgError)
            @test contains(ex.msg, "Cannot clone Example from notarealprotocol://github.com/JuliaLang/Example.jl.git")
        end
    end

    Pkg.setprotocol!("")
    @test Pkg.Cache.rewrite_url_to === nothing
    Pkg.setprotocol!("https")
    Pkg.add("Example")
    @test [keys(Pkg.installed())...] == ["Example"]
    iob = IOBuffer()
    Pkg.checkout("Example")
    Pkg.status("Example", iob)
    str = chomp(String(take!(iob)))
    @test startswith(str, " - Example")
    @test endswith(str, "master")
    Pkg.free("Example")
    Pkg.status("Example", iob)
    str = chomp(String(take!(iob)))
    @test endswith(str, string(Pkg.installed("Example")))
    Pkg.checkout("Example")
    Pkg.free(("Example",))
    Pkg.status("Example", iob)
    str = chomp(String(take!(iob)))
    @test endswith(str, string(Pkg.installed("Example")))
    Pkg.rm("Example")
    @test isempty(Pkg.installed())
    @test !isempty(Pkg.available("Example"))
    @test !in("Example", keys(Pkg.installed()))
    Pkg.rm("Example")
    @test isempty(Pkg.installed())
    @test !isempty(Pkg.available("Example"))
    @test !in("Example", keys(Pkg.installed()))
    Pkg.clone("https://github.com/JuliaLang/Example.jl.git")
    @test [keys(Pkg.installed())...] == ["Example"]
    Pkg.status("Example", iob)
    str = chomp(String(take!(iob)))
    @test startswith(str, " - Example")
    @test endswith(str, "master")
    Pkg.free("Example")
    Pkg.status("Example", iob)
    str = chomp(String(take!(iob)))
    @test endswith(str, string(Pkg.installed("Example")))
    Pkg.checkout("Example")
    Pkg.free(("Example",))
    Pkg.status("Example", iob)
    str = chomp(String(take!(iob)))
    @test endswith(str, string(Pkg.installed("Example")))

    # 17364 - a, Pkg.checkout with specific local branch
    let branch_name = "test-branch-1",
        branch_commit = "ba3888212e30a7974ac6803a89e64c7098f4865e"

        # create a branch in Example package
        LibGit2.with(LibGit2.GitRepo, Pkg.dir("Example")) do repo
            LibGit2.branch!(repo, branch_name, branch_commit, set_head=false)
        end

        Pkg.clone(Pkg.dir("Example"), "Example2")
        Pkg.clone(Pkg.dir("Example"), "Example3")
        open(Pkg.dir("Example3", "README.md"), "w") do f
            println(f, "overwritten")
        end
        LibGit2.with(LibGit2.GitRepo, Pkg.dir("Example3")) do repo
            LibGit2.add!(repo, "README.md")
            test_sig = LibGit2.Signature("TEST", "TEST@TEST.COM", round(time(), 0), 0)
            LibGit2.commit(repo, "testmsg"; author=test_sig, committer=test_sig)
        end

        Pkg.checkout("Example2", branch_name)
        Pkg.checkout("Example3", branch_name)

        LibGit2.with(LibGit2.GitRepo, Pkg.dir("Example2")) do repo
            @test LibGit2.head_oid(repo) == LibGit2.GitHash(branch_commit)
        end
        LibGit2.with(LibGit2.GitRepo, Pkg.dir("Example3")) do repo
            @test LibGit2.head_oid(repo) == LibGit2.GitHash(branch_commit)
        end
    end

    # 17364 - b, remote off-tree branch
    let branch_name = "test-branch-2",
        branch_commit = "ba3888212e30a7974ac6803a89e64c7098f4865e"

        # create a branch in Example package
        LibGit2.with(LibGit2.GitRepo, Pkg.dir("Example")) do repo
            LibGit2.branch!(repo, branch_name, branch_commit, set_head=true)
        end

        # Make changes to local branch
        open(Pkg.dir("Example", "README.md"), "w") do f
            println(f, "overwritten")
        end

        test_commit = LibGit2.with(LibGit2.GitRepo, Pkg.dir("Example")) do repo
            LibGit2.add!(repo, "README.md")
            test_sig = LibGit2.Signature("TEST", "TEST@TEST.COM", round(time(), 0), 0)
            LibGit2.commit(repo, "testmsg"; author=test_sig, committer=test_sig)
        end
        Pkg.checkout("Example")

        Pkg.clone(Pkg.dir("Example"), "Example4")
        Pkg.checkout("Example4", branch_name)

        LibGit2.with(LibGit2.GitRepo, Pkg.dir("Example4")) do repo
            @test LibGit2.head_oid(repo) == test_commit
        end
    end

    # adding a package with unsatisfiable julia version requirements (REPL.jl) errors
    try
        Pkg.add("REPL")
        error("unexpected")
    catch err
        @test isa(err.exceptions[1].ex, PkgError)
        @test err.exceptions[1].ex.msg == "REPL can't be installed because " *
            "it has no versions that support $VERSION of julia. You may " *
            "need to update METADATA by running `Pkg.update()`"
    end

    # trying to add, check availability, or pin a nonexistent package errors
    try
        Pkg.add("NonexistentPackage")
        error("unexpected")
    catch err
        @test isa(err.exceptions[1].ex, PkgError)
        @test err.exceptions[1].ex.msg == "unknown package NonexistentPackage"
    end
    try
        Pkg.available("NonexistentPackage")
        error("unexpected")
    catch err
        @test isa(err, PkgError)
        @test err.msg == "NonexistentPackage is not a package (not registered or installed)"
    end
    try
        Pkg.pin("NonexistentPackage", v"1.0.0")
        error("unexpected")
    catch err
        @test isa(err, PkgError)
        @test err.msg == "NonexistentPackage is not a git repo"
    end

    # trying to pin a git repo under Pkg.dir that is not an installed package errors
    try
        Pkg.pin("METADATA", v"1.0.0")
        error("unexpected")
    catch err
        @test isa(err, PkgError)
        @test err.msg == "METADATA cannot be pinned – not an installed package"
    end

    # trying to pin an installed, registered package to an unregistered version errors
    try
        Pkg.pin("Example", v"2147483647.0.0")
        error("unexpected")
    catch err
        @test isa(err, PkgError)
        @test err.msg == "Example – 2147483647.0.0 is not a registered version"
    end

    # PR #13572, handling of versions with untagged detached heads
    LibGit2.with(LibGit2.GitRepo, Pkg.dir("Example")) do repo
        LibGit2.checkout!(repo, "72f09c7d0099793378c645929a9961155faae6d2")
    end
    @test Pkg.installed()["Example"] > v"0.0.0"

    # issue #13583
    begin
        try
            Pkg.test("IDoNotExist")
            error("unexpected")
        catch ex
            @test isa(ex,Pkg.Entry.PkgTestError)
            @test ex.msg == "IDoNotExist is not an installed package"
        end

        try
            Pkg.test("IDoNotExist1", "IDoNotExist2")
            error("unexpected")
        catch ex
            @test isa(ex,Pkg.Entry.PkgTestError)
            @test ex.msg == "IDoNotExist1 and IDoNotExist2 are not installed packages"
        end
    end

    # Various pin/free/re-pin/change-pin patterns (issue #17176)
    @test "" == capture_stdout() do
        nothingtodomsg = (:info,"No packages to install, update or remove")

        @test_logs((:info,"Freeing Example"),
                   nothingtodomsg,
                   Pkg.free("Example"))

        @test_logs (:info,r"^Creating Example branch pinned\.[0-9a-f]{8}\.tmp$") Pkg.pin("Example")
        vers = Pkg.installed("Example")
        branch = LibGit2.with(LibGit2.GitRepo, Pkg.dir("Example")) do repo
            LibGit2.branch(repo)
        end

        @test_logs((:info,"Freeing Example"),
                   nothingtodomsg,
                   Pkg.free("Example"))

        @test_logs((:info,"Creating Example branch pinned.b1990792.tmp"),
                   nothingtodomsg, Pkg.pin("Example", v"0.4.0"))
        @test Pkg.installed("Example") == v"0.4.0"

        @test_logs (:info,r"^Package Example is already pinned to the selected commit$") Pkg.pin("Example", v"0.4.0")
        @test Pkg.installed("Example") == v"0.4.0"

        @test_logs (:info,r"^Package Example is already pinned$") Pkg.pin("Example")
        @test Pkg.installed("Example") == v"0.4.0"

        @test_logs (:info,"Package Example: skipping update (pinned)...") match_mode=:any Pkg.update()
        @test Pkg.installed("Example") == v"0.4.0"

        @test_logs((:info,"Creating Example branch pinned.d1ef7b00.tmp"),
                   nothingtodomsg, Pkg.pin("Example", v"0.3.1"))
        @test Pkg.installed("Example") == v"0.3.1"

        @test_logs((:info,"Package Example: checking out existing branch pinned.b1990792.tmp"),
                   nothingtodomsg, Pkg.pin("Example", v"0.4.0"))
        @test Pkg.installed("Example") == v"0.4.0"

        @test_logs((:info,"Freeing Example"),
                   nothingtodomsg, Pkg.free("Example"))
        @test Pkg.installed("Example") == vers

        @test_logs (:info,Regex("^Package Example: checking out existing branch $branch\$")) Pkg.pin("Example")
        @test Pkg.installed("Example") == vers

        @test_logs((:info,"Freeing Example"),
                   nothingtodomsg, Pkg.free("Example"))
        @test Pkg.installed("Example") == vers
    end

    begin
        # bug identified in #16850, Base.url \ vs / for non-Base methods
        include(Pkg.dir("Example","src","Example.jl"))
        meth = first(methods(Example.domath))
        fname = string(meth.file)
        @test ('\\' in fname) == Sys.iswindows()
        @test startswith(Base.url(meth), "https://github.com/JuliaLang/Example.jl/tree")
    end

    # add a directory that is not a git repository
    begin
        mkdir(joinpath(Pkg.dir(), "NOTGIT"))
        Pkg.installed("NOTGIT") == typemin(VersionNumber)
        Pkg.installed()["NOTGIT"] == typemin(VersionNumber)
    end

    begin
        # don't bork when a Pkg repo is bare, issue #13804
        pth = joinpath(Pkg.dir(), "BAREGIT")
        mkdir(pth)
        # create a bare repo (isbare = true)
        repo = LibGit2.init(pth, true)
        @test repo.ptr != C_NULL
        close(repo)
        Pkg.update()
    end

    #test PkgDev redirects
    begin
        try
            Pkg.register("IDoNotExist")
            error("unexpected")
        catch ex
            @test ex.msg == "Pkg.register(pkg,[url]) has been moved to the package PkgDev.jl.\nRun Pkg.add(\"PkgDev\") to install PkgDev on Julia v0.5-"
        end

        try
            Pkg.tag("IDoNotExist")
            error("unexpected")
        catch ex
            @test ex.msg == "Pkg.tag(pkg, [ver, [commit]]) has been moved to the package PkgDev.jl.\nRun Pkg.add(\"PkgDev\") to install PkgDev on Julia v0.5-"
        end

        try
            Pkg.generate("IDoNotExist","MIT")
            error("unexpected")
        catch ex
            @test ex.msg == "Pkg.generate(pkg, license) has been moved to the package PkgDev.jl.\nRun Pkg.add(\"PkgDev\") to install PkgDev on Julia v0.5-"
        end

        try
            Pkg.publish()
            error("unexpected")
        catch ex
            @test ex.msg == "Pkg.publish() has been moved to the package PkgDev.jl.\nRun Pkg.add(\"PkgDev\") to install PkgDev on Julia v0.5-"
        end

        try
            Pkg.license()
            error("unexpected")
        catch ex
            @test ex.msg == "Pkg.license([lic]) has been moved to the package PkgDev.jl.\nRun Pkg.add(\"PkgDev\") to install PkgDev on Julia v0.5-"
        end
        try
            Pkg.submit("IDoNotExist")
            error("unexpected")
        catch ex
            @test ex.msg == "Pkg.submit(pkg[, commit]) has been moved to the package PkgDev.jl.\nRun Pkg.add(\"PkgDev\") to install PkgDev on Julia v0.5-"
        end
        try
            Pkg.submit("IDoNotExist", "nonexistentcommit")
            error("unexpected")
        catch ex
            @test ex.msg == "Pkg.submit(pkg[, commit]) has been moved to the package PkgDev.jl.\nRun Pkg.add(\"PkgDev\") to install PkgDev on Julia v0.5-"
        end
    end

    # Test Pkg.Read.url works
    @test Pkg.Read.url("Example") == "git://github.com/JuliaLang/Example.jl.git"

    # issue #15789, build failure warning are printed correctly.
    # Also makes sure `Pkg.build()` works even for non-git repo
    begin
        pth = joinpath(Pkg.dir(), "BuildFail")
        mkdir(pth)
        depspath = joinpath(pth, "deps")
        mkdir(depspath)
        depsbuild = joinpath(depspath, "build.jl")
        touch(depsbuild)
        # Pkg.build works without the src directory now
        # but it's probably fine to require it.
        msg = read(`$(Base.julia_cmd()) --startup-file=no -e 'redirect_stderr(stdout); using Logging; global_logger(SimpleLogger(stdout)); import Pkg; Pkg.build("BuildFail")'`, String)
        @test contains(msg, "Building BuildFail")
        @test !contains(msg, "Build failed for BuildFail")
        open(depsbuild, "w") do fd
            println(fd, "error(\"Throw build error\")")
        end
        msg = read(`$(Base.julia_cmd()) --startup-file=no -e 'redirect_stderr(stdout); using Logging; global_logger(SimpleLogger(stdout)); import Pkg; Pkg.build("BuildFail")'`, String)
        @test contains(msg, "Building BuildFail")
        @test contains(msg, "Build failed for BuildFail")
        @test contains(msg, "Pkg.build(\"BuildFail\")")
        @test contains(msg, "Throw build error")
    end

    # issue #15948
    let package = "Example"
        Pkg.rm(package)  # Remove package if installed
        @test Pkg.installed(package) === nothing  # Registered with METADATA but not installed
        msg = read(ignorestatus(`$(Base.julia_cmd()) --startup-file=no -e "redirect_stderr(stdout); using Logging; global_logger(SimpleLogger(stdout)); import Pkg; Pkg.build(\"$package\")"`), String)
        @test contains(msg, "$package is not an installed package")
        @test !contains(msg, "signal (15)")
    end

    # issue #20695
    Pkg.cd() do
        @test Pkg.Entry.url_and_pkg("Example") == ("git://github.com/JuliaLang/Example.jl.git", "Example")
        for url = [
            "https://github.com/Org/Nonsense",
            "git@github.com:Org/Nonsense",
            "file:///home/user/Nonsense",
            "/home/user/Nonsense",
        ]
            @test Pkg.Entry.url_and_pkg(url) == (url, "Nonsense")
            @test Pkg.Entry.url_and_pkg("$url.jl") == ("$url.jl", "Nonsense")
            @test Pkg.Entry.url_and_pkg("$url.git") == ("$url.git", "Nonsense")
            @test Pkg.Entry.url_and_pkg("$url.jl.git") == ("$url.jl.git", "Nonsense")
        end
        pkg = randstring(20)
        @test Pkg.Entry.url_and_pkg(pkg) == (pkg, pkg)
    end

    # partial Pkg.update
    @test "" == capture_stdout() do
        nothingtodomsg = (:info,"No packages to install, update or remove")

        @test_logs (:info,r"Installing Example v") match_mode=:any begin
            Pkg.rm("Example")
            Pkg.add("Example")
        end

        @test_logs nothingtodomsg match_mode=:any Pkg.update("Example")

        @test_logs (:info,"Installing Example v0.4.0") match_mode=:any begin
            Pkg.rm("Example")
            Pkg.add("Example", v"0", v"0.4.1-") # force version to be < 0.4.1
        end

        @test_logs((:info,r"""Package Example was set to version 0\.4\.0, but a higher version \d+\.\d+\.\d+\S* exists.
                              The update is prevented by explicit requirements constraints. Edit your REQUIRE file to change this."""),
                   nothingtodomsg, match_mode=:any, Pkg.update("Example"))

        @test_logs (:info,r"Installing Example") match_mode=:any begin
            Pkg.rm("Example")
            Pkg.add("Example")
            Pkg.pin("Example", v"0.4.0")
        end

        @test_logs((:info,"Package Example: skipping update (pinned)..."),
                   (:info,r"""Package Example was set to version 0\.4\.0, but a higher version \d+\.\d+\.\d+\S* exists.
                              The package is fixed. You can try using `Pkg.free\("Example"\)` to update it."""),
                   nothingtodomsg, match_mode=:any, Pkg.update("Example"))

        metadata_dir = Pkg.dir("METADATA")
        old_commit = "313bfaafa301e82d40574a778720e893c559a7e2"

        # Force a METADATA rollback to an old version, so that we will install some
        # outdated versions of some packages and then update some of those
        # (note that the following Pkg.update calls will update METADATA to the
        # latest version even though they don't update all packages)
        LibGit2.with(LibGit2.GitRepo, metadata_dir) do repo
            LibGit2.reset!(repo, LibGit2.GitHash(old_commit), LibGit2.Consts.RESET_HARD)
        end

        # run these at an old metadata commit where it's guaranteed no
        # packages depend on Example.jl
        @test isempty(Pkg.dependents("Example"))
        @test isempty(Pkg.dependents("Example.jl"))

        logs,_ = Test.collect_test_logs() do
            Pkg.add("Iterators")
            Pkg.update("Iterators")
        end
        @test all(!contains(l.message,"updated but were already imported") for l in logs)

        # Do it again, because the above Iterators test will update things prematurely
        LibGit2.with(LibGit2.GitRepo, metadata_dir) do repo
            LibGit2.reset!(repo, LibGit2.GitHash(old_commit), LibGit2.Consts.RESET_HARD)
        end

        @test_logs((:info,"Installing Colors v0.6.4"),
                   (:info,"Installing ColorTypes v0.2.2"),
                   (:info,"Installing FixedPointNumbers v0.1.3"),
                   (:info,"Installing Compat v0.7.18"),
                   (:info,"Installing Reexport v0.0.3"), match_mode=:any, Pkg.add("Colors"))

        logs,_ = Test.collect_test_logs() do
            Pkg.update("ColorTypes")
        end
        @test any(contains(l, (:info,r"Upgrading ColorTypes: v0\.2\.2 => v\d+\.\d+\.\d+")) for l in logs)
        @test any(contains(l, (:info,r"Upgrading Compat: v0\.7\.18 => v\d+\.\d+\.\d+")) for l in logs)
        @test !any(contains(l, (:info,r"Upgrading Colors")) for l in logs)

        @test Pkg.installed("Colors") == v"0.6.4"

        @test_logs nothingtodomsg match_mode=:any Pkg.update("FixedPointNumbers")
    end

    # issue #18239
    let package = "Example"
        Pkg.free(package)
        Pkg.rm(package)  # Remove package if installed

        metadata_dir = Pkg.dir("METADATA")
        old_commit = "83ff7116e51fc9cdbd7e67affbd344b9f5c9dbf2"

        # Reset METADATA to the second to last update of Example.jl
        LibGit2.with(LibGit2.GitRepo, metadata_dir) do repo
            LibGit2.reset!(repo, LibGit2.GitHash(old_commit), LibGit2.Consts.RESET_HARD)
        end

        Pkg.add(package)
        msg = read(ignorestatus(`$(Base.julia_cmd()) --startup-file=no -e
            "redirect_stderr(stdout); using Logging; global_logger(SimpleLogger(stdout)); using Example; import Pkg; Pkg.update(\"$package\")"`), String)
        @test contains(msg, Regex("- $package.*Restart Julia to use the updated versions","s"))
    end

    # Verify that the --startup-file flag is respected by Pkg.build / Pkg.test
    let package = "StartupFile"
        content = """
            @info "JULIA_RC_LOADED defined \$(isdefined(@__MODULE__, :JULIA_RC_LOADED))"
            @info "Main.JULIA_RC_LOADED defined \$(isdefined(Main, :JULIA_RC_LOADED))"
            """

        write_build(package, content)

        test_filename = Pkg.dir(package, "test", "runtests.jl")
        mkpath(dirname(test_filename))
        write(test_filename, content)

        # Make a ~/.julia/config/startup.jl
        home = Pkg.dir(".home")
        mkpath(joinpath(home, ".julia", "config"))
        write(joinpath(home, ".julia", "config", "startup.jl"),
            "const JULIA_RC_LOADED = true")

        withenv((Sys.iswindows() ? "USERPROFILE" : "HOME") => home) do
            code = "redirect_stderr(stdout); using Logging; global_logger(SimpleLogger(stdout)); import Pkg; Pkg.build(\"$package\")"
            msg = read(`$(Base.julia_cmd()) --startup-file=no -e $code`, String)
            @test contains(msg, "JULIA_RC_LOADED defined false")
            @test contains(msg, "Main.JULIA_RC_LOADED defined false")

            msg = read(`$(Base.julia_cmd()) --startup-file=yes -e $code`, String)
            @test contains(msg, "JULIA_RC_LOADED defined false")
            @test contains(msg, "Main.JULIA_RC_LOADED defined true")

            code = "redirect_stderr(stdout); using Logging; global_logger(SimpleLogger(stdout)); import Pkg; Pkg.test(\"$package\")"

            msg = read(`$(Base.julia_cmd()) --startup-file=no -e $code`, String)
            @test contains(msg, "JULIA_RC_LOADED defined false")
            @test contains(msg, "Main.JULIA_RC_LOADED defined false")

            # Note: Since both the startup-file and "runtests.jl" are run in the Main
            # module any global variables created in the startup file can be referenced.
            msg = read(`$(Base.julia_cmd()) --startup-file=yes -e $code`, String)
            @test contains(msg, "JULIA_RC_LOADED defined true")
            @test contains(msg, "Main.JULIA_RC_LOADED defined true")
        end
    end

    let package = "Output"
        stdout_file = Pkg.dir(package, "stdout.txt")
        stderr_file = Pkg.dir(package, "stderr.txt")
        content = """
            println(stdout, "stdout")
            println(stderr, "stderr")
            """
        write_build(package, content)

        code = "import Pkg; Pkg.build(\"$package\")"
        msg = run(pipeline(
            `$(Base.julia_cmd()) --startup-file=no -e $code`,
            stdout=stdout_file, stderr=stderr_file))
        @test last(readlines(stdout_file)) == "stdout"
        @test last(readlines(stderr_file)) == "stderr"
    end
end

@testset "Pkg functions with .jl extension" begin
    temp_pkg_dir() do
        @test Pkg.installed("Example.jl") === nothing
        Pkg.add("Example.jl")
        @test [keys(Pkg.installed())...] == ["Example"]
        iob = IOBuffer()
        Pkg.update("Example.jl")
        Pkg.checkout("Example.jl")
        Pkg.status("Example.jl", iob)
        str = chomp(String(take!(iob)))
        @test startswith(str, " - Example")
        @test endswith(str, "master")
        Pkg.free("Example.jl")
        Pkg.status("Example.jl", iob)
        str = chomp(String(take!(iob)))
        @test endswith(str, string(Pkg.installed("Example.jl")))
        Pkg.checkout("Example.jl")
        Pkg.free(("Example.jl",))
        Pkg.status("Example.jl", iob)
        str = chomp(String(take!(iob)))
        @test endswith(str, string(Pkg.installed("Example.jl")))
        Pkg.rm("Example.jl")
        @test isempty(Pkg.installed())
        @test !isempty(Pkg.available("Example.jl"))
        @test !in("Example", keys(Pkg.installed()))
        Pkg.rm("Example.jl")
        @test isempty(Pkg.installed())
        @test !isempty(Pkg.available("Example.jl"))
        @test !in("Example", keys(Pkg.installed()))
        Pkg.clone("https://github.com/JuliaLang/Example.jl.git")
        @test [keys(Pkg.installed())...] == ["Example"]
        Pkg.status("Example.jl", iob)
        str = chomp(String(take!(iob)))
        @test startswith(str, " - Example")
        @test endswith(str, "master")
        Pkg.free("Example.jl")
        Pkg.status("Example.jl", iob)
        str = chomp(String(take!(iob)))
        @test endswith(str, string(Pkg.installed("Example.jl")))
        Pkg.checkout("Example.jl")
        Pkg.free(("Example.jl",))
        Pkg.status("Example.jl", iob)
        str = chomp(String(take!(iob)))
        @test endswith(str, string(Pkg.installed("Example.jl")))
    end
end

let io = IOBuffer()
    Base.showerror(io, Pkg.Entry.PkgTestError("ppp"), backtrace())
    @test !contains(String(take!(io)), "backtrace()")
end

@testset "Relative path operations" begin
    cd(tempdir()) do
        temp_pkg_dir(randstring()) do
            Pkg.add("Example")
            @test [keys(Pkg.installed())...] == ["Example"]
        end
    end
end

temp_pkg_dir(initialize=false) do
    write_build("Normal", "")
    write_build("Error", "error(\"An error has occurred while building a package\")")
    write_build("Exit", "exit()")

    cd(Pkg.dir()) do
        errors = Dict()

        # Log forwarding TODO - port these to @test_logs
        empty!(errors)
        @test_warn ("Building Error",
                    "Building Normal") Pkg.Entry.build!(["Error", "Normal"], errors)

        empty!(errors)
        @test_warn ("Building Exit",
                    "Building Normal") Pkg.Entry.build!(["Exit", "Normal"], errors)

        empty!(errors)
        @test_warn ("Building Exit",
                    "Building Normal",
                    "Building Exit",
                    "Building Normal") Pkg.Entry.build!(["Exit", "Normal", "Exit", "Normal"], errors)
    end
end

# VersionSet tests
import Pkg.Types: VersionSet, VersionInterval

function chkint(a::VersionSet)
    ints = a.intervals
    for k = 1:length(ints)
        ints[k].lower < ints[k].upper || return false
        k < length(ints) && (ints[k].upper < ints[k+1].lower || return false)
    end
    return true
end

const empty_versionset = VersionSet(VersionInterval[])
@test isempty(empty_versionset)

# VersionSet intersections and unions
@test empty_versionset ∩ empty_versionset == empty_versionset
@test empty_versionset ∪ empty_versionset == empty_versionset
for t = 1:1_000
    a = VersionSet(sort!(map(v->VersionNumber(v...), [(rand(0:8),rand(0:3)) for i = 1:rand(0:10)]))...)
    b = VersionSet(sort!(map(v->VersionNumber(v...), [(rand(0:8),rand(0:3)) for i = 1:rand(0:10)]))...)
    @assert chkint(a)
    @assert chkint(b)
    u = a ∪ b
    @test chkint(u)
    i = a ∩ b
    @test chkint(i)
    for vM = 0:9, vm = 0:5
        v = VersionNumber(vM, vm)
        @test (v ∈ a || v ∈ b) ? (v ∈ u) : (v ∉ u)
        @test (v ∈ a && v ∈ b) ? (v ∈ i) : (v ∉ i)
    end
end

end # module

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
module FEM

import Base.push!
import Base.show
import Base.copy

using Base.LinAlg

using Compat
using FixedSizeArrays # 0.5 seconds
using Devectorize # 0.8 seconds
using FastAnonymous # 0.04 seconds
using Requires # 0.08 seconds
using Logging # 0.20 seconds

macro lintpragma(s) end

export NodeSet, gennodeset, ElementSet, DirichletBC, NodeLoad, Dof
export FENode2, FENode3
export GeoMesh, GeoTrig, GeoQTrig, GeoQuad, GeoTetra, GeoNode3, GeoNode2

# Elements, interpolators, materials, materialstatuses
export AbstractFElement, LinTrig, LinQuad, QuadTrig, GradTrig, activedofs
export AbstractInterpolator, LinTrigInterp, LinQuadInterp, QuadTrigInterp
export AbstractMaterial, LinearIsotropic, GradMekh
export AbstractElemStorage

export FESection, MaterialSection, ElementSection
export FEProblem
export Stress, Strain, InvFp, KAlpha, VonMises, Kappa, τ
export Solver, NRSolver, solve
export meshquad, gencook, read_mphtxt
export create_feproblem
export AbstractDataExporter, AbstractField

@Logging.configure(level=DEBUG, output = open("FEMlog.log", "w"))

t = @elapsed include("geomesh.jl")
@debug("Time to include geomesh.jl: $t s")

t = @elapsed include("core.jl")
@debug("Time to include core.jl: $t s")

t = @elapsed include("fields.jl")
@debug("Time to include fields.jl: $t s")

t = @elapsed include("materials/material.jl")
@debug("Time to include material.jl: $t s")

t = @elapsed include("interpolators/interpolator.jl")
@debug("Time to include interpolator.jl: $t s")

t = @elapsed include("elements/element.jl")
@debug("Time to include element.jl: $t s")

t = @elapsed include("sections.jl")
@debug("Time to include sections.jl: $t s")

t = @elapsed include("fe_problem.jl")
@debug("Time to include fe_problem.jl: $t s")

include("sparse_tools.jl")

include("export/export.jl")
include("solver.jl")
include("mesh/quad_mesher.jl")
include("mesh/read_mphtxt.jl")
include("mesh/mesh_conversions.jl")

end


text:
function eval_metric(::Val{:mse}, p::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where {T <: AbstractFloat}
    eval = zero(T)
    @inbounds for i in eachindex(y)
        eval += (p[1,i] - y[i])^2
    end
    eval /= length(y)
    return eval
end

function eval_metric(::Val{:rmse}, p::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where {T <: AbstractFloat}
    eval = zero(T)
    @inbounds for i in eachindex(y)
        eval += (p[1,i] - y[i])^2
    end
    eval = sqrt(eval / length(y))
    return eval
end

function eval_metric(::Val{:mae}, p::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where {T <: AbstractFloat}
    eval = zero(T)
    @inbounds for i in eachindex(y)
        eval += abs(p[1,i] - y[i])
    end
    eval /= length(y)
    return eval
end

function eval_metric(::Val{:logloss}, p::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where {T <: AbstractFloat}
    eval = zero(T)
    @inbounds for i in eachindex(y)
        pred = sigmoid(p[1,i])
        eval -= y[i] * log(pred) + (1 - y[i]) * log(1 - pred)
    end
    eval /= length(y)
    return eval
end

function eval_metric(::Val{:mlogloss}, p::AbstractMatrix{T}, y::AbstractVector{S}, α=0.0) where {T <: AbstractFloat,S <: Integer}
    eval = zero(T)
    p_prob = exp.(p) ./ sum(exp.(p), dims=1)
    @inbounds for i in eachindex(y)
        eval -= log(p_prob[y[i], i])
    end
    eval /= length(y)
    return eval
end

function eval_metric(::Val{:poisson}, p::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where {T <: AbstractFloat}
    eval = zero(T)
    @inbounds for i in eachindex(y)
        eval += exp(p[1,i]) * (1 - y[i]) + log(factorial(y[i]))
    end
    eval /= length(y)
    return eval
end

function eval_metric(::Val{:gaussian}, p::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where {T <: AbstractFloat}
    eval = zero(T)
    @inbounds for i in eachindex(y)
        eval += p[2,i] + (y[i] - p[1,i])^2 / (2 * exp(2 * p[2,i]))
    end
    eval /= length(y)
    return eval
end

function eval_metric(::Val{:quantile}, p::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where {T <: AbstractFloat}
    eval = zero(T)
    for i in eachindex(y)
        eval += α * max(y[i] - p[1,i], zero(T)) + (1 - α) * max(p[1,i] - y[i], zero(T))
    end
    eval /= length(y)
    return eval
end

function gini_raw(y::T, p::S) where {T,S}
    if length(y) < 2
        return 0.0
    end
    random = (1:length(p)) ./ length(p)
    l_sort = y[sortperm(p)]
    l_cum_w = cumsum(l_sort) ./ sum(y)
    gini = sum(l_cum_w .- random)
    return gini
end

function gini_norm(y::T, p::S) where {T,S}
    if length(y) < 2
        return 0.0
    end
    return gini_raw(y, p) / gini_raw(y, y)
end

function eval_metric(::Val{:gini}, p::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where {T <: AbstractFloat}
    return -gini_norm(y, view(p, 1, :))
end


text:
#!/usr/bin/env julia
using PkgServer, Sockets, Logging, LoggingExtras, Dates

# Accept optional environment-based arguments
# Eventually, do this via Pkg preferences
pkgserver = get(ENV, "JULIA_PKG_SERVER", "http://0.0.0.0:8000")
try
    m = match(r"(https?://)?(.+):(\d+)", pkgserver)
    global host = m.captures[2]
    global port = parse(Int, m.captures[3])
catch
    @warn("Invalid JULIA_PKG_SERVER setting, ignoring and using default of 0.0.0.0:8000!")
    global host = "0.0.0.0"
    global port = 8000
end

storage_root = get(ENV, "JULIA_PKG_SERVER_STORAGE_ROOT", "/tmp/pkgserver")
storage_servers = strip.(split(get(ENV, "JULIA_PKG_SERVER_STORAGE_SERVERS", "https://us-east.storage.juliahub.com,https://kr.storage.juliahub.com"), ","))
log_dir = get(ENV, "JULIA_PKG_SERVER_LOGS_DIR", joinpath(storage_root, "logs"))

mkpath(storage_root)
mkpath(log_dir)

# Set up logging
const date_format = "yyyy-mm-dd HH:MM:SS"
timestamp_logger(logger) = TransformerLogger(logger) do log
    merge(log, (; message = "[$(Dates.format(now(), date_format))] $(log.message)"))
end

global_logger(TeeLogger(
    timestamp_logger(
        MinLevelLogger(
            DatetimeRotatingFileLogger(
                log_dir,
                string(raw"yyyy-mm-dd-\p\k\g\s\e\r\v\e\r.\l\o\g"),
            ),
            Logging.Info,
        ),
    ),
    current_logger(),
))

PkgServer.start(;
    listen_addr=Sockets.InetAddr(host, port),
    storage_root=storage_root,
    storage_servers=storage_servers,
)


text:
# Linear Discriminant Analysis

#### Type to represent a linear discriminant functional

abstract type Discriminant{T} end

struct LinearDiscriminant{T<:Real} <: Discriminant{T}
    w::Vector{T}
    b::T
end

length(f::LinearDiscriminant) = length(f.w)

evaluate(f::LinearDiscriminant, x::AbstractVector) = dot(f.w, x) + f.b

function evaluate(f::LinearDiscriminant, X::AbstractMatrix)
    R = transpose(X) * f.w
    if f.b != 0
        broadcast!(+, R, R, f.b)
    end
    return R
end

predict(f::Discriminant, x::AbstractVector) = evaluate(f, x) > 0

predict(f::Discriminant, X::AbstractMatrix) = (Y = evaluate(f, X); Bool[y > 0 for y in Y])

#### function to solve linear discriminant

function ldacov(C::DenseMatrix{T},
                μp::DenseVector{T},
                μn::DenseVector{T}) where T<:Real

    w = cholesky(C) \ (μp - μn)
    ap = w ⋅ μp
    an = w ⋅ μn
    c = 2 / (ap - an)
    LinearDiscriminant(rmul!(w, c), 1 - c * ap)
end

ldacov(Cp::DenseMatrix{T},
       Cn::DenseMatrix{T},
       μp::DenseVector{T},
       μn::DenseVector{T}) where T<:Real = ldacov(Cp + Cn, μp, μn)

#### interface functions

function fit(::Type{LinearDiscriminant}, Xp::DenseMatrix{T}, Xn::DenseMatrix{T};
             covestimator::CovarianceEstimator = SimpleCovariance()) where T<:Real
    μp = vec(mean(Xp, dims=2))
    μn = vec(mean(Xn, dims=2))
    Zp = Xp .- μp
    Zn = Xn .- μn
    Cp = calcscattermat(covestimator, Zp)
    Cn = calcscattermat(covestimator, Zn)
    ldacov(Cp, Cn, μp, μn)
end

#### Multiclass LDA Stats

mutable struct MulticlassLDAStats{T<:Real, M<:AbstractMatrix{T}, N<:AbstractMatrix{T}}
    dim::Int              # sample dimensions
    nclasses::Int         # number of classes
    cweights::Vector{T}   # class weights
    tweight::T            # total sample weight
    mean::Vector{T}       # overall sample mean
    cmeans::Matrix{T}     # class-specific means
    Sw::M                 # within-class scatter matrix
    Sb::N                 # between-class scatter matrix
end

mean(S::MulticlassLDAStats) = S.mean
classweights(S::MulticlassLDAStats) = S.cweights
classmeans(S::MulticlassLDAStats) = S.cmeans

withclass_scatter(S::MulticlassLDAStats) = S.Sw
betweenclass_scatter(S::MulticlassLDAStats) = S.Sb

function MulticlassLDAStats(cweights::Vector{T},
                            mean::Vector{T},
                            cmeans::Matrix{T},
                            Sw::AbstractMatrix{T},
                            Sb::AbstractMatrix{T}) where T<:Real
    d, nc = size(cmeans)
    length(mean) == d || throw(DimensionMismatch("Incorrect length of mean"))
    length(cweights) == nc || throw(DimensionMismatch("Incorrect length of cweights"))
    tw = sum(cweights)
    size(Sw) == (d, d) || throw(DimensionMismatch("Incorrect size of Sw"))
    size(Sb) == (d, d) || throw(DimensionMismatch("Incorrect size of Sb"))
    MulticlassLDAStats(d, nc, cweights, tw, mean, cmeans, Sw, Sb)
end

function multiclass_lda_stats(nc::Int, X::AbstractMatrix{T}, y::AbstractVector{Int};
                              covestimator_within::CovarianceEstimator=SimpleCovariance(),
                              covestimator_between::CovarianceEstimator=SimpleCovariance()) where T<:Real
    # check sizes
    d = size(X, 1)
    n = size(X, 2)
    n ≥ nc || throw(ArgumentError("The number of samples is less than the number of classes"))
    length(y) == n || throw(DimensionMismatch("Inconsistent array sizes."))

    # compute class-specific weights and means
    cmeans, cweights, Z = center(X, y, nc)

    Sw = calcscattermat(covestimator_within, Z)

    # compute between-class scattering
    mean = cmeans * (cweights ./ T(n))
    U = rmul!(cmeans .- mean, Diagonal(sqrt.(cweights)))
    Sb = calcscattermat(covestimator_between, U)

    return MulticlassLDAStats(Vector{T}(cweights), mean, cmeans, Sw, Sb)
end

#### Multiclass LDA

mutable struct MulticlassLDA{T<:Real}
    proj::Matrix{T}
    pmeans::Matrix{T}
    stats::MulticlassLDAStats{T}
end

indim(M::MulticlassLDA) = size(M.proj, 1)
outdim(M::MulticlassLDA) = size(M.proj, 2)

projection(M::MulticlassLDA) = M.proj

mean(M::MulticlassLDA) = mean(M.stats)
classmeans(M::MulticlassLDA) = classmeans(M.stats)
classweights(M::MulticlassLDA) = classweights(M.stats)

withclass_scatter(M::MulticlassLDA) = withclass_scatter(M.stats)
betweenclass_scatter(M::MulticlassLDA) = betweenclass_scatter(M.stats)

transform(M::MulticlassLDA, x::AbstractVecOrMat{<:Real}) = M.proj'x

function fit(::Type{MulticlassLDA}, nc::Int, X::DenseMatrix{T}, y::AbstractVector{Int};
             method::Symbol=:gevd,
             outdim::Int=min(size(X,1), nc-1),
             regcoef::T=T(1.0e-6),
             covestimator_within::CovarianceEstimator=SimpleCovariance(),
             covestimator_between::CovarianceEstimator=SimpleCovariance()) where T<:Real

    multiclass_lda(multiclass_lda_stats(nc, X, y;
                                        covestimator_within=covestimator_within,
                                        covestimator_between=covestimator_between);
                   method=method,
                   regcoef=regcoef,
                   outdim=outdim)
end

function multiclass_lda(S::MulticlassLDAStats{T};
                        method::Symbol=:gevd,
                        outdim::Int=min(size(X,1), S.nclasses-1),
                        regcoef::T=T(1.0e-6)) where T<:Real

    P = mclda_solve(S.Sb, S.Sw, method, outdim, regcoef)
    MulticlassLDA(P, P'S.cmeans, S)
end

mclda_solve(Sb::AbstractMatrix{T}, Sw::AbstractMatrix{T}, method::Symbol, p::Int, regcoef::T) where T<:Real =
    mclda_solve!(copy(Sb), copy(Sw), method, p, regcoef)

function mclda_solve!(Sb::AbstractMatrix{T},
                      Sw::AbstractMatrix{T},
                      method::Symbol, p::Int, regcoef::T) where T<:Real

    p <= size(Sb, 1) || throw(ArgumentError("p cannot exceed sample dimension."))

    if method == :gevd
        regularize_symmat!(Sw, regcoef)
        E = eigen!(Symmetric(Sb), Symmetric(Sw))
        ord = sortperm(E.values; rev=true)
        P = E.vectors[:, ord[1:p]]

    elseif method == :whiten
        W = _lda_whitening!(Sw, regcoef)
        wSb = transpose(W) * (Sb * W)
        Eb = eigen!(Symmetric(wSb))
        ord = sortperm(Eb.values; rev=true)
        P = W * Eb.vectors[:, ord[1:p]]

    else
        throw(ArgumentError("Invalid method name $(method)"))
    end
    return P::Matrix{T}
end

function _lda_whitening!(C::AbstractMatrix{T}, regcoef::T) where T<:Real
    n = size(C,1)
    E = eigen!(Symmetric(C))
    v = E.values
    a = regcoef * maximum(v)
    for i = 1:n
        @inbounds v[i] = 1.0 / sqrt(v[i] + a)
    end
    return rmul!(E.vectors,  Diagonal(v))
end

#### SubspaceLDA

# When the dimensionality is much higher than the number of samples,
# it makes more sense to perform LDA on the space spanned by the
# within-group scatter.

struct SubspaceLDA{T<:Real}
    projw::Matrix{T}
    projLDA::Matrix{T}
    λ::Vector{T}
    cmeans::Matrix{T}
    cweights::Vector{Int}
end

indim(M::SubspaceLDA) = size(M.projw,1)
outdim(M::SubspaceLDA) = size(M.projLDA, 2)

projection(M::SubspaceLDA) = M.projw * M.projLDA

mean(M::SubspaceLDA) = vec(sum(M.cmeans * Diagonal(M.cweights / sum(M.cweights)), dims=2))
classmeans(M::SubspaceLDA) = M.cmeans
classweights(M::SubspaceLDA) = M.cweights

transform(M::SubspaceLDA, x) = M.projLDA' * (M.projw' * x)

fit(::Type{F}, X::AbstractMatrix{T}, nc::Int, label::AbstractVector{Int}) where {T<:Real, F<:SubspaceLDA} =
    fit(F, X, label, nc)

function fit(::Type{F}, X::AbstractMatrix{T},
             label::AbstractVector{Int},
             nc=maximum(label);
             normalize::Bool=false) where {T<:Real, F<:SubspaceLDA}
    d, n = size(X, 1), size(X, 2)
    n ≥ nc || throw(ArgumentError("The number of samples is less than the number of classes"))
    length(label) == n || throw(DimensionMismatch("Inconsistent array sizes."))
    # Compute centroids, class weights, and deviation from centroids
    # Note Sb = Hb*Hb', Sw = Hw*Hw'
    cmeans, cweights, Hw = center(X, label, nc)
    dmeans = cmeans .- (normalize ? mean(cmeans, dims=2) : cmeans * (cweights / T(n)))
    Hb = normalize ? dmeans : dmeans * Diagonal(convert(Vector{T}, sqrt.(cweights)))
    if normalize
        Hw /= T(sqrt(n))
    end
    # Project to the subspace spanned by the within-class scatter
    # (essentially, PCA before LDA)
    Uw, Σw, _ = svd(Hw, full=false)
    keep = Σw .> sqrt(eps(T)) * maximum(Σw)
    projw = Uw[:,keep]
    pHb = projw' * Hb
    pHw = projw' * Hw
    λ, G = lda_gsvd(pHb, pHw, cweights)
    SubspaceLDA(projw, G, λ, cmeans, cweights)
end

# Reference: Howland & Park (2006), "Generalizing discriminant analysis
# using the generalized singular value decomposition", IEEE
# Trans. Patt. Anal. & Mach. Int., 26: 995-1006.
function lda_gsvd(Hb::AbstractMatrix{T}, Hw::AbstractMatrix{T}, cweights::AbstractVector{Int}) where T<:Real
    nc = length(cweights)
    K = vcat(Hb', Hw')
    P, R, Q = svd(K, full=false)
    keep = R .> sqrt(eps(T))*maximum(R)
    R = R[keep]
    Pk = P[1:nc, keep]
    U, ΣA, W = svd(Pk)
    ncnz = sum(cweights .> 0)
    G = Q[:,keep]*(Diagonal(1 ./ R) * W[:,1:ncnz-1])
    # Normalize
    Gw = G' * Hw
    nrm = Gw * Gw'
    G = G ./ reshape(sqrt.(diag(nrm)), 1, ncnz-1)
    # Also get the eigenvalues
    Gw = G' * Hw
    Gb = G' * Hb
    λ = diag(Gb * Gb')./diag(Gw * Gw')
    λ, G
end

function center(X::AbstractMatrix{T}, label::AbstractVector{Int}, nc=maximum(label)) where T<:Real
    d, n = size(X,1), size(X,2)
    # Calculate the class weights and means
    cmeans = zeros(T, d, nc)
    cweights = zeros(Int, nc)
    for j = 1:n
        k = label[j]
        for i = 1:d
            cmeans[i,k] += X[i,j]
        end
        cweights[k] += 1
    end
    for j = 1:nc
        cw = cweights[j]
        cw == 0 && continue
        for i = 1:d
            cmeans[i,j] /= cw
        end
    end
    # Compute differences from the means
    dX = Matrix{T}(undef, d, n)
    for j = 1:n
        k = label[j]
        for i = 1:d
            dX[i,j] = X[i,j] - cmeans[i,k]
        end
    end
    cmeans, cweights, dX
end


text:
using Test

# run tests on Travis CI in parallel
const TRIXI_TEST = get(ENV, "TRIXI_TEST", "all")
const ON_APPVEYOR = lowercase(get(ENV, "APPVEYOR", "false")) == "true"

@time @testset "Trixi.jl tests" begin
  @time if TRIXI_TEST == "all" || TRIXI_TEST == "2D"
    include("test_examples_2d.jl")
  end

  @time if (TRIXI_TEST == "all" && !ON_APPVEYOR) || TRIXI_TEST == "3D"
    include("test_examples_3d.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc"
    include("test_manual.jl")
    include("test_elixirs.jl")
  end

  @time if (TRIXI_TEST == "all" && !ON_APPVEYOR) || TRIXI_TEST == "paper-self-gravitating-gas-dynamics"
    include("test_paper-self-gravitating-gas-dynamics.jl")
  end
end