text:
import Base.LinAlg:
    A_mul_B!, At_mul_B, A_ldiv_B!

import Tulip:
    Model, PrimalDualPoint

"""
    solve(model, tol, verbose)

Solve model m using an infeasible predictor-corrector Interior-Point algorithm.

# Arguments
- `model::Model`: the optimization model
- `tol::Float64`: numerical tolerance
- `verbose::Int`: 0 means no output, 1 displays log at each iteration
"""
function solve!(
    model::Model;
    tol::Float64 = 10.0^-8,
    verbose::Int = 0
)

    model.status = :Built

    N_ITER_MAX = 100  # maximum number of IP iterations
    niter = 0  # number of IP iterations

    # TODO: pre-optimization stuff
    F = symbolic_cholesky(model.A)  # Symbolic factorization
    θ = zeros(model.sol.x)

    # X = Array{PrimalDualPoint, 1}()
    # TODO: check stopping criterion for possible early termination

    # compute starting point
    # TODO: decide which starting point
    compute_starting_point!(
        model.A,
        F,
        model.sol,
        model.b, model.c, model.uind, model.uval
    )

    # IPM log
    if verbose == 1
        println(" Itn      Primal Obj        Dual Obj    Prim Inf Dual Inf UBnd Inf\n")
    end

    # main loop
    # push!(X, copy(model.sol))

    while niter < N_ITER_MAX

        # I. Form and factor Newton System
        compute_newton!(
            model.A,
            model.sol.x,
            model.sol.s,
            model.sol.w,
            model.sol.z,
            model.uind,
            θ,
            F
        )

        # II. Compute and take step
        compute_next_iterate!(model, F)

        # III. Book-keeping + display log
        # compute residuals
        rb = model.A * model.sol.x - model.b
        rc = At_mul_B(model.A, model.sol.y) + model.sol.s - model.c
        spxpay!(-1.0, rc, model.uind, model.sol.z)

        ru = model.sol.x[model.uind] + model.sol.w - model.uval

        obj_primal = dot(model.sol.x, model.c)
        obj_dual = dot(model.b, model.sol.y) - dot(model.uval, model.sol.z)

        niter += 1

        # check stopping criterion
        eps_p = (norm(rb)) / (1.0 + norm(model.b))
        eps_d = (norm(rc)) / (1.0 + norm(model.c))
        eps_u = (norm(ru)) / (1.0 + norm(model.uval))
        eps_g = abs(obj_primal - obj_dual) / (1.0 + abs(obj_primal))
        # if verbose == 1
        #     print("\teps_p=")
        #     print(@sprintf("%9.2e", eps_p))
        #     print("\teps_d=")
        #     print(@sprintf("%9.2e", eps_d))
        #     print("\teps_u=")
        #     print(@sprintf("%9.2e", eps_u))
        #     print("\teps_g=")
        #     print(@sprintf("%9.2e", eps_g))
        #     println("\n")
        # end
        if verbose == 1
            print(@sprintf("%4d", niter))  # iteration count
            print(@sprintf("%+18.7e", obj_primal))  # primal objective
            print(@sprintf("%+16.7e", obj_dual))  # dual objective
            print(@sprintf("%10.2e", norm(rb, Inf)))  # primal infeas
            print(@sprintf("%9.2e", norm(rc, Inf)))  # dual infeas
            print(@sprintf("%9.2e", norm(ru, Inf)))  # upper bound infeas
            print(@sprintf("%9.2e", abs(obj_primal - obj_dual) / (model.nvars + size(model.uind, 1))))
            print("\n")
        end

        if (eps_p < tol) && (eps_u < tol) && (eps_d < tol) && (eps_g < tol)
            model.status = :Optimal
        end

        # push!(X, copy(model.sol))

        # check status
        if model.status == :Optimal
            if verbose == 1
                println()
                println("Optimal solution found.")
            end
            return model.status
        end

    end

    #
    return model.status

end

"""
    compute_starting_point!
    Compute a starting point

# Arguments
-`model::Model`
-`F::Factorization`: Cholesky factor of A*A', where A is the constraint matrix
    of the optimization problem.
"""
function compute_starting_point!(
    A::AbstractMatrix{Tv},
    F::Factorization{Tv},
    Λ::Tulip.PrimalDualPoint{Tv},
    b::StridedVector{Tv},
    c::StridedVector{Tv},
    uind::StridedVector{Ti},
    uval::StridedVector{Tv}
    ) where{Tv<:Real, Ti<:Integer}

    (m, n) = size(A)
    p = size(uind, 1)

    rhs = - 2 * b
    u_ = zeros(n)
    spxpay!(1.0, u_, uind, uval)
    rhs += A* u_

    #==============================#
    #
    #   I. Solve two QPs
    #
    #==============================#

    # Compute x0
    v = F \ rhs

    copy!(Λ.x, -0.5 * At_mul_B(A, v))
    spxpay!(0.5, Λ.x, uind, uval)

    # Compute w0
    @inbounds for i in 1:p
        j = uind[i]
        Λ.w[i] = uval[i] - Λ.x[j]
    end

    # Compute y0
    copy!(Λ.y, F \ (A*c))

    # Compute s0
    copy!(Λ.s, 0.5 * (At_mul_B(A, Λ.y) - c))

    # Compute z0
    @inbounds for i in 1:p
        j = uind[i]
        Λ.z[i] = - Λ.s[j]
    end

    #==============================#
    #
    #   II. Compute correction
    #
    #==============================#

    dp = zero(Tv)
    dd = zero(Tv)

    @inbounds for i in 1:n
        tmp = -1.5 * Λ.x[i]
        if tmp > dp
            dp = tmp
        end
    end

    @inbounds for i in 1:p
        tmp = -1.5 * Λ.w[i]
        if tmp > dp
            dp = tmp
        end
    end

    @inbounds for i in 1:n
        tmp = -1.5 * Λ.s[i]
        if tmp > dd
            dd = tmp
        end
    end

    @inbounds for i in 1:p
        tmp = -1.5 * Λ.z[i]
        if tmp > dd
            dd = tmp
        end
    end

    tmp = dot(Λ.x + dp, Λ.s + dd) + dot(Λ.w + dp, Λ.z + dd)

    dp += 0.5 * tmp / (sum(Λ.s + dd) + sum(Λ.z + dd))
    dd += 0.5 * tmp / (sum(Λ.x + dp) + sum(Λ.w + dp))

    #==============================#
    #
    #   III. Apply correction
    #
    #==============================#

    @inbounds for i in 1:n
        Λ.x[i] += dp
    end

    @inbounds for i in 1:n
        Λ.s[i] += dd
    end

    @inbounds for i in 1:p
        Λ.w[i] += dp
    end

    @inbounds for i in 1:p
        Λ.z[i] += dd
    end

    # Done
    return Λ
end

function compute_next_iterate!(model::Model, F::Factorization)
    (x, y, s, w, z) = (model.sol.x, model.sol.y, model.sol.s, model.sol.w, model.sol.z)
    (m, n, p) = model.nconstr, model.nvars, size(model.uind, 1)

    d_aff = copy(model.sol)
    d_cc = copy(model.sol)

    # compute residuals
    μ = (
        (dot(x, s) + dot(w, z))
        / (n + p)
    )

    rb = (model.A * x) - model.b
    rc = At_mul_B(model.A, y) + s - model.c
    spxpay!(-1.0, rc, model.uind, model.sol.z)

    ru = x[model.uind] + w - model.uval
    rxs = x .* s
    rwz = w .* z

    θ = x ./ s
    update_theta!(θ, x, s, z, w, model.uind)

    # compute predictor
    solve_newton!(
        model.A,
        θ,
        F,
        model.sol,
        d_aff,
        model.uind,
        -rb,
        -rc,
        -ru,
        -rxs,
        -rwz
    )

    # compute step length
    (α_pa, α_da) = compute_stepsize(model.sol, d_aff)

    # update centrality parameter
    μ_aff = (
        (
            dot(x + α_pa * d_aff.x, s + α_da * d_aff.s)
            + dot(w + α_pa * d_aff.w, z + α_da * d_aff.z)
        ) / (n + p)
    )

    σ = clamp((μ_aff / μ)^3, 10.0^-12, 1.0 - 10.0^-12)  # clamped for numerical stability
    # compute corrector
    solve_newton!(
        model.A,
        θ,
        F,
        model.sol,
        d_cc,
        model.uind,
        zeros(m),
        zeros(n),
        zeros(p),
        σ*μ*ones(n) - d_aff.x .* d_aff.s,
        σ*μ*ones(p) - d_aff.w .* d_aff.z
    )

    # final step size
    d = d_aff + d_cc
    (α_p, α_d) = compute_stepsize(model.sol, d, damp=0.99995)

    # take step
    model.sol.x += α_p * d.x
    model.sol.y += α_d * d.y
    model.sol.s += α_d * d.s
    model.sol.w += α_p * d.w
    model.sol.z += α_d * d.z

    return model.sol
end

"""
    symbolic_cholesky
    Compute Cholesky factorization of A*A'
"""
function symbolic_cholesky(A::AbstractMatrix{T}) where {T<:Real}

    F = Cholesky.cholesky(A, ones(A.n))
    return F

end

"""
    compute_newton!
    Form and factorize the Newton system, using the normal equations.
"""
function compute_newton!(
    A::AbstractMatrix{Ta},
    x::AbstractVector{Tx},
    s::AbstractVector{Ts},
    w::AbstractVector{Tw},
    z::AbstractVector{Tz},
    uind::AbstractVector{Ti},
    θ::AbstractVector{T},
    F::Factorization{Ta}
    ) where {Ta<:Real, Tx<:Real, Ts<:Real, Tw<:Real, Tz<:Real, Ti<:Integer, T<:Real}

    # Compute Θ = (X^{-1} S + W^{-1} Z)^{-1}
    θ = x ./ s
    for i in 1:size(uind, 1)
        j = uind[i]
        θ[j] = 1.0 / (s[j] / x[j] + z[i] / w[i])
    end

    # Form the normal equations matrix and compute its factorization
    Cholesky.cholesky!(A, θ, F)

    return θ
end

"""
    solve_newton
    Solve Newton system with the given right-hand side.
    Overwrites the input d
"""
function solve_newton!(
    A::AbstractMatrix{Ta},
    θ::AbstractVector{T1},
    F::Factorization{Ta},
    Λ::PrimalDualPoint,
    d::PrimalDualPoint,
    uind::AbstractVector{Ti},
    ξ_b::AbstractVector{T2},
    ξ_c::AbstractVector{T3},
    ξ_u::AbstractVector{T4},
    ξ_xs::AbstractVector{T5},
    ξ_wz::AbstractVector{T6},
) where {Ta<:Real, T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real, T6<:Real, Ti<:Integer}

    ξ_tmp = ξ_c - (ξ_xs ./ Λ.x)
    ξ_tmp[uind] += (ξ_wz - (Λ.z .* ξ_u)) ./ Λ.w

    d.y = F \ (ξ_b + A * (θ .* ξ_tmp))

    d.x = θ .* (At_mul_B(A, d.y) - ξ_tmp)

    d.z = (Λ.z .* (-ξ_u + d.x[uind]) + ξ_wz) ./ Λ.w
    d.s = (ξ_xs - Λ.s .* d.x) ./ Λ.x
    d.w = (ξ_wz - Λ.w .* d.z) ./ Λ.z

    # # check if system is solved correctly
    # rb = ξ_b - A*d.x
    # rc = ξ_c - (A'*d.y + d.s)
    # rc[uind] += d.z
    # ru = ξ_u - (d.x[uind] + d.w)
    # rxs = ξ_xs - (Λ.s .* d.x + Λ.x .* d.s)
    # rwz = ξ_wz - (Λ.z .* d.w + Λ.w .* d.z)

    # println("Residuals\t(normal eqs)")
    # println("||rb||   \t", @sprintf("%.6e", maximum(abs.(rb))))
    # println("||rc||   \t", @sprintf("%.6e", maximum(abs.(rc))))
    # println("||ru||   \t", @sprintf("%.6e", maximum(abs.(ru))))
    # println("||rxs||  \t", @sprintf("%.6e", maximum(abs.(rxs))))
    # println("||rwz||  \t", @sprintf("%.6e", maximum(abs.(rwz))))
    # println()
    return d
end

function compute_stepsize(
    tx::AbstractVector{T}, tw::AbstractVector{T}, ts::AbstractVector{T}, tz::AbstractVector{T},
    dx::AbstractVector{T}, dw::AbstractVector{T}, ds::AbstractVector{T}, dz::AbstractVector{T};
    damp=1.0
) where T<:Real

    n = size(tx, 1)
    p = size(tw, 1)
    n == size(ts, 1) || throw(DimensionMismatch("t.s is wrong size"))
    p == size(tz, 1) || throw(DimensionMismatch("t.z is wrong size"))

    n == size(dx, 1) || throw(DimensionMismatch("d.x is wrong size"))
    n == size(ds, 1) || throw(DimensionMismatch("d.s is wrong size"))
    p == size(dw, 1) || throw(DimensionMismatch("d.w is wrong size"))
    p == size(dz, 1) || throw(DimensionMismatch("d.z is wrong size"))

    ap, ad = -1.0, -1.0

    @inbounds for i in 1:n
        if dx[i] < 0.0
            if (tx[i] / dx[i]) > ap
                ap = (tx[i] / dx[i])
            end
        end
    end

    @inbounds for i in 1:n
        if ds[i] < 0.0
            if (ts[i] / ds[i]) > ad
                ad = (ts[i] / ds[i])
            end
        end
    end

    @inbounds for j in 1:p
        if dw[j] < 0.0
            if (tw[j] / dw[j]) > ap
                ap = (tw[j] / dw[j])
            end
        end
    end

    @inbounds for j in 1:p
        if dz[j] < 0.0
            if (tz[j] / dz[j]) > ad
                ad = (tz[j] / dz[j])
            end
        end
    end

    ap = - damp * ap
    ad = - damp * ad

    # println("\t(ap, ad) = ", (ap, ad))

    return (ap, ad)

end

function compute_stepsize(t::PrimalDualPoint{T}, d::PrimalDualPoint{T}; damp=1.0) where T<:Real
    (ap, ad) = compute_stepsize(t.x, t.w, t.s, t.z, d.x, d.w, d.s, d.z, damp=damp)
    return (ap, ad)
end

function update_theta!(θ, x, s, z, w, colind)
    # only called from within the optimization, so bounds were checked before
    for i in 1:size(colind, 1)
        j = colind[i]
        θ[j] = 1.0 / (s[j] / x[j] + z[i] / w[i])
    end
    return nothing
end

"""
    spxpay!(α, x, y_ind, y_val)

In-place computation of x += α * y, where y = sparse(y_ind, y_val)

# Arguments
"""
function spxpay!(α::Tv, x::AbstractVector{Tv}, y_ind::AbstractVector{Ti}, y_val::AbstractVector{Tv}) where{Ti<:Integer, Tv<:Real}
    for i in 1:size(y_ind, 1)
        j = y_ind[i]
        x[j] = x[j] + α * y_val[i]
    end
    return nothing
end


text:
using GLM, RDatasets, DataFrames, Distributions, Random, LinearAlgebra
Random.seed!(0)

df = dataset("MASS", "cpus")
n = size(df)[1]
df = df[shuffle(1:n),:]

pTest = 0.2
lastTindex = Int(floor(n*(1-pTest)))
numTest = n - lastTindex

train = df[1:lastTindex,:]
test = df[lastTindex+1:n,:]

form = @formula(Perf~CycT+MMin+MMax+Cach+ChMin+ChMax)
model1 = glm(form, train, Normal(),  IdentityLink())
model2 = glm(form, train, Poisson(), LogLink())
model3 = glm(form, train, Gamma(),  InverseLink())

invIdenityLink(x) = x
invLogLink(x) = exp(x)
invInverseLink(x) = 1/x

A = [ones(numTest) test.CycT test.MMin test.MMax test.Cach test.ChMin test.ChMax]
pred1 = invIdenityLink.(A*coef(model1))
pred2 = invLogLink.(A*coef(model2))
pred3 = invInverseLink.(A*coef(model3))

actual = test.Perf
lossModel1 = norm(pred1 - actual)
lossModel2 = norm(pred2 - actual)
lossModel3 = norm(pred3 - actual)

println("Model 1: ", coef(model1))
println("Model 2: ", coef(model2))
println("Model 3: ", coef(model3))
println("\nLoss of models 1,2,3: ",(lossModel1 ,lossModel2, lossModel3))

ext: jl

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# [[file:../fresa.org::testrosenbrock][testrosenbrock]]
using Fresa
nx = 2
x0 = 0.5*ones(nx)
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x)))
rosenbrock(x) = ([(1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2], 0)
# f = x -> ((x[1]-3)^2+(x[2]-5)^2+8, 0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,rosenbrock)]
# now invoke Fresa to solve the problem
best, pop = Fresa.solve(rosenbrock, p0, domain; ngen=1000, tolerance=1e-8)
println("Population at end: $pop")
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")
# testrosenbrock ends here


text:
lookup, _, lines... = readlines() .|> s -> (Int.(collect(s) .== '#') |> join)
default = '0'
for iteration in 1:50
    sy, sx = length(lines), length(lines[1])
    atmap(y, x) = (y ∈ 1:sy && x ∈ 1:sx) ? lines[y][x] : default
    get_index(Y, X) = parse(Int, join(atmap(Y+y, X+x) for x=-1:1, y=-1:1), base=2) + 1
    global lines = [[lookup[get_index(Y, X)] for X=0:sx+1] for Y=0:sy+1]
    global default = (default == '0') ? lookup[1] : lookup[end]
    iteration ∈ [2, 50] && println(lines .|> (line -> line .== '1') |> sum |> sum)
end

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
type PRacos
    rc::RacosCommon
    evaluation_server_num
    sample_set
    result_set
    asyn_result
    history
    is_finish

    function PRacos(ncomputer)
        new(RacosCommon(), ncomputer, RemoteChannel(()->Channel(ncomputer)),
        RemoteChannel(()->Channel(ncomputer)), RemoteChannel(()->Channel(1)),
        RemoteChannel(()->Channel(1)), false)
    end
end


text:
using Mosek

escapename(n::String) =
    replace(n,r"(^#)","_#")

fmtcof(c::Float64) =
    if     c < 0.0
        if     c < -1.0 || c > -1.0
             "- $(-c)"
        else
            "-"
        end
    elseif c > 0.0
        if c < 1.0 || c > 1.0
            "+ $c"
        else
            "+"
        end
    else
        "0"
    end

write(t::Mosek.Task, filename::String) =
    open(filename,"w") do f
        write(t,f)
    end

function write(t::Mosek.Task, f::IO)
    numvar    = getnumvar(t)
    numbarvar = getnumbarvar(t)
    numcon    = getnumcon(t)
    numcone   = getnumcone(t)
    numanz    = getnumanz(t)

    varnames = String[ escapename(getvarname(t,i)) for i in 1:numvar]
    barvarnames = String[ escapename(getbarvarname(t,i)) for i in 1:numbarvar]
    connames = String[ escapename(getconname(t,i)) for i in 1:numcon]
    conenames = String[ escapename(getconename(t,i)) for i in 1:numcone]

    for i in 1:numvar
        if length(varnames[i]) == 0
            varnames[i] = "#x$i"
        end
    end

    for i in 1:numcon
        if length(connames[i]) == 0
            connames[i] = "#c$i"
        end
    end

    for i in 1:numcone
        if length(conenames[i]) == 0
            conenames[i] = "#k$i"
        end
    end

    for i in 1:numbarvar
        if length(barvarnames[i]) == 0
            barvarnames[i] = "#X̂$i"
        end
    end

    maxvarnamelen = max(map(n -> length(n),varnames)...)
    maxconnamelen = max(map(n -> length(n),connames)...)

    taskname = gettaskname(t)

    println(f,"Task $(escape_string(taskname))")

    _,barcidx = getbarcsparsity(t)

    _,baraidx = getbarasparsity(t)
    baraptrb = zeros(Int,numcon+1)
    for idx in baraidx
        (barai,baraj) = getbaraidxij(t,idx)
        baraptrb[barai+1] += 1
    end
    for i in 1:numcon
        baraptrb[i+1] += baraptrb[i]
    end
    baraptrb += 1

    baraptre = copy(baraptrb)
    baraidx_idx = Vector{Int}(length(baraidx))
    for idx in baraidx
        (barai,baraj) = getbaraidxij(t,idx)
        baraidx_idx[baraptre[barai]] = idx
        baraptre[barai] += 1
    end

    # Objective
    objsense = getobjsense(t)
    objname  = getobjname(t)

    if objsense == MSK_OBJECTIVE_SENSE_MAXIMIZE
        println(f,"Maximize $(escape_string(objname))")
    else
        println(f,"Minimze $(escape_string(objname))")
    end
    c = getc(t)
    print(f,"    ")
    for j in 1:numvar
        if c[j] < 0.0 || c[j] > 0.0
            print(f,"$(fmtcof(c[j])) $(varnames[j]) ")
        end
    end

    # Bar entries go here

    cfix = getcfix(t)
    if cfix < 0 || cfix > 0
        print(f,fmtcof(cfix))
    end
    println()

    # Constraints
    if numcon+numcone > 0
        println(f,"Subject to")

        for i in 1:numcon
            (nzi,subi,vali) = getarow(t,Int32(i))
            (bk,bl,bu) = getconbound(t,Int32(i))

            print(f,"    $(connames[i]): ")
            if bk == MSK_BK_LO || bk == MSK_BK_RA
                print(f,"$bl < ")
            end
            for k in 1:nzi
                print(f,"$(fmtcof(vali[k])) $(varnames[subi[k]]) ")
            end

            for k in baraptrb[i]:baraptrb[i+1]-1
                idx = baraidx_idx[k]
                (barai,baraj,num,sub,w) = getbaraidx(t,idx)
                if num == 1
                    print(f,"$(fmtcof(w[1])) #MX$(sub[1]) ⋅ $(barvarnames[baraj]) ")
                elseif num > 1
                    print(f,"(")
                    printf(f,"$(fmtcof(w[1])) #MX$(sub[1])")
                    print(f,") ⋅ $(barvarnames[baraj]) ")
                end
            end

            if bk == MSK_BK_UP || bk == MSK_BK_RA
                println(f,"< $bu")
            elseif bk == MSK_BK_FX
                println(f,"= $bu")
            end

        end

        for k in 1:numcone
            (ct,conepar,nummem,submem) = getcone(t,k)
            dom =
                if ct == MSK_CT_QUAD "Q"
                elseif ct == MSK_CT_RQUAD "Q_r"
                else "?"
                end

            print(f,"    $(conenames[k]): ($(varnames[submem[1]])")
            for i in 2:nummem
                print(f,",$(varnames[submem[i]])")
            end
            println(f,") ∈ $dom(nummem)")
        end
    end

    # Variable bounds
    if numvar > 0
        println(f,"Variables")
        for j in 1:numvar
            bk,bl,bu = getvarbound(t,j)

            if bk == MSK_BK_FX
                println(f,"    $(varnames[j]) = $bl")
            else
                blstr =
                    if     bk == MSK_BK_LO || bk == MSK_BK_RA "[$bl"
                    else   "]-Inf"
                    end
                bustr =
                    if     bk == MSK_BK_UP || bk == MSK_BK_RA "$bu]"
                    else   "+Inf["
                    end

                println(f,"    $(varnames[j]) ∈ $blstr;$bustr")
            end
        end

        for j in 1:numbarvar
            dim = getdimbarvarj(t,j)
            println(f,"    $(barvarnames[j]) ∈ S($dim)")
        end

    end

    # Matsto
    numsymmat = getnumsymmat(t)
    if numsymmat > 0
        println(f,"Symmetric matrixes")
        for i in 1:numsymmat
            dim,nz,tp = getsymmatinfo(t,i)
            (msubi,msubj,mvalij) = getsparsesymmat(t,i)
            print(f,"    #MX$i ($dim): [ " )
            if nz > 0
                print(f,"($(msubi[1]),$(msubj[1]),$(mvalij[1]))")
                for i in 2:nz
                    print(f,",($(msubi[i]),$(msubj[i]),$(mvalij[i]))")
                end
            end
            println(f," ]")
        end
    end
end

maketask() do t
    readdata(t,ARGS[1])
    write(t,STDOUT)
end


text:
export cg
using LinearAlgebra

"""
`cg(A, b [, x0];  maxiter=10, tol=1e-3, log=false) -> x [, history_x, history_r]`

Solve the system `Ax = b` of a symmetric, positive definite linear map given by `A` using the conjugate gradients method.

Parameters
----------
* `A`

Linear map, which must have two methods defined: multiplication by `AbstractVector` and transpose multiplication by `AbstractVector`.

* `b` : `AbstractVector`

RHS of equation `Ax = b`.

* `x0` : `AbstractVector`

Starting solution.

* `maxiter` : Int, optional

Maximum number of iterations. Defaults to 10.

* `tol` : Real, optional

Stopping criterion tolerance. Stops if the square root of the ration between the previous and current residuals is smaller than `tol`. Defaults to 1e-3.

* `log` : Bool, optional

If `true`, returns also model and residuals at each iterations.
"""
function cg(A, b::AbstractVector, x0::AbstractVector; maxiter::Int=10, tol::Real=1e-3, log::Bool=false)
    x = x0[:]
    r = b - A * x
    p = r
    rsold = dot(r, r)
    r1 = rsold
    if log
        history_x = [x]
        history_r = [rsold]
    end

    for i = 1:minimum([length(b), maxiter])
        Ap = A * p
        alpha = rsold / dot(p, Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        rsnew = dot(r, r)
        p = r + (rsnew / rsold) * p
        if log
            push!(history_x, x)
            push!(history_r, rsnew)
        end
        if sqrt(rsnew/rsold) < tol || p ≈ zero(p)
            break
        end
        rsold = rsnew
    end
    if log
        return x, history_x, history_r
    else
        return x
    end
end
cg(A, b; kwargs...) = cg(A, b, zeros(length(b)); kwargs...)


text:
function corner(r)
    Flg = true
    while !isborder(r,Sud) || Flg
        if !isborder(r,Sud)
            move!(r,Sud)
        end
        if isborder(r,Sud)
            Flg =  ob(r,2)
        end
    end
    while !isborder(r,West)
        move!(r,West)
    end
    if !isborder(r,Nord)
        move!(r,Nord)
    else
        ob(r,0)
    end
end

function ob(r,Side)
    x = 0
    Flg = true
    while isborder(r,HorizonSide(Side)) && !isborder(r,HorizonSide((Side + 3) % 4))
        move!(r,HorizonSide((Side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 3) % 4))
        Flg = false
    end
    if  !isborder(r,HorizonSide(Side))
        move!(r,HorizonSide(Side))
    end
    while !isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 1) % 4))
        move!(r,HorizonSide(Side))
    end
    for i = 1:x
        move!(r,HorizonSide((Side+1) % 4))
    end
    return Flg
end

function counter(r)
    c1 = 0
    c2 = 0
    FlagCount = true
    while  !isborder(r,Ost) || FlagCount
        if isborder(r,Sud)
            c1 += 1
        else
            if c1 != 0
                c2 += 1
            end
            c1 = 0
        end
        if !isborder(r,Ost)
            move!(r,Ost)
        else
            FlagCount = false
            if c1 != 0
                c2 += 1
            end
        end
    end
    while !isborder(r,West)
        move!(r,West)
    end
    return(c2)
end

function task_20(r)
    corner(r)
    Flg2 = true
    sum = 0
    while !isborder(r,Nord)
        sum += counter(r)
        if !isborder(r,Nord)
            move!(r,Nord)
        end
    end
    sum += counter(r)
    return sum
end
task_20(r)


text:
"""
    determinant(mat)

Given a non singluar matrix, calculate its determinant using LU decomposition.

L and U are lower triangular and upper triangular matrices respectively such that

A = L*U

If we want to find the determinant, then

det(A) = det(LU) = det(L)*det(U)

Determinant of triangualar matrices is the product of their diagonal entries. Hence, makes finding the determinant easy.
"""
function determinant(mat)
    n, m = size(mat)
    if n != m
        DomainError(mat, "The matrix should be a square matrix.")
    end
    L, U = lu_decompose(mat)
    l_prod = prod([L[i, i] for i in 1:n])
    u_prod = prod([U[i, i] for i in 1:n])

    return l_prod * u_prod
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
type outOutput
    outid
    keepid
    outval
    nout
    method
end

type indtOutput
    method
    dstat
    pval_d
    wstat
    pval_w
    flag
end

type idealfOutput
    lower_quartile
    upper_quartile
    method
end

type testOutput
    method
    h0
    n
    df
    estimate
    se
    statistic
    crit
    ci
    p
end

type indirectTestOutput
    nboot
    n
    regfit
    sobel
    bootest
    bootci
    bootse
    p
end

indirectTestOutput()=indirectTestOutput(nothing,
                                        nothing,
                                        nothing,
                                        nothing,
                                        nothing,
                                        nothing,
                                        nothing,
                                        nothing)

testOutput()=testOutput(nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing,
                        nothing)


text:
"""
    ArrayOutput <: Output

    ArrayOutput(init; tspan::AbstractRange, [aux, mask, padval])

A simple output that stores each step of the simulation in a vector of arrays.

# Arguments

- `init`: initialisation `Array` or `NamedTuple` of `Array`

# Keywords

- `tspan`: `AbstractRange` timespan for the simulation
- `aux`: NamedTuple of arbitrary input data. Use `get(data, Aux(:key), I...)`
    to access from a `Rule` in a type-stable way.
- `mask`: `BitArray` for defining cells that will/will not be run.
- `padval`: padding value for grids with neighborhood rules. The default is `zero(eltype(init))`.
"""
mutable struct ArrayOutput{T,F<:AbstractVector{T},E} <: Output{T,F}
    frames::F
    running::Bool
    extent::E
end
function ArrayOutput(; frames, running, extent, kw...)
    append!(frames, _zerogrids(init(extent), length(tspan(extent))-1))
    ArrayOutput(frames, running, extent)
end

"""
    ResultOutput <: Output

    ResultOutput(init; tspan::AbstractRange, kw...)

A simple output that only stores the final result, not intermediate frames.

# Arguments

- `init`: initialisation `Array` or `NamedTuple` of `Array`

# Keywords

- `tspan`: `AbstractRange` timespan for the simulation
- `aux`: NamedTuple of arbitrary input data. Use `get(data, Aux(:key), I...)`
    to access from a `Rule` in a type-stable way.
- `mask`: `BitArray` for defining cells that will/will not be run.
- `padval`: padding value for grids with neighborhood rules. The default is `zero(eltype(init))`.
"""
mutable struct ResultOutput{T,F<:AbstractVector{T},E} <: Output{T,F}
    frames::F
    running::Bool
    extent::E
end
ResultOutput(; frames, running, extent, kw...) = ResultOutput(frames, running, extent)

isstored(o::ResultOutput) = false
storeframe!(o::ResultOutput, data::AbstractSimData) = nothing

function finalise!(o::ResultOutput, data::AbstractSimData)
    _storeframe!(eltype(o), o, data)
end


text:
"""
UnweightedParticleFilter

A particle filter that does not use any reweighting, but only keeps particles if the observation matches the true observation exactly. This does not require obs_weight, but it will not work well in real-world situations.
"""
struct UnweightedParticleFilter{M, RNG<:AbstractRNG} <: Updater
    model::M
    n::Int
    rng::RNG
end

function UnweightedParticleFilter(model, n::Integer; rng=Base.GLOBAL_RNG)
    return UnweightedParticleFilter(model, n, rng)
end

function update(up::UnweightedParticleFilter, b::ParticleCollection, a, o)
    new = Random.gentype(b)[]
    i = 1
    while i <= up.n
        s = particle(b, mod1(i, n_particles(b)))
        sp, o_gen = @gen(:sp, :o)(up.model, s, a, up.rng)
        if o_gen == o
            push!(new, sp)
        end
        i += 1
    end
    if isempty(new)
        warn("""
             Particle Depletion!

             The UnweightedParticleFilter generated no particles consistent with observation $o. Consider upgrading to a BootstrapFilter or a BasicParticleFilter or creating your own domain-specific updater.
             """
            )
    end
    return ParticleCollection(new)
end

function update(up::UnweightedParticleFilter, b, a, o)
    return update(up, initialize_belief(up, b), a, o)
end

function initialize_belief(up::UnweightedParticleFilter, b)
    return ParticleCollection(collect(rand(up.rng, b) for i in 1:up.n))
end


text:
using PeeEmTee
using Test
using DelimitedFiles

const chargedist_file = joinpath(@__DIR__, "data", "chargedist_example.txt")

@testset "make_qfunc()" begin
    data = readdlm(chargedist_file)
    params = [-0.0009971235435351372,
              0.01266241070136633,
              0.1999583888727496,
              0.0996819367396509,
              0.09768461062509437,
              499.5820279060795]
    qfunc = make_qfunc(pmtresp, data[:, 1], data[:, 2])
    @test qfunc(params) == 66.02414426599655
end

@testset "pre_fit()" begin
    data = readdlm(chargedist_file)
    chargedist = ChargeDist(data[:, 1], data[:, 2])
    prefit_target = PreFit(-0.0009439575776608711,
                           0.012781642902921997,
                           455.31645300021285,
                           0.20214663928517476,
                           0.10553835465661239,
                           46.63963250021425)
    prefit_actual = pre_fit(chargedist)
    @test isapprox([getfield(prefit_target, field) for field in fieldnames(PreFit)],
                   [getfield(prefit_actual, field) for field in fieldnames(PreFit)])
end

@testset "pmtresp_fit()" begin
    data = readdlm(chargedist_file)
    chargedist = ChargeDist(data[:, 1], data[:, 2])
    prffit_target = PMTRespFit(-0.0009971235435351222,
                               0.0126624107013663,
                               0.1999583888727531,
                               0.09968193673964637,
                               0.0976846106250966,
                               499.58202790610073)
    prefit = pre_fit(chargedist)
    prffit_actual = pmtresp_fit(chargedist, prefit)
    @test isapprox([getfield(prffit_target, field) for field in fieldnames(PMTRespFit)],
                   [getfield(prffit_actual, field) for field in fieldnames(PMTRespFit)])
    prfuapfit_target = PMTRespUapFit(-0.0010866013489182497,
                                     0.012573335066167545,
                                     0.20255116463024672,
                                     0.09761207500313138,
                                     0.09719044076647844,
                                     495.05266162357447,
                                     0.012985248350779485,
                                     0.021308630425982703,
                                     4.592234142351302)
    prfuapfit_actual = pmtresp_fit(chargedist, prefit; mod=:uap)
    @test isapprox([getfield(prfuapfit_target, field) for field in fieldnames(PMTRespUapFit)],
                   [getfield(prfuapfit_actual, field) for field in fieldnames(PMTRespUapFit)])
end


text:
"""
    vertex4(para::GenericPara,
        extK = [DiagTree.getK(para.totalLoopNum, 1), DiagTree.getK(para.totalLoopNum, 2), DiagTree.getK(para.totalLoopNum, 3)],
        chan::AbstractVector = [PHr, PHEr, PPr, Alli],
        subdiagram = false;
        level = 1, name = :none, resetuid = false,
        phi_toplevel = ParquetBlocks().phi, ppi_toplevel = ParquetBlocks().ppi, Γ4_toplevel = ParquetBlocks().Γ4)

    Generate 4-vertex diagrams using Parquet Algorithm

# Arguments
- `para`            : parameters. It should provide internalLoopNum, interactionTauNum, firstTauIdx
- `extK`            : basis of external loops as a vector [left in, left out, right in, right out].
- `chan`            : vector of channels of the current 4-vertex.
- `subdiagram`      : a sub-vertex or not
- `name`            : name of the vertex
- `level`           : level in the diagram tree
- `resetuid`        : restart uid count from 1
- `phi_toplevel`    : channels of left sub-vertex for the particle-hole and particle-hole-exchange of the bubble at level one.
- `ppi_toplevel`    : channels of left sub-vertex for the particle-particle bubble at level one
- `Γ4_toplevel`     : channels of right sub-vertex for all all bubbles at level one

# Output
- A DataFrame with fields :response, :type, :extT, :diagram, :hash
"""
function vertex4(para::GenericPara,
    extK=[DiagTree.getK(para.totalLoopNum, 1), DiagTree.getK(para.totalLoopNum, 2), DiagTree.getK(para.totalLoopNum, 3)],
    chan::AbstractVector=[PHr, PHEr, PPr, Alli], subdiagram=false;
    level=1, name=:none, resetuid=false,
    phi_toplevel=ParquetBlocks().phi, ppi_toplevel=ParquetBlocks().ppi, Γ4_toplevel=ParquetBlocks().Γ4,
    subchannel::Symbol=:All #:All, :W, :Lver3, :Rver3, :RPA
)

    for k in extK
        @assert length(k) >= para.totalLoopNum "expect dim of extK>=$(para.totalLoopNum), got $(length(k))"
    end
    extK = [k[1:para.totalLoopNum] for k in extK]
    legK = extK

    resetuid && uidreset()

    if (para.extra isa ParquetBlocks) == false
        para = reconstruct(para, extra=ParquetBlocks())
    end

    @assert para.extra isa ParquetBlocks
    @assert para.totalTauNum >= maxVer4TauIdx(para) "Increase totalTauNum!\n$para"
    @assert para.totalLoopNum >= maxVer4LoopIdx(para) "Increase totalLoopNum\n$para"

    phi, ppi = para.extra.phi, para.extra.ppi

    @assert (PHr in phi) == false "PHi vertex is particle-hole irreducible, so that PHr channel is not allowed in $phi"
    @assert (PPr in ppi) == false "PPi vertex is particle-particle irreducible, so that PPr channel is not allowed in $ppi"
    @assert (PHr in phi_toplevel) == false "PHi vertex is particle-hole irreducible, so that PHr channel is not allowed in $phi_toplevel"
    @assert (PPr in ppi_toplevel) == false "PPi vertex is particle-particle irreducible, so that PPr channel is not allowed in $ppi_toplevel"

    @assert length(legK[1]) == length(legK[2]) == length(legK[3]) == para.totalLoopNum

    KinL, KoutL, KinR = legK[1], legK[2], legK[3]
    KoutR = (length(legK) > 3) ? legK[4] : KinL + KinR - KoutL
    @assert KoutR ≈ KinL + KinR - KoutL
    legK = [KinL, KoutL, KinR, KoutR]

    loopNum = para.innerLoopNum
    @assert loopNum >= 0

    diags = Diagram{para.weightType}[]

    if loopNum == 0
        if subchannel == :W || subchannel == :RPA || subchannel == :LVer3 || subchannel == :RVer3
            permutation = [Di,]
        else
            permutation = [Di, Ex]
        end
        append!(diags, bareVer4(para, legK, permutation))
    else # loopNum>0
        for c in chan
            if c == Alli
                continue
            end

            partition = orderedPartition(loopNum - 1, 4, 0)

            for p in partition

                if c == PHr || c == PHEr || c == PPr
                    bub = bubble(para, legK, c, p, level, name, phi_toplevel, ppi_toplevel, Γ4_toplevel, 1.0, subchannel)
                    append!(diags, bub)
                    if (NoBubble in para.filter) && (loopNum == 1) && (c == PHr || c == PHEr)
                        #add bubble counter-diagram to remove the bubble
                        cbub = bubble(para, legK, c, p, level, Symbol("$(name)_counter"), phi_toplevel, ppi_toplevel, Γ4_toplevel, -1.0, :RPA)
                        append!(diags, cbub)
                    end
                    # println(bub)
                end
            end
        end
        # # TODO: add envolpe diagrams
    end
    @assert all(x -> x.id isa Ver4Id, diags) "not all id are Ver4Id! $diags"
    @assert all(x -> x.id.extK ≈ legK, diags) "not all extK are the same! $diags"

    # @assert isempty(diags) == false "got empty ver4! $chan with\n $para\n"
    if isempty(diags)
        return DataFrame(response=[], type=[], extT=[], diagram=[])
    end

    df = toDataFrame(diags, expand=true)
    # println(df[:, [:response, :type, :extT, :diagram]])
    groups = mergeby(df, [:response, :type, :extT], name=name,
        getid=g -> Ver4Id(para, g[1, :response], g[1, :type], k=legK, t=g[1, :extT]) #generate id from the dataframe
    )
    return groups
end

function bubble(para::GenericPara, legK, chan::TwoBodyChannel, partition::Vector{Int}, level::Int, name::Symbol,
    phi_toplevel, ppi_toplevel, Γ4_toplevel, extrafactor=1.0, subchannel=:All)

    diag = Diagram{para.weightType}[]

    TauNum = para.interactionTauNum # maximum tau number for each bare interaction
    oL, oG0, oR, oGx = partition[1], partition[2], partition[3], partition[4]
    if isValidG(para.filter, oG0) == false || isValidG(para.filter, oGx) == false
        return diag
    end

    #the first loop idx is the inner loop of the bubble!
    LoopIdx = para.firstLoopIdx
    idx, maxLoop = findFirstLoopIdx(partition, LoopIdx + 1)
    LfirstLoopIdx, G0firstLoopIdx, RfirstLoopIdx, GxfirstLoopIdx = idx
    @assert maxLoop == maxVer4LoopIdx(para)

    diagType = [Ver4Diag, GreenDiag, Ver4Diag, GreenDiag]
    idx, maxTau = findFirstTauIdx(partition, diagType, para.firstTauIdx, TauNum)
    LfirstTauIdx, G0firstTauIdx, RfirstTauIdx, GxfirstTauIdx = idx
    @assert maxTau == maxVer4TauIdx(para) "Partition $partition with tauNum configuration $idx. maxTau = $maxTau, yet $(maxTauIdx(para)) is expected!"

    lPara = reconstruct(para, diagType=Ver4Diag, innerLoopNum=oL, firstLoopIdx=LfirstLoopIdx, firstTauIdx=LfirstTauIdx)
    rPara = reconstruct(para, diagType=Ver4Diag, innerLoopNum=oR, firstLoopIdx=RfirstLoopIdx, firstTauIdx=RfirstTauIdx)
    gxPara = reconstruct(para, diagType=GreenDiag, innerLoopNum=oGx, firstLoopIdx=GxfirstLoopIdx, firstTauIdx=GxfirstTauIdx)
    g0Para = reconstruct(para, diagType=GreenDiag, innerLoopNum=oG0, firstLoopIdx=G0firstLoopIdx, firstTauIdx=G0firstTauIdx)

    phi, ppi, Γ4 = para.extra.phi, para.extra.ppi, para.extra.Γ4
    if chan == PHr || chan == PHEr
        Γi = (level == 1) ? phi_toplevel : phi
        Γf = (level == 1) ? Γ4_toplevel : Γ4
    elseif chan == PPr
        Γi = (level == 1) ? ppi_toplevel : ppi
        Γf = (level == 1) ? Γ4_toplevel : Γ4
    else
        error("chan $chan isn't implemented!")
    end

    LLegK, K, RLegK, Kx = legBasis(chan, legK, LoopIdx)
    # println(K, ", ", Kx)

    ls, rs = subChannel(subchannel)

    Lver = vertex4(lPara, LLegK, Γi, true; level=level + 1, name=:Γi, subchannel=ls)
    isempty(Lver) && return diag
    # println("Γf: ", Γf)
    Rver = vertex4(rPara, RLegK, Γf, true; level=level + 1, name=:Γf, subchannel=rs)
    isempty(Rver) && return diag

    for ldiag in Lver.diagram
        for rdiag in Rver.diagram
            extT, G0T, GxT = tauBasis(chan, ldiag.id.extT, rdiag.id.extT)
            g0 = green(g0Para, K, G0T, true, name=:G0)
            gx = green(gxPara, Kx, GxT, true, name=:Gx)
            @assert g0 isa Diagram && gx isa Diagram
            append!(diag, bubble2diag(para, chan, ldiag, rdiag, legK, g0, gx, extrafactor))
        end
    end
    return diag
end

function bubble2diag(para, chan, ldiag, rdiag, extK, g0, gx, extrafactor)
    lid, rid = ldiag.id, rdiag.id
    ln, rn = lid.response, rid.response
    lo, ro = lid.para.innerLoopNum, rid.para.innerLoopNum
    vtype = typeMap(lid.type, rid.type)

    extT, G0T, GxT = tauBasis(chan, lid.extT, rid.extT)
    Factor = factor(para, chan) * extrafactor
    spin(response) = (response == UpUp ? "↑↑" : "↑↓")

    diag = Diagram{para.weightType}[]

    function add(Lresponse::Response, Rresponse::Response, Vresponse::Response, factor=1.0)
        if ln == Lresponse && rn == Rresponse
            nodeName = Symbol("$(spin(Lresponse))x$(spin(Rresponse)) → $chan,")
            id = Ver4Id(para, Vresponse, vtype, k=extK, t=extT, chan=chan)
            push!(diag, Diagram(id, Prod(), [g0, gx, ldiag, rdiag], factor=factor * Factor, name=nodeName))
        end
    end

    if chan == PHr
        add(UpUp, UpUp, UpUp, 1.0)
        add(UpDown, UpDown, UpUp, 1.0)
        add(UpUp, UpDown, UpDown, 1.0)
        add(UpDown, UpUp, UpDown, 1.0)
    elseif chan == PHEr
        add(UpUp, UpUp, UpUp, 1.0)
        add(UpUp, UpUp, UpDown, 1.0)
        add(UpDown, UpDown, UpUp, 1.0)
        add(UpDown, UpDown, UpDown, 1.0)
        #! the sign here is from the spin symmetry, not from the fermionic statistics
        add(UpUp, UpDown, UpDown, -1.0)
        #! the sign here is from the spin symmetry, not from the fermionic statistics
        add(UpDown, UpUp, UpDown, -1.0)
    elseif chan == PPr
        add(UpUp, UpUp, UpUp, 1.0)
        #! the sign here is from the spin symmetry, not from the fermionic statistics
        add(UpDown, UpDown, UpDown, -2.0)
        add(UpUp, UpDown, UpDown, 1.0)
        add(UpDown, UpUp, UpDown, 1.0)
    else
        error("chan $chan isn't implemented!")
    end

    return diag
end

function bareVer4(para::GenericPara, legK, diex::Vector{Permutation}=[Di, Ex])
    # @assert para.diagType == Ver4Diag

    KinL, KoutL, KinR = legK[1], legK[2], legK[3]
    t0 = para.firstTauIdx

    nodes = Diagram{para.weightType}[]

    q = [KinL - KoutL, KinR - KoutL]

    if para.hasTau
        extT_ins = [(t0, t0, t0, t0), (t0, t0, t0, t0)]
        extT_dyn = [(t0, t0, t0 + 1, t0 + 1), (t0, t0 + 1, t0 + 1, t0)]
        innerT_ins = [(1, 1), (1, 1)]
        innerT_dyn = [(t0, t0 + 1), (t0, t0 + 1)]
    else
        extT_ins = [(1, 1, 1, 1), (1, 1, 1, 1)]
        extT_dyn = extT_ins
        innerT_ins = [(1, 1), (1, 1)]
        innerT_dyn = innerT_ins
    end

    function bare(response::Response, type::AnalyticProperty, _diex::Permutation, _innerT, _q, _factor=1.0)
        @assert _diex == Di || _diex == Ex

        # there is an overall sign coming from Taylor expansion of exp(-S) depsite the statistics
        if _diex == Di
            sign = -1.0
        elseif _diex == Ex
            sign = para.isFermi ? 1.0 : -1.0
        else
            error("not implemented!")
        end

        if notProper(para, _q) == false && _diex in diex
            #create new bare ver4 only if _diex is required in the diex table
            vid = BareInteractionId(para, response, type, k=_q, t=_innerT, permu=_diex)
            return Diagram(vid, factor=sign * _factor)
        else
            return nothing
        end
    end

    function addver4!(response::Response, type, _extT, vd, ve)
        id_di = Ver4Id(para, response, type, k=legK, t=_extT[DI])
        (isnothing(vd) == false) && push!(nodes, Diagram(id_di, Sum(), [vd,]))

        id_ex = Ver4Id(para, response, type, k=legK, t=_extT[EX])
        (isnothing(ve) == false) && push!(nodes, Diagram(id_ex, Sum(), [ve,]))
        # end
    end

    function addresponse!(response::Response, type, _extT, _innerT)
        if response == UpUp
            vd = bare(response, type, Di, _innerT[DI], q[DI])
            ve = bare(response, type, Ex, _innerT[EX], q[EX])
            addver4!(UpUp, type, _extT, vd, ve)
        elseif response == UpDown
            vd = bare(UpDown, type, Di, _innerT[DI], q[DI])
            ve = nothing
            addver4!(UpDown, type, _extT, vd, ve)
        elseif response == ChargeCharge
            # UpUp channel
            vuud = bare(ChargeCharge, type, Di, _innerT[DI], q[DI])
            vuue = bare(ChargeCharge, type, Ex, _innerT[EX], q[EX])
            addver4!(UpUp, type, _extT, vuud, vuue)

            # UpDown channel
            vupd = bare(ChargeCharge, type, Di, _innerT[DI], q[DI])
            vupe = nothing
            # UpDown, exchange channel doesn't exist for the charge-charge interaction
            addver4!(UpDown, type, _extT, vupd, vupe)
        elseif response == SpinSpin
            # see manual/interaction.md for more details

            # UpUp channel
            vuud = bare(SpinSpin, type, Di, _innerT[DI], q[DI])
            vuue = bare(SpinSpin, type, Ex, _innerT[EX], q[EX])
            addver4!(UpUp, type, _extT, vuud, vuue)

            # UpDown channel
            vupd = bare(SpinSpin, type, Di, _innerT[DI], q[DI], -1.0)
            vupe = bare(SpinSpin, type, Ex, _innerT[EX], q[EX], 2.0)
            addver4!(UpDown, type, _extT, vupd, vupe)
        else
            error("not implemented!")
        end
    end

    for interaction in para.interaction
        response = interaction.response
        typeVec = interaction.type

        if Instant ∈ typeVec && Dynamic ∉ typeVec
            addresponse!(response, Instant, extT_ins, innerT_ins)
        elseif Instant ∉ typeVec && Dynamic ∈ typeVec
            addresponse!(response, Dynamic, extT_dyn, innerT_dyn)
        elseif Instant ∈ typeVec && Dynamic ∈ typeVec
            #if hasTau, instant interaction has an additional fake tau variable, making it similar to the dynamic interaction
            addresponse!(response, Instant, extT_ins, innerT_dyn)
            addresponse!(response, Dynamic, extT_dyn, innerT_dyn)
        end

        if D_Instant ∈ typeVec && D_Dynamic ∉ typeVec
            addresponse!(response, D_Instant, extT_ins, innerT_ins)
        elseif D_Instant ∉ typeVec && D_Dynamic ∈ typeVec
            addresponse!(response, D_Dynamic, extT_dyn, innerT_dyn)
        elseif D_Instant ∈ typeVec && D_Dynamic ∈ typeVec
            #if hasTau, instant interaction has an additional fake tau variable, making it similar to the dynamic interaction
            addresponse!(response, D_Instant, extT_ins, innerT_dyn)
            addresponse!(response, D_Dynamic, extT_dyn, innerT_dyn)
        end
    end

    return nodes
end

######################### utility functions ############################
maxVer4TauIdx(para) = (para.innerLoopNum + 1) * para.interactionTauNum + para.firstTauIdx - 1
maxVer4LoopIdx(para) = para.firstLoopIdx + para.innerLoopNum - 1

function legBasis(chan::TwoBodyChannel, legK, loopIdx)
    KinL, KoutL, KinR, KoutR = legK[1], legK[2], legK[3], legK[4]
    K = zero(KinL)
    K[loopIdx] = 1
    if chan == PHr
        Kx = KoutL + K - KinL
        LLegK = [KinL, KoutL, Kx, K]
        RLegK = [K, Kx, KinR, KoutR]
    elseif chan == PHEr
        Kx = KoutR + K - KinL
        LLegK = [KinL, KoutR, Kx, K]
        RLegK = [K, Kx, KinR, KoutL]
    elseif chan == PPr
        Kx = KinL + KinR - K
        LLegK = [KinL, Kx, KinR, K]
        RLegK = [K, KoutL, Kx, KoutR]
    else
        error("not implemented!")
    end

    # check conservation and momentum assignment
    @assert LLegK[INL] ≈ KinL
    @assert LLegK[INL] + LLegK[INR] ≈ LLegK[OUTL] + LLegK[OUTR]
    @assert RLegK[INL] + RLegK[INR] ≈ RLegK[OUTL] + RLegK[OUTR]

    return LLegK, K, RLegK, Kx
end

function tauBasis(chan::TwoBodyChannel, LvT, RvT)
    G0T = (LvT[OUTR], RvT[INL])
    if chan == PHr
        extT = (LvT[INL], LvT[OUTL], RvT[INR], RvT[OUTR])
        GxT = (RvT[OUTL], LvT[INR])
    elseif chan == PHEr
        extT = (LvT[INL], RvT[OUTR], RvT[INR], LvT[OUTL])
        GxT = (RvT[OUTL], LvT[INR])
    elseif chan == PPr
        extT = (LvT[INL], RvT[OUTL], LvT[INR], RvT[OUTR])
        GxT = (LvT[OUTL], RvT[INR])
    else
        error("not implemented!")
    end

    # make sure all tidx are used once and only once
    t1 = sort(vcat(collect(G0T), collect(GxT), collect(extT)))
    t2 = sort(vcat(collect(LvT), collect(RvT)))
    @assert t1 == t2 "chan $(chan): G0=$G0T, Gx=$GxT, external=$extT don't match with Lver4 $LvT and Rver4 $RvT"
    @assert extT[INL] == LvT[INL]
    return extT, G0T, GxT
end

function factor(para, chan)
    Factor = SymFactor[Int(chan)] / (2π)^para.loopDim
    if para.isFermi == false
        Factor = abs(Factor)
    end
    return Factor
end

function typeMap(ltype, rtype)
    if (ltype == Instant || ltype == Dynamic) && (rtype == Instant || rtype == Dynamic)
        return Dynamic
    elseif (ltype == D_Instant || ltype == D_Dynamic) && (rtype == Instant || rtype == Dynamic)
        return D_Dynamic
    elseif (ltype == Instant || ltype == Dynamic) && (rtype == D_Instant || rtype == D_Dynamic)
        return D_Dynamic
    else
        return nothing
    end
end

function subChannel(subchan)
    if subchan == :RPA
        return :RPA, :RPA
    elseif subchan == :LVer3
        return :LVer3, :All
    elseif subchan == :RVer3
        return :All, :RVer3
    elseif subchan == :W
        return :LVer3, :RVer3
    elseif subchan == :All
        return :All, :All
    else
        error("not implemented!")
    end
end


text:
function nextindex(k,m,A)
    return (k+m) % length(A)
end

function sdvig(A,m)
    for i in 1:length(A)
        A[i], A[nextindex(i,m,A)] = A[nextindex(i,m,A)], A[i]
    end
    return A
end


text:
using GraphNeuralNetworks
using Statistics:mean
using NNlib:gather

struct GlobalConcatenation{F} <: GNNLayer
  aggr::F
end

function (l::GlobalConcatenation)(g::GNNGraph, x::AbstractArray)
  d = reduce_nodes(l.aggr, g, x)
  d = gather(d, graph_indicator(g))
  return vcat(x, d)
end

(l::GlobalConcatenation)(g::GNNGraph) = GNNGraph(g, gdata=l(g, node_features(g)))

"""
    GinHP

Hyperparameters for the gin architecture.

"""
@kwdef struct GinHP
  depth_common :: Int = 12
  depth_phead :: Int = 3
  depth_vhead :: Int = 3
  hidden_size :: Int = 64
end

"""
    Gin <: TwoHeadGraphNeuralNetwork

A simple two-headed GNN architecture with only dense layers.
"""
mutable struct Gin <: TwoHeadGraphNeuralNetwork
  gspec
  hyper
  common
  vhead
  phead
end

function Gin(gspec::AbstractGameSpec, hyper::GinHP)
  Dense_layers(size, depth) = [Dense(size, size, relu) for _ in 1:depth]
  f(in, out) = Chain(Dense(in, out, relu),
                    Dense_layers(out, 2)...,
                    BatchNorm(out, relu))
  GIN_layers(in, out, depth) = [GNNChain(GINConv(f(i==1 ? in : out, out), 0)) for i in 1:depth]

  indim, _ = GI.state_dim(gspec)
  common = GNNChain(GIN_layers(indim, hyper.hidden_size, 3)...)
  vhead = GNNChain(GlobalPool(mean),
                   Dense_layers(hyper.hidden_size, 2)...,
                   Dense(hyper.hidden_size, 1, σ))
  phead = GNNChain(GlobalConcatenation(mean),
                   GIN_layers(2*hyper.hidden_size, 32, 2)...,
                   Dense(32, 1))
  return Gin(gspec, hyper, common, vhead, phead)
end

Network.HyperParams(::Type{Gin}) = GinHP

function Base.copy(nn::Gin)
  return Gin(
    nn.gspec,
    nn.hyper,
    deepcopy(nn.common),
    deepcopy(nn.vhead),
    deepcopy(nn.phead)
  )
end


text:
# Autogenerated wrapper script for Deno_jll for aarch64-apple-darwin
export deno

JLLWrappers.@generate_wrapper_header("Deno")
JLLWrappers.@declare_executable_product(deno)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        deno,
        "bin/deno",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
Learning repo for testing various CI and benchmarking tools.

Mainly benchmarking my naive ForwardDiff implementation and Dormand-Prince solver on models I find interesting.

Working through `Mathematical Physiology` by James Keener and `Biomolecular Feedback Systems` by Murray.
"""
module DemoActions

using StaticArrays

include("mysolvers.jl")
include("models.jl")
include("fwd.jl")

# solver interfaces
export MyProbInplace, mysolve, mysolve_inplace

# solver algorithms
export eulers_inplace_f!,
    midpoint!,
    midpoint_inplace_f!,
    rk4!,
    rk4_inplace_f!,
    rk4_inplace_f_prealloc_k!,
    dp5,
    dp5!,
    dp5_inplace!

# mini model-zoo
export lotka!, transcriptional_component!, bistable_gene_circuit!, negative_autoregulation!, phosphorylation_cycle!

end


text:
module ACTRSimulators
    using Reexport
    @reexport using DiscreteEventsLite, ACTRModels
    using ACTRModels, DiscreteEventsLite, DataStructures, CodeTracking
    import DiscreteEventsLite: run!, last_event!, is_running, print_event, register!
    import DiscreteEventsLite: Now, At, After, Every, AbstractScheduler
    import ACTRModels: AbstractACTR, utility_match
    export run!, next_event!, register!, vo_to_chunk, add_to_visicon!, clear_buffer!, add_to_buffer!, get_time, attending!
    export attend!, retrieving!, retrieve!, responding!, respond!
    export encoding!, encode!, all_match, AbstractTask, ACTRScheduler, setup_window, import_gui
    export remove_visual_object!, clear_visicon!, move_vo!
    export repaint!, draw_object!, clear!, draw_attention!, why_not

    include("simulator.jl")
    include("buffers.jl")
    include("declarative.jl")
    include("imaginal.jl")
    include("visual.jl")
    include("motor.jl")
    include("procedural.jl")
    include("utilities.jl")
end


text:
#=

# 3: Water collapse (implicit)

```@raw html
	<img src='../assets/collapse_exp.png' width="50%" height="50%" alt='missing' /><br>
```

Simulation of a water column collapsing under its own weight onto dry bottom.
Here with strictly incompressible approach (Projection method).
=#

module collapse_dry_implicit

using Printf
using SmoothedParticles
using LinearAlgebra
using IterativeSolvers
using IncompleteLU
using Parameters

#=
Declare constant parameters
=#

##kernel functions
const kernel = spline23
const Dkernel = Dspline23
const rDkernel = rDspline23

const dr = 2.0e-2         #average particle distance (decrease to make finer simulation)
const h = 2.8*dr          #size of kernel support
const rho0 = 1000.0       #fluid density
const g = -9.8*VECY       #gravitational acceleration
const mu = 8.4e-4     #dynamic viscosity
const m = dr^2*rho0       #particle mass

##geometry parameters
const water_column_width = 1.0
const water_column_height = 2.0
const box_height = 3.0
const box_width = 4.0
const nlayers = 3.5 #number of wall layers
const wall_width = nlayers*dr

##temporal parameters
const dt = h/40.0
const t_end = 1.0
const dt_frame = t_end/50

##labels for particle types
const FLUID = 0.
const  WALL = 1.
const DUMMY = 2.

#=
Declare variables to be stored in a Particle
=#

@with_kw mutable struct Particle <: AbstractParticle
	x::RealVector = VEC0 #position
	v::RealVector = VEC0 #velocity
	a::RealVector = VEC0 #acceleration
	P::Float64 = 0. #pressure
	div::Float64 = 0. #divergence of velocity
	L::Float64 = 0. #free surface identier
	L0::Float64 = 0.
	type::Float64 #particle type
end

#=
Define geometry and create particles
=#
function make_system()
	grid = Grid(dr, :hexagonal)
	box = Rectangle(0., 0., box_width, box_height)
	fluid = Rectangle(0., 0., water_column_width, water_column_height)
	walls = BoundaryLayer(box, grid, 1.2*dr)
	walls = Specification(walls, x -> (x[2] < box_height))
	dummy = BoundaryLayer(box, grid, nlayers*dr) - walls
	dummy = Specification(dummy, x -> (x[2] < box_height))
	domain = Rectangle(-nlayers*dr, -nlayers*dr, 2*box_width, 3*box_height)
	sys = ParticleSystem(Particle, domain, h)
	generate_particles!(sys, grid, fluid, x -> Particle(x=x, type=FLUID))
	generate_particles!(sys, grid, walls, x -> Particle(x=x, type=WALL))
	generate_particles!(sys, grid, dummy, x -> Particle(x=x, type=DUMMY))
	create_cell_list!(sys)
	apply!(sys, find_L0!)
	Lmax = maximum(p -> p.L0, sys.particles)
	ParticleField(sys, :L0) .= Lmax
	return sys
end

#=
Particle interactions
=#

function initialize!(p::Particle)
	if p.type == FLUID
		p.x += dt*p.v
		p.v += dt*g
	end
	p.div = 0.
	p.L = 0.
end

@inbounds function viscous_force!(p::Particle, q::Particle, r::Float64)
	p.a += 2.0*m*mu*rDkernel(h,r)/rho0^2*(p.v - q.v)
end

@inbounds function find_div_and_L!(p::Particle, q::Particle, r::Float64)
	p.div += -SmoothedParticles.dot(p.x - q.x, p.v - q.v)*m*rDkernel(h,r)/rho0
	p.L += -2.0*m*rDkernel(h,r)/rho0^2
end

@inbounds function find_L0!(p::Particle, q::Particle, r::Float64)
	p.L0 += -2.0*m*rDkernel(h,r)/rho0^2
end

@inbounds function internal_force!(p::Particle, q::Particle, r::Float64)
	p.a -= m*rDkernel(h,r)*(p.P + q.P)/rho0^2*(p.x - q.x)
end

function accelerate!(p::Particle)
	if p.type == FLUID
		p.v += dt*p.a
	end
	p.a = VEC0
end

#=
Functions to build the linear system
=#

function minus_laplace(p::Particle, q::Particle, r::Float64)::Float64
	if p == q
		return p.type == FLUID ? max(p.L, p.L0) : p.L
	end
	return 2.0*m*rDkernel(h,r)/rho0^2
end

function rhs(p::Particle)::Float64
	return -p.div/dt
end

#=
Time iteration
=#
function main()
	sys = make_system()
	out_pvd = new_pvd_file("results/collapse_dry_implicit")
	out_txt = open("results/collapse_dry_implicit/data.txt", "w")
	P = ParticleField(sys, :P)
	for k = 0 : Int64(round(t_end/dt))
		if (k %  Int64(round(dt_frame/dt)) == 0)
			@printf("t = %.6e\n", k*dt)
			println("# of particles = ", length(sys.particles))
			save_frame!(out_pvd, sys, :v, :P, :type)
			dimless_time = string(k*dt*sqrt(-g[2]/water_column_height))
			leading_edge = maximum(p -> (p.type == FLUID ? p.x[1] - water_column_width : 0.), sys.particles)/water_column_height
			write(out_txt, string(dimless_time)*" "*string(leading_edge)*"\n")
		end
		apply!(sys, initialize!)
		create_cell_list!(sys)
		apply!(sys, viscous_force!)

		##assemble linear system and solve for pressure
		apply!(sys, find_div_and_L!)
		A = assemble_matrix(sys, minus_laplace)
		b = assemble_vector(sys, rhs)
		try
			P .= cg(A, b; Pl = ilu(A; τ = 5.0))
		catch
			save_pvd_file(out_pvd)
			close(out_txt)
			error("Unable to solve linear system.")
		end
		apply!(sys, internal_force!)
		apply!(sys, accelerate!)
	end
	save_pvd_file(out_pvd)
	close(out_txt)
end

end


text:
#=
Track chunks begin with four bytes spelling out "MTrk", followed by the length
in bytes of the track (see readvariablelength in util.jl), followed by a sequence
of events.
=#
type MIDITrack
    events::Array{TrackEvent, 1}

    MIDITrack() = new(TrackEvent[])
    MIDITrack(events) = new(events)
end

function readtrack(f::IO)
    mtrk = join(map(Char, read(f, UInt8, 4)))
    if mtrk != MTRK
        error("Not a valid MIDI file. Expected MTrk, got $(mtrk) starting at byte $(hex(position(f)-4, 2))")
    end
    track = MIDITrack()

    # Get the length in bytes of the track
    tracklength = ntoh(read(f, UInt32))

    trackstart = position(f)
    bytesread = 0
    laststatus = UInt8(0) # Keeps track of running status
    while bytesread < tracklength
        # A track is made up of events. All events start with a variable length
        # value indicating the number of ticks (time) since the last event.

        # Get the time offset of the next event
        dT = readvariablelength(f)

        # Figure out the event type
        # Remember, endianness is byte order, not bit order. No need to ntoh here.
        event_start = read(f, UInt8)
        skip(f, -1)

        local event
        if isMIDIevent(event_start)
            event = readMIDIevent(dT, f, laststatus)
            laststatus = event.status
        elseif issysexevent(event_start)
            event = readsysexevent(dT, f)
            laststatus = UInt8(0)
        elseif ismetaevent(event_start)
            event = readmetaevent(dT, f)
            laststatus = UInt8(0)
        else
            error("Unrecognized event $(hex(event_start,2))")
        end
        push!(track.events, event)

        bytesread = position(f) - trackstart
    end

    # Validate that the track ends with a track end event
    lastevent = track.events[length(track.events)]
    if !isa(lastevent, MetaEvent) || lastevent.metatype != METATRACKEND
        error("Invalid track - does not end with track metaevent")
    else
        # strip the track end event - we don't need to worry about manipulating it
        track.events = track.events[1:length(track.events)-1]
    end

    track
end

function writetrack(f::IO, track::MIDITrack)
    write(f, convert(Array{UInt8, 1}, MTRK)) # Track identifier

    writingMIDI = false
    previous_status = UInt8(0)

    event_buffer = IOBuffer()

    for event in track.events
        if isa(event, MIDIEvent) && previous_status != 0 && previous_status == event.status
            writeevent(event_buffer, event, false)
        elseif isa(event, MIDIEvent)
            writeevent(event_buffer, event)
            previous_status = event.status
        else
            writeevent(event_buffer, event)
            previous_status = UInt8(0)
        end
    end

    # Write the track end event
    writeevent(event_buffer, MetaEvent(0, METATRACKEND, UInt8[]))

    bytes = takebuf_array(event_buffer)

    write(f, hton(UInt32(length(bytes))))
    write(f, bytes)
end

# Adds an event to a track, with an absolute time
function addevent(track::MIDITrack, time::Integer, newevent::TrackEvent)
    tracktime = 0
    addedevent = false
    for (i, event) in enumerate(track.events)
        if tracktime + event.dT > time
            # Add to track at position
            newdt = time - tracktime
            newevent.dT = newdt

            insert!(track.events, i, newevent)
            addedevent = true

            nextevent = track.events[i+1]
            nextevent.dT -= newdt

            break
        else
            tracktime += event.dT
        end
    end

    if !addedevent
        newdt = time - tracktime
        newevent.dT = newdt
        push!(track.events, newevent)
    end
end

# Adds a note to a track
function addnote(track::MIDITrack, note::Note)
    for (status, position) in [(NOTEON, note.position), (NOTEOFF, note.position + note.duration)]
        addevent(track, position, MIDIEvent(0, status | note.channel, UInt8[note.value, note.velocity]))
    end
end

# Adds a series of notes to a track
function addnotes(track::MIDITrack, notes::Array{Note, 1})
    for note in notes
        addnote(track, note)
    end
end

# Gets all of the notes on a track
function getnotes(track::MIDITrack)
    # Read through events until a noteon is found
    notes = Note[]
    tracktime = UInt64(0)
    for (i, event) in enumerate(track.events)
        tracktime += event.dT
        if isa(event, MIDIEvent) && event.status & 0xF0 == NOTEON
            duration = UInt64(0)
            for event2 in track.events[i+1:length(track.events)]
                duration += event2.dT
                # If we have a MIDI event & it's a noteoff, and it's for the same note as the first event we found, make a note
                if isa(event2, MIDIEvent) && event2.status & 0xF0 == NOTEOFF && event.data[1] == event2.data[1]
                    push!(notes, Note(event.data[1], duration, tracktime, event.status & 0x0F, event.data[2]))
                    break
                end
            end
        end
    end
    sort!(notes, lt=((x, y)->x.position<y.position))
end

# Change the program (instrument) on the given channel. Time is absolute, not relative to the last event.
function programchange(track::MIDITrack, time::Integer, channel::UInt8, program::UInt8)
    program = program - 1 # Program changes are typically given in range 1-128, but represented internally as 0-127.
    addevent(track, time, MIDIEvent(0, PROGRAMCHANGE | channel, UInt8[program]))
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function triangle(points, degree)
    sethue(cols[degree])
    poly(points, :fill)
end

function sierpinski(points, degree)
    triangle(points, degree)
    if degree > 1
        p1, p2, p3 = points
        sierpinski([p1, midpoint(p1, p2),
                        midpoint(p1, p3)], degree-1)
        sierpinski([p2, midpoint(p1, p2),
                        midpoint(p2, p3)], degree-1)
        sierpinski([p3, midpoint(p3, p2),
                        midpoint(p1, p3)], degree-1)
    end
end

function draw(fname, n)
  Drawing(200, 200, fname)
  origin()
  background("ivory")
  circle(O, 75, :clip)
  my_points = ngon(O, 150, 3, -pi/2, vertices=true)
  depth = 8 # 12 is ok, 20 is right out
  sierpinski(my_points, n)
  @test finish() == true
end

cols = distinguishable_colors(8)
fname = "sierpinski.pdf"
draw(fname, 8)
println("...finished test: output in $(fname)")


text:
# This file is a part of UpROOT.jl, licensed under the MIT License (MIT).

const kByteCountMask = UInt32(0x40000000)
const kByteCountMask_UInt16 = UInt16(kByteCountMask >> 16)

struct ROOTIOBuffer{T}
    buffer::Base.GenericIOBuffer{T}
end

ROOTIOBuffer(data::AbstractVector{UInt8}) = ROOTIOBuffer(IOBuffer(data))

Base.eof(io::ROOTIOBuffer) = eof(io.buffer)

function Base.read(io::ROOTIOBuffer, ::Type{String})
    short_len = read(io.buffer, UInt8)
    len = (short_len < 0xff) ? UInt32(short_len) : read(io.buffer, UInt32)
    String(read(io.buffer, len))
end

function Base.read(io::ROOTIOBuffer, ::Type{T}) where {T<:Real}
    x::T = read(io.buffer, T)
    if sizeof(T) > 1
        x = ntoh(x)
    end
    x
end

function Base.read(io::ROOTIOBuffer, ::Type{Vector{T}}) where {T<:Real}
    len = read(io, UInt32)
    A = Vector{T}(undef, len)
    read!(io.buffer, A)
    if sizeof(T) > 1
        A .= ntoh.(A)
    end
    A
end

struct ROOTClassHeader
    classname::String
end

function Base.read(io::ROOTIOBuffer, ::Type{ROOTClassHeader})
    classname = read(io, String)
    padding = read(io, UInt8)
    ROOTClassHeader(classname)
end

struct ROOTClassVersion
    nbytes::UInt32
    version::Int16
end

function Base.read(io::ROOTIOBuffer, ::Type{ROOTClassVersion})
    a = read(io, UInt16)

    if a & kByteCountMask_UInt16 > 0
        # byte count available
        b = read(io, UInt32)
        cnt = (UInt32(a & ~kByteCountMask_UInt16) << 16 | b >> 16)
        version = b % Int16
        bytecount = cnt - sizeof(version)
        ROOTClassVersion(bytecount, version)
    else
        ROOTClassVersion(0, signed(a))
    end
end

struct TObjectContent
    fUniqueID::UInt32
    fBits::UInt32
end

function Base.read(io::ROOTIOBuffer, ::Type{TObjectContent})
    ver = read(io, ROOTClassVersion)
    @assert ver.nbytes == 0
    @assert ver.version == 1
    fUniqueID = read(io, UInt32)
    fBits = read(io, UInt32)
    TObjectContent(fUniqueID, fBits)
end

struct TNamedContent{S<:AbstractString}
    tobject::TObjectContent
    fName::S
    fTitle::S
end

function Base.read(io::ROOTIOBuffer, ::Type{TNamedContent})
    ver = read(io, ROOTClassVersion)
    @assert ver.version == 1
    tobject = read(io, TObjectContent)
    fName = read(io, String)
    fTitle = read(io, String)
    TNamedContent(tobject, fName, fTitle)
end

struct TClonesArrayHeader{S<:AbstractString}
    tobject::TObjectContent
    fName::S
    fClass::S
    nobjects::Int
    fLowerBound::Int
end

function Base.read(io::ROOTIOBuffer, ::Type{TClonesArrayHeader})
    ver = read(io, ROOTClassVersion)
    @assert ver.version == 4
    tobject = read(io, TObjectContent)
    fName = read(io, String)
    fClass = read(io, String)
    nobjects = read(io, Int32)
    fLowerBound = read(io, Int32)

    # assuming CanBypassStreamer() is false:
    nch = read(io, Int8)
    TClonesArrayHeader(tobject, fName, fClass, Int(nobjects), Int(fLowerBound))
end


text:
using Pkg
Pkg.activate(".")

using LPSjl.BayesianNetworks
using LPSjl.ParameterLearning
using SparseArrays
using DataFrames

"""
    Neapolitan, R. Learning Bayesian Networks. Ch. 1, Pg. 36
    Figure 1.5
"""
node_C = DiscreteNode("C", ["c1", "c2"])
set_parameters!(node_C, sparse([9/13 4/13]))

node_V = DiscreteNode("V", ["v1", "v2"])
add_parent!(node_V, node_C)
set_parameters!(node_V, sparse([[1/3 2/3]; [1/2 1/2]]))

node_S = DiscreteNode("S", ["s1", "s2"])
add_parent!(node_S, node_C)
set_parameters!(node_S, sparse([[2/3 1/3]; [1/2 1/2]]))


text:
export ParameterServer, update!

using Flux: Params

struct ParameterServer
    rwl::ReaderCountRWLock
    params::Params
    ParameterServer(ps) = new(ReaderCountRWLock(), ps)
end

function update!(f, ps::ParameterServer)
    lock(ps.rwl) do
        f(ps.params)
    end
end

function Base.get(f, ps::ParameterServer)
    read_lock(ps.rwl) do
        f(ps.params)
    end
end


text:
module aoc

include("day.jl")
include("day1.jl")
include("day2.jl")
include("day3.jl")
include("day4.jl")
include("day5.jl")
include("day6.jl")
include("day7.jl")

include("day10.jl")
aoc.
end


text:
module GPAR

using Random, Stheno

include("util.jl")
include("gpar.jl")

end # module


text:
using Unitful

const NPS5 = [0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 30]
const S5i = [18, 23.4, 30.1, 38.9, 45, 57, 68.78, 84.68, 97.38, 110.08, 135.76, 162.76, 213.56, 266.2, 315.88, 347.68, 398.02, 448.62, 498.44, 549.44, 598.92, 749.3]u"mm"
const S5o = [21.3, 26.7, 33.4, 42.2, 48.3, 60.3, 73, 88.9, 101.6, 114.3, 141.3, 168.3, 219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610, 762]u"mm"
const S5t = [1.65, 1.65, 1.65, 1.65, 1.65, 1.65, 2.11, 2.11, 2.11, 2.11, 2.77, 2.77, 2.77, 3.4, 3.96, 3.96, 4.19, 4.19, 4.78, 4.78, 5.54, 6.35]u"mm"

const NPS10 = [0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36]
const S10i = [7.82, 10.4, 13.8, 17.08, 22.48, 27.86, 36.66, 42.76, 54.76, 66.9, 82.8, 95.5, 108.2, 134.5, 161.5, 211.58, 264.62, 314.66, 342.9, 393.7, 444.3, 495.3, 546.3, 597.3, 644.16, 695.16, 746.16, 797.16, 848.16, 898.16]u"mm"
const S10o = [10.3, 13.7, 17.1, 21.3, 26.7, 33.4, 42.2, 48.3, 60.3, 73, 88.9, 101.6, 114.3, 141.3, 168.3, 219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610, 660, 711, 762, 813, 864, 914]u"mm"
const S10t = [1.24, 1.65, 1.65, 2.11, 2.11, 2.77, 2.77, 2.77, 2.77, 3.05, 3.05, 3.05, 3.05, 3.4, 3.4, 3.76, 4.19, 4.57, 6.35, 6.35, 6.35, 6.35, 6.35, 6.35, 7.92, 7.92, 7.92, 7.92, 7.92, 7.92]u"mm"

const NPS20 = [8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36]
const S20i = [206.4, 260.3, 311.1, 339.76, 390.56, 441.16, 488.94, 539.94, 590.94, 634.6, 685.6, 736.6, 787.6, 838.6, 888.6]u"mm"
const S20o = [219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610, 660, 711, 762, 813, 864, 914]u"mm"
const S20t = [6.35, 6.35, 6.35, 7.92, 7.92, 7.92, 9.53, 9.53, 9.53, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7]u"mm"

const NPS30 = [0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 8, 10, 12, 14, 16, 18, 20, 22, 24, 28, 30, 32, 34, 36]
const S30i = [7.4, 10, 13.4, 16.48, 21.88, 27.6, 36.26, 41.94, 53.94, 63.44, 79.34, 92.04, 104.74, 205.02, 257.4, 307.04, 336.54, 387.34, 434.74, 482.6, 533.6, 581.46, 679.24, 730.24, 781.24, 832.24, 882.24]u"mm"
const S30o = [10.3, 13.7, 17.1, 21.3, 26.7, 33.4, 42.2, 48.3, 60.3, 73, 88.9, 101.6, 114.3, 219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610, 711, 762, 813, 864, 914]u"mm"
const S30t = [1.45, 1.85, 1.85, 2.41, 2.41, 2.9, 2.97, 3.18, 3.18, 4.78, 4.78, 4.78, 4.78, 7.04, 7.8, 8.38, 9.53, 9.53, 11.13, 12.7, 12.7, 14.27, 15.88, 15.88, 15.88, 15.88, 15.88]u"mm"

const NPS40 = [0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 24, 32, 34, 36]
const S40i = [6.84, 9.22, 12.48, 15.76, 20.96, 26.64, 35.08, 40.94, 52.48, 62.68, 77.92, 90.12, 102.26, 128.2, 154.08, 202.74, 254.46, 303.18, 333.34, 381, 428.46, 477.82, 575.04, 778.04, 829.04, 875.9]u"mm"
const S40o = [10.3, 13.7, 17.1, 21.3, 26.7, 33.4, 42.2, 48.3, 60.3, 73, 88.9, 101.6, 114.3, 141.3, 168.3, 219.1, 273, 323.8, 355.6, 406.4, 457, 508, 610, 813, 864, 914]u"mm"
const S40t = [1.73, 2.24, 2.31, 2.77, 2.87, 3.38, 3.56, 3.68, 3.91, 5.16, 5.49, 5.74, 6.02, 6.55, 7.11, 8.18, 9.27, 10.31, 11.13, 12.7, 14.27, 15.09, 17.48, 17.48, 17.48, 19.05]u"mm"

const NPS60 = [8, 10, 12, 14, 16, 18, 20, 22, 24]
const S60i = [198.48, 247.6, 295.26, 325.42, 373.08, 418.9, 466.76, 514.54, 560.78]u"mm"
const S60o = [219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610]u"mm"
const S60t = [10.31, 12.7, 14.27, 15.09, 16.66, 19.05, 20.62, 22.23, 24.61]u"mm"

const NPS80 = [0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]
const S80i = [5.48, 7.66, 10.7, 13.84, 18.88, 24.3, 32.5, 38.14, 49.22, 58.98, 73.66, 85.44, 97.18, 122.24, 146.36, 193.7, 242.82, 288.84, 317.5, 363.52, 409.34, 455.62, 501.84, 548.08]u"mm"
const S80o = [10.3, 13.7, 17.1, 21.3, 26.7, 33.4, 42.2, 48.3, 60.3, 73, 88.9, 101.6, 114.3, 141.3, 168.3, 219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610]u"mm"
const S80t = [2.41, 3.02, 3.2, 3.73, 3.91, 4.55, 4.85, 5.08, 5.54, 7.01, 7.62, 8.08, 8.56, 9.53, 10.97, 12.7, 15.09, 17.48, 19.05, 21.44, 23.83, 26.19, 28.58, 30.96]u"mm"

const NPS100 = [8, 10, 12, 14, 16, 18, 20, 22, 24]
const S100i = [188.92, 236.48, 280.92, 307.94, 354.02, 398.28, 442.92, 489.14, 532.22]u"mm"
const S100o = [219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610]u"mm"
const S100t = [15.09, 18.26, 21.44, 23.83, 26.19, 29.36, 32.54, 34.93, 38.89]u"mm"

const NPS120 = [4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]
const S120i = [92.04, 115.9, 139.76, 182.58, 230.12, 273, 300.02, 344.48, 387.14, 431.8, 476.44, 517.96]u"mm"
const S120o = [114.3, 141.3, 168.3, 219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610]u"mm"
const S120t = [11.13, 12.7, 14.27, 18.26, 21.44, 25.4, 27.79, 30.96, 34.93, 38.1, 41.28, 46.02]u"mm"

const NPS140 = [8, 10, 12, 14, 16, 18, 20, 22, 24]
const S140i = [177.86, 222.2, 266.64, 292.1, 333.34, 377.66, 419.1, 463.74, 505.26]u"mm"
const S140o = [219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610]u"mm"
const S140t = [20.62, 25.4, 28.58, 31.75, 36.53, 39.67, 44.45, 47.63, 52.37]u"mm"

const NPS160 = [0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24]
const S160i = [11.74, 15.58, 20.7, 29.5, 34.02, 42.82, 53.94, 66.64, 87.32, 109.54, 131.78, 173.08, 215.84, 257.16, 284.18, 325.42, 366.52, 407.98, 451.04, 490.92]u"mm"
const S160o = [21.3, 26.7, 33.4, 42.2, 48.3, 60.3, 73, 88.9, 114.3, 141.3, 168.3, 219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610]u"mm"
const S160t = [4.78, 5.56, 6.35, 6.35, 7.14, 8.74, 9.53, 11.13, 13.49, 15.88, 18.26, 23.01, 28.58, 33.32, 35.71, 40.49, 45.24, 50.01, 53.98, 59.54]u"mm"

# Schedules designated STD, XS, and XXS
const NPSSTD = [0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48]
const STDi = [6.84, 9.22, 12.48, 15.76, 20.96, 26.64, 35.08, 40.94, 52.48, 62.68, 77.92, 90.12, 102.26, 128.2, 154.08, 202.74, 254.46, 304.74, 336.54, 387.34, 437.94, 488.94, 539.94, 590.94, 640.94, 691.94, 742.94, 793.94, 844.94, 894.94, 945.94, 996.94, 1047.94, 1098.94, 1148.94, 1199.94]u"mm"
const STDo = [10.3, 13.7, 17.1, 21.3, 26.7, 33.4, 42.2, 48.3, 60.3, 73, 88.9, 101.6, 114.3, 141.3, 168.3, 219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610, 660, 711, 762, 813, 864, 914, 965, 1016, 1067, 1118, 1168, 1219]u"mm"
const STDt = [1.73, 2.24, 2.31, 2.77, 2.87, 3.38, 3.56, 3.68, 3.91, 5.16, 5.49, 5.74, 6.02, 6.55, 7.11, 8.18, 9.27, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53, 9.53]u"mm"

const NPSXS = [0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48]
const XSi = [5.48, 7.66, 10.7, 13.84, 18.88, 24.3, 32.5, 38.14, 49.22, 58.98, 73.66, 85.44, 97.18, 122.24, 146.36, 193.7, 247.6, 298.4, 330.2, 381, 431.6, 482.6, 533.6, 584.6, 634.6, 685.6, 736.6, 787.6, 838.6, 888.6, 939.6, 990.6, 1041.6, 1092.6, 1142.6, 1193.6]u"mm"
const XSo = [10.3, 13.7, 17.1, 21.3, 26.7, 33.4, 42.2, 48.3, 60.3, 73, 88.9, 101.6, 114.3, 141.3, 168.3, 219.1, 273, 323.8, 355.6, 406.4, 457, 508, 559, 610, 660, 711, 762, 813, 864, 914, 965, 1016, 1067, 1118, 1168, 1219]u"mm"
const XSt = [2.41, 3.02, 3.2, 3.73, 3.91, 4.55, 4.85, 5.08, 5.54, 7.01, 7.62, 8.08, 8.56, 9.53, 10.97, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7]u"mm"

const NPSXXS = [0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 4, 5, 6, 8, 10, 12]
const XXSi = [6.36, 11.06, 15.22, 22.8, 28, 38.16, 44.96, 58.42, 80.06, 103.2, 124.4, 174.64, 222.2, 273]u"mm"
const XXSo = [21.3, 26.7, 33.4, 42.2, 48.3, 60.3, 73, 88.9, 114.3, 141.3, 168.3, 219.1, 273, 323.8]u"mm"
const XXSt = [7.47, 7.82, 9.09, 9.7, 10.15, 11.07, 14.02, 15.24, 17.12, 19.05, 21.95, 22.23, 25.4, 25.4]u"mm"

const schedule_lookup = Dict(     "5"=> (NPS5, S5i, S5o, S5t),
                            "10"=> (NPS10, S10i, S10o, S10t),
                            "20"=> (NPS20, S20i, S20o, S20t),
                            "30"=> (NPS30, S30i, S30o, S30t),
                            "40"=> (NPS40, S40i, S40o, S40t),
                            "60"=> (NPS60, S60i, S60o, S60t),
                            "80"=> (NPS80, S80i, S80o, S80t),
                            "100"=> (NPS100, S100i, S100o, S100t),
                            "120"=> (NPS120, S120i, S120o, S120t),
                            "140"=> (NPS140, S140i, S140o, S140t),
                            "160"=> (NPS160, S160i, S160o, S160t),
                            "STD"=> (NPSSTD, STDi, STDo, STDt),
                            "XS"=> (NPSXS, XSi, XSo, XSt),
                            "XXS"=> (NPSXXS, XXSi, XXSo, XXSt)
)

function Di_lookup(Di, NPSes, Dis, Dos, ts)
    if Dis[end] < Di
        return error("No hay tubería con diámetro interno igual o mayor a este diámetro")
    else
        for i=1:length(Dis) # Go up ascending list; once larger than specified, return

            if Dis[i] >= Di
                _nps, _di, _do, _t = NPSes[i], Dis[i], Dos[i], ts[i]
                return (_nps, _di, _do, _t)
            end
        end
    end
end

function Do_lookup(Do, NPSes, Dis, Dos, ts)
    if Dos[end] < Do
        return error("No hay tubería con diámetro externo igual o mayor a este diámetro")
    else
        for i=1:length(Dos) # Go up ascending list; once larger than specified, return

            if Dos[i] >= Do
                _nps, _di, _do, _t = NPSes[i], Dis[i], Dos[i], ts[i]
                return (_nps, _di, _do, _t)
            end
        end
    end
end

function NPS_lookup(NPS, NPSes, Dis, Dos, ts)
    i = findfirst(isequal(NPS),NPSes)
    if i==nothing
        error("No hay esa cedula")
    end
    _nps, _di, _do, _t = NPSes[i], Dis[i], Dos[i], ts[i]
    return (_nps, _di, _do, _t)
end

function tuberia(schedule;Do=nothing, Di=nothing, NPS=nothing)
    NPSes, Dis, Dos, ts = schedule_lookup[schedule]
    if Di != nothing
        Di_lookup(Di, NPSes, Dis, Dos, ts)
    elseif Do != nothing
        Do_lookup(Do, NPSes, Dis, Dos, ts)
    elseif NPS != nothing
        NPS_lookup(NPS, NPSes, Dis, Dos, ts)
    end
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(concrete_term),Term,Vector{String},Nothing})   # time: 0.4439514
    Base.precompile(Tuple{typeof(modelcols),MatrixTerm{Tuple{CategoricalTerm{DummyCoding, String, 3}}},NamedTuple{(:X, :Z, :W), Tuple{Vector{String}, Vector{String}, Vector{String}}}})   # time: 0.0605238
    Base.precompile(Tuple{typeof(+),Tuple{Term, Term},Term})   # time: 0.0308432
    Base.precompile(Tuple{Type{ContrastsMatrix},Matrix{Float64},Vector{String},Vector{String},DummyCoding})   # time: 0.0124063
    Base.precompile(Tuple{typeof(schema),Vector{Term},NamedTuple{(:X, :Z, :W), Tuple{Vector{String}, Vector{String}, Vector{String}}},Dict{Symbol, Any}})   # time: 0.0070532
    Base.precompile(Tuple{typeof(getindex),ContrastsMatrix{DummyCoding, String, String},Vector{String},Function})   # time: 0.0040177
    Base.precompile(Tuple{typeof(contrasts_matrix),DummyCoding,Int64,Int64})   # time: 0.0032101
    Base.precompile(Tuple{Type{CategoricalTerm},Symbol,ContrastsMatrix{DummyCoding, String, String}})   # time: 0.0029878
    Base.precompile(Tuple{typeof(~),ConstantTerm{Int64},Tuple{Term, Term, Term}})   # time: 0.0017909
    Base.precompile(Tuple{typeof(termnames),DummyCoding,Vector{String},Int64})   # time: 0.0012678
    Base.precompile(Tuple{typeof(collect_matrix_terms),CategoricalTerm{DummyCoding, String, 3}})   # time: 0.0012
end


text:
const R  = 6372.1  # Radius of the Earth [km]

function haversine_distance(lat1, lon1, lat2, lon2)
    d = 0
    # Your code goes here!
    return d
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
export compose_rotation

################################################################################
#                              Compose Rotations
################################################################################

"""
    compose_rotation(R1, [, R2, R3, R4, R5, ...])

Compute a composed rotation using the rotations `R1`, `R2`, `R3`, `R4`, ..., in
the following order:

     First rotation
     |
     |
    R1 => R2 => R3 => R4 => ...
           |
           |
           Second rotation

The rotations can be described by:

* A direction cosina matrix (`DCM`);
* An Euler angle and axis (`EulerAngleAxis`);
* A set of Euler anlges (`EulerAngles`); or
* A quaternion (`Quaternion`).

Notice, however, that all rotations **must be** of the same type (DCM or
quaternion).

The output will have the same type as the inputs.

# Example

```julia-repl
julia> D1 = angle_to_dcm(+pi/3,+pi/4,+pi/5,:ZYX);

julia> D2 = angle_to_dcm(-pi/5,-pi/4,-pi/3,:XYZ);

julia> compose_rotation(D1,D2)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 1.0          0.0          5.55112e-17
 0.0          1.0          5.55112e-17
 5.55112e-17  5.55112e-17  1.0

julia> ea1 = EulerAngleAxis(30*pi/180, [0;1;0]);

julia> ea2 = EulerAngleAxis(45*pi/180, [0;1;0]);

julia> compose_rotation(ea1,ea2)
EulerAngleAxis{Float64}:
  Euler angle:   1.3090 rad ( 75.0000 deg)
   Euler axis: [  0.0000,   1.0000,   0.0000]

julia> Θ1 = EulerAngles(1,2,3,:ZYX);

julia> Θ2 = EulerAngles(-3,-2,-1,:XYZ);

julia> compose_rotation(Θ1, Θ2)
EulerAngles{Float64}:
  R(X):  -0.0000 rad (  -0.0000 deg)
  R(Y):   0.0000 rad (   0.0000 deg)
  R(Z):  -0.0000 rad (  -0.0000 deg)

julia> q1 = angle_to_quat(+pi/3,+pi/4,+pi/5,:ZYX);

julia> q2 = angle_to_quat(-pi/5,-pi/4,-pi/3,:XYZ);

julia> compose_rotation(q1,q2)
Quaternion{Float64}:
  + 1.0 + 0.0.i + 2.0816681711721685e-17.j + 5.551115123125783e-17.k
```
"""
@inline compose_rotation(D::DCM) = D
@inline compose_rotation(D::DCM, Ds::DCM...) = compose_rotation(Ds...)*D

@inline compose_rotation(ea::EulerAngleAxis) = ea
@inline compose_rotation(ea::EulerAngleAxis, eas::EulerAngleAxis...) =
    compose_rotation(eas...)*ea

@inline compose_rotation(Θ::EulerAngles) = Θ
@inline compose_rotation(Θ::EulerAngles, Θs::EulerAngles...) =
    compose_rotation(Θs...)*Θ

@inline compose_rotation(q::Quaternion) = q
@inline compose_rotation(q::Quaternion, qs::Quaternion...) =
    q*compose_rotation(qs...)

# This algorithm was proposed by @Per in
#
#   https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/24


text:
using JLD
using BH_ED
using Gadfly
using DataFrames
using Cairo

#]add JLD, Gadfly, DataFrames, Cairo    # 安装依赖

# Tuple 格式数据
# result = Tuple{Int, Int, Float64, Float64, Float64, Float64, Float64, Float64, Float64}[]

# DataFrame格式数据
result_df = DataFrame(N = Int[], Lp = Int[], μ = Float64[], U = Float64[], J = Float64[], V = Float64[], ε = Float64[], E = Float64[], np = Float64[])

for N_ in 6, Lp_ in 3, μ_ in 1.0, U_ in 1.0, J_ in 1.0, V_ in 0.1:0.1:2.0, ε_ in 1.0
# for N_ in 2:4, Lp_ in 1:4, μ_ in 0.1:0.2:1.0, U_ in 0.1:0.2:1.0, J_ in 0.1:0.2:1.0, V_ in 0.1:0.2:1.0, ε_ in 0.1:0.2:1.0
    sol = BH_solver(N = N_, Lp = Lp_, μ = μ_, U = U_, J = J_, V = V_, ε = ε_, screen_p = false)
    # push!( result, (N_, Lp_, μ_, U_, J_, V_, ε_, sol[1][1], sol[3])) )
    push!(result_df, ( N_, Lp_, μ_, U_, J_, V_, ε_, sol[1][1], sol[3]) )
end

# 保存数据
# save("data.jld", "result", result_df)

# 绘图
mypl1 = plot(result_df, x = :V, y = :E, Geom.point, Geom.line)
# 保存图片
# img = SVG("N = 6, Lp = 3.svg", 12inch, 8inch) # SVG 格式
# img = PDF("N = 6, Lp = 3.pdf", 12inch, 8inch) # PDF 格式
# img = SVG("N = 6, Lp = 3.png", 12inch, 8inch) # PNG 格式
# draw(img, mypl1)  # 输出文件


text:
@testset "Random" begin
    @testset for k in 1:64
        if k <= 32
            for _ in 1:10
                kmer = rand(DNAMer{k})
                @test isa(kmer, DNAMer{k})
                kmer = rand(RNAMer{k})
                @test isa(kmer, RNAMer{k})
            end
            for size in [0, 1, 2, 5, 10, 100]
                @test length(rand(DNAMer{k}, size)) == size
                @test length(rand(RNAMer{k}, size)) == size
            end
            kmers = rand(DNAMer{k}, 10_000)
            for i in 1:k
                a = sum([kmer[i] for kmer in kmers] .== DNA_A)
                c = sum([kmer[i] for kmer in kmers] .== DNA_C)
                g = sum([kmer[i] for kmer in kmers] .== DNA_G)
                t = sum([kmer[i] for kmer in kmers] .== DNA_T)
                @test 2200 ≤ a ≤ 2800
                @test 2200 ≤ c ≤ 2800
                @test 2200 ≤ g ≤ 2800
                @test 2200 ≤ t ≤ 2800
                @test a + c + g + t == 10_000
            end
        end

        for _ in 1:10
            kmer = rand(BigDNAMer{k})
            @test isa(kmer, BigDNAMer{k})
            kmer = rand(BigRNAMer{k})
            @test isa(kmer, BigRNAMer{k})
        end
        for size in [0, 1, 2, 5, 10, 100]
            @test length(rand(BigDNAMer{k}, size)) == size
            @test length(rand(BigRNAMer{k}, size)) == size
        end
        kmers = rand(BigDNAMer{k}, 10_000)
        for i in 1:k
            a = sum([kmer[i] for kmer in kmers] .== DNA_A)
            c = sum([kmer[i] for kmer in kmers] .== DNA_C)
            g = sum([kmer[i] for kmer in kmers] .== DNA_G)
            t = sum([kmer[i] for kmer in kmers] .== DNA_T)
            @test 2200 ≤ a ≤ 2800
            @test 2200 ≤ c ≤ 2800
            @test 2200 ≤ g ≤ 2800
            @test 2200 ≤ t ≤ 2800
            @test a + c + g + t == 10_000
        end
    end
end


text:
using DataStructures
include("intcode.jl")

"""Paint the hull, putting current tile colour on input and reading output to move robot"""
function paint_hull_io(hull, input, output)
    current_pos = [0 0]
    current_angle = 0
    try
        while true
            #print("First: $(something(hull[current_pos], 0))")
            put!(input, something(hull[current_pos], 0))
            # Paint the hull the appropriate colour based on the first command
            hull[current_pos] = take!(output)
            # Turn the robot the appropriate direction based on the second command
            current_angle = (current_angle + (take!(output) == 0 ? -90 : 90)) % 360
            # Move forward one panel
            current_pos += [Int32(sind(current_angle)) -Int32(cosd(current_angle))]
        end
    catch InvalidStateException end
end

"""Run the hull-painting robot prog on the hull defined as a DefaultDict of coords to colours"""
function paint_hull(hull, prog)
    @sync begin
        input = Channel(2)
        output = Channel(2)
        @async paint_hull_io(hull, input, output)
        @async begin
            run_intcode!(copy(prog), input, output)
            close(input)
            close(output)
        end
    end
end

# Part 1: run the robot and count how many panels it paints
prog = open(f -> parse.(Int64, split(readline(f), ",")), "day11_input.txt")
hull = DefaultDict(nothing)
paint_hull(hull, prog)
println("Number of panels painted: $(length(filter(x -> !isnothing(x), collect(values(hull)))))")

# Part 2: start on a white panel and show the result
hull = DefaultDict(nothing, Dict{Any, Any}([0 0] => 1))
paint_hull(hull, prog)
x_range = (:)(extrema(first.(keys(hull)))...)
y_range = (:)(extrema(last.(keys(hull)))...)
image = zeros(Bool, length(x_range), length(y_range))
for (pos, colour) in hull image[(pos - [first(x_range) first(y_range)] + [1 1])...] = something(colour, 0) end
println.(join.(eachrow(map(x -> x ? "█" : " ", image'))))

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
using Distributed
addprocs(4)
@everywhere begin
    push!(LOAD_PATH,"/home/wanglab/Programs/WhiskerTracking.jl/src")
end
@everywhere using WhiskerTracking, Distributed, Knet

Knet.cuallocator()=false
@time myhandles=make_gui();
@time WhiskerTracking.add_callbacks(myhandles.b,myhandles)


text:
module Stat

const cosmic_fl = "../data/CosmicCompleteExport.tsv"

const data_dir = "../data"

function heredgenes()
    s = Set(readcsv("../data/58gene.csv",ASCIIString))
    push!(s,"BRCA1")
    push!(s,"BRCA2")
    s
end

const panel_fl = "/home/guo/haplox/Bone/data/panel/Cancer_gene.list"

function genespanel()
    genespanel(panel_fl)
end

@doc """ Get sorted genes list from a panel file .list
""" ->
function genespanel(fl)
    genes = Set{ASCIIString}()
    open(fl) do file
        header = readline(file)
        while !eof(file)
            line = readline(file)
            gene = strip(line,'\n')
            if !in(gene,genes)
                push!(genes,gene)
            end
        end
    end
    filter!(gene->gene!="KAT6A" && gene!="MEK1",genes)
    push!(genes,"MYST3")
    push!(genes,"MAP2K1")
    sort(collect(genes))
end

function check_heredgene()
    hdgenes = heredgenes()
    cgenes  = cosmicgenes()
    sygns = Set{ASCIIString}()
    for gn in hdgenes
        if in(gn,cgenes)
            continue
        end
        push!(sygns,gn)
    end
    sygns
end

function check_pnlhdg()
    data = collect(setdiff(Set(genespanel()), heredgenes()))
    writecsv("panel_genes_not_hereditary.csv",sort(data))
end
function cosmicgenes()
    genes = Set{ASCIIString}()
    open(cosmic_fl) do file
        header = readline(file)
        while !eof(file)
            line = readline(file)
            row = split(line,"\t")
            gene = row[1]
            if contains(gene,"_")
                gene = split(gene,"_")[1]
            end
            gene = convert(ASCIIString,gene)
            push!(genes,gene)
        end
    end
    genes
end
function check_pnlgenes()
    pnlgenes = genespanel()
    cosgenes = cosmicgenes()
    for gene in pnlgenes
        if !in(gene,cosgenes)
            @show gene
        end
    end
end

function cancernames()
    cancers = Set{ASCIIString}()
#    cancer_gene = Dict{ASCIIString,Dict{ASCIIString,Arr}}()
    open(cosmic_fl) do file
        header = readline(file)
#        @show split(header,"\t")[17]
#        @show header
        while !eof(file)
            line = readline(file)
            row = split(line,"\t")
 #           gene = convert(ASCIIString, row[1])
            cancer = convert(ASCIIString, row[8])
 #           sample_name = convert(ASCIIString,row[5])
 #           mut = length(row[17]) != 0
 #           @show gene,cancer,sample_name,mut
 #           break
            push!(cancers,cancer)
        end
    end
    collect(cancers)
end

function genesforcancer(cancer)
    gene_samples = Dict{ASCIIString,Dict{ASCIIString,Bool}}()
    samples = Set{ASCIIString}()
    hdgenes = heredgenes()
    pnlgenes = genespanel()
    open(cosmic_fl) do file
        header = readline(file)
        while !eof(file)
            line = readline(file)
            row = split(line,"\t")
            cancer_name = convert(ASCIIString, row[8])
            if cancer_name != cancer
                continue
            end
            sample_name = convert(ASCIIString, row[5])
            if !in(sample_name, samples)
                push!(samples, sample_name)
            end
            gene = convert(ASCIIString, row[1])
            if contains(gene,"_")
                gene = convert(ASCIIString, split(gene,"_")[1])
            end
            if in(gene,hdgenes) || !in(gene, pnlgenes)
                continue
            end
            mut = length(row[17]) != 0
            if !in(gene, keys(gene_samples))
                gene_samples[gene] = Dict{ASCIIString,Bool}(sample_name=>mut)
            else
                if in(sample_name, keys(gene_samples[gene]))
                    gene_samples[gene][sample_name] |= mut
                else
                    gene_samples[gene][sample_name] = mut
                end
            end
        end
    end
    genes = collect(keys(gene_samples)) #first column genename
    samplemuts = map(collect, collect(values(gene_samples))) # [[sp1=>0,sp2=>1],[],[]]
    numsamples = map(length, samplemuts) # third column tol num samples
    nummuts = map(samplemuts) do spmut # second column num mutations
        length(filter(x->x[2]==1,spmut))
    end

    ratios = nummuts ./ numsamples # fourth column of mutation rates

    #=
    samplemuts = map(length,collect(values(gene_samples)))
    data = sort(collect(Dict(zip(genes,counts))),by=x->x[2],rev=true)
    gene = map(x->x[1],data)
    len  = length(samples)
    count = map(x->x[2],data)
    ratio = map(x->x[2]/float(len),data)
    lens  = repmat([len],size(data,1),1)
    =#

    header = ["Gene name","The number of muated samples ",
              "The total number of samples","mutation rates"]

    data = hcat(genes, nummuts, numsamples, ratios)
    data = sortrows(data,by=x->(x[3],x[2]),rev=true)
    data = data[(data[:,4] .>= 0.1) & (data[:,3] .>= 10),:]
    view = vcat(header', data)
    writecsv(joinpath(data_dir, "$(cancer).csv"), view)
end

function genesampleofeachcancer()
    cancers = cancernames()
    pmap(genesforcancer, cancers)
    #run(`zip ../data/StarGene.zip ../data/*.csv`)
end

end


text:
#=
Campo Minado

Calcular o número de minas das células
em volta para todas as células do campo.

Fomato de entrada:

  - Colunas são separadas por espaços
  - Linhas são separadas por("\n")
  - "X" marca uma mina,
    "O"(ó maiúsculo) marca um campo livre de minas

Formato de saída:

  - Campo com "X" são deixados como eles são.
  - Para cada campo contendo "O",
    deve-se calcular o número de minas em volta daquele campo.

    Calculate the number of mines in the surrounding cells for every cell in the field.

Entrada de Exemplo:
O O O O X O O O O O
X X O O O O O O X O
O O O O O O O O O O
O O O O O O O O O O
O O O O O X O O O O

Saída de Exemplo:
2 2 1 1 X 1 0 1 1 1
X X 1 1 1 1 0 1 X 1
2 2 1 0 0 0 0 1 1 1
0 0 0 0 1 1 1 0 0 0
0 0 0 0 1 X 1 0 0 0

=#
include("..\\src\\campo_minado.jl")
include("error_handler_pt_br.jl")

module CampoMinadoTeste
end

using Base.Test

import CampoMinado

Test.with_handler(custom_handler) do
  @test ehValido("O\n") == false
  @test ehValido("X\n") == false
  @test ehValido("O \n") == true
  @test ehValido("X \n") == true
  @test ehValido("")  == true
  @test ehValido("O \nX \n") == true
  @test ehValido("O O O\nX X\n") == false
  @test ehValido("zcakjsfajbsclkak") == false
  @test ehValido("O OOXXXX") == false

  @test paraArray("O X \n") == Vector[Char['O','X']];
  @test paraArray("O \nX \n") == Vector[Char['O'],Char['X']];

  @test_throws ArgumentError campoMinado("asndnasdasddkj")
  @test campoMinado("O \n") == "0 \n"
  @test campoMinado("O X \n") == "1 X \n"
end

max_stars_repo_name: raphaelpanta/julia-lang-exemplos

max_issues_repo_name: raphaelpanta/julia-lang-exemplos

max_forks_repo_name: raphaelpanta/julia-lang-exemplos


text:
module NDReducibles

export ndreducible

using Base:
    StridedReinterpretArray,
    StridedReshapedArray,
    tuple_type_head,
    tuple_type_tail
using Base.Broadcast: Broadcasted

using LinearAlgebra: Adjoint, Transpose
using Referenceables: ReferenceableArray, parenttype
using Transducers:
    @next,
    @return_if_reduced,
    @simd_if,
    Foldable,
    Reduced,
    Transducers,
    complete,
    foldlargs,
    ifunreduced,
    reduced,
    unreduced

include("base.jl")
include("planning.jl")
include("reduce.jl")
include("show.jl")

end # module