text:
@testset "Trading metrics from transactions blotter"
  @testset "Transaction blotter transformations"
    blotter = TradingLogic.emptyblotter()

    # long than short with exit
    p1, p2, p3 = 108.0, 85.0, 74.0
    along, ashort = 50, 70
    tlong = Dates.DateTime(2015,1,1,3,8,17)
    tshort = Dates.DateTime(2015,1,1,5,2,4)
    texit = Dates.DateTime(2015,1,3,15,12,48)

    # enter transactions not in chronological order
    blotter[tshort] = (-(along + ashort), p2)
    blotter[tlong] = (along, p1)
    blotter[texit] = (ashort, p3)
    TradingLogic.printblotter(STDOUT, blotter)

    # sorted vectors
    @test TradingLogic.vtblotter(blotter) == [tlong, tshort, texit]
    va, vp = TradingLogic.vapblotter(blotter)
    @test va == [along, -(along + ashort), ashort]
    @test vp == roughly([p1, p2, p3])

    # print to file
    dtfmt = "yyyy-mm-dd HH:MM:SS" # not default
    fnm = "writetest.csv"
    TradingLogic.writeblotter(fnm, blotter, dtformat = dtfmt)
    rcsv = readcsv(fnm)
    run(`rm $fnm`)
    @test size(rcsv) == (4,3)
    @test rcsv[1,2] == "Amount"
    @test rcsv[2,1] == Dates.format(tlong, dtfmt)
    @test rcsv[4,1] == Dates.format(texit, dtfmt)
    ##@test int(rcsv[3,2]) == -(along + ashort)
    @test round(Int64, rcsv[3,2]) == -(along + ashort)
    @test float(rcsv[4,3]) == roughly(p3)
  end
  @testset "Long enter with partial fill"
    blotter = TradingLogic.emptyblotter()
    # fill prices and exit price
    p1, p2, pexit = 120.0, 135.0, 150.0
    # amounts
    a1, a2 = 20, 30

    blotter[DateTime(2015,1,1,0)] = (a1, p1)
    blotter[DateTime(2015,1,1,1)] = (a2, p2)
    blotter[DateTime(2015,1,2)] = (-(a1 + a2), pexit)
    pnlfin = a1*(p2 - p1) + (a1 + a2)*(pexit - p2)

    metr = [:DDown]
    vt, perfm = TradingLogic.tradeperf(blotter, metr)
    ntrans = length(blotter)
    @test length(vt) == ntrans

    @test length(perfm[:PnL]) == ntrans

    @test perfm[:Qty] == [a1, a2, -(a1 + a2)]
    @test perfm[:FillPrice] == roughly([p1, p2, pexit])

    @test perfm[:PnL] == roughly([0.0, a1*(p2 - p1), pnlfin])
    @test TradingLogic.tradepnlfinal(blotter) == roughly(pnlfin)
  end
  @testset "Short enter and exit"
    blotter = TradingLogic.emptyblotter()
    # short enter and exit price
    pent, pexit = 120.0, 115.2
    # short position held
    a = 100

    blotter[DateTime(2015,1,1,0)] = (-a, pent)
    blotter[DateTime(2015,1,1,1)] = (a, pexit)
    pnlfin = -a*(pexit - pent)

    metr = [:DDown]
    vt, perfm = TradingLogic.tradeperf(blotter, metr)
    ntrans = length(blotter)
    @test length(vt) == ntrans

    @test length(perfm[:PnL]) == ntrans

    @test perfm[:Qty] == [-a, a]
    @test perfm[:FillPrice] == roughly([pent, pexit])

    @test perfm[:PnL] == roughly([0.0, pnlfin])
    @test TradingLogic.tradepnlfinal(blotter) == roughly(pnlfin)
  end
  @testset "Long then short"
    blotter = TradingLogic.emptyblotter()
    # prices and amounts
    p1, p2, p3 = 120.0, 85.0, 70.0
    along, ashort = 50, 70

    losslong = along*(p2 - p1)
    pnlfin = losslong + ashort*(p2 - p3)

    blotter[DateTime(2015,1,1,0)] = (along, p1)
    blotter[DateTime(2015,1,1,1)] = (-(along + ashort), p2)

    # not yet final PnL with incomplete blotter
    @test abs(TradingLogic.tradepnlfinal(blotter) - pnlfin) > 10.0 == true
    # test tradepnlfinal method with pnow: final PnL
    #  doing this before adding the final (exit) transaction to blotter
    @test TradingLogic.tradepnlfinal(blotter, p3) == roughly(pnlfin)

    # final transaction (e.g. stop backtest point)
    blotter[DateTime(2015,1,2)] = (0, p3)

    metr = [:DDown]
    vt, perfm = TradingLogic.tradeperf(blotter, metr)
    #println(vt, perfm)
    ntrans = length(blotter)
    @test length(vt) == ntrans
    @test length(perfm[:PnL]) == ntrans

    @test perfm[:Qty] == [along, -(along + ashort), 0]
    @test perfm[:FillPrice] == roughly([p1, p2, p3])

    @test perfm[:PnL] == roughly([0.0, losslong, pnlfin])
    @test TradingLogic.tradepnlfinal(blotter) == roughly(pnlfin)

    # return-based drawdown
    ddmax = abs(losslong)
    @test perfm[:DDown] == roughly([0.0, -ddmax, -abs(pnlfin)])
  end
end

@testset "Running performance metrics"
  @testset "Current maximum PnL and drawdown over trading history"
    pnl_dd_max = TradingLogic.tradeperffold((10.0, 0.0), (true, 15.0))
    @test pnl_dd_max[1] == roughly(15.0)
    @test pnl_dd_max[2] == roughly(0.0)

    pnl_dd_max = TradingLogic.tradeperffold((10.0, 20.0), (true, 15.0))
    @test pnl_dd_max[1] == roughly(15.0)
    @test pnl_dd_max[2] == roughly(20.0)

    pnl_dd_max = TradingLogic.tradeperffold((10.0, 0.0), (true, 5.0))
    @test pnl_dd_max[1] == roughly(10.0)
    @test pnl_dd_max[2] == roughly(5.0)

    pnl_dd_max = TradingLogic.tradeperffold((10.0, 50.0), (true, 5.0))
    @test pnl_dd_max[1] == roughly(10.0)
    @test pnl_dd_max[2] == roughly(50.0)
  end
end

@testset "Trades and corresponding metrics from transactions blotter"
  blotter = TradingLogic.emptyblotter()
  dthour(h::Int64) = Dates.DateTime(2015,6,1,h)
  blotter[dthour(1)] = (25, 85.0)
  blotter[dthour(2)] = (25, 95.0)
  blotter[dthour(3)] = (-60, 100.0)
  blotter[dthour(4)] = (20, 80.0)
  blotter[dthour(5)] = (-10, 75.0)
  blotter[dthour(6)] = (10, 85.0)
  blotter[dthour(7)] = (-10, 70.0)
  blotter[dthour(8)] = (-50, 65.0)
  blotter[dthour(11)] = (50, 25.0)
  # open trade
  blotter[dthour(12)] = (20, 35.0)
  #TradingLogic.printblotter(STDOUT, blotter)
  @test TradingLogic.tradepnlfinal(blotter) == roughly(2500.0)

  @testset "Trades PnL vector with selected metrics from blotter"
    vtrpnl, ntrpos, avwin, ntrneg, avloss = TradingLogic.vtradespnl(blotter)
    @test length(vtrpnl) == 5
    @test sum(vtrpnl) == roughly(2500.0)
    @test ntrpos == 3
    @test avwin == roughly(900.0)
    @test ntrneg == 2
    @test avloss == roughly(100.0)
  end

  @testset "Trades-based performance metrics: PROM"
    marg = 5e3

    # no transactions
    @test TradingLogic.perf_prom(
      TradingLogic.emptyblotter()) == roughly(0.0)

    # from blotter, using all completed trades
    prom0 = ((3.0 - sqrt(3.0))*900.0 - (2.0 + sqrt(2.0))*100.0)/marg
    @test TradingLogic.perf_prom(blotter, marg=marg) == roughly(prom0)

    # removing specific number of best wins
    prom1 = ((2.0 - sqrt(2.0))*350.0 - (2.0 + sqrt(2.0))*100.0)/marg
    @test TradingLogic.perf_prom(blotter, marg=marg,
                                 nbest_remove = 1) == roughly(prom1)

    # from trades PnL vector
    vtrpnl, ntrpos, avwin, ntrneg, avloss = TradingLogic.vtradespnl(blotter)
    @test TradingLogic.perf_prom(vtrpnl, marg=marg,
                                 nbest_remove = 1) == roughly(prom1)
    prom2 = -(2.0 + sqrt(2.0))*100.0/marg
    @test TradingLogic.perf_prom(vtrpnl, marg=marg,
                                 nbest_remove = 2) == roughly(prom2)
    # if asking to remove more best wins than the total number of wins
    prom = TradingLogic.perf_prom(vtrpnl, marg=marg, nbest_remove = 5)
    # -Inf
    @test prom < 0.0 == true
    @test isfinite(prom) == false
    prom = TradingLogic.perf_prom(vtrpnl, marg=marg, nbest_remove = 6)
    # -Inf
    @test prom < 0.0 == true
    @test isfinite(prom) == false
  end

  @testset "Trades-based performance metrics: PROR"
    pror0 = (3.0 - sqrt(3.0))*900.0 / ((2.0 + sqrt(2.0))*100.0)
    @test TradingLogic.perf_prom(blotter, pror=true) == roughly(pror0)
    pror1 = (2.0 - sqrt(2.0))*350.0 / ((2.0 + sqrt(2.0))*100.0)
    @test TradingLogic.perf_prom(blotter, pror=true,
                                 nbest_remove = 1) == roughly(pror1)

    # PROR without profits
    @test TradingLogic.perf_prom(TradingLogic.emptyblotter(),
                                 pror=true) == roughly(0.0)
    @test TradingLogic.perf_prom(Array(Float64,0), pror=true) == roughly(0.0)
    @test TradingLogic.perf_prom([-2.0], pror=true) == roughly(0.0)
    pror = TradingLogic.perf_prom([-2.0], pror=true, nbest_remove = 1)
    # -Inf
    @test pror < 0.0 == true
    @test isfinite(pror) == false

    # PROR without profits or losses: 1 - sqrt(1)
    @test TradingLogic.perf_prom([5.0, -2.0], pror=true) == roughly(0.0)
    @test TradingLogic.perf_prom([5.0, -2.0], pror=true,
                                 nbest_remove = 1) == roughly(0.0)
    # PROR without losses
    pror = TradingLogic.perf_prom([5.0, 2.0], pror=true, nbest_remove = 0)
    # +Inf
    @test pror < 0.0 == false
    @test isfinite(pror) == false
    @test TradingLogic.perf_prom([5.0, 2.0], pror=true,
                                 nbest_remove = 1) == roughly(0.0)
    pror = TradingLogic.perf_prom([5.0, 2.0], pror=true, nbest_remove = 2)
    # -Inf
    @test pror < 0.0 == true
    @test isfinite(pror) == false
  end
end


text:
# __precompile__()

module SubMatrixSelectionSVD

using LinearAlgebra
using Statistics
using Distributed

export
    smssvd,
    projectionscore,
    projectionscorefiltered

include("projectionscore.jl")
include("smssvdimpl.jl")
include("precompile.jl")
# _precompile_()

end

ext: jl


text:
using DoubleBLAS
using LinearAlgebra
using DoubleFloats: Double32, Double64
using Test, Random

let nt = Threads.nthreads()
    nts = Sys.CPU_THREADS
    @info "running tests with $nt threads of $nts"
end

Random.seed!(1234)

include("level1.jl")
include("givens.jl")

function gemmcheck(T,m,n,k,tol)
    A = rand(T,m,k)
    B = rand(T,k,n)
    C = rand(T,m,n)
    Ab = big.(A)
    Bb = big.(B)
    C = A * B
    Cb = Ab * Bb
    e = (T <: Complex) ? eps(real(T)) : eps(T)
    @test opnorm(Cb - big.(C),1) / (opnorm(A,1) * opnorm(B,1)) < tol * k * e

    Ar = Matrix(reshape(A,k,m))
    Abr = Matrix(reshape(Ab,k,m))
    C = Ar' * B
    Cb = Abr' * Bb
    @test opnorm(Cb - big.(C),1) / (opnorm(A,1) * opnorm(B,1)) < tol * k * e
    Br = Matrix(reshape(B,n,k))
    Bbr = Matrix(reshape(Bb,n,k))
    C = Ar' * Br'
    Cb = Abr' * Bbr'
    @test opnorm(Cb - big.(C),1) / (opnorm(A,1) * opnorm(B,1)) < tol * k * e

    #FIXME: needs cases w/ simple transpose for complex types
end

@testset "gemm $T" for T in (Double64, Double32)
    # make sure to exercise the clean-up loops
    mA, nA, nB = 129, 67, 33
    gemmcheck(T,mA,nB,nA,1)
    # also run below MT threshold
    if Threads.nthreads() > 1
        t = DoubleBLAS.get_mt_threshold(:gemm)
        DoubleBLAS.set_mt_threshold(1.0e12,:gemm)
        gemmcheck(T,mA,nB,nA,1)
        DoubleBLAS.set_mt_threshold(t,:gemm)
    end
end

@testset "gemm $T" for T in (Complex{Double64}, Complex{Double32})
    # make sure to exercise the clean-up loops
    mA, nA, nB = 129, 67, 33
    gemmcheck(T,mA,nB,nA,1)
    # also run below MT threshold
    if Threads.nthreads() > 1
        t = DoubleBLAS.get_mt_threshold(:gemm)
        DoubleBLAS.set_mt_threshold(1.0e12,:gemm)
        gemmcheck(T,mA,nB,nA,1)
        DoubleBLAS.set_mt_threshold(t,:gemm)
    end
end

function gemvcheck(T,m,n,k,tol)
    A = rand(T,m,k)
    B = rand(T,k)
    C = rand(T,m)
    Ab = big.(A)
    Bb = big.(B)
    C = A * B
    Cb = Ab * Bb
    e = (T <: Complex) ? eps(real(T)) : eps(T)
    @test norm(Cb - big.(C),1) / (opnorm(A,1) * norm(B,1)) < tol * k * e

    Ar = Matrix(reshape(A,k,m))
    Abr = Matrix(reshape(Ab,k,m))
    C = Ar' * B
    Cb = Abr' * Bb
    @test norm(Cb - big.(C),1) / (opnorm(A,1) * norm(B,1)) < tol * k * e
end

@testset "gemv $T" for T in (Double64, Double32)
    # make sure to exercise the clean-up loops
    mA, nA, nB = 129, 67, 33
    gemvcheck(T,mA,nB,nA,1)
    # also run below MT threshold
    if Threads.nthreads() > 1
        t = DoubleBLAS.get_mt_threshold(:gemv)
        DoubleBLAS.set_mt_threshold(1.0e12,:gemv)
        gemvcheck(T,mA,nB,nA,1)
        DoubleBLAS.set_mt_threshold(t,:gemv)
    end
end

function lucheck(T,k,tol)
    A = rand(T,k,k)
    if (T <: Complex)
        κ = cond(ComplexF64.(A),Inf)
    else
        κ = cond(Float64.(A),Inf)
    end
    F = lu(A)
    xb = rand(T,k)
    b = A * xb
    x = F \ b
    e = eps(real(T))
    @test norm(x - xb, Inf) / norm(xb, Inf) < tol * κ * e

    Ai = inv(A)
    # CHECKME: why did I pick this norm?
    @test opnorm(Ai * A - I, 1) / opnorm(A, 1) < tol * k * e
end

@testset "lu $T" for T in (Double64, Double32, Complex{Double64})
    # make sure to exercise the clean-up loops
    nA = 67
    lucheck(T,nA,10)
    # also run below MT threshold
    if Threads.nthreads() > 1
        t = DoubleBLAS.get_mt_threshold(:lu)
        DoubleBLAS.set_mt_threshold(1.0e12,:lu)
        lucheck(T,nA,10)
        DoubleBLAS.set_mt_threshold(t,:lu)
    end
end

squareQ(Q::LinearAlgebra.AbstractQ) = (sq = size(Q.factors, 1); lmul!(Q, Matrix{eltype(Q)}(I, sq, sq)))

rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q)

function qrcheck(T,k,m,tol)
    A = rand(T,k,k)
    qra = qr(A)
    q, r = qra.Q, qra.R
    e = eps(real(T))
    @test norm(q' * squareQ(q) - I) < tol * k * e
    @test norm(A - q * r) / norm(A) < tol * k * e
    b = rand(T,k)
    @test norm(A * (qra \ b) - b) / norm(b) < tol * k * e
    A = rand(T,k,m)
    qra = qr(A)
    q, r = qra.Q, qra.R
    e = eps(real(T))
    @test norm(q' * squareQ(q) - I) < tol * k * e
    @test norm(q' * rectangularQ(q) - I) < tol * k * e
    @test norm(A - q * r) / norm(A) < tol * k * e
end

@testset "qr $T" for T in (Double64, Double32, Complex{Double64})
    # make sure to exercise the clean-up loops
    nA = 67
    mA = 49
    qrcheck(T,mA,nA,10)
end

function cholcheck(T,k,tol)
    A = rand(T,k,k)
    A = A * adjoint(A)
    # noncommutative multiplication means we need to be extra thorough
    A = A + adjoint(A)
    if (T <: Complex)
        κ = cond(ComplexF64.(A),Inf)
    else
        κ = cond(Float64.(A),Inf)
    end
    F = cholesky(A)
    xb = rand(T,k)
    b = A * xb
    x = F \ b
    e = eps(real(T))
    @test norm(x - xb,Inf) / norm(xb,Inf) < tol * κ * e

    C, info = LinearAlgebra._chol!(copy(A), LowerTriangular)
    F = LinearAlgebra.Cholesky(C.data,'L',info)
    x = F \ b
    @test norm(x - xb,Inf) / norm(xb,Inf) < tol * κ * e
end

# put this last since it depends on patches to DoubleFloats
@testset "chol $T" for T in (Double64, Double32, Complex{Double64})
    # make sure to exercise the clean-up loops
    nA = 67
    tol = 10
    cholcheck(T,nA,tol)
    # also run below MT threshold
    if Threads.nthreads() > 1
        t = DoubleBLAS.get_mt_threshold(:chol)
        DoubleBLAS.set_mt_threshold(1.0e12,:chol)
        cholcheck(T,nA,tol)
        DoubleBLAS.set_mt_threshold(t,:chol)
    end
end

include("refine.jl")


text:
# # Distribution of first return times for the one-dimensional random walk

# Here we verify that the tail of the distribution of the time of first return to the origin
# decays like $n^{-3/2}$

# First, we load packages that we will use
using RandomWalks, EmpiricalCDFs

# Set the number of samples (trajectories) to collect.
nsamples = 10^4;

# The following `AbstractActor` stops the walk on the first return, recording the number of steps taken.
first_return_actor = FirstReturnActor();

# `ECDFValueActor` builds an `EmpiricalCDF` from another (child) actor.
# In this case, we want to build an `EmpiricalCDF` of the first return time.
# If the child actor satisfies a simple interface, and `FirstReturnActor` does,
# then the entire specification of the action is simply
ecdf_actor = ECDFValueActor(first_return_actor);

# The probability that the walker returns to the origin is $1$, but the mean return time is infinite.
# Fluctuations are large and some walkers make a large number of steps before returning.
# So, we want to specify two criteria for terminating a walk: the first return and a limit on the number of steps.
# We put these `Actors` (and any others we may want) in an `ActorSet`.
plan_actor_set = ActorSet(StepLimitActor(10^9), first_return_actor);

# Now we run the trial, using a basic one-dimensional random walker `WalkB()`.

@time trial!(WalkPlan(WalkB(), plan_actor_set), SampleLoopActor(nsamples, ecdf_actor))

# Note that we have used `first_return_actor` in two places.
# First in the generation of the sample, and again in collecting data after each sample.

# Now let's assume (correctly) the tail of the first-return distribution follows a power law
# and verify that the exponent is close to $-3/2$.
using MaximumLikelihoodPower
mle(ecdf_actor[end-200:end]) # mle returns the estimate and error estimate

# ### Ensuring accurate statistics

# There is a potential problem with the example above.
# `StepLimitActor` breaks the sample loop silently, as it should.
# In the case at hand, with only `10^4` samples, the chance of one walker making $10^{9}$ steps before
# returning is very small. (How small ?)
# However, if we want to be sure that none of the walks were cut short,
# We can replace `StepLimitActor(10^9)` above by `ErrorOnExitActor(StepLimitActor(10^9)`.
# This replaces the silent exit with a fatal error.

# Another option is to use `CountActor` to count the number of times the step limit was reached.

step_limit_actor = StepLimitActor(10^6)
plan_actor_set = ActorSet(step_limit_actor, first_return_actor);
count_actor = CountActor(step_limit_actor)
trial_actor = SampleLoopActor(nsamples, ActorSet(ecdf_actor, count_actor))
@time trial!(WalkPlan(WalkB(), plan_actor_set), trial_actor)

get_actor_value(count_actor)

# Note that in the last example, we bound `step_limit_actor` to the instance of `StepLimitActor`
# because we need to refer to it in two places.


text:
using StochasticSearch, FactCheck

facts("[Configuration]") do
    context("constructors") do
        c = Configuration("test")
        p = [StringParameter("value",  "test1"),
             IntegerParameter(1, 5, 4, "test2")]
        c = Configuration(p, "test")
        @fact (c["test1"] == p[1])                   --> true
        @fact (c["test2"] == p[2])                   --> true
        p = Dict{String, Parameter}()
        p["test1"] = StringParameter("value",  "test")
        p["test2"] = IntegerParameter(1, 5, 4, "test")
        c = Configuration(p, "test")
        @fact (c["test1"] == p["test1"])             --> true
        @fact (c["test2"] == p["test2"])             --> true
        p = [StringParameter("valuea",  "test1"),
             StringParameter("valueb", "test2"),
             IntegerParameter(1, 44, 3, "test3")]
        c = Configuration(p, "test")
        @fact (c["test1"] == p[1])                   --> true
        @fact (c["test2"] == p[2])                   --> true
    end
    context("perturb!") do
        p = [StringParameter("valuea",  "test1"),
             StringParameter("valueb", "test2"),
             IntegerParameter(1, 44, 3, "test3")]
        t = Configuration(p, "test1")
        l = [IntegerParameter(1, 38, 3, "i0"),
             IntegerParameter(4, 66, 55, "i1"),
             StringParameter("value", "i2"),
             FloatParameter(2.33, 85.33, 22.2, "i3"),
             t]
        c = Configuration(l, "test")
        r = perturb!(c)
        r = perturb!(c)
        r = perturb!(c)
        r = perturb!(c)
        @fact (typeof(r["test1"]) <: Configuration)              --> true
        @fact (r["test1"]["test1"].value == "valuea")            --> true
        @fact (r["test1"]["test2"].value == "valueb")            --> true
        @fact (r["i2"].value             == "value")             --> true
        @fact (typeof(r["test1"]["test3"]) <: NumberParameter)   --> true
        @fact (typeof(r["i0"])             <: NumberParameter)   --> true
        @fact (typeof(r["i1"])             <: NumberParameter)   --> true
        @fact (typeof(r["i3"])             <: NumberParameter)   --> true
        d = Dict{String, Any}()
        d["i0"]             = 3
        d["i1"]             = 4
        d["i3"]             = 2.3
        d["test1"]          = Dict{String, Any}()
        d["test1"]["test3"] = 2
        @fact (typeof(perturb!(c, d)) <: Configuration)          --> true
        r = perturb!(c, d)
        r = perturb!(c, d)
        r = perturb!(c, d)
        r = perturb!(c, d)
        @fact (r["i0"].value == c["i0"].value)                   --> true
        @fact (r["i1"].value == c["i1"].value)                   --> true
        @fact (r["i3"].value == c["i3"].value)                   --> true
        @fact (r["test1"]["test3"] == c["test1"]["test3"])       --> true
        d["test1"]["test1"] = 23
        @fact_throws MethodError perturb!(c, d)
    end
    context("neighbor!") do
        l = [IntegerParameter(1, 38, 3, "i0"),
             IntegerParameter(4, 66, 55, "i1"),
             StringParameter("value", "i2"),
             FloatParameter(2.33, 85.33, 22.2, "i3")]
        d = Dict{String, Any}()
        c = Configuration(l, "test")
        println(c)
        d["i0"] = 5
        v = c["i0"].value
        neighbor!(c, d)
        @fact v                                              --> not(exactly(c["i0"].value))
        v = c["i0"].value
        neighbor!(c, d, distance = 9)
        d = Dict{String, Any}()
        d["i2"] = 5
        @fact_throws MethodError neighbor!(c, d)
        @fact_throws MethodError neighbor!(c, d, 2)
        l = [IntegerParameter(1, 38, 3, "i0")]
        c = Configuration(l, "test2")
        v = c["i0"].value
        neighbor!(c)
        @fact v                                              --> not(exactly(c["i0"].value))
    end
    context("update! and convert") do
        c     = Configuration("test")
        c["test"] = IntegerParameter(1, 2323, 2, "i0")
        a     = convert(Array{Parameter}, c)
        @fact (a[1] == c["test"])                            --> true
        push!(a, StringParameter("name", "i2"))
        update!(c, a)
        @fact (c["i2"] != a[1])                              --> true
        @fact (c["i2"].value == a[2].value)                  --> true
    end
end


text:
using Rankers.Storages
using Rankers.ProTeams
using LOLTools.SummonerV4: SummonerV4, http_action
using Calendars.Times

api_key = get(ENV, "RIOT_TOKEN", "")

region = "kr"
summonerName = "오뚜기3분미트볼"
summonerdto = SummonerV4.summoner_by_name(api_key, region, summonerName; action=store_with(region, http_action, debug=true))
@info millisecond2datetime(summonerdto.revisionDate)

for summoner in summonersof(players(DRX))
    summonerdto = SummonerV4.summoner_by_name(api_key, summoner.region, summoner.name; action=store_with(summoner.region, http_action, debug=true))
    @info millisecond2datetime(summonerdto.revisionDate)
end

for summoner in (summonersof ∘ players)(C9, ADC)
    summonerdto = SummonerV4.summoner_by_name(api_key, summoner.region, summoner.name; action=store_with(summoner.region, http_action, debug=true))
    @info millisecond2datetime(summonerdto.revisionDate)
end


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

using Base.Test

using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

srand(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2
a2real = randn(n,n)/2
a2img  = randn(n,n)/2

@testset for eltya in (Float32, Float64, Complex64, Complex128, Int)
    aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    aa2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
    asym = aa'+aa                  # symmetric indefinite
    apd  = aa'*aa                 # symmetric positive-definite
    @testset for atype in ("Array", "SubArray")
        if atype == "Array"
            a = aa
            a2 = aa2
        else
            a = view(aa, 1:n, 1:n)
            a2 = view(aa2, 1:n, 1:n)
        end
        ε = εa = eps(abs(float(one(eltya))))

        usv = svdfact(a)
        @testset "singular value decomposition" begin
            @test usv[:S] === svdvals(usv)
            @test usv[:U] * (Diagonal(usv[:S]) * usv[:Vt]) ≈ a
            @test AbstractArray(usv) ≈ a
            @test usv[:Vt]' ≈ usv[:V]
            @test_throws KeyError usv[:Z]
            b = rand(eltya,n)
            @test usv\b ≈ a\b

            if eltya <: BlasFloat
                svdz = svdfact!(ones(eltya,0,0))
                @test svdz[:U] ≈ eye(eltya,0,0)
                @test svdz[:S] ≈ real(zeros(eltya,0))
                @test svdz[:Vt] ≈ eye(eltya,0,0)
            end
        end
        @testset "Generalized svd" begin
            a_svd = a[1:n1, :]
            gsvd = svdfact(a,a_svd)
            @test gsvd[:U]*gsvd[:D1]*gsvd[:R]*gsvd[:Q]' ≈ a
            @test gsvd[:V]*gsvd[:D2]*gsvd[:R]*gsvd[:Q]' ≈ a_svd
            @test usv[:Vt]' ≈ usv[:V]
            @test_throws KeyError usv[:Z]
            @test_throws KeyError gsvd[:Z]
            @test gsvd[:vals] ≈ svdvals(a,a_svd)
            α = eltya == Int ? -1 : rand(eltya)
            β = svdfact(α)
            @test β[:S] == [abs(α)]
            @test svdvals(α) == abs(α)
            u,v,q,d1,d2,r0 = svd(a,a_svd)
            @test u ≈ gsvd[:U]
            @test v ≈ gsvd[:V]
            @test d1 ≈ gsvd[:D1]
            @test d2 ≈ gsvd[:D2]
            @test q ≈ gsvd[:Q]
            @test gsvd[:a].^2 + gsvd[:b].^2 ≈ ones(eltya,length(gsvd[:a]))

            #testing the other layout for D1 & D2
            b = rand(eltya,n,2*n)
            c = rand(eltya,n,2*n)
            gsvd = svdfact(b,c)
            @test gsvd[:U]*gsvd[:D1]*gsvd[:R]*gsvd[:Q]' ≈ b
            @test gsvd[:V]*gsvd[:D2]*gsvd[:R]*gsvd[:Q]' ≈ c
        end
    end
end


text:
module CCP

using Random: shuffle!
using Statistics: median
using StatsBase: sample, weights
using Serialization: serialize
using GZip
using JLD
using ArgParse
using Parameters

nbgens = 10000
nbruns = 4
nbsteps = 5000
popsize = 100
const beta = 1.
const a = 5.
const b = 3.
const maxcoop = 20
logdir = "juliatest2"

mutable struct Agent
    id::Int64
    fitness::Float64

    # Inputs
    seen_coop::Float64
    seen_punish::Float64
    new_partner::Bool

    # Outputs
    coop::Float64
    punish::Float64
    leave::Float64

    # inside
    coef::Float64

    # Layers
    layers::Array{Matrix{Float64}}
end

@with_kw struct Conf
    a::Float64 = 5
    b::Float64 = 3
    beta::Float64 = 1
    nbgens::Int = 20000
    nbruns::Int = 4
    nbsteps::Int = 5000
    popsize::Int = 100
    fakeprop::Int = 2
    maxcoop::Float64 = 30
    nbmut::Int = 1
    logdir::String = "."
    partnercontrol::Bool = false
    nointeract::Bool = false
end

Agent(id::Int) = Agent(id, 0, 0, 0, true, 0, 0, 0, 1, init_layers(-5., 5., [5]))
Agent(id::Int, layers::Array{Matrix{Float64}}) = Agent(id, 0, 0, 0, true, 0, 0, 0, 1, deepcopy(layers))
Agent(id::Int, parent::Agent) = Agent(id, parent.layers)

function Base.show(io::Base.IO, a::Agent)
    print(io, "Agent ", a.id)
end

function init_layers(low::Float64, high::Float64, hidden_dim::Array{Int64})::Array{Matrix{Float64}}
    dim::Array{Int64} = [4 hidden_dim 3]
    layers = Matrix{Float64}[]
    for (i, j) in zip(dim, dim[2:end])
        push!(layers, rand(Float64, (j, i)) .* (high - low) .+ low)
    end
    return layers
end

function decision(agent::Agent)::Agent
    output::Array{Float64} = [1.; agent.seen_coop / maxcoop;
                              agent.seen_punish / maxcoop;
                              convert(Float64, agent.new_partner)]
    for layer in agent.layers
        output = tanh.(layer * output)
    end
    agent.leave = output[1]
    agent.coop = (output[2] + 1) / 2 * maxcoop
    agent.punish = (output[3] + 1) / 2 * maxcoop
    return agent
end

function clear_inputs!(agent::Agent)::Agent
    agent.seen_coop = 0.
    agent.seen_punish = 0.
    agent.new_partner = true
    return agent
end

function create_pairs(pool, beta::Real)::Tuple{Array{Array{Agent}}, Array{Agent}}
    new_pairs = Array{Agent}[]
    alone_pool = Agent[]
    shuffle!(pool)
    for i = 1:2:length(pool)
        if rand(Float64) < beta
            push!(new_pairs, [pool[i], pool[i+1]])
        else
            push!(alone_pool, pool[i], pool[i+1])
        end
    end
    return new_pairs, alone_pool
end

function mutate!(agent::Agent)::Agent
    il = rand(1:length(agent.layers))
    i = rand(1:size(agent.layers[il])[1])
    j = rand(1:size(agent.layers[il])[2])
    agent.layers[il][i, j] = rand(Float32) * 10 - 5  # [-5, 5[
    return agent
end

isfake(agent::Agent)::Bool = agent.coef != 1
getcoop(agent::Agent)::Float64 = clamp(agent.coef * agent.coop, 0, maxcoop)
doleave(agent::Agent)::Bool = agent.leave > 0

function payoff(agent::Agent, pair, a::Real, b::Real)::Real
    x = getcoop(agent)
    x0 = sum(getcoop.(a for a in pair)) - x
    n = length(pair)
    return (a * (x + x0) + b * x0)/n - 1/2 * x*x
end

function main(conf::Conf)
    nbgens = conf.nbgens
    nbruns = conf.nbruns
    nbsteps = conf.nbsteps
    beta = conf.beta
    fakeprop = conf.fakeprop
    a = conf.a
    b = conf.b
    logdir = conf.logdir

    @assert (nbruns % fakeprop == 0) "nbruns should be a multiple of fakeprop so that each robot is fake the same amount of runs"

    mkpath(logdir)

    population = Agent[]
    nbfakes = div(popsize, fakeprop)
    GZip.open("$logdir/fitnesslog.txt.gz", "w") do ffit
            for agent in population
                write(ffit, "gen, agent, fitness\n")
            end
        end

    for i in 1:popsize
        push!(population, Agent(i))
    end
    for igen in 1:nbgens
        log = igen % 1000 == 0
        println(igen)
        if log
            JLD.save("$logdir/population_$igen.jld", "population", population)
            println("genome written")
            runlog = GZip.open("$logdir/jrunlog_$igen.txt.gz", "w")
            write(runlog, "gen, run, step, agent, pair, coef, seencoop, seenpunish, newpartner, coop, punish, leave\n")
        end
        fakes = zeros(Int, popsize)
        for irun in 1:nbruns
            if all(fakes .== 2)
                fakes = cat(zeros(Int, popsize - nbfakes), ones(Int, nbfakes), dims=1)
            else
                to_do = sample([i for i in 1:popsize if fakes[i] == 0], nbfakes, replace=false)
                fakes[to_do] .= 1
            end
            nbfakesdone = 0
            fakelistiter = collect(zip(population, fakes))
            shuffle!(fakelistiter)
            for (agent, fake) in fakelistiter
                if fake == 1
                    agent.coef = nbfakesdone / (nbfakes - 1) * 2
                    nbfakesdone += 1
                else
                    agent.coef = 1
                end
            end

            alone_pool = copy(population)
            pairs = Array{Agent}[]
            clear_inputs!.(population)
            for istep in 1:nbsteps
                newpairs, alone_pool = create_pairs(alone_pool, beta)
                push!(pairs, newpairs...)
                next_pairs = Array{Agent}[]
                leavers = Agent[]
                for agent in alone_pool
                    write(runlog, "$igen, $irun, $istep, -1, $(agent.id), $(agent.coef)," *
                                  "$(agent.seen_coop), $(agent.seen_punish), $(agent.new_partner)," *
                                  "$(agent.coop), $(agent.punish), $(agent.leave)\n")
                end
                if conf.nointeract
                    clear_inputs!.(population)
                end
                for (i_pair, pair) in enumerate(pairs)
                    for agent in pair
                        decision(agent)
                        if log
                            write(runlog, "$igen, $irun, $istep, $i_pair, $(agent.id), $(agent.coef)," *
                                  "$(agent.seen_coop), $(agent.seen_punish), $(agent.new_partner)," *
                                  "$(agent.coop), $(agent.punish), $(agent.leave)\n")
                        end
                    end
                    if !conf.partnercontrol && any([doleave(agent) for agent in pair])
                        clear_inputs!.([agent for agent in pair])
                        push!(leavers, pair...)
                    else
                        push!(next_pairs, pair)
                        total_coop::Float64 = 0.
                        total_spite::Float64 = 0.
                        for agent in pair
                            if !isfake(agent)
                                agent.fitness += payoff(agent, pair, a, b)
                                agent.fitness -= agent.punish
                            end
                            total_coop += getcoop(agent)
                            total_spite += agent.punish
                            for o_agent in pair
                                if o_agent != agent && !isfake(o_agent)
                                    o_agent.fitness -= 3 * agent.punish
                                end
                            end
                        end
                        for agent in pair
                            agent.new_partner = false
                            agent.seen_coop = total_coop - getcoop(agent)
                            agent.seen_punish = total_spite - agent.punish
                        end
                    end
                end
                pairs = next_pairs
                if length(leavers) > 0
                    push!(alone_pool, leavers...)
                end
            end # step
            fakes[fakes .== 1] .= 2
        end # run
        @assert all(fakes .== 2) "not all robots has been faked at the end"
        fitnesses = Float64[agent.fitness for agent in population]
        GZip.open("$logdir/fitnesslog.txt.gz", "a") do ffit
            for agent in population
                write(ffit, "$igen, $(agent.id), $(agent.fitness)\n")
            end
        end
        println(minimum(fitnesses), " ",  median(fitnesses), " ", maximum(fitnesses))
        minfit = minimum(fitnesses)
        if minfit <= 0
            wfit = weights(fitnesses .- (minfit .- 1))
        else
            wfit = weights(fitnesses)
        end
        parents = sample(population, wfit, popsize, replace=true)
        population::Array{Agent} = [Agent(i, parent) for (i, parent) in enumerate(parents)]
        for mutant in sample(population, conf.nbmut)
            mutate!(mutant)
        end
        if log
            close(runlog)
        end
    end  # gen
end

s = ArgParseSettings()
@add_arg_table s begin
    "--logdir"
        arg_type = String
        default = "."
    "--partnercontrol"
        action = :store_true
    "--nointeract"
        action = :store_true
    "--nbgens"
        arg_type = Int
        default = 10000
    "--nbruns"
        arg_type = Int
        default = 4
    "--nbsteps"
        arg_type = Int
        default = 5000
    "--popsize"
        arg_type = Int
        default = 100
    "--nbmut"
        arg_type = Int
        default = 1
    "--fakeprop"
        arg_type = Int
        default = 2
end

parsed_args = parse_args(ARGS, s)

conf = Conf(nbgens=parsed_args["nbgens"],
            logdir=parsed_args["logdir"],
            partnercontrol=parsed_args["partnercontrol"],
            nointeract=parsed_args["nointeract"],
            nbruns=parsed_args["nbruns"],
            nbsteps=parsed_args["nbsteps"],
            popsize=parsed_args["popsize"],
            nbmut=parsed_args["nbmut"],
            fakeprop=parsed_args["fakeprop"])
main(conf)

end

max_stars_repo_name: PaulEcoffet/controlchoicepunish

max_issues_repo_name: PaulEcoffet/controlchoicepunish

max_forks_repo_name: PaulEcoffet/controlchoicepunish


text:
using Plots, LaTeXStrings, Statistics, ProgressBars, JLD
cd(dirname(@__FILE__))
Plots.default(titlefontsize=12, tickfontsize=10, labelfontsize=10,
    fontfamily="Computer Modern")

function Deposition(; len, tot_time, dep_rate, color_step)
    surface = InitialSurf(len)
    AllLens = []
    for n in 0:tot_time
        randsurf = rand(1:len, dep_rate)
        for i in randsurf
            surface = check_and_depose(surface, i, len, (n%color_step) + 1)
        end
        append!(AllLens, Check_Dlen(surface))
    end
    return AllLens
end

function Check_Dlen(surf)
    Dlen = []
    for row in surf
        yindex = []
        for num in 1:length(row)
            if row[num] != 0
                append!(yindex, num)
            end
        end
        append!(Dlen, max(yindex...) - min(yindex...))
    end
    return max(Dlen...)
end

function check_and_depose(surf, index, len, addval) #checking down and sides and first layer
    i1 , i2 = sides(index, len)
    if surf[end][index] != 0
        newsurf = [0 for j = 1:len]
        newsurf[index] = addval
        push!(surf, newsurf)
    end
    for num in 0:length(surf) -1
        if surf[end - num][index] == 0 && (surf[end - num][i1] != 0 || surf[end - num][i2] != 0)
            surf[end - num][index] = addval
            break
        elseif length(surf) > num + 1 && surf[end - 1 - num][index] != 0
            surf[end - num][index] = addval
            break
        end
    end
    return surf
end

function sides(n, L)
    if n == L
        return n-1 , 1
    elseif n == 1
        return L , n+1
    else
        return n-1, n+1
    end
end

function InitialSurf(len)
    pos1 = [0 for i=1:len]
    pos1[floor(Int, len/2)] = 1
    return [pos1]
end

function Linear_fit(;Time, AvgList, steps)
    A = [hcat(Time) reshape(ones(steps), steps, 1)]
    b = reshape(AvgList, steps, 1)
    line = (A \ b)
    return line
end

Parameters = Dict(:len => 200,
                    :tot_time => 403,
                        :dep_rate => 100,
                            :color_step => 100)

up = 5.35; low = 1; steps = 100; SelctedSpan = 56:95; SelectedSteps = length(SelctedSpan);
Time = round.(Int,exp.(low:(up-low)/(steps-1):up))[SelctedSpan]
RawTime = hcat(low:(up-low)/(steps-1):up)[SelctedSpan]

iternum = 1000
allVal = [[0.0 for i ∈ 1:SelectedSteps] for j ∈ 1:iternum]
meanVal = [0.0 for i ∈ 1:SelectedSteps]
vars = [0.0 for i ∈ 1:SelectedSteps]
for i ∈ ProgressBar(1:iternum)
    VarList = Deposition(;Parameters...)[Time]
    allVal[i] = VarList
    meanVal += VarList
end
meanVal /= iternum
for i ∈ 1:SelectedSteps
    vars[i] = std(hcat(allVal...)[i,:])
end

Paraline = Dict(:Time => log.(Time),
                :AvgList => log.(meanVal),
                :steps => SelectedSteps)
Line = Linear_fit(;Paraline...)
X = RawTime
Y = X .* Line[1] .+ Line[2]

plot(X,Y,c = :black,label = L"y = %$(round(Line[1],digits= 3))x + %$(round(Line[2],digits= 3))")
scatter!(RawTime, log.(meanVal),
    c = :steelblue,
    xlabel= L"\mathrm{ln}(Time)",
    ylabel= L"\mathrm{ln}({L_C}_{(t)})",
    title= L"Log-Log\ Plot\ of\ L_C\ (L = %$(Parameters[:len]),\ %$iternum\ runs, %$(Time[end]*Parameters[:dep_rate])\ Particles)",
    label = L"Data\ point",
    ribbon = log.(vars),
    yerror = log.(vars),
    legend = 150)

savefig("../../Figs/Q3/Q3-L(t).pdf")
save("../../Data/Q3/Q3-TotData.jld","Data", allVal)


text:
module DataFramesMeta

using Reexport

@reexport using DataFrames

# Basics:
export @with, @where, @orderby, @transform, @by, @combine, @select, @eachrow,
       @transform!, @select!,
       @byrow, @byrow!, @based_on # deprecated

global const DATAFRAMES_GEQ_22 = isdefined(DataFrames, :pretty_table) ? true : false

include("linqmacro.jl")
include("eachrow.jl")

##############################################################################
##
## @with
##
##############################################################################

function addkey!(membernames, nam)
    if !haskey(membernames, nam)
        membernames[nam] = gensym()
    end
    membernames[nam]
end

onearg(e::Expr, f) = e.head == :call && length(e.args) == 2 && e.args[1] == f
onearg(e, f) = false

mapexpr(f, e) = Expr(e.head, map(f, e.args)...)

replace_syms!(x, membernames) = x
replace_syms!(q::QuoteNode, membernames) =
    replace_syms!(Meta.quot(q.value), membernames)
replace_syms!(e::Expr, membernames) =
    if onearg(e, :^)
        e.args[2]
    elseif onearg(e, :_I_)
        @warn "_I_() for escaping variables is deprecated, use cols() instead"
        addkey!(membernames, :($(e.args[2])))
    elseif onearg(e, :cols)
        addkey!(membernames, :($(e.args[2])))
    elseif e.head == :quote
        addkey!(membernames, Meta.quot(e.args[1]) )
    elseif e.head == :.
        replace_dotted!(e, membernames)
    else
        mapexpr(x -> replace_syms!(x, membernames), e)
    end

is_simple_non_broadcast_call(x) = false
function is_simple_non_broadcast_call(expr::Expr)
    expr.head == :call &&
        length(expr.args) >= 2 &&
        expr.args[1] isa Symbol &&
        all(x -> x isa QuoteNode || onearg(x, :cols), expr.args[2:end])
end

is_simple_broadcast_call(x) = false
function is_simple_broadcast_call(expr::Expr)
    expr.head == :. &&
        length(expr.args) == 2 &&
        expr.args[1] isa Symbol &&
        expr.args[2] isa Expr &&
        expr.args[2].head == :tuple &&
        all(x -> x isa QuoteNode || onearg(x, :cols), expr.args[2].args)
end

function args_to_selectors(v)
    t = map(v) do arg
        if arg isa QuoteNode
            arg
        elseif onearg(arg, :cols)
            arg.args[2]
        else
            throw(ArgumentError("This path should not be reached, arg: $(arg)"))
        end
    end

    :(DataFramesMeta.make_source_concrete($(Expr(:vect, t...))))
end

"""
    get_source_fun(function_expr)

Given an expression that may contain `QuoteNode`s (`:x`)
and items wrapped in `cols`, return a function
that is equivalent to that expression where the
`QuoteNode`s and `cols` items are the inputs
to the function.

For fast compilation `get_source_fun` returns
the name of a called function where possible.

* `f(:x, :y)` will return `f`
* `f.(:x, :y)` will return `ByRow(f)`
* `:x .+ :y` will return `.+`

`get_source_fun` also returns an expression
representing the vector of inputs that will be
used as the `src` in the `src => fun => dest`
call later on.

### Examples

julia> using MacroTools

julia> ex = :(:x + :y)

julia> DataFramesMeta.get_source_fun(ex)
(:(DataFramesMeta.make_source_concrete([:x, :y])), :+)

julia> ex = quote
           :x .+ 1 .* :y
       end |> MacroTools.prettify

julia> src, fun = DataFramesMeta.get_source_fun(ex);

julia> MacroTools.prettify(fun)
:((mammoth, goat)->mammoth .+ 1 .* goat)

"""
function get_source_fun(function_expr)
    # recursive step for begin :a + :b end
    if function_expr isa Expr &&
        function_expr.head == :block &&
        length(function_expr.args) == 2 # omitting the line number node

        return get_source_fun(function_expr.args[2])
    elseif is_simple_non_broadcast_call(function_expr)
        source = args_to_selectors(function_expr.args[2:end])
        fun_t = function_expr.args[1]

        # .+ to +
        if startswith(string(fun_t), '.')
            f_sym_without_dot = Symbol(chop(string(fun_t), head = 1, tail = 0))
            fun = :(DataFrames.ByRow($f_sym_without_dot))
        else
            fun = fun_t
        end

        return source, fun
    elseif is_simple_broadcast_call(function_expr)
        # extract source symbols from quotenodes
        source = args_to_selectors(function_expr.args[2].args)
        fun_t = function_expr.args[1]
        fun = :(DataFrames.ByRow($fun_t))

        return source, fun
    else
        membernames = Dict{Any, Symbol}()

        body = replace_syms!(function_expr, membernames)

        source = :(DataFramesMeta.make_source_concrete($(Expr(:vect, keys(membernames)...))))
        inputargs = Expr(:tuple, values(membernames)...)

        fun = quote
            $inputargs -> begin
                $body
            end
        end

        return source, fun
    end
end

"""
    @col(kw)

`@col` transforms an expression of the form `z = :x + :y` into it's equivalent in
DataFrames's `source => fun => destination` syntax.

### Details

Parsing follows the same convention as other DataFramesMeta.jl macros, such as `@with`. All
terms in the expression that are `Symbol`s are treated as columns in the data frame, except
`Symbol`s wrapped in `^`. To use a variable representing a column name, wrap the variable
in `cols`.

`@col` constructs an anonymous function `fun` based on the given expression. It then creates
a `source => fun => destination` pair that is suitable for the `select`, `transform`, and
`combine` functions in DataFrames.jl.

### Examples

```julia
julia> @col z = :x + :y
[:x, :y] => (##595 => :z)
```

In the above example, `##595` is an anonymous function equivalent to the following

```julia
(_x, _y) -> _x + _y
```

```jldoctest
julia> using DataFramesMeta;

julia> df = DataFrame(x = [1, 2], y = [3, 4]);

julia> import DataFramesMeta: @col;

julia> DataFrames.transform(df, @col z = :x .* :y)
2×3 DataFrame
│ Row │ x     │ y     │ z     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 3     │ 3     │
│ 2   │ 2     │ 4     │ 8     │

```
"""
macro col(kw)
    esc(fun_to_vec(kw))
end

# `nolhs` needs to be `true` when we have syntax of the form
# `@combine(gd, fun(:x, :y))` where `fun` returns a `table` object.
# We don't create the "new name" pair because new names are
# given by the table.
function fun_to_vec(ex::Expr; nolhs::Bool = false, gensym_names::Bool = false)
    # classify the type of expression
    # :x # handled via dispatch
    # cols(:x) # handled as though above
    # f(:x) # nohls == true, re-write as simple call
    # (; a = :x, ) # nolhs == true, complicated call
    # y = :x # :x is a QuoteNode
    # y = cols(:x) # use cols on RHS
    # cols(:y) = :x # RHS in :block
    # cols(:y) = cols(:x) #
    # y = f(:x) # re-write as simple call
    # y = f(cols(:x)) # re-write as simple call, use cols
    # y = :x + 1 # re-write as complicated call
    # y = cols(:x) + 1 # re-write as complicated call, with cols
    # cols(:y) = f(:x) # re-write as simple call, but RHS is :block
    # cols(:y) = f(cols(:x)) # re-write as simple call, RHS is block, use cols
    # cols(y) = :x + 1 # re-write as complicated col, but RHS is :block
    # cols(:y) = cols(:x) + 1 # re-write as complicated call, RHS is block, use cols

    if gensym_names
        ex = Expr(:kw, gensym(), ex)
    end

    nokw = (ex.head !== :(=)) && (ex.head !== :kw) && nolhs

    # :x
    # handled below via dispatch on ::QuoteNode

    # cols(:x)
    if onearg(ex, :cols)
        return ex.args[2]
    end

    # The above cases are the only ones allowed
    # if you don't have nolhs explicitely stated
    # or are just `:x` or `cols(x)`
    if !(ex.head === :(=) || ex.head === :kw || nolhs)
        throw(ArgumentError("Expressions not of the form `y = f(:x)` are currently disallowed."))
    end

    # f(:x) # it's assumed this returns a Table
    # (; a = :x, ) # something more explicit we might see
    if nokw
        source, fun = get_source_fun(ex)

        return quote
            $source => $fun => AsTable
        end
    end

    if !nokw
        lhs = ex.args[1]
        rhs_t = ex.args[2]
        # if lhs is a cols(y) then the rhs gets parsed as a block
        if onearg(lhs, :cols) && rhs_t.head === :block && length(rhs_t.args) == 2
            rhs = rhs_t.args[2]
        else
            rhs = rhs_t
        end
    else
        throw(ArgumentError("This path should not be reached"))
    end

    # y = :x
    if lhs isa Symbol && rhs isa QuoteNode
        source = rhs
        dest = QuoteNode(lhs)

        return quote
            $source => $dest
        end
    end

    # y = cols(:x)
    if lhs isa Symbol && onearg(rhs, :cols)
        source = rhs.args[2]
        dest = QuoteNode(lhs)

        return quote
            $source => $dest
        end
    end

    # cols(:y) = :x
    if onearg(lhs, :cols) && rhs isa QuoteNode
        source = rhs
        dest = lhs.args[2]

        return quote
            $source => $dest
        end
    end

    # cols(:y) = cols(:x)
    if onearg(lhs, :cols) && onearg(rhs, :cols)
        source = rhs.args[2]
        dest = lhs.args[2]

        return quote
            $source => $dest
        end
    end

    # y = f(:x)
    # y = f(cols(:x))
    # y = :x + 1
    # y = cols(:x) + 1
    if lhs isa Symbol
        source, fun = get_source_fun(rhs)
        dest = QuoteNode(lhs)

        return quote
            $source => $fun => $dest
        end
    end

    # cols(:y) = f(:x)
    if onearg(lhs, :cols)
        source, fun = get_source_fun(rhs)
        dest = lhs.args[2]

        return quote
            $source => $fun => $dest
        end
    end

    throw(ArgumentError("This path should not be reached"))
end
fun_to_vec(ex::QuoteNode; nolhs::Bool = false, gensym_names::Bool = false) = ex

function make_source_concrete(x::AbstractVector)
    if isempty(x) || isconcretetype(eltype(x))
        return x
    elseif all(t -> t isa Union{AbstractString, Symbol}, x)
        return Symbol.(x)
    else
        throw(ArgumentError("Column references must be either all the same " *
                            "type or a a combination of `Symbol`s and strings"))
    end
end

protect_replace_syms!(e, membernames) = e
function protect_replace_syms!(e::Expr, membernames)
    if e.head == :quote
        e
    else
        replace_syms!(e, membernames)
    end
end

function replace_dotted!(e, membernames)
    x_new = replace_syms!(e.args[1], membernames)
    y_new = protect_replace_syms!(e.args[2], membernames)
    Expr(:., x_new, y_new)
end

exec(df, v, fun) = fun(map(c -> DataFramesMeta.getsinglecolumn(df, c), v)...)

getsinglecolumn(df, s::DataFrames.ColumnIndex) = df[!, s]
getsinglecolumn(df, s) = throw(ArgumentError("Only indexing with Symbols, strings and integers " *
    "is currently allowed with cols"))

function with_helper(d, body)
    source, fun = get_source_fun(body)
    :(DataFramesMeta.exec($d, $source, $fun))
end

"""
    @with(d, expr)

`@with` allows DataFrame columns keys to be referenced as symbols.

### Arguments

* `d` : an AbstractDataFrame type
* `expr` : the expression to evaluate in `d`

### Details

`@with` works by parsing the expression body for all columns indicated
by symbols (e.g. `:colA`). Then, a function is created that wraps the
body and passes the columns as function arguments. This function is
then called. Operations are efficient because:

- A pseudo-anonymous function is defined, so types are stable.
- Columns are passed as references, eliminating DataFrame indexing.

The following

```julia
@with(d, :a .+ :b .+ 1)
```

becomes

```julia
tempfun(a, b) = a .+ b .+ 1
tempfun(d[!, :a], d[!, :b])
```

If an expression is wrapped in `^(expr)`, `expr` gets passed through untouched.
If an expression is wrapped in  `cols(expr)`, the column is referenced by the
variable `expr` rather than a symbol.

### Examples

```jldoctest
julia> using DataFramesMeta

julia> y = 3;

julia> df = DataFrame(x = 1:3, y = [2, 1, 2]);

julia> x = [2, 1, 0];

julia> @with(df, :y .+ 1)
3-element Array{Int64,1}:
 3
 2
 3

julia> @with(df, :x + x)
3-element Array{Int64,1}:
 3
 3
 3

julia> @with df begin
            res = 0.0
            for i in 1:length(:x)
                res += :x[i] * :y[i]
            end
            res
        end
10.0

julia> @with(df, df[:x .> 1, ^(:y)]) # The ^ means leave the :y alone
2-element Array{Int64,1}:
 1
 2

julia> colref = :x;

julia> @with(df, :y + cols(colref)) # Equivalent to df[!, :y] + df[!, colref]
3-element Array{Int64,1}:
 3
 3
 5
```

!!! note
    `@with` creates a function, so the scope within `@with` is a local scope.
    Variables in the parent can be read. Writing to variables in the parent scope
    differs depending on the type of scope of the parent. If the parent scope is a
    global scope, then a variable cannot be assigned without using the `global` keyword.
    If the parent scope is a local scope (inside a function or let block for example),
    the `global` keyword is not needed to assign to that parent scope.
"""
macro with(d, body)
    esc(with_helper(d, body))
end

##############################################################################
##
## @where - select row subsets
##
##############################################################################

function where_helper(x, args...)
    t = (fun_to_vec(arg; nolhs = true, gensym_names = true) for arg in args)
    quote
        $where($x, $(t...))
    end
end

function df_to_bool(res::AbstractDataFrame)
    if any(t -> !(t isa AbstractVector{<:Union{Missing, Bool}}), eachcol(res))
        throw(ArgumentError("All arguments in @where must return an " *
                            "AbstractVector{<:Union{Missing, Bool}}"))
    end

    return reduce((x, y) -> x .& y, eachcol(res)) .=== true
end

function where(df::AbstractDataFrame, @nospecialize(args...))
    res = DataFrames.select(df, args...; copycols = false)
    tokeep = df_to_bool(res)
    df[tokeep, :]
end

function where(gd::GroupedDataFrame, @nospecialize(args...))
    res = DataFrames.select(gd, args...; copycols = false, keepkeys = false)
    tokeep = df_to_bool(res)
    parent(gd)[tokeep, :]
end

function where(df::SubDataFrame, @nospecialize(args...))
    res = DataFrames.select(df, args...)
    tokeep = df_to_bool(res)
    df[tokeep, :]
end

"""
    @where(d, i...)

Select row subsets in `AbstractDataFrame`s and `GroupedDataFrame`s.

### Arguments

* `d` : an AbstractDataFrame or GroupedDataFrame
* `i...` : expression for selecting rows

Multiple `i` expressions are "and-ed" together.

If given a `GroupedDataFrame`, `@where` applies transformations by
group, and returns a fresh `DataFrame` containing the rows
for which the generated values are all `true`.

!!! note
    `@where` treats `missing` values as `false` when filtering rows.
    Unlike `DataFrames.filter` and other boolean operations with
    `missing`, `@where` will *not* error on missing values, and
    will only keep `true` values.

### Examples

```jldoctest
julia> using DataFramesMeta, Statistics

julia> df = DataFrame(x = 1:3, y = [2, 1, 2]);

julia> globalvar = [2, 1, 0];

julia> @where(df, :x .> 1)
2×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 2     │ 1     │
│ 2   │ 3     │ 2     │

julia> @where(df, :x .> globalvar)
2×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 2     │ 1     │
│ 2   │ 3     │ 2     │

julia> @where(df, :x .> globalvar, :y .== 3)
0×2 DataFrame

julia> d = DataFrame(n = 1:20, x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1,
                                    2, 1, 1, 2, 2, 2, 3, 1, 1, 2]);

julia> g = groupby(d, :x);

julia> @where(g, :n .> mean(:n))
8×2 DataFrame
│ Row │ n     │ x     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 12    │ 1     │
│ 2   │ 13    │ 1     │
│ 3   │ 15    │ 2     │
│ 4   │ 16    │ 2     │
│ 5   │ 17    │ 3     │
│ 6   │ 18    │ 1     │
│ 7   │ 19    │ 1     │
│ 8   │ 20    │ 2     │

julia> @where(g, :n .== first(:n))
3×2 DataFrame
│ Row │ n     │ x     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 3     │
│ 2   │ 5     │ 1     │
│ 3   │ 8     │ 2     │

julia> d = DataFrame(a = [1, 2, missing], b = ["x", "y", missing]);

julia> @where(d, :a .== 1)
1×2 DataFrame
│ Row │ a      │ b       │
│     │ Int64? │ String? │
├─────┼────────┼─────────┤
│ 1   │ 1      │ x       │
```
"""
macro where(x, args...)
    esc(where_helper(x, args...))
end

##############################################################################
##
## @orderby
##
##############################################################################

function orderby_helper(x, args...)
    t = (fun_to_vec(arg; nolhs = true, gensym_names = true) for arg in args)
    quote
        $DataFramesMeta.orderby($x, $(t...))
    end
end

function orderby(x::AbstractDataFrame, @nospecialize(args...))
    t = DataFrames.select(x, args...; copycols = false)
    x[sortperm(t), :]
end

function orderby(x::GroupedDataFrame, @nospecialize(args...))
    throw(ArgumentError("@orderby with a GroupedDataFrame is reserved"))
end

function orderby(x::SubDataFrame, @nospecialize(args...))
    t = DataFrames.select(x, args...)
    x[sortperm(t), :]
end

"""
    @orderby(d, i...)

Sort rows by values in one of several columns or a transformation of columns.
Always returns a fresh `DataFrame`. Does not accept a `GroupedDataFrame`.

When given a `DataFrame`, `@orderby` applies the transformation
given by its arguments (but does not create new columns) and sorts
the given `DataFrame` on the result, returning a new `DataFrame`.

### Arguments

* `d` : an AbstractDataFrame
* `i...` : expression for sorting

### Examples

```jldoctest
julia> using DataFramesMeta, Statistics

julia> d = DataFrame(x = [3, 3, 3, 2, 1, 1, 1, 2, 1, 1], n = 1:10,
                     c = ["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"]);

julia> @orderby(d, -1 .* :n)
10×3 DataFrame
 Row │ x      n      c
     │ Int64  Int64  String
─────┼──────────────────────
   1 │     1     10  j
   2 │     1      9  i
   3 │     2      8  h
   4 │     1      7  g
   5 │     1      6  f
   6 │     1      5  e
   7 │     2      4  d
   8 │     3      3  c
   9 │     3      2  b
  10 │     3      1  a

julia> @orderby(d, sortperm(:c, rev = true))
10×3 DataFrame
 Row │ x      n      c
     │ Int64  Int64  String
─────┼──────────────────────
   1 │     1     10  j
   2 │     1      9  i
   3 │     2      8  h
   4 │     1      7  g
   5 │     1      6  f
   6 │     1      5  e
   7 │     2      4  d
   8 │     3      3  c
   9 │     3      2  b
  10 │     3      1  a

julia> @orderby(d, :x, abs.(:n .- mean(:n)))
10×3 DataFrame
 Row │ x      n      c
     │ Int64  Int64  String
─────┼──────────────────────
   1 │     1      5  e
   2 │     1      6  f
   3 │     1      7  g
   4 │     1      9  i
   5 │     1     10  j
   6 │     2      4  d
   7 │     2      8  h
   8 │     3      3  c
   9 │     3      2  b
  10 │     3      1  a

"""
macro orderby(d, args...)
    esc(orderby_helper(d, args...))
end

##############################################################################
##
## transform & @transform
##
##############################################################################

function transform_helper(x, args...)

    t = (fun_to_vec(arg) for arg in args)

    quote
        $DataFrames.transform($x, $(t...))
    end
end

"""
    @transform(d, i...)

Add additional columns or keys based on keyword arguments.

### Arguments

* `d` : an `AbstractDataFrame`, or `GroupedDataFrame`
* `i...` : keyword arguments defining new columns or keys

### Returns

* `::AbstractDataFrame` or `::GroupedDataFrame`

### Examples

```jldoctest
julia> using DataFramesMeta

julia> df = DataFrame(A = 1:3, B = [2, 1, 2]);

julia> @transform(df, a = 2 * :A, x = :A .+ :B)
3×4 DataFrame
│ Row │ A     │ B     │ a     │ x     │
│     │ Int64 │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 2     │ 3     │
│ 2   │ 2     │ 1     │ 4     │ 3     │
│ 3   │ 3     │ 2     │ 6     │ 5     │
```
"""
macro transform(x, args...)
    esc(transform_helper(x, args...))
end

##############################################################################
##
## transform! & @transform!
##
##############################################################################

function transform!_helper(x, args...)

    t = (fun_to_vec(arg) for arg in args)

    quote
        $DataFrames.transform!($x, $(t...))
    end
end

"""
    @transform!(d, i...)

Mutate `d` inplace to add additional columns or keys based on keyword arguments and return it.
No copies of existing columns are made.

### Arguments

* `d` : an `AbstractDataFrame`, or `GroupedDataFrame`
* `i...` : keyword arguments defining new columns or keys

### Returns

* `::DataFrame`

### Examples

```jldoctest
julia> using DataFramesMeta

julia> df = DataFrame(A = 1:3, B = [2, 1, 2]);

julia> df2 = @transform!(df, a = 2 * :A, x = :A .+ :B)
3×4 DataFrame
│ Row │ A     │ B     │ a     │ x     │
│     │ Int64 │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┼───────┤
│ 1   │ 1     │ 2     │ 2     │ 3     │
│ 2   │ 2     │ 1     │ 4     │ 3     │
│ 3   │ 3     │ 2     │ 6     │ 5     │

julia> df === df2
true
```
"""
macro transform!(x, args...)
    esc(transform!_helper(x, args...))
end

##############################################################################
##
## @combine - summarize a grouping operation
##
##############################################################################

function combine_helper(x, args...; deprecation_warning = false)
    deprecation_warning && @warn "`@based_on` is deprecated. Use `@combine` instead."

    # Only allow one argument when returning a Table object
    if length(args) == 1 &&
        !(first(args) isa QuoteNode) &&
        !(first(args).head == :(=) || first(args).head == :kw)

        t = fun_to_vec(first(args); nolhs = true)
        # 0.22: No pair as first arg, needs AsTable in other args to return table
        if DATAFRAMES_GEQ_22
            quote
                $DataFrames.combine($x, $t)
            end
        # 0.21: Pair as first arg, other args can't return table
        else
            quote
                $DataFrames.combine($t, $x)
            end
        end
    else
        t = (fun_to_vec(arg) for arg in args)
        quote
            $DataFrames.combine($x, $(t...))
        end
    end
end

"""
    @combine(x, args...)

Summarize a grouping operation

### Arguments

* `x` : a `GroupedDataFrame` or `AbstractDataFrame`
* `args...` : keyword arguments defining new columns

### Examples

```jldoctest
julia> using DataFramesMeta

julia> d = DataFrame(
            n = 1:20,
            x = [3, 3, 3, 3, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 3, 1, 1, 2]);

julia> g = groupby(d, :x);

julia> @combine(g, nsum = sum(:n))
3×2 DataFrame
│ Row │ x     │ nsum  │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 3     │ 27    │
│ 2   │ 1     │ 99    │
│ 3   │ 2     │ 84    │

julia> @combine(g, x2 = 2 * :x, nsum = sum(:n))
20×3 DataFrame
│ Row │ x     │ x2    │ nsum  │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 3     │ 6     │ 27    │
│ 2   │ 3     │ 6     │ 27    │
│ 3   │ 3     │ 6     │ 27    │
│ 4   │ 3     │ 6     │ 27    │
│ 5   │ 3     │ 6     │ 27    │
│ 6   │ 1     │ 2     │ 99    │
│ 7   │ 1     │ 2     │ 99    │
⋮
│ 13  │ 1     │ 2     │ 99    │
│ 14  │ 1     │ 2     │ 99    │
│ 15  │ 2     │ 4     │ 84    │
│ 16  │ 2     │ 4     │ 84    │
│ 17  │ 2     │ 4     │ 84    │
│ 18  │ 2     │ 4     │ 84    │
│ 19  │ 2     │ 4     │ 84    │
│ 20  │ 2     │ 4     │ 84    │

```
"""
macro combine(x, args...)
    esc(combine_helper(x, args...))
end

"""
    @based_on(d, args...)

Deprecated version of `@combine`, see: [`@combine`](@ref)
"""
macro based_on(x, args...)
    esc(combine_helper(x, args...; deprecation_warning = true))
end

##############################################################################
##
## @by - grouping
##
##############################################################################

function by_helper(x, what, args...)
    # Only allow one argument when returning a Table object
    if length(args) == 1 &&
        !(first(args) isa QuoteNode) &&
        !(first(args).head == :(=) || first(args).head == :kw)

        t = fun_to_vec(first(args); nolhs = true)
        # 0.22: No pair as first arg, needs AsTable in other args to return table
        if DATAFRAMES_GEQ_22
            quote
                $DataFrames.combine($groupby($x, $what), $t)
            end
        # 0.21: Pair as first arg, other args can't return table
        else
            quote
                $DataFrames.combine($t, $groupby($x, $what))
            end
        end
    else
        t = (fun_to_vec(arg) for arg in args)
        quote
            $DataFrames.combine($groupby($x, $what), $(t...))
        end
    end
end

"""
    @by(d::AbstractDataFrame, cols, e...)

Split-apply-combine in one step.

### Arguments

* `d` : an AbstractDataFrame
* `cols` : a column indicator (Symbol, Int, Vector{Symbol}, etc.)
* `e` :  keyword arguments specifying new columns in terms of column groupings

### Returns

* `::DataFrame`

### Examples

```jldoctest
julia> using DataFramesMeta, Statistics

julia> df = DataFrame(
            a = repeat(1:4, outer = 2),
            b = repeat(2:-1:1, outer = 4),
            c = 1:8);

julia> @by(df, :a, d = sum(:c))
4×2 DataFrame
│ Row │ a     │ d     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 6     │
│ 2   │ 2     │ 8     │
│ 3   │ 3     │ 10    │
│ 4   │ 4     │ 12    │

julia> @by(df, :a, d = 2 * :c)
8×2 DataFrame
│ Row │ a     │ d     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 2     │
│ 2   │ 1     │ 10    │
│ 3   │ 2     │ 4     │
│ 4   │ 2     │ 12    │
│ 5   │ 3     │ 6     │
│ 6   │ 3     │ 14    │
│ 7   │ 4     │ 8     │
│ 8   │ 4     │ 16    │

julia> @by(df, :a, c_sum = sum(:c), c_mean = mean(:c))
4×3 DataFrame
│ Row │ a     │ c_sum │ c_mean  │
│     │ Int64 │ Int64 │ Float64 │
├─────┼───────┼───────┼─────────┤
│ 1   │ 1     │ 6     │ 3.0     │
│ 2   │ 2     │ 8     │ 4.0     │
│ 3   │ 3     │ 10    │ 5.0     │
│ 4   │ 4     │ 12    │ 6.0     │

julia> @by(df, :a, c = :c, c_mean = mean(:c))
8×3 DataFrame
│ Row │ a     │ c     │ c_mean  │
│     │ Int64 │ Int64 │ Float64 │
├─────┼───────┼───────┼─────────┤
│ 1   │ 1     │ 1     │ 3.0     │
│ 2   │ 1     │ 5     │ 3.0     │
│ 3   │ 2     │ 2     │ 4.0     │
│ 4   │ 2     │ 6     │ 4.0     │
│ 5   │ 3     │ 3     │ 5.0     │
│ 6   │ 3     │ 7     │ 5.0     │
│ 7   │ 4     │ 4     │ 6.0     │
│ 8   │ 4     │ 8     │ 6.0     │

```
"""
macro by(x, what, args...)
    esc(by_helper(x, what, args...))
end

##############################################################################
##
## @select - select and transform columns
##
##############################################################################

function select_helper(x, args...)
    t = (fun_to_vec(arg) for arg in args)

    quote
        $DataFrames.select($x, $(t...))
    end
end

"""
    @select(d, e...)

Select and transform columns.

### Arguments

* `d` : an `AbstractDataFrame` or `GroupedDataFrame`
* `e` :  keyword arguments specifying new columns in terms of existing columns
  or symbols to specify existing columns

### Returns

* `::AbstractDataFrame`

### Examples

```jldoctest
julia> using DataFrames, DataFramesMeta

julia> df = DataFrame(a = repeat(1:4, outer = 2), b = repeat(2:-1:1, outer = 4), c = 1:8);

julia> @select(df, :c, :a)
8×2 DataFrame
│ Row │ c     │ a     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 1     │
│ 2   │ 2     │ 2     │
│ 3   │ 3     │ 3     │
│ 4   │ 4     │ 4     │
│ 5   │ 5     │ 1     │
│ 6   │ 6     │ 2     │
│ 7   │ 7     │ 3     │
│ 8   │ 8     │ 4     │

julia> @select(df, :c, x = :b + :c)
8×2 DataFrame
│ Row │ c     │ x     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 3     │
│ 2   │ 2     │ 3     │
│ 3   │ 3     │ 5     │
│ 4   │ 4     │ 5     │
│ 5   │ 5     │ 7     │
│ 6   │ 6     │ 7     │
│ 7   │ 7     │ 9     │
│ 8   │ 8     │ 9     │
```
"""
macro select(x, args...)
    esc(select_helper(x, args...))
end

##############################################################################
##
## @select! - in-place select and transform columns
##
##############################################################################

function select!_helper(x, args...)
    t = (fun_to_vec(arg) for arg in args)

    quote
        $DataFrames.select!($x, $(t...))
    end
end

"""
    @select!(d, e...)

Mutate `d` in-place to retain only columns or transformations specified by `e` and return it. No copies of existing columns are made.

### Arguments

* `d` : an AbstractDataFrame
* `e` :  keyword arguments specifying new columns in terms of existing columns
  or symbols to specify existing columns

### Returns

* `::DataFrame`

### Examples

```jldoctest
julia> using DataFrames, DataFramesMeta

julia> df = DataFrame(a = repeat(1:4, outer = 2), b = repeat(2:-1:1, outer = 4), c = 1:8);

julia> df2 = @select!(df, :c, :a)
8×2 DataFrame
│ Row │ c     │ a     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 1     │
│ 2   │ 2     │ 2     │
│ 3   │ 3     │ 3     │
│ 4   │ 4     │ 4     │
│ 5   │ 5     │ 1     │
│ 6   │ 6     │ 2     │
│ 7   │ 7     │ 3     │
│ 8   │ 8     │ 4     │

julia> df === df2
true

julia> df = DataFrame(a = repeat(1:4, outer = 2), b = repeat(2:-1:1, outer = 4), c = 1:8);

julia> df2 = @select!(df, :c, x = :b + :c)
8×2 DataFrame
│ Row │ c     │ x     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 3     │
│ 2   │ 2     │ 3     │
│ 3   │ 3     │ 5     │
│ 4   │ 4     │ 5     │
│ 5   │ 5     │ 7     │
│ 6   │ 6     │ 7     │
│ 7   │ 7     │ 9     │
│ 8   │ 8     │ 9     │

julia> df === df2
true
```
"""
macro select!(x, args...)
    esc(select!_helper(x, args...))
end

end # module


text:
# This is a Node in the Binning Analysis tree that averages two values. There
# is one of these for each binning level. When two values should be
# compressed, this is done immediately, so that only one value needs to be saved.
# switch indicates whether value should be written to or averaging should happen.
mutable struct Compressor{T}
    value::T
    switch::Bool
end

struct LogBinner{T, N}
    # list of Compressors, one per level
    compressors::NTuple{N, Compressor{T}}

    # sum(x) for all values on a given lvl
    x_sum::Vector{T}
    # sum(x.^2) for all values on a given lvl
    x2_sum::Vector{T}
    # number of values that are summed on a given lvl
    count::Vector{Int64}
end

# Overload some basic Base functions
Base.eltype(B::LogBinner{T,N}) where {T,N} = T
Base.length(B::LogBinner) = B.count[1]
Base.ndims(B::LogBinner{T,N}) where {T,N} = ndims(eltype(B))
Base.isempty(B::LogBinner) = length(B) == 0
Base.:(==)(a::T, b::T) where {T <: Compressor} = (a.value == b.value) && (a.switch == b.switch)
Base.:(!=)(a::T, b::T) where {T <: Compressor} = !(a == b)

function Base.:(==)(a::LogBinner{T, N}, b::LogBinner{T, M}) where {T, N, M}
    # Switch order so that we can deal with just N ≤ M here
    (N > M) && (return b == a)

    # Does every level match?
    for i in 1:N
        (a.compressors[i] == b.compressors[i]) || return false
        (a.x_sum[i] == b.x_sum[i]) || return false
        (a.x2_sum[i] == b.x2_sum[i]) || return false
        (a.count[i] == b.count[i]) || return false
    end

    # Are extra levels empty?
    for i in N+1:M
        (b.count[i] == 0) || return false
    end

    return true
end
Base.:(!=)(a::LogBinner{T, N}, b::LogBinner{T, M}) where {T, N, M} = !(a == b)

function _print_header(io::IO, B::LogBinner{T,N}) where {T,N}
    print(io, "LogBinner{$(T),$(N)}")
    nothing
end

function _println_body(io::IO, B::LogBinner{T,N}) where {T,N}
    n = length(B)
    println(io)
    print(io, "| Count: ", n)
    if n > 0 && ndims(B) == 0
        print(io, "\n| Mean: ", round.(mean(B), digits=5))
        print(io, "\n| StdError: ", round.(std_error(B), digits=5))
    end
    nothing
end

# short version (shows up in arrays etc.)
Base.show(io::IO, B::LogBinner{T,N}) where {T,N} = print(io, "LogBinner{$(T),$(N)}()")
# verbose version (shows up in the REPL)
Base.show(io::IO, m::MIME"text/plain", B::LogBinner) = (_print_header(io, B); _println_body(io, B))

"""
    empty!(B::LogBinner)

Clear the binner, i.e. reset it to its inital state.
"""
function Base.empty!(B::LogBinner)
    !isempty(B) || return

    # get zero
    z = zero(eltype(eltype(B)))

    # reset x_sum and x2_sum to all zeros
    @inbounds for i in eachindex(B.x_sum)
        if ndims(B) == 0 # compile time
            # numbers
            B.x_sum[i] = z
            B.x2_sum[i] = z
        else
            # arrays
            fill!(B.x_sum[i], z)
            fill!(B.x2_sum[i], z)
        end
    end

    # reset counts
    fill!(B.count, 0)

    # reset compressors
    for c in B.compressors
        if ndims(B) == 0 # compile time
            c.value = z
        else
            fill!(c.value, z)
        end
        c.switch = false
    end

    nothing
end

_nlvls2capacity(N::Int) = 2^N - 1
_capacity2nlvls(capacity::Int) = ceil(Int, log(2, capacity + 1))

"""
    capacity(B)

Capacity of the binner, i.e. how many values can be handled before overflowing.
"""
capacity(B::LogBinner{T,N}) where {T,N} = _nlvls2capacity(N)
nlevels(B::LogBinner{T,N}) where {T,N} = N

"""
    LogBinner([::Type{T}; capacity::Int])

Creates a `LogBinner` which can handle (at least) `capacity` many values of type `T`.

The default is `T = Float64` and `capacity = 2^32-1 ≈ 4e9`.
"""
LogBinner(::Type{T} = Float64; kw...) where T = LogBinner(zero(T); kw...)

# TODO
# Currently, the Binning Analysis requires a "zero" to initialize x_sum and
# x2_sum. It's not necessary for the Compressors.
# Since Arrays do not have a static size, we cannot generate a fitting zero
# automatically. So currently we let the user supply it.
# Other possibilities:
# > generate x_sum, x2_sum with #undef
# bad: requires many checks (if isassigned ... else ...)
# > force StaticArrays (and/or tuples)
# good: faster for small arrays
# bad: requires user to use another package
# > initialize is first push!
# bad: also requires frequent checks (if first push ... else ...)
# ...?
"""
    LogBinner(zero_element::T[; capacity::Int])

Creates a new `LogBinner` which can take (at least) `capacity` many values of type `T`. The type
and the size are inherited from the given `zero_element`, which must exclusively contain zeros.

Values can be added using `push!` and `append!`.

---

    LogBinner(timeseries::AbstractVector{T})

Creates a new `LogBinner` and adds all elements from the given timeseries.
"""
function LogBinner(x::T;
        capacity::Int64 = _nlvls2capacity(32)
        ) where {T <: Union{Number, AbstractArray}}

    # check keyword args
    capacity <= 0 && throw(ArgumentError("`capacity` must be finite and positive."))

    # got_timeseries = didn't receive a zero && is a vector
    got_timeseries = count(!iszero, x) > 0 && ndims(T) == 1

    if got_timeseries
        # x = time series
        N = _capacity2nlvls(length(x))
        S = _sum_type_heuristic(eltype(T), ndims(x[1]))
        el = zero(x[1])
    else
        # x = zero_element
        N = _capacity2nlvls(capacity)
        S = _sum_type_heuristic(T, ndims(x))
        el = x
    end

    B = LogBinner{S, N}(
        tuple([Compressor{S}(copy(el), false) for i in 1:N]...),
        [copy(el) for _ in 1:N],
        [copy(el) for _ in 1:N],
        zeros(Int64, N)
    )

    got_timeseries && append!(B, x)

    return B
end

function _sum_type_heuristic(::Type{T}, elndims::Integer) where T
    # heuristic to set sum type (#2)
    S = if eltype(T)<:Real
        elndims > 0 ? Array{Float64, elndims} : Float64
    else
        elndims > 0 ? Array{ComplexF64, elndims} : ComplexF64
    end
    return S
end

"""
    LogBinner(B::LogBinner[; capacity::Int])

Creates a new `LogBinner` from an existing LogBinner, copying the data inside.
The new LogBinner may be larger or smaller than the given one.
"""
function LogBinner(B::LogBinner{S, M}; capacity::Int64 = _nlvls2capacity(32)) where {S, M}
    N = _capacity2nlvls(capacity)
    B.count[min(M, N)] > 1 && throw(OverflowError(
        "The new LogBinner is too small to reconstruct the given LogBinner. " *
        "New capacity = $capacity   Old capacity = $(B.count[1])"
    ))
    el = zero(B.x_sum[1])

    LogBinner{S, N}(
        tuple([i > M ? Compressor{S}(copy(el), false) : deepcopy(B.compressors[i]) for i in 1:N]...),
        [i > M ? copy(el) : copy(B.x_sum[i]) for i in 1:N],
        [i > M ? copy(el) : copy(B.x2_sum[i]) for i in 1:N],
        [i > M ? 0 : copy(B.count[i]) for i in 1:N]
    )
end

# TODO typing?
"""
    append!(LogBinner, values)

Adds an array of values to the binner by `push!`ing each element.
"""
function Base.append!(B::LogBinner, values::AbstractArray)
    @inbounds for i in eachindex(values)
        _push!(B, 1, values[i])
    end
    nothing
end

"""
    push!(LogBinner, value)

Pushes a new value into the Binning Analysis.
"""
function Base.push!(B::LogBinner{T,N}, value::S) where {N, T, S}
    ndims(T) == ndims(S) || throw(DimensionMismatch("Expected $(ndims(T)) dimensions but got $(ndims(S))."))

    _push!(B, 1, value)
end

@inline _square(x) = x^2
@inline _square(x::Complex) = Complex(real(x)^2, imag(x)^2)
@inline _square(x::AbstractArray) = _square.(x)

# recursion, back-end function
@inline function _push!(B::LogBinner{T,N}, lvl::Int64, value::S) where {N, T <: Number, S}
    C = B.compressors[lvl]

    # any value propagating through this function is new to lvl. Therefore we
    # add it to the sums. Note that values pushed to the output arrays are not
    # added here until the array drops to the next level. (New compressors are
    # added)
    B.x_sum[lvl] += value
    B.x2_sum[lvl] += _square(value)
    B.count[lvl] += 1

    if !C.switch
        # Compressor has space -> save value
        C.value = value
        C.switch = true
        return nothing
    else
        # Do averaging
        if lvl == N
            # No more propagation possible -> throw error
            throw(OverflowError("The Binning Analysis has exceeded its maximum capacity."))
        else
            # propagate to next lvl
            C.switch = false
            _push!(B, lvl+1, 0.5 * (C.value + value))
            return nothing
        end
    end
    return nothing
end

function _push!(
        B::LogBinner{T,N},
        lvl::Int64,
        value::S
    ) where {N, T <: AbstractArray, S}

    C = B.compressors[lvl]
    B.x_sum[lvl] .+= value
    B.x2_sum[lvl] .+= _square(value)
    B.count[lvl] += 1

    if !C.switch
        C.value .= value
        C.switch = true
        return nothing
    else
        if lvl == N
            throw(OverflowError("The Binning Analysis has exceeded its maximum capacity."))
        else
            C.switch = false
            _push!(B, lvl+1, 0.5 * (C.value .+ value))
            return nothing
        end
    end
    return nothing
end


text:
module GridapGmsh

using Libdl
using Gridap
using Gridap.CellValuesGallery
using StaticArrays
using UnstructuredGrids.Kernels

export GmshDiscreteModel
export GmshCLagrangianFESpace
export GmshDLagrangianFESpace

using Gridap.DiscreteModels: DiscreteModelFromData

deps_jl = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(deps_jl)
  s = """
  Package GridapGmsh not installed properly.
  Run Pkg.build(\"GridapGmsh\"), restart Julia and try again
  """
  error(s)
end

include(deps_jl)

include(gmsh_jl)

# Hack taken from MPI.jl
function __init__()
    @static if Sys.isunix()
        Libdl.dlopen(gmsh.lib, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
    end
end

include("GmshDiscreteModels.jl")

include("LagrangianFESpaces.jl")

end # module


text:
# Partial JSON API implementation.
# Author: Chuan-Zheng Lee

module JsonAPI

using JSON

export jsonapidict, writejsonapi, jsonapi

typealias JsonDict Dict{AbstractString,Any}

makename(s::Symbol) = lowercase(string(s))

makevalue(x::Integer) = x
makevalue(x::AbstractFloat) = x
makevalue(x::AbstractString) = x
makevalue(x::Void) = x
makevalue(x::Enum) = Integer(x)

addfield!(res::JsonDict, k::Symbol, v) = addfield!(res, makename(k), v)

addfield!(res::JsonDict, k::AbstractString, v::Integer) = addattribute!(res, k, v)
addfield!(res::JsonDict, k::AbstractString, v::AbstractFloat) = addattribute!(res, k, v)
addfield!(res::JsonDict, k::AbstractString, v::AbstractString) = addattribute!(res, k, v)
addfield!(res::JsonDict, k::AbstractString, v::Void) = addattribute!(res, k, v)
addfield!(res::JsonDict, k::AbstractString, v::Enum) = addattribute!(res, k, v)

addfield!{T<:Integer}(res::JsonDict, k::AbstractString, v::Array{T}) = addattributearray!(res, k, v)
addfield!{T<:AbstractFloat}(res::JsonDict, k::AbstractString, v::Array{T}) = addattributearray!(res, k, v)
addfield!{T<:AbstractString}(res::JsonDict, k::AbstractString, v::Array{T}) = addattributearray!(res, k, v)
addfield!{T<:Void}(res::JsonDict, k::AbstractString, v::Array{T}) = addattributearray!(res, k, v)
addfield!{T<:Enum}(res::JsonDict, k::AbstractString, v::Array{T}) = addattributearray!(res, k, v)

function addfield!(res::JsonDict, k::AbstractString, v)
    if :id ∉ fieldnames(v)
        error("Not sure what to do with field $k, type $(typeof(v)) has no field id")
    end
    addrelationship!(res, k, v)
end

function addfield!(res::JsonDict, k::AbstractString, v::Array)
    if :id ∉ fieldnames(eltype(v))
        error("Not sure what to do with field $k, type $(eltype(v)) has no field id")
    end
    addrelationshiparray!(res, k, v)
end

function addattribute!(res::JsonDict, k::AbstractString, v)
    attributes = get!(JsonDict, res, "attributes")
    attributes[k] = makevalue(v)
end

function addattributearray!(res::JsonDict, k::AbstractString, v::Array)
    attributes = get!(JsonDict, res, "attributes")
    attributes[k] = map(makevalue, v)
end

function addrelationship!(res::JsonDict, k::AbstractString, v)
    relationships = get!(JsonDict, res, "relationships")
    resource = resourceid(v)
    linkage = JsonDict("data"=>resource)
    relationships[k] = linkage
end

function addrelationshiparray!(res::JsonDict, k::AbstractString, v)
    relationships = get!(JsonDict, res, "relationships")
    resources = map(resourceid, v)
    linkage = JsonDict("data"=>resources)
    relationships[k] = linkage
end

function resourceid{T}(obj::T)
    res = JsonDict()
    res["type"] = makename(T.name.name)
    res["id"] = obj.id
    return res
end

function resource{T}(obj::T, id::Int)
    res = resource(obj)
    if !haskey(res, "id")
        res["id"] = id
    end
    return res
end

function resource{T}(obj::T)
    res = JsonDict()
    res["type"] = makename(T.name.name)
    if :id ∈ fieldnames(T)
        res["id"] = obj.id
    end
    for field in fieldnames(T)
        if field == :id
            continue
        end
        fvalue = getfield(obj, field)
        addfield!(res, field, fvalue)
    end
    return res
end

function primarydata(a::Array)
    primarydata = Array{JsonDict}(length(a))
    for (i, item) in enumerate(a)
        primarydata[i] = resource(item, i)
    end
    return primarydata
end

function primarydata(obj)
    resource(obj, 1)
end

function jsonapidict(a)
    d = JsonDict()
    d["data"] = primarydata(a)
    return d
end

function writejsonapi(io::IO, obj)
    d = jsonapidict(obj)
    JSON.print(io, d)
end

function jsonapi(obj)
    d = jsonapidict(obj)
    JSON.json(d)
end

end # module


text:
"""
`Label` is the root abstract type used to provide unique identifiers for variables within
`UncertainValues` and `LabeledValues`. See [BasicLabel]
"""
abstract type Label end

"""
`BasicLabel{T}` is a mechanism for creating `Label` around other objects.
For example, BasicLabel{String} would create a `Label` using a `String`
to provide a unique identifier.
"""
struct BasicLabel{T} <: Label
    value::T
    hc::UInt # cache the hash to ensure Dict and Set searches are fast

    BasicLabel(v::T) where {T} = new{T}(v, hash(v, UInt(0xEA7F00DBADF00D)))
end

label(item::T) where {T} = BasicLabel(item)

macro nl_str(str::AbstractString)
    BasicLabel(str)
end

Base.show(io::IO, bl::BasicLabel) = print(io, "$(bl.value)")
Base.isequal(bl1::BasicLabel{T}, bl2::BasicLabel{T}) where {T} =
    isequal(bl1.hc, bl2.hc) && isequal(bl1.value, bl2.value)
Base.isequal(l1::Label, l2::Label) = false
Base.hash(bl::BasicLabel, h::UInt) = hash(bl.hc, h)

"""
    labelsByType(ty::Type{<:Label}, labels::AbstractVector{<:Label})
    labelsByType(types::AbstractVector{DataType}, labels::AbstractVector{<:Label})

Extracts all of a specific type of `Label` from a list of `Label`s. The second version
extracts multiple different types in a single call.
"""
labelsByType(ty::Type{<:Label}, labels::AbstractVector{<:Label}) =
    filter(lbl -> typeof(lbl) == ty, labels)

labelsByType(types::AbstractVector{DataType}, labels::AbstractVector{<:Label}) =
    mapreduce(ty -> labelsByType(ty, labels), append!, types)


text:
using Test

using Dice
using Dice: CuddMgr, DistBool

@testset "DistBool Tests" begin

    mgr = CuddMgr()

    t = DistBool(mgr, true)
    f = DistBool(mgr, false)
    x = DistBool(mgr, 0.5)

    @test Dice.isvalid(prob_equals(x & t, x))
    @test Dice.isvalid(prob_equals(x & f, f))
    @test Dice.isvalid(prob_equals(x | t, t))
    @test Dice.isvalid(prob_equals(x | f, x))

    @test Dice.isvalid(prob_equals(x & true, x))
    @test Dice.isvalid(prob_equals(true & x, x))
    @test Dice.isvalid(prob_equals(x & false, f))
    @test Dice.isvalid(prob_equals(false & x, f))
    @test Dice.isvalid(prob_equals(x | true, t))
    @test Dice.isvalid(prob_equals(true | x, t))
    @test Dice.isvalid(prob_equals(x | false, x))
    @test Dice.isvalid(prob_equals(false | x, x))

    @test Dice.isvalid(prob_equals(x & x, x))
    @test Dice.isvalid(prob_equals(x | x, x))

end


text:
using Polynomials

hasPositiveRealPart(x::Number) = (real(x) > 0)

pp(x::Number) = max(x,0) # positive part

function findRootMonotoneFunction(f)
  # find root t of f such that t >= 0
  # assume f is a monotonically increasing function such that f(0) < 0

  if f(0.0) >= 0
    error("monotonically increasing function with f(0) >= 0 has no positive root.")
  end
  t0 = 0.0

  prec = 1.0e-10

  # find t2 such that p(t2) > 0
  # may cause infinite loop if p not increasing
  t2 = 1.0
  while (f(t2) <= 0)
    t2 *= 2
  end
  t1 = t2/2

  while abs(f(t1)) > prec
    if f(t1) < 0
      t0 = t1
    else
      t2 = t1
    end
    t1 = t0 + (t2-t0)/2
  end

  return t1

end

function findRootMonotonePolynomial(p)
  # find root t of p such that t >= 0
  # assume f is a monotonically increasing polynomial such that p(0) < 0

  theroots = filter(isreal,roots(p))
  if isempty(r)
    error("No root found")
  else
    maximum(real(theroots))
  end

end

function findRootPolyMaxZero(p::Polynomials.Poly, val::Real)
  # for m(t) = (p(t))^+, find t such that \int_0^t m(s) \ d s = val
  # assume p(t) is monotonically increasing polynomial for t >= 0
  if p(0) >= 0
    # println("using direct method, p(0) >= 0, p'(0) = ", polyder(p)(0), ", value = ", val)
    # return findRootMonotonePolynomial(polyint(p) - val)
    return findRootMonotoneFunction(polyint(p) - val)
  else
    println("shifting")
    t0 = findRootMonotoneFunction(p) # p(t) >= 0 for t >= t0 by assumption
    P_prelim = polyint(p)
    f = t -> P_prelim(t+t0) - P_prelim(t0) - val
    return findRootMonotoneFunction(f) + t0
    # P = P_prelim - P_prelim(t0) - val
    # println("P(x) = ", P)
    # return findRootMonotonePolynomial(P)
  end
end

function equidistantSamples(T::Vector{Float64},X::Array{Float64,2}, n_samples::Int = length(T))

  n_skeletonpoints = length(T)
  dim = size(X)[1]
  t_max = T[end]
  dt = t_max / (n_samples-1)
  samples = Array{Float64,2}(dim,n_samples)
  t_current = dt
  t0 = T[1]
  x0 = X[:,1]
  n_sampled = 0
  for i = 1:n_skeletonpoints-1
    x1 = X[:,i+1]
    t1 = T[i+1]
    while (t_current < t1 && n_sampled < n_samples)
      samples[:,n_sampled + 1] = x0 + (x1-x0) *(t_current - t0)/(t1-t0)
      n_sampled += 1
      t_current += dt
    end
    x0 = x1
    t0 = t1
  end
  # t_travelled = 0
  # current_skeleton_index = 1
  # delta_skeleton_time = T[current_skeleton_index + 1] - T[current_skeleton_index]
  # for i = 1:n_samples
  #   while ((current_skeleton_index + 1) * delta_skeleton_time < stepsize * i)
  #     t_travelled = (current_skeleton_index + 1) * delta_skeleton_time
  #     current_skeleton_index += 1
  #     if current_skeleton_index < n_skeletonpoints
  #       delta_skeleton_time = T[current_skeleton_index + 1] - T[current_skeleton_index]
  #     else
  #       # reached the end, now what?
  #     end
  #   end # of while
  #
  #   samples[i] = X[:,current_skeleton_index] + (stepsize * i - t_travelled)/(delta_skeleton_time) * (X[:,current_skeleton_index+1] - X[:,current_skeleton_index])
  #   t_travelled += stepsize * i
  # end
  #
  #
  #
  return samples[:,1:n_sampled]

end

switchingtime = function (a,b,u)
# generate switching time for rate of the form max(0, a + b s)
  if (b > 0)
    if (a < 0)
      return -a/b + switchingtime(0, b, u);
    else # a >= 0
      return -a/b + sqrt(a^2/b^2 - 2 * log(u)/b);
    end
  elseif (b == 0) # degenerate case
    if (a < 0)
      return Inf;
    else # a >= 0
      return -log(u)/a;
    end
  else # b <= 0
    if (a <= 0)
      return Inf;
    else # a > 0
      y = -log(u); t1=-a/b;
      if (y >= a * t1 + b *t1^2/2)
        return Inf;
      else
        return -a/b - sqrt(a^2/b^2 + 2 * y /b);
      end
    end
  end
end

max_stars_repo_name: jbierkens/zigzag-experiments

max_issues_repo_name: jbierkens/zigzag-experiments

max_forks_repo_name: jbierkens/zigzag-experiments


text:
module YetAnotherHelperStickyJellyfish

using ..Ahorn, Maple

@mapdef Entity "YetAnotherHelper/StickyJellyfish" StickyJellyfish(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Sticky Jellyfish (Yet Another Helper)" => Ahorn.EntityPlacement(
        StickyJellyfish,
        "point"
    )
)

sprite = "ahorn/YetAnotherHelper/stickyJellyfish"

function Ahorn.selection(entity::StickyJellyfish)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y - 5)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::StickyJellyfish, room::Maple.Room)
    Ahorn.drawSprite(ctx, sprite, 0, -5)
end

end

max_stars_repo_name: YetAnotherHelper/YetAnotherHelper

max_issues_repo_name: YetAnotherHelper/YetAnotherHelper

max_forks_repo_name: YetAnotherHelper/YetAnotherHelper


text:
@testset "Capillary pressure" begin
    f = reshape(collect(1.0:25),5,5)
    f_float = reshape(collect(1.0f0:25.0f0),5,5)
    res = zeros(5,5)
    @testset "No contact angle" begin
        Swalbe.filmpressure!(res, f, 1.0, 0.0, 3, 2, 0.1, 0.1)
        sol = [-30.0 -5.0 -5.0 -5.0 20;
               -25.0 0.0 0.0 0.0 25.0;
               -25.0 0.0 0.0 0.0 25.0;
               -25.0 0.0 0.0 0.0 25.0;
               -20.0 5.0 5.0 5.0 30.0;
               ]
        for i in eachindex(sol)
            @test res[i] .≈ sol[i] atol=1e-10
        end
        Swalbe.filmpressure!(res, f, 0.0)
        for i in eachindex(sol)
            @test res[i] .≈ 0.01 * sol[i] atol=1e-10
        end
    end
    @testset "No height gradient" begin
        nograd = ones(5,5)
        Swalbe.filmpressure!(res, nograd, 1.0, 1/2, 3, 2, 0.1, 0.0)
        for i in eachindex(res)
            @test res[i] .≈ -2(0.1^2-0.1) atol=1e-10
        end
    end
    @testset "Gradient and contact angle" begin
        Swalbe.filmpressure!(res, f, 1.0, 1/2, 3, 2, 0.1, 0.0)
        sol = [30.0 5.0 5.0 5.0 -20;
               25.0 0.0 0.0 0.0 -25.0;
               25.0 0.0 0.0 0.0 -25.0;
               25.0 0.0 0.0 0.0 -25.0;
               20.0 -5.0 -5.0 -5.0 -30.0]
        for i in eachindex(res)
            # Now the result comprisses of two components the disjoining potential and the laplace term.
            @test res[i] .≈ -1(sol[i] + 20((0.1/f[i])^3-(0.1/f[i])^2)) atol=1e-10
        end
    end
    @testset "Gradient and contact angle Float32" begin
        Swalbe.filmpressure!(res, f_float, 1.0f0, 0.5f0, 3, 2, 0.1f0, 0.0f0)
        sol = [30.0f0 5.0f0 5.0f0 5.0f0 -20f0;
               25.0f0 0.0f0 0.0f0 0.0f0 -25.0f0;
               25.0f0 0.0f0 0.0f0 0.0f0 -25.0f0;
               25.0f0 0.0f0 0.0f0 0.0f0 -25.0f0;
               20.0f0 -5.0f0 -5.0f0 -5.0f0 -30.0f0]

        for i in eachindex(res)
            # Now the result comprisses of two components the disjoining potential and the laplace term.
            @test res[i] .≈ -1(sol[i] + 20((0.1f0/f_float[i])^3-(0.1f0/f_float[i])^2)) atol=1e-6
        end
    end
    dgrad = zeros(5,5,8)
    @testset "No contact angle circshift!" begin
        Swalbe.filmpressure!(res, f, dgrad, 1.0, 0.0, 3, 2, 0.1, 0.1)
        sol = [-30.0 -5.0 -5.0 -5.0 20;
               -25.0 0.0 0.0 0.0 25.0;
               -25.0 0.0 0.0 0.0 25.0;
               -25.0 0.0 0.0 0.0 25.0;
               -20.0 5.0 5.0 5.0 30.0;
               ]
        for i in eachindex(sol)
            @test res[i] .≈ sol[i] atol=1e-10
        end
    end
end

@testset "Capillary pressure 1D" begin
    f = collect(1.0:30)
    f_float = collect(1.0f0:30.0f0)
    res = zeros(30)
    dummy = zeros(30,3)
    @testset "No contact angle" begin
        Swalbe.filmpressure!(res, f, dummy, 1.0, 0.0, 3, 2, 0.1, 0.1)
        # println("My result: $res")
        sol = zeros(30)
        sol[1] = 30
        sol[end] = -30
        @test all(res .== -sol)
    end
    @testset "No height gradient" begin
        nograd = ones(30)
        Swalbe.filmpressure!(res, nograd, dummy, 1.0, 1/2, 3, 2, 0.1, 0.0)
        for i in eachindex(res)
            @test res[i] .≈ -2(0.1^2-0.1) atol=1e-10
        end
    end
    @testset "Gradient and contact angle" begin
        Swalbe.filmpressure!(res, f, dummy, 1.0, 1/2, 3, 2, 0.1, 0.0)
        sol = zeros(30)
        sol[1] = 30
        sol[end] = -30
        for i in eachindex(res)
            # Now the result comprisses of two components the disjoining potential and the laplace term.
            @test res[i] .≈ -1(sol[i] + 20((0.1/f[i])^3-(0.1/f[i])^2)) atol=1e-10
        end
    end
end

@testset "Capillary pressure 1D with Gamma" begin
    f = collect(1.0:30)
    res = zeros(30)
    rho = fill(0.1, 30)
    dummy = zeros(30,2)
    @testset "No Gamma no theta" begin
        Swalbe.filmpressure!(res, f, dummy, rho, 1.0, 0.0, 3, 2, 0.1, 0.1)
        # println("My result: $res")
        sol = zeros(30)
        sol[1] = 30
        sol[end] = -30
        @test all(res .== -sol)
    end
    @testset "No Gamma no laplacian" begin
        nograd = ones(30)
        Swalbe.filmpressure!(res, nograd, dummy, rho, 1.0, 1/2, 3, 2, 0.1, 0.0)
        for i in eachindex(res)
            @test res[i] .≈ -2(0.1^2-0.1) atol=1e-10
        end
    end
    @testset "Gamma but rho=0 no laplacian" begin
        nograd = ones(30)
        Swalbe.filmpressure!(res, nograd, dummy, zeros(30), 1.0, 1/2, 3, 2, 0.1, 0.0)
        for i in eachindex(res)
            @test res[i] .≈ -2(0.1^2-0.1) atol=1e-10
        end
    end
    @testset "Gamma, rho, gradient and contact angle" begin
        rho = fill(0.1, 30)
        Swalbe.filmpressure!(res, f, dummy, rho, 1.0, 1/2, 3, 2, 0.1, 0.0, Gamma=0.1)
        sol = zeros(30)
        sol[1] = 30
        sol[end] = -30
        for i in eachindex(res)
            # Now the result comprisses of two components the disjoining potential and the laplace term.
            # The prefactor (1 + 0.01) -> (gamma + Gamma * rho)
            @test res[i] .≈ -(1 + 0.01) * (sol[i] + 20((0.1/f[i])^3-(0.1/f[i])^2)) atol=1e-10
        end
    end
end

@testset "Power function" begin
    arg = 2
    argsimple = 5.0
    argsimplefp32 = 5.0f0
    argvec = [2.0 3.0 4.0]
    n = [3 6]
    r = [8 64]
    r1 = [125.0 15625.0]
    r2 = [125.0f0 15625.0f0]
    r3 = [8.0 27.0 64.0; 64.0 729.0 4096.0]
    for i in 1:2
        res0 = Swalbe.power_broad(arg, n[i])
        @test res0 == r[i]
        res0 = Swalbe.power_broad(argsimple, n[i])
        @test res0 == r1[i]
        res1 = Swalbe.power_broad(argsimplefp32, n[i])
        @test res1 == r2[i]
        res2 = Swalbe.power_broad.(argvec, n[i])
        for j in 1:3
            @test res2[j] .== r3[i,j]
        end
    end
end