text:
using Arrow
using Test, Random, Dates
using CategoricalArrays

const ≅ = isequal

const SEED = 999

const N_OUTER = 4
const N_PAD_TESTS = 32
const N_WRITE_TESTS = 16
const N_DATE_TESTS = 32
const N_IDX_CHECK = 32
const MAX_IDX_LEN = 32
const MAX_VECTOR_LENGTH = 256
const MAX_STRING_LENGTH = 32

const PRIMITIVE_ELTYPES = [Float32, Float64, Int32, Int64, UInt16]
const OFFSET_ELTYPES = [Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, UInt64]

Random.seed!(SEED)

include("testutils.jl")

include("indexing.jl")
include("writing.jl")

# used in locate testing
struct LocateTest1 end
struct LocateTest2 end
struct LocateTest3 end

@testset "locate" begin
    v = rand(Float32, rand(1:MAX_VECTOR_LENGTH))
    p = Primitive(v)

    Locate.length(::LocateTest1) = length(v)
    Locate.values(::LocateTest1) = 1
    Locate.valueslength(::LocateTest1) = length(v)

    data = rawpadded(p, values)

    pp = locate(data, Float32, LocateTest1())
    @test typeof(pp) == Primitive{Float32}
    @test p[:] == pp[:]

    v = Union{Int64,Missing}[2,3,missing,7,11]
    p = NullablePrimitive(v)

    Locate.length(::LocateTest2) = length(v)
    Locate.bitmask(::LocateTest2) = 1
    Locate.values(::LocateTest2) = 9
    Locate.valueslength(::LocateTest2) = length(v)

    data = rawpadded(p, bitmask, values)

    pp = locate(data, Union{Int64,Missing}, LocateTest2())
    @test typeof(pp) == NullablePrimitive{Int64}
    @test p[:] ≅ pp[:]

    v = String["fire", "walk", "with", "me"]
    p = List(v)

    Locate.length(::LocateTest3) = length(v)
    Locate.offsets(::LocateTest3) = 1
    Locate.values(::LocateTest3) = padding((length(v)+1)*sizeof(Arrow.DefaultOffset)) + 1
    Locate.valueslength(::LocateTest3) = 14

    data = rawpadded(p, offsets, values)

    pp = locate(data, String, LocateTest3())
    @test typeof(pp) == List{String,Arrow.DefaultOffset,Primitive{UInt8}}
    @test p[:] == pp[:]
end

@testset "DateTime" begin
    for i ∈ 1:N_DATE_TESTS
        d = Date(rand(0:4000), rand(1:12), rand(1:20))
        ad = convert(Arrow.Datestamp, d)
        @test convert(Date, ad) == d

        dt = DateTime(rand(0:4000), rand(1:12), rand(1:20), rand(0:23), rand(0:59), rand(0:59))
        adt = convert(Arrow.Timestamp, dt)
        @test convert(DateTime, adt) == dt

        t = Time(rand(0:23), rand(0:59), rand(0:59))
        at = convert(Arrow.TimeOfDay, t)
        @test convert(Time, at) == t
    end
end

@testset "arrowformat_construct" begin
    len = rand(1:MAX_VECTOR_LENGTH)
    randstring_() = randstring(rand(1:MAX_STRING_LENGTH))

    p = arrowformat(rand(Float64, len))
    @test typeof(p) == Primitive{Float64}

    mask = rand(Bool, len)

    v = Union{Float32,Missing}[mask[i] ? rand(Float32) : missing for i ∈ 1:len]
    p = arrowformat(v)
    @test typeof(p) == NullablePrimitive{Float32}

    v = String[randstring_() for i ∈ 1:len]
    p = arrowformat(v)
    @test typeof(p) == List{String,Arrow.DefaultOffset,Primitive{UInt8}}

    v = Union{String,Missing}[mask[i] ? randstring_() : missing for i ∈ 1:len]
    p = arrowformat(v)
    @test typeof(p) == NullableList{String,Arrow.DefaultOffset,Primitive{UInt8}}

    v = rand(Bool, len)
    p = arrowformat(v)
    @test typeof(p) == BitPrimitive

    v = Union{Bool,Missing}[mask[i] ? rand(Bool) : missing for i ∈ 1:len]
    p = arrowformat(v)
    @test typeof(p) == NullableBitPrimitive

    v = Date[Date(1), Date(2)]
    p = arrowformat(v)
    @test typeof(p) == Primitive{Arrow.Datestamp}

    v = DateTime[DateTime(0), DateTime(1)]
    p = arrowformat(v)
    @test typeof(p) == Primitive{Arrow.Timestamp{Millisecond}}

    v = Time[Time(0), Time(1)]
    p = arrowformat(v)
    @test typeof(p) == Primitive{Arrow.TimeOfDay{Nanosecond,Int64}}

    v = Union{Date,Missing}[Date(1), missing]
    p = arrowformat(v)
    @test typeof(p) == NullablePrimitive{Arrow.Datestamp}

    v = Union{DateTime,Missing}[DateTime(0), missing]
    p = arrowformat(v)
    @test typeof(p) == NullablePrimitive{Arrow.Timestamp{Millisecond}}

    v = Union{Time,Missing}[Time(0), missing]
    p = arrowformat(v)
    @test typeof(p) == NullablePrimitive{Arrow.TimeOfDay{Nanosecond,Int64}}

    v = categorical(["a", "b", "c"])
    p = arrowformat(v)
    @test typeof(p) == DictEncoding{String,Primitive{Int32},List{String,Arrow.DefaultOffset,Primitive{UInt8}}}
end


text:
using Mandelbrot
using Test

@testset "Mandelbrot.jl" begin
    # Write your tests here.
end

ext: jl


text:
"""
    v = RecursiveVec(vecs)

Create a new vector `v` from an existing (homogeneous or heterogeneous) list of vectors
`vecs` with one or more elements, represented as a `Tuple` or `AbstractVector`. The elements
of `vecs` can be any type of vectors that are supported by KrylovKit. For a heterogeneous
list, it is best to use a tuple for reasons of type stability, while for a homogeneous list,
either a `Tuple` or a `Vector` can be used. From a mathematical perspectve, `v` represents
the direct sum of the vectors in `vecs`. Scalar multiplication and addition of vectors `v`
acts simultaneously on all elements of `v.vecs`. The inner product corresponds to the sum
of the inner products of the individual vectors in the list `v.vecs`.

The vector `v` also adheres to the iteration syntax, but where it will just produce the
individual vectors in `v.vecs`. Hence, `length(v) = length(v.vecs)`. It can also be indexed,
so that `v[i] = v.vecs[i]`, which can be useful in writing a linear map that acts on `v`.
"""
struct RecursiveVec{T<:Union{Tuple,AbstractVector}}
    vecs::T
end
function RecursiveVec(arg1::AbstractVector{T}) where {T}
    if isbitstype(T)
        return RecursiveVec((arg1,))
    else
        return RecursiveVec{typeof(arg1)}(arg1)
    end
end
RecursiveVec(arg1, args...) = RecursiveVec((arg1, args...))

Base.getindex(v::RecursiveVec, i) = v.vecs[i]

Base.iterate(v::RecursiveVec, args...) = iterate(v.vecs, args...)

Base.IteratorEltype(::Type{RecursiveVec{T}}) where {T} = Base.IteratorEltype(T)
Base.IteratorSize(::Type{RecursiveVec{T}}) where {T} = Base.IteratorSize(T)

Base.eltype(v::RecursiveVec) = eltype(v.vecs)
Base.size(v::RecursiveVec) = size(v.vecs)
Base.length(v::RecursiveVec) = length(v.vecs)

Base.first(v::RecursiveVec) = first(v.vecs)
Base.last(v::RecursiveVec) = last(v.vecs)

Base.:-(v::RecursiveVec) = RecursiveVec(map(-, v.vecs))
Base.:+(v::RecursiveVec, w::RecursiveVec) = RecursiveVec(map(+, v.vecs, w.vecs))
Base.:-(v::RecursiveVec, w::RecursiveVec) = RecursiveVec(map(-, v.vecs, w.vecs))
Base.:*(v::RecursiveVec, a::Number) = RecursiveVec(map(x -> x * a, v.vecs))
Base.:*(a::Number, v::RecursiveVec) = RecursiveVec(map(x -> a * x, v.vecs))
Base.:/(v::RecursiveVec, a::Number) = RecursiveVec(map(x -> x / a, v.vecs))
Base.:\(a::Number, v::RecursiveVec) = RecursiveVec(map(x -> a \ x, v.vecs))

function Base.similar(v::RecursiveVec)
    return RecursiveVec(similar.(v.vecs))
end

function Base.copy!(w::RecursiveVec, v::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i in 1:length(w.vecs)
        copyto!(w.vecs[i], v.vecs[i])
    end
    return w
end

function LinearAlgebra.mul!(w::RecursiveVec, a::Number, v::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i in 1:length(w.vecs)
        mul!(w.vecs[i], a, v.vecs[i])
    end
    return w
end

function LinearAlgebra.mul!(w::RecursiveVec, v::RecursiveVec, a::Number)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i in 1:length(w.vecs)
        mul!(w.vecs[i], v.vecs[i], a)
    end
    return w
end

function LinearAlgebra.rmul!(v::RecursiveVec, a::Number)
    for x in v.vecs
        rmul!(x, a)
    end
    return v
end

function LinearAlgebra.axpy!(a::Number, v::RecursiveVec, w::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i in 1:length(w.vecs)
        axpy!(a, v.vecs[i], w.vecs[i])
    end
    return w
end
function LinearAlgebra.axpby!(a::Number, v::RecursiveVec, b, w::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i in 1:length(w.vecs)
        axpby!(a, v.vecs[i], b, w.vecs[i])
    end
    return w
end

LinearAlgebra.dot(v::RecursiveVec{T}, w::RecursiveVec{T}) where {T} =
    sum(dot.(v.vecs, w.vecs))
LinearAlgebra.norm(v::RecursiveVec) = norm(norm.(v.vecs))


text:
module X4J

function xmlaConnect()
	print("xmlaConnect")
end

function xmlaClose()
	print("xmlaClose")
end

function xmlaExecute()
	print("xmlaExecute")
end

function xmlaDiscover()
	print("xmlaDiscover")
end

end # module


text:
using Plots
gr()                   # pyplot()
x = LinRange(0,2π,100) # with 100 points in the interval.
plot(x, sin.(x))
savefig("figSin0.pdf")

max_stars_repo_name: lazarusA/helpfulJulia

max_issues_repo_name: lazarusA/helpfulJulia

max_forks_repo_name: lazarusA/helpfulJulia


text:
@testset "RotatingTimeValue" begin
    rt0 = RotatingTimeValue(1, Date(1))
    rt1 = RotatingTimeValue(1.0, Date(1))
    @test RotatingTimeValue(typeof(rt0), 1.0, Date(1)) === rt0

    @test rotatingtime(1.0, Date(1)) == rt1
    rot = repeat(5:-1:1, inner=5)
    time = repeat(1:5, outer=5)
    rt = rotatingtime(rot, time)
    @test eltype(rt) == RotatingTimeValue{Int, Int}
    @test getfield.(rt, :rotation) == rot
    @test getfield.(rt, :time) == time

    @test isless(rt[1], rt[2])
    @test isless(rt[6], rt[1])
    @test isless(rt[1], rt[7])
    @test isless(rt[1], 2.0)
    @test isless(1, rt[2])
    @test isless(rt[1], missing)
    @test !isless(missing, rt[1])

    @test isequal(rt[1], RotatingTimeValue(5, 1))

    @test rt[1] == RotatingTimeValue(Int32(5), Int32(1))
    @test isequal(rt[1] == RotatingTimeValue(missing, 1), missing)
    @test isequal(rt[1] == RotatingTimeValue(5, missing), missing)
    @test isequal(rt[1] == RotatingTimeValue(1, missing), missing)
    @test rt[1] == 1
    @test 1 == rt[1]
    @test isequal(rt[1] == missing, missing)
    @test isequal(missing == rt[1], missing)

    @test zero(typeof(rt[1])) === RotatingTimeValue(0, 0)
    @test iszero(RotatingTimeValue(1, 0))
    @test !iszero(RotatingTimeValue(0, 1))
    @test one(typeof(rt[1])) === RotatingTimeValue(1, 1)
    @test isone(RotatingTimeValue(1, 1))
    @test !isone(RotatingTimeValue(1, 0))

    @test convert(typeof(rt[1]), RotatingTimeValue(Int32(5), Int32(1))) ===
        RotatingTimeValue(5, 1)

    @test nonmissingtype(RotatingTimeValue{Union{Int,Missing}, Union{Int,Missing}}) ==
        RotatingTimeValue{Int, Int}

    @test iterate(rt1) == (rt1, nothing)
    @test iterate(rt1, 1) === nothing
    @test length(rt1) == 1

    @test sprint(show, rt[1]) == "5_1"
    w = VERSION < v"1.6.0" ? "" : " "
    @test sprint(show, MIME("text/plain"), rt[1]) == """
        RotatingTimeValue{Int64,$(w)Int64}:
          rotation: 5
          time:     1"""
end

@testset "RotatingTimeArray" begin
    rot = collect(1:5)
    time = collect(1.0:5.0)
    a = RotatingTimeArray(rot, time)
    @test eltype(a) == RotatingTimeValue{Int, Float64}
    @test size(a) == (5,)
    @test IndexStyle(typeof(a)) == IndexLinear()

    @test a[1] == RotatingTimeValue(1, 1.0)
    @test a[5:-1:1] == RotatingTimeArray(rot[5:-1:1], time[5:-1:1])

    a[1] = RotatingTimeValue(2, 2.0)
    @test a[1] == RotatingTimeValue(2, 2.0)
    a[1:2] = a[3:4]
    @test a[1:2] == a[3:4]

    v = view(a, 2:4)
    @test v isa RotatingTimeArray
    @test v == a[2:4]

    @test typeof(similar(a)) == typeof(a)
    @test size(similar(a, 3)) == (3,)

    @test v.rotation == rot[2:4]
    @test v.time == time[2:4]

    @test refarray(a).rotation == rot
end


text:
function set_wdot(Qbase, cellxmax, cellymax, Rhat, nch, nre, nval)

    wdot = zeros(cellxmax, cellymax, nval)
    npre = nval - nch

    mw = set_mw()
    avgdro  = set_avgdro_const()

    fmol = set_fmol()
    bmol = set_bmol()

    for i in 2:cellxmax-1
        for j in 2:cellymax-1

            # Initial setting
            total = 0.0
            rho_mw = zeros(nch)
            for s in 1:nch
                rho_mw[s] = Qbase[i,j,npre+s] * Qbase[i,j,1] / mw[s]   # [mol/m3] = [kg/m3] / [kg/mol]
                total     = total + rho_mw[s]
            end
            total = log10(total*avgdro)

            if total <= 20.0
                nsb = 1
            elseif total <= 21.0
                nsb = 2
            elseif total <= 22.0
                nsb = 3
            elseif total <= 23.0
                nsb = 4
            elseif total <= 24.0
                nsb = 5
            else
                nsb = 6
            end

            T   = Qbase[i,j,npre]/(Qbase[i,j,1]*Rhat[i,j])

            keq = k_eq(T, nsb, nre)
            kf  = freaction_rate(T, nre)
            kb  = breaction_rate(kf, keq, nre)

            # Forward and backward rate for reaction ir
            lf = ones(nre)
            lb = ones(nre)
            for r in 1:nre
                for s in 1:nch
                    lf[r] = lf[r]*1.0e-6*rho_mw[s]^fmol[r][s]
                    lb[r] = lb[r]*1.0e-6*rho_mw[s]^bmol[r][s]
                end
                lf[r] = lf[r]*kf[r]
                lb[r] = lb[r]*kb[r]
            end

            for s in 1:nch
                wdot[i,j,npre+s] = 0.0
                for r in 1:nre
                    temp = (bmol[r][s] - fmol[r][s])*(lf[r] - lb[r])
                    wdot[i,j,npre+s] += temp
                end
                wdot[i,j,npre+s] = 1.0e6 * wdot[i,j,npre+s] * mw[s]
            end
        end
    end
    return wdot
end

function set_wdot_for_implicit(Qbase, cellxmax, cellymax, Rhat, nch, nre, nval, volume)

    wdot  = zeros(cellxmax, cellymax, nval)
    H_hat = zeros(cellxmax, cellymax, nval, nval) # source term for point implicit
    npre = nval - nch
    epsion = zeros(nval)
    for s in 1:nch
        epsion[npre+s] = 1.0
    end

    mw = set_mw()
    avgdro  = set_avgdro_const()
    rhos_inv = ones(nval)
    R   = set_gas_const()
    Rs  = zeros(nch)
    for s in 1:nch
        Rs[s] = R/mw[s]
    end

    fmol = set_fmol()
    bmol = set_bmol()

    for i in 2:cellxmax-1
        for j in 2:cellymax-1

            # Initial setting
            total = 0.0
            rho_mw = zeros(nch)
            for s in 1:nch
                rho_mw[s] = Qbase[i,j,npre+s] * Qbase[i,j,1] / mw[s]   # [mol/m3] = [kg/m3] / [kg/mol]
                total     = total + rho_mw[s]
            end
            total = log10(total*avgdro)

            if total <= 20.0
                nsb = 1
            elseif total <= 21.0
                nsb = 2
            elseif total <= 22.0
                nsb = 3
            elseif total <= 23.0
                nsb = 4
            elseif total <= 24.0
                nsb = 5
            else
                nsb = 6
            end

            T   = Qbase[i,j,npre]/(Qbase[i,j,1]*Rhat[i,j])

            keq = k_eq(T, nsb, nre)
            kf  = freaction_rate(T, nre)
            kb  = breaction_rate(kf, keq, nre)

            # Forward and backward rate for reaction ir
            lf = ones(nre)
            lb = ones(nre)
            for r in 1:nre
                for s in 1:nch
                    lf[r] = lf[r]*1.0e-6*rho_mw[s]^fmol[r][s]
                    lb[r] = lb[r]*1.0e-6*rho_mw[s]^bmol[r][s]
                end
                lf[r] = lf[r]*kf[r]
                lb[r] = lb[r]*kb[r]
            end

            for s in 1:nch
                wdot[i,j,npre+s] = 0.0
                for r in 1:nre
                    temp = (bmol[r][s] - fmol[r][s])*(lf[r] - lb[r])
                    wdot[i,j,npre+s] += temp
                end
                wdot[i,j,npre+s] = 1.0e6 * wdot[i,j,npre+s] * mw[s]
            end

            for s in 1:nch
                if Qbase[i,j,4+s] < 1e-20
                    rhos_inv[4+s] = 0.0
                else
                    rhos_inv[4+s] = Qbase[i,j,1] * Qbase[i,j,4+s]
                end
            end

            pTpq = pT_pq(Qbase[i,j,:], T, Rs, nval, nch)
            AZ   = set_AZ(T, nsb)
            sr   = set_sr()
            tr   = set_tr()

            for s in 1:nch
                for l in 1:nval
                    pwpq = 0.0
                    for r in 1:nre
                        tempf = epsion[s]*fmol[r][s]*rhos_inv[s] + (sr[r]*T-tr[r])/T^2 * pTpq[s]
                        tempb = epsion[s]*bmol[r][s]*rhos_inv[s] + ((sr[r]*T-tr[r])/T^2 - AZ/T) * pTpq[s]
                        temp  = (bmol[r][s] - fmol[r][s])*(tempf*lf[r] - tempb*lb[r])
                        pwpq += temp
                    end
                    pwpq = 1.0e6 * wdot[i,j,npre+s] * mw[s]
                    #H_hat[i,j,4+s,l] = pwpq / volume[i,j]
                    H_hat[i,j,4+s,l] = pwpq
                end
            end
        end
    end

    return wdot, H_hat
end

function pT_pq(Qbase_cell, T, Rs, nval, nch)
    pTpq = zeros(nval)

    rhos = zeros(nch)
    rho = Qbase_cell[1]
    u   = Qbase_cell[2]
    v   = Qbase_cell[3]
    for s in 1:nch
        rhos[s] = Qbase_cell[4+s]*rho
    end

    q2         = u^2 + v^2
    delta_hs_0 = set_delta_hs_0()
    Cvhat      = 0.0
    for s in 1:nch
        Cvhat += 1.5*rhos[s]*Rs[s]/rho
    end

    pTpq[1] = 0.5 * q2 / (rho*Cvhat)
    pTpq[2] = -u / (rho*Cvhat)
    pTpq[3] = -v / (rho*Cvhat)
    pTpq[4] =  1 / (rho*Cvhat)
    for s in 1:nch
        pTpq[4+s] = -(delta_hs_0[s] + T*1.5*Rs[s]) / (rho*Cvhat)
    end

    return pTpq
end

function set_AZ(T, nsb)
    ar1, ar2, ar3, ar4, ar5 = set_ar()
    Z  = 10^4/T

    AZ = ar1[nsb]/Z - ar3[nsb] - ar4[nsb]*Z - 2*ar5[nsb]*Z^2
    return AZ
end

function chemical_value(Qbase, cellxmax, cellymax, Rhat, nval, nch)

    chmu        = zeros(cellxmax, cellymax)
    chlambda_tr = zeros(cellxmax, cellymax)
    chD         = zeros(cellxmax, cellymax, nch)
    chmolef         = zeros(cellxmax, cellymax, nch)
    # N2 - N2
    # N2 - N
    # N  - N
    num_nncross = binomial(BigInt(2*nch-1), BigInt(nch)) # 重複組み合わせ nHr = n+r-1Cr
    npre = nval - nch

    mw = set_mw()  # 原子量[kg/mol]
    avgdro  = set_avgdro_const()    # アボガドロ定数
    kb = set_boltzmann_const()      # ボルツマン定数
    pi = 3.1415              # pi

    ms = zeros(nch)          # 原子一個当たりの質量[kg]
    for s in 1:nch
        ms[s] = mw[s] / avgdro
    end

    for i in 1:cellxmax
        for j in 1:cellymax

            total = 0.0
            for s in 1:nch
                total += Qbase[i,j,npre+s]/mw[s]
            end
            for s in 1:nch
                chmolef[i,j,s] = (Qbase[i,j,npre+s]/mw[s]) /total
            end

            T = Qbase[i,j,npre]/(Qbase[i,j,1]*Rhat[i,j])
            p = Qbase[i,j,npre]

            # 対象行列に格納
            # [ N2-N2   N2-N ]
            # [   0.0    N-N ]
            #
            ps11 = collision_cross_section_pisigma11(T, nch)
            ps22 = collision_cross_section_pisigma22(T, nch)

            deltaij1 = zeros(nch, nch)
            deltaij2 = zeros(nch, nch)

            # 以下では行列のすべてを計算しているが、下三角は0なのでキャンセルされる
            for si in 1:nch
                for sj in 1:nch
                    mi = ms[si]
                    mj = ms[sj]

                    temp = 2*mi*mj / (pi*kb*T*(mi+mj))
                    deltaij1[si,sj] = 8/3 * temp^0.5 * ps11[si,sj]
                    deltaij2[si,sj] = 16/5 * temp^0.5 * ps22[si,sj]
                end
            end

            chmu[i,j] = 0.0
            temp      = 0.0
            for si in 1:nch
                for sj in 1:nch
                    temp += chmolef[i,j,sj] * deltaij2[si,sj]
                end
                chmu[i,j] += ms[si]*chmolef[i,j,si]/temp
            end

            chlambda_tr[i,j] = 0.0
            temp =0.0
            for si in 1:nch
                for sj in 1:nch
                    alphaij = 1 + (1-ms[si]/ms[sj]) * (0.45-2.54*ms[si]/ms[sj]) / (1+ms[si]/ms[sj])^2
                    temp += alphaij * chmolef[i,j,sj] * deltaij2[si, sj]
                end
                chlambda_tr[i,j] += chmolef[i,j,si]/temp
            end

            # 対象の行と列を足して(si,si)を引く
            # [ N2-N2   N2-N ]
            # [   0.0    N-N ]
            #
            #=
            (例) N2
                sum[1,si] + sum[si,1] - [si,si]*2 = [ N2-N2  N2-N ] + [ N2-N2 0.0 ] - [N2-N2]*2
                                                  = [ N2-N ]

            Dijを導入すると0割りが発生するため、inverse Dijを計算している
            =#
            inv_Dij = zeros(nch,nch)
            for si in 1:nch
                for sj in 1:nch
                    inv_Dij[si,sj] = (p*deltaij1[si,sj])/kb*T
                end
            end

            for si in 1:nch
                tempsum = 0
                for sj in 1:nch
                    tempsum += chmolef[i,j,sj]*inv_Dij[si,sj]
                    tempsum += chmolef[i,j,sj]*inv_Dij[sj,si]
                end
                tempsum -= chmolef[i,j,si]*inv_Dij[si,si]*2
                rho  = Qbase[i,j,1]
                rhos = Qbase[i,j,npre+si]
                Cs   = rhos / rho
                if chmolef[i,j,si] == 0.0 || tempsum == 0.0
                    chD[i,j,si] = 0.0
                else
                    chD[i,j,si] = (1-Cs)/tempsum * Cs/chmolef[i,j,si]
                end
            end
        end
    end

    return chmu, chlambda_tr, chD, chmolef
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
# Arrays that hold grades and corresponding point values
grades = ["A+", "A", "A-", "B+", "B", "B-", "C+", "C", "C-", "D+", "D", "D-", "F"]
gradepoints = [4.33, 4.0, 3.67, 3.33, 3.0, 2.67, 2.33, 2.0, 1.67, 1.33, 1.0, 0.67, 0.0]
simplegrades = ["A", "B", "C", "D", "E", "F"]
simplegradepoints = [4.0, 3.0, 2.0, 1.0, 0.0]

# Converts a grade to a numerical value
function gradeconvert(grade; expanded=true)
    # The expanded variable determines wheter pluses and minuses
    # are accounted for, or discounted.

    # Stores a translation dictionary
    translation = nothing

    if expanded
        # Expanded Translation
        translation = Dict(
            "A+"    => 4.33,
            "A"     => 4.0,
            "A-"    => 3.67,
            "B+"    => 3.33,
            "B"     => 3.0,
            "B-"    => 2.67,
            "C+"    => 2.33,
            "C"     => 2.0,
            "C-"    => 1.67,
            "D+"    => 1.33,
            "D"     => 1.0,
            "D-"    => 0.67,
            "F"     => 0.0
        )
    else
        # Truncated Translation
        translation = Dict(
            "A+"    => 4.0,
            "A"     => 4.0,
            "A-"    => 3.0,
            "B+"    => 3.0,
            "B"     => 3.0,
            "B-"    => 2.0,
            "C+"    => 2.0,
            "C"     => 2.0,
            "C-"    => 1.0,
            "D+"    => 1.0,
            "D"     => 1.0,
            "D-"    => 0.0,
            "F"     => 0.0
        )
    end

    if in(grade, grades)
        # Checks for a standard grade and tralsate.
        return translation[grade]
    else
        # Otherwise, return 0.0
        return 0.0
    end
end

# Converts a numerical value to a grade
function valueconvert(value; expanded=true)
    # The expanded variable determines wheter pluses and minuses
    # are accounted for, or discounted.

    # Variables to hold grade and point arrays depending on the
    # expanded variable
    g = nothing; p = nothing

    if expanded
        g = grades
        p = gradepoints
    else
        g = simplegrades
        p = simplegradepoints
    end

    # Find point value closests to passed in value
    diffs = abs(p .- value)
    ind = indmin(diffs)

    # Return corresponding grade
    return g[ind]
end

# Converts a numerical value to a grade
function nearestvalue(value; expanded=true)
    # The expanded variable determines wheter pluses and minuses
    # are accounted for, or discounted.

    # Variables to hold grade and point arrays depending on the
    # expanded variable
    g = nothing; p = nothing

    if expanded
        g = grades
        p = gradepoints
    else
        g = simplegrades
        p = simplegradepoints
    end

    # Find point value closests to passed in value
    diffs = abs(p .- value)
    ind = indmin(diffs)

    # Return corresponding grade
    return p[ind]
end

max_stars_repo_name: M-S-H/CASL

max_issues_repo_name: M-S-H/CASL

max_forks_repo_name: M-S-H/CASL


text:
using JLD
traintest = load("CTSV.jld", "simdata")
thetas = traintest[:,1:10]
Zs = traintest[:,11:end]
trainsize = 450000
mZs = mean(Zs,1)
sZs = std(Zs,1)
Zs = (Zs .- mZs) ./ sZs
# get least squares fit using only the training data
x = [ones(size(Zs,1)) Zs]
beta = x\thetas
# least squares errors, both train and test, using coefficients only from train
preprocess = x*beta
errors = thetas - preprocess
# standardize the inputs using stds of training data
sErrors = std(errors,1)
errors = errors ./ sErrors
Y = errors[1:trainsize,:]'
X = Zs[1:trainsize,:]'
YT = errors[trainsize+1:end,:]'
XT = Zs[trainsize+1:end,:]'
preprocess = preprocess[trainsize+1:end,:]
thetas = thetas[trainsize+1:end,:]


text:
tests = ["dists"]

for t in tests
	fn = joinpath("test", "test_$t.jl")
	println("$fn ...")
	include(fn)
end


text:
function rotate(n::Int, c::Char)
    if c in 'a':'z'
        c = 'a' + (c - 'a' + n) % 26
    elseif c in 'A':'Z'
        c = 'A' + (c - 'A' + n) % 26
    end
    return c
end

rotate(n::Int, s::String) = join(rotate(n, c) for c in s)

for n in 0:26
    # current_module() is necessary due to constraints of runtests.jl
    # This is not required for the user solving the exercise
    Core.eval(@__MODULE__, :(macro $(Symbol(:R, n, :_str))(s::String)
        :(rotate($$n, $s))
    end))
end


text:
using MathOptInterface
const MOI = MathOptInterface
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

const MOIU = MOI.Utilities

const SF = Union{MOI.SingleVariable, MOI.ScalarAffineFunction{Float64}, MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}}
const SS = Union{MOI.EqualTo{Float64}, MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone, MOI.ExponentialCone}

struct Solution
    ret_val::Int
    primal::Vector{Float64}
    dual_eq::Vector{Float64}
    dual_ineq::Vector{Float64}
    slack::Vector{Float64}
    objval::Float64
    objbnd::Float64
end
const OPTIMIZE_NOT_CALLED = -1
Solution() = Solution(OPTIMIZE_NOT_CALLED, Float64[], Float64[], Float64[],
                      Float64[], NaN, NaN)

# Used to build the data with allocate-load during `copy_to`.
# When `optimize!` is called, a the data is used to build `ECOSMatrix`
# and the `ModelData` struct is discarded
mutable struct ModelData
    m::Int # Number of rows/constraints
    n::Int # Number of cols/variables
    IA::Vector{Int} # List of conic rows
    JA::Vector{Int} # List of conic cols
    VA::Vector{Float64} # List of conic coefficients
    b::Vector{Float64} # List of conic coefficients
    IG::Vector{Int} # List of equality rows
    JG::Vector{Int} # List of equality cols
    VG::Vector{Float64} # List of equality coefficients
    h::Vector{Float64} # List of equality coefficients
    objconstant::Float64 # The objective is min c'x + objconstant
    c::Vector{Float64}
end

# This is tied to ECOS's internal representation
mutable struct ConeData
    f::Int # number of linear equality constraints
    l::Int # length of LP cone
    q::Int # length of SOC cone
    qa::Vector{Int} # array of second-order cone constraints
    ep::Int # number of primal exponential cone triples
    # The following four field store model information needed to compute `ConstraintPrimal` and `ConstraintDual`
    eqsetconstant::Dict{Int, Float64}   # For the constant of EqualTo
    eqnrows::Dict{Int, Int}             # The number of rows of Zeros
    ineqsetconstant::Dict{Int, Float64} # For the constant of LessThan and GreaterThan
    ineqnrows::Dict{Int, Int}           # The number of rows of each vector sets except Zeros
    function ConeData()
        new(0, 0, 0, Int[], 0,
            Dict{Int, Float64}(),
            Dict{Int, UnitRange{Int}}(),
            Dict{Int, Float64}(),
            Dict{Int, UnitRange{Int}}())
    end
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    cone::ConeData
    maxsense::Bool
    data::Union{Nothing, ModelData} # only non-Nothing between MOI.copy_to and MOI.optimize!
    sol::Solution
    options
    function Optimizer(; kwargs...)
        new(ConeData(), false, nothing, Solution(), kwargs)
    end
end

MOI.get(::Optimizer, ::MOI.SolverName) = "ECOS"

function MOI.is_empty(instance::Optimizer)
    !instance.maxsense && instance.data === nothing
end

function MOI.empty!(instance::Optimizer)
    instance.maxsense = false
    instance.data = nothing # It should already be nothing except if an error is thrown inside copy_to
    instance.sol = Solution()
end

MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names

function MOI.supports(::Optimizer,
                      ::Union{MOI.ObjectiveSense,
                              MOI.ObjectiveFunction{MOI.SingleVariable},
                              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}})
    return true
end

MOI.supports_constraint(::Optimizer, ::Type{<:SF}, ::Type{<:SS}) = true

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

using Compat.SparseArrays

const ZeroCones = Union{MOI.EqualTo, MOI.Zeros}
const LPCones = Union{MOI.GreaterThan, MOI.LessThan, MOI.Nonnegatives, MOI.Nonpositives}

# Computes cone dimensions
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, <:ZeroCones}) = ci.value
function _allocate_constraint(cone::ConeData, f, s::ZeroCones)
    ci = cone.f
    cone.f += MOI.dimension(s)
    ci
end
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, <:LPCones}) = ci.value
function _allocate_constraint(cone::ConeData, f, s::LPCones)
    ci = cone.l
    cone.l += MOI.dimension(s)
    ci
end
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, <:MOI.SecondOrderCone}) = cone.l + ci.value
function _allocate_constraint(cone::ConeData, f, s::MOI.SecondOrderCone)
    push!(cone.qa, s.dimension)
    ci = cone.q
    cone.q += MOI.dimension(s)
    ci
end
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, <:MOI.ExponentialCone}) = cone.l + cone.q + ci.value
function _allocate_constraint(cone::ConeData, f, s::MOI.ExponentialCone)
    ci = 3cone.ep
    cone.ep += 1
    ci
end
constroffset(instance::Optimizer, ci::CI) = constroffset(instance.cone, ci::CI)
function MOIU.allocate_constraint(instance::Optimizer, f::F, s::S) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    CI{F, S}(_allocate_constraint(instance.cone, f, s))
end

# Build constraint matrix
scalecoef(rows, coef, minus, s) = minus ? -coef : coef
scalecoef(rows, coef, minus, s::Union{MOI.LessThan, Type{<:MOI.LessThan}, MOI.Nonpositives, Type{MOI.Nonpositives}}) = minus ? coef : -coef
output_index(t::MOI.VectorAffineTerm) = t.output_index
variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value
variable_index_value(t::MOI.VectorAffineTerm) = variable_index_value(t.scalar_term)
coefficient(t::MOI.ScalarAffineTerm) = t.coefficient
coefficient(t::MOI.VectorAffineTerm) = coefficient(t.scalar_term)
constrrows(::MOI.AbstractScalarSet) = 1
constrrows(s::MOI.AbstractVectorSet) = 1:MOI.dimension(s)
constrrows(instance::Optimizer, ci::CI{<:MOI.AbstractScalarFunction, <:MOI.AbstractScalarSet}) = 1
constrrows(instance::Optimizer, ci::CI{<:MOI.AbstractVectorFunction, MOI.Zeros}) = 1:instance.cone.eqnrows[constroffset(instance, ci)]
constrrows(instance::Optimizer, ci::CI{<:MOI.AbstractVectorFunction, <:MOI.AbstractVectorSet}) = 1:instance.cone.ineqnrows[constroffset(instance, ci)]
matrix(data::ModelData, s::ZeroCones) = data.b, data.IA, data.JA, data.VA
matrix(data::ModelData, s::Union{LPCones, MOI.SecondOrderCone, MOI.ExponentialCone}) = data.h, data.IG, data.JG, data.VG
matrix(instance::Optimizer, s) = matrix(instance.data, s)
MOIU.load_constraint(instance::Optimizer, ci, f::MOI.SingleVariable, s) = MOIU.load_constraint(instance, ci, MOI.ScalarAffineFunction{Float64}(f), s)
function MOIU.load_constraint(instance::Optimizer, ci, f::MOI.ScalarAffineFunction, s::MOI.AbstractScalarSet)
    a = sparsevec(variable_index_value.(f.terms), coefficient.(f.terms))
    # sparsevec combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(a)
    offset = constroffset(instance, ci)
    row = constrrows(s)
    i = offset + row
    # The ECOS format is b - Ax ∈ cone
    # so minus=false for b and minus=true for A
    setconstant = MOIU.getconstant(s)
    if s isa MOI.EqualTo
        instance.cone.eqsetconstant[offset] = setconstant
    else
        instance.cone.ineqsetconstant[offset] = setconstant
    end
    constant = f.constant - setconstant
    b, I, J, V = matrix(instance, s)
    b[i] = scalecoef(row, constant, false, s)
    append!(I, fill(i, length(a.nzind)))
    append!(J, a.nzind)
    append!(V, scalecoef(row, a.nzval, true, s))
end
MOIU.load_constraint(instance::Optimizer, ci, f::MOI.VectorOfVariables, s) = MOIU.load_constraint(instance, ci, MOI.VectorAffineFunction{Float64}(f), s)
# ECOS orders differently than MOI the second and third dimension of the exponential cone
orderval(val, s) = val
function orderval(val, s::Union{MOI.ExponentialCone, Type{MOI.ExponentialCone}})
    val[[1, 3, 2]]
end
orderidx(idx, s) = idx
expmap(i) = (1, 3, 2)[i]
function orderidx(idx, s::MOI.ExponentialCone)
    expmap.(idx)
end
function MOIU.load_constraint(instance::Optimizer, ci, f::MOI.VectorAffineFunction, s::MOI.AbstractVectorSet)
    A = sparse(output_index.(f.terms), variable_index_value.(f.terms), coefficient.(f.terms))
    # sparse combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(A)
    I, J, V = findnz(A)
    offset = constroffset(instance, ci)
    rows = constrrows(s)
    if s isa MOI.Zeros
        instance.cone.eqnrows[offset] = length(rows)
    else
        instance.cone.ineqnrows[offset] = length(rows)
    end
    i = offset .+ rows
    # The ECOS format is b - Ax ∈ cone
    # so minus=false for b and minus=true for A
    b, Is, Js, Vs = matrix(instance, s)
    b[i] .= scalecoef(rows, orderval(f.constants, s), false, s)
    append!(Is, offset .+ orderidx(I, s))
    append!(Js, J)
    append!(Vs, scalecoef(I, V, true, s))
end

function MOIU.allocate_variables(instance::Optimizer, nvars::Integer)
    instance.cone = ConeData()
    VI.(1:nvars)
end

function MOIU.load_variables(instance::Optimizer, nvars::Integer)
    cone = instance.cone
    m = cone.l + cone.q + 3cone.ep
    IA = Int[]
    JA = Int[]
    VA = Float64[]
    b = zeros(cone.f)
    IG = Int[]
    JG = Int[]
    VG = Float64[]
    h = zeros(m)
    c = zeros(nvars)
    instance.data = ModelData(m, nvars, IA, JA, VA, b, IG, JG, VG, h, 0., c)
end

function MOIU.allocate(instance::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    instance.maxsense = sense == MOI.MaxSense
end
function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction,
                       ::MOI.Union{MOI.SingleVariable,
                                   MOI.ScalarAffineFunction{Float64}})
end

function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.SingleVariable)
    MOIU.load(optimizer,
              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
              MOI.ScalarAffineFunction{Float64}(f))
end
function MOIU.load(instance::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.ScalarAffineFunction)
    c0 = Vector(sparsevec(variable_index_value.(f.terms), coefficient.(f.terms),
                          instance.data.n))
    instance.data.objconstant = f.constant
    instance.data.c = instance.maxsense ? -c0 : c0
    return nothing
end

function MOI.optimize!(instance::Optimizer)
    if instance.data === nothing
        # optimize! has already been called and no new model has been copied
        return
    end
    cone = instance.cone
    m = instance.data.m
    n = instance.data.n
    A = ECOS.ECOSMatrix(sparse(instance.data.IA, instance.data.JA, instance.data.VA, cone.f, n))
    b = instance.data.b
    G = ECOS.ECOSMatrix(sparse(instance.data.IG, instance.data.JG, instance.data.VG, m, n))
    h = instance.data.h
    objconstant = instance.data.objconstant
    c = instance.data.c
    instance.data = nothing # Allows GC to free instance.data before A is loaded to ECOS
    ecos_prob_ptr = ECOS.setup(n, m, cone.f, cone.l, length(cone.qa), cone.qa,
                               cone.ep, G, A, c, h, b; instance.options...)
    ret_val = ECOS.solve(ecos_prob_ptr)
    ecos_prob = unsafe_wrap(Array, ecos_prob_ptr, 1)[1]
    primal    = unsafe_wrap(Array, ecos_prob.x, n)[:]
    dual_eq   = unsafe_wrap(Array, ecos_prob.y, cone.f)[:]
    dual_ineq = unsafe_wrap(Array, ecos_prob.z, m)[:]
    slack     = unsafe_wrap(Array, ecos_prob.s, m)[:]
    ECOS.cleanup(ecos_prob_ptr, 0)
    objval = (instance.maxsense ? -1 : 1) * dot(c, primal)
    if ret_val != ECOS.ECOS_DINF
        objval += objconstant
    end
    objbnd = -(dot(b, dual_eq) + dot(h, dual_ineq))
    if ret_val != ECOS.ECOS_PINF
        objbnd += objconstant
    end
    instance.sol = Solution(ret_val, primal, dual_eq, dual_ineq, slack, objval,
                            objbnd)
end

# Implements getter for result value and statuses
function MOI.get(instance::Optimizer, ::MOI.TerminationStatus)
    flag = instance.sol.ret_val
    if flag == OPTIMIZE_NOT_CALLED
        return MOI.OptimizeNotCalled
    elseif flag == ECOS.ECOS_OPTIMAL
        return MOI.Optimal
    elseif flag == ECOS.ECOS_PINF
        return MOI.Infeasible
    elseif flag == ECOS.ECOS_DINF
        return MOI.DualInfeasible
    elseif flag == ECOS.ECOS_MAXIT
        return MOI.IterationLimit
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        return MOI.AlmostOptimal
    elseif flag == ECOS.ECOS_PINF + ECOS.ECOS_INACC_OFFSET
        return MOI.AlmostInfeasible
    else
        return MOI.OtherError
    end
    # TODO: AlmostDualInfeasible for ECOS.ECOS_DINF + ECOS.ECOS_INACC_OFFSET
    # https://github.com/JuliaOpt/MathOptInterface.jl/issues/601
end

MOI.get(instance::Optimizer, ::MOI.ObjectiveValue) = instance.sol.objval
MOI.get(instance::Optimizer, ::MOI.ObjectiveBound) = instance.sol.objbnd

function MOI.get(instance::Optimizer, ::MOI.PrimalStatus)
    flag = instance.sol.ret_val
    if flag == ECOS.ECOS_OPTIMAL
        return MOI.FeasiblePoint
    elseif flag == ECOS.ECOS_PINF
        return MOI.InfeasiblePoint
    elseif flag == ECOS.ECOS_DINF
        return MOI.InfeasibilityCertificate
    elseif flag == ECOS.ECOS_MAXIT
        return MOI.UnknownResultStatus
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        return MOI.NearlyFeasiblePoint
    elseif flag == ECOS.ECOS_PINF + ECOS.ECOS_INACC_OFFSET
        return MOI.InfeasiblePoint
    elseif flag == ECOS.ECOS_DINF + ECOS.ECOS_INACC_OFFSET
        return MOI.NearlyInfeasibilityCertificate
    else
        return MOI.OtherResultStatus
    end
end
# Swapping indices 2 <-> 3 is an involution (it is its own inverse)
const reorderval = orderval
function MOI.get(instance::Optimizer, ::MOI.VariablePrimal, vi::VI)
    instance.sol.primal[vi.value]
end
MOI.get(instance::Optimizer, a::MOI.VariablePrimal, vi::Vector{VI}) = MOI.get.(instance, Ref(a), vi)
# setconstant: Retrieve set constant stored in `ConeData` during `copy_to`
setconstant(instance::Optimizer, offset, ::CI{<:MOI.AbstractFunction, <:MOI.EqualTo}) = instance.cone.eqsetconstant[offset]
setconstant(instance::Optimizer, offset, ::CI) = instance.cone.ineqsetconstant[offset]
_unshift(instance::Optimizer, offset, value, ::CI) = value
_unshift(instance::Optimizer, offset, value, ci::CI{<:MOI.AbstractScalarFunction, <:MOI.AbstractScalarSet}) = value + setconstant(instance, offset, ci)
function MOI.get(instance::Optimizer, ::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, MOI.Zeros})
    rows = constrrows(instance, ci)
    zeros(length(rows))
end
function MOI.get(instance::Optimizer, ::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, <:MOI.EqualTo})
    offset = constroffset(instance, ci)
    setconstant(instance, offset, ci)
end
function MOI.get(instance::Optimizer, ::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(instance, ci)
    rows = constrrows(instance, ci)
    _unshift(instance, offset, scalecoef(rows, reorderval(instance.sol.slack[offset .+ rows], S), false, S), ci)
end

function MOI.get(instance::Optimizer, ::MOI.DualStatus)
    flag = instance.sol.ret_val
    if flag == ECOS.ECOS_OPTIMAL
        return MOI.FeasiblePoint
    elseif flag == ECOS.ECOS_PINF
        return MOI.InfeasibilityCertificate
    elseif flag == ECOS.ECOS_DINF
        return MOI.InfeasiblePoint
    elseif flag == ECOS.ECOS_MAXIT
        return MOI.UnknownResultStatus
    elseif flag == ECOS.ECOS_OPTIMAL + ECOS.ECOS_INACC_OFFSET
        return MOI.NearlyFeasiblePoint
    elseif flag == ECOS.ECOS_PINF + ECOS.ECOS_INACC_OFFSET
        return MOI.NearlyInfeasibilityCertificate
    elseif flag == ECOS.ECOS_DINF + ECOS.ECOS_INACC_OFFSET
        return MOI.InfeasiblePoint
    else
        m.solve_stat = MOI.OtherResultStatus
    end
end
_dual(instance, ci::CI{<:MOI.AbstractFunction, <:ZeroCones}) = instance.sol.dual_eq
_dual(instance, ci::CI) = instance.sol.dual_ineq
function MOI.get(instance::Optimizer, ::MOI.ConstraintDual, ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(instance, ci)
    rows = constrrows(instance, ci)
    scalecoef(rows, reorderval(_dual(instance, ci)[offset .+ rows], S), false, S)
end

MOI.get(instance::Optimizer, ::MOI.ResultCount) = 1


text:
for (@var i = 0; i < triangles; i ++)
    positions.push(Math.random() - 0.5 )
    positions.push(Math.random() - 0.5 )
    positions.push(Math.random() - 0.5 )
    colors.push( Math.random() * 255 )
    colors.push( Math.random() * 255 )
    colors.push( Math.random() * 255 )
    colors.push( Math.random() * 255 )
end

@var triangles = 500
positions = Observable(threejs, "positions", rand(triangles * 3) .- 0.5
colors = Observable(threejs, "colors", rand(triangles * 4) * 255

vert_shader = Observable(threejs, "vert_shader", """
precision mediump float;
precision mediump int;

uniform mat4 modelViewMatrix; # optional
uniform mat4 projectionMatrix; # optional

attribute vec3 position;
attribute vec4 color;

varying vec3 vPosition;
varying vec4 vColor;

void main(){
    vPosition = position;
    vColor = color;
    gl_Position = projectionMatrix * modelViewMatrix * vec4(position, 1.0);
}
""")

frag_shader = Observable(threejs, "frag_shader", """
precision mediump float;
precision mediump int;
uniform float time;
varying vec3 vPosition;
varying vec4 vColor;
void main(){
    vec4 color = vec4(vColor);
    color.r += sin(vPosition.x * 10.0 + time) * 0.5;
    gl_FragColor = color;
}
""")

@var container, stats;
@var camera, scene, renderer;
function onWindowResize(event)
    camera.aspect = window.innerWidth / window.innerHeight
    camera.updateProjectionMatrix()
    renderer.setSize(window.innerWidth, window.innerHeight)
end
function animate()
    requestAnimationFrame(animate)
    render()
    stats.update()
end
function render()
    @var time = performance.now()
    @var object = scene.children[0]
    object.rotation.y = time * 0.0005
    object.material.uniforms.time.value = time * 0.005
    renderer.render(scene, camera)
end

function init()
    @var container = document.querySelector("#container");
    camera = @new THREE.PerspectiveCamera(50, window.innerWidth / window.innerHeight, 1, 10 );
    camera.position.z = 2;
    scene = @new THREE.Scene();
    scene.background = @new THREE.Color(0x101010)
    # geometry
    @var triangles = 500
    @var geometry = @new THREE.BufferGeometry()
    @var positions = []
    @var colors = []

    @var positionAttribute = @new THREE.Float32BufferAttribute(positions, 3);
    @var colorAttribute = @new THREE.Uint8BufferAttribute(colors, 4);
    colorAttribute.normalized = true; # this will map the buffer values to 0.0f - +1.0f in the shader
    geometry.addAttribute("position", positionAttribute);
    geometry.addAttribute("color", colorAttribute);
    # material
    @var material = @new THREE.RawShaderMaterial(d(
        uniforms = d(
            time = d(value = 1.0)
        ),
        vertexShader = $vert_shader[],
        fragmentShader = $frag_shader[],
        side = THREE.DoubleSide,
        transparent = true
    ))
    @var mesh = @new THREE.Mesh(geometry, material)
    scene.add(mesh)
    renderer = @new THREE.WebGLRenderer()
    renderer.setPixelRatio( window.devicePixelRatio )
    renderer.setSize( window.innerWidth, window.innerHeight )
    container.appendChild( renderer.domElement )
    stats = @new Stats()
    container.appendChild(stats.dom)
    window.addEventListener("resize", onWindowResize, false)
end
init()
animate()


text:
using AssociativeWindowAggregation
using OnlineStatsBase
using Plots
using Statistics

struct WindowedStat{Stat} <: OnlineStat{Stat}
    state::FixedWindowAssociativeOp{Stat,merge,merge!}

    function WindowedStat{Stat}(window::Int) where {Stat<:OnlineStat}
        return new{Stat}(FixedWindowAssociativeOp{Stat,merge,merge!}(window))
    end
end

function OnlineStatsBase.fit!(x::WindowedStat{Stat}, v) where {Stat}
    # Create a single-observation instance of the statistic.
    wrapped = Stat()
    fit!(wrapped, v)
    update_state!(x.state, wrapped)
    return x
end

OnlineStatsBase.value(x::WindowedStat) = window_value(x.state)
OnlineStatsBase.nobs(x::WindowedStat) = window_size(x.state)

function Base.show(io::IO, o::WindowedStat{Stat}) where {Stat}
    nobs(o) == 0 && return print(io, "WindowStat{$Stat}()")
    return invoke(show, Tuple{IO,OnlineStat}, io, o)
end

x = WindowedStat{Variance}(100)
values = rand(1000) .- 0.5
output = [std(value(fit!(x, v))) for v in values]

begin
    plot(values; alpha=0.4, label="Input")
    plot!(output; label="std (window=100)")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ fb37684a-0e82-11ec-227f-4595b4e5fa74
begin

	using Pkg
	 # activate the shared project environment
    Pkg.activate("/home/graham_s/julia/vs/ScottishTaxBenefitMode/")
	# Pkg.activate(Base.current_project())
	# Pkg.add( "https://github.com/grahamstark/ScottishTaxBenefitModel.jl" )
	using ScottishTaxBenefitModel
	using .Definitions
	using .RunSettings: Settings
	using CSV,DataFrames,GLM,RegressionTables

end

# ╔═╡ f00633e2-af9c-46cf-b4e6-2bf403d10392
begin

const settings = Settings()

frshh = CSV.File("$(MODEL_DATA_DIR)/model_households.tab" ) |> DataFrame
frspeople = CSV.File("$(MODEL_DATA_DIR)/model_people.tab") |> DataFrame

end

# ╔═╡ dc101061-f980-4fba-b223-8bbb88a7218f
begin

fm = innerjoin( frshh, frspeople, on=[:data_year, :hid ], makeunique=true )
fm.age_sq = fm.age.^2
fm.deaf_blind=fm.registered_blind .| fm.registered_deaf .| fm.registered_partially_sighted
fm.yr = fm.data_year .- 2014
fm.any_dis = (
    fm.disability_vision .|
    fm.disability_hearing .|
    fm.disability_mobility .|
    fm.disability_dexterity .|
    fm.disability_learning .|
    fm.disability_memory .|
    fm.disability_other_difficulty .|
    fm.disability_mental_health .|
    fm.disability_stamina .|
    fm.disability_socially )
fm.adls_bad=fm.adls_are_reduced.==1
fm.adls_mid=fm.adls_are_reduced.==2
fm.rec_carers = fm.income_carers_allowance .> 0
fm.rec_dla = ( fm.income_dlamobility.>0.0) .| ( fm.income_dlaself_care .>0.0 )
fm.rec_dla_care = ( fm.income_dlaself_care .>0.0 )
fm.rec_dla_mob = ( fm.income_dlamobility.>0.0 )
fm.rec_pip = ( fm.income_personal_independence_payment_mobility.>0.0) .| ( fm.income_personal_independence_payment_daily_living .>0.0 )
fm.rec_pip_care = ( fm.income_personal_independence_payment_daily_living .>0.0 )
fm.rec_pip_mob = ( fm.income_personal_independence_payment_mobility.>0.0)
fm.rec_esa = ( fm.income_employment_and_support_allowance.>0.0)
fm.rec_aa = ( fm.income_attendance_allowance.>0.0)
fm.rec_carers = ( fm.income_carers_allowance.>0.0)
fm_rec_aa = ( fm.income_attendance_allowance.>0.0)
fm.scotland = fm.region .== 299999999
fm.male = fm.sex .== 1

# subsets
fm_working_age = fm[((fm.age .< 65).&(fm.from_child_record.!= 1)),:]
fm_pension_age = fm[(fm.age .>= 65),:]

sum( fm.scotland )
	size( fm )

end

# ╔═╡ 5d54d384-2b8e-4729-ab4a-dba5858c6f10
sum(fm_pension_age.rec_aa)

# ╔═╡ 187d3bbe-489b-485b-9db6-b3a44ef22c49
begin

bens = [
  :rec_carers,
  :rec_pip,
  :rec_pip_mob,
  :rec_pip_care,
  :rec_esa,
  :rec_dla,
  :rec_dla_mob,
  :rec_dla_care,
  :rec_aa ]

datasets = Dict(
  :rec_carers=>fm,
  :rec_pip=>fm_working_age,
  :rec_pip_mob=>fm_working_age,
  :rec_pip_care=>fm_working_age,
  :rec_esa=>fm_working_age,
  :rec_dla=>fm_working_age,
  :rec_dla_mob=>fm_working_age,
  :rec_dla_care=>fm_working_age,
  :rec_aa=>fm_pension_age )

	main_out=[]

	for target in bens
		  v = glm(
				Term( target ) ~
				  Term( :scotland ) +
				  Term( :data_year )+
				  Term( :male )+
				  Term( :age )+
				  Term( :age_sq ) +
				  Term( :has_long_standing_illness )+
				  Term( :adls_bad )+
				  Term( :adls_mid )+
				  Term( :registered_blind )+
				  Term( :registered_partially_sighted )+
				  Term( :registered_deaf )+

				  Term( :disability_dexterity )+
				  Term( :disability_learning )+
				  Term( :disability_memory )+
				  Term( :disability_mental_health )+
				  Term( :disability_mobility )+
				  Term( :disability_socially )+
				  Term( :disability_stamina )+
				  Term( :disability_other_difficulty ),
			  datasets[target ],
			  Binomial(),
			  ProbitLink(), contrasts=Dict( :data_year=>DummyCoding()))
		push!( main_out, v )
	end
	main_out
end

# ╔═╡ faa4a768-baf0-4c79-a2e3-1c18b2ce903b
dla_1=glm(@formula(rec_dla ~ scotland + male+age+age^2+has_long_standing_illness+data_year+adls_bad+adls_mid+registered_blind+registered_partially_sighted+registered_deaf+disability_dexterity+disability_memory+disability_mental_health +disability_mobility+disability_socially+disability_stamina+disability_learning+disability_other_difficulty), fm_working_age, Binomial(), ProbitLink(), contrasts=Dict( :data_year=>DummyCoding()) )

# ╔═╡ Cell order:
# ╠═fb37684a-0e82-11ec-227f-4595b4e5fa74
# ╠═f00633e2-af9c-46cf-b4e6-2bf403d10392
# ╠═dc101061-f980-4fba-b223-8bbb88a7218f
# ╠═5d54d384-2b8e-4729-ab4a-dba5858c6f10
# ╠═187d3bbe-489b-485b-9db6-b3a44ef22c49
# ╠═faa4a768-baf0-4c79-a2e3-1c18b2ce903b

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function SphericalTriangle{T}(
    domain::SphericalTriangle,
    parent::ArbField,
) where {T<:Union{fmpq,arb}}
    SphericalTriangle(domain.angles, parent)
end

function Base.show(io::IO, domain::SphericalTriangle{fmpq})
    angles_string = ["$(numerator(a))π/$(denominator(a))" for a in domain.angles]
    print(
        io,
        "Spherical triangle with angles ($(angles_string[1]), $(angles_string[2]), $(angles_string[3]))",
    )
end

function Base.show(io::IO, domain::SphericalTriangle{arb})
    print(io, "Spherical triangle with angles $(domain.angles)")
end

"""
    angles(domain::SphericalTriangle)
> Return the angles of the spherical triangle.
"""
function angles(domain::SphericalTriangle{fmpq})
    domain.parent(π) .* domain.angles
end

function angles(domain::SphericalTriangle{arb})
    domain.angles
end

"""
    angle(domain::SphericalTriangle, i::Integer)
> Return the angle for vertex i of the spherical triangle.
"""
function angle(domain::SphericalTriangle{fmpq}, i::Integer)
    domain.parent(π) * domain.angles[i]
end

function angle(domain::SphericalTriangle{arb}, i::Integer)
    domain.angles[i]
end

"""
    vertex(domain::SphericalTriangle, i::Integer)
> Compute vertex i of the spherical triangle.
"""
function vertex(domain::SphericalTriangle, i::Integer)
    RR = domain.parent
    if i == 1
        return RR.(SVector(0, 0, 1))
    elseif i == 2
        α, β, γ = angles(domain)
        S = (α + β + γ) / 2
        θ = 2asin(sqrt(-cos(S) * cos(S - γ) / (sin(α) * sin(β))))
        return SVector(sin(θ), RR(0), cos(θ))
    elseif i == 3
        α, β, γ = angles(domain)
        S = (α + β + γ) / 2
        θ = 2asin(sqrt(-cos(S) * cos(S - β) / (sin(α) * sin(γ))))
        return SVector(sin(θ) * cos(α), sin(θ) * sin(α), cos(θ))
    else
        throw(ErrorException("attempt to get vertex $i from a spherical triangle"))
    end
end

function area(domain::SphericalTriangle)
    sum(angles(domain)) - domain.parent(π)
end

"""
    center(domain::SphericalTriangle)
> Compute the center of the spherical triangle, given by the
  normalised sum of the vertices.
"""
function center(domain::SphericalTriangle)
    xyz = sum(vertex(domain, i) for i = 1:3)
    r = sqrt(sum(xyz .^ 2))
    xyz ./ r
end

"""
    greatcircleplane(domain::SphericalTriangle)
> Compute coefficients a, b, c giving the plane ax + by + cz = 0 which
  is parallel to the great circle that defines the bottom boundary of
  the spherical triangle.
"""
function greatcircleplane(domain::SphericalTriangle)
    # Compute the normal vector of the plane
    w = cross(vertex(domain, 2), vertex(domain, 3))
    # w = [a, b, c]
    return (w[1], w[2], w[3])
end

"""
    greatcircle(ϕ::arb, a::arb, b::arb, c::arb)
> Compute the θ value corresponding to the give ϕ value on the great
circle defined by the plane ax + by + cz = 0.
"""
function greatcircle(ϕ::arb, a::arb, b::arb, c::arb)
    θ = atan(-c, a * cos(ϕ) + b * sin(ϕ))
    if θ < 0
        θ += ϕ.parent(π)
    end
    return θ
end

"""
    theta_bound(domain::SphericalTriangle, type = :lower)
> Compute a bound of θ on the boundary opposite of the north pole. By
  default a lower bound is computed, an upper bound is computed with
  `type = :upper`.
"""
function theta_bound(domain::SphericalTriangle, type = :lower)
    if type == :lower
        cmp = min
    elseif type == :upper
        cmp = max
    else
        throw(ArgumentError("type must be :lower or :upper, got $type"))
    end

    # Compute the critical point with respect to θ on the boundary
    i = 1
    u = vertex(domain, 2)
    v = vertex(domain, 3)
    w = u .- v

    critical_point = (
        (v[3] * sum(v .* w) - w[3] * sum(v .* v)) /
        (w[3] * sum(v .* w) - v[3] * sum(w .* w))
    )

    if !isfinite(critical_point)
        @warn "Problem computing critical point in theta_bound, got $critical_point"
    end

    a, b, c = greatcircleplane(domain)
    # Compute θ for the two endpoints
    θ = cmp(greatcircle(spherical(u).ϕ, a, b, c), greatcircle(spherical(v).ϕ, a, b, c))

    # If the critical point lies in the interval (0, 1) compute the
    # angle corresponding to it
    if !isnonpositive(critical_point) && !isnonnegative(critical_point - 1)
        θ = cmp(θ, greatcircle(spherical(v .+ critical_point .* w).ϕ, a, b, c))
    end

    return θ
end

function boundary_parameterization(t, domain::SphericalTriangle, i::Integer)
    i >= 1 && i <= 3 ||
        throw(ErrorException("attempt to use vertex number $i from a spherical triangle"))
    v = vertex(domain, mod1(i + 1, 3))
    w = vertex(domain, mod1(i + 2, 3))

    normalize(v .+ t .* (w - v))
end

"""
    boundary_points(domain::SphericalTriangle,
                    i::Integer,
                    n::Integer;
                    distribution = :linear,
                    half::Bool = true)
> Return n points on the boundary opposite of vertex number i.

  If `half` is true then only take points from the first half of the
  boundary.
"""
function boundary_points(
    domain::SphericalTriangle,
    i::Integer,
    n::Integer;
    distribution = :linear,
    half::Bool = false,
)
    i >= 1 && i <= 3 ||
        throw(ErrorException("attempt to use vertex number $i from a spherical triangle"))
    points = Vector{SVector{3,arb}}(undef, n)

    v = vertex(domain, mod1(i + 1, 3))
    w = vertex(domain, mod1(i + 2, 3))

    for j = 1:n
        t = ifelse(half, domain.parent(j // (2n + 1)), domain.parent(j // (n + 1)))

        points[j] = normalize(v .+ t .* (w - v))
    end

    return points, fill(i, n)
end

"""
    boundary_points(domain::SphericalTriangle,
                    eigenfunction::SphericalVertexEigenfunction,
                    n::Integer)
> Return n points on the boundary opposite of the vertex used for the
  eigenfunction.
"""
function boundary_points(
    domain::SphericalTriangle,
    eigenfunction::SphericalVertexEigenfunction,
    n::Integer,
)
    boundary_points(domain, eigenfunction.vertex, n, half = eigenfunction.stride == 2)
end

"""
    boundary_points(domain::SphericalTriangle,
                    eigenfunction::KrewerasEigenfunction,
                    n::Integer)
> Return n points from half of boundary one.
"""
function boundary_points(
    domain::SphericalTriangle,
    eigenfunction::KrewerasEigenfunction,
    n::Integer,
)
    boundary_points(domain, 1, n, half = true)
end

function interior_points(domain::SphericalTriangle, n::Integer; rng = MersenneTwister(42))
    # The points are generated by uniformly choosing a value for ϕ and
    # for θ we uniformly choose an u ∈ [-1, 1] and let θ = acos(1 -
    # u), this gives a uniform distribution on the sphere. To ensure
    # that the generated point is indeed inside the triangle we
    # restrict the possible values for ϕ, given by the angle at the
    # north pole. For θ we generate it and then check that it's inside
    # the triangle by checking that it's above the great circle
    # determining the lower boundary of the spherical triangle, if
    # it's outside we generate ϕ and θ again until we get something
    # which is inside.
    RR = domain.parent

    points = Vector{SVector{3,arb}}(undef, n)

    # Maximum value of ϕ
    A = angle(domain, 1)

    # Compute coefficients for plane determining the great circle
    # giving the lower boundary of the spherical triangle.
    (d, e, f) = greatcircleplane(domain)

    for i = 1:n
        # Randomly pick θ and ϕ
        θ = acos(1 - 2RR(rand(rng)))
        ϕ = rand(rng) * A

        # To determine if the point (θ, ϕ) is inside the triangle we
        # check if θ is above the θ value of great circle going
        # between the two vertices of the triangle not on the north
        # pole at the value ϕ.
        while θ > greatcircle(ϕ, d, e, f)
            θ = acos(1 - 2RR(rand(rng)))
            ϕ = rand(rng) * A
        end

        points[i] = cartesian(θ, ϕ)
    end

    points
end

"""
    anglesfromvertices(a, b, c)
> Compute the angles for the spherical triangle having the vertices
  `a`, `b` and `c`.
"""
function anglesfromvertices(a, b, c)
    v = cross(a, b)
    w = cross(a, c)
    α = acos(dot(v, w) / LinearAlgebra.norm(v) / LinearAlgebra.norm(w))

    v = cross(b, a)
    w = cross(b, c)
    β = acos(dot(v, w) / LinearAlgebra.norm(v) / LinearAlgebra.norm(w))

    v = cross(c, a)
    w = cross(c, b)
    γ = acos(dot(v, w) / LinearAlgebra.norm(v) / LinearAlgebra.norm(w))

    (α, β, γ)
end

"""
    subtriangle(domain::SphericalTriangle, ratio = 0.5)
> Return vertices for a spherical triangle in the interior of the
  domain. The size of the triangle is determined by `ratio`.
"""
function subtriangle(domain::SphericalTriangle; ratio = 0.5)
    v = center(domain)

    a = normalize((1 - ratio) * v + ratio * vertex(domain, 1))
    b = normalize((1 - ratio) * v + ratio * vertex(domain, 2))
    c = normalize((1 - ratio) * v + ratio * vertex(domain, 3))

    return (a, b, c)
end

"""
    partitiontriangle(a, b, c; iterations::Int = 1)
> Given the vertices of a spherical triangle return vertices for the
  four spherical triangle obtained by partitioning the original
  triangle into smaller ones.

  If `iterations` is not equal to 1 then recursively do this the give
  number of times, in case `iterations = 0` then return the original
  triangle.
"""
function partitiontriangle(a, b, c; iterations::Int = 1)
    if iterations < 1
        return [(a, b, c)]
    end
    x = normalize(a + b)
    y = normalize(b + c)
    z = normalize(c + a)

    triangles = [(a, x, z), (x, b, y), (z, y, c), (x, y, z)]
    if iterations == 1
        return triangles
    else
        return collect(
            Iterators.flatten([
                partitiontriangle(i, j, k, iterations = iterations - 1) for
                (i, j, k) in triangles
            ]),
        )
    end
end


text:
using Test
using MPI

MPI.Init()

try
    comm = MPI.COMM_WORLD
    size = MPI.Comm_size(comm)
    rank = MPI.Comm_rank(comm)

    # Tests with MPI.UNWEIGHTED
    sources      = Cint.(0:rank)
    destinations = Cint.(rank:(size-1))
    graph_comm = MPI.Dist_graph_create_adjacent(comm, sources, destinations)

    indegree, outdegree, weighted = MPI.Dist_graph_neighbors_count(graph_comm)
    @test indegree == (rank+1)
    @test outdegree == (size-rank)
    @test weighted == false

    incoming_edges = Vector{Cint}(undef,indegree)
    outgoing_edges = Vector{Cint}(undef,outdegree)
    MPI.Dist_graph_neighbors!(graph_comm, incoming_edges, outgoing_edges)

    @test sort(incoming_edges) == collect(0:rank)
    @test sort(outgoing_edges) == collect(rank:(size-1))

    # Tests with MPI.WEIGHTS_EMPTY
    sources      = Cint[]
    destinations = Cint[]
    graph_comm = MPI.Dist_graph_create_adjacent(comm, sources, destinations, source_weights=MPI.WEIGHTS_EMPTY,destination_weights=MPI.WEIGHTS_EMPTY)
    indegree, outdegree, weighted = MPI.Dist_graph_neighbors_count(graph_comm)
    @test indegree == 0
    @test outdegree == 0
    @test weighted == true
catch e
    if isa(e, MPI.FeatureLevelError)
        @test_broken e.min_version <= MPI.Get_version()
    else
        rethrow(e)
    end
end

MPI.Finalize()
@test MPI.Finalized()


text:
export make_continuous!

"""
`make_continuous!(vals,mod)` assumes that the list of numbers in `vals`
should represent a continuous stream of numbers, but with ambiguity modulo `mod`.
For example, they are the angles of a continuous list of complex values, but there
is an ambiguity modulo 2π. This function adjusts the values so they
appear continuous.

For finer control, use `make_continuous!(vals,mod,thresh)` where `thresh` is
the maximum allowable difference between consecutive entries in `val`.
"""
function make_continuous!(vals::Array{T,1}, mod::Real, thresh::Real) where T
    @assert mod>0 "mod must be positive, not $mod"
    @assert thresh>0 "threshold must be positive, not $thresh"
    @assert thresh<mod/2 "threshold must be less than half the modulus, not $thresh"

    nvals = length(vals)
    for j=2:nvals
        while abs(vals[j-1]-vals[j]) > thresh
            if vals[j] < vals[j-1]
                vals[j] += mod
            else
                vals[j] -= mod
            end
        end
    end
    nothing
end
make_continuous!(vals::Array{T,1}, mod::Real) where T = make_continuous!(vals, mod, mod/5)


text:
double_inc(x) = 2x + 1

@testset "rrule utils" begin

    rr = make_rrule(double_inc, 2.0)
    val, pb = rr(double_inc, 3.0)
    @test val == 7
    @test pb(1.0) == (ZeroTangent(), 2.0)

    config = YotaRuleConfig()
    val, pb = rrule_via_ad(config, double_inc, 3.0)
    @test val == 7
    @test pb(1.0) == (ZeroTangent(), 2.0)
end


text:
import Base
using DataStructures

"""Node of the state tree to explore"""
mutable struct AStarNode{TState, TCost <: Number}
  data::TState
  g::TCost
  f::TCost
  parent::Union{AStarNode{TState, TCost}, Nothing}
end

"order by f = g + h"
Base.isless(n1::AStarNode, n2::AStarNode) = Base.isless(n1.f, n2.f)

"Results structure"
struct AStarResult{TState, TCost <: Number}
  status::Symbol
  path::Vector{TState}
  cost::TCost
  closedsetsize::Int64
  opensetsize::Int64
end

function astar_compatibility_warn(; kwargs...)
  if "maxdepth" in keys(kwargs)
    @warn "'maxdepth' is no longer supported in the astar function, use 'maxcost' instead. Ignoring it..."
  end
end

mutable struct AStarSearchState{TState, TCost <: Number, THash}
  openheap::Vector{AStarNode{TState, TCost}}
  opennodedict::Dict{THash, AStarNode{TState, TCost}}
  closedset::Set{THash}
  start_time::Float64
  best_node::AStarNode{TState, TCost}
  best_heuristic::TCost

  function AStarSearchState(
    start_node::AStarNode{TState, TCost},
    start_hash::THash,
    start_heuristic::TCost,
  ) where {TState, TCost, THash}
    start_time = time()
    closedset = Set{THash}()
    openheap = [start_node]
    opennodedict = Dict(start_hash => start_node)
    return new{TState, TCost, THash}(
      openheap,
      opennodedict,
      closedset,
      start_time,
      start_node,
      start_heuristic,
    )
  end
end

function _astar!(
  astar_state::AStarSearchState{TState, TCost, THash},
  neighbours,
  goal,
  heuristic,
  cost,
  isgoal,
  hashfn,
  timeout,
  maxcost,
  enable_closedset,
) where {TState, TCost <: Number, THash}
  while !isempty(astar_state.openheap)
    node = heappop!(astar_state.openheap)

    if isgoal(node.data, goal)
      return AStarResult{TState, TCost}(
        :success,
        reconstructpath(node),
        node.g,
        length(astar_state.closedset),
        length(astar_state.openheap),
      )
    end

    nodehash = hashfn(node.data)
    delete!(astar_state.opennodedict, nodehash)

    if timeout < Inf && time() - astar_state.start_time >= timeout
      return AStarResult{TState, TCost}(
        :timeout,
        reconstructpath(astar_state.best_node),
        astar_state.best_node.g,
        length(astar_state.closedset),
        length(astar_state.openheap),
      )
    end

    if enable_closedset
      push!(astar_state.closedset, nodehash)
    end

    nodeheuristic = node.f - node.g
    if nodeheuristic < astar_state.best_heuristic
      astar_state.best_heuristic = nodeheuristic
      astar_state.best_node = node
    end

    neighbour_states = neighbours(node.data)

    for neighbour in neighbour_states
      neighbourhash = hashfn(neighbour)
      if enable_closedset && neighbourhash in astar_state.closedset
        continue
      end

      gfromthisnode = node.g + cost(node.data, neighbour)

      if gfromthisnode > maxcost
        continue
      end

      if neighbourhash in keys(astar_state.opennodedict)
        neighbournode = astar_state.opennodedict[neighbourhash]
        if gfromthisnode < neighbournode.g
          neighbourheuristic = neighbournode.f - neighbournode.g
          neighbournode.g = gfromthisnode
          neighbournode.f = gfromthisnode + neighbourheuristic
          neighbournode.parent = node
          heapify!(astar_state.openheap)
        end
      else
        neighbourheuristic = heuristic(neighbour, goal)
        neighbournode = AStarNode{TState, TCost}(
          neighbour,
          gfromthisnode,
          gfromthisnode + neighbourheuristic,
          node,
        )
        heappush!(astar_state.openheap, neighbournode)
        push!(astar_state.opennodedict, neighbourhash => neighbournode)
      end
    end
  end

  return AStarResult{TState, TCost}(
    :nopath,
    reconstructpath(astar_state.best_node),
    astar_state.best_node.g,
    length(astar_state.closedset),
    length(astar_state.openheap),
  )
end

"""
  astar(neighbours, start, goal;
        heuristic=defaultheuristic, cost=defaultcost, isgoal=defaultisgoal, hashfn=hash, timeout=Inf, maxcost=Inf)

Execute the A* algorithm to get the best path from the start state to reach a goal condition.
Only the first 3 arguments are mandatory, all the others are optional.

It returns a structure in which the `status` field is a Symbol that can be either:
- `:success`: the algorithm found a path from start to goal
- `:timeout`: the algorithm timed out, a partial path to the best state is returned in the `path` field
- `:nopath`: the algorithm didn't find any path to a goal, the path to the best state is still returned

The other fields are:
- `path`: an array of states from the start state to the goal or the best found state
- `cost`: the cost of the returned path
- `closedsetsize`: how many states the algorithm tested if they were a goal (size of the closed set)
- `opensetsize`: how many states were still in the open set when the algorithm ended

# Arguments
- `neighbours`: a function that takes a state and returns the neighbour states as an array (or iterable)
- `start`: the starting state, the type of the state is completely unrestricted
- `goal`: the goal state, the type is unrestricted, usually it's the same as the start
- `heuristic`: a function that given a state and the goal returns an estimate of the cost to reach goal. This estimate should be optimistic if you want to be sure to get the best path. Notice that the best path could be very expensive to find, so if you want a good but not guaranteed optimal path, you could multiply your heuristic by a constant, the algorithm will usually be much faster
- `cost`: a function that takes the current state and a neighbour and returns the cost to do that state transition. By default all transitions cost 1
- `isgoal`: a function that takes a state and the goal and evaluates if the goal is reached (by default ==)
- `hashfn`: a function that takes a state and returns a compact representation to use as dictionary key (usually one of UInt, Int, String), by default it is the base hash function. This is a very important field for composite states in order to avoid duplications. *WARNING* states with arrays as fields might return a different hash every time! If this is the case, please pass an hashfn that always returns the same value for the same state!
- `timeout`: timeout in number of seconds after which the algorithm stops returning the best partial path to the state with the lowest heuristic, by default it is unrestricted. Please notice that the algorithm wil run _AT LEAST_ the specified time.
- `maxcost`: a maximum bound of the accumulated cost of the path, this can result in a :nopath result even if a path to the goal (with a greater cost) exists. By default it is Inf
- `enable_closedset`: keep track of already visited nodes to avoid visiting them again, you might want to disable this if you know there isn't any loop in the state space graph (by default true)
"""
function astar(
  neighbours,
  start,
  goal;
  heuristic = defaultheuristic,
  cost = defaultcost,
  isgoal = defaultisgoal,
  hashfn = hash,
  timeout = Inf,
  maxcost = Inf,
  enable_closedset = true,
  kwargs...,
)
  astar_compatibility_warn(; kwargs...)

  start_heuristic = heuristic(start, goal)
  start_cost = zero(start_heuristic)
  start_node = AStarNode(start, start_cost, start_heuristic, nothing)
  start_hash = hashfn(start)

  astar_state = AStarSearchState(start_node, start_hash, start_heuristic)

  return _astar!(
    astar_state,
    neighbours,
    goal,
    heuristic,
    cost,
    isgoal,
    hashfn,
    timeout,
    maxcost,
    enable_closedset,
  )
end


text:
using FastMarching
using Test

function testMultipleStarts(tsize::Int=300, npoints::Int=20)
    SpeedImage = ones(tsize, tsize)
    SourcePoint = rand(2, npoints) .* tsize .+ 1

    t1 = FastMarching.msfm(SpeedImage, SourcePoint, true, true)
    t2 = FastMarching.msfm(SpeedImage, SourcePoint, true, true)

    @test norm(t1 - t2) < 1e-10 # Result is deterministic
end  # function testMultipleStarts

@testset "testMultipleStarts is deterministic" begin
    @testset "npoints=$npoints" for npoints in [10, 15]
        testMultipleStarts(300, npoints)
    end
end


text:
using DataFrames
using CSV
using StatsPlots

df = DataFrame(CSV.File("../mods/df_dac.csv"))
print(df)

remain_abs = df.sF + df.a2 - df.vAbs
remain_reg = df.sS + df.aR1 - df.vReg

groupedbar([remain_abs df.vAbs df.a1 remain_reg df.vReg], bar_position=:stack)
xlabel!("Hour")
savefig("bars.pdf")

max_stars_repo_name: dthierry/princetonDacLti

max_issues_repo_name: dthierry/princetonDacLti

max_forks_repo_name: dthierry/princetonDacLti


text:
#############################
# Partition type references #
#############################

const PartitionTypeReference = Tuple{ValueId,Integer}

############################
# Partitioning constraints #
############################

pt_ref_to_jl(pt_ref) =
    if pt_ref isa Tuple
        (convert(Future, pt_ref[1]).value_id, pt_ref[2] - 1)
    else
        (convert(Future, pt_ref).value_id, 0)
    end

pt_refs_to_jl(refs) =
    [pt_ref_to_jl(ref) for ref in refs]

struct PartitioningConstraintOverGroup
    type::String
    args::Vector{PartitionTypeReference}
end

struct PartitioningConstraintOverGroups
    type::String
    args::Vector{Vector{PartitionTypeReference}}
end

const PartitioningConstraint = Union{PartitioningConstraintOverGroup, PartitioningConstraintOverGroups}

to_jl(pc::PartitioningConstraint) = Dict("type" => pc.type, "args" => pc.args)

arg_to_jl_for_co(arg) =
    if arg isa Vector
        pt_refs_to_jl(arg)
    else
        [pt_ref_to_jl(arg)]
    end

function constraint_for_co(args)::PartitioningConstraintOverGroups
    if any(arg isa Vector for arg in args)
        args = [arg_to_jl_for_co(arg) for arg in args]
        PartitioningConstraintOverGroups("CO_GROUP", args)
    else
        PartitioningConstraintOverGroups("CO", pt_refs_to_jl(args))
    end
end

# TODO: Support Ordered
Co(args...) = constraint_for_co(args)
Cross(args...) = PartitioningConstraintOverGroup("CROSS", pt_refs_to_jl(args))
Equal(args...) = PartitioningConstraintOverGroup("EQUAL", pt_refs_to_jl(args))
Sequential(args...) =
    PartitioningConstraintOverGroup("SEQUENTIAL", pt_refs_to_jl(args))
# Note that a Match constraint will result in the same PT-level constraints
# being applied to all matching PTs
Match(args...) = PartitioningConstraintOverGroup("MATCH", pt_refs_to_jl(args))
MatchOn(on, args...) =
    PartitioningConstraintOverGroup(
        "MATCH_ON=" * string(on),
        pt_refs_to_jl(args),
    )
AtMost(npartitions, args...) =
    PartitioningConstraintOverGroup(
        "AT_MOST=$npartitions",
        pt_refs_to_jl(args)
    )
ScaleBy(arg, factor::Real = 1.0, relative_to...) =
    PartitioningConstraintOverGroup(
        "SCALE_BY=$factor",
        pt_refs_to_jl([arg; relative_to...])
    )

# Co, Cross, Equal, Sequential, Match, MatchOn, AtMost are PA-level constraints.
# AtMost, ScaleBy are PT-level constraints.
# For all constraints, contradictions result in the PA not being used. The
# exception is ScaleBy constraints where multiple PAs can specify a ScaleBy
# for the same PT and these are not fused in any way.

# TODO: Make the above constraint constructors produce dictionaries that have
# fields that make sense and are specialized for each one. This will reduce
# a significant amount of messy and hard-to-read code here (e.g., what in the
# world isa PartitioningConstraintOverGroup vs. a
# PartitioningConstraintOverGroups).

# NOTE: ScaleBy constraints accept PT references but only the values of PT
# references where the index is 1 are taken because we only scale relative
# to the memory usage that is split by the first PT in the PT compositions
# referenced

# NOTE: If you require a constraint for a particular PT, the onus is on you to
# ensure that whereever you use a value with that PT assigned, you always
# have the PTs that are referenced by the constraint. For example, if you use
# an AtMost constraint which references both PTs from a PT composition for a
# value where the first PT splits across workers and the second across batches,
# you need to ensure that anywhere you use the value, you actually do have a
# PT composition of length 2.

# NOTE: Currently, only AtMost and ScaleBy are supported as PT-level
# constraints (meaning they are included as part of a PT so that the PT
# cannot be applied to a variable unless the constraints are also be enforced)
# while ScaleBy may not be used as PA-level constraints (constraints that are
# applicable only for a single code region annotated with a PA)

mutable struct PartitioningConstraints
    constraints::Vector{Union{PartitioningConstraint, Function}}
end

PartitioningConstraints() = PartitioningConstraints([])

function to_jl(constraints::PartitioningConstraints)
    return Dict(
        "constraints" =>
            [to_jl(constraint) for constraint in constraints.constraints],
    )
end

###################
# Partition types #
###################

const PartitionTypeParameters = Dict{String, Any}

mutable struct PartitionType
    parameters::PartitionTypeParameters
    constraints::PartitioningConstraints

    PartitionType(
        parameters::Dict{String, <:Any} = PartitionTypeParameters(),
        constraints::PartitioningConstraints = PartitioningConstraints(),
    ) = new(parameters, constraints)
    PartitionType(s::String) = new(Dict("name" => s), PartitioningConstraints())
    PartitionType(parameters::PartitionTypeParameters) = new(parameters, PartitioningConstraints())

    function PartitionType(args::Union{String, Pair{String,<:Any}, PartitioningConstraint, Function}...)
        parameters = Dict()
        constraints = PartitioningConstraints()

        # Construct parameters and constraints from arguments
        for arg in args
            if arg isa String
                parameters["name"] = arg
            elseif arg isa Pair
                parameters[first(arg)] = last(arg)
            elseif arg isa PartitioningConstraint || arg isa Function
                push!(constraints.constraints, arg)
            else
                throw(ArgumentError("Expected either a partition type parameter or constraint"))
            end
        end

        new(parameters, constraints)
    end
end

# We probably need this so we can iterate over PTs produced by Grouped and then
# check the key property
function Base.getproperty(pt::PartitionType, name::Symbol)
    if hasfield(PartitionType, name)
        return getfield(pt, name)
    end

    n = string(name)
    if haskey(pt.parameters, n)
        return pt.parameters[n]
    end
    error("$name not found in partition type parameters")
end

function to_jl(pt::PartitionType)
    # Interpret bangs as random IDs
    for (k, v) in pt.parameters
        if v == "!"
            pt.parameters[k] = generate_bang_value()
        end
    end

    # Construct dictionary
    Dict("parameters" => pt.parameters, "constraints" => to_jl(pt.constraints))
end

##############################
# Partition type composition #
##############################

# This is mutable so that we can append PTs
mutable struct PartitionTypeComposition
    pts::Vector{PartitionType}
end

to_jl(ptc::PartitionTypeComposition) = [to_jl(pt) for pt in ptc.pts]

##############################
# Partition type combinators #
##############################

const PTOrPTUnion = Union{PartitionType,Vector{PartitionType}}

Base.:&(a::PartitionType, b::PartitionType) =
    if all(
        a.parameters[param_name] == b.parameters[param_name] for
        param_name in keys(a.parameters) if param_name in keys(b.parameters)
    )
        [PartitionType(
            merge(a.parameters, b.parameters),
            PartitioningConstraints(
                [a.constraints.constraints; b.constraints.constraints]
            )
        )]
    else
        PartitionType[]
    end

Base.:&(a::Vector{PartitionType}, b::PartitionType)::Vector{PartitionType} =
    vcat([pt & b for pt in a]...)
Base.:&(a::PartitionType, b::Vector{PartitionType}) = b & a
Base.:&(a::Vector{PartitionType}, b::Vector{PartitionType})::Vector{PartitionType} =
    vcat([aa & bb for aa in a for bb in b]...)
Base.:|(a::PTOrPTUnion, b::PTOrPTUnion) = [a; b]

#########################
# Partition annotations #
#########################

# TODO: Rename Partitions to PartitionTypeBinding and keep Partitioning as is
mutable struct Partitions
    # TODO: Only use either PT stack or PT composition to be consistent in
    # terminology
    pt_stacks::Dict{ValueId,PartitionTypeComposition}
end

Partitions() = Partitions(Dict())

function to_jl(p::Partitions)
    # NOTE: This assumes that the PT compositions in `p.pt_stacks` are _not_
    # delayed
    return Dict(
        "pt_stacks" =>
            Dict(v => ptc |> to_jl for (v, ptc) in p.pt_stacks),
    )
end

mutable struct PartitionAnnotation
    partitions::Partitions
    constraints::PartitioningConstraints
end

PartitionAnnotation() = PartitionAnnotation(Partitions(), PartitioningConstraints())

function to_jl(pa::PartitionAnnotation)
    return Dict(
        "partitions" => to_jl(pa.partitions),
        "constraints" => to_jl(pa.constraints),
    )
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
@testset "interval" begin
    @testset "radius" begin
        x = Arb(2)
        y = Arb(2)
        Arblib.set!(Arblib.radref(y), one(Mag))

        @test radius(x) == radius(Mag, x) == radius(Arf, x) == radius(Arb, x) == zero(Mag)
        @test radius(y) == radius(Mag, y) == radius(Arf, y) == radius(Arb, y) == one(Mag)

        @test radius(x) isa Mag
        @test radius(Mag, x) isa Mag
        @test radius(Arf, x) isa Arf
        @test radius(Arb, x) isa Arb

        @test precision(radius(Arf, Arb(prec = 80))) == 80
        @test precision(radius(Arb, Arb(prec = 80))) == 80
    end

    @testset "midpoint" begin
        x = zero(Arb)
        y = one(Arb)

        @test midpoint(x) == midpoint(Arf, x) == midpoint(Arb, x) == zero(Arf)
        @test midpoint(y) == midpoint(Arf, y) == midpoint(Arb, y) == one(Arf)
        @test midpoint(x) isa Arf
        @test midpoint(Arf, x) isa Arf
        @test midpoint(Arb, x) isa Arb

        @test precision(midpoint(Arb(prec = 80))) == 80
        @test precision(midpoint(Arf, Arb(prec = 80))) == 80
        @test precision(midpoint(Arb, Arb(prec = 80))) == 80
    end

    @testset "lbound/ubound" begin
        x = Arb(2)
        y = Arb(0)
        Arblib.set!(Arblib.radref(y), 1)

        @test lbound(x) == lbound(Arf, x) == lbound(Arb, x) == Arf(2)
        @test lbound(y) == lbound(Arf, y) == lbound(Arb, y) == Arf(-1)
        @test ubound(x) == ubound(Arf, x) == ubound(Arb, x) == Arf(2)
        @test ubound(y) == ubound(Arf, y) == ubound(Arb, y) == Arf(1)

        @test typeof(lbound(x)) == typeof(ubound(x)) == Arf
        @test typeof(lbound(Arf, x)) == typeof(ubound(Arf, x)) == Arf
        @test typeof(lbound(Arb, x)) == typeof(ubound(Arb, x)) == Arb

        @test precision(lbound(Arb(prec = 80))) == 80
        @test precision(lbound(Arf, Arb(prec = 80))) == 80
        @test precision(lbound(Arb, Arb(prec = 80))) == 80
        @test precision(ubound(Arb(prec = 80))) == 80
        @test precision(ubound(Arf, Arb(prec = 80))) == 80
        @test precision(ubound(Arb, Arb(prec = 80))) == 80
    end

    @testset "abs_lbound/abs_ubound" begin
        x = Arb(-2)
        y = Arb(1)
        Arblib.set!(Arblib.radref(x), 1)
        Arblib.set!(Arblib.radref(y), 2)
        xc = Acb(0, x)
        yc = Acb(0, y)

        @test abs_lbound(x) ==
              abs_lbound(Arf, x) ==
              abs_lbound(Arb, x) ==
              abs_lbound(xc) ==
              abs_lbound(Arf, xc) ==
              abs_lbound(Arb, xc) ==
              Arf(1)
        @test abs_lbound(y) ==
              abs_lbound(Arf, y) ==
              abs_lbound(Arb, y) ==
              abs_lbound(yc) ==
              abs_lbound(Arf, yc) ==
              abs_lbound(Arb, yc) ==
              Arf(0)
        @test abs_ubound(x) ==
              abs_ubound(Arf, x) ==
              abs_ubound(Arb, x) ==
              abs_ubound(xc) ==
              abs_ubound(Arf, xc) ==
              abs_ubound(Arb, xc) ==
              Arf(3)
        @test abs_ubound(y) ==
              abs_ubound(Arf, y) ==
              abs_ubound(Arb, y) ==
              abs_ubound(yc) ==
              abs_ubound(Arf, yc) ==
              abs_ubound(Arb, yc) ==
              Arf(3)

        @test typeof(abs_lbound(x)) ==
              typeof(abs_ubound(x)) ==
              typeof(abs_lbound(xc)) ==
              typeof(abs_ubound(xc)) ==
              Arf
        @test typeof(abs_lbound(Arf, x)) ==
              typeof(abs_ubound(Arf, x)) ==
              typeof(abs_lbound(Arf, xc)) ==
              typeof(abs_ubound(Arf, xc)) ==
              Arf
        @test typeof(abs_lbound(Arb, x)) ==
              typeof(abs_ubound(Arb, x)) ==
              typeof(abs_lbound(Arb, xc)) ==
              typeof(abs_ubound(Arb, xc)) ==
              Arb

        @test precision(abs_lbound(Arb(prec = 80))) ==
              precision(abs_lbound(Acb(prec = 80))) ==
              80
        @test precision(abs_lbound(Arf, Arb(prec = 80))) ==
              precision(abs_lbound(Arf, Acb(prec = 80))) ==
              80
        @test precision(abs_lbound(Arb, Arb(prec = 80))) ==
              precision(abs_lbound(Arb, Acb(prec = 80))) ==
              80
        @test precision(abs_ubound(Arb(prec = 80))) ==
              precision(abs_ubound(Acb(prec = 80))) ==
              80
        @test precision(abs_ubound(Arf, Arb(prec = 80))) ==
              precision(abs_ubound(Arf, Acb(prec = 80))) ==
              80
        @test precision(abs_ubound(Arb, Arb(prec = 80))) ==
              precision(abs_ubound(Arb, Acb(prec = 80))) ==
              80
    end

    @testset "getinterval" begin
        x = one(Arb)
        y = one(Arb)
        Arblib.set!(Arblib.radref(y), 1)

        @test getinterval(x) ==
              getinterval(Arf, x) ==
              getinterval(BigFloat, x) ==
              getinterval(Arb, x) ==
              (one(Arf), one(Arf))
        @test getinterval(y) ==
              getinterval(Arf, y) ==
              getinterval(BigFloat, y) ==
              getinterval(Arb, y) ==
              (Arf(0), Arf(2))

        @test getinterval(x) isa NTuple{2,Arf}
        @test getinterval(Arf, x) isa NTuple{2,Arf}
        @test getinterval(BigFloat, x) isa NTuple{2,BigFloat}
        @test getinterval(Arb, x) isa NTuple{2,Arb}

        @test precision.(getinterval(Arb(prec = 80))) == (80, 80)
        @test precision.(getinterval(Arf, Arb(prec = 80))) == (80, 80)
        @test precision.(getinterval(BigFloat, Arb(prec = 80))) == (80, 80)
        @test precision.(getinterval(Arb, Arb(prec = 80))) == (80, 80)
    end

    @testset "getball" begin
        x = one(Arb)
        y = one(Arb)
        Arblib.set!(Arblib.radref(y), 1)

        @test getball(x) == getball(Arb, x) == (one(Arf), zero(Mag))
        @test getball(y) == getball(Arb, y) == (one(Arf), one(Mag))

        @test getball(x) isa Tuple{Arf,Mag}
        @test getball(Arb, x) isa Tuple{Arb,Arb}

        @test precision(getball(Arb(prec = 80))[1]) == 80
        @test precision.(getball(Arb, Arb(prec = 80))) == (80, 80)
    end

    @testset "union" begin
        xs = [Arb(i) for i in vcat(1:10, 10:-1:1)]
        xsc = Acb.(xs)

        @test Arblib.contains(union(zero(Arb), one(Arb)), Arb(0.5))
        @test Arblib.contains(union(zero(Acb), Acb(1, 1)), Acb(0.5, 0.5))
        @test all(i -> Arblib.contains(union(xs...), i), 1:10)
        @test all(i -> Arblib.contains(union(xsc...), i), Acb.(1:10))

        @test precision(union(Arb(prec = 80), Arb(prec = 90))) == 90
        @test precision(union(Acb(prec = 80), Acb(prec = 90))) == 90
        @test precision(union([Arb(prec = p) for p = 70:10:100]...)) == 100
        @test precision(union([Acb(prec = p) for p = 70:10:100]...)) == 100
    end

    @testset "intersect" begin
        xs = [Arb((0, i)) for i in vcat(1:10, 10:-1:1)]

        @test Arblib.contains(intersect(Arb((0, 2)), Arb((1, 3))), Arb((1, 2)))
        @test Arblib.contains(intersect(xs...), Arb((0, 1)))
        @test !Arblib.contains(intersect(xs...), Arb(2))

        @test precision(intersect(Arb(prec = 80), Arb(prec = 90))) == 90
        @test precision(intersect([Arb(prec = p) for p = 70:10:100]...)) == 100

        @test_throws ArgumentError intersect(Arb(1), Arb(2))
        @test_throws ArgumentError intersect([xs; Arb(2)]...)
    end
end


text:
module PkgCacher

import Pkg.TOML
using Pkg
using Pkg.Types: UUID, PackageSpec, VersionSpec, SHA1

import LibGit2
import BinaryProvider

const WORKDIR = Ref(@__DIR__)
function setdir(path)
    WORKDIR[] = path
end

registries(path = depot_path()) = map(readdir(joinpath(path, "registries"))) do r
    abspath(path, "registries", r)
end
depot_path(path = WORKDIR[]) = joinpath(path, "depot")

function read_all_pkgs(filter = Set{UUID}())
    pkgs = Tuple{PackageSpec, SHA1}[]
    for path in registries()
        tomlpath = joinpath(path, "Registry.toml")
        isfile(tomlpath) || error("$path does not exist")
        registry = TOML.parsefile(tomlpath)
        haskey(registry, "packages") || error("Registry at $path has no packages")

        packages = registry["packages"]

        for (uuid, pkg) in packages
            name = pkg["name"]
            uuid = UUID(uuid)
            if uuid ∉ filter
                continue
            end
            name == "julia" && continue
            versions_file = abspath(path, pkg["path"], "Versions.toml")
            if !isfile(versions_file)
                @warn "$name has no Versions.toml file, skipping"
                continue
            end

            versions = TOML.parsefile(versions_file)

            for (version, hash) in versions
                push!(pkgs, (PackageSpec(name, uuid, VersionSpec(version)), SHA1(hash)))
            end
        end
    end
    pkgs
end

function read_project(path)
    pkgs = Tuple{PackageSpec, SHA1}[]
    isfile(path) || error("$path does not exist")
    project = TOML.parsefile(path)
    haskey(project, "deps") || error("Project at $path has no deps")

    uuids = Set{UUID}()
    for (name, uuid) in project["deps"]
        push!(uuids, UUID(uuid))
    end
    uuids
end

function cache(pkg)
    try
        # If a pkg got deleted this will hang.
        Pkg.add(pkg)
        # TODO: Does adding automatically rebuild?
        #       CUDAnative et.al could get stale and would need to rebuild
        #       regularly.
        Pkg.rm(pkg)
    catch
    end
end

function cache_all(filter_uuids)
    old_depot = copy(Base.DEPOT_PATH)
    empty!(Base.DEPOT_PATH)
    push!(Base.DEPOT_PATH, depot_path())

    ctx = Pkg.API.Context()
    Pkg.API.update_registry(ctx)
    BinaryProvider.probe_platform_engines!()

    pkgs=read_all_pkgs(filter_uuids)

    pkgs_to_install = Tuple{PackageSpec, SHA1, String}[]
    for (pkg, hash) in pkgs
        pkg.uuid in keys(ctx.stdlibs) && continue
        pkg.path == nothing || continue
        pkg.repo == nothing || continue
        path = Pkg.Operations.find_installed(pkg.name, pkg.uuid, hash)
        if !ispath(path)
            push!(pkgs_to_install, (pkg, hash, path))
        end
    end

    widths = [textwidth(pkg.name) for (pkg, _) in pkgs_to_install]
    max_name = length(widths) == 0 ? 0 : maximum(widths)

    ########################################
    # Install from archives asynchronously #
    ########################################
    jobs = Channel(ctx.num_concurrent_downloads);
    results = Channel(ctx.num_concurrent_downloads);
    @async begin
        for pkg in pkgs_to_install
            put!(jobs, pkg)
        end
    end

    for i in 1:ctx.num_concurrent_downloads
        @async begin
            for (pkg, hash, path) in jobs
                try
                    success = Pkg.Operations.install_archive(urls[pkg.uuid], hash, path)
                    if ctx.use_only_tarballs_for_downloads && !success
                        pkgerror("failed to get tarball from $(urls[pkg.uuid])")
                    end
                    put!(results, (pkg, success, path))
                catch err
                    put!(results, (pkg, err, catch_backtrace()))
                end
            end
        end
    end

    missed_packages = Tuple{PackageSpec, String}[]
    for i in 1:length(pkgs_to_install)
        pkg, exc_or_success, bt_or_path = take!(results)
        exc_or_success isa Exception && pkgerror("Error when installing packages:\n", sprint(Base.showerror, exc_or_success, bt_or_path))
        success, path = exc_or_success, bt_or_path
        if success
            vstr = pkg.version != nothing ? "v$(pkg.version)" : "[$h]"
            printpkgstyle(ctx, :Installed, string(rpad(pkg.name * " ", max_name + 2, "─"), " ", vstr))
        else
            push!(missed_packages, (pkg, path))
        end
    end

    ##################################################
    # Use LibGit2 to download any remaining packages #
    ##################################################
    for (pkg, path) in missed_packages
        uuid = pkg.uuid
        if !ctx.preview
            install_git(ctx, pkg.uuid, pkg.name, hashes[uuid], urls[uuid], pkg.version::VersionNumber, path)
        end
        vstr = pkg.version != nothing ? "v$(pkg.version)" : "[$h]"
        @info "Installed $(rpad(pkg.name * " ", max_name + 2, "─")) $vstr"
    end

    # Cleanup
    rm(joinpath(depot_path(), "environments"), recursive=true)

    empty!(Base.DEPOT_PATH)
    append!(Base.DEPOT_PATH, old_depot)
end

function install(target)
    rmdir(target, recursive=true)
    cp(depot_path(), target, recursive=true)
    rm(joinpath(target, "environments"), recursive=true)
    # Fix a Julia 1.0.0 bug that we attempt to update global registries
    for path in registries(target)
        rm(joinpath(path, ".git"), recursive=true)
    end
end

end

max_stars_repo_name: JuliaTagBot/PkgCacher

max_issues_repo_name: JuliaTagBot/PkgCacher

max_forks_repo_name: JuliaTagBot/PkgCacher


text:
seed!(RNG, 1337)
x = randn(RNG, 10000)

struct Scale{T}
    r::T
end
(f::Scale)(x) = f.r * x  # functor

@testset "default params" begin
    p = @inferred histogram(x)
    test_ref("references/histogram/default.txt", @print_col(p))
    test_ref("references/histogram/default_nocolor.txt", @print_nocol(p))
    p = @inferred histogram(x, closed = :left)
    test_ref("references/histogram/default.txt", @print_col(p))
    hist = fit(Histogram, x, closed = :left)
    p = @inferred histogram(hist)
    test_ref("references/histogram/default.txt", @print_col(p))
    p = @inferred histogram(x * 100)
    test_ref("references/histogram/default_1e2.txt", @print_col(p))
    p = @inferred histogram(x * 0.01)
    test_ref("references/histogram/default_1e-2.txt", @print_col(p))
    p = @inferred histogram(x, xscale = :log10)
    test_ref("references/histogram/log10.txt", @print_col(p))
    p = @inferred histogram(x, xscale = Scale(π))
    test_ref("references/histogram/functor.txt", @print_col(p))
    p = @inferred histogram(x, xlabel = "custom label", xscale = :log10)
    test_ref("references/histogram/log10_label.txt", @print_col(p))
    p = @inferred histogram([0.1f0, 0.1f0, 0.0f0])
    test_ref("references/histogram/float32.txt", @print_col(p))
    p = @inferred histogram(Histogram([0.0, 0.1, 1.0, 10.0, 100.0], [1, 2, 3, 4]))
    test_ref("references/histogram/nonuniformbins.txt", @print_col(p))
    x2 = copy(reshape(x, (1, length(x), 1, 1)))
    p = @inferred Plot histogram(x2)
    test_ref("references/histogram/default.txt", @print_col(p))
    p = @inferred histogram(Histogram(1:15.0, vcat(collect(1:13), 400)))
    test_ref("references/histogram/fraction.txt", @print_col(p))
end

@testset "hist params" begin
    hist = fit(Histogram, x, nbins = 5, closed = :right)
    p = @inferred histogram(hist)
    test_ref("references/histogram/hist_params.txt", @print_col(p))
    p = @inferred histogram(x, nbins = 5, closed = :right)
    test_ref("references/histogram/hist_params.txt", @print_col(p))
end

@testset "keyword arguments" begin
    p = @inferred histogram(
        x,
        title = "My Histogram",
        xlabel = "Absolute Frequency",
        color = :blue,
        margin = 7,
        padding = 3,
    )
    test_ref("references/histogram/parameters1.txt", @print_col(p))
    p = @inferred histogram(
        x,
        title = "My Histogram",
        xlabel = "Absolute Frequency",
        color = :blue,
        margin = 7,
        padding = 3,
        labels = false,
    )
    test_ref("references/histogram/parameters1_nolabels.txt", @print_col(p))

    p = @inferred histogram(
        x,
        title = "My Histogram",
        xlabel = "Absolute Frequency",
        color = :yellow,
        border = :solid,
        symbols = ["="],
        width = 50,
    )
    test_ref("references/histogram/parameters2.txt", @print_col(p))
    # same but with Char as symbols

    p = @inferred histogram(
        x,
        title = "My Histogram",
        xlabel = "Absolute Frequency",
        color = :yellow,
        border = :solid,
        symbols = ['='],
        width = 50,
    )
    test_ref("references/histogram/parameters2.txt", @print_col(p))

    # colors
    p = @inferred histogram(x, title = "Gray color", color = 240)
    test_ref("references/histogram/col1.txt", @print_col(p))

    p = @inferred histogram(x, title = "Green color", color = (0, 135, 95))
    test_ref("references/histogram/col2.txt", @print_col(p))
end


text:
const REPO_URL = "000000000000.dkr.ecr.us-east-1.amazonaws.com/foobar"
const RESPONSE = Response(400, "foobar")
const STATUS_ERROR = StatusError(400, "foo", "bar", RESPONSE)

describe_repository_exists_patch = @patch function JuliaOnLambda.ECR.describe_repositories(
    params
)
    return Dict("repositories" => [Dict("repositoryUri" => REPO_URL)])
end

describe_repository_dne_patch = @patch function JuliaOnLambda.ECR.describe_repositories(
    params
)
    return throw(
        AWSException(
            JuliaOnLambda.REPO_NOT_FOUND_EXCEPTION, "foobar", nothing, STATUS_ERROR, nothing
        ),
    )
end

create_repository_patch = @patch function JuliaOnLambda.ECR.create_repository(
    repository_name
)
    return Dict("repository" => Dict("repositoryUri" => REPO_URL))
end

create_repository_failure_patch = @patch function JuliaOnLambda.ECR.create_repository(
    repository_name
)
    return throw(
        AWSException("Alternative Error", "foobar", nothing, STATUS_ERROR, nothing)
    )
end


text:
using Documenter
# Pkg.activate("E:\\Programs\\blockpop\\TSSOS")
using TSSOS

makedocs(
# root = "E:\\Programs\\blockpop\\TSSOS\\docs",
sitename = "TSSOS",
pages = ["Home" => "index.md",
         "Polynomial Optimization" => "spop.md",
         "Noncommutative Polynomial Optimization" => "ncpop.md"],
modules = [TSSOS],
format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true")
)

deploydocs(
    repo = "github.com/wangjie212/TSSOS.git"
)