text:
function expand(v1,v2)
    v = vec([])
    for i in 1:size(v1)[1]
        for j in 1:size(v2)[1]
            push!(v, v1[i]*v2[j])
        end
    end
    return v
end

function calc(x)
    v = vec([])
    fin = vec([])
    for i in 1:size(x)[1]
        push!(v,join(sort(collect(x[i]))))
    end
    t = Set(v)
    temp = 0
    for j in t
        temp = count(i->(i == j),v)
        if temp != 1
            push!(fin,string(temp)*j)
        else
            push!(fin,j)
        end
    end
    return fin
end

function poly_expand(m,n)
    alpha = "abcdefghijklmnopqrstuvwxyz"
    x = vec([])
    for i in 1:m
        push!(x,alpha[i])
    end
    t = copy(x)
    for i in 1:(n-1)
        t = copy(expand(x,t))
    end
    return calc(t)
end

max_stars_repo_name: Arnav-Bista/Poly-Expansion

max_issues_repo_name: Arnav-Bista/Poly-Expansion

max_forks_repo_name: Arnav-Bista/Poly-Expansion


text:
@testset "FiniteFields $p" for p in [2, 3, 7]

    @test GF{p}(0) isa GF{p}
    @test_throws AssertionError GF{-p}(0)
    @test_throws AssertionError GF{p^2}(0)

    @test Int(zero(GF{p})) == 0
    @test Int(one(GF{p})) == 1
    @test characteristic(GF{p}) == p
    @test characteristic(zero(GF{p})) == p

    @test GF{p}(0) == zero(GF{p})
    @test zero(GF{p}) != one(GF{p})
    @test iszero(GF{p}(0))
    @test !isone(GF{p}(0))
    @test isone(GF{p}(1))
    @test !isone(GF{p}(2))

    @test GF{p}(2p-1) == GF{p}(p-1) == GF{p}(-1)

    @test hash(GF{p}(1)) == hash(GF{p}(p+1)) == hash(one(GF{p}))
    @test hash(GF{p}(0)) != hash(GF{p}(1))

    p != 23 && @test_throws DomainError GF{p}(1) == GF{23}(1)

    @test -GF{p}(23) == GF{p}(-(23 % p))

    @test inv(GF{p}(23)) isa GF{p}
    @test isone(GF{p}(23)*inv(GF{p}(23)))
    @test_throws DomainError inv(GF{p}(p))

    @test all(GF{p}(i) + GF{p}(j) == GF{p}(i + j) for i in 0:p-1 for j in 0:p-1)
    @test all(GF{p}(i) - GF{p}(j) == GF{p}(i - j) for i in 0:p-1 for j in 0:p-1)
    @test all(GF{p}(i) * GF{p}(j) == GF{p}(i * j) for i in 0:p-1 for j in 0:p-1)

    @test all(iszero, (-GF{p}(i) + GF{p}(i) for i in 0:p-1))
    @test all(isone, (GF{p}(i)*inv(GF{p}(i)) for i in 1:p-1))
    @test all(GF{p}(i)^4 == GF{p}(i^4) for i in 0:p-1)
    @test all(GF{p}(23)/GF{p}(i) isa GF{p} for i in 1:p-1)

    @test all(GF{p}(i) + j == GF{p}(i + j) for i in 0:p-1 for j in 0:p-1)
    @test all(GF{p}(i) - j == GF{p}(i - j) for i in 0:p-1 for j in 0:p-1)
    @test all(GF{p}(i) * j == GF{p}(i * j) for i in 0:p-1 for j in 0:p-1)
    @test all(i/GF{p}(j) == GF{p}(i)/GF{p}(j) for i in 0:p-1 for j in 1:p-1)

    p!= 23 && @test_throws DomainError GF{p}(1) + GF{23}(2)
    p!= 23 && @test_throws DomainError GF{p}(1) - GF{23}(2)
    p!= 23 && @test_throws DomainError GF{p}(1) * GF{23}(2)
    p!= 23 && @test_throws DomainError GF{p}(1) / GF{23}(2)

    for i in 0:p
        x = GF{p}(i)
        try
            w = sqrt(x)
            @test w*w == x
            @test FiniteFields.issquare(x)
        catch er
            if er isa DomainError
                @test FiniteFields.legendresymbol(Int(x), p) == -1
                @test !FiniteFields.issquare(x)
            else
                rethrow(er)
            end
        end
    end

    @test FiniteFields.generator(GF{p}) isa GF{p}
    g = FiniteFields.generator(GF{p})
    @test !iszero(g)
    @test sort(Int.(g^i for i in 1:p-1)) == 1:p-1

    @test collect(GF{p}) isa Vector{GF{p}}
    F = collect(GF{p})
    @test Int.(F) == collect(0:p-1)
    @test all(characteristic.(F) .== p)
    @test iszero(first(F))
    @test all(F .== F .+ F[1])
    @test all(F .* F[1] .== zeros(GF{p}, p))

    for f in F
        @test iszero(sum(f for _ in 1:characteristic(GF{p})))
        @test iszero(f*characteristic(GF{p}))
    end

    @test isone(F[2])
    @test all(F .== F .* F[2])

    @test string(one(GF{p})) == "1"*FiniteFields.subscriptify(p)
end

@testset "rootofunity" begin
    p = 31
    FFs = SymbolicWedderburn.FiniteFields
    F = FFs.GF{p}
    @test isone(FFs.generator(F)^30)
    @test isone(FFs.rootofunity(F, 5)^5)
    x = FFs.rootofunity(F, 5)
    @test !any(isone, (x^i for i in 1:4))
    @test_throws DomainError FFs.rootofunity(F, 7)
    y = FFs.rootofunity(F, 6)
    @test isone(y^6)
    @test !any(isone, (y^i for i in 1:5))
    @test !any(isone, ((x*y)^i for i in 1:29))
    @test isone((x*y)^30)
end

ext: jl


text:
using Juno
using Flux.Tracker: back!

tocb(f) = f
tocb(fs::AbstractVector) = () -> foreach(call, fs)

"""
    train!(loss, data, opt; cb = () -> ())

For each datapoint `d` in `data` computes the gradient of `loss(d...)` through
backpropagation and calls the optimizer `opt` and the callback `cb`
(i.e. `opt()` and `cb()`).

Multiple callbacks can be passed to `cb` as an array.
"""
function train!(loss, data, opt; cb = () -> ())
  cb = tocb(cb)
  @progress for d in data
    l = loss(d...)
    isinf(l.data[]) && error("Loss is Inf")
    isnan(l.data[]) && error("Loss is NaN")
    back!(l)
    opt()
    cb()
  end
end


text:
using PkgPage
using Test
import Franklin
const P = PkgPage
const F = Franklin

include("testutils.jl")

include("utils.jl")

include("coms/basic.jl")
include("coms/begin-end.jl")


text:
@testset "GMC_NS_Ensemble tests..." begin
    sample_dist=Normal(100.,5.)
    samples=rand(sample_dist, 10000)

    n_chains=10

    max_μ=1000
    min_μ=eps()
    max_λ=100
    min_λ=1e-4

    ln_max_μ=log(max_μ^2/sqrt(max_μ^2 + inv(max_λ)))
    ln_min_μ=log(min_μ^2/sqrt(min_μ^2 + inv(min_λ)))
    ln_max_λ=1/log(1+(inv(max_λ)/max_μ^2))
    ln_min_λ=1/log(1+(inv(min_λ)/min_μ^2))

    n_priors=[Uniform(min_μ,max_μ),Uniform(min_λ,max_λ)]
    n_box=[min_μ max_μ;min_λ max_λ]
    ln_priors=[Uniform(ln_min_μ,ln_max_μ),Uniform(ln_min_λ,ln_max_λ)]
    ln_box=[ln_min_μ ln_max_μ;ln_min_λ ln_max_λ]

    gmc=GMC_DEFAULTS
    gmc[1]=3
    gmc[2]=1e-5
    gmc[end]=1e4

    ne=Normal_Ensemble("NE_test", n_chains, samples, n_priors, n_box, gmc...)
    lne=LogNormal_Ensemble("LNE_test", n_chains, samples, ln_priors, ln_box, gmc...)

    uds=Vector{Vector{Function}}()

    lds=Vector{Vector{Function}}()

    ne_ev = converge_ensemble!(ne,backup=(true,100),upper_displays=uds, lower_displays=lds)
    lne_ev = converge_ensemble!(lne,backup=(true,100),upper_displays=uds, lower_displays=lds)

    @test ne_ev > lne_ev

    for e in [ne, lne]
        @info "Testing detailed balance at $(e.path)..."

        @test issorted(e.log_Li)
        @test all(e.log_Li.<0.)
        @test issorted(e.log_Xi,rev=true)
        @test all(e.log_Xi.<=0.)
        @test all(e.log_wi.<0.)
        @test all(e.log_Liwi.<0.)
        @test all(e.log_Zi.<0.)

        Nvec=Vector{Int64}()
        for i in 2:length(e.log_Li)
            Xi=e.log_Xi[i]
            N=round(-(i-1)/Xi)
            push!(Nvec, N)
        end

        @test all(e.GMC_Nmin.<=Nvec.<=n_chains)

        for chain in 1:e.t_counter-1
            ch_pfx=e.path*'/'*string(chain)*'.'
            total_spls=0
            while isfile(ch_pfx*string(total_spls+1))
                total_spls+=1
            end

            lastm=deserialize(ch_pfx*"1")

            for spl in 2:total_spls
                splm=deserialize(ch_pfx*string(spl))

                @test !(splm.θ!=lastm.θ && splm.v!=lastm.v)
                @test !(splm.θ==lastm.θ && splm.v==lastm.v)
                if splm.θ==lastm.θ
                    @test splm.v!=lastm.v
                    @test splm.log_Li==lastm.log_Li
                elseif splm.v==lastm.v
                    @test splm.θ!=lastm.θ
                    @test splm.log_Li >= lastm.log_Li
                else
                    @error "Galilean trajectory generation is broken!"
                end
                lastm=splm
            end
        end
    end

    rm("NE_test", recursive=true)
    rm("LNE_test", recursive=true)
end


text:
using Dierckx,LinearAlgebra

"Tristimulus values of digital color, based on spectral measurement and matching functions"
function matchcolors(C,λ,I,cmf;f=nothing)
    n = length(C);MC=Vector{Vector{Float64}}(undef,n)
    @views r = cmf[:,1]
    if isnothing(f)
        f=1;dr=1
    else
        dr = r[2]-r[1]
    end
    for i in 1:n
        @views s = Spline1D(λ[i],I[i],k=3,bc="zero")(r)'
        @views MC[i] = f*dr*[s*cmf[:,2],s*cmf[:,3],s*cmf[:,4]]
    end
    return MC,C
end
newcmf(cmf,m) = hcat(cmf[:,1],(m*cmf[:,2:end]')')
divsum(m)=m./sum(m,dims=1)
"Color matching Tristimulus values of single Wavelength unit spectral"
function matchlambda(λ,cmf)
    fs = [Spline1D(cmf[:,1],cmf[:,i+1],k=3,bc="zero") for i in 1:3]
    return map(f->f(λ),fs)
end
"Cone activations of digital colors, based on spectral measurement and cone fundamentals"
function lmsresponse(C,λ,I;observer=10)
    conef = observer==10 ? sscone10le : sscone2le
    matchcolors(C,λ,I,conef)
end
"""
CIE XYZ Tristimulus values of digital colors, based on following function:

XYZ = 683∫s(λ)x̄ȳz̄(λ)dλ

where s(λ) is the power spectrum, x̄ȳz̄(λ) are the matching functions.
Y is the luminance(cd/m²)
"""
function xyzresponse(C,λ,I;observer=10)
    if observer == 10
        conef = sscone10le
        m = LMSToXYZ10
    else
        conef = sscone2le
        m = LMSToXYZ2
    end
    matchcolors(C,λ,I,newcmf(conef,m),f=683)
end
"Michelson Contrast, where ``michelson(Lmax, Lmin) = weber(Lmax, Lmean), Lmean = (Lmax+Lmin)/2``"
contrast_michelson(Lmax,Lmin) = (Lmax-Lmin)/(Lmax+Lmin)
"Weber Contrast, where ``weber(L, Lb) = michelson(L, 2Lb-L)``"
contrast_weber(L,Lb) = (L/Lb)-1

TranslateXYZMatrix(;x=0,y=0,z=0) = TranslateXYZMatrix([x,y,z])
function TranslateXYZMatrix(t::Vector)
    tm =  [1.0 0.0 0.0   t[1];
           0.0 1.0 0.0   t[2];
           0.0 0.0 1.0   t[3];
           0.0 0.0 0.0 1.0]
end
function RotateXYZMatrix(rad::Real;dims=1)
    rads = [0.0,0.0,0.0]
    if 1<=dims<=3
        rads[dims]=rad
    end
    RotateXYZMatrix(rads)
end
function RotateXYZMatrix(rads::Vector)
    axis=nothing
    ri = findall(!iszero,rads)
    if length(ri) == 1 && ri[1] <= 3
        axis = ri[1]
    end
    isnothing(axis) && return Matrix{Float64}(I,4,4)
    c = cos(rads[axis])
    s = sin(rads[axis])
    if axis==1
      rm = [1.0 0.0 0.0 0.0;
            0.0   c  -s 0.0;
            0.0   s   c 0.0;
            0.0 0.0 0.0 1.0]
    elseif axis==2
      rm = [  c 0.0   s 0.0;
            0.0 1.0 0.0 0.0;
             -s 0.0   c 0.0;
            0.0 0.0 0.0 1.0]
    else
      rm = [  c  -s 0.0 0.0;
              s   c 0.0 0.0;
            0.0 0.0 1.0 0.0;
            0.0 0.0 0.0 1.0]
    end
end
quavectors(vs) = vcat(vs,ones(size(vs,2))')
trivectors(vs) = vs[1:3,:]
quamatrix(m) = vcat(hcat(m,zeros(3)),[0.0 0.0 0.0 1.0])
trimatrix(m) = m[1:3,1:3]
"Converting Matrix between digital RGB and LMS Space, based on spectral measurement and cone fundamentals"
function RGBLMSMatrix(C,λ,I;observer=10)
    TC,C = lmsresponse(C,λ,I,observer=observer)
    # Solve System of Linear Equations in form of Space Converting Matrix
    mc=hcat(C...);mtc=hcat(TC...)
    isqua = size(mc,1)==4
    if isqua
        mc = trivectors(mc)
    end
    RGBToLMS = mtc/mc
    LMSToRGB = inv(RGBToLMS)
    if isqua
        RGBToLMS = quamatrix(RGBToLMS)
        LMSToRGB = quamatrix(LMSToRGB)
    end
    return RGBToLMS,LMSToRGB
end
"Converting Matrix between digital RGB and CIE XYZ Space, based on spectral measurement and xyz matching functions"
function RGBXYZMatrix(C,λ,I;observer=10)
    TC,C = xyzresponse(C,λ,I,observer=observer)
    # Solve System of Linear Equations in form of Space Converting Matrix
    mc=hcat(C...);mtc=hcat(TC...)
    isqua = size(mc,1)==4
    if isqua
        mc = trivectors(mc)
    end
    RGBToXYZ = mtc/mc
    XYZToRGB = inv(RGBToXYZ)
    if isqua
        RGBToXYZ = quamatrix(RGBToXYZ)
        XYZToRGB = quamatrix(XYZToRGB)
    end
    return RGBToXYZ,XYZToRGB
end
"""
Converting Matrix between LMS and Cone Contrast(Weber) Space
(DH Brainard, Cone contrast and opponent modulation color spaces, human color vision, 1996)
"""
function LMSContrastMatrix(bg)
    # translate origin to bg to get differential cone activation
    t = [1 0 0 -bg[1];
         0 1 0 -bg[2];
         0 0 1 -bg[3];
         0 0 0 1]
    # scale relative to bg to get cone contrast
    s = [inv(bg[1]) 0 0 0;
         0 inv(bg[2]) 0 0;
         0 0 inv(bg[3]) 0;
         0 0 0 1]
    LMSToContrast = s*t
    ContrastToLMS = inv(LMSToContrast)
    return LMSToContrast,ContrastToLMS
end
"""
Converting Matrix between differential LMS and DKL[L+M, L-M, S-(L+M)] Space
(DH Brainard, Cone contrast and opponent modulation color spaces, human color vision, 1996)
"""
function dLMSDKLMatrix(bg;cone=nothing,v=nothing,isnorm=true)
    wl = 1;wm = 1
    if !isnothing(cone) && !isnothing(v)
        wl,wm = v[:,2]'/cone[:,2:3]'
    end
    dLMSToDKL = [wl  wm           0;
                1 -bg[1]/bg[2] 0;
               -wl -wm (wl*bg[1] + wm*bg[2])/bg[3]]
    if isnorm
        # Each column of the inverse of dLMSToDKL is the differential LMS relative to bg that isolating each DKL mechanism
        dlms_dkliso = inv(dLMSToDKL)
        # Cone contrast relative to bg
        cc = dlms_dkliso./bg[1:3]
        # Pooled cone contrast of each differential LMS relative to bg that isolating each DKL mechanism
        pcc = [norm(cc[:,i]) for i in 1:3]
        # Scale differential LMS vector by its pooled cone contrast
        udlms_dkliso = dlms_dkliso./pcc'
        # Rescale dLMSToDKL so that differential LMS which isolating DKL mechanism and having unit pooled cone contrast will result unit DKL response
        dLMSToDKL = inv(dLMSToDKL*udlms_dkliso)*dLMSToDKL
    end
    DKLTodLMS = inv(dLMSToDKL)
    return dLMSToDKL,DKLTodLMS
end
"Converting Matrix between LMS and DKL[L+M, L-M, S-(L+M)] Space"
function LMSDKLMatrix(bg;observer=10,isnorm=true)
    if observer == 10
        conef = sscone10le
        v = v10le
    else
        conef = sscone2le
        v = v2le
    end
    dLMSToDKL,DKLTodLMS = dLMSDKLMatrix(bg,cone=conef,v=v,isnorm=isnorm)
    t = [1 0 0 -bg[1];
         0 1 0 -bg[2];
         0 0 1 -bg[3];
         0 0 0 1]
    LMSToDKL = quamatrix(dLMSToDKL)*t
    DKLToLMS = inv(LMSToDKL)
    return LMSToDKL,DKLToLMS
end
"""
Intersection point of a line and a plane.
points of a line are defined as a direction(Dₗ) through a point(Pₗ): P = Pₗ + λDₗ , where λ is a scaler
points of a plane are defined as a plane through a point(Pₚ) and with normal vector(Nₚ) : Nₚᵀ(P - Pₚ) = 0 , where Nᵀ is the transpose of N

return point of intersection and if it's on direction
"""
function intersectlineplane(Pₗ,Dₗ,Pₚ,Nₚ)
    NₚᵀDₗ = Nₚ'*Dₗ
    NₚᵀDₗ == 0 && return nothing
    λ = Nₚ'*(Pₚ - Pₗ) / NₚᵀDₗ
    return Pₗ + λ*Dₗ, λ >=0
end
"""
Intersection point of a line and the six faces of the unit cube with origin as a vertex and three axies as edges.
points of a line are defined as a direction(Dₗ) through a point(Pₗ)

return intersection point on direction
"""
function intersectlineunitorigincube(Pₗ,Dₗ)
    ps = [zeros(3,3) ones(3,3)]
    ns = [Matrix{Float64}(I,3,3) Matrix{Float64}(I,3,3)]
    for i in 1:6
        p,d = intersectlineplane(Pₗ,Dₗ,ps[:,i],ns[:,i])
        if !isnothing(p) && all(i->-eps(1.0)<=i<=1+eps(1.0),p) && d
            return p
        end
    end
    return nothing
end
"""
Points of a line segment defined by two points P₀ and P₁
d: points density of unit line length
"""
function linepoints(P₀,P₁;d=10)
    d = max(1,d)
    Dₗ = P₁-P₀
    l =norm(Dₗ)
    hcat([P₀.+λ*Dₗ for λ in 0:1/(d*l):1]...)
end
function XYZ2xyY(m)
    xyz = divsum(trivectors(m))
    xyz[3,:]=m[2,:]
    replace!(xyz,NaN=>0.0)
    return xyz
end
function xyY2XYZ(m)
    s = m[3:3,:]./m[2:2,:]
    XYZ = s.*vcat(m[1:2,:], 1 .- sum(m[1:2,:],dims=1))
    replace!(XYZ,NaN=>0.0)
    return XYZ
end
function desaturate2gamut!(vs)
    for i in 1:size(vs,2)
        minc = minimum(vs[:,i])
        if minc<0
            vs[:,i].-=minc
        end
        maxc = maximum(vs[:,i])
        if maxc>0
            vs[:,i]./=maxc
        end
    end
    return vs
end
function cam16uniquehueindex(h′)
    h′ᵢ = zeros(Int,size(h′))
    for j in 1:length(h′)
        for i in 1:4
            if cam16uniquehue[1,i] <= h′[j] < cam16uniquehue[1,i+1]
                h′ᵢ[j]=i
                break
            end
        end
    end
    return h′ᵢ
end
function cam16uniquehuecomposition(H,i)
    Pl = round(Int,cam16uniquehue[3,i+1]-H)
    Pr = round(Int,H-cam16uniquehue[3,i])
    return (cam16uniquehue[4,i]=>Pl,cam16uniquehue[4,i+1]=>Pr)
end
function cam16nonlinearcompression(x,Fl;checksign=true)
    if checksign
        p = (Fl*abs(x)/100.0)^0.42
        sign(x)*400*p/(27.13 + p) + 0.1
    else
        p = (Fl*x/100.0)^0.42
        400*p/(27.13 + p) + 0.1
    end
end
"""
CAM16 Viewing Conditions

W: White in test illuminant [Xw, Yw, Zw]
Yb: Background in test conditions
La: Luminance of test adapting field (cd/m2)
Surround: Surround condition {Average, Dim, Dark}, Nc and F are functions of c, and their values can be linearly interpolated.
          Sr = Lsw / Ldw, where Lsw is the luminance of reference white in surround and Ldw in the display area.
          Sr == 0: Dark
          0 < Sr < 0.2: Dim
          Sr >= 0.2: Average
"""
function cam16view(;W=100*WP_D65,Surround=:Average,La=40,Yb=20)
    F,c,Nc = cam16surround[cam16surround[:,1].==Surround,2:end]
    k = 1.0/(5*La + 1)
    Fl = La*(k^4) + 0.1 * ((1 - k^4)^2) * ((5*La)^(1/3))

    Yw = W[2]
    n = Yb/Yw
    Nbb = Ncb = 0.725*(1/n)^0.2
    z = 1.48 + sqrt(n)
    D = clamp(F * (1 - (1/3.6)*exp(-(La+42)/92)), 0, 1)

    RGBw = cat16*W
    Da = Diagonal([D*Yw/RGBw[i] + 1 - D for i in 1:3])
    RGBwc = Da*RGBw
    RGBaw = cam16nonlinearcompression.(RGBwc,Fl,checksign=false)
    Aw = Nbb*(([2 1 1/20]*RGBaw)[1] - 0.305)

    return (c=c,Nc=Nc,Nbb=Nbb,Ncb=Ncb,Fl=Fl,n=n,z=z,Da=Da,Aw=Aw)
end
"""
Convert CIE XYZ to CIE Color Appearance Model 2016
Li, C., Li, Z., Wang, Z., Xu, Y., Luo, M.R., Cui, G., Melgosa, M., Brill, M.H., and Pointer, M. (2017). Comprehensive color solutions: CAM16, CAT16, and CAM16‐UCS.

return:
J is the lightness
C is the chroma
h is the hue angle
Q is the brightness
M is the colourfulness
s is the saturation
"""
function XYZ2CAM16(XYZ;Fl,Nc,Ncb,Nbb,n,c,z,Da,Aw)
    # Chromatic Adaptation
    RGBc = Da*cat16*XYZ
    # Non Linear Response Compression
    RGBa = cam16nonlinearcompression.(RGBc,Fl,checksign=true)
    # Perceptual Attribute Correlates
    a = ([1 -12/11 1/11]*RGBa)[:]
    b = ([1/9 1/9 -2/9]*RGBa)[:]
    h = mod.(2pi .+ atan.(b,a), 2pi)

    h′ = rad2deg.(h)
    i = h′ .< cam16uniquehue[1,1]
    h′[i] = h′[i] .+ 360

    h′ᵢ = cam16uniquehueindex(h′)
    H = cam16uniquehue[3,h′ᵢ] .+ (100*(h′.-cam16uniquehue[1,h′ᵢ])./cam16uniquehue[2,h′ᵢ]) ./ ( (h′.-cam16uniquehue[1,h′ᵢ])./cam16uniquehue[2,h′ᵢ] .+ (cam16uniquehue[1,h′ᵢ.+1] .- h′)./cam16uniquehue[2,h′ᵢ.+1] )
    Hc = cam16uniquehuecomposition.(H,h′ᵢ)

    et = 0.25 .* (cos.(deg2rad.(h′) .+ 2) .+ 3.8)
    t = (50000/13)*Nc*Ncb .* et .* sqrt.(a.^2 .+ b.^2) ./ (([1 1 21/20]*RGBa)[:])
    A = (Nbb.*([2 1 1/20]*RGBa .- 0.305))[:]

    J = 100 .* (A./Aw).^(c*z)
    Q = (4/c) .* sqrt.(J./100) .* (Aw + 4)*Fl^0.25
    C = t.^0.9 .* sqrt.(J./100) .* (1.64 - 0.29^n)^0.73
    M = (Fl^0.25).*C
    s = 100 .* sqrt.(M./Q)

    return (J=J,Q=Q,C=C,M=M,h=h,s=s,H=H,Hc=Hc)
end
"""
CAM16 Uniform Color Space in polar[J′, M′, h] or cartesian[J′, a′, b′]
Luo, M.R., Cui, G., and Li, C. (2006). Uniform colour spaces based on CIECAM02 colour appearance model. Color Research & Application 31, 320–330.
"""
function CAM16UCS(J,M,h;c1=0.007,c2=0.0228,form=:cartesian)
    J′ = (1 + 100c1).*J ./ (1 .+ c1.*J)
    M′ = log.(1 .+ c2.*M) ./ c2
    a′ = M′.*cos.(h)
    b′ = M′.*sin.(h)
    return form == :cartesian ? hcat(J′,a′,b′)' : hcat(J′,M′,h)'
end
function CAM16UCSinv(UCS;c1=0.007,c2=0.0228,form=:cartesian)
    J′ = UCS[1,:]
    if form == :cartesian
        M′ = sqrt.(UCS[2,:].^2 .+ UCS[3,:].^2)
        h = mod.(2pi .+ atan.(UCS[3,:],UCS[2,:]), 2pi)
    else
        M′ = UCS[2,:]
        h = UCS[3,:]
    end
    M = (exp.(c2.*M′) .- 1) ./ c2
    J = J′ ./ (1 .+ 100c1 .- c1.*J′)
    return (J=J,M=M,h=h)
end
function cam16nonlinearcompressioninv(x,Fl)
    p = x-0.1
    sign(p)*(100/Fl)*(27.13*abs(p) / (400-abs(p)))^(1/0.42)
end
"""
Convert CIE Color Appearance Model 2016 to CIE XYZ
"""
function CAM162XYZ(J,M,h;Fl,Nc,Ncb,Nbb,n,c,z,Da,Aw)
    C = M ./ (Fl^0.25)

    t = (C ./ (sqrt.(J./100) .* (1.64 - 0.29^n)^0.73)).^(1/0.9)
    et = 0.25 .* (cos.(h .+ 2) .+ 3.8)
    A = Aw.*(J./100).^(1/(c*z))
    p2 = A./Nbb .+ 0.305
    p3 = 21/20
    sh = sin.(h)
    ch = cos.(h)

    a=b=zeros(length(h))
    for i in 1:length(h)
        if t[i] !=0
            p1 = 50000/13 * Nc*Ncb * et[i] * (1 / t[i])
            if abs(sh[i]) >= abs(ch[i])
                p4 = p1/sh[i]
                b[i] = p2[i]*(2+p3)*(460/1403) / ( p4+(2+p3)*(220/1403)*(ch[i]/sh[i]) - (27/1403) + p3*(6300/1403) )
                a[i] = b[i]*ch[i]/sh[i]
            else
                p5 = p1/ch[i]
                a[i] = p2[i]*(2+p3)*(460/1403) / ( p5+(2+p3)*(220/1403) - (27/1403 - p3*6300/1403)*sh[i]/ch[i] )
                b[i] = a[i]*sh[i]/ch[i]
            end
        end
    end

    RGBa = [460/1403 451/1403 288/1403;
            460/1403 -891/1403 -261/1403;
            460/1403 -220/1403 -6300/1403]*([p2 a b]')
    RGBc = cam16nonlinearcompressioninv.(RGBa,Fl)
    XYZ = cat16\Da\RGBc
end

max_stars_repo_name: Experica/ColorLab

max_issues_repo_name: Experica/ColorLab

max_forks_repo_name: Experica/ColorLab


text:
module SingleFactorANOVA

using Distributions

export anova, AnovaResult, tukey_kramer, TukeyKramerResult

"""
Contains the result of calling `anova`.

# Fields
- SSB: Between treatment variation
- SSE: Error variation
- DFB: Between treatment degrees of freedom
- DFE: Error degrees of freedom
- MSB: Between treatment mean squares
- MSE: Error mean squares
- F: F-Statistic
- p: p-value for F-statistic

"""
struct AnovaResult
    SSB::Real
	SSE::Real
	DFB::Real
	DFE::Real
	MSB::Real
	MSE::Real
	F::Real
	p::Real
end

function ssb(X::AbstractArray)
    X̄ = mean(vcat(X...))
	nⱼ = [length(X[i]) for i in 1:length(X)]
	X̄ⱼ = [mean(X[i]) for i in 1:length(X)]
	sum(nⱼ .* (X̄ⱼ .- X̄).^2)
end

function sse(X::AbstractArray)
	X̄ⱼ = [mean(X[i]) for i in 1:length(X)]
	sse = []
	for j in 1:length(X)
		for i in 1:length(X[j])
			push!(sse, (X[j][i] - X̄ⱼ[j])^2)
		end
	end
	sum(sse)
end

"""
    anova(X::AbstractArray)

Perform ANOVA on `X`. `X` is an `AbstractArray` of `AbstractArray`s (typically a `Vector` of `Vectors`). Each vector
within `X` is treated as a different group and `anova` tests whether the means of the values contained in each vector
are equal or if at least two of them are different. In other words:

    H₀: μ₁ = μ₂... = μₖ
    H₁: means are not equal

where `k` is the length of `X`.

# Example
```julia-repl
julia> result = anova([[1, 2, 5, 9], [2, 6, 4, 2, 3, 8], [15, 6, 26]])
AnovaResult(303.44230769230774, 268.25, 2, 10, 151.72115384615387, 26.825, 5.655961000788588, 0.022745050729447377)
```

To get the p-value for the F statistic:

```julia-repl
julia> result.p
0.022745050729447377
```

To view the fields of `AnovaResult`:

```julia-repl
help?> AnovaResult
search: AnovaResult

  Contains the result of calling anova.

  Fields
  ≡≡≡≡≡≡≡≡

    •  SSB: Between treatment variation
    •  SSE: Error variation
    •  DFB: Between treatment degrees of freedom
    •  DFE: Error degrees of freedom
    •  MSB: Between treatment mean squares
    •  MSE: Error mean squares
    •  F: F-Statistic
    •  p: p-value for F-statistic
```
"""
function anova(X::AbstractArray)
	N = length(vcat(X...))
	k = length(X)
	SSB = ssb(X)
	SSE = sse(X)
	DFB = k - 1
	DFE = N - k
	MSB = SSB / DFB
	MSE = SSE / DFE
	F = MSB / MSE
	fdist = FDist(DFB,DFE)
	AnovaResult(
		SSB,
		SSE,
		DFB,
		DFE,
		MSB,
		MSE,
		F,
		ccdf(fdist, F)
	)
end

"""
Contains the result of calling `tukey-kramer`.

# Fields
- mean_difference: Pairwise differences between means
- q_crit: Pairwise critical q values
- significant: Boolean values indicating whether pairwise mean differences are significant
"""
struct TukeyKramerResult
	mean_differences::NamedTuple
	q_crit::NamedTuple
	significant::NamedTuple
end

"""
    tukey_kramer(X::AbstractArray, anova_result::AnovaResult, α=0.05)

Perform Tukey-Kramer test on `X`. `X` is the same `AbstractArray` that you passed to the `anova` function. This test
is only performed if you have rejected the null hypothesis after having called `anova`. It returns a `TukeyKramerResult` which
shows which pairs of means are significantly different from one another.

# Example
```julia-repl
julia> X = [[1, 2, 5, 9], [2, 6, 4, 2, 3, 8], [15, 6, 26]]
3-element Vector{Vector{Int64}}:
 [1, 2, 5, 9]
 [2, 6, 4, 2, 3, 8]
 [15, 6, 26]

julia> result = anova(X)
AnovaResult(303.44230769230774, 268.25, 2, 10, 151.72115384615387, 26.825, 5.655961000788588, 0.022745050729447377)

julia> julia> tk_result = tukey_kramer(X, result, 0.05)
TukeyKramerResult((|x1 - x2| = 0.08333333333333304, |x1 - x3| = 11.416666666666666, |x2 - x3| = 11.5), (|x1 - x2| = 9.164737679909036, |x1 - x3| = 10.843863861104227, |x2 - x3| = 10.03946712180748), (|x1 - x2| = false, |x1 - x3| = true, |x2 - x3| = true))
```

To view the fields of `TukeyKramerResult`:

```julia-repl
help?> TukeyKramerResult
search: TukeyKramerResult

  Contains the result of calling tukey-kramer.

  Fields
  ≡≡≡≡≡≡≡≡

    •  mean_difference: Pairwise differences between means

    •  q_crit: Pairwise critical q values

    •  significant: Boolean values indicating whether pairwise mean differences are significant
```

To just see the `significant` field:

```julia-repl
julia> tk_result.significant
(|x1 - x2| = false, |x1 - x3| = true, |x2 - x3| = true)
```
"""
function tukey_kramer(X, anova_result, α=0.05)
	x̄ = [abs(mean(X[i]) - mean(X[j])) for i in 1:length(X) for j in i+1:length(X)]
	d = StudentizedRange(anova_result.DFE, length(X))
	q_crit = [quantile(d, 1 - α) * sqrt((anova_result.MSE/2)*((1/length(X[i]))+(1/length(X[j])))) for i in 1:length(X) for j in i+1:length(X)]
	comparisons = Tuple(Symbol("|x$i - x$j|") for i in 1:length(X) for j in i+1:length(X))
	TukeyKramerResult(
		NamedTuple{comparisons}(x̄),
		NamedTuple{comparisons}(q_crit),
		NamedTuple{comparisons}(x̄ .> q_crit)
	)
end

end


text:
export Shamanskii

function Shamanskii(
    nlp_stop,
    PData::PDataFact,
    dₙ::Vector,
    gt::Vector;
    λfact2::Bool = false,
)
    # printstyled("on est dans Shamanskii \n", color = :bold)
    nlp_at_x = nlp_stop.current_state
    x = nlp_at_x.x
    n = length(x)
    xdemi = x + dₙ
    gtemp = grad(nlp_stop.pb, xdemi)

    ϵ2 = sqrt(eps(eltype(gt)))
    Γ = max.(abs.(PData.Δ), ϵ2)

    if (PData.λ == 0) && !(λfact2)
        d̃ = PData.L \ gtemp[PData.pp]
        d̂ = PData.L' \ (PData.Q * (PData.Q' * d̃ ./ Γ))
        d = -d̂[invperm(PData.pp)]
    else
        # println("on est ici 🐵↓")
        # @show (nlp_at_x.Hx + PData.λ * Matrix(1.0I, n, n))
        L, D, pp, ρ, ncomp = ldlt_symm((nlp_at_x.Hx + PData.λ * Matrix(1.0I, n, n)))
        # println("On a L")
        D = Hermitian(D)
        X = eigen(D)
        Δ = X.values
        Q = X.vectors
        d̃ = L \ gtemp[pp]
        d̂ = L' \ (Q * (Q' * d̃ ./ Γ))
        d = -d̂[invperm(pp)]
    end

    dₕₒ = dₙ .+ d

    return dₕₒ
end


text:
"""
    AbstractMillNode

Supertype for any structure representing a data node.
"""
abstract type AbstractMillNode end

"""
    AbstractProductNode <: AbstractMillNode

Supertype for any structure representing a data node implementing a Cartesian product of data in subtrees.
"""
abstract type AbstractProductNode <: AbstractMillNode end

"""
    AbstractBagNode <: AbstractMillNode

Supertype for any data node structure representing a multi-instance learning problem.
"""
abstract type AbstractBagNode <: AbstractMillNode end

"""
    Mill.data(n::AbstractMillNode)

Return data stored in node `n`.

# Examples
```jldoctest
julia> Mill.data(ArrayNode([1 2; 3 4], "metadata"))
2×2 Matrix{Int64}:
 1  2
 3  4

julia> Mill.data(BagNode(ArrayNode([1 2; 3 4]), bags([1:3, 4:4]), "metadata"))
2×2 ArrayNode{Matrix{Int64}, Nothing}:
 1  2
 3  4
```

See also: [`Mill.metadata`](@ref)
"""
data(n::AbstractMillNode) = getfield(n, :data)

"""
    Mill.metadata(n::AbstractMillNode)

Return metadata stored in node `n`.

# Examples
```jldoctest
julia> Mill.metadata(ArrayNode([1 2; 3 4], "metadata"))
"metadata"

julia> Mill.metadata(BagNode(ArrayNode([1 2; 3 4]), bags([1:3, 4:4]), "metadata"))
"metadata"
```

See also: [`Mill.data`](@ref)
"""
metadata(x::AbstractMillNode) = getfield(x, :metadata)

"""
    catobs(ns...)

Merge multiple nodes storing samples (observations) into one suitably promoting in the process if possible.

Similar to `Base.cat` but concatenates along the abstract \"axis\" where samples are stored.

In case of repeated calls with varying number of arguments or argument types, use `reduce(catobs, [ns...])`
to save compilation time.

# Examples
```jldoctest
julia> catobs(ArrayNode(zeros(2, 2)), ArrayNode([1 2; 3 4]))
2×4 ArrayNode{Matrix{Float64}, Nothing}:
 0.0  0.0  1.0  2.0
 0.0  0.0  3.0  4.0

julia> n = ProductNode((t1=ArrayNode(randn(2, 3)), t2=BagNode(ArrayNode(randn(3, 8)), bags([1:3, 4:5, 6:8]))))
ProductNode with 3 obs
  ├── t1: ArrayNode(2×3 Array with Float64 elements)
  └── t2: BagNode with 3 obs
            └── ArrayNode(3×8 Array with Float64 elements)

julia> catobs(n[1], n[3])
ProductNode with 2 obs
  ├── t1: ArrayNode(2×2 Array with Float64 elements)
  └── t2: BagNode with 2 obs
            └── ArrayNode(3×6 Array with Float64 elements)
```

See also: [`Mill.subset`](@ref).
"""
function catobs end

"""
    subset(n, i)

Extract a subset `i` of samples (observations) stored in node `n`.

Similar to `Base.getindex` or `MLDataPattern.getobs` but defined for all `Mill.jl` compatible data as well.

# Examples
```jldoctest
julia> Mill.subset(ArrayNode(NGramMatrix(["Hello", "world"])), 2)
2053×1 ArrayNode{NGramMatrix{String, Vector{String}, Int64}, Nothing}:
 "world"

julia> Mill.subset(BagNode(ArrayNode(randn(2, 8)), [1:2, 3:3, 4:7, 8:8]), 1:3)
BagNode with 3 obs
  └── ArrayNode(2×7 Array with Float64 elements)
```

See also: [`catobs`](@ref).
"""
function subset end

"""
    removeinstances(n::AbstractBagNode, mask)

Remove instances from `n` using `mask` and remap bag indices accordingly.

# Examples
```jldoctest
julia> b1 = BagNode(ArrayNode([1 2 3; 4 5 6]), bags([1:2, 0:-1, 3:3]))
BagNode with 3 obs
  └── ArrayNode(2×3 Array with Int64 elements)

julia> b2 = removeinstances(b1, [false, true, true])
BagNode with 3 obs
  └── ArrayNode(2×2 Array with Int64 elements)

julia> b2.data
2×2 ArrayNode{Matrix{Int64}, Nothing}:
 2  3
 5  6

julia> b2.bags
AlignedBags{Int64}(UnitRange{Int64}[1:1, 0:-1, 2:2])
```
"""
function removeinstances end

"""
    dropmeta(n:AbstractMillNode)

Drop metadata stored in data node `n`.

# Examples
```jldoctest
julia> n1 = ArrayNode(NGramMatrix(["foo", "bar"]), ["metafoo", "metabar"])
2053×2 ArrayNode{NGramMatrix{String, Vector{String}, Int64}, Vector{String}}:
 "foo"
 "bar"

julia> n2 = dropmeta(n1)
2053×2 ArrayNode{NGramMatrix{String, Vector{String}, Int64}, Nothing}:
 "foo"
 "bar"

julia> isnothing(Mill.metadata(n2))
true
```

See also: [`Mill.metadata`](@ref).
"""
function dropmeta end

"""
    mapdata(f, x)

Recursively apply `f` to data in all leaves of `x`.

# Examples
```jldoctest
julia> n1 = ProductNode((a=ArrayNode(zeros(2,2)), b=ArrayNode(ones(2,2))))
ProductNode with 2 obs
  ├── a: ArrayNode(2×2 Array with Float64 elements)
  └── b: ArrayNode(2×2 Array with Float64 elements)

julia> n2 = Mill.mapdata(x -> x .+ 1, n1)
ProductNode with 2 obs
  ├── a: ArrayNode(2×2 Array with Float64 elements)
  └── b: ArrayNode(2×2 Array with Float64 elements)

julia> Mill.data(n2).a
2×2 ArrayNode{Matrix{Float64}, Nothing}:
 1.0  1.0
 1.0  1.0

julia> Mill.data(n2).b
2×2 ArrayNode{Matrix{Float64}, Nothing}:
 2.0  2.0
 2.0  2.0
```
"""
mapdata(f, x) = f(x)

# functions to make datanodes compatible with getindex and with MLDataPattern
Base.getindex(x::T, i::BitArray{1}) where T <: AbstractMillNode = x[findall(i)]
Base.getindex(x::T, i::Vector{Bool}) where T <: AbstractMillNode = x[findall(i)]
Base.getindex(x::AbstractMillNode, i::Int) = x[i:i]
Base.lastindex(ds::AbstractMillNode) = nobs(ds)
MLDataPattern.getobs(x::AbstractMillNode, i) = x[i]
MLDataPattern.getobs(x::AbstractMillNode, i, ::ObsDim.Undefined) = x[i]
MLDataPattern.getobs(x::AbstractMillNode, i, ::ObsDim.Last) = x[i]

subset(x::AbstractMatrix, i) = x[:, i]
subset(x::AbstractVector, i) = x[i]
subset(x::AbstractMillNode, i) = x[i]
subset(x::DataFrame, i) = x[i, :]
subset(::Missing, i) = missing
subset(::Nothing, i) = nothing
subset(xs::Union{Tuple, NamedTuple}, i) = map(x -> x[i], xs)

include("arraynode.jl")

StatsBase.nobs(::Missing) = nothing

include("bagnode.jl")
include("weighted_bagnode.jl")
StatsBase.nobs(a::AbstractBagNode) = length(a.bags)
StatsBase.nobs(a::AbstractBagNode, ::Type{ObsDim.Last}) = nobs(a)
Base.ndims(x::AbstractBagNode) = Colon()

include("productnode.jl")

Base.ndims(x::AbstractProductNode) = Colon()
StatsBase.nobs(a::AbstractProductNode) = nobs(a.data[1], ObsDim.Last)
StatsBase.nobs(a::AbstractProductNode, ::Type{ObsDim.Last}) = nobs(a)

include("lazynode.jl")

catobs(as::Maybe{ArrayNode}...) = reduce(catobs, collect(as))
catobs(as::Maybe{BagNode}...) = reduce(catobs, collect(as))
catobs(as::Maybe{WeightedBagNode}...) = reduce(catobs, collect(as))
catobs(as::Maybe{ProductNode}...) = reduce(catobs, collect(as))
catobs(as::Maybe{LazyNode}...) = reduce(catobs, collect(as))

Base.cat(as::AbstractMillNode...; dims=:) = catobs(as...)

# reduction of common datatypes the way we like it
Base.reduce(::typeof(catobs), as::Vector{<:DataFrame}) = reduce(vcat, as)
Base.reduce(::typeof(catobs), as::Vector{<:AbstractMatrix}) = reduce(hcat, as)
Base.reduce(::typeof(catobs), as::Vector{<:AbstractVector}) = reduce(vcat, as)
Base.reduce(::typeof(catobs), as::Vector{Missing}) = missing
Base.reduce(::typeof(catobs), as::Vector{Nothing}) = nothing
Base.reduce(::typeof(catobs), as::Vector{Union{Missing, Nothing}}) = nothing

function Base.reduce(::typeof(catobs), as::Vector{Maybe{T}}) where T <: AbstractMillNode
    reduce(catobs, skipmissing(as) |> collect)
end

function Base.show(io::IO, @nospecialize(n::AbstractMillNode))
    print(io, nameof(typeof(n)))
    if !get(io, :compact, false)
        _show_data(IOContext(io, :compact => true), n)
        if !(n isa ArrayNode)
            print(io, " with ", nobs(n), " obs")
        end
    end
end

_show_data(io, n::LazyNode{Name}) where {Name} = print(io, "{", Name, "}")
_show_data(io, _) = print(io)


text:
# Packages required for testing
using Test
using Logging

# Package Under Test
using PGFMapPlots

# Set logging level
global_logger(SimpleLogger(stderr, Logging.Debug))

# Define package tests
@time @testset "PGFMapPlots Package Tests" begin
    testdir = joinpath(dirname(@__DIR__), "test")
    @time @testset "PGFMapPlots" begin
        include(joinpath(testdir, "test_plots.jl"))
    end
end


text:
# This file was generated, do not modify it. # hide
@load XGBoostClassifier
dtc = machine(XGBoostClassifier(), Xtrain, ytrain)
fit!(dtc)
ŷ = predict(dtc, Xtrain)
cross_entropy(ŷ, ytrain) |> mean


text:
using HPAT

@acc hpat function cumsum_test(file_name)
    df = DataSource(DataTable{:id=Int64,:x=Float64,:y=Float64}, HDF5, file_name)
    cs = cumsum(df[:x])
    return sum(cs)
end

cumsum_test(ENV["HOME"]*"/tmp/cumsum.hdf5")

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
using Gen
import Random
import MacroTools

@gen function model()
    n = @trace(geometric(0.5), :n)
    total = 0.
    for i=1:n
        total += @trace(normal(0, 1), (:x, i))
    end
    @trace(normal(total, 1.), :y)
    total
end

@gen function random_walk_proposal(trace, i::Int)
    @trace(normal(trace[(:x, i)], 0.1), (:x, i))
end

@gen function add_remove_proposal(trace)
    n = trace[:n]
    total = get_retval(trace)
    add = (n == 0) || @trace(bernoulli(0.5), :add)
    if add
        @trace(normal(trace[:y] - total, 1.), :new_x)
    end
    (n, add)
end

function add_remove_involution(trace, fwd_choices, ret, args)
    (n, add) = ret
    bwd_choices = choicemap()
    new_n = add ? n + 1 : n - 1
    constraints = choicemap((:n, new_n))
    if add
        bwd_choices[:add] = false
        constraints[(:x, new_n)] = fwd_choices[:new_x]
    else
        bwd_choices[:new_x] = trace[(:x, n)]
        (new_n > 0) && (bwd_choices[:add] = true)
    end
    new_trace, weight, = update(trace, (), (), constraints)
    (new_trace, bwd_choices, weight)
end

@pkern function k2(trace; check=false, observations=EmptyChoiceMap())
    mh(
        trace, add_remove_proposal, (), add_remove_involution,
        check=check, observations=observations)
end

@rkern k2 : k2

@pkern function k3(trace; check=false, observations=EmptyChoiceMap())
    perm = Random.randperm(trace[:n])
    constraints = choicemap()
    for (i, j) in enumerate(perm)
        constraints[(:x, i)] = trace[(:x, j)]
        constraints[(:x, j)] = trace[(:x, i)]
    end
    trace, = update(trace, (), (), constraints)
    trace, nothing # return trace and metadata..
end

@rkern k3 : k3

max_n_add_remove = 10 # to test that we have escaped the body of the kern properly

ex = quote
@kern function my_kernel(trace)

    # cycle through the x's and do a random walk update on each one
    for i in 1:trace[:n]
        trace ~ mh(trace, random_walk_proposal, (i,))
    end

    # repeatedly pick a random x and do a random walk update on it
    if trace[:n] > 0
        for rep in 1:10
            let i ~ uniform_discrete(1, trace[:n])
                trace ~ mh(trace, random_walk_proposal, (i,))
            end
        end
    end

    # remove the last x, or add a new one, a random number of times
    let n_add_remove_reps ~ uniform_discrete(0, max_n_add_remove)
        for rep in 1:n_add_remove_reps
            trace ~ k2(trace)
        end
    end

    # permute the x's
    trace ~ k3(trace)
end
end

println(MacroTools.striplines(macroexpand(Main, MacroTools.postwalk(MacroTools.unblock, ex))))
eval(ex)

function run_dsl_kernel(n::Int, iters::Int, check)
    obs = choicemap((:y, 10))
    trace, = generate(model, (), merge(obs, choicemap((:n, n))))
    for i=1:iters
        trace, = my_kernel(trace, check=check, observations=obs)
    end
end

function run_reversal_dsl_kernel(n::Int, iters::Int, check)
    obs = choicemap((:y, 10))
    trace, = generate(model, (), merge(obs, choicemap((:n, n))))
    for i=1:iters
        trace, = reversal(my_kernel)(trace, check=check, observations=obs)
    end
end

# checks disabled
@time run_dsl_kernel(100, 100, false)
@time run_dsl_kernel(100, 100, false)

# checks enabled
@time run_dsl_kernel(100, 100, true)
@time run_dsl_kernel(100, 100, true)

# run reverse DSL with checks enabled
@time run_reversal_dsl_kernel(100, 100, true)
@time run_reversal_dsl_kernel(100, 100, true)

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
### Day 3 - Problem 1
## Squares With Three Sides
## Author: Thanasis Georgiou

# Read triangle data
input = readdlm("triangles.csv")

# For each row
count = 0
for i in 1:size(input, 1)
    # Get a row and sort it
    triangle = sort(view(input, i, :))

    # Is this row a valid triangle?
    count += triangle[1] + triangle[2] > triangle[3]
end

# Print result
println("Possible Triangles: $(count)")


text:
julia> julia> Inf32 == Inf64 == Inf16 == Inf
true

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# string escaping & unescaping
cx = {
    0x00000000      '\0'        "\\0"
    0x00000001      '\x01'      "\\x01"
    0x00000006      '\x06'      "\\x06"
    0x00000007      '\a'        "\\a"
    0x00000008      '\b'        "\\b"
    0x00000009      '\t'        "\\t"
    0x0000000a      '\n'        "\\n"
    0x0000000b      '\v'        "\\v"
    0x0000000c      '\f'        "\\f"
    0x0000000d      '\r'        "\\r"
    0x0000000e      '\x0e'      "\\x0e"
    0x0000001a      '\x1a'      "\\x1a"
    0x0000001b      '\e'        "\\e"
    0x0000001c      '\x1c'      "\\x1c"
    0x0000001f      '\x1f'      "\\x1f"
    0x00000020      ' '         " "
    0x0000002f      '/'         "/"
    0x00000030      '0'         "0"
    0x00000039      '9'         "9"
    0x0000003a      ':'         ":"
    0x00000040      '@'         "@"
    0x00000041      'A'         "A"
    0x0000005a      'Z'         "Z"
    0x0000005b      '['         "["
    0x00000060      '`'         "`"
    0x00000061      'a'         "a"
    0x0000007a      'z'         "z"
    0x0000007b      '{'         "{"
    0x0000007e      '~'         "~"
    0x0000007f      '\x7f'      "\\x7f"
    0x000000bf      '\ubf'      "\\ubf"
    0x000000ff      '\uff'      "\\uff"
    0x00000100      '\u100'     "\\u100"
    0x000001ff      '\u1ff'     "\\u1ff"
    0x00000fff      '\ufff'     "\\ufff"
    0x00001000      '\u1000'    "\\u1000"
    0x00001fff      '\u1fff'    "\\u1fff"
    0x0000ffff      '\uffff'    "\\uffff"
    0x00010000      '\U10000'   "\\U10000"
    0x0001ffff      '\U1ffff'   "\\U1ffff"
    0x0002ffff      '\U2ffff'   "\\U2ffff"
    0x00030000      '\U30000'   "\\U30000"
    0x000dffff      '\Udffff'   "\\Udffff"
    0x000e0000      '\Ue0000'   "\\Ue0000"
    0x000effff      '\Ueffff'   "\\Ueffff"
    0x000f0000      '\Uf0000'   "\\Uf0000"
    0x000fffff      '\Ufffff'   "\\Ufffff"
    0x00100000      '\U100000'  "\\U100000"
    0x0010ffff      '\U10ffff'  "\\U10ffff"
}

for i = 1:size(cx,1)
    @test cx[i,1] == cx[i,2]
    @test string(cx[i,2]) == unescape_string(cx[i,3])
    if isascii(cx[i,2]) || !isprint(cx[i,2])
        @test cx[i,3] == escape_string(string(cx[i,2]))
    end
    for j = 1:size(cx,1)
        str = string(cx[i,2], cx[j,2])
        @test str == unescape_string(escape_string(str))
    end
end

for i = 0:0x7f, p = {"","\0","x","xxx","\x7f","\uFF","\uFFF",
                     "\uFFFF","\U10000","\U10FFF","\U10FFFF"}
    c = char(i)
    cp = string(c,p)
    op = string(char(div(i,8)), oct(i%8), p)
    hp = string(char(div(i,16)), hex(i%16), p)
    @test string(unescape_string(string("\\",oct(i,1),p))) == cp
    @test string(unescape_string(string("\\",oct(i,2),p))) == cp
    @test string(unescape_string(string("\\",oct(i,3),p))) == cp
    @test string(unescape_string(string("\\",oct(i,4),p))) == op
    @test string(unescape_string(string("\\x",hex(i,1),p))) == cp
    @test string(unescape_string(string("\\x",hex(i,2),p))) == cp
    @test string(unescape_string(string("\\x",hex(i,3),p))) == hp
end

@test "\z" == unescape_string("\z") == "z"
@test "\X" == unescape_string("\X") == "X"
@test "\AbC" == unescape_string("\AbC") == "AbC"

@test "\0" == unescape_string("\\0")
@test "\1" == unescape_string("\\1")
@test "\7" == unescape_string("\\7")
@test "\0x" == unescape_string("\\0x")
@test "\1x" == unescape_string("\\1x")
@test "\7x" == unescape_string("\\7x")
@test "\00" == unescape_string("\\00")
@test "\01" == unescape_string("\\01")
@test "\07" == unescape_string("\\07")
@test "\70" == unescape_string("\\70")
@test "\71" == unescape_string("\\71")
@test "\77" == unescape_string("\\77")
@test "\00x" == unescape_string("\\00x")
@test "\01x" == unescape_string("\\01x")
@test "\07x" == unescape_string("\\07x")
@test "\70x" == unescape_string("\\70x")
@test "\71x" == unescape_string("\\71x")
@test "\77x" == unescape_string("\\77x")
@test "\000" == unescape_string("\\000")
@test "\001" == unescape_string("\\001")
@test "\007" == unescape_string("\\007")
@test "\070" == unescape_string("\\070")
@test "\071" == unescape_string("\\071")
@test "\077" == unescape_string("\\077")
@test "\170" == unescape_string("\\170")
@test "\171" == unescape_string("\\171")
@test "\177" == unescape_string("\\177")
@test "\0001" == unescape_string("\\0001")
@test "\0011" == unescape_string("\\0011")
@test "\0071" == unescape_string("\\0071")
@test "\0701" == unescape_string("\\0701")
@test "\0711" == unescape_string("\\0711")
@test "\0771" == unescape_string("\\0771")
@test "\1701" == unescape_string("\\1701")
@test "\1711" == unescape_string("\\1711")
@test "\1771" == unescape_string("\\1771")

@test "\x0" == unescape_string("\\x0")
@test "\x1" == unescape_string("\\x1")
@test "\xf" == unescape_string("\\xf")
@test "\xF" == unescape_string("\\xF")
@test "\x0x" == unescape_string("\\x0x")
@test "\x1x" == unescape_string("\\x1x")
@test "\xfx" == unescape_string("\\xfx")
@test "\xFx" == unescape_string("\\xFx")
@test "\x00" == unescape_string("\\x00")
@test "\x01" == unescape_string("\\x01")
@test "\x0f" == unescape_string("\\x0f")
@test "\x0F" == unescape_string("\\x0F")

# integer parsing
@test is(parseint(Int32,"0",36),int32(0))
@test is(parseint(Int32,"1",36),int32(1))
@test is(parseint(Int32,"9",36),int32(9))
@test is(parseint(Int32,"A",36),int32(10))
@test is(parseint(Int32,"a",36),int32(10))
@test is(parseint(Int32,"B",36),int32(11))
@test is(parseint(Int32,"b",36),int32(11))
@test is(parseint(Int32,"F",36),int32(15))
@test is(parseint(Int32,"f",36),int32(15))
@test is(parseint(Int32,"Z",36),int32(35))
@test is(parseint(Int32,"z",36),int32(35))

@test parseint("0") == 0
@test parseint("-0") == 0
@test parseint("1") == 1
@test parseint("-1") == -1
@test parseint("9") == 9
@test parseint("-9") == -9
@test parseint("10") == 10
@test parseint("-10") == -10
@test parseint(Int64,"3830974272") == 3830974272
@test parseint(Int64,"-3830974272") == -3830974272
@test parseint('3') == 3
@test parseint('3', 8) == 3

parsebin(s) = parseint(s,2)
parseoct(s) = parseint(s,8)
parsehex(s) = parseint(s,16)

@test parsebin("0") == 0
@test parsebin("-0") == 0
@test parsebin("1") == 1
@test parsebin("-1") == -1
@test parsebin("10") == 2
@test parsebin("-10") == -2
@test parsebin("11") == 3
@test parsebin("-11") == -3
@test parsebin("1111000011110000111100001111") == 252645135
@test parsebin("-1111000011110000111100001111") == -252645135

@test parseoct("0") == 0
@test parseoct("-0") == 0
@test parseoct("1") == 1
@test parseoct("-1") == -1
@test parseoct("7") == 7
@test parseoct("-7") == -7
@test parseoct("10") == 8
@test parseoct("-10") == -8
@test parseoct("11") == 9
@test parseoct("-11") == -9
@test parseoct("72") == 58
@test parseoct("-72") == -58
@test parseoct("3172207320") == 434704080
@test parseoct("-3172207320") == -434704080

@test parsehex("0") == 0
@test parsehex("-0") == 0
@test parsehex("1") == 1
@test parsehex("-1") == -1
@test parsehex("9") == 9
@test parsehex("-9") == -9
@test parsehex("a") == 10
@test parsehex("-a") == -10
@test parsehex("f") == 15
@test parsehex("-f") == -15
@test parsehex("10") == 16
@test parsehex("-10") == -16
@test parsehex("0BADF00D") == 195948557
@test parsehex("-0BADF00D") == -195948557
@test parseint(Int64,"BADCAB1E",16) == 3135023902
@test parseint(Int64,"-BADCAB1E",16) == -3135023902
@test parseint(Int64,"CafeBabe",16) == 3405691582
@test parseint(Int64,"-CafeBabe",16) == -3405691582
@test parseint(Int64,"DeadBeef",16) == 3735928559
@test parseint(Int64,"-DeadBeef",16) == -3735928559

@test parseint("2\n") == 2
@test parseint("   2 \n ") == 2
@test parseint(" 2 ") == 2
@test parseint("2 ") == 2
@test parseint(" 2") == 2
@test parseint("+2\n") == 2
@test parseint("-2") == -2
@test_throws parseint("   2 \n 0")
@test_throws parseint("2x")
@test_throws parseint("-")

@test parseint("1234") == 1234
@test parseint("0x1234") == 0x1234
@test parseint("0o1234") == 0o1234
@test parseint("0b1011") == 0b1011
@test parseint("-1234") == -1234
@test parseint("-0x1234") == -int(0x1234)
@test parseint("-0o1234") == -int(0o1234)
@test parseint("-0b1011") == -int(0b1011)

for T in (Int8,Uint8,Int16,Uint16,Int32,Uint32,Int64,Uint64)#,Int128,Uint128) ## FIXME: #4905
    @test parseint(T,string(typemin(T))) == typemin(T)
    @test parseint(T,string(typemax(T))) == typemax(T)
    @test_throws parseint(T,string(big(typemin(T))-1))
    @test_throws parseint(T,string(big(typemax(T))+1))
end

# string manipulation
@test strip("\t  hi   \n") == "hi"

# some test strings
astr = "Hello, world.\n"
u8str = "∀ ε > 0, ∃ δ > 0: |x-y| < δ ⇒ |f(x)-f(y)| < ε"

# ascii search
for str in {astr, Base.GenericString(astr)}
    @test search(str, 'x') == 0
    @test search(str, '\0') == 0
    @test search(str, '\u80') == 0
    @test search(str, '∀') == 0
    @test search(str, 'H') == 1
    @test search(str, 'l') == 3
    @test search(str, 'l', 4) == 4
    @test search(str, 'l', 5) == 11
    @test search(str, 'l', 12) == 0
    @test search(str, ',') == 6
    @test search(str, ',', 7) == 0
    @test search(str, '\n') == 14
    @test search(str, '\n', 15) == 0
end

# ascii rsearch
for str in {astr}
    @test rsearch(str, 'x') == 0
    @test rsearch(str, '\0') == 0
    @test rsearch(str, '\u80') == 0
    @test rsearch(str, '∀') == 0
    @test rsearch(str, 'H') == 1
    @test rsearch(str, 'H', 0) == 0
    @test rsearch(str, 'l') == 11
    @test rsearch(str, 'l', 5) == 4
    @test rsearch(str, 'l', 4) == 4
    @test rsearch(str, 'l', 3) == 3
    @test rsearch(str, 'l', 2) == 0
    @test rsearch(str, ',') == 6
    @test rsearch(str, ',', 5) == 0
    @test rsearch(str, '\n') == 14
end

# utf-8 search
for str in {u8str, Base.GenericString(u8str)}
    @test search(str, 'z') == 0
    @test search(str, '\0') == 0
    @test search(str, '\u80') == 0
    @test search(str, '∄') == 0
    @test search(str, '∀') == 1
    @test search(str, '∀', 2) == 0
    @test search(str, '∃') == 13
    @test search(str, '∃', 14) == 0
    @test search(str, 'x') == 26
    @test search(str, 'x', 27) == 43
    @test search(str, 'x', 44) == 0
    @test search(str, 'δ') == 17
    @test search(str, 'δ', 18) == 33
    @test search(str, 'δ', 34) == 0
    @test search(str, 'ε') == 5
    @test search(str, 'ε', 6) == 54
    @test search(str, 'ε', 55) == 0
end

# utf-8 rsearch
for str in {u8str}
    @test rsearch(str, 'z') == 0
    @test rsearch(str, '\0') == 0
    @test rsearch(str, '\u80') == 0
    @test rsearch(str, '∄') == 0
    @test rsearch(str, '∀') == 1
    @test rsearch(str, '∀', 0) == 0
    @test rsearch(str, '∃') == 13
    @test rsearch(str, '∃', 14) == 13
    @test rsearch(str, '∃', 13) == 13
    @test rsearch(str, '∃', 12) == 0
    @test rsearch(str, 'x') == 43
    @test rsearch(str, 'x', 42) == 26
    @test rsearch(str, 'x', 25) == 0
    @test rsearch(str, 'δ') == 33
    @test rsearch(str, 'δ', 32) == 17
    @test rsearch(str, 'δ', 16) == 0
    @test rsearch(str, 'ε') == 54
    @test rsearch(str, 'ε', 53) == 5
    @test rsearch(str, 'ε', 4) == 0
end

# string search with a single-char string
@test search(astr, "x") == 0:-1
@test search(astr, "H") == 1:1
@test search(astr, "H", 2) == 0:-1
@test search(astr, "l") == 3:3
@test search(astr, "l", 4) == 4:4
@test search(astr, "l", 5) == 11:11
@test search(astr, "l", 12) == 0:-1
@test search(astr, "\n") == 14:14
@test search(astr, "\n", 15) == 0:-1

@test search(u8str, "z") == 0:-1
@test search(u8str, "∄") == 0:-1
@test search(u8str, "∀") == 1:1
@test search(u8str, "∀", 4) == 0:-1
@test search(u8str, "∃") == 13:13
@test search(u8str, "∃", 16) == 0:-1
@test search(u8str, "x") == 26:26
@test search(u8str, "x", 27) == 43:43
@test search(u8str, "x", 44) == 0:-1
@test search(u8str, "ε") == 5:5
@test search(u8str, "ε", 7) == 54:54
@test search(u8str, "ε", 56) == 0:-1

# string rsearch with a single-char string
@test rsearch(astr, "x") == 0:-1
@test rsearch(astr, "H") == 1:1
@test rsearch(astr, "H", 2) == 1:1
@test rsearch(astr, "H", 0) == 0:-1
@test rsearch(astr, "l") == 11:11
@test rsearch(astr, "l", 10) == 4:4
@test rsearch(astr, "l", 4) == 4:4
@test rsearch(astr, "l", 3) == 3:3
@test rsearch(astr, "l", 2) == 0:-1
@test rsearch(astr, "\n") == 14:14
@test rsearch(astr, "\n", 13) == 0:-1

@test rsearch(u8str, "z") == 0:-1
@test rsearch(u8str, "∄") == 0:-1
@test rsearch(u8str, "∀") == 1:1
@test rsearch(u8str, "∀", 0) == 0:-1
#TODO: setting the limit in the middle of a wide char
#      makes search fail but rsearch succeed.
#      Should rsearch fail as well?
#@test rsearch(u8str, "∀", 2) == 0:-1 # gives 1:3
@test rsearch(u8str, "∃") == 13:13
@test rsearch(u8str, "∃", 12) == 0:-1
@test rsearch(u8str, "x") == 43:43
@test rsearch(u8str, "x", 42) == 26:26
@test rsearch(u8str, "x", 25) == 0:-1
@test rsearch(u8str, "ε") == 54:54
@test rsearch(u8str, "ε", 53) == 5:5
@test rsearch(u8str, "ε", 4) == 0:-1

# string search with a single-char regex
@test search(astr, r"x") == 0:-1
@test search(astr, r"H") == 1:1
@test search(astr, r"H", 2) == 0:-1
@test search(astr, r"l") == 3:3
@test search(astr, r"l", 4) == 4:4
@test search(astr, r"l", 5) == 11:11
@test search(astr, r"l", 12) == 0:-1
@test search(astr, r"\n") == 14:14
@test search(astr, r"\n", 15) == 0:-1
@test search(u8str, r"z") == 0:-1
@test search(u8str, r"∄") == 0:-1
@test search(u8str, r"∀") == 1:1
@test search(u8str, r"∀", 4) == 0:-1
@test search(u8str, r"∀") == search(u8str, r"\u2200")
@test search(u8str, r"∀", 4) == search(u8str, r"\u2200", 4)
@test search(u8str, r"∃") == 13:13
@test search(u8str, r"∃", 16) == 0:-1
@test search(u8str, r"x") == 26:26
@test search(u8str, r"x", 27) == 43:43
@test search(u8str, r"x", 44) == 0:-1
@test search(u8str, r"ε") == 5:5
@test search(u8str, r"ε", 7) == 54:54
@test search(u8str, r"ε", 56) == 0:-1
for i = 1:endof(astr)
    @test search(astr, r"."s, i) == i:i
end
for i = 1:endof(u8str)
    if isvalid(u8str,i)
        @test search(u8str, r"."s, i) == i:i
    end
end

# string search with a zero-char string
for i = 1:endof(astr)
    @test search(astr, "", i) == i:i-1
end
for i = 1:endof(u8str)
    @test search(u8str, "", i) == i:i-1
end
@test search("", "") == 1:0

# string rsearch with a zero-char string
for i = 1:endof(astr)
    @test rsearch(astr, "", i) == i:i-1
end
for i = 1:endof(u8str)
    @test rsearch(u8str, "", i) == i:i-1
end
@test rsearch("", "") == 1:0

# string search with a zero-char regex
for i = 1:endof(astr)
    @test search(astr, r"", i) == i:i-1
end
for i = 1:endof(u8str)
    # TODO: should regex search fast-forward invalid indices?
    if isvalid(u8str,i)
        @test search(u8str, r""s, i) == i:i-1
    end
end

# string search with a two-char string literal
@test search("foo,bar,baz", "xx") == 0:-1
@test search("foo,bar,baz", "fo") == 1:2
@test search("foo,bar,baz", "fo", 3) == 0:-1
@test search("foo,bar,baz", "oo") == 2:3
@test search("foo,bar,baz", "oo", 4) == 0:-1
@test search("foo,bar,baz", "o,") == 3:4
@test search("foo,bar,baz", "o,", 5) == 0:-1
@test search("foo,bar,baz", ",b") == 4:5
@test search("foo,bar,baz", ",b", 6) == 8:9
@test search("foo,bar,baz", ",b", 10) == 0:-1
@test search("foo,bar,baz", "az") == 10:11
@test search("foo,bar,baz", "az", 12) == 0:-1

# string rsearch with a two-char string literal
@test rsearch("foo,bar,baz", "xx") == 0:-1
@test rsearch("foo,bar,baz", "fo") == 1:2
@test rsearch("foo,bar,baz", "fo", 1) == 0:-1
@test rsearch("foo,bar,baz", "oo") == 2:3
@test rsearch("foo,bar,baz", "oo", 2) == 0:-1
@test rsearch("foo,bar,baz", "o,") == 3:4
@test rsearch("foo,bar,baz", "o,", 1) == 0:-1
@test rsearch("foo,bar,baz", ",b") == 8:9
@test rsearch("foo,bar,baz", ",b", 6) == 4:5
@test rsearch("foo,bar,baz", ",b", 3) == 0:-1
@test rsearch("foo,bar,baz", "az") == 10:11
@test rsearch("foo,bar,baz", "az", 10) == 0:-1

# array rsearch
@test rsearch(Uint8[1,2,3],Uint8[2,3],3) == 2:3
@test rsearch(Uint8[1,2,3],Uint8[2,3],1) == 0:-1

# string search with a two-char regex
@test search("foo,bar,baz", r"xx") == 0:-1
@test search("foo,bar,baz", r"fo") == 1:2
@test search("foo,bar,baz", r"fo", 3) == 0:-1
@test search("foo,bar,baz", r"oo") == 2:3
@test search("foo,bar,baz", r"oo", 4) == 0:-1
@test search("foo,bar,baz", r"o,") == 3:4
@test search("foo,bar,baz", r"o,", 5) == 0:-1
@test search("foo,bar,baz", r",b") == 4:5
@test search("foo,bar,baz", r",b", 6) == 8:9
@test search("foo,bar,baz", r",b", 10) == 0:-1
@test search("foo,bar,baz", r"az") == 10:11
@test search("foo,bar,baz", r"az", 12) == 0:-1

# split
@test isequal(split("foo,bar,baz", 'x'), ["foo,bar,baz"])
@test isequal(split("foo,bar,baz", ','), ["foo","bar","baz"])
@test isequal(split("foo,bar,baz", ","), ["foo","bar","baz"])
@test isequal(split("foo,bar,baz", r","), ["foo","bar","baz"])
@test isequal(split("foo,bar,baz", ',', 0), ["foo","bar","baz"])
@test isequal(split("foo,bar,baz", ',', 1), ["foo,bar,baz"])
@test isequal(split("foo,bar,baz", ',', 2), ["foo","bar,baz"])
@test isequal(split("foo,bar,baz", ',', 3), ["foo","bar","baz"])
@test isequal(split("foo,bar", "o,b"), ["fo","ar"])

@test isequal(split("", ','), [""])
@test isequal(split(",", ','), ["",""])
@test isequal(split(",,", ','), ["","",""])
@test isequal(split("", ',', false), [])
@test isequal(split(",", ',', false), [])
@test isequal(split(",,", ',', false), [])

@test isequal(split("a b c"), ["a","b","c"])
@test isequal(split("a  b \t c\n"), ["a","b","c"])

@test isequal(rsplit("foo,bar,baz", 'x'), ["foo,bar,baz"])
@test isequal(rsplit("foo,bar,baz", ','), ["foo","bar","baz"])
@test isequal(rsplit("foo,bar,baz", ","), ["foo","bar","baz"])
@test isequal(rsplit("foo,bar,baz", ',', 0), ["foo","bar","baz"])
@test isequal(rsplit("foo,bar,baz", ',', 1), ["foo,bar,baz"])
@test isequal(rsplit("foo,bar,baz", ',', 2), ["foo,bar","baz"])
@test isequal(rsplit("foo,bar,baz", ',', 3), ["foo","bar","baz"])
@test isequal(rsplit("foo,bar", "o,b"), ["fo","ar"])

@test isequal(rsplit("", ','), [""])
@test isequal(rsplit(",", ','), ["",""])
@test isequal(rsplit(",,", ','), ["","",""])
@test isequal(rsplit(",,", ',', 2), [",",""])
@test isequal(rsplit("", ',', false), [])
@test isequal(rsplit(",", ',', false), [])
@test isequal(rsplit(",,", ',', false), [])

#@test isequal(rsplit("a b c"), ["a","b","c"])
#@test isequal(rsplit("a  b \t c\n"), ["a","b","c"])

let str = "a.:.ba..:..cba.:.:.dcba.:."
@test isequal(split(str, ".:."), ["a","ba.",".cba",":.dcba",""])
@test isequal(split(str, ".:.", false), ["a","ba.",".cba",":.dcba"])
@test isequal(split(str, ".:."), ["a","ba.",".cba",":.dcba",""])
@test isequal(split(str, r"\.(:\.)+"), ["a","ba.",".cba","dcba",""])
@test isequal(split(str, r"\.(:\.)+", false), ["a","ba.",".cba","dcba"])
@test isequal(split(str, r"\.+:\.+"), ["a","ba","cba",":.dcba",""])
@test isequal(split(str, r"\.+:\.+", false), ["a","ba","cba",":.dcba"])

@test isequal(rsplit(str, ".:."), ["a","ba.",".cba.:","dcba",""])
@test isequal(rsplit(str, ".:.", false), ["a","ba.",".cba.:","dcba"])
@test isequal(rsplit(str, ".:.", 2), ["a.:.ba..:..cba.:.:.dcba", ""])
@test isequal(rsplit(str, ".:.", 3), ["a.:.ba..:..cba.:", "dcba", ""])
@test isequal(rsplit(str, ".:.", 4), ["a.:.ba.", ".cba.:", "dcba", ""])
@test isequal(rsplit(str, ".:.", 5), ["a", "ba.", ".cba.:", "dcba", ""])
@test isequal(rsplit(str, ".:.", 6), ["a", "ba.", ".cba.:", "dcba", ""])
end

# zero-width splits
@test isequal(rsplit("", ""), [""])

@test isequal(split("", ""), [""])
@test isequal(split("", r""), [""])
@test isequal(split("abc", ""), ["a","b","c"])
@test isequal(split("abc", r""), ["a","b","c"])
@test isequal(split("abcd", r"b?"), ["a","c","d"])
@test isequal(split("abcd", r"b*"), ["a","c","d"])
@test isequal(split("abcd", r"b+"), ["a","cd"])
@test isequal(split("abcd", r"b?c?"), ["a","d"])
@test isequal(split("abcd", r"[bc]?"), ["a","","d"])
@test isequal(split("abcd", r"a*"), ["","b","c","d"])
@test isequal(split("abcd", r"a+"), ["","bcd"])
@test isequal(split("abcd", r"d*"), ["a","b","c",""])
@test isequal(split("abcd", r"d+"), ["abc",""])
@test isequal(split("abcd", r"[ad]?"), ["","b","c",""])

# replace
@test replace("foobar", 'o', '0') == "f00bar"
@test replace("foobar", 'o', '0', 1) == "f0obar"
@test replace("foobar", 'o', "") == "fbar"
@test replace("foobar", 'o', "", 1) == "fobar"
@test replace("foobar", 'f', 'F') == "Foobar"
@test replace("foobar", 'r', 'R') == "foobaR"

@test replace("foofoofoo", "foo", "bar") == "barbarbar"
@test replace("foobarfoo", "foo", "baz") == "bazbarbaz"
@test replace("barfoofoo", "foo", "baz") == "barbazbaz"

@test replace("", "", "") == ""
@test replace("", "", "x") == "x"
@test replace("", "x", "y") == ""

@test replace("abcd", "", "^") == "^a^b^c^d^"
@test replace("abcd", "b", "^") == "a^cd"
@test replace("abcd", r"b?", "^") == "^a^c^d^"
@test replace("abcd", r"b+", "^") == "a^cd"
@test replace("abcd", r"b?c?", "^") == "^a^d^"
@test replace("abcd", r"[bc]?", "^") == "^a^^d^"

@test replace("foobarfoo", r"(fo|ba)", "xx") == "xxoxxrxxo"
@test replace("foobarfoo", r"(foo|ba)", "bar") == "barbarrbar"

@test replace("foobar", 'o', 'ø') == "føøbar"
@test replace("foobar", 'o', 'ø', 1) == "føobar"
@test replace("føøbar", 'ø', 'o') == "foobar"
@test replace("føøbar", 'ø', 'o', 1) == "foøbar"
@test replace("føøbar", 'ø', 'ö') == "fööbar"
@test replace("føøbar", 'ø', 'ö', 1) == "föøbar"
@test replace("føøbar", 'ø', "") == "fbar"
@test replace("føøbar", 'ø', "", 1) == "føbar"
@test replace("føøbar", 'f', 'F') == "Føøbar"
@test replace("ḟøøbar", 'ḟ', 'F') == "Føøbar"
@test replace("føøbar", 'f', 'Ḟ') == "Ḟøøbar"
@test replace("ḟøøbar", 'ḟ', 'Ḟ') == "Ḟøøbar"
@test replace("føøbar", 'r', 'R') == "føøbaR"
@test replace("føøbaṙ", 'ṙ', 'R') == "føøbaR"
@test replace("føøbar", 'r', 'Ṙ') == "føøbaṘ"
@test replace("føøbaṙ", 'ṙ', 'Ṙ') == "føøbaṘ"

@test replace("ḟøøḟøøḟøø", "ḟøø", "bar") == "barbarbar"
@test replace("ḟøøbarḟøø", "ḟøø", "baz") == "bazbarbaz"
@test replace("barḟøøḟøø", "ḟøø", "baz") == "barbazbaz"

@test replace("foofoofoo", "foo", "ƀäṙ") == "ƀäṙƀäṙƀäṙ"
@test replace("fooƀäṙfoo", "foo", "baz") == "bazƀäṙbaz"
@test replace("ƀäṙfoofoo", "foo", "baz") == "ƀäṙbazbaz"

@test replace("foofoofoo", "foo", "bar") == "barbarbar"
@test replace("foobarfoo", "foo", "ƀäż") == "ƀäżbarƀäż"
@test replace("barfoofoo", "foo", "ƀäż") == "barƀäżƀäż"

@test replace("ḟøøḟøøḟøø", "ḟøø", "ƀäṙ") == "ƀäṙƀäṙƀäṙ"
@test replace("ḟøøƀäṙḟøø", "ḟøø", "baz") == "bazƀäṙbaz"
@test replace("ƀäṙḟøøḟøø", "ḟøø", "baz") == "ƀäṙbazbaz"

@test replace("ḟøøḟøøḟøø", "ḟøø", "bar") == "barbarbar"
@test replace("ḟøøbarḟøø", "ḟøø", "ƀäż") == "ƀäżbarƀäż"
@test replace("barḟøøḟøø", "ḟøø", "ƀäż") == "barƀäżƀäż"

@test replace("ḟøøḟøøḟøø", "ḟøø", "ƀäṙ") == "ƀäṙƀäṙƀäṙ"
@test replace("ḟøøƀäṙḟøø", "ḟøø", "ƀäż") == "ƀäżƀäṙƀäż"
@test replace("ƀäṙḟøøḟøø", "ḟøø", "ƀäż") == "ƀäṙƀäżƀäż"

@test replace("", "", "ẍ") == "ẍ"
@test replace("", "ẍ", "ÿ") == ""

@test replace("äƀçđ", "", "π") == "πäπƀπçπđπ"
@test replace("äƀçđ", "ƀ", "π") == "äπçđ"
@test replace("äƀçđ", r"ƀ?", "π") == "πäπçπđπ"
@test replace("äƀçđ", r"ƀ+", "π") == "äπçđ"
@test replace("äƀçđ", r"ƀ?ç?", "π") == "πäπđπ"
@test replace("äƀçđ", r"[ƀç]?", "π") == "πäππđπ"

@test replace("foobarfoo", r"(fo|ba)", "ẍẍ") == "ẍẍoẍẍrẍẍo"

@test replace("ḟøøbarḟøø", r"(ḟø|ba)", "xx") == "xxøxxrxxø"
@test replace("ḟøøbarḟøø", r"(ḟøø|ba)", "bar") == "barbarrbar"

@test replace("fooƀäṙfoo", r"(fo|ƀä)", "xx") == "xxoxxṙxxo"
@test replace("fooƀäṙfoo", r"(foo|ƀä)", "ƀäṙ") == "ƀäṙƀäṙṙƀäṙ"

@test replace("ḟøøƀäṙḟøø", r"(ḟø|ƀä)", "xx") == "xxøxxṙxxø"
@test replace("ḟøøƀäṙḟøø", r"(ḟøø|ƀä)", "ƀäṙ") == "ƀäṙƀäṙṙƀäṙ"

# {begins,ends}with
@test beginswith("abcd", 'a')
@test beginswith("abcd", "a")
@test beginswith("abcd", "ab")
@test !beginswith("ab", "abcd")
@test !beginswith("abcd", "bc")
@test endswith("abcd", 'd')
@test endswith("abcd", "d")
@test endswith("abcd", "cd")
@test !endswith("abcd", "dc")
@test !endswith("cd", "abcd")

# RepStrings and SubStrings
u8str2 = u8str^2
len_u8str = length(u8str)
slen_u8str = length(u8str)
len_u8str2 = length(u8str2)
slen_u8str2 = length(u8str2)

@test len_u8str2 == 2 * len_u8str
@test slen_u8str2 == 2 * slen_u8str

u8str2plain = utf8(u8str2)

for i1 = 1:length(u8str2)
    if !isvalid(u8str2, i1); continue; end
    for i2 = i1:length(u8str2)
        if !isvalid(u8str2, i2); continue; end
        @test length(u8str2[i1:i2]) == length(u8str2plain[i1:i2])
        @test length(u8str2[i1:i2]) == length(u8str2plain[i1:i2])
        @test u8str2[i1:i2] == u8str2plain[i1:i2]
    end
end

str = "aa\u2200\u2222bb"
u = SubString(str, 3, 6)
@test length(u)==2
b = IOBuffer()
write(b, u)
@test takebuf_string(b) == "\u2200\u2222"

str = "føøbar"
u = SubString(str, 4, 3)
@test length(u)==0
b = IOBuffer()
write(b, u)
@test takebuf_string(b) == ""

str = "føøbar"
u = SubString(str, 10, 10)
@test length(u)==0
b = IOBuffer()
write(b, u)
@test takebuf_string(b) == ""

@test replace("\u2202", '*', '\0') == "\u2202"

# search and SubString (issue #5679)
str = "Hello, world!"
u = SubString(str, 1, 5)
@test rsearch(u, "World") == 0:-1
@test rsearch(u, 'z') == 0
@test rsearch(u, "ll") == 3:4

# quotes + interpolation (issue #455)
@test "$("string")" == "string"
arr = ["a","b","c"]
@test "[$(join(arr, " - "))]" == "[a - b - c]"

# string iteration, and issue #1454
str = "é"
str_a = [str...]
@test length(str_a)==1
@test str_a[1] == str[1]

str = "s\u2200"
@test str[1:end] == str

# triple-quote delimited strings
@test """abc""" == "abc"
@test """ab"c""" == "ab\"c"
@test """ab""c""" == "ab\"\"c"
@test """ab"\"c""" == "ab\"\"c"
@test """abc\"""" == "abc\""
n = 3
@test """$n\n""" == "$n\n"
@test """$(n)""" == "3"
@test """$(2n)""" == "6"
@test """$(n+4)""" == "7"
@test """$("string")""" == "string"
a = [3,1,2]
@test """$(a[2])""" == "1"
@test """$(a[3]+7)""" == "9"
@test """$(ifloor(4.5))""" == "4"
nl = "
"
@test """
     a
     b

     c
     """ == "a$(nl)b$(nl)$(nl)c$(nl)"
@test """
      """ == ""
@test """x
     a
    """ == "x$(nl) a$(nl)"
@test """
     $n
   """ == "  $n$(nl)"
@test """
      a
     b
       c""" == " a$(nl)b$(nl)  c"
# note tab/space mixing
@test """
	a
     b
     """ == "   a$(nl)b$(nl)"
@test """
      a
       """ == "a$(nl)"
s = "   p"
@test """
      $s""" == "$s"
@test """
       $s
      """ == " $s$(nl)"

# bytes2hex and hex2bytes
hex_str = "d7a8fbb307d7809469ca9abcb0082e4f8d5651e46d3cdb762d02d0bf37c9e592"
bin_val = hex2bytes(hex_str)

@test div(length(hex_str), 2) == length(bin_val)
@test hex_str == bytes2hex(bin_val)

bin_val = hex2bytes("07bf")
@test bin_val[1] == 7
@test bin_val[2] == 191
@test typeof(bin_val) == Array{Uint8, 1}
@test length(bin_val) == 2

# all valid hex chars
@test "0123456789abcdefabcdef" == bytes2hex(hex2bytes("0123456789abcdefABCDEF"))

# odd size
@test_throws hex2bytes("0123456789abcdefABCDEF0")

#non-hex characters
@test_throws hex2bytes("0123456789abcdefABCDEFGH")

# sizeof
@test sizeof("abc") == 3
@test sizeof("\u2222") == 3
@test sizeof(SubString("abc\u2222def",4,4)) == 3
@test sizeof(RopeString("abc","def")) == 6

# issue #3597
@test string(UTF32String(['T', 'e', 's', 't'])[1:1], "X") == "TX"

# issue #3710
@test prevind(SubString("{var}",2,4),4) == 3

# printf
# int
@test (@sprintf "%d" typemax(Int64)) == "9223372036854775807"
@test (@sprintf "%i" 42) == "42"
@test (@sprintf "%u" 42) == "42"
@test (@sprintf "Test: %i" 42) == "Test: 42"
@test (@sprintf "%#x" 42) == "0x2a"
@test (@sprintf "%#o" 42) == "052"
@test (@sprintf "%X" 42) == "2A"
@test (@sprintf "%X" 42) == "2A"
@test (@sprintf "% i" 42) == " 42"
@test (@sprintf "%+i" 42) == "+42"
@test (@sprintf "%4i" 42) == "  42"
@test (@sprintf "%-4i" 42) == "42  "
# float
@test (@sprintf "%7.2f" 1.2345) == "   1.23"
@test (@sprintf "%-7.2f" 1.2345) == "1.23   "
@test (@sprintf "%07.2f" 1.2345) == "0001.23"
@test (@sprintf "%.0f" 1.2345) == "1"
@test (@sprintf "%#.0f" 1.2345) == "1."
# Inf / NaN handling
@test (@sprintf "%f" Inf) == "Inf"
@test (@sprintf "%f" NaN) == "NaN"
# scientific notation
@test (@sprintf "%.4e" 1.2345) == "1.2345e+00"
@test (@sprintf "%.0e" 3e142) == "3e+142"
@test (@sprintf "%#.0e" 3e142) == "3.e+142"
# chars
@test (@sprintf "%c" 65) == "A"
@test (@sprintf "%c" 'A') == "A"
@test (@sprintf "%c" 248) == "ø"
@test (@sprintf "%c" 'ø') == "ø"
# strings
@test (@sprintf "%s" "test") == "test"
@test (@sprintf "%s" "tést") == "tést"
# reasonably complex
@test (@sprintf "Test: %s%c%C%c%#-.0f." "t" 65 66 67 -42) == "Test: tABC-42.."

# issue #4183
@test split(SubString(ascii("x"), 2, 0), "y") == String[""]
@test split(SubString(utf8("x"), 2, 0), "y") == String[""]

# issue #4586
@test rsplit(RevString("ailuj"),'l') == {"ju","ia"}
@test float64(RevString("64")) === 46.0

for T = (Uint8,Int8,Uint16,Int16,Uint32,Int32,Uint64,Int64,Uint128,Int128,BigInt),
    b = 2:62, _ = 1:10
    n = T != BigInt ? rand(T) : BigInt(rand(Int128))
    @test parseint(T,base(b,n),b) == n
end

# normalize_string (Unicode normalization etc.):
@test normalize_string("\u006e\u0303", :NFC) == "\u00f1"
@test "\u006e\u0303" == normalize_string("\u00f1", :NFD)
@test normalize_string("\ufb00", :NFC) != "ff"
@test normalize_string("\ufb00", :NFKC) == "ff"
@test normalize_string("\u006e\u0303\ufb00", :NFKC) == "\u00f1"*"ff"
@test normalize_string("\u00f1\ufb00", :NFKD) == "\u006e\u0303"*"ff"
@test normalize_string("\u006e\u0303", compose=true) == "\u00f1"
@test "\u006e\u0303" == normalize_string("\u00f1", decompose=true)
@test normalize_string("\u006e\u0303\u00b5",compat=true) == "\u00f1\u03bc"
@test normalize_string("Σσς",casefold=true) == "σσσ"
@test normalize_string("∕⁄", lump=true) == "//"
@test normalize_string("\ua\n\r\r\ua", newline2lf=true) == "\ua\ua\ua\ua"
@test normalize_string("\ua\n\r\r\ua", newline2ls=true) == "\u2028\u2028\u2028\u2028"
@test normalize_string("\ua\n\r\r\ua", newline2ps=true) == "\u2029\u2029\u2029\u2029"
@test normalize_string("\u00f1", stripmark=true) == "n"
@test isempty(normalize_string("\u00ad", stripignore=true))
@test normalize_string("\t\r", stripcc=true) == "  "
@test normalize_string("\t\r", stripcc=true, newline2ls=true) == " \u2028"

#Tests from Unicode SA#15, "Unicode normalization forms"
#http://www.unicode.org/reports/tr15/

#1. Canonical equivalence
let ==(a::Array{Char},b::Array{Char}) = normalize_string(string(a...), :NFC)==normalize_string(string(b...), :NFC)
    ==(a,b) = Base.(:(==))(a,b)
    @test ['C', '̧'] == ['Ç']
    @test ['q', '̇', '̣'] == ['q', '̣', '̇']
    @test ['가'] == ['ᄀ', 'ᅡ']
    @test ['Ω'] == ['Ω']
end

#2. Compatibility Equivalence
let ==(a::Array{Char},b::Array{Char}) = normalize_string(string(a...), :NFKC)==normalize_string(string(b...), :NFKC)
    ==(a,b) = Base.(:(==))(a,b)
    @test ['ℌ'] == ['ℍ'] == ['H']
    @test ['ﻨ'] == ['ﻧ'] == ['ﻦ'] == ['ﻥ']
    @test ['①'] == ['1']
    @test ['ｶ'] == ['カ']
    @test ['︷'] == ['{']
    @test ['⁹'] == ['₉']
    @test ['㌀'] == ['ア', 'パ', 'ー', 'ト']
    @test ['¼'] == ['1', '⁄', '4']
    @test ['ǆ'] == ['d', 'ž']
end

#3. Singletons
@test normalize_string("\U212b", :NFD) == "A\U030a"
@test normalize_string("\U212b", :NFC) == "\U00c5"
@test normalize_string("\U2126", :NFC) == normalize_string("\U2126", :NFD) == "\U03a9"

#4. Canonical Composites
@test normalize_string("\U00c5", :NFC) == "\U00c5"
@test normalize_string("\U00c5", :NFD) == "A\U030a"
@test normalize_string("\U00f4", :NFC) == "\U00f4"
@test normalize_string("\U00f4", :NFD) == "o\U0302"

#5. Multiple Combining Marks
@test normalize_string("\U1e69", :NFD) == "s\U0323\U0307"
@test normalize_string("\U1e69", :NFC) == "\U1e69"
@test normalize_string("\U1e0b\U0323", :NFD) == "d\U0323\U0307"
@test normalize_string("\U1e0b\U0323", :NFC) == "\U1e0d\U0307"
@test normalize_string("q\U0307\U0323", :NFC) == "q\U0323\U0307"
@test normalize_string("q\U0307\U0323", :NFD) == "q\U0323\U0307"

#6. Compatibility Composites
@test normalize_string("\Ufb01", :NFD) == normalize_string("\Ufb01", :NFC) == "\Ufb01"
@test normalize_string("\Ufb01", :NFKD) == normalize_string("\Ufb01", :NFKC) == "fi"
@test normalize_string("2\U2075", :NFD) == normalize_string("2\U2075", :NFC) == "2\U2075"
@test normalize_string("2\U2075", :NFKD) == normalize_string("2\U2075", :NFKC) == "25"
@test normalize_string("\U1e9b\U0323", :NFD) == "\U017f\U0323\U0307"
@test normalize_string("\U1e9b\U0323", :NFC) == "\U1e9b\U0323"
@test normalize_string("\U1e9b\U0323", :NFKD) == "s\U0323\U0307"
@test normalize_string("\U1e9b\U0323", :NFKC) == "\U1e69"

# issue #5870
@test !ismatch(Regex("aa"), SubString("",1,0))
@test ismatch(Regex(""), SubString("",1,0))

# issue #6027
let
    # make symbol with invalid char
    sym = symbol(char(0xdcdb))
    @test string(sym) == "\udcdb"
    @test expand(sym) === sym
    @test parse("\udcdb = 1",1,raise=false)[1] == Expr(:error, "error normalizing identifier \udcdb: Invalid UTF-8 string")
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


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

using Markdown
using InteractiveUtils

# ╔═╡ 7c62a2c4-9159-11eb-1b54-452017ea90d9
using LinearAlgebra, SparseArrays

# ╔═╡ 7ad16586-9162-11eb-1d94-5188eb8d4c5d
using DifferentialEquations, Plots

# ╔═╡ 44df6004-9158-11eb-1a63-651d7dca957b
md"# Heat equation with Robin boundary conditions

Consider a one-dimensional rod $\Omega = (0,L)$ and time $t \in (0,T).$
"

# ╔═╡ aff71d6e-9158-11eb-0559-cb7334f4a624
L = 0.1 # Length of rod

# ╔═╡ 3771dab8-9159-11eb-3a69-cd9b94a10399
Tf = 60.0 # Final simulation time

# ╔═╡ 44dbf81c-9159-11eb-2920-27ee60a7db9a
md"The rod has the physical properties
- thermal conductivity $\lambda$,
- specific heat capacity $c$,
- mass density $\rho$.

The diffusivity constant is given by $\alpha = \frac{\lambda}{c \rho}$. The heat equation is noted as

$\dot{\vartheta}(t,x) = \alpha \frac{\partial^2}{\partial x^2} \vartheta(t,x)$

for $(t,x) \in (0,T) \times \Omega$ with initial condition

$\vartheta(0,x) = \vartheta_{0}(x)$

for $x \in \overline{\Omega}$."

# ╔═╡ 86a40dc8-9158-11eb-16a3-1f776d6bde4e
begin
	λ = 45.0;    # Thermal conductivity
	ρ = 7800.0;  # Mass density
	cap = 480.0; # Specific heat capacitivity
end

# ╔═╡ 9e8ee0a2-9158-11eb-0210-4513221b7823
α = λ/(cap * ρ)   # Diffusivity

# ╔═╡ a4f1e1e2-9158-11eb-0754-a70399a341ba
md"## Boundary conditions

On both sides of the rod, heat transfer

$\phi_{t}(t,x) = - h ~ (\vartheta(t,x) - \vartheta_{amb})$

and heat radiation

$\phi_{r}(t,x) = - \epsilon \varrho ~ (\vartheta(t,x)^4 - \vartheta_{amb}^4)$

is assumed to operate. Parameter $\vartheta_{amb}$ describes the ambient temperature, $h>0$ is called the heat transfer coefficient, $\epsilon \in (0,1)$ is called emissivity and $\varrho\approx 5.67 \cdot 10^{-8}$ is known as Stefan-Boltzmann constant. For simplicity, $k=\epsilon \cdot \varrho$.

The sum of $\phi_{t}$ and $\phi_{r}$

$\phi_{out}(t,x) = \phi_{t}(t,x) + \phi_{r}(t,x) = - h ~ (\vartheta(t,x) - \vartheta_{amb}) - k ~ (\vartheta(t,x)^4 - \vartheta_{amb}^4)$

is the flux from the rod to the environment. Thus, the natural (Robin) boundary conditions are noted as

$\lambda ~ \left.\frac{\partial}{\partial x} \vartheta(\cdot, x) \right\rvert_{x = 0} \cdot \vec{n} ~=~ \lambda ~  \left. \frac{\partial}{\partial x} \vartheta(\cdot, x) \right\rvert_{x = L} \cdot \vec{n} ~=~ \phi_{out}(t,x)$

with outer [normal vector](https://en.wikipedia.org/wiki/Normal_(geometry)) $\vec{n}$ on the left or right boundary. Here, the normal vector $\vec{n} = -1$ on the left boundary at $x=0$, and  $\vec{n} = 1$ on the right boundary at $x=L$."

# ╔═╡ 7c9e9a4a-9159-11eb-3dda-59d2aa15f32a
begin
	h = 10.0; # Heat transfer coefficient
	ϵ = 0.6;  # Emissivity
	sb = 5.67*10^(-8) # Stefan-Boltzmann constant
	k = ϵ * sb; 	  # Radiation coefficient
end

# ╔═╡ c5da798e-9160-11eb-1d16-75ce792ef2fe
θamb = 298.0 # Ambient temperature in Kelvin

# ╔═╡ 7c80218c-9159-11eb-301a-95a428fb76d9
md"### Spatial approximation

The one-dimensional rod is discretized as a one-dimensional grid with N points and the finite discretization is noted by

$\Delta x = \frac{L}{N-1}.$

The resulting grid points have the position $x^{0} = 0$, $x^{n} = n~\Delta x$ and $x^{N-1} = L$. The second order derivative is approximated using the [Taylor series](https://en.wikipedia.org/wiki/Taylor_series) as

$\left. \frac{\partial^2 f(x)}{\partial x^2} \right\rvert_{x = \tilde{x}} \approx \frac{1}{\Delta x^2} \left[ f(\tilde{x} - \Delta x) - 2 f(\tilde{x}) + f(\tilde{x} + \Delta x)  \right].$

This [Finite Difference](https://en.wikipedia.org/wiki/Finite_difference_method) scheme is written in matrix-vector form as

$\frac{\partial^2 \vartheta(t,x)}{\partial x^2} \approx \frac{1}{\Delta x^2}
\begin{pmatrix}
-2 & 1 & 0 & \cdots &  & 0 \\
1 & -2 & 1 & 0  & \cdots & \vdots \\
0 & 1 & -2 & 1 & \ddots &  \\
 & & \ddots & \ddots & \ddots &  & \\
& & & 1 & -2 & 1 \\
& & &  & 1 & -2 \\
\end{pmatrix}
~
\begin{pmatrix}
\vartheta(t, x^{0}) \\
\vartheta(t, x^{1}) \\
\vdots \\
\vartheta(t, x^{N-1}) \\
\end{pmatrix}.$
"

# ╔═╡ 3487893a-915c-11eb-3a04-1fdd31009413
N = 101 # Number of grid elements

# ╔═╡ 39663578-915c-11eb-108c-2718b706729e
Δx = L/(N-1) # Finite discretization

# ╔═╡ 3468092a-915c-11eb-1719-d1b9b04b68e0
# Diffusion matrix
M = spdiagm(-1 => ones(N-1), 0 => -2*ones(N), 1 => ones(N-1));

# ╔═╡ 344c1bb6-915c-11eb-2be7-35f340f4caa8
Matrix(M)[1:5,1:5]

# ╔═╡ 342fe676-915c-11eb-3f22-1d20d5b502ed
md"### Approximated boundary conditions

The spatial approximation of the heat equation

$\frac{\partial^2 \vartheta(t,x)}{\partial x^2} \approx \frac{1}{\Delta x^2} \left( \vartheta(\cdot, x^{n-1}) - 2 ~ \vartheta(\cdot, x^{n}) + \vartheta(\cdot, x^{n+1}) \right)$

can not be evaluated directly at the grid points $x^{0}$ and $x^{N-1}$ because they depend on values from the not-existing grid points $x^{-1}$ and $x^{N}$. Therefore, the Neumann boundary condition is used to find $\vartheta(\cdot, x^{-1})$ and $\vartheta(\cdot, x^{N})$.

The Robin boundary condition states that the heat flux at both sides is described by

$\lambda \left. \frac{\partial \vartheta(\cdot,x)}{\partial x}  \right\rvert_{x = 0} \cdot \vec{n} = -1 \cdot \lambda \frac{\vartheta(\cdot, x^{1}) - \vartheta(\cdot, x^{-1})}{2 \Delta x}   = \phi_{out}(t,x^{0})$

on the left boundary and

$\lambda \left. \frac{\partial \vartheta(\cdot,x)}{\partial x}  \right\rvert_{x = L} \cdot \vec{n} = \lambda \frac{\vartheta(\cdot, x^{N}) - \vartheta(\cdot, x^{N-2})}{2 \Delta x}  = \phi_{out}(t,x^{N-1})$

on the right boundary. Reformulating both equations lead to

$\vartheta(\cdot, x^{-1}) = \vartheta(\cdot, x^{1}) + 2 \Delta x ~ \phi_{out}(t,x^{0})$

and

$\vartheta(\cdot, x^{N}) = \vartheta(\cdot, x^{N-2})  + 2 \Delta x ~ \phi_{out}(t,x^{N-1}).$

On the left boundary one yields

$\left. \frac{\partial^2 \vartheta(t,x)}{\partial x^2} \right\rvert_{x = 0} \approx \frac{1}{\Delta x^2} \left( - 2 ~ \vartheta(\cdot, x^{0}) + 2 ~ \vartheta(\cdot, x^{1}) \right) + \frac{2}{\Delta x} \phi_{out}(t,x^{0})$

and analog on the right boundary

$\left. \frac{\partial^2 \vartheta(t,x)}{\partial x^2} \right\rvert_{x = 0} \approx \frac{1}{\Delta x^2} \left(2 ~ \vartheta(\cdot, x^{N-2}) - 2 ~ \vartheta(\cdot, x^{N-1}) \right) + \frac{2}{\Delta x} \phi_{out}(t,x^{N-1}).$

The diffusion matrix and the approximated flux are noted as

$M =
\begin{pmatrix}
-2 & 2 & 0 & \cdots &  & 0 \\
1 & -2 & 1 & 0  & \cdots & \vdots \\
0 & 1 & -2 & 1 & \ddots &  \\
 & & \ddots & \ddots & \ddots &  & \\
& & & 1 & -2 & 1 \\
& & &  & 2 & -2 \\
\end{pmatrix}, \quad

\Phi_{out}(t) =
\begin{pmatrix}
\phi_{out}(t, x^{0}) \\
0 \\
\vdots \\
\\
0 \\
\phi_{out}(t, x^{N-1}) \\
\end{pmatrix}.$
"

# ╔═╡ 3414b70c-915c-11eb-2606-4d1ad425d3ab
# First row
M[1,2] = 2;

# ╔═╡ 33cd51a0-915c-11eb-3c62-efc367e08126
# Last row
M[end,end-1] = 2;

# ╔═╡ 33af0164-915c-11eb-30f0-313d926c2dae
md"### Heat equation as ODE

After the spatial approximation the heat equation has a form of

$\dot{\theta}(t) =  \frac{\alpha}{\Delta x^2} ~ M ~ \theta(t) + \frac{2 \alpha}{\Delta x} ~ \Phi_{out}(t)$

with M as diffusion matrix, temperature vector $\quad \theta(t) = \left( \vartheta(t, x^{0}), \cdots, \vartheta(t, x^{N-1}) \right)^{\top}$. This is an ordinary differential equation (ODE) that can be solved with common solvers like forward Euler method or Runge-Kutta scheme."

# ╔═╡ 328a1ed0-9160-11eb-0892-7d0a9d05a54d
# Heat Equation as ODE
function heat_eq!(dθ, θ, p, t)

	N = size(θ)
	Φout = zeros(N);

	Φout[1] = -h * (θ[1] - θamb) - k*(θ[1]^4 - θamb^4)

	Φout[end] = -h * (θ[end] - θamb) - k*(θ[end]^4 - θamb^4)

    dθ .= (1/Δx^2) * α * M * θ + (2/Δx) * α * Φout
end

# ╔═╡ 66d5aee4-9161-11eb-1ef1-2193d9c2f6a2
md"### Initial Conditions

For simplicity, the initial data of the heat equation is assumed as $\vartheta_{0}(x)= 1000$ Kelvin for all $x \in \left[0,L\right]$. This function is approximated to gain the initial conditions of the ODE

$\theta(0) = \left( 10^3, \cdots , 10^3 \right)^{\top}.$
"

# ╔═╡ 91a272be-9162-11eb-1046-8dd9e33653e5
θ₀ = 10^3 * ones(N)

# ╔═╡ 1571f3f2-9162-11eb-0b63-ff4e25a07572
md"The upper limit of the sampling time has to

$\Delta t < \frac{1}{2} \frac{\Delta x^2}{\alpha}.$

to guarantee numerical stability."

# ╔═╡ f85a2d4a-9161-11eb-0e0f-851c0750832a
Δt_ul = (0.5*Δx^2)/α # upper limit of sampling time

# ╔═╡ 5aa78f74-9162-11eb-1bd4-c96f2b90e586
Δt = 10^(-2) # Sampling time

# ╔═╡ 66ba8b22-9162-11eb-13b0-b14dfcead641
md"## Simulation"

# ╔═╡ 861018de-9162-11eb-2a01-2354acd868e2
tspan = (0.0, Tf)

# ╔═╡ e0e033e8-9162-11eb-3fe6-05565991ce77
# 1-dimensional grid
xspan = 0 : Δx : L

# ╔═╡ 01b373f0-9163-11eb-3f35-3f645825720f
prob = ODEProblem( heat_eq!, θ₀, tspan ) # ODE Problem

# ╔═╡ 0103888c-9163-11eb-16ae-371cd6d04c45
sol = solve(prob,Euler(),dt=Δt,progress=true, save_everystep=false, save_start=true) # Solving the ODE

# ╔═╡ 08dcb90e-9163-11eb-2ff0-7129add62934
plot(xspan, sol.u[2], xlabel = "Position x", ylabel="Temperature", legend=false)

# ╔═╡ Cell order:
# ╟─44df6004-9158-11eb-1a63-651d7dca957b
# ╠═aff71d6e-9158-11eb-0559-cb7334f4a624
# ╠═3771dab8-9159-11eb-3a69-cd9b94a10399
# ╟─44dbf81c-9159-11eb-2920-27ee60a7db9a
# ╠═86a40dc8-9158-11eb-16a3-1f776d6bde4e
# ╠═9e8ee0a2-9158-11eb-0210-4513221b7823
# ╟─a4f1e1e2-9158-11eb-0754-a70399a341ba
# ╠═7c9e9a4a-9159-11eb-3dda-59d2aa15f32a
# ╠═c5da798e-9160-11eb-1d16-75ce792ef2fe
# ╟─7c80218c-9159-11eb-301a-95a428fb76d9
# ╠═7c62a2c4-9159-11eb-1b54-452017ea90d9
# ╠═3487893a-915c-11eb-3a04-1fdd31009413
# ╠═39663578-915c-11eb-108c-2718b706729e
# ╠═3468092a-915c-11eb-1719-d1b9b04b68e0
# ╠═344c1bb6-915c-11eb-2be7-35f340f4caa8
# ╟─342fe676-915c-11eb-3f22-1d20d5b502ed
# ╠═3414b70c-915c-11eb-2606-4d1ad425d3ab
# ╠═33cd51a0-915c-11eb-3c62-efc367e08126
# ╟─33af0164-915c-11eb-30f0-313d926c2dae
# ╠═328a1ed0-9160-11eb-0892-7d0a9d05a54d
# ╟─66d5aee4-9161-11eb-1ef1-2193d9c2f6a2
# ╠═91a272be-9162-11eb-1046-8dd9e33653e5
# ╟─1571f3f2-9162-11eb-0b63-ff4e25a07572
# ╠═f85a2d4a-9161-11eb-0e0f-851c0750832a
# ╠═5aa78f74-9162-11eb-1bd4-c96f2b90e586
# ╠═66ba8b22-9162-11eb-13b0-b14dfcead641
# ╠═7ad16586-9162-11eb-1d94-5188eb8d4c5d
# ╠═861018de-9162-11eb-2a01-2354acd868e2
# ╠═e0e033e8-9162-11eb-3fe6-05565991ce77
# ╠═01b373f0-9163-11eb-3f35-3f645825720f
# ╠═0103888c-9163-11eb-16ae-371cd6d04c45
# ╠═08dcb90e-9163-11eb-2ff0-7129add62934


text:
module ConditionalEquations

println("\nDemonstrating conditional equations")

using Modia

# Desired:
#   using ModiaMath: plot
#
# In order that these packages need not to be defined in the user environment, they are included via Modia:
using Modia.ModiaMath: plot

@model Conditional begin
    x = Var(start=1.0)
    u = Var()
    steadyState = Boolean(false, variability=parameter)
    cond = false
    y = Var()
@equations begin
    u = sin(time)
    if ! steadyState
      der(x) + 2x = u
    else
      0 + 2x = u
    end
    if ! cond
      y = 1
    else
      y = 2
    end
  end
end

result = simulate(Conditional, 2, logTranslation=true, removeSingularities=false)
plot(result, ("x", "y"), heading="Conditional", figure=1)

@model ConditionalInstance1 begin
    @extends Conditional(steadyState=true, cond=true)
end

result = simulate(ConditionalInstance1, 2, logTranslation=true, removeSingularities=false)
plot(result, ("x", "y"), heading="ConditionalInstance", figure=2)

@model ConditionalInstance2 begin
    @extends Conditional(steadyState=false, cond=false)
end

result = simulate(ConditionalInstance2, 2, logTranslation=true, removeSingularities=false)
plot(result, ("x", "y"), heading="ConditionalInstance", figure=3)

@model Conditional2 begin
    x = Var(start=1.0)
    u = Var()
    steadyState = Boolean(false, variability=parameter)
    cond = false
    y = if cond; Var() end
@equations begin
    u = sin(time)
    if ! steadyState #&& cond
      der(x) + 2x = u
    else
      0 + 2x = u
    end
    if cond
      y = 1
    end
  end
end

result = simulate(Conditional2, 2, logTranslation=false, removeSingularities=false)
plot(result, ("x"), heading="Conditional2", figure=1)

end


text:
# This file provides optimal integer-programming-powered algorithms for
# tree layout. JuMP is not a hard dependency of GraphLayout.jl, but if
# a user wishes to use these methods then they will need JuMP and an
# integer programming solver installed, such as GLPK, Cbc, or one of the
# commercial solvers like Gurobi or CPLEX
# Ideally, only when one of these methods is used will JuMP will be
# loaded. However, one of the phases only has a JuMP-dependent version
# so we need to always load JuMP for now.
########################################################################

using JuMP

"""
    Given a layer assignment, decide a permutation for each layer
    that minimizes edge crossings using integer programming.

    Based on the IP described in
      M. Junger, E. K. Lee, P. Mutzel, and T. Odenthal.
      A polyhedral approach to the multi-layer crossing minimization problem.
      In G. Di Battista, editor, Graph Drawing: 5th International Symposium,
      GD ’97, volume 1353 of Lecture Notes in Computer Science, pages 13–24,
      Rome, Italy, September 1997. Springer-Verlag.

    Arguments:
    adj_list        Directed graph in adjacency list format
    layers          Assignment of vertices
    layer_verts     Dictionary of layer => vertices (initial perm.)

    Returns:
    new_layer_verts An improved dictionary of layer => vertices (opt. perm.)
"""
function _ordering_ip{T}(adj_list::AdjList{T}, layers, layer_verts)
    num_layers = maximum(layers)

    m = Model()

    # Define crossing binary variables
    @variable(m, c[L=1:num_layers,        # for each layer
                 i=layer_verts[L],      # for each vertex in this layer
                 j=adj_list[i],         # and vertex in the next layer
                 k=layer_verts[L],      # for each vertex in this layer
                 l=adj_list[k]], Bin)   # and vertex in the next layer

    # Objective: minimize crossings
    @objective(m, Min, sum(c))

    # Define permutation variables for each layer
    # We'll define for both (i,j) and (j,i), and ensure they consistency
    # by adding constraints. Presolve in the IP solver will simplify
    # the problem for us by removing one of the variables.
    @variable(m, x[L=1:num_layers, i=layer_verts[L], j=layer_verts[L]], Bin)
    for L in 1:num_layers
        for i in layer_verts[L]
            for j in layer_verts[L]
                j <= i && continue  # Don't double-add
                # Ensure x[i,j] and x[j,i] are consistent
                @constraint(m, x[L,i,j] == 1 - x[L,j,i])
                # And ensure that triples are consistent
                for k in layer_verts[L]
                    k <= j && continue
                    @constraint(m, 0 <= x[L,i,j] + x[L,j,k] - x[L,i,k] <= 1)
                end
            end
        end
    end

    # Link permutations to crossings
    for L in 1:num_layers-1
        # For all (i,j)
        for i in layer_verts[L]
        for j in adj_list[i]
            # For all (k,l)
            for k in layer_verts[L]
            k == i && continue  # Can't cross if starting from same vertex!
            for l in adj_list[k]
                @constraint(m, -c[L,i,j,k,l] <= x[L+1,j,l] - x[L,i,k])
                @constraint(m,  c[L,i,j,k,l] >= x[L+1,j,l] - x[L,i,k])
            end
            end
        end
        end
    end

    # Solve the IP
    solve(m)

    # Extract permutation from solution
    x_sol = getvalue(x)
    new_layer_verts = [L => Int[] for L in 1:num_layers]
    for L in 1:num_layers
        old_perm = layer_verts[L]
        # For each vertex, count the number of times it is "in front"
        # of the other vertices. The higher this number, the earlier
        # the vertex appears in the layer.
        scores = zeros(length(old_perm))
        for (p,i) in enumerate(old_perm)
            for j in old_perm
                i == j && continue
                if round(Integer, x_sol[L,i,j]) == 1
                    # i appears before j
                    scores[p] += 1
                end
            end
        end
        new_layer_verts[L] = old_perm[sortperm(scores,rev=true)]
    end

    return new_layer_verts
end

########################################################################

"""
    Given a layer assignment and permutation, decide the coordinates for
    each vertex. The objective is to encourage straight edges, especially
    for longer edges. This function uses an integer program to decide the
    coordinates (although it is solved as a linear program), as described in
      Gansner, Emden R., et al.
      A technique for drawing directed graphs.
      Software Engineering, IEEE Transactions on 19.3 (1993): 214-230.

    Arguments:
    adj_list        Directed graph in adjacency list format
    layers          Assignment of vertices
    layer_verts     Dictionary of layer => vertices (final perm.)
    orig_n          Number of original (non-dummy) vertices
    widths          Width of each vertex
    xsep            Minimum seperation between each vertex

    Returns:
    layer_coords    For each layer and vertex, the x-coord
"""
function _coord_ip{T}(adj_list::AdjList{T}, layers, layer_verts, orig_n, widths, xsep)
    num_layers = maximum(layers)

    m = Model()

    # One variable for each vertex
    @variable(m, x[L=1:num_layers, i=layer_verts[L]] >= 0)

    # Constraint: must respect permutation, and spacign constraint
    for L in 1:num_layers
        for i in 1:length(layer_verts[L])-1
            a = layer_verts[L][i]
            b = layer_verts[L][i+1]
            @constraint(m, x[L,b] - x[L,a] >=
                (widths[a] + widths[b])/2 + xsep)
        end
    end

    # Objective: minimize total misalignment
    # Use the weights from the Ganser paper:
    #   1 if both nodes "real"
    #   2 if one of the nodes is "real"
    #   8 if neither node is "real"
    # We use absolute distance in the objective so we'll need
    # auxilary variables for each pair of edges
    obj = AffExpr()
    @variable(m, absdiff[L=1:num_layers-1,
                        i=layer_verts[L], j=adj_list[i]] >= 0)
    for L in 1:num_layers-1
        for i in layer_verts[L]
            for j in adj_list[i]
                @constraint(m, absdiff[L,i,j] >= x[L,i] - x[L+1,j])
                @constraint(m, absdiff[L,i,j] >= x[L+1,j] - x[L,i])
                if i > orig_n && j > orig_n
                    # Both dummy vertices
                    obj += 8*absdiff[L,i,j]
                elseif (i <= orig_n && j >  orig_n) ||
                       (i >  orig_n && j <= orig_n)
                    # Only one dummy vertix
                    obj += 2*absdiff[L,i,j]
                else
                    # Both real
                    obj += absdiff[L,i,j]
                end
            end
        end
    end
    @objective(m, Min, obj)

    # Solve it...
    solve(m)

    # ... and mangle the solution into shape
    x_sol = getvalue(x)
    locs_x = zeros(length(layers))
    for L in 1:num_layers
        for i in layer_verts[L]
            locs_x[i] = x_sol[L,i]
        end
    end
    return locs_x
end


text:
using OpenQuantumBase, Test

coupling = ConstantCouplings(["Z"], unit=:ħ)
cfun(t₁, t₂) = 1.0
cfun(τ) = 1.0
# TODO: add test for unitary using StaticArrays
# const Sx = SMatrix{2,2}(σx)
unitary(t) = exp(-1.0im * t * σx)
tf = 5.0
u0 = PauliVec[1][1]
ρ = u0 * u0'
kernels = [(((1, 1),), coupling, OpenQuantumBase.SingleFunctionMatrix(cfun))]
redfield = OpenQuantumBase.RedfieldLiouvillian(kernels, unitary, tf, 1e-8, 1e-6)
p = ODEParams(nothing, 5.0, (tf, t) -> t / tf)

Λ = OpenQuantumBase.quadgk((x) -> unitary(x)' * σz * unitary(x), 0, 5)[1]
dρ = zero(ρ)
redfield(dρ, ρ, p, 5.0)
@test dρ ≈ -(σz * (Λ * ρ - ρ * Λ') - (Λ * ρ - ρ * Λ') * σz) atol = 1e-6 rtol =
    1e-6

A = zero(ρ ⊗ σi)
update_vectorized_cache!(A, redfield, p, 5.0)
@test A * ρ[:] ≈ -(σz * (Λ * ρ - ρ * Λ') - (Λ * ρ - ρ * Λ') * σz)[:]

# Λ = OpenQuantumBase.quadgk((x) -> unitary(x)' * σz * unitary(x), 0, 2.5)[1]
# dρ = zero(ρ)
# redfield(dρ, ρ, p, 2.5)
# @test dρ ≈ -(σz * (Λ * ρ - ρ * Λ') - (Λ * ρ - ρ * Λ') * σz) atol = 1e-6 rtol =
#     1e-6
#
# A = zero(ρ ⊗ σi)
# update_vectorized_cache!(A, redfield, UnitTime(5.0), 2.5)
# @test A * ρ[:] ≈ -(σz*(Λ*ρ-ρ*Λ')-(Λ*ρ-ρ*Λ')*σz)[:]

# test for CustomCouplings
coupling = CustomCouplings([(s) -> σz], unit=:ħ)
bath = CustomBath(correlation=(τ) -> 1.0)
interactions = InteractionSet(Interaction(coupling, bath))
redfield = OpenQuantumBase.redfield_from_interactions(interactions, unitary, tf, 1e-8, 1e-6)[1]

A = zero(ρ ⊗ σi)
Λ = OpenQuantumBase.quadgk((x) -> unitary(x)' * σz * unitary(x), 0, 2.5)[1]
update_vectorized_cache!(A, redfield, p, 2.5)
@test A * ρ[:] ≈ -(σz * (Λ * ρ - ρ * Λ') - (Λ * ρ - ρ * Λ') * σz)[:]

# =============== CGME Test ===================
kernels = [(((1, 1),), coupling, OpenQuantumBase.SingleFunctionMatrix(cfun), 1)]
cgop = OpenQuantumBase.CGLiouvillian(kernels, unitary, 1e-8, 1e-6)
dρ = zero(ρ)
function integrand(x)
    a1 = unitary(x[1])' * σz * unitary(x[1])
    a2 = unitary(x[2])' * σz * unitary(x[2])
    a1 * ρ * a2 - 0.5 * (a2 * a1 * ρ + ρ * a2 * a1)
end
exp_res, err = OpenQuantumBase.hcubature(integrand, [-0.5, -0.5], [0.5, 0.5])
cgop(dρ, ρ, p, 2.5)
@test dρ ≈ exp_res


text:
import JSON
using VegaLite

function vega_header(title::String="Plot", other::String="")
    return """
      <head>
        <title>$title</title>
        <meta charset="UTF-8">
        <script>$(read(VegaLite.asset("vega.min.js"), String))</script>
        <script>$(read(VegaLite.asset("vega-lite.min.js"), String))</script>
        <script>$(read(VegaLite.asset("vega-embed.min.js"), String))</script>
        $other
      </head>
    """
end

function vega_style()
    return """
    <style media="screen">
      .vega-actions a {
        margin-right: 10px;
        font-family: sans-serif;
        font-size: x-small;
        font-style: italic;
      }
    </style>
    """
end

function vega_plot_html(specs::Dict{String, VegaLite.VLSpec})
    res = """
      <script type="text/javascript">
        var opt = {
          mode: "vega-lite",
          renderer: "$(VegaLite.Vega.RENDERER)",
          actions: $(VegaLite.Vega.ACTIONSLINKS)
        }
    """
    for (i, (divid,spec)) in enumerate(specs)
        res *= "var spec$i = $(JSON.json(VegaLite.add_encoding_types(VegaLite.Vega.getparams(spec))))\nvegaEmbed('#$divid', spec$i, opt);\n\n"
    end

    return res * "</script>"
end

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
# libpurex:
# Tracking

# Tracking helper types

struct CTracking
    sumw;
    CTracking(N) = new(
        Float64.(N)
    ); # makes a copy of the starting situation
end

abbrev(_::Type{CTracking}) = "C";

# Add a weight vector and track it
function track(t::CTracking, N, w)
    @assert all(N .≤ t.sumw.+1) "N $N  sumw $(t.sumw)";
    @assert sum(N) ≈ sum(t.sumw);

    t.sumw .+= w;
    argmin(N .- t.sumw);
end

struct DTracking
    DTracking(N) = new();
end

abbrev(_::Type{DTracking}) = "D";

function track(t::DTracking, N, w)
    argmin(N .- sum(N).*w);
end

# Wrapper to add forced exploration to a tracking rule
struct ForcedExploration
    t;
end

function track(fe::ForcedExploration, N, w)
    t = sum(N);
    K = length(N);
    undersampled = N .≤ sqrt(t) .- K/2;
    if any(undersampled)
        track(fe.t, N, undersampled/sum(undersampled));
    else
        track(fe.t, N, w);
    end
end


text:
import LinearAlgebra.norm

# deprecate these soon
norm_inf(x::AbstractExpr) = maximum(abs(x))
norm_1(x::AbstractExpr) = sum(abs(x))
norm_fro(x::AbstractExpr) = norm2(vec(x))

"""
    norm(x::AbstractExpr, p::Real=2)

Computes the `p`-norm `‖x‖ₚ = (∑ᵢ |xᵢ|^p)^(1/p)` of a vector expression `x`.

This function uses specialized methods for `p=1, 2, Inf`. For `p > 1` otherwise,
this function uses the procedure documented at
[`rational_to_socp.pdf`](https://github.com/jump-dev/Convex.jl/raw/master/docs/supplementary/rational_to_socp.pdf),
based on the paper "Second-order cone programming" by F. Alizadeh and D. Goldfarb,
Mathematical Programming, Series B, 95:3-51, 2001.

!!! warning
    For versions of Convex.jl prior to v0.14.0, `norm` on a matrix expression returned
    the operator norm ([`opnorm`](@ref)), which matches Julia v0.6 behavior. This functionality
    was deprecated since Convex.jl v0.8.0, and has been removed. In the future,
    `norm(x, p)` will return `‖vec(x)‖ₚ`, matching the behavior of [`norm`](@ref)
    for numeric matrices.
"""
function LinearAlgebra.norm(x::AbstractExpr, p::Real=2)
    if length(size(x)) <= 1 || minimum(size(x))==1
        # x is a vector
        if p == 1
            return norm_1(x)
        elseif p == 2
            return norm2(x)
        elseif p == Inf
            return norm_inf(x)
        elseif p > 1
            # TODO: allow tolerance in the rationalize step
            return rationalnorm(x, rationalize(Int, float(p)))
        else
            error("vector p-norms not defined for p < 1")
        end
    else
        error("In Convex.jl v0.13 and below, `norm(x, p)` meant `opnorm(x, p)` (but was deprecated since v0.8.0). In the future, `norm(x,p)` for matrices will be equivalent to `norm(vec(x),p)`. This is currently an error to ensure you update your code!")
    end
end