text:
push!(LOAD_PATH, "Modules")
using LinearAlgebra
using Colors
using Revise
using JuliRay
using Statistics
RealVector=JuliRay.RealVector
RealMatrix=JuliRay.RealMatrix
Base.@irrational ° 0.0174532925199432957692369076848861271344 (big(pi)/big(180))
φ=MathConstants.φ
include("Polyhedra.jl")

##
ε=10^(-1.2)

function f(x)
    ifelse(x>0,exp(-1/x),0)
end
function g(x)
    f(2x/√3)/(f(2x/√3)+f(2/√3-2x/√3))
end
function Smooth(a,b,x)
    g((x-a)/(b-a))
end

function ℝ³⭢S³(p::Array{T,1}) where T<:Real
    if length(p)≠3
        error("No Point on ℝ³")
    end
    return [2p[1]/(1+p[1]^2+p[2]^2+p[3]^2),2p[2]/(1+p[1]^2+p[2]^2+p[3]^2),2p[3]/(1+p[1]^2+p[2]^2+p[3]^2),(-1+p[1]^2+p[2]^2+p[3]^2)/(1+p[1]^2+p[2]^2+p[3]^2)]
end
function ℝ³⭢S³(p::Array{T,1},θ::Real) where T<:Real
    if length(p)≠3
        error("No Point on ℝ³")
    elseif θ<ε
        return [p[1:3]...,0.0]
    else
        c=[0,0,0,cot(θ)]
        r=1/sin(θ)
        P=p*cot((π/2-θ)/2+π/4)
        return ℝ³⭢S³(P)*r+c
    end
end
function S³⭢ℝ³(q::Array{T,1}) where T<:Real
    if length(q)≠4
        error("No Point on S³")
    elseif !(norm(q) ≈ 1.0)
        print(norm(q))
        error("No Point on unit S³, radius: $(norm(q))")
    end
    return [q[1]/(1-q[4]),q[2]/(1-q[4]),q[3]/(1-q[4])]
end
function S³⭢ℝ³(q::Array{T,1},θ::Real) where T<:Real
    if length(q)≠4
        error("No Point on S³")
    elseif θ<ε
        return q[1:3]
    else
        c=[0,0,0,cot(θ)]
        r=1/sin(θ)
        return S³⭢ℝ³((q-c)/r)*tan((π/2-θ)/2+π/4)
    end
end
function ℝ⁴⭢S³(p::Array{T,1},θ::Real) where T<:Real
    if θ<ε
        return [p[1:3]...,0.0]
    else
        R=1/sin(θ)
        O=[0,0,0,cot(θ)]
        return O+R*normalize(p-O)
    end
end

function NormalVector(p₁::RealVector,p₂::RealVector,p₃::RealVector,p₄::RealVector)
    A=hcat(p₁-p₄,p₂-p₄,p₃-p₄)
    return normalize([(-1)^i*det(A[deleteat!(collect(1:4),i),:]) for i ∈ 1:4])
end

function Mirror(q::RealVector,p₁::RealVector,p₂::RealVector,p₃::RealVector)
    𝒏=JuliRay.NormalVector(p₁,p₂,p₃)
    return q-2*dot(𝒏,q-p₁)*𝒏
end
function Mirror(q::RealVector,p₁::RealVector,p₂::RealVector,p₃::RealVector,p₄::RealVector)
    𝒏=JuliRay.NormalVector(p₁,p₂,p₃,p₄)
    return q-2*dot(𝒏,q-p₁)*𝒏
end
function Mirror(q::RealVector,p₁::RealVector,p₂::RealVector,p₃::RealVector,p₄::RealVector, θ)
    𝒏=NormalVector(p₁,p₂,p₃,p₄)
    return ℝ⁴⭢S³(q-2*dot(𝒏,q-p₁)*𝒏,θ)
end

function PickFace(cell::CELL, v::RealVector,POINTS³)
    if norm(v) ≈ 0
        error("vector v must be non-zero")
    end
    faces=vertices.(cell)
    fpts=(i->POINTS³[i]).(faces)
    return cell[findmax([dot(v,+(pts...)) for pts ∈ fpts])[2]]
end

function NewCell(cell::CELL,v::RealVector,θ,POINTS³, POINTS⁴)
    n=length(POINTS³)
    c=copy(cell)
    face=PickFace(cell,v,POINTS³)
    IND_face=vertices(face)
    IND_cell=Int[]
    for f ∈ c for e ∈ f for v ∈ e push!(IND_cell,v) end end end
    IND_cell=union(IND_cell)
    IND_cell2=[
        if i ∈ IND_face
            i
        else
            push!(POINTS³,Mirror(POINTS³[i],POINTS³[IND_face[1]],POINTS³[IND_face[2]],POINTS³[IND_face[3]]));
            if θ<ε
                push!(POINTS⁴,[POINTS³[end]...,0]);
            else
                O=[0,0,0,cot(θ)]
                push!(POINTS⁴,Mirror(POINTS⁴[i],O,POINTS⁴[IND_face[1]],POINTS⁴[IND_face[2]],POINTS⁴[IND_face[3]],θ));
            end
            n=n+1
        end
        for i ∈ IND_cell]
    return [[[IND_cell2[findfirst(w->w==v,IND_cell)] for v ∈ e] for e ∈ f] for f ∈ c]
end

function SphericalSphere(v,r::Real,θ::Real) where T<:RealVector
    V=S³⭢ℝ³(v,θ)
    return Sphere(V,r)
end

function SphericalCylinder(v₁,v₂,r::Real,θ::Real) where T<:RealVector
    w₁=ℝ⁴⭢S³((v₁+v₂)/2,θ)
    V₁=S³⭢ℝ³(v₁,θ)
    V₂=S³⭢ℝ³(v₂,θ)
    W₁=S³⭢ℝ³(w₁,θ)
    return Arc(V₁,W₁,V₂,r)
end

function strictplane(ps::Array{Array{Float64,1}})
    m = length(ps) # m-gon
    c = mean(ps) # center
    δs = [ps[mod(i+2,m)+1] - ps[i] for i in 1:m] # edge
    qs = [ps[i]-c for i in 1:m]
    ns = [cross(qs[i],δs[i]) for i in 1:m] # normal
    n = normalize(mean(ns))
    ps′ = [qs[i]-dot(qs[i],n)*n+c for i in 1:m]
    return ps′
end

function SphericalPolygon(v::Array{T,1},θ::Real) where T<:RealVector
    n=length(v)
    u=ℝ⁴⭢S³(+(v...)/n,θ)
    v₁=v[1]
    v₂=v[2]
    v₃=v[3]
    V=(q->S³⭢ℝ³(q,θ)).(v)
    U=S³⭢ℝ³(u,θ)
    w=[ℝ⁴⭢S³((v[i]+v[mod(i,length(v))+1])/2,θ) for i ∈ 1:length(v)]
    W=(q->S³⭢ℝ³(q,θ)).(w)
    if θ<ε
    # if θ<10^(-3)
        return Polygon(V)
    elseif rank(hcat(V...),atol=1.0e-4)==2
    # elseif rank(hcat(V...),atol=1.0e-12)==2
        m=4
        vw=copy(v)
        # m-gon
        for _ ∈ 1:m
            l=length(vw)
            vw=[ℝ⁴⭢S³((vw[(i+1)÷2]+vw[mod(i÷2,l)+1])/2,θ) for i ∈ 1:2l]
        end
        VW=(q->S³⭢ℝ³(q,θ)).(vw)
        # return rgbftColor(csgUnion(Polygon(strictplane(VW))),RGB(1,0,0),FT(0.3,0.3))
        return csgUnion(Polygon(strictplane(VW)))
    else
        O=JuliRay.Circumcenter(U,V[1],V[2],V[3])
        sphere=Sphere(O,norm(U-O))
        N=JuliRay.NormalVector(V[1],W[1],V[2]);
        direction=sign(dot(N,O-U))

        cylinders=csgIntersection(
                (V₁=V[i];
                V₂=V[mod(i,n)+1];
                V₃=V[mod(i+1,n)+1];
                C=JuliRay.Circumcenter(V₁,W[i],V₂);
                N=JuliRay.NormalVector(V₁,W[i],V₂);
                direction=sign(dot(N,U-W[i]));
                cylinder=Cylinder(C,C+2*direction*norm(U-O)*N,norm(U-O)))
                for i ∈ 1:n)
        return csgClip(sphere,cylinders)
    end
end

function approxin(α,A)
    for a in A
        if norm(a - α) < 0.05
            return true
        end
    end
    return false
end

function Cells2Object(cells::Array{CELL,1},θ,POINTS⁴;rᵥ=0.05,rₑ=0.025,color=RGB(1,1,1))
    cs=copy(cells)
    # POINTS⁴=copy(POINTS⁴)
    fs=vcat(cs...)
    nf = length(fs)

    Cs = [S³⭢ℝ³(ℝ⁴⭢S³(mean([POINTS⁴[i] for i ∈ vertices(f)]),θ),θ) for f in fs]

    # get unique center points
    index_double_unique = []
    Cs_double = []
    Cs_dummy = []
    for i in 1:nf
        if approxin(Cs[i], Cs_dummy)
            push!(Cs_double,Cs[i])
            push!(index_double_unique,i)
        end
        push!(Cs_dummy,Cs[i])
    end

    # get double indices
    index_double = []
    for i in 1:nf
        if approxin(Cs[i], Cs_double)
            push!(index_double,i)
        end
    end

    # get single indices
    index_single = setdiff(1:nf, index_double)

    fs_double = fs[index_double_unique]
    fs_single = fs[index_single]

    # fs_b
    es=unique(vcat(fs...))
    vs=unique(vcat(es...))
    V=rgbColor(csgUnion([SphericalSphere(POINTS⁴[v],rᵥ,θ) for v ∈ vs]),RGB(0.1,0.1,0.1))
    E=rgbColor(csgUnion([SphericalCylinder(POINTS⁴[e[1]],POINTS⁴[e[2]],rₑ,θ) for e ∈ es]),RGB(0.2,0.2,0.2))
    # F=rgbftColor(csgUnion([SphericalPolygon([POINTS⁴[i] for i ∈ vertices(f)],θ) for f ∈ fs]),color,FT(0,0))
    # F=rgbftColor(csgUnion([SphericalPolygon([POINTS⁴[i] for i ∈ vertices(f)],θ) for f ∈ fs]),color,FT(0.4,0.2))
    F_double=rgbftColor2(csgUnion([SphericalPolygon([POINTS⁴[i] for i ∈ vertices(f)],θ) for f ∈ fs_double]),color,FT(0.6,0.3))
    F_single=rgbftColor(csgUnion([SphericalPolygon([POINTS⁴[i] for i ∈ vertices(f)],θ) for f ∈ fs_single]),color,FT(0.6,0.3))
    return csgUnion(V,E,F_double, F_single)
end

##%

function C₈(θ)
    h=1/2
    r=√(1-h^2)

    R=1/sin(θ)
    O=[0,0,0,cot(θ)]
    N=O+[0,0,0,R]
    H=O-[0,0,0,√(R^2-r^2)]

    POINTS, c₁=F₆
    POINTS=r*copy(POINTS)
    POINTS³=copy(POINTS)
    POINTS⁴=[[𝒑...,0]+H for 𝒑 ∈ POINTS]
    c₂=NewCell(c₁,[-1,0,0],θ,POINTS³, POINTS⁴)
    c₃=NewCell(c₂,[-1,0,0],θ,POINTS³, POINTS⁴)
    c₄=NewCell(c₁,[1,0,0],θ,POINTS³, POINTS⁴)
    c₅=NewCell(c₁,[0,0,-1],θ,POINTS³, POINTS⁴)
    c₆=NewCell(c₅,[0,-1,0],θ,POINTS³, POINTS⁴)
    c₇=NewCell(c₆,[0,-1,0],θ,POINTS³, POINTS⁴)
    c₈=NewCell(c₅,[0,1,0],θ,POINTS³, POINTS⁴)

    cells=[c₁,c₂,c₃,c₄,c₅,c₆,c₇,c₈]
    return cells, POINTS³, POINTS⁴
end

##%
M=120
for i ∈ 0:2M-1
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₈(θ)
    render(Cells2Object(cells,θ,POINTS⁴,color=RGB(1.0,0.5,0.1),rₑ=0.012,rᵥ=0.03),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.15,width=1200,height=900),name="8-Cell",index=i+1)
end

##%

function C₁₆(θ)
    h=1/2
    r=√(1-h^2)

    R=1/sin(θ)
    O=[0,0,0,cot(θ)]
    N=O+[0,0,0,R]
    H=O-[0,0,0,√(R^2-r^2)]

    POINTS, c₁=F₄
    POINTS=r*copy(POINTS)
    POINTS³=copy(POINTS)
    POINTS⁴=[[𝒑...,0]+H for 𝒑 ∈ POINTS]
    V=[1,-2,0]
    c₂=NewCell(c₁,V,θ,POINTS³, POINTS⁴)
    c₃=NewCell(c₂,V,θ,POINTS³, POINTS⁴)
    c₄=NewCell(c₃,V,θ,POINTS³, POINTS⁴)
    c₅=NewCell(c₁,-V,θ,POINTS³, POINTS⁴)
    c₆=NewCell(c₅,-V,θ,POINTS³, POINTS⁴)
    c₇=NewCell(c₆,-V,θ,POINTS³, POINTS⁴)
    c₈=NewCell(c₇,-V,θ,POINTS³, POINTS⁴)
    c₉=NewCell(c₁,[-2,-1,0],θ,POINTS³, POINTS⁴)
    W=[-1,-0.5,1]
    c₁₀=NewCell(c₉,W,θ,POINTS³, POINTS⁴)
    c₁₁=NewCell(c₁₀,W,θ,POINTS³, POINTS⁴)
    c₁₂=NewCell(c₁₁,W,θ,POINTS³, POINTS⁴)
    c₁₃=NewCell(c₉,-W,θ,POINTS³, POINTS⁴)
    c₁₄=NewCell(c₁₃,-W,θ,POINTS³, POINTS⁴)
    c₁₅=NewCell(c₁₄,-W,θ,POINTS³, POINTS⁴)
    c₁₆=NewCell(c₁₅,-W,θ,POINTS³, POINTS⁴)

    cells=[c₁,c₂,c₃,c₄,c₅,c₆,c₇,c₈,c₉,c₁₀,c₁₁,c₁₂,c₁₃,c₁₄,c₁₅,c₁₆]
    return cells, POINTS³, POINTS⁴
end

##%

function C₂₄(θ)
    h=1/√2
    r=√(1-h^2)

    R=1/sin(θ)
    O=[0,0,0,cot(θ)]
    N=O+[0,0,0,R]
    H=O-[0,0,0,√(R^2-r^2)]

    POINTS, c₁=F₈
    POINTS=r*copy(POINTS)
    POINTS³=copy(POINTS)
    POINTS⁴=[[𝒑...,0]+H for 𝒑 ∈ POINTS]
    c₂=NewCell(c₁,[1,1,1],θ,POINTS³, POINTS⁴)
    c₃=NewCell(c₂,[1,1,1],θ,POINTS³, POINTS⁴)
    c₄=NewCell(c₃,[1,1,1],θ,POINTS³, POINTS⁴)
    c₅=NewCell(c₁,-[1,1,1],θ,POINTS³, POINTS⁴)
    c₆=NewCell(c₅,-[1,1,1],θ,POINTS³, POINTS⁴)
    c₇=NewCell(c₁,[1,-1,-1],θ,POINTS³, POINTS⁴)
    c₈=NewCell(c₁,[-1,1,-1],θ,POINTS³, POINTS⁴)
    c₉=NewCell(c₁,[-1,-1,1],θ,POINTS³, POINTS⁴)
    c₁₀=NewCell(c₂,[-1,0,0],θ,POINTS³, POINTS⁴)
    c₁₁=NewCell(c₂,[0,-1,0],θ,POINTS³, POINTS⁴)
    c₁₂=NewCell(c₂,[0,0,-1],θ,POINTS³, POINTS⁴)
    c₁₃=NewCell(c₃,[1,-1,-1],θ,POINTS³, POINTS⁴)
    c₁₄=NewCell(c₃,[-1,1,-1],θ,POINTS³, POINTS⁴)
    c₁₅=NewCell(c₃,[-1,-1,1],θ,POINTS³, POINTS⁴)
    c₁₆=NewCell(c₄,[-1,0,0],θ,POINTS³, POINTS⁴)
    c₁₇=NewCell(c₄,[0,-1,0],θ,POINTS³, POINTS⁴)
    c₁₈=NewCell(c₄,[0,0,-1],θ,POINTS³, POINTS⁴)
    c₁₉=NewCell(c₆,[1,-1,-1],θ,POINTS³, POINTS⁴)
    c₂₀=NewCell(c₆,[-1,1,-1],θ,POINTS³, POINTS⁴)
    c₂₁=NewCell(c₆,[-1,-1,1],θ,POINTS³, POINTS⁴)
    c₂₂=NewCell(c₅,[-1,0,0],θ,POINTS³, POINTS⁴)
    c₂₃=NewCell(c₅,[0,-1,0],θ,POINTS³, POINTS⁴)
    c₂₄=NewCell(c₅,[0,0,-1],θ,POINTS³, POINTS⁴)

    cells=[c₁,c₂,c₃,c₄,c₅,c₆,c₇,c₈,c₉,c₁₀,c₁₁,c₁₂,c₁₃,c₁₄,c₁₅,c₁₆,c₁₇,c₁₈,c₁₉,c₂₀,c₂₁,c₂₂,c₂₃,c₂₄]
    return cells, POINTS³, POINTS⁴
end

##%

function C₁₂₀(θ)
    h=norm([φ^3,0,0,0])/norm([φ^3,1,1,1])
    r=√(1-h^2)

    R=1/sin(θ)
    O=[0,0,0,cot(θ)]
    N=O+[0,0,0,R]
    H=O-[0,0,0,√(R^2-r^2)]

    POINTS, c₁=F₁₂
    POINTS=r*copy(POINTS)
    POINTS³=copy(POINTS)
    POINTS⁴=[[𝒑...,0]+H for 𝒑 ∈ POINTS]
    Va,_=F₂₀
    Vb=[v-2*Va[1]*dot(Va[1],v) for v ∈ Va]
    c₂=NewCell(c₁,Va[1],θ,POINTS³, POINTS⁴)
    c₃=NewCell(c₂,Va[1],θ,POINTS³, POINTS⁴)
    c₄=NewCell(c₃,Va[1],θ,POINTS³, POINTS⁴)
    c₅=NewCell(c₄,Va[1],θ,POINTS³, POINTS⁴)
    c₆=NewCell(c₅,Va[1],θ,POINTS³, POINTS⁴)
    c₇=NewCell(c₁,-Va[1],θ,POINTS³, POINTS⁴)
    c₈=NewCell(c₇,-Va[1],θ,POINTS³, POINTS⁴)
    c₉=NewCell(c₈,-Va[1],θ,POINTS³, POINTS⁴)
    c₁₀=NewCell(c₉,-Va[1],θ,POINTS³, POINTS⁴)

    cells=[c₁,c₂,c₃,c₄,c₅,c₆,c₇,c₈,c₉,c₁₀]

    push!(cells,NewCell(cells[1],Va[3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[1],Va[7],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[1],Va[8],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[1],Va[10],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[1],Va[12],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[2],Vb[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[2],Vb[5],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[2],Vb[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[2],Vb[9],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[2],Vb[11],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[3],Va[3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[3],Va[7],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[3],Va[8],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[3],Va[10],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[3],Va[12],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[4],Vb[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[4],Vb[5],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[4],Vb[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[4],Vb[9],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[4],Vb[11],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[5],Va[3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[5],Va[7],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[5],Va[8],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[5],Va[10],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[5],Va[12],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[6],Vb[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[6],Vb[5],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[6],Vb[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[6],Vb[9],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[6],Vb[11],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[7],Vb[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[7],Vb[5],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[7],Vb[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[7],Vb[9],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[7],Vb[11],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[8],Va[3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[8],Va[7],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[8],Va[8],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[8],Va[10],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[8],Va[12],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[9],Vb[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[9],Vb[5],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[9],Vb[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[9],Vb[9],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[9],Vb[11],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[10],Va[3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[10],Va[7],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[10],Va[8],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[10],Va[10],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[10],Va[12],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[18],Va[12],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end],Va[5],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end],-Va[11],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end],-Va[11],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end],-Va[11],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end],-Va[11],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end-4],Va[11],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end],Va[11],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end],Va[11],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end],Va[11],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[end],Va[11],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[17+45],Va[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[17+45],Va[4],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[17+45],Va[12],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[17+45],Vb[12],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[18+45],Va[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[18+45],Va[12],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[18+45],Va[1],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[18+45],Va[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[18+45],Va[8],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[19+45],[-1,-2,3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[19+45],[1,-1,3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[19+45],Vb[12],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[19+45],Vb[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[19+45],[-1,-1,-1],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[20+45],Va[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[20+45],Va[12],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[20+45],Va[1],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[20+45],Va[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[20+45],Va[8],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[21+45],[-1,-2,3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[21+45],[1,-1,3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[21+45],Vb[12],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[21+45],Vb[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[21+45],[-1,-1,-1],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[22+45],Va[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[22+45],Va[12],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[22+45],Va[1],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[22+45],Va[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[22+45],Va[8],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[23+45],[-1,-2,3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[23+45],[1,-1,3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[23+45],[2,-2,-1],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[23+45],Vb[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[23+45],[-1,-1,-1],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[24+45],Va[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[24+45],Va[12],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[24+45],Va[1],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[24+45],Va[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[24+45],Va[8],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[25+45],[-1,-2,3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[25+45],[1,-1,3],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[25+45],[2,-2,-1],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[25+45],Vb[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[25+45],[-1,-1,-1],θ,POINTS³, POINTS⁴))

    push!(cells,NewCell(cells[26+45],Va[2],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[26+45],Va[12],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[26+45],Va[1],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[26+45],Va[6],θ,POINTS³, POINTS⁴))
    push!(cells,NewCell(cells[26+45],Va[8],θ,POINTS³, POINTS⁴))

    return cells, POINTS³, POINTS⁴
end

##
M=240
for i ∈ 0:2M-1
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₁₆(θ)
    render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.1,1.0,0.1),rₑ=0.012,rᵥ=0.03),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.15,width=1200,height=900),name="16-Cell",index=i+1)
end

##

M=240
for i ∈ 0:2M-1
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₈(θ)
    render(Cells2Object(cells,θ,POINTS⁴,color=RGB(1.0,0.5,0.1),rₑ=0.012,rᵥ=0.03),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.15,width=1200,height=900),name="8-Cell",index=i+1)
end

##

M=240
for i ∈ 0:2M-1
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₂₄(θ)
    render(Cells2Object(cells,θ,POINTS⁴,color=RGB(1.0,0.1,1.0),rₑ=0.012,rᵥ=0.03),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.12,width=1200,height=900),name="24-Cell",index=i+1)
end

##

M=240
for i ∈ 0:2M-1
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₁₂₀(θ)
    try
        render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.1,0.1,1.0),rₑ=0.012,rᵥ=0.03),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.12,width=1200,height=900),name="120-Cell",index=i+1)
    catch
    end
end

dropped = Int[]
M=150
for i ∈ 0:2M-1
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₁₂₀(θ)
    try
        render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.1,0.3,1.0),rₑ=0.012,rᵥ=0.03),camera=LngLatCamera(lng=180°+180°*i/M,lat=25°,pers=0.2,zoom=0.075,width=500,height=500),name="120e-Cell",index=i+1)
    catch
        push!(dropped,i)
    end
end

dropped

for i ∈ 0:39
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₁₂₀(θ)
    try
        # render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.1,0.3,1.0),rₑ=0.012,rᵥ=0.03),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.075,width=500,height=500),name="120a-Cell",index=i+1)
        render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.1,0.3,1.0),rₑ=0.012,rᵥ=0.03),camera=LngLatCamera(lng=180°,lat=25°,pers=0.2,zoom=0.075,width=500,height=500),name="120b-Cell",index=i+1)
    catch
    end
end

try
    [1][3]
catch
end

## %
M=120
for i ∈ 0:2M-1
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₈(θ)
    render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.2,1,1),rₑ=0.012,rᵥ=0.03),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.15,width=1200,height=900),name="8-Cell",index=i+1)
end

M=120
for i ∈ 80
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₁₆(θ)
    render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.2,1,1)),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.12,width=600,height=450),name="16-Cell",index=i+1)
end

M=120
for i ∈ 50
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₂₄(θ)
    render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.2,1,1)),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.12,width=600,height=450),name="24-Cell",index=i+1)
end

M=120
i=24
θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
cells, POINTS³, POINTS⁴=C₁₂₀(θ)
render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.2,1,1)),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.12,width=600,height=450),name="120-CellB",index=i+1)

M=120
for i ∈ 0:2M
    θ=π/2*(Smooth(0,1,i/M)-Smooth(1,2,i/M))
    cells, POINTS³, POINTS⁴=C₁₂₀(θ)
    try
        render(Cells2Object(cells,θ,POINTS⁴,color=RGB(0.2,1,1)),camera=LngLatCamera(lng=180°+360°*i/M,lat=25°,pers=0.2,zoom=0.07,width=500,height=500),name="120-CellD",index=i+1)
    catch
        println(i)
    end
end

3

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
# a Hamiltonian
struct RGHamiltonian
    # terms that anticommute with some other terms
    A::Vector{Term}
    # terms that commute with all other terms
    C::Vector{Term}
end

# constructor
RGHamiltonian(A::Vector{Term}) = RGHamiltonian(A, typeof(A)[])

# retrieve largest anticommuting term
function popH0!(H::RGHamiltonian)
    # init values to be overridden
    iH0 = -1
    h0 = 0.
    # run over all terms and compare
    for (i, T) in enumerate(H.A)
        h = norm(T)
        if h > h0
            # this term is the largest so far
            iH0 = i
            h0 = h
        end
    end
    # return H0 and remove it from H
    H0 = H.A[iH0]
    deleteat!(H.A, iH0)
    H0
end

# retrieve terms that anticommute with H0
function popS!(H::RGHamiltonian, H0::Term)
    # terms that anticommute with H0
    iS = Int[]
    for (i, T) in enumerate(H.A)
        if anticommute(T, H0)
            push!(iS, i)
        end
    end
    S = H.A[iS]
    deleteat!(H.A, iS)
    S
end

# compute the renormalized H0 ("on diagonal" Schrieffer-Wolff term)
function renormalizeH0(H0::Term, S::Vector{Term})
    # sum of h^2 over terms in S
    sumh2 = 0.
    for T in S
        sumh2 += T.h^2
    end
    # renormalized h0
    h0r = H0.h + sumh2 / (2 * H0.h)
    # renormalized H0
    h0r * H0.O
end

# compute the "off diagonal" Schrieffer-Wolff terms
function schriefferwolff(H0::Term, S::Vector{Term})
    # "off diagonal" terms in the SW correction
    SW = Term[]
    lS = length(S)
    for i in 1:lS
        for j in i+1:lS
            if commute(S[i], S[j])
                push!(SW, ((1 / H0.h) * H0.O) * S[i] * S[j])
            end
        end
    end
    # combine terms by adding them
    combine!(SW)
    # drop any zero terms
    dropzeros!(SW)
    SW
end

# perform an RG step
function step!(H::RGHamiltonian, maxrate::Int = 2)
    # retrieve the largest term
    H0 = popH0!(H)
    # retrieve anticommuting terms
    S = popS!(H, H0)
    # the renormalized H0 term
    rH0 = renormalizeH0(H0, S)
    # add the renormalized H0 to H
    push!(H.C, rH0)
    # the SW terms
    SW = schriefferwolff(H0, S)
    # truncate SW terms
    maxlen = maxrate * length(S)
    if maxlen < length(SW)
        sort!(SW, by=norm, rev=true)
        resize!(SW, maxlen)
    end
    # add SW terms to H
    append!(H.A, SW)
    # add terms together if possible and drop zeros
    combine!(H.A)
    dropzeros!(H.A)
    # return the G value of this step
    isempty(S) ? nothing : log(maximum(norm, S)) - log(norm(H0))
end

ext: jl

max_stars_repo_name: aormorningstar/SBRG

max_issues_repo_name: aormorningstar/SBRG

max_forks_repo_name: aormorningstar/SBRG


text:
const LEAP_YEAR = 4
const CENTURIAL_YEAR = 100
const LEAP_CYCLE_LENGTH = LEAP_YEAR * CENTURIAL_YEAR

function is_leap_year(year::Int)
    is_leap(year) && (is_non_centurial(year) || is_leap_cycle(year))
end

function is_leap(year::Int)
    year % LEAP_YEAR == 0
end

function is_non_centurial(year::Int)
    year % CENTURIAL_YEAR != 0
end

function is_leap_cycle(year::Int)
    year % LEAP_CYCLE_LENGTH == 0
end

max_stars_repo_name: paulfioravanti/exercism

max_issues_repo_name: paulfioravanti/exercism

max_forks_repo_name: paulfioravanti/exercism


text:
filename = "ff_collapse.d"
infile = open(filename)
ρ      = 1.0E4
av     = 10.0E0
temp   = 10.0E0

outfile = open("isochoric_$filename","w")
for ln in eachline(infile)
#  println(ln)
  newln = ln
  ρ_str = @sprintf("%14.6E",ρ)
  av_str = @sprintf("%14.6E",av)
  temp_str = @sprintf("%14.6E",temp)
  println(ρ_str,av_str,temp_str)
  outline = "$(newln[1:14])$(ρ_str)$(av_str)$(temp_str)\n"
  write(outfile,outline)
end

close(infile)
close(outfile)


text:
"In place version of `select`."
function select!(
    op::SelectUnion,
    C::GBVecOrMat,
    A::GBArray,
    thunk::Union{GBScalar, Nothing, Number} = nothing;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    desc = _handledescriptor(desc; in1=A)
    thunk === nothing && (thunk = C_NULL)
    accum = getaccum(accum, eltype(C))
    if thunk isa Number
        thunk = GBScalar(thunk)
    end
    libgb.GxB_Matrix_select(C, mask, accum, op, parent(A), thunk, desc)
    return C
end

function select!(
    op::Function,
    C::GBVecOrMat,
    A::GBArray,
    thunk::Union{GBScalar, Nothing, Number} = nothing;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    return select!(SelectOp(op), C, A, thunk; mask, accum, desc)
end
"""
    select(op::SelectUnion, A::GBArray; kwargs...)::GBArray
    select(op::SelectUnion, A::GBArray, thunk; kwargs...)::GBArray

Return a `GBArray` whose elements satisfy the predicate defined by `op`.
Some SelectOps may require an additional argument `thunk`, for use in comparison operations
such as `C[i,j] = A[i,j] >= thunk ? A[i,j] : nothing`, which maps to
`select(GT, A, thunk)`.

# Arguments
- `op::SelectUnion`: A select operator from the SelectOps submodule.
- `A::GBArray`: GBVector or optionally transposed GBMatrix.
- `thunk::Union{GBScalar, nothing, valid_union}`: Optional value used to evaluate `op`.

# Keywords
- `mask::Union{Nothing, GBMatrix} = nothing`: optional mask which determines the output
    pattern.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: optional binary accumulator
    operation where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc = nothing`

# Returns
- `GBArray`: The output matrix whose `eltype` is determined by `A`.
"""
function select(
    op::SelectUnion,
    A::GBArray,
    thunk::Union{GBScalar, Nothing, valid_union} = nothing;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    mask, accum = _handlenothings(mask, accum)
    C = similar(A)
    select!(op, C, A, thunk; accum, mask, desc)
    return C
end
function select(
    op::Function, A::GBArray, thunk;
    mask = nothing, accum = nothing, desc = nothing
)
    select(SelectOp(op), A, thunk; mask, accum, desc)
end

function select(
    op::Function,
    A::GBArray,
    thunk::Union{GBScalar, Nothing, Number} = nothing;
    mask = nothing,
    accum = nothing,
    desc = nothing
)
    return select(SelectOp(op), A, thunk; mask, accum, desc)
end

LinearAlgebra.tril(A::GBArray, k::Integer = 0) = select(tril, A, k)
LinearAlgebra.triu(A::GBArray, k::Integer = 0) = select(triu, A, k)
SparseArrays.dropzeros(A::GBArray) = select(nonzeros, A)


text:
function hide(x) x end

for n = 1 : hide("foo")
    println(n)
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
@testset "Histogram Construction" begin
    # Consider an image where each intensity occurs only once and vary the number
    # of bins used in the histogram in powers of two. With the exception of the
    # first bin (with index 0), all other bins should have equal counts.
    expected_counts = [2^i for i = 0:7]
    bins = [2^i for i = 8:-1:1]
    for i = 1:length(bins)
        for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
            edges, counts  = build_histogram(T.(collect(0:1/255:1)),bins[i], minval = 0, maxval = 1)
            @test length(edges) == length(counts) - 1
            @test all(counts[1:end] .== expected_counts[i]) && counts[0] == 0
            @test axes(counts) == (0:length(edges),)
        end

        # Verify that the function can also take a color image as an input.
        for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
            imgg = collect(0:1/255:1)
            img = colorview(RGB,imgg,imgg,imgg)
            edges, counts  = build_histogram(T.(img),bins[i], minval = 0, maxval = 1)
            @test length(edges) == length(counts) - 1
            @test all(counts[1:end] .== expected_counts[i]) && counts[0] == 0
            @test axes(counts) == (0:length(edges),)
        end

        # Consider also integer-valued images.
        edges, counts  = build_histogram(0:1:255,bins[i], minval = 0, maxval = 255)
        @test length(edges) == length(counts) - 1
        @test all(counts[1:end] .== expected_counts[i]) && counts[0] == 0
        @test axes(counts) == (0:length(edges),)
    end

    # Consider truncated intervals.
    for T in (Int, Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        if T == Int
            edges, counts  = build_histogram(0:1:255,4, minval = 128, maxval = 192)
            @test length(edges) == length(counts) - 1
            @test collect(counts) == [128; 16; 16; 16; 80]
            @test axes(counts) == (0:length(edges),)
        else
            img = collect(0:1/255:1)
            edges, counts  = build_histogram(T.(img),4, minval = 128/255, maxval = 192/255)
            @test length(edges) == length(counts) - 1
            @test collect(counts) == [128; 16; 16; 16; 80]
            @test axes(counts) == (0:length(edges),)
        end

        if T == Int
            edges, counts  = build_histogram(0:1:255,4, minval = 120,maxval = 140)
            @test length(edges) == length(counts) - 1
            @test collect(counts) == [120, 5, 5, 5, 121]
            @test axes(counts) == (0:length(edges),)
        else
            img = collect(0:1/255:1)
            edges, counts  = build_histogram(T.(img),4,minval = 120/255, maxval = 140/255)
            @test length(edges) == length(counts) - 1
            @test axes(counts) == (0:length(edges),)
            # Due to roundoff errors the bins are not the same as in the
            # integer case above.
            @test all([120, 4, 4, 4, 121] .<= collect(counts) .<= [120, 6, 6, 6, 121])
            @test sum(counts) == length(img)
        end
    end

    # Consider the case where the minimum and maximum values are not the start and
    # end points of the dynamic range. Because of numerical precision, the
    # results will be slightly different depending on the Image type.
    for T in (Int, Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        if T == Int
            edges, counts  = build_histogram(200:1:240,4, minval = 200, maxval = 240)

            @test length(edges) == length(counts) - 1
            @test collect(counts) == [0, 10, 10, 10, 11]
            @test axes(counts) == (0:length(edges),)

            edges, counts  = build_histogram(200:1:240,4)
            @test length(edges) == length(counts) - 1
            @test collect(counts) == [0, 10, 10, 10, 11]
            @test axes(counts) == (0:length(edges),)
        else
            img = 200/255:1/255:240/255
            edges, counts  = build_histogram(T.(img),4, minval = 200/255, maxval = 240/255)
            @test length(edges) == length(counts) - 1
            @test all([0, 9, 9, 9, 11] .<= collect(counts) .<= [0, 11, 11, 11, 11])
            @test sum(counts) == length(img)
            @test axes(counts) == (0:length(edges),)

            edges, counts  = build_histogram(T.(img),4)
            @test axes(counts) == (0:length(edges),)
            @test length(edges) == length(counts) - 1
            @test all([0, 9, 9, 9, 10] .<= collect(counts) .<= [0, 11, 11, 11, 11])
            @test sum(counts) == length(img)
        end
    end

    # Consider the effect of NaN on the histogram counts.
    for j = 2:255
        for T in (Gray{Float32}, Gray{Float64})
            img = collect(0:1/255:1)
            img[j] = NaN
            edges, counts  = build_histogram(T.(img),256, minval = 0, maxval = 1)
            target = [1 for k = 1:length(counts)-1]
            target[j] -= 1
            @test length(edges) == length(counts) - 1
            @test counts[1:end] == target && counts[0] == 0
            @test axes(counts) == (0:length(edges),)

            # Verify that the minimum value (0) and maximum value (1) is
            # determined automatically even when there are NaN values
            # between these two bounds.
            img = collect(0:1/255:1)
            img[j] = NaN
            edges, counts  = build_histogram(T.(img),256)
            target = [1 for k = 1:length(counts)-1]
            target[j] -= 1 # Because of the NaN value in the image
            @test length(edges) == length(counts) - 1
            @test counts[1:end] == target && counts[0] == 0
            @test axes(counts) == (0:length(edges),)
        end
    end
end


text:
# This file was generated, do not modify it. # hide
using PyPlot

figure(figsize=(8,6))
plot(y, ls="none", marker="o")

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Index", fontsize=14), ylabel("Salary", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-6-g1.svg")) # hide

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
#Save a (trained) neural network to an HDF5 file
function nnSave(nn::neuralNet,fileName::String)
  f = h5open(fileName,"w");
  write(f,"nNeurons",nn.nNeurons);
  write(f,"w",nn.w);
  write(f,"b",nn.b);
  close(f);
  println("Saved to: $fileName.");
end

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet


text:
using DLLBFGS
# regular pkgs
using LinearAlgebra
# JSO
using JSOSolvers, SolverTools, Logging, NLPModels, Test, LinearOperators

function tests()
    problems = [(x->2*x[1]^2 - 1.05*x[1]^4 + (x[1]^6)/6 + x[1]*x[2] + x[2]^2),
                (x->(1 - x[1])^2 + 100(x[2] - x[1]^2)^2),
                (x->(x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2),
                (x->0.26(x[1]^2 + x[2]^2) - 0.48*x[1]*x[2])]

    nlps = [ADNLPModel(i, [2.0, 3.0]) for i in problems]
    solvers = [ dllbfgs]
    @testset "DL-LBFGS" begin
      for nlp in nlps
        with_logger(NullLogger()) do
          for s in solvers
            out = s(nlp)
            @test -1e-4 ≤ out.objective ≤ 1e-4
          end
        end
      end
    end
end

tests()


text:
using LinearAlgebra

# degenerate cases (x₁ + x₂ + ⋯ + xₘ)¹
@test [multiexponents(1,1)...] == [[1]]
@test [multiexponents(2,1)...] == [[1,0],[0,1]]
@test [multiexponents(3,1)...] == [[1,0,0],[0,1,0],[0,0,1]]
@test hcat([multiexponents(10,1)...]...) == Matrix{Float64}(I, 10, 10)

# degenerate cases (x₁ + x₂ + ⋯ + xₘ)⁰
@test [multiexponents(1,0)...] == [[0]]
@test [multiexponents(2,0)...] == [[0,0]]
@test [multiexponents(3,0)...] == [[0,0,0]]
@test [multiexponents(10,0)...] == [zeros(Int, 10)]

# degenerate cases (x₁)ⁿ
@test [multiexponents(1,1)...] == [[1]]
@test [multiexponents(1,2)...] == [[2]]
@test [multiexponents(1,3)...] == [[3]]
@test [multiexponents(1,10)...] == [[10]]

# general cases
@test [multiexponents(3,2)...] == [[2,0,0],[1,1,0],[1,0,1],[0,2,0],[0,1,1],[0,0,2]]
@test [multiexponents(2,3)...] == [[3,0],[2,1],[1,2],[0,3]]
@test [multiexponents(2,2)...] == [[2,0],[1,1],[0,2]]
@test [multiexponents(3,3)...] == [[3,0,0],[2,1,0],[2,0,1],[1,2,0],[1,1,1],[1,0,2],[0,3,0],[0,2,1],[0,1,2],[0,0,3]]


text:
# simulation: DataFrame with all additional incident information created during the simulation
function prepare_simulation_data(incidents::DataFrame, districts::DataFrame)
    sim_data = DataFrame(incidentid                 = Int[],
                         location_responsible       = Int[],
                         location_dispatched_first  = Union{Missing, Int}[],
                         incident_minute            = Int[],
                         dispatch_minute_first      = Union{Missing, Int}[],
                         dispatch_minute_all        = Union{Missing, Int}[],
                         arrival_minute_first       = Union{Missing, Int}[],
                         arrival_minute_all         = Union{Missing, Int}[],
                         cars_dispatched            = Union{Missing, Int}[],
                         cars_missing               = Union{Missing, Int}[],
                         cars_location_responsible  = Union{Missing, Int}[])
    for i = 1:nrow(incidents)
        push!(sim_data, (incidentid                 = incidents[i, :incidentid],
                         location_responsible       = districts[incidents[i, :location], :location],
                         location_dispatched_first  = missing,
                         incident_minute            = floor(incidents[i, :epoch]/60),
                         dispatch_minute_first      = missing,
                         dispatch_minute_all        = missing,
                         arrival_minute_first       = missing,
                         arrival_minute_all         = missing,
                         cars_dispatched            = 0,
                         cars_missing               = incidents[i, :cars],
                         cars_location_responsible  = 0))
    end
    return sim_data
end


text:
module TestCopy

using Folds
using Folds.Testing: eval_call_data, test_with_sequential
using Test
using Transducers: DistributedEx, PreferParallel, SequentialEx, ThreadedEx

preamble = quote
    using DataFrames: DataFrame
    using StructArrays: StructVector
end

rawdata = """
copy(StructVector, ((x = x, y = x^2) for x in 1:10))
copy(StructVector, ((x = x, y = x^2) for x in 1:10 if isodd(x)))
copy(DataFrame, ((x = x, y = x^2) for x in 1:10))
copy(DataFrame, ((x = x, y = x^2) for x in 1:10 if isodd(x)))
"""

tests = eval_call_data(rawdata, preamble)
executors = [
    SequentialEx(),
    ThreadedEx(),
    ThreadedEx(basesize = 1),
    DistributedEx(),
    DistributedEx(basesize = 1),
    PreferParallel(),
    PreferParallel(basesize = 1),
]

test_with_sequential(tests, executors)

end  # module


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module LE

using ExtractMacro
using ..Interface
using ..Common

using ...RRRMC # this is silly but it's required for correct cross-linking in docstrings, apparently

export GraphLE, GraphLocalEntropy, LEenergies

import ..Interface: energy, delta_energy, neighbors, allΔE,
                    update_cache!, delta_energy_residual,
                    cenergy, distances

mutable struct GraphLE{M,γT} <: DiscrGraph{Float64}
    N::Int
    Nk::Int
    cache::LocalFields{Int}

    function GraphLE{M,γT}(N::Integer) where {M,γT}
        isa(M, Int) || throw(ArgumentError("invalid parameter M, expected Int, given: $(typeof(M))"))
        M > 2 || throw(ArgumentError("M must be greater than 2, given: $M"))
        isa(γT, Float64) || throw(ArgumentError("invalid parameter γT, expected Float64, given: $(typeof(γT))"))
        N % (M + 1) == 0 || throw(ArgumentError("N must be divisible by M+1, given: N=$N M=$M"))
        Nk = N ÷ (M+1)

        cache = LocalFields{Int}(N)
        return new{M,γT}(N, Nk, cache)
    end
end

@doc """
    GraphLE{M,γT}(N::Integer) <: DiscrGraph

An auxiliary `DiscrGraph` used to implement the interactions in the
Local Entropy Ensemble.

It is only useful when implementing other graph types; see [`GraphLocalEntropy`](@ref).
""" -> GraphLE{M,γT}(N::Integer)

GraphLE(X::GraphLE{M,oldγ}, newγ::Float64) where {M,oldγ} = GraphLE{M,newγ}(X.N)

# @generated function ΔElist{M,γT}(::Type{GraphLE{M,γT}})
#     Expr(:tuple, ntuple(d->fk(2*(d - 1 - (M-1) >>> 0x1) - iseven(M), γT), M)...)
# end
# lstind(μ̄::Int, M::Int) = (μ̄ + M-1) >>> 0x1 + 1
#
# function getk{M,γT}(X::GraphLE{M,γT}, μ̄::Int)
#     @inbounds k = ΔElist(GraphLE{M,γT})[lstind(μ̄, M)]
#     return k
# end

function energy(X::GraphLE{M,γT}, C::Config) where {M,γT}
    # @assert X.N == C.N
    @extract X : Nk cache
    @extract cache : lfields lfields_last
    @extract C : s

    n = 0
    j = 0
    for i = 1:Nk
        j += 1
        # @assert j == 1 + (i-1) * (M+1)
        jc = j
        σc = 2s[jc] - 1

        μ = 0
        for k = 2:(M+1)
            j += 1
            # @assert j == k + (i-1) * (M+1)
            σj = 2s[j] - 1
            lfields[j] = σc * σj
            μ += σj
        end
        f = σc * μ
        lfields[jc] = f
        n -= f
    end
    cache.move_last = 0
    fill!(lfields_last, 0.0)
    return n * γT
end

function kinterval(move::Integer, M::Integer)
    j0 = move - ((move-1) % (M+1))
    j1 = j0 + M
    return j0:j1
end

function update_cache!(X::GraphLE{M,γT}, C::Config, move::Int) where {M,γT}
    # @assert X.N == C.N
    # @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : Nk cache

    k = mod1(move, M+1)
    if k ≠ 1
        @inbounds σx = 2s[move] - 1
        i = (move-1) ÷ (M+1) + 1
        jc = 1 + (i-1) * (M+1)
    else
        Ux = kinterval(move, M)
    end

    @extract cache : lfields lfields_last move_last
    if move_last == move
        @inbounds begin
            #for y in neighbors(X, move)
            if k ≠ 1
                lfields[jc], lfields_last[jc] = lfields_last[jc], lfields[jc]
                lfields[move], lfields_last[move] = lfields_last[move], lfields[move]
            else
                for y in Ux
                    lfields[y], lfields_last[y] = lfields_last[y], lfields[y]
                end
            end
            #lfields[move] = -lfields[move]
            #lfields_last[move] = -lfields_last[move]
        end
        return
    end

    @inbounds begin
        if k ≠ 1
            #jc = 1 + (i-1) * (M+1)
            σc = 2s[jc] - 1
            lfc = lfields[jc]
            lfm = lfields[move]
            lfields_last[jc] = lfc
            lfields_last[move] = lfm
            lfields[jc] = lfc + 2 * (σc * σx)
            lfields[move] = -lfm
        else
            for y in Ux
                lfy = lfields[y]
                lfields_last[y] = lfy
                lfields[y] = -lfy
            end
        end
        #lfm = lfields[move]
        #lfields_last[move] = lfm
        #lfields[move] = -lfm
    end
    cache.move_last = move

    # lfields_bk = copy(lfields)
    # energy(X, C)
    # lfields_bk ≠ lfields && @show move hcat(lfields,lfields_bk) find(lfields .≠ lfields_bk)
    # @assert lfields_bk == lfields

    return
end

@inline function delta_energy(X::GraphLE{M,γT}, C::Config, move::Int) where {M,γT}
    # @assert X.N == C.N
    # @assert 1 ≤ move ≤ C.N
    @extract X : cache
    @extract cache : lfields

    @inbounds Δ = 2γT * lfields[move]
    return Δ
end

@inline function neighbors(X::GraphLE{M}, j::Int) where {M}
    r = (j-1) % (M+1)
    if r == 0
        return (j+1):(j+M)
    else
        j0 = j - r
        return j0:j0
    end
end

@generated function allΔE(::Type{GraphLE{M,γT}}) where {M,γT}
    iseven(M) ? Expr(:tuple, insert!([4*d*abs(γT) for d = 0:(M÷2)], 2, 2abs(γT))...) :
                Expr(:tuple, [2*(2d-1)*abs(γT) for d = 1:((M+1)÷2)]...)
end

# Replicate an existsing graph

mutable struct GraphLocalEntropy{M,γT,G<:AbstractGraph} <: DoubleGraph{DiscrGraph{Float64},Float64}
    N::Int
    Nk::Int
    X0::GraphLE{M,γT}
    Xc::G
    X1::Vector{G}
    Cc::Config
    C1::Vector{Config}
    function GraphLocalEntropy{M,γT,G}(N::Integer, g0::G, Gconstr, args...) where {M,γT,G}
        X0 = GraphLE{M,γT}(N)
        Nk = X0.Nk
        X1 = Array{G}(undef, M)
        Xc = g0
        for k = 1:M
            X1[k] = Gconstr(args...)
        end
        Cc = Config(Nk, init=false)
        C1 = [Config(Nk, init=false) for k = 1:M]
        return new{M,γT,G}(N, Nk, X0, Xc, X1, Cc, C1)
    end
end

"""
    GraphLocalEntropy(N::Integer, M::Integer, γ::Float64, β::Float64, Gconstr, args...) <: DoubleGraph

A `DoubleGraph` which implements a "Local Entropy" model, given any other Ising model previously defined.
This simulates a replicated system in which each replica interacts with an extra "reference"
configuration.
This kind of graph can be simulated efficiently with [`rrrMC`](@ref).

`N` is the number of spins, `M` the number of replicas, `γ` the interaction strength,
`β` the inverse temperature. `GConstr` is the (original) graph constructor, and `args` the
arguments to the contructor.

This is similar to [`GraphRobustEnsemble`](@ref RRRMC.GraphRobustEnsemble), but the reference is simulated explicitly.

See also [`LEenergies`](@ref) and [`cenergy`](@ref).
"""
function GraphLocalEntropy(Nk::Integer, M::Integer, γ::Float64, β::Float64, Gconstr, args...)
    g0 = Gconstr(args...)
    G = typeof(g0)
    return GraphLocalEntropy{M,γ/β,G}(Nk * (M+1), g0, Gconstr, args...)
end

function update_cache!(X::GraphLocalEntropy{M}, C::Config, move::Int) where {M}
    @extract X : X0 Xc X1 Cc C1
    k = mod1(move, M+1)
    i = (move - 1) ÷ (M+1) + 1

    if k == 1
        #spinflip!(Xc, Cc, i)
        spinflip!(Cc, i) # Xc cache not updated!
    else
        spinflip!(X1[k-1], C1[k-1], i)
    end

    update_cache!(X0, C, move)
end

function energy(X::GraphLocalEntropy{M}, C::Config) where {M}
    # @assert X.N == C.N
    @extract X : Nk X0 Xc X1 Cc C1
    @extract C : s

    E = energy(X0, C)

    Cc.s .= s[1:(M+1):(1 + (M+1) * (Nk-1))]
    energy(Xc, Cc)

    for k = 1:M
        C1[k].s .= s[(k+1):(M+1):(k+1 + (M+1) * (Nk-1))]
        E += energy(X1[k], C1[k])
    end

    return E
end

"""
    LEenergies(X::GraphLocalEntropy)

Returns a Vector with the individual energy (as defined by the original model)
of each replica in a [`GraphLocalEntropy`](@ref) graph.
"""
function LEenergies(X::GraphLocalEntropy{M}) where {M}
    @extract X : X1 C1
    return Float64[energy(X1[k], C1[k]) for k = 1:M]
end

function cenergy(X::GraphLocalEntropy{M}) where {M}
    @extract X : Xc Cc
    return energy(Xc, Cc)
end

function delta_energy_residual(X::GraphLocalEntropy{M}, C::Config, move::Int) where {M}
    @extract X : X1 C1

    k = mod1(move, M+1)
    k == 1 && return 0.0

    i = (move - 1) ÷ (M+1) + 1

    return delta_energy(X1[k-1], C1[k-1], i)
end

function delta_energy(X::GraphLocalEntropy, C::Config, move::Int)
    return delta_energy(X.X0, C, move) +
           delta_energy_residual(X, C, move)
end

# This is type-unstable and inefficient. On the other hand, it is
# basically only written for testing purposes...
function neighbors(X::GraphLocalEntropy{M}, i::Int) where {M}
    @extract X : X0 X1
    jts = neighbors(X0, i)

    k = mod1(i, M+1)

    k == 1 && return jts

    i1 = (i - 1) ÷ (M+1) + 1

    is = neighbors(X1[k-1], i1)

    return tuple(jts..., map(i->((i-1)*(M+1)+k), is)...)
end

function distances(X::GraphLocalEntropy{M}) where {M}
    @extract X : C1
    # @extract C : s
    # @assert Nk == C.N ÷ (M+1)

    # rs = BitVector[s[(k-1)*(M+1) + (2:(M+1):end)] for k = 1:M]

    # return [sum(s1 .⊻ s2) for s1 in rs, s2 in rs]
    return [sum(C1[k1].s .⊻ C1[k2].s) for k1 = 1:M, k2 = 1:M]
end

end # module


text:
"""
    symplectic_integrate(x₀, p₀, Λ, U, δUδx, N=50, ϵ=0.1, progress=false)

Do a symplectic integration of the potential energy `U` (with gradient
`δUδx`) starting from point `x₀` with momentum `p₀` and mass matrix
`Λ`. The number of steps is `N` and the step size `ϵ`.

Returns `ΔH, xᵢ, pᵢ` corresponding to change in Hamiltonian, and final
position and momenta. If `history_keys` is specified a history of
requested variables throughout each step is also returned.

"""
function symplectic_integrate(
    x₀::AbstractVector{T}, p₀, Λ, U, δUδx=x->gradient(U,x)[1];
    N=50, ϵ=T(0.1), progress=false, history_keys=nothing
) where {T}

    xᵢ, pᵢ = x₀, p₀
    δUδxᵢ = δUδx(xᵢ)
    H(x,p) = U(x) - p⋅(Λ\p)/2

    history = []

    @showprogress (progress ? 1 : Inf) "Symplectic Integration: " for i=1:N
        xᵢ₊₁    = xᵢ - T(ϵ) * (Λ \ (pᵢ - T(ϵ)/2 * δUδxᵢ))
        δUδxᵢ₊₁ = δUδx(xᵢ₊₁)
        pᵢ₊₁    = pᵢ - T(ϵ)/2 * (δUδxᵢ₊₁ + δUδxᵢ)
        xᵢ, pᵢ, δUδxᵢ = xᵢ₊₁, pᵢ₊₁, δUδxᵢ₊₁

        if !isnothing(history_keys)
            historyᵢ = (;i, x=xᵢ, p=pᵢ, δUδx=δUδxᵢ₊₁, H=(haskey(history_keys,:H) ? H(xᵢ,pᵢ) : nothing))
            push!(history, select(historyᵢ, history_keys))
        end

    end

    ΔH = H(xᵢ,pᵢ) - H(x₀,p₀)

    if isnothing(history)
        return ΔH, xᵢ, pᵢ
    else
        return ΔH, xᵢ, pᵢ, history
    end

end

@doc doc"""
    grid_and_sample(lnP::Function; range::NamedTuple; progress=false, nsamples=1)

Interpolate the log pdf `lnP` with support on `range`, and return  the
integrated log pdf as well `nsamples` samples (drawn via inverse transform
sampling)

`lnP` should either accept a NamedTuple argument and `range` should be a
NamedTuple mapping those same names to `range` objects specifying where to
evaluate `lnP`, e.g.:

```julia
grid_and_sample(nt->-(nt.x^2+nt.y^2)/2, (x=range(-3,3,length=100),y=range(-3,3,length=100)))
```

or `lnP` should accept a single scalar argument and `range` should be directly
the range for this variable:

```julia
grid_and_sample(x->-x^2/2, range(-3,3,length=100))
```

The return value is `(lnP, samples, Px)` where `lnP` is an interpolated/smoothed
log PDF which can be evaluated anywhere within the original range, `Px` are
sampled points of the original PDF, and `samples` is a NamedTuple giving the
Monte-Carlo samples of each of the parameters.

(Note: only 1D sampling is currently implemented, but 2D like in the example
above is planned)
"""
function grid_and_sample(lnP::Function, range::AbstractVector; progress=false, kwargs...)
    lnPs = @showprogress (progress ? 1 : Inf) "Grid Sample: " map(lnP, range)
    grid_and_sample(lnPs, range; progress=progress, kwargs...)
end

function grid_and_sample(lnPs::Vector{<:BatchedReal}, xs::AbstractVector; kwargs...)
    batches = [grid_and_sample(batch_index.(lnPs,i), xs; kwargs...) for i=1:batch_length(lnPs[1])]
    ((batch(getindex.(batches,i)) for i=1:3)...,)
end

function grid_and_sample(lnPs::Vector, xs::AbstractVector; progress=false, nsamples=1, span=0.25, require_convex=false)

    # trim leading/trailing zero-probability regions
    support = findnext(isfinite,lnPs,1):findprev(isfinite,lnPs,length(lnPs))
    xs = xs[support]
    lnPs = lnPs[support]

    if require_convex
        support = longest_run_of_trues(finite_second_derivative(lnPs) .< 0)
        xs = xs[support]
        lnPs = lnPs[support]
    end

    # interpolate PDF
    xmin, xmax = first(xs), last(xs)
    lnPs = lnPs .- maximum(lnPs)
    ilnP = loess(xs, lnPs, span=span)

    # normalize the PDF. note the smoothing is done of the log PDF.
    cdf(x) = quadgk(nan2zero∘exp∘ilnP,xmin,x,rtol=1e-3)[1]
    logA = nan2zero(log(cdf(xmax)))
    lnPs = (ilnP.ys .-= logA)
    ilnP.bs[:,1] .-= logA

    # draw samples via inverse transform sampling
    θsamples = @showprogress (progress ? 1 : Inf) map(1:nsamples) do i
        r = rand()
        if (cdf(xmin)-r)*(cdf(xmax)-r) >= 0
            first(lnPs) > last(lnPs) ? xmin : xmax
        else
            fzero(x->cdf(x)-r, xmin, xmax, xatol=(xmax-xmin)*1e-3)
        end
    end

    (nsamples==1 ? θsamples[1] : θsamples), ilnP, lnPs

end

function grid_and_sample(lnP::Function, range::NamedTuple{S,<:NTuple{1}}; kwargs...) where {S}
    NamedTuple{S}.(Ref.(grid_and_sample(x -> lnP(NamedTuple{S}(x)), first(range); kwargs...)))
end

# allow more convenient evaluation of Loess-interpolated functions
(m::Loess.LoessModel)(x) = Loess.predict(m,x)

@doc doc"""
    sample_joint(ds::DataSet; kwargs...)

Sample the joint posterior, $\mathcal{P}(f,\phi,\theta\,|\,d)$.

Keyword arguments:

* `nsamps_per_chain` — The number of samples per chain.
* `nchains = 1` — Number of chains in parallel.
* `nchunk = 1` — Number of steps between parallel chain communication.
* `nsavemaps = 1` — Number of steps in between saving maps into chain.
* `nburnin_always_accept = 0` — Number of steps at the beginning of
  the chain to always accept HMC steps regardless of integration
  error.
* `nburnin_fixθ = 0` — Number of steps at the beginning of the chain
  before starting to sample `θ`.
* `Nϕ = :qe` — Noise to use in the initial approximation to the
  Hessian. Can give `:qe` to use the quadratic estimate noise.
* `chains = nothing` — `nothing` to start a new chain; the return
  value from a previous call to `sample_joint` to resume those chains;
  `:resume` to resume chains from a file given by `filename`
* `θrange` — Range and density to grid sample parameters as a
  NamedTuple, e.g. `(Aϕ=range(0.7,1.3,length=20),)`.
* `θstart` — Starting values of parameters as a NamedTuple, e.g.
  `(Aϕ=1.2,)`, or nothing to randomly sample from θrange
* `ϕstart` — Starting `ϕ`, either a `Field` object, `:quasi_sample`,
  or `:best_fit`
* `metadata` — Does nothing, but is saved into the chain file
* `nhmc = 1` — Number of HMC passes per `ϕ` Gibbs step.
* `symp_kwargs = fill((N=25, ϵ=0.01), nhmc)` — an array of NamedTupe
  kwargs to pass to [`symplectic_integrate`](@ref). E.g.
  `[(N=50,ϵ=0.1),(N=25,ϵ=0.01)]` would do 50 large steps then 25
  smaller steps per each Gibbs pass. If specified, `nhmc` is ignored.
* `wf_kwargs` — Keyword arguments to pass to [`argmaxf_lnP`](@ref) in
  the Wiener Filter Gibbs step.
* `MAP_kwargs` — Keyword arguments to pass to [`MAP_joint`](@ref) when
  computing the starting point.
"""
function sample_joint(
    ds :: DataSet;
    gibbs_initializers = [
        gibbs_initialize_θ!,
        gibbs_initialize_ϕ!,
        gibbs_initialize_f!
    ],
    gibbs_samplers = [
        gibbs_sample_f!,
        gibbs_mix!,
        gibbs_sample_ϕ!,
        gibbs_unmix!,
        gibbs_postprocess!
    ],
    nsamps_per_chain,
    nchains = nworkers(),
    nfilewrite = 5,
    nsavemaps = 1,
    filename = nothing,
    ϕstart = :prior,
    θstart = :prior,
    θrange = NamedTuple(),
    pmap = (myid() in workers() ? map : pmap),
    pool = default_worker_pool(),
    conjgrad_kwargs = (tol=1e-1, nsteps=500),
    preconditioner = :diag,
    nhmc = 1,
    symp_kwargs = fill((N=25, ϵ=0.01), nhmc),
    MAP_kwargs = (nsteps=40,),
    metadata = nothing,
    progress = false,
    storage = basetype(fieldinfo(ds.d).M),
    grid_and_sample_kwargs = (;),
    kwargs...
)

    # rundat is a Dict with all the args and kwargs minus a few removed ones
    rundat = merge!(foldl(delete!, (:ds, :gibbs_initializers, :gibbs_samplers, :kwargs, :pmap, :pool), init=Base.@locals()), kwargs)
    rundat[:Nbatch] = batch_length(ds.d)

    # dont adapt things passed in kwargs when we adapt the state dict
    _adapt(storage, state) = Dict(k => (haskey(rundat,k) ? v : adapt(storage, v)) for (k,v) in state)

    function filter_for_saving(state, step)
        Dict(k=>v for (k,v) in state if (!(haskey(rundat,k)) && !(k in (:pbar_dict, :timer)) && (step == 1 || (step % nsavemaps) == 0 || !isa(v,Field))))
    end

    # validate arguments
    if !(progress in [false,:summary,:verbose])
        error("`progress` should be one of [false, :summary, :verbose]")
    end
    if (filename!=nothing && splitext(filename)[2]!=".jld2")
        error("Chain filename '$filename' should have '.jld2' extension.")
    end

    # seed
    @everywhere @eval CMBLensing seed!(global_rng_for($storage))

    # distributed dataset object to workers once
    set_distributed_dataset(ds)

    # initialize chains
    if (filename != nothing) && isfile(filename)

        @info "Resuming chain at $filename"
        local chunks_index, prev_chunks
        jldopen(filename,"r") do io
            chunks_index = maximum([parse(Int,k[8:end]) for k in keys(io) if startswith(k,"chunks_")]) + 1
            states = last.(read(io, "chunks_$(chunks_index-1)"))
            merge!.(states, Ref(rundat))
            @unpack step = states[1]
            chain_chunks = map(copy, repeated([], nchains))
        end

    else

        chunks_index = step = 1

        states = pmap(map(copy,repeated(rundat,nchains))) do state
            state = _adapt(storage, state)
            for gibbs_initialize! in gibbs_initializers
                gibbs_initialize!(state, get_distributed_dataset())
            end
            @pack! state = step
            _adapt(Array, state)
        end
        chain_chunks = [Any[filter_for_saving(state, step)] for state in states]

        if filename != nothing
            save(filename, "rundat", cpu(rundat))
        end

    end

    # setup progressbar
    setindex!.(states, copy.(Ref(OrderedDict{String,Any}("step"=>step))), :pbar_dict)
    if progress == :summary
        pbar = Progress(nsamps_per_chain-step, dt=0, desc="Gibbs chain: ")
        ProgressMeter.update!(pbar, showvalues=[("step", step)])
        @everywhere $reset_timer!()
    end

    # start sampling
    for step in (step+1):nsamps_per_chain

        setindex!.(states, step, :step)

        state₁, = states = pmap(pool, states) do state

            state = @⌛ _adapt(storage, state)
            @unpack step, pbar_dict = state

            timing = @⌛ "Gibbs passes" map(gibbs_samplers) do gibbs_sample!
                @elapsed gibbs_sample!(state, get_distributed_dataset())
            end

            state = @⌛ _adapt(Array, state)

            timer = get_defaulttimer()
            @pack! state = timing, timer

            state

        end

        push!.(chain_chunks, filter_for_saving.(states, step))

        if (filename != nothing) && ((step % nfilewrite) == 0)
            @⌛ jldopen(filename,"a+") do io
                wsession = JLDWriteSession()
                write(io, "chunks_$chunks_index", chain_chunks, wsession)
            end
            chunks_index += 1
            empty!.(chain_chunks)
        end

        if progress == :summary
            print("\033[2J")
            next!(pbar, showvalues = [("step",step), ("timing",state₁[:timing]), collect(state₁[:pbar_dict])...])
            merge!(state₁[:timer].inner_timers, get_defaulttimer().inner_timers)
            print(state₁[:timer])
            flush(stdout)
        end

    end

    Chains(chain_chunks)

end

lnP_arg_names(p::Union{Int,Symbol}, ds::DataSet) = lnP_arg_names(Val(p), ds)
lnP_arg_names(::Val{0},    ds::DataSet) = (:f,  :ϕ,  :θ)
lnP_arg_names(::Val{:mix}, ds::DataSet) = (:f°, :ϕ°, :θ)

## initialization

function gibbs_initialize_θ!(state, ds::DataSet)
    @unpack θstart, θrange, nchains, Nbatch = state
    θ = @match θstart begin
        :prior => map(range->batch((first(range) .+ rand(Nbatch) .* (last(range) - first(range)))...), θrange)
        (_::NamedTuple) => θstart
        _ => throw(ArgumentError(θstart))
    end
    lnPθ = map(_->missing, θrange)
    @pack! state = θ, lnPθ
end

function gibbs_initialize_f!(state, ds::DataSet)
    f = missing
    @pack! state = f
end

function gibbs_initialize_ϕ!(state, ds::DataSet)
    @unpack ϕstart, θ, nchains, Nbatch = state
    ϕ = @match ϕstart begin
        :prior     => simulate(ds.Cϕ(θ); Nbatch)
        0          => zero(diag(ds.Cϕ)) * batch(ones(Int,Nbatch)...)
        (_::Field) => ϕstart
        (:quasi_sample|:best_fit) => MAP_joint(
            adapt(storage, ds(θstart)),
            progress = (progress==:verbose ? :summary : false),
            Nϕ = adapt(storage,Nϕ),
            quasi_sample = (ϕstart==:quasi_sample); MAP_kwargs...
        ).ϕ
        _ => throw(ArgumentError(ϕstart))
    end
    @pack! state = ϕ
end

## gibbs passes

@⌛ function gibbs_sample_f!(state, ds::DataSet)
    @unpack f, progress, preconditioner, conjgrad_kwargs = state
    z_other = delete(select(state, lnP_arg_names(0, ds)), :f)
    f = argmaxf_lnP(
        z_other..., ds;
        which = :sample,
        guess = ismissing(f) ? nothing : f,
        preconditioner = preconditioner,
        conjgrad_kwargs = (progress=(progress==:verbose), conjgrad_kwargs...)
    )
    @pack! state = f
end

@⌛ function gibbs_sample_ϕ!(state, ds::DataSet)
    @unpack f°, ϕ°, θ, symp_kwargs, progress, step, nburnin_always_accept = state
    U = ϕ° -> lnP(:mix, f°, ϕ°, θ, ds)
    ϕ°, ΔH, accept = hmc_step(U, ϕ°, mass_matrix_ϕ(θ,ds); symp_kwargs, progress, always_accept=(step<nburnin_always_accept))
    @pack! state = ϕ°, ΔH, accept
end

function hmc_step(U::Function, x, Λ; symp_kwargs, progress, always_accept)
    local ΔH, accept
    for kwargs in symp_kwargs
        p = simulate(Λ)
        (ΔH, xtest) = symplectic_integrate(
            x, p, Λ, U;
            progress = (progress==:verbose),
            kwargs...
        )
        accept = batch(@. always_accept | (log(rand()) < $unbatch(ΔH)))
        @. x = accept * xtest + (1 - accept) * x
    end
    x, ΔH, accept
end

@⌛ function mass_matrix_ϕ(θ, ds)
    @unpack G, Cϕ, Nϕ = ds(θ)
    pinv(G)^2 * (pinv(Cϕ) + pinv(Nϕ))
end

function gibbs_sample_slice_θ!(k::Symbol)
    @⌛ function gibbs_sample_slice_θ!(state, ds::DataSet)
        @unpack θ, θrange, lnPθ, pbar_dict, progress, grid_and_sample_kwargs = state
        z_other = delete(select(state, lnP_arg_names(:mix, ds)), :θ)
        θₖ, lnPθₖ = grid_and_sample(θₖ -> lnP(:mix, z_other..., @set(θ[k]=θₖ), ds), cpu(θrange[k]); progress=(progress==:verbose), grid_and_sample_kwargs...)
        @set! θ[k] = θₖ
        @set! lnPθ[k] = lnPθₖ
        pbar_dict[string(k)] = string_trunc(θₖ)
        @pack! state = θ, lnPθ
    end
end

@⌛ function gibbs_mix!(state, ds::DataSet)
    @unpack f, ϕ, θ = state
    f°, ϕ° = mix(f, ϕ, θ, ds)
    @pack! state = f°, ϕ°
end

@⌛ function gibbs_unmix!(state, ds::DataSet)
    @unpack f°, ϕ°, θ = state
    f, ϕ = unmix(f°, ϕ°, θ, ds)
    @pack! state = f, ϕ
end

## postprocessing

@⌛ function gibbs_postprocess!(state, ds::DataSet)
    @unpack f, ϕ, θ, pbar_dict, ΔH = state
    lnP = pbar_dict["lnP"] = CMBLensing.lnP(0, select(state, lnP_arg_names(0, ds))..., ds)
    pbar_dict["ΔH"] = ΔH
    f̃ = ds.L(ϕ) * f
    @pack! state = f̃, lnP
end

## util

function once_every(n)
    function (gibbs_sample!)
        function (state, ds)
            if iszero(state[:step] % n)
                gibbs_sample!(state, ds)
            end
        end
    end
end

function start_after_burnin(n)
    function (gibbs_sample!)
        function (state, ds)
            if state[:step] > n
                gibbs_sample!(state, ds)
            end
        end
    end
end


text:
@unix_only libpath = replace(@__FILE__, "src/trace/taskcopy.jl", "deps/")
@windows_only libpath = replace(@__FILE__, "src\\trace\\taskcopy.jl", "deps\\")

if !(libpath in Base.DL_LOAD_PATH)
  push!(Base.DL_LOAD_PATH, libpath)
end
# println(libpath)

function sweepandmark(t::Task)
  s :: ObjectIdDict = t.storage
  for k in keys(s)
    # All copy-on-write family data structures stores data in task_local_storage in
    # the form of (taskid::Task, data::Any). Setting taskid to nothing will triger
    #  copy-on-write for the parent task (which is copied).
    if isa(s[k], Tuple{Union{Void,Task}, Any})
      _, d = s[k]
      s[k] = (nothing, d)
    end
  end
end

function Base.copy(t::Task)
  t.state != :runnable && t.state != :done &&
    error("Only runnable or finished tasks can be copied.")
  newt = ccall((:jl_clone_task, "libtask"), Any, (Any,), t)::Task
  if t.storage != nothing
    newt.storage = copy(t.storage)
    sweepandmark(t)
  else
    newt.storage = nothing
  end
  newt.code = t.code
  newt.state = t.state
  newt.result = t.result
  newt.parent = t.parent
  newt
end

olderr = STDERR
# supress warning message for re-defining Base.produce
redirect_stderr()
function Base.produce(v)
  #### un-optimized version
  #q = current_task().consumers
  #t = shift!(q.waitq)
  #empty = isempty(q.waitq)
  ct = current_task()
  local empty, t, q
  while true
    q = ct.consumers
    if isa(q,Task)
      t = q
      ct.consumers = nothing
      empty = true
      break
    elseif isa(q,Condition) && !isempty(q.waitq)
      t = shift!(q.waitq)
      empty = isempty(q.waitq)
      break
    end
    wait()
  end

  t.state = :runnable
  if empty
    if isempty(Base.Workqueue)
      yieldto(t, v)
    else
      schedule_and_wait(t, v)
    end
    ct = current_task() # When a task is copied, ct should be updated to new task ID.
    while true
      # wait until there are more consumers
      q = ct.consumers
      if isa(q,Task)
        return q.result
      elseif isa(q,Condition) && !isempty(q.waitq)
        return q.waitq[1].result
      end
      wait()
    end
  else
    schedule(t, v)
    # make sure `t` runs before us. otherwise, the producer might
    # finish before `t` runs again, causing it to see the producer
    # as done, causing done(::Task, _) to miss the value `v`.
    # see issue #7727
    yield()
    return q.waitq[1].result
  end
end
redirect_stderr(olderr)


text:
@testset "True and False" begin
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Int}()))
    config = MOIT.Config()

    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> (MOIU.mock_optimize!(mock, [1])),
        (mock::MOIU.MockOptimizer) -> (MOIU.mock_optimize!(mock, [0])),
    )
    COIDT.truefalsetest(mock, config)
end

function true_test(model::MOI.ModelLike, config::MOIT.Config)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(model, MOI.VectorOfVariables, CP.True)

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(
        model,
        _vov([x1, x1]),
        CP.IfThenElse(CP.True(), MOI.EqualTo(1), MOI.EqualTo(0)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, c1)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
    end
end

function false_test(model::MOI.ModelLike, config::MOIT.Config)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(model, MOI.VectorOfVariables, CP.False)

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(
        model,
        _vov([x1, x1]),
        CP.IfThenElse(CP.False(), MOI.EqualTo(1), MOI.EqualTo(0)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, c1)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 0
    end
end

const truefalsetests = Dict("true" => true_test, "false" => false_test)

MOIT.@moitestset truefalse


text: #  ArgParse

max_stars_repo_name: Greg-R/spiavalonfpga

max_issues_repo_name: Greg-R/spiavalonfpga

max_forks_repo_name: Greg-R/spiavalonfpga


text:
module Grids
using ..Topologies

export DiscontinuousSpectralElementGrid, AbstractGrid
export dofs_per_element, arraytype, dimensionality, polynomialorder
import Canary

abstract type AbstractGrid{FloatType, dim, polynomialorder, numberofDOFs,
                         DeviceArray} end

dofs_per_element(::AbstractGrid{T, D, N, Np}) where {T, D, N, Np} = Np

polynomialorder(::AbstractGrid{T, dim, N}) where {T, dim, N} = N

dimensionality(::AbstractGrid{T, dim}) where {T, dim} = dim

arraytype(::AbstractGrid{T, D, N, Np, DA}) where {T, D, N, Np, DA} = DA

# {{{
const _nvgeo = 14
const _ξx, _ηx, _ζx, _ξy, _ηy, _ζy, _ξz, _ηz, _ζz, _M, _MI,
       _x, _y, _z = 1:_nvgeo
const vgeoid = (ξxid = _ξx, ηxid = _ηx, ζxid = _ζx,
                ξyid = _ξy, ηyid = _ηy, ζyid = _ζy,
                ξzid = _ξz, ηzid = _ηz, ζzid = _ζz,
                Mid  = _M , MIid = _MI,
                xid  = _x , yid  = _y , zid  = _z)

const _nsgeo = 5
const _nx, _ny, _nz, _sM, _vMI = 1:_nsgeo
const sgeoid = (nxid = _nx, nyid = _ny, nzid = _nz, sMid = _sM,
                vMIid = _vMI)
# }}}

"""
    DiscontinuousSpectralElementGrid(topology; FloatType, DeviceArray,
                                     polynomialorder,
                                     meshwarp = (x...)->identity(x))

Generate a discontinuous spectral element (tensor product,
Legendre-Gauss-Lobatto) grid/mesh from a `topology`, where the order of the
elements is given by `polynomialorder`. `DeviceArray` gives the array type used
to store the data (`CuArray` or `Array`), and the coordinate points will be of
`FloatType`.

The optional `meshwarp` function allows the coordinate points to be warped after
the mesh is created; the mesh degrees of freedom are orginally assigned using a
trilinear blend of the element corner locations.
"""
struct DiscontinuousSpectralElementGrid{T, dim, N, Np, DA,
                                        DAT2, DAT3, DAT4, DAI1, DAI2, DAI3,
                                        TOP
                                       } <: AbstractGrid{T, dim, N, Np, DA }
  "mesh topology"
  topology::TOP

  "volume metric terms"
  vgeo::DAT3

  "surface metric terms"
  sgeo::DAT4

  "element to boundary condition map"
  elemtobndy::DAI2

  "volume DOF to element minus side map"
  vmapM::DAI3

  "volume DOF to element plus side map"
  vmapP::DAI3

  "list of elements that need to be communicated (in neighbors order)"
  sendelems::DAI1

  "1-D derivative operator on the device"
  D::DAT2

  function DiscontinuousSpectralElementGrid(topology::AbstractTopology{dim};
                                            FloatType = nothing,
                                            DeviceArray = nothing,
                                            polynomialorder = nothing,
                                            meshwarp::Function =
                                            (x...)->identity(x)) where dim
    @assert FloatType != nothing
    @assert DeviceArray != nothing
    @assert polynomialorder != nothing

    N = polynomialorder
    (ξ, ω) = Canary.lglpoints(FloatType, N)
    D = Canary.spectralderivative(ξ)

    (vmapM, vmapP) = mappings(N, topology.elemtoelem, topology.elemtoface,
                              topology.elemtoordr)

    (vgeo, sgeo) = computegeometry(topology, D, ξ, ω, meshwarp, vmapM)
    Np = (N+1)^dim
    @assert Np == size(vgeo, 1)

     # Create arrays on the device
     vgeo = DeviceArray(vgeo)
     sgeo = DeviceArray(sgeo)
     elemtobndy = DeviceArray(topology.elemtobndy)
     vmapM = DeviceArray(vmapM)
     vmapP = DeviceArray(vmapP)
     sendelems = DeviceArray(topology.sendelems)
     D = DeviceArray(D)

     # FIXME: There has got to be a better way!
     DAT2 = typeof(D)
     DAT3 = typeof(vgeo)
     DAT4 = typeof(sgeo)
     DAI1 = typeof(sendelems)
     DAI2 = typeof(elemtobndy)
     DAI3 = typeof(vmapM)
     TOP = typeof(topology)

    new{FloatType, dim, N, Np, DeviceArray, DAT2, DAT3, DAT4, DAI1, DAI2, DAI3,
        TOP
       }(topology, vgeo, sgeo, elemtobndy, vmapM, vmapP, sendelems, D)
  end
end

function Base.getproperty(G::DiscontinuousSpectralElementGrid, s::Symbol)
  if s ∈ keys(vgeoid)
    vgeoid[s]
  elseif s ∈ keys(sgeoid)
    sgeoid[s]
  else
    getfield(G, s)
  end
end

function Base.propertynames(G::DiscontinuousSpectralElementGrid)
  (fieldnames(DiscontinuousSpectralElementGrid)...,
   keys(vgeoid)..., keys(sgeoid)...)
end

# {{{ mappings
"""
    mappings(N, elemtoelem, elemtoface, elemtoordr)

This function takes in a polynomial order `N` and parts of a topology (as
returned from `connectmesh`) and returns index mappings for the element surface
flux computation.  The returned `Tuple` contains:

 - `vmapM` an array of linear indices into the volume degrees of freedom where
   `vmapM[:,f,e]` are the degrees of freedom indices for face `f` of element
    `e`.

 - `vmapP` an array of linear indices into the volume degrees of freedom where
   `vmapP[:,f,e]` are the degrees of freedom indices for the face neighboring
   face `f` of element `e`.
"""
function mappings(N, elemtoelem, elemtoface, elemtoordr)
  nface, nelem = size(elemtoelem)

  d = div(nface, 2)
  Np, Nfp = (N+1)^d, (N+1)^(d-1)

  p = reshape(1:Np, ntuple(j->N+1, d))
  fd(f) =   div(f-1,2)+1
  fe(f) = N*mod(f-1,2)+1
  fmask = hcat((p[ntuple(j->(j==fd(f)) ? (fe(f):fe(f)) : (:), d)...][:]
                for f=1:nface)...)
  inds = LinearIndices(ntuple(j->N+1, d-1))

  vmapM = similar(elemtoelem, Nfp, nface, nelem)
  vmapP = similar(elemtoelem, Nfp, nface, nelem)

  for e1 = 1:nelem, f1 = 1:nface
    e2 = elemtoelem[f1,e1]
    f2 = elemtoface[f1,e1]
    o2 = elemtoordr[f1,e1]

    vmapM[:,f1,e1] .= Np*(e1-1) .+ fmask[:,f1]

    if o2 == 1
      vmapP[:,f1,e1] .= Np*(e2-1) .+ fmask[:,f2]
    elseif d == 3 && o2 == 3
      n = 1
      @inbounds for j = 1:N+1, i = N+1:-1:1
        vmapP[n,f1,e1] = Np*(e2-1) + fmask[inds[i,j],f2]
        n+=1
      end
    else
      error("Orientation '$o2' with dim '$d' not supported yet")
    end
  end

  (vmapM, vmapP)
end
# }}}

# {{{ compute geometry
function computegeometry(topology::AbstractTopology{dim}, D, ξ, ω, meshwarp,
                         vmapM) where {dim}
  # Compute metric terms
  Nq = size(D, 1)
  DFloat = eltype(D)

  (nface, nelem) = size(topology.elemtoelem)

  # crd = creategrid(Val(dim), elemtocoord(topology), ξ)

  vgeo = zeros(DFloat, Nq^dim, _nvgeo, nelem)
  sgeo = zeros(DFloat, _nsgeo, Nq^(dim-1), nface, nelem)

  (ξx, ηx, ζx, ξy, ηy, ζy, ξz, ηz, ζz, MJ, MJI, x, y, z) =
      ntuple(j->(@view vgeo[:, j, :]), _nvgeo)
  J = similar(x)
  (nx, ny, nz, sMJ, vMJI) = ntuple(j->(@view sgeo[ j, :, :, :]), _nsgeo)
  sJ = similar(sMJ)

  X = ntuple(j->(@view vgeo[:, _x+j-1, :]), dim)
  Canary.creategrid!(X..., topology.elemtocoord, ξ)

  @inbounds for j = 1:length(x)
    (x[j], y[j], z[j]) = meshwarp(x[j], y[j], z[j])
  end

  # Compute the metric terms
  if dim == 1
    Canary.computemetric!(x, J, ξx, sJ, nx, D)
  elseif dim == 2
    Canary.computemetric!(x, y, J, ξx, ηx, ξy, ηy, sJ, nx, ny, D)
  elseif dim == 3
    Canary.computemetric!(x, y, z, J, ξx, ηx, ζx, ξy, ηy, ζy, ξz, ηz, ζz, sJ,
                   nx, ny, nz, D)
  end

  M = kron(1, ntuple(j->ω, dim)...)
  MJ .= M .* J
  MJI .= 1 ./ MJ
  vMJI .= MJI[vmapM]

  sM = dim > 1 ? kron(1, ntuple(j->ω, dim-1)...) : one(DFloat)
  sMJ .= sM .* sJ

  (vgeo, sgeo)
end
# }}}

end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
# import Particles package
using Particles

# parameters
Δt = 0.1    # time step size
nt = 100    # number of time steps
ns = 5      # number of samples

# random initial conditions
x₀ = hcat(rand(ns) .* 2π, rand(ns) .* 2 .- 1)'

# vector field
function pendulum!(ẋ, x)
    ẋ[1] = x[2]
    ẋ[2] = - sin(x[1])
end

# create an Equation instance
equ = Equation(pendulum!, x₀)

# create an Integrator instance
int = ExplicitEuler(equ, Δt)

# create Simulation instance
sim = Simulation(equ, int, nt)

# run simulation
x = run!(sim)

# load Plots package
using Plots

# set plot ranges
xlim = (-π, +π)
ylim = (-3, +3)

# plot energy contour lines
plt = contour(
    LinRange(xlim..., 100),
    LinRange(ylim..., 100),
    (x,v) -> v^2 / 2  + 1 - cos(x),
    xlim = xlim,
    ylim = ylim,
    title = "Pendulum",
    legend = false,
    size = (800, 600)
)

# plot solutions
for i in axes(x,2)
    scatter!(plt, (x[1,i,:] .+ π) .% 2π .- π, x[2,i,:], marker=5)
end

# save figure to file
savefig("pendulum.png")

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
using ForwardDiff
using ForwardDiff

import NLPModels.obj, NLPModels.grad, NLPModels.grad!, NLPModels.cons, NLPModels.cons!
import NLPModels.jac, NLPModels.jprod, NLPModels.jprod!, NLPModels.jtprod
import NLPModels.jtprod!, NLPModels.hess, NLPModels.hprod, NLPModels.hprod!
import NLPModels.increment!

export HOADNLPModel, obj, grad, grad!, cons, cons!, jac, jprod,
       jprod!, jtprod, jtprod!, hess, hprod, hprod!, ∇f³xuv, ∇f³xuv2

"""HOADNLPModel is an AbstractNLPModel using ForwardDiff to compute the
derivatives. HOHOADNLPModel uses ForwardDiff to compute 3rd order information.
This interface is very similar to HOADNLPModel with few additional functions which
are documented here.
For further documentation, refer to the HOADNLPModel.
"""
mutable struct HOADNLPModel <: AbstractNLPModel
  meta     :: NLPModelMeta

  counters :: Counters

  # Functions
  f :: Function
  c :: Function

  # AD Packages
  AD :: Any
end

function HOADNLPModel(f::Function, x0::AbstractVector; y0::AbstractVector = eltype(x0)[],
                      lvar::AbstractVector = eltype(x0)[], uvar::AbstractVector = eltype(x0)[],
                      lcon::AbstractVector = eltype(x0)[], ucon::AbstractVector = eltype(x0)[],
                      c::Function = (args...)->throw(NotImplementedError("cons")),
                      name::String = "Generic", lin::AbstractVector{Int}=Int[],
                      AD :: Any = ForwardDiff)

  nvar = length(x0)
  length(lvar) == 0 && (lvar = -Inf*ones(nvar))
  length(uvar) == 0 && (uvar =  Inf*ones(nvar))
  ncon = maximum([length(lcon); length(ucon); length(y0)])

  A = AD.hessian(f, x0)
  for i = 1:ncon
    A += AD.hessian(x->c(x)[i], x0) * (-1)^i
  end
  nnzh = nvar * (nvar + 1) / 2
  nnzj = 0

  if ncon > 0
    length(lcon) == 0 && (lcon = -Inf*ones(ncon))
    length(ucon) == 0 && (ucon =  Inf*ones(ncon))
    length(y0) == 0   && (y0 = zeros(ncon))
    nnzj = nvar * ncon
  end
  nln = setdiff(1:ncon, lin)

  meta = NLPModelMeta(nvar, x0=x0, lvar=lvar, uvar=uvar, ncon=ncon, y0=y0,
    lcon=lcon, ucon=ucon, nnzj=nnzj, nnzh=nnzh, lin=lin, nln=nln, minimize=true,
    islp=false, name=name)

  return HOADNLPModel(meta, Counters(), f, c, AD)
end

function obj(nlp :: HOADNLPModel, x :: AbstractVector)
  # increment!(nlp, :neval_obj)
  return nlp.f(x)
end

function grad(nlp :: HOADNLPModel, x :: AbstractVector)
  # increment!(nlp, :neval_grad)
  return ForwardDiff.gradient(nlp.f, x)
end

function grad!(nlp :: HOADNLPModel, x :: AbstractVector, g :: AbstractVector)
  # increment!(nlp, :neval_grad)
  ForwardDiff.gradient!(view(g, 1:length(x)), nlp.f, x)
  return g
end

function cons(nlp :: HOADNLPModel, x :: AbstractVector)
  # increment!(nlp, :neval_cons)
  return nlp.c(x)
end

function cons!(nlp :: HOADNLPModel, x :: AbstractVector, c :: AbstractVector)
  increment!(nlp, :neval_cons)
  c[1:nlp.meta.ncon] = nlp.c(x)
  return c
end

function jac(nlp :: HOADNLPModel, x :: AbstractVector)
  # increment!(nlp, :neval_jac)
  return ForwardDiff.jacobian(nlp.c, x)
end

function jac_coord(nlp :: HOADNLPModel, x :: AbstractVector)
  Jx = jac(nlp, x)
  m, n = nlp.meta.ncon, nlp.meta.nvar
  if VERSION < v"1.0"
    I = [(i,j) for i = 1:m, j = 1:n]
  else
    I = ((i,j) for i = 1:m, j = 1:n)
  end
  return (getindex.(I, 1)[:], getindex.(I, 2)[:], Jx[:])
end

function jprod(nlp :: HOADNLPModel, x :: AbstractVector, v :: AbstractVector)
  # increment!(nlp, :neval_jprod)
  return ForwardDiff.jacobian(nlp.c, x) * v
end

function jprod!(nlp :: HOADNLPModel, x :: AbstractVector, v :: AbstractVector, Jv :: AbstractVector)
  # increment!(nlp, :neval_jprod)
  Jv[1:nlp.meta.ncon] = ForwardDiff.jacobian(nlp.c, x) * v
  return Jv
end

function jtprod(nlp :: HOADNLPModel, x :: AbstractVector, v :: AbstractVector)
  # increment!(nlp, :neval_jtprod)
  return ForwardDiff.jacobian(nlp.c, x)' * v
end

function jtprod!(nlp :: HOADNLPModel, x :: AbstractVector, v :: AbstractVector, Jtv :: AbstractVector)
  # increment!(nlp, :neval_jtprod)
  Jtv[1:nlp.meta.nvar] = ForwardDiff.jacobian(nlp.c, x)' * v
  return Jtv
end

function hess(nlp :: HOADNLPModel, x :: AbstractVector; obj_weight :: Real = one(eltype(x)), y :: AbstractVector = eltype(x)[])
  increment!(nlp, :neval_hess)
  Hx = obj_weight == 0.0 ? spzeros(nlp.meta.nvar, nlp.meta.nvar) :
       ForwardDiff.hessian(nlp.f, x) * obj_weight
  for i = 1:min(length(y), nlp.meta.ncon)
    if y[i] != 0.0
      Hx += ForwardDiff.hessian(x->nlp.c(x)[i], x) * y[i]
    end
  end
  return tril(Hx)
end

function hess_coord(nlp :: HOADNLPModel, x :: AbstractVector; obj_weight :: Real = one(eltype(x)), y :: AbstractVector = eltype(x)[])
  Hx = hess(nlp, x, obj_weight=obj_weight, y=y)
  n = nlp.meta.nvar
  if VERSION < v"1.0"
    I = [(i,j,Hx[i,j]) for i = 1:n, j = 1:n if i ≥ j]
  else
    I = ((i,j,Hx[i,j]) for i = 1:n, j = 1:n if i ≥ j)
  end
  return (getindex.(I, 1), getindex.(I, 2), getindex.(I, 3))
end

function hprod(nlp :: HOADNLPModel, x :: AbstractVector, v :: AbstractVector;
               obj_weight :: Real = one(eltype(x)), y :: AbstractVector = eltype(x)[])
  Hv = zeros(nlp.meta.nvar)
  return hprod!(nlp, x, v, Hv, obj_weight=obj_weight, y=y)
end

function hprod!(nlp :: HOADNLPModel, x :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector;
                obj_weight :: Real = one(eltype(x)), y :: AbstractVector = eltype(x)[])
  increment!(nlp, :neval_hprod)
  n = nlp.meta.nvar
  Hv[1:n] = obj_weight == 0.0 ? zeros(nlp.meta.nvar) :
          ForwardDiff.hessian(nlp.f, x) * v * obj_weight
  for i = 1:min(length(y), nlp.meta.ncon)
    if y[i] != 0.0
      Hv[1:n] += ForwardDiff.hessian(x->nlp.c(x)[i], x) * v * y[i]
    end
  end
  return Hv
end

function ∇f³xuv(nlp :: HOADNLPModel, x :: AbstractVector, u :: AbstractVector, v :: AbstractVector)
    ∇f(x) = ForwardDiff.gradient(nlp.f, x)
    # println("On  a ∇f")
    ϕᵤ(x) = ∇f(x)'*u
    # println("On  a ϕᵤ")
    ∇ϕᵤ(x) = ForwardDiff.gradient(ϕᵤ, x)
    # println("On  a ∇ϕᵤ")
    ψᵥ(x) =  ∇ϕᵤ(x)' * v
    # println("On  a ψᵥ")
    ∇ψᵥ(x) = ForwardDiff.gradient(ψᵥ, x)
    # println("On  a ∇ψᵥ")
    # increment!(nlp, :neval_Tuvprod)
    return ∇ψᵥ(x)
end

function ∇f³xuv2(nlp :: HOADNLPModel, x :: AbstractVector, u :: AbstractVector, v :: AbstractVector)
    ∇²fu(x) = ForwardDiff.hessian(nlp.f, x) * u
    ∇²fuv(x) = ∇²fu(x)' * v
    return ForwardDiff.gradient(∇²fuv, x)
end


text:
# For why we use Base.unsafe_trunc instead of floor below see:
# https://github.com/CliMA/Oceananigans.jl/issues/828
# https://github.com/CliMA/Oceananigans.jl/pull/997

"""
    permute_index(i, N)

Permute `i` such that, for example, `i ∈ 1:N` becomes

    [1, 2, 3, 4, 5, 6, 7, 8] -> [1, 8, 2, 7, 3, 6, 4, 5]

    [1, 2, 3, 4, 5, 6, 7, 8, 9] -> [1, 9, 2, 8, 3, 7, 4, 6, 5]

for `N=8` and `N=9` respectively.

See equation (20) of [Makhoul80](@cite).
"""
@inline permute_index(i, N)::Int = ifelse(isodd(i),
                                          Base.unsafe_trunc(Int, i/2) + 1,
                                          N - Base.unsafe_trunc(Int, (i-1)/2))

"""
    unpermute_index(i, N)

Permute `i` in the opposite manner as `_permute_index`, such that,
for example, `i ∈ 1:N` becomes

   [1, 2, 3, 4, 5, 6, 7, 8] -> [1, 3, 5, 7, 8, 6, 4, 2]

   [1, 2, 3, 4, 5, 6, 7, 8, 9] -> [1, 3, 5, 7, 9, 8, 6, 4, 2]

for `N=8` and `N=9` respectively.

See equation (20) of [Makhoul80](@cite).
"""
@inline unpermute_index(i, N) = ifelse(i <= ceil(N/2), 2i-1, 2(N-i+1))

@kernel function permute_x_indices!(dst, src, grid)
    i, j, k = @index(Global, NTuple)
    i′ = permute_index(i, grid.Nx)
    @inbounds dst[i′, j, k] = src[i, j, k]
end

@kernel function permute_y_indices!(dst, src, grid)
    i, j, k = @index(Global, NTuple)
    j′ = permute_index(j, grid.Ny)
    @inbounds dst[i, j′, k] = src[i, j, k]
end

@kernel function permute_z_indices!(dst, src, grid)
    i, j, k = @index(Global, NTuple)
    k′ = permute_index(k, grid.Nz)
    @inbounds dst[i, j, k′] = src[i, j, k]
end

@kernel function unpermute_x_indices!(dst, src, grid)
    i, j, k = @index(Global, NTuple)
    i′ = unpermute_index(i, grid.Nx)
    @inbounds dst[i′, j, k] = src[i, j, k]
end

@kernel function unpermute_y_indices!(dst, src, grid)
    i, j, k = @index(Global, NTuple)
    j′ = unpermute_index(j, grid.Ny)
    @inbounds dst[i, j′, k] = src[i, j, k]
end

@kernel function unpermute_z_indices!(dst, src, grid)
    i, j, k = @index(Global, NTuple)
    k′ = unpermute_index(k, grid.Nz)
    @inbounds dst[i, j, k′] = src[i, j, k]
end

permute_kernel! = Dict(
    1 => permute_x_indices!,
    2 => permute_y_indices!,
    3 => permute_z_indices!
)

unpermute_kernel! = Dict(
    1 => unpermute_x_indices!,
    2 => unpermute_y_indices!,
    3 => unpermute_z_indices!
)

function permute_indices!(dst, src, arch, grid, dim)
    event = launch!(arch, grid, :xyz, permute_kernel![dim], dst, src, grid, dependencies=Event(device(arch)))
    wait(device(arch), event)
end

function unpermute_indices!(dst, src, arch, grid, dim)
    event = launch!(arch, grid, :xyz, unpermute_kernel![dim], dst, src, grid, dependencies=Event(device(arch)))
    wait(device(arch), event)
end