text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack
using StaticArrays
using DelimitedFiles
using Polyester
using MuladdMacro

push!(LOAD_PATH, "./src")

using CommonUtils
using Basis2DTri
using UniformTriMesh
using NodesAndModes

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
include("../../EntropyStableEuler.jl/src/logmean.jl")
include("../SBP_quad_data.jl")
using EntropyStableEuler
using EntropyStableEuler.Fluxes2D

function init_reference_tri_sbp_GQ(N, qnode_choice)
    include("SBP_quad_data.jl")
    # initialize a new reference element data struct
    rd = RefElemData()

    fv = tri_face_vertices() # set faces for triangle
    Nfaces = length(fv)
    @pack! rd = fv, Nfaces

    # Construct matrices on reference elements
    r, s = nodes(Tri(),N)
    VDM = vandermonde(Tri(),N, r, s)
    Vr, Vs = grad_vandermonde(Tri(),N, r, s)
    Dr = Vr/VDM
    Ds = Vs/VDM
    @pack! rd = r,s,VDM,Dr,Ds

    # low order interpolation nodes
    r1,s1 = nodes(Tri(),1)
    V1 = vandermonde(Tri(),1,r,s)/vandermonde(Tri(),1,r1,s1)
    @pack! rd = V1

    #Nodes on faces, and face node coordinate
    if qnode_choice == "GQ"
        r1D, w1D = gauss_quad(0,0,N)
    elseif qnode_choice == "GL" || qnode_choice == "tri_diage"
        r1D, w1D = gauss_lobatto_quad(0,0,N+1)
    end
    Nfp = length(r1D) # number of points per face
    e = ones(Nfp) # vector of all ones
    z = zeros(Nfp) # vector of all zeros
    rf = [r1D; -r1D; -e];
    sf = [-e; r1D; -r1D];
    wf = vec(repeat(w1D,3,1));
    nrJ = [z; e; -e]
    nsJ = [-e; e; z]
    @pack! rd = rf,sf,wf,nrJ,nsJ

    if qnode_choice == "GQ"
        rq,sq,wq = GQ_SBP[N];
    elseif qnode_choice == "GL"
        rq,sq,wq = GL_SBP[N];
    elseif qnode_choice == "tri_diage"
        rq,sq,wq = Tri_diage[N];
    end
    # rq,sq,wq = GQ_SBP[N]
    # rq,sq,wq = GL_SBP[N]
    # rq,sq,wq = Tri_diage[N]
    Vq = vandermonde(Tri(),N,rq,sq)/VDM
    M = Vq'*diagm(wq)*Vq
    Pq = M\(Vq'*diagm(wq))
    @pack! rd = rq,sq,wq,Vq,M,Pq

    Vf = vandermonde(Tri(),N,rf,sf)/VDM # interpolates from nodes to face nodes
    LIFT = M\(Vf'*diagm(wf)) # lift matrix used in rhs evaluation
    @pack! rd = Vf,LIFT

    # plotting nodes
    rp, sp = equi_nodes(Tri(),10)
    Vp = vandermonde(Tri(),N,rp,sp)/VDM
    @pack! rd = rp,sp,Vp

    return rd
end

function build_meshfree_sbp(rq,sq,wq,rf,sf,wf,nrJ,nsJ,α)
    # [-1,1,0], [-1,-1,sqrt(4/3)]
    equilateral_map(r,s) = (@. .5*(2*r+1*s+1), @. sqrt(3)*(1+s)/2 - 1/sqrt(3) )
    req,seq = equilateral_map(rq,sq)
    ref,sef = equilateral_map(rf,sf)
    barycentric_coords(r,s) = ((@. (1+r)/2), (@. (1+s)/2), (@. -(r+s)/2))
    λ1,λ2,λ3 = barycentric_coords(rq,sq)
    λ1f,λ2f,λ3f = barycentric_coords(rf,sf)

    Br = diagm(nrJ.*wf)
    Bs = diagm(nsJ.*wf)

    # build extrapolation matrix
    E = zeros(length(rf),length(rq))
    for i = 1:length(rf)
        # d = @. (λ1 - λ1f[i])^2 + (λ2 - λ2f[i])^2 + (λ3 - λ3f[i])^2
        d2 = @. (req-ref[i])^2 + (seq-sef[i])^2
        p = sortperm(d2)
        h2 = (wf[i]/sum(wf))*2/pi # set so that h = radius of circle with area w_i = face weight
        nnbrs = min(4,max(3,count(d2[p] .< h2))) # find 3 closest points
        p = p[1:nnbrs]
        Ei = vandermonde_2D(1,[rf[i]],[sf[i]])/vandermonde_2D(1,rq[p],sq[p])
        E[i,p] = Ei
    end
    E = Matrix(droptol!(sparse(E),1e-13))

    # build stencil
    A = spzeros(length(req),length(req))
    for i = 1:length(req)
        d2 = @. (req-req[i])^2 + (seq-seq[i])^2
        p = sortperm(d2)

        # h^2 = wq[i]/pi = radius of circle with area wq[i]
        # h2 =     (sqrt(3)/sum(wq))*wq[i]/pi
        h2 = α^2*(sqrt(3)/sum(wq))*wq[i]/pi

        nnbrs = count(d2[p] .< h2)
        nbrs = p[1:nnbrs]
        A[i,nbrs] .= one(eltype(A))
    end
    A = (A+A')
    A.nzval .= one(eltype(A)) # bool-ish

    # build graph Laplacian
    L1 = (A-diagm(diag(A))) # ignore
    L1 -= diagm(vec(sum(L1,dims=2)))

    b1r = -sum(.5*E'*Br*E,dims=2)
    b1s = -sum(.5*E'*Bs*E,dims=2)
    ψ1r = pinv(L1)*b1r
    ψ1s = pinv(L1)*b1s

    function fillQ(adj,ψ)
        Np = length(ψ)
        S = zeros(Np,Np)
        for i = 1:Np
            for j = 1:Np
                if adj[i,j] != 0
                        S[i,j] += (ψ[j]-ψ[i])
                end
            end
        end
        return S
    end

    S1r,S1s = fillQ.((A,A),(ψ1r,ψ1s))
    Qr = Matrix(droptol!(sparse(S1r + .5*E'*Br*E),1e-14))
    Qs = Matrix(droptol!(sparse(S1s + .5*E'*Bs*E),1e-14))

    return Qr,Qs,E,Br,Bs,A
end

function construct_SBP(N)
    # Construct matrices on reference elements
    rd = init_reference_tri_sbp_GQ(N, "GQ")
    @unpack r,s,rf,sf,wf,rq,sq,wq,nrJ,nsJ = rd
    @unpack VDM,V1,Vq,Vf,Dr,Ds,M,Pq,LIFT = rd

    V = vandermonde_2D(N, r, s)
    Vr, Vs = grad_vandermonde_2D(N, r, s)
    M = inv(V*V')
    Dr_ES = Vr/V; Ds = Vs/V;
    M = Vq'*diagm(wq)*Vq

    Qr_ES = M*Dr;
    Qs_ES = M*Ds;
    Pq = M\(Vq'*diagm(wq));

    # diff. matrices redefined in terms of quadrature points
    Qr_ES = Pq'*Qr_ES*Pq;
    Qs_ES = Pq'*Qs_ES*Pq;
    E_ES = Vf*Pq;

    return Qr_ES,Qs_ES,E_ES
end

const TOL = 5e-16
const POSTOL = 1e-14
const WALLPT = 1.0/6.0
const Nc = 4 # number of components
"Approximation parameters"
const N = 2
const K1D = 10
const T = 1.0
const dt0 = 1e-2
const XLENGTH = 2.0
const CFL = 0.5
const NUM_THREADS = Threads.nthreads()
qnode_choice = "GQ" #"GQ" "GL" "tri_diage"

# Initial condition 2D shocktube
const γ = 1.4

# Mesh related variables
Kx = Int(round(XLENGTH*K1D))
Ky = K1D
VX, VY, EToV = uniform_tri_mesh(Kx,Ky)

# @. VX = 15*(1+VX)/2
# @. VY = 5*VY
@. VX = 10*(1+VX)/2
@. VY = 5*(1+VY)/2
FToF = CommonUtils.connect_mesh(EToV,tri_face_vertices())
Nfaces,K = size(FToF)

# Construct matrices on reference elements
rd = init_reference_tri_sbp_GQ(N, qnode_choice)
@unpack r,s,rf,sf,wf,rq,sq,wq,nrJ,nsJ = rd
@unpack VDM,V1,Vq,Vf,Dr,Ds,M,Pq,LIFT = rd

V = vandermonde_2D(N, r, s)
Vr, Vs = grad_vandermonde_2D(N, r, s)
M = inv(V*V')
Dr_ES = Vr/V; Ds = Vs/V;
M = Vq'*diagm(wq)*Vq

Qr_ES = M*Dr;
Qs_ES = M*Ds;
Pq = M\(Vq'*diagm(wq));

# diff. matrices redefined in terms of quadrature points
Qr_ES = Pq'*Qr_ES*Pq;
Qs_ES = Pq'*Qs_ES*Pq;
E_ES = Vf*Pq;

# Need to choose α so that Qr, Qs have zero row sums (and maybe a minimum number of neighbors)
# α = 4 # for N=1
# α = 2.5 #for N=2
α = 3.5 # for N=3
Qr_ID,Qs_ID,E,Br,Bs,A = build_meshfree_sbp(rq,sq,wq,rf,sf,wf,nrJ,nsJ,α)
if (norm(sum(Qr_ID,dims=2)) > 1e-10) | (norm(sum(Qs_ID,dims=2)) > 1e-10)
    error("Qr_ID or Qs_ID doesn't sum to zero for α = $α")
end
Qr_ID = Matrix(droptol!(sparse(Qr_ID),1e-15))
Qs_ID = Matrix(droptol!(sparse(Qs_ID),1e-15))
Qrskew_ID = .5*(Qr_ID-transpose(Qr_ID))
Qsskew_ID = .5*(Qs_ID-transpose(Qs_ID))

# "Construct global coordinates"
x = V1*VX[transpose(EToV)]
y = V1*VY[transpose(EToV)]

# "Connectivity maps"
xf,yf = (x->Vf*x).((x,y))
mapM,mapP,mapB = build_node_maps((xf,yf),FToF)
mapM = reshape(mapM,length(wf),K)
global mapP = reshape(mapP,length(wf),K)

mapM = reshape(collect(1:K*length(wf)),length(wf),K)
mapP = reshape(collect(1:K*length(wf)),length(wf),K)

for k = 1:K
    Npf = div(length(wf),Nfaces)
    k1 = Int(k+2*(mod(k,2)-.5))
    # k2 = Int(mod1(k+2*3*(mod(k-1,2)-.5),2*Kx)+2*Kx*div(k-1,2*Kx))
    k2 = Int(mod1(k+2*3*(mod(k,2)-.5),2*Kx)+2*Kx*div(k-1,2*Kx))
    k3 = Int(mod1(k+2*(2*Kx-1)*(mod(k-1,2)-.5),K))
    for i = 1:Npf
        mapP[i,k]       = mapM[Npf-i+1,k1]
        mapP[i+Npf,k]   = mapM[2*Npf-i+1,k2]
        mapP[i+2*Npf,k] = mapM[3*Npf-i+1,k3]
    end
end

# "Geometric factors and surface normals"
rxJ, sxJ, ryJ, syJ, J = CommonUtils.geometric_factors(x, y, Dr, Ds)
rxJ = Matrix(droptol!(sparse(rxJ),1e-14)); sxJ = Matrix(droptol!(sparse(sxJ),1e-14))
ryJ = Matrix(droptol!(sparse(ryJ),1e-14)); syJ = Matrix(droptol!(sparse(syJ),1e-14))
nxJ = (Vf*rxJ).*nrJ + (Vf*sxJ).*nsJ;
nyJ = (Vf*ryJ).*nrJ + (Vf*syJ).*nsJ;
sJ = @. sqrt(nxJ^2 + nyJ^2)
nx = nxJ./sJ; ny = nyJ./sJ;
rxJ,sxJ,ryJ,syJ,J = (x->Vq*x).((rxJ,sxJ,ryJ,syJ,J))

# TODO:
J = -J
rxJ = -rxJ
ryJ = -ryJ
sxJ = -sxJ
syJ = -syJ

QNr = [Qr_ES - .5*E_ES'*Br*E_ES .5*E_ES'*Br;
    -.5*Br*E_ES .5*Br];
QNs = [Qs_ES - .5*E_ES'*Bs*E_ES .5*E_ES'*Bs;
    -.5*Bs*E_ES .5*Bs];

VN_sbp = [eye(length(wq));E];
QNr_sbp = VN_sbp'*QNr*VN_sbp;
@show norm(QNr_sbp+QNr_sbp' - E'*diagm(wf.*nrJ)*E)
QNs_sbp = VN_sbp'*QNs*VN_sbp;
@show norm(QNs_sbp+QNs_sbp' - E'*diagm(wf.*nsJ)*E)

Qrskew_ES = .5*(QNr_sbp-QNr_sbp');
Qsskew_ES = .5*(QNs_sbp-QNs_sbp');

Sr = Qrskew_ES
Ss = Qsskew_ES
S0r = Qrskew_ID
S0s = Qsskew_ID
M   = diagm(wq)
Minv   = diagm(1 ./ wq)

Fmask = zeros(Int64,size(Vf,1))
for i = 1:size(Vf,1)
    tmparr = findall(abs.(E[i,:] .- 1) .< 1e-12)
    Fmask[i] = tmparr[1]
end

xq = Vq*x
yq = Vq*y

gr(aspect_ratio=:equal,legend=false,
   markerstrokewidth=0,markersize=2)
VDM = vandermonde(Tri(),N, rq, sq)
rp,sp = equi_nodes_2D(10)
Vp = vandermonde_2D(N,rp,sp)/VDM

@inline function pfun(rho,rhou,E)
    return (γ-1)*(E-.5*rhou^2/rho)
end

@inline function pfun(rho,rhou,rhov,E)
    return (γ-1)*(E-.5*(rhou^2+rhov^2)/rho)
end

@inline function Efun(rho,u,v,p)
    return p/(γ-1) + .5*rho*(u^2+v^2)
end

@inline function wavespeed_1D(rho,rhou,E)
    p = pfun(rho,rhou,E)
    return abs(rhou/rho) + sqrt(γ*p/rho)
end

@inline function logmean(aL,aR,logL,logR)

    # "from: Entropy stable num. approx. for the isothermal and polytropic Euler"

    da = aR-aL;
    aavg = .5*(aR+aL);
    f = da/aavg;
    v = f^2;
    if abs(f)<1e-4
        # numbers assume the specific value γ = 1.4
        return aavg*(1 + v*(-.2-v*(.0512 - v*0.026038857142857)))
    else
        return -da/(logL-logR)
    end
end

@inline function euler_fluxes_2D(rhoL,uL,vL,betaL,rhologL,betalogL,
                                 rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog  = logmean(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = .5*(rhoL+rhoR)
    uavg   = .5*(uL+uR)
    vavg   = .5*(vL+vR)

    unorm = uL*uR + vL*vR
    pa    = rhoavg/(betaL+betaR)
    f4aux = rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm

    FxS1 = rholog*uavg
    FxS2 = FxS1*uavg + pa
    FxS3 = FxS1*vavg
    FxS4 = f4aux*uavg

    FyS1 = rholog*vavg
    FyS2 = FxS3
    FyS3 = FyS1*vavg + pa
    FyS4 = f4aux*vavg

    return FxS1,FxS2,FxS3,FxS4,FyS1,FyS2,FyS3,FyS4
end

@inline function inviscid_flux_prim(rho,u,v,p)
    E = Efun(rho,u,v,p)

    rhou  = rho*u
    rhov  = rho*v
    rhouv = rho*u*v
    Ep    = E+p

    fx1 = rhou
    fx2 = rhou*u+p
    fx3 = rhouv
    fx4 = u*Ep

    fy1 = rhov
    fy2 = rhouv
    fy3 = rhov*v+p
    fy4 = v*Ep

    return fx1,fx2,fx3,fx4,fy1,fy2,fy3,fy4
end

@inline function get_infoP(mapP,Fmask,i,k)
    gP = mapP[i,k]           # exterior global face node number
    kP = fld1(gP,Nfp)        # exterior element number
    iP = Fmask[mod1(gP,Nfp)] # exterior node number
    return iP,kP
end

function rhs_ESDG!(U,rhsU,ops,geom)
    Sr,Ss,S0r,S0s,Br,Bs,Minv,E = ops
    mapP,Fmask,rxJ,ryJ,sxJ,syJ,J = geom

    FHx = zeros(Nc,Np,Np)
    FHy = zeros(Nc,Np,Np)
    FHPx = zeros(Nc,Nfp)
    FHPy = zeros(Nc,Nfp)
    Sx  = zeros(Np,Np)
    Sy  = zeros(Np,Np)
    EBx  = zeros(Np,Np)
    EBy  = zeros(Np,Np)
    rhsvec = zeros(Nc,Np)

    Br = diagm(wf.*nrJ)
    Bs = diagm(wf.*nsJ)
    for k = 1:K
        for j = 1:Np
            for i = 1:Np
                rhoM     = U[1,i,k]
                rhouM    = U[2,i,k]
                rhovM    = U[3,i,k]
                EM       = U[4,i,k]
                uM       = rhouM/rhoM
                vM       = rhovM/rhoM
                pM       = pfun(rhoM,rhouM,rhovM,EM)
                betaM    = rhoM/(2*pM)
                rhologM  = log(rhoM)
                betalogM = log(betaM)
                rhoP     = U[1,j,k]
                rhouP    = U[2,j,k]
                rhovP    = U[3,j,k]
                EP       = U[4,j,k]
                uP       = rhouP/rhoP
                vP       = rhovP/rhoP
                pP       = pfun(rhoP,rhouP,rhovP,EP)
                betaP    = rhoP/(2*pP)
                rhologP  = log(rhoP)
                betalogP = log(betaP)

                Fx1,Fx2,Fx3,Fx4,Fy1,Fy2,Fy3,Fy4 = euler_fluxes_2D(rhoM,uM,vM,betaM,rhologM,betalogM,
                                                                  rhoP,uP,vP,betaP,rhologP,betalogP)

                FHx[1,i,j] = Fx1
                FHx[2,i,j] = Fx2
                FHx[3,i,j] = Fx3
                FHx[4,i,j] = Fx4
                FHy[1,i,j] = Fy1
                FHy[2,i,j] = Fy2
                FHy[3,i,j] = Fy3
                FHy[4,i,j] = Fy4
            end
        end

        for i = 1:Nfp
            iM       = Fmask[i]
            rhoM     = U[1,iM,k]
            rhouM    = U[2,iM,k]
            rhovM    = U[3,iM,k]
            EM       = U[4,iM,k]
            uM       = rhouM/rhoM
            vM       = rhovM/rhoM
            pM       = pfun(rhoM,rhouM,rhovM,EM)
            betaM    = rhoM/(2*pM)
            rhologM  = log(rhoM)
            betalogM = log(betaM)

            iP,kP = get_infoP(mapP,Fmask,i,k)
            rhoP     = U[1,iP,kP]
            rhouP    = U[2,iP,kP]
            rhovP    = U[3,iP,kP]
            EP       = U[4,iP,kP]
            uP       = rhouP/rhoP
            vP       = rhovP/rhoP
            pP       = pfun(rhoP,rhouP,rhovP,EP)
            betaP    = rhoP/(2*pP)
            rhologP  = log(rhoP)
            betalogP = log(betaP)

            Fx1,Fx2,Fx3,Fx4,Fy1,Fy2,Fy3,Fy4 = euler_fluxes_2D(rhoM,uM,vM,betaM,rhologM,betalogM,
                                                              rhoP,uP,vP,betaP,rhologP,betalogP)

            FHPx[1,i] = Fx1
            FHPx[2,i] = Fx2
            FHPx[3,i] = Fx3
            FHPx[4,i] = Fx4
            FHPy[1,i] = Fy1
            FHPy[2,i] = Fy2
            FHPy[3,i] = Fy3
            FHPy[4,i] = Fy4
        end

        Sx  =     rxJ[1,k]*Sr + sxJ[1,k]*Ss
        Sy  =     ryJ[1,k]*Sr + syJ[1,k]*Ss
        EBx = E'*(rxJ[1,k]*Br + sxJ[1,k]*Bs)
        EBy = E'*(ryJ[1,k]*Br + syJ[1,k]*Bs)

        for c = 1:Nc
            rhsvec[c,:] = 2*(Sx .* FHx[c,:,:])*ones(Np) + 2*(Sy .* FHy[c,:,:])*ones(Np) + EBx*FHPx[c,:] + EBy*FHPy[c,:]
            rhsvec[c,:] = -(1/J[1,k]*Minv)*rhsvec[c,:]
        end

        for c = 1:Nc
            rhsU[c,:,k] = rhsvec[c,:]
        end
    end
end

function vortex_sol(x,y,t)
    # x0 = 5
    # y0 = 0
    x0 = 4.5
    y0 = 2.5
    beta = 5
    r2 = @. (x-x0-t)^2 + (y-y0)^2

    u = @. 1 - beta*exp(1-r2)*(y-y0)/(2*pi)
    v = @. beta*exp(1-r2)*(x-x0-t)/(2*pi)
    rho = @. 1 - (1/(8*γ*pi^2))*(γ-1)/2*(beta*exp(1-r2))^2
    rho = @. rho^(1/(γ-1))
    p = @. rho^γ

    return (rho, u, v, p)
end

# const Np  = size(x,1)
const Np  = size(xq,1)
const Nfp = size(Vf,1)
U = zeros(Nc,Np,K)
for k = 1:K
    for i = 1:Np
        # rho,u,v,p = vortex_sol(x[i,k],y[i,k],0.0)
        rho,u,v,p = vortex_sol(xq[i,k],yq[i,k],0.0)
        U[1,i,k] = rho
        U[2,i,k] = rho*u
        U[3,i,k] = rho*v
        U[4,i,k] = Efun(rho,u,v,p)
    end
end

ops      = (Sr,Ss,S0r,S0s,Br,Bs,Minv,E)
geom     = (mapP,Fmask,rxJ,ryJ,sxJ,syJ,J)

# TODO: previous Time stepping
"Time integration"
t = 0.0
U = collect(U)
rhsU = zeros(size(U))
resW = zeros(size(U))

dt_hist = []
i = 1

@time while t < T
    @show t
    # SSPRK(3,3)
    # dt = min(dt0,T-t)
    dt = dt0
    rhs_ESDG!(U,rhsU,ops,geom)
    dt = min(CFL*dt,T-t)
    @. resW = U + dt*rhsU
    rhs_ESDG!(resW,rhsU,ops,geom)
    @. resW = resW+dt*rhsU
    @. resW = 3/4*U+1/4*resW
    rhs_ESDG!(resW,rhsU,ops,geom)
    @. resW = resW+dt*rhsU
    @. U = 1/3*U+2/3*resW

    push!(dt_hist,dt)
    global t = t + dt
    println("Current time $t with time step size $dt, and final time $T, at step $i")
    flush(stdout)
    global i = i + 1
end

exact_U = @. vortex_sol.(x,y,T)
exact_rho = [x[1] for x in exact_U]
exact_u = [x[2] for x in exact_U]
exact_v = [x[3] for x in exact_U]
exact_rhou = exact_rho .* exact_u
exact_rhov = exact_rho .* exact_v
exact_p = [x[4] for x in exact_U]
exact_E = Efun.(exact_rho,exact_u,exact_v,exact_p)

rho = U[1,:,:]
u = U[2,:,:]./U[1,:,:]
v = U[3,:,:]./U[1,:,:]
rhou = U[2,:,:]
rhov = U[3,:,:]
# E = U[4,:,:]

vv = Vp*U[1,:,:]
scatter(Vp*xq,Vp*yq,vv,zcolor=vv,camera=(0,90))
savefig("~/Desktop/tmptri.png")

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
using CompositeStructs, Test

@testset "CompositeStructs" begin

    # tests work by explicilty redefing the @composite structs, which
    # will be an error if its not exactly the same

    # docstring example
    @test_nowarn @eval module $(gensym())
        using CompositeStructs
        struct Foo{X,Y}
            x :: X
            y :: Y
        end
        @composite struct Bar{X,Y,Z}
            Foo{X,Y}...
            z :: Z
        end
        struct Bar{X,Y,Z}
            x :: X
            y :: Y
            z :: Z
        end
    end

    # composite type in Foo, swap Y/X
    @test_nowarn @eval module $(gensym())
        using CompositeStructs
        struct Foo{X<:Real,Y}
            x :: X
            y :: Complex{Y}
        end
        @composite struct Bar{X,Y,Z}
            Foo{Y,X}...
            z :: Z
        end
        struct Bar{X,Y,Z}
            x :: Y
            y :: Complex{X}
            z :: Z
        end
    end

    # multiple splices / order
    @test_nowarn @eval module $(gensym())
        using CompositeStructs
        struct X x end
        struct Z z end
        @composite struct Bar
            X...
            y
            Z...
        end
        struct Bar
            x
            y
            z
        end
    end

    # non-concrete field
    @test_nowarn @eval module $(gensym())
        using CompositeStructs
        struct Foo{T}
            x :: Array{T,N} where {N}
        end
        @composite struct Bar{T}
            Foo{T}...
            y
        end
        struct Bar{T}
            x :: Array{T,N} where {N}
            y
        end
    end

    # mutable structs, with supertype
    @test_nowarn @eval module $(gensym())
        using CompositeStructs
        mutable struct Foo <: Real
            x
        end
        @composite mutable struct Bar <: Number
            Foo...
            y
        end
        mutable struct Bar <: Number
            x
            y
        end
    end

    # kwdef
    @test_nowarn @eval module $(gensym())
        using CompositeStructs, Test

        Base.@kwdef struct Child1{T} <: Real
            x::T = 1
        end

        Base.@kwdef struct Child2
            y
        end

        @composite Base.@kwdef struct ParentInner{T<:Any} <: Number
            Child1{T}...
            Child2...
            z = 3
            w
        end

        @test ParentInner(y=2, w=4) == ParentInner{Int64}(1,2,3,4)
        @test ParentInner{Float64}(y=2, w=4) == ParentInner{Float64}(1.0,2,3,4)
        @test ParentInner(x="x", y="y", z="z", w="w") == ParentInner{String}("x", "y", "z", "w")

        Base.@kwdef @composite struct ParentOuter{T<:Any} <: Number
            Child1{T}...
            Child2...
            z = 3
            w
        end

        @test ParentOuter(x=1, y=2, w=4) == ParentOuter{Int64}(1,2,3,4)

    end

end

ext: jl


text:
using ODEInterfaceDiffEq, DiffEqProblemLibrary, DiffEqBase, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear,
                            prob_ode_2Dlinear, prob_ode_vanderpol

prob = prob_ode_linear
sol =solve(prob,dopri5(),dt=1//2^(4))

sol =solve(prob,dop853();dt=1//2^(4))

sol =solve(prob,odex();dt=1//2^(4))

sol =solve(prob,seulex();dt=1//2^(4))

sol =solve(prob,radau();dt=1//2^(4))

sol =solve(prob,radau5();dt=1//2^(4))

sol =solve(prob,rodas();dt=1//2^(4))

sol =solve(prob,ddeabm();dt=1//2^(4))

sol =solve(prob,ddebdf();dt=1//2^(4))

prob = prob_ode_2Dlinear

sol =solve(prob,dopri5(),dt=1//2^4)

sol =solve(prob,dop853();dt=1//2^(4))

sol =solve(prob,odex();dt=1//2^(4))

sol =solve(prob,seulex();dt=1//2^(4))

sol =solve(prob,radau();dt=1//2^(4))

sol =solve(prob,radau5();dt=1//2^(4))

sol =solve(prob,rodas();dt=1//2^(4))

sol =solve(prob,ddeabm();dt=1//2^(4))

sol =solve(prob,ddebdf();dt=1//2^(4))

prob = prob_ode_vanderpol

sol =solve(prob,dopri5(),dt=1//2^4)

sol =solve(prob,dop853();dt=1//2^(4))

sol =solve(prob,odex();dt=1//2^(4))

sol =solve(prob,seulex();dt=1//2^(4))

sol =solve(prob,radau();dt=1//2^(4))

sol =solve(prob,radau5();dt=1//2^(4))

sol =solve(prob,rodas();dt=1//2^(4))

sol =solve(prob,ddeabm();dt=1//2^(4))

sol =solve(prob,ddebdf();dt=1//2^(4))


text:
####
# Helpers
###
# Contains statistical and other routines used by data-driven sets
using Distributions, Roots, JuMP, Optim

export boot, calcMeansT, calcSigsBoot
export boot_mu, boot_sigma, bootDY_mu, bootDY_sigma, calc_ab_thresh

###Bootstrapping code
# ideally this should all be moved to some base level function
function boot(data::Vector, fun::Function, prob::Float64, numBoots::Int, f_args...)
	const N = size(data, 1)
	dist = DiscreteUniform(1, N)
	out = zeros(Float64, numBoots)
	indices = collect(1:N)
	for i = 1:numBoots
		rand!(dist, indices)
		out[i] = fun(data[indices], f_args...)
	end
	quantile(out, prob)
end

function boot(data::Matrix, fun::Function, prob::Float64, numBoots::Int, f_args...)
	const N = size(data, 1)
	dist = DiscreteUniform(1, N)
	out = zeros(Float64, numBoots)
	indices = collect(1:N)
	for i = 1:numBoots
		rand!(dist, indices)
		out[i] = fun(data[indices, :], f_args...)
	end
	quantile(out, prob)
end

#search for a valid bracket multiplicatively
#primarily used to prep for root-finding
function bracket(fun, guess; max_iter=32, ratio=2., trace=false)
    low = guess/ratio
    hi  = guess * ratio

    f_low  = fun(low)
    f_hi = fun(hi)
    iter = 0
    trace && println("Bracketing:")
    while((f_low * f_hi >= 0) && (iter < max_iter))
        trace && println("$iter \t $low \t $hi \t $f_low \t $f_hi")
        low /= ratio
        hi  *= ratio
        f_low  = fun(low)
        f_hi = fun(hi)
        iter += 1
    end
    iter == max_iter && error("Bracketing Failed")

    if f_hi * fun(hi/ratio) <= 0
        return hi/ratio, hi
    elseif f_low * fun(low*ratio) <= 0
        return low, low*ratio
    else
        error("iter $iter \t $f_low \t $f_hi")
    end
end

#Assumes both, joint
function calcSigsBoot(data, alpha, numBoots)
    sigf = calcSigsBoot(data, alpha/2., numBoots, :Fwd)
    sigb = calcSigsBoot(data, alpha/2., numBoots, :Back)
    sigf, sigb
end

#####
# b = maximum(data) if x > 0, b= minimum(data) otherwse
function f_sig(x, mu_adj, data_b, b)
    expdata = exp(x*data_b)
    logmeanexp = log(mean(expdata))
    sqrt(-2mu_adj/x + 2/x^2 * logmeanexp)
end

#approximate deriv or rootfinding
#muadj and data_b are already adjusted for numerical stability
function df_sig(x, mu_adj, data_b)
    expdata = exp(x*data_b)
    logmeanexp = log(mean(expdata))
    wght_mean = dot(expdata, data_b) / sum(expdata)
    (x*(wght_mean + mu_adj) -2logmeanexp)/x^3
end

######

#Sgn of guess encodes whether we are looking at Fwd or Backwd
function sigFwdBack(data, guess, trace=false)
    #adjust things for stability
    const b = guess > 0 ? maximum(data) : minimum(data)
    const mu_adj = mean(data)-b
    data_b = data-b

    df_sig_(x) = df_sig(x, mu_adj, data_b)
    #repair bracket if necessary
    low, hi = bracket(df_sig_, guess, trace=trace)

    trace && println("Bracketing Success:\t $low \t $hi")

    #rootfind
    xstar = fzero(df_sig_, low, hi)
    #use root to calc sig
    f_sig(xstar, mu_adj, data_b, b)
end

#Takes a Case Specification
function calcSigsBoot(data, alpha, numBoots, CASE)
    if CASE == :Fwd
        sigFwdBack_(data_) = sigFwdBack(data_, 1.)
    elseif CASE == :Back
        sigFwdBack_(data_) = sigFwdBack(data_, -1.)
    else
        error("Case must be either :Fwd or :Back $CASE")
    end
    boot(data, sigFwdBack_, 1-alpha, numBoots)
end

#calculates the means via t approx
# if joint, bounds hold (jointly) simultaneously at level 1-alpha_
# o.w. bounds hold individually at level 1-alpha_
function calcMeansT(data, alpha_; joint=true)
    const N   = length(data)
    const sig_rt_N = std(data)/sqrt(N)
    dist      = TDist(N-1)
    alpha = joint ? alpha_/2 : alpha_
    mean(data) + quantile(dist, alpha)*sig_rt_N, mean(data) + quantile(dist, 1-alpha)*sig_rt_N
end

#Currently computed using Stephens Approximation
#Journaly of Royal Statistical Society 1970
function KSGamma(alpha, N)
       const sqrt_N = sqrt(N)
       const num = sqrt(.5 * log(2/alpha))
       const denom = sqrt_N + .12 + .11/sqrt_N
       num/denom
end

kappa(eps_) = sqrt(1./eps_ - 1.)

function boot_mu(data, alpha, numBoots)
    const muhat = mean(data, 1)
    myfun(data_b) = norm(mean(data_b, 1) - muhat)
    boot(data, myfun, 1-alpha, numBoots)
end

function boot_sigma(data, alpha, numBoots)
    const covhat = cov(data)
    myfun(data_b) = vecnorm(cov(data_b) - covhat)
    boot(data, myfun, 1-alpha, numBoots)
end

function bootDY_mu(data, alpha, numBoots)
    muhat = mean(data, 1)
    myfun(data_b) = (mu = mean(data_b, 1); ((mu-muhat) * inv(cov(data)) * (mu-muhat)')[1])
    boot(data, myfun, 1-alpha, numBoots)
end

function bootDY_sigma(data, alpha, numBoots)
    mu0  = mean(data, 1)
    sig0 = cov(data)

    function myfun(data_b)
        muhat = mean(data_b, 1)
        sighat = cov(data_b)
        mudiff = (mu0 - muhat)'*(mu0-muhat)
        f(g2) = eigmin(g2*sighat - sig0 -mudiff)

        #If the dist of sighat is degenerate, return Inf
        if eigmin(sighat) <= 0
            println("Infinite 2nd moment matrix")
            return Inf
        end

        #solve
        try
            Roots.fzero(f, [.01, 500])
        catch e
            show(e); println()
            println("Using Manual Bracketing: f(1) $(f(1))  f(2) $(f(2))")
            #do your own bracketing
            if f(1) < 0
                lb = 1; ub = 2
                iter = 0
                const MAX_ITER = 100
                while f(ub) < 0 && iter < MAX_ITER
                    println("$iter \t $(f(ub))")
                    lb = ub
                    ub = ub * 2
                    iter = iter + 1
                end
            else
                lb = .5; ub = 1
                iter = 0
                const MAX_ITER = 100
                while f(lb) > 0 && iter < MAX_ITER
                    ub = lb
                    lb = lb * .5
                    iter = iter + 1
                end
            end
            try
                fzero(f, [lb, ub])
            catch e2
                println("Manual bracketing failed")
                println("Debug Sequence")
                for i = 10.0.^[-2:5]
                    println(f(i))
                end
            end
        end

   end
   boot(data, myfun, 1-alpha, numBoots)
end

### Used by UM and UIOracle
function sort_data_cols(data)
    data_sort = zeros(eltype(data), size(data))
    const d = size(data, 2)
    for i = 1:d
        data_sort[:, i] = sort(data[:, i])
    end
    data_sort
end

########
# LCX Stuff
####

########################
# Bootstrapping computation

#Returns indx of last instance of val, 0 if not found
#Assumes sorted_list, and sort_list[start] >= val
function findlast_sort(val::Float64, sort_list::Vector; TOL::Float64=1e-10, start::Int64=1)
    i::Int64 = 1
    for i = start:length(sort_list)
        if abs(val - sort_list[i]) > TOL
            break
        end
    end
    i-1
end

## Makes a single pass through zetas/zetahats to solve
# 1/N * max_b { sum_i [zeta_i - b] ^+ - sum_i[zetahat_i - b]^+ }
function singlepass!(zetas::Vector, zetahats::Vector)
    sort!(zetas)
    sort!(zetahats)

    vstar::Float64 = mean(zetas) - zetas[1]
    vb::Float64    = mean(zetahats) - zetas[1]

    Gamma::Float64 = vstar - vb
    const N::Int64 = length(zetas)
    pbar::Float64    = 1.0
    hat_indx::Int64  = 1
    hat_indx_::Int64 = 0

    for k = 2:length(zetas)
        vstar += (zetas[k-1] - zetas[k]) * (N-k+1)/N
        hat_indx = findlast_sort(zetas[k-1], zetahats, start=hat_indx_ + 1)
        pbar -=  (hat_indx - hat_indx_)/N
        hat_indx_ = hat_indx
        vb  += (zetas[k-1] - zetas[k]) * pbar
        Gamma = max(Gamma, vstar - vb)
    end
    Gamma
end

#a::Vector, sgns::Vector is used as storage between bootstraps
function f2(boot_sample::Matrix, data::Matrix, numSamples::Int, a::Vector, sgns::Vector)
    Gamma::Float64 = 0.
    for i = 1:numSamples
        randL1!(a, sgns)
        Gamma_ = singlepass!(data*a, boot_sample*a)
        Gamma = max(Gamma, Gamma_)
    end
    Gamma
end

function randL1!(a, sgns)
    a = rand!(a)
    a = a/ sum(a)
    rand!(Bernoulli(), sgns)
    a = a .* (2*sgns-1)
end

#Approximates the threshold by sampling a bunch of abs for each bootstrap rep
function calc_ab_thresh(data::Matrix, alpha::Float64, numBoots::Int, numSamples::Int)
    const N = size(data, 1)
    const d = size(data, 2)
    a::Vector{Float64} = zeros(Float64, d)
    sgns::Vector{Int}  = zeros(Int64, d)
    boot(data, f2, 1-alpha, numBoots, data, numSamples, a, sgns)
end

function compute_cs(boot_indx)
    cs = ones(length(boot_indx))
    for i = boot_indx
        cs[i] = cs[i] - 1
    end
    return cs/length(boot_indx)
end


text:
function parse_cli_args()

  settings = ArgParseSettings(prog="sialia",
    description="Command line tools for the Sialia framework",
    add_version=true)

  @add_arg_table! settings begin
    "generate", "gen"
      help = "Create a new Sialia project, service, or infrastructure definition"
      action = :command
    "build"
      help = "Build the project"
      action = :command
    "run"
      help = "Run the project"
      action = :command
    "deploy"
      help = "Deploy the project"
      action = :command
  end

  @add_arg_table! settings["generate"] begin
    "project"
      help = "Create a new Sialia project"
      action = :command
    "service"
      help = "Create a service in an existing Sialia project"
      action = :command
    "infrastructure"
      help = "Create a new infrastructure definition in an existing Sialia project"
      action = :command
  end

  @add_arg_table! settings["build"] begin
    # Empty for now
  end

  @add_arg_table! settings["run"] begin
    # Empty for now
  end

  @add_arg_table! settings["deploy"] begin
    # Empty for now
  end

  return parse_args(settings, as_symbols=true)

end


text:
println("Parsing cross-sections")
cs = parse_cs(creek)

println("Loading cross-section data")
csd = load_data(cs)


text:
using Documenter
using CosmoCentral
using Plots
using FisherPlot
using PlotThemes
using CosmoCentralDocs

ENV["GKSwstype"] = "100"

push!(LOAD_PATH,"../src/")

makedocs(
    modules = [CosmoCentralDocs, CosmoCentral],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
    sidebar_sitename=false),
    sitename = "CosmoCentralDocs.jl",
    authors  = "Marco Bonici",
    pages = [
        "Home" => "index.md",
        "Cosmology" => "Cosmology.md",
        "Background Universe" => "BackgroundUniverse.md",
        "Source Density" => "SourceDensity.md",
        "Bias" => "Bias.md",
        "Boltzmann Solver" => "BoltzmannSolver.md",
        "Weight Functions" => "WeightFunction.md",
        "Angular Coefficients" => "AngularCoefficients.md",
        "Derivatives" => "Derivatives.md",
        "Covariance Matrix" => "Covariance.md",
        "Math Utils" => "MathUtils.md",
	"Fisher" => "Fisher.md"
    ]
)

deploydocs(
    repo = "github.com/marcobonici/CosmoCentralDocs.jl.git",
    devbranch = "main"
)


text:
"""
Providing *naive generalized generated* functions.
"""
module NGG
export to_type, from_type, TypeLevel
export compress, decompress
export RuntimeFn, Unset, Argument, mkngg, rmlines
using DataStructures
const List = LinkedList

rmlines(ex::Expr) = begin
    hd = ex.head
    tl = map(rmlines, filter(!islinenumbernode, ex.args))
    Expr(hd, tl...)
end
rmlines(@nospecialize(a)) = a
islinenumbernode(@nospecialize(x)) = x isa LineNumberNode

include("typeable.jl")
include("runtime_fns.jl")

"""
julia> using .NGG
julia> mkngg(
           :f, #fname
           [
               Argument(:a, nothing, Unset()),
               Argument(:b, nothing, Unset())
           ],  # args
           Argument[], # kwargs
           :(a + b) #expression
       )
f = (a, b;) -> a + b

julia> ans(1, 2)
3
"""
function mkngg(
    name::Symbol,
    args::Vector{Argument},
    kwargs::Vector{Argument},
    @nospecialize(ex)
)
    Args = to_type(_typed_list(Argument, args...))
    Kwargs = to_type(_typed_list(Argument, kwargs...))
    Ex = to_type(ex)
    RuntimeFn{Args,Kwargs,Ex,name}()
end

# const f = mkngg(
#     :f, #fname
#     [
#         Argument(:a, nothing, Unset()),
#         Argument(:b, nothing, Unset())
#     ],  # args
#     Argument[], # kwargs
#     :(a + b) #expression
# )

# println(f)
# println(typeof(f))
# using BenchmarkTools

# @btime f(1, 2)

end


text:
using Interpolations: LinearInterpolation, Extrapolation
using DataStructures: OrderedDict

struct Interpolate{V} <: State{V}
    #TODO: check performance of non-concrete type
    value::Extrapolation{V}
end

Interpolate(; unit, knotunit, reverse, _value, _type, _...) = begin
    l = if _value isa Extrapolation
        i = _value.itp
        if reverse
            #HACK: reverse interpolation
            zip(i.coefs, i.knots[1])
        else
            zip(i.knots[1], i.coefs)
        end
    else
        if _value isa Matrix
            zip(_value[:, 1], _value[:, 2])
        else
            _value
        end
    end
    d = OrderedDict(l)
    sort!(d)
    K = unitfy(collect(keys(d)), value(knotunit))
    V = unitfy(collect(values(d)), value(unit))
    v = LinearInterpolation(K, V)
    #HACK: pick up unitfy-ed valuetype
    V = typeof(v).parameters[1]
    Interpolate{V}(v)
end

constructortags(::Val{:Interpolate}) = (:unit, :knotunit, :reverse)

updatetags!(d, ::Val{:Interpolate}; _...) = begin
    !haskey(d, :reverse) && (d[:reverse] = false)
    !haskey(d, :knotunit) && (d[:knotunit] = nothing)
end

genvartype(v::VarInfo, ::Val{:Interpolate}; V, U, _...) = @q Interpolate{$V}

gendefault(v::VarInfo, ::Val{:Interpolate}) = gendefaultvalue(v, parameter=true, unitfy=false)

genupdate(v::VarInfo, ::Val{:Interpolate}, ::MainStep; kw...) = nothing


text:
export HomotopyNullCache,
    nvariables,
    cache,
    evaluate!,
    evaluate,
    jacobian!, jacobian,
    evaluate_and_jacobian!, evaluate_and_jacobian,
    dt!, dt,
    jacobian_and_dt!, jacobian_and_dt,
    basehomotopy

# Cache

"""
    HomotopyNullCache

The default `AbstractHomotopyCache` containing nothing.
"""
struct HomotopyNullCache <: AbstractHomotopyCache end
const HNC = HomotopyNullCache

"""
    cache(H::AbstractHomotopy, x, t)::AbstractHomotopyCache

Create a cache for the evaluation (incl. Jacobian) of `F` with elements of the type
of `x`. The default implementation returns [`HomotopyNullCache`](@ref).
"""
cache(H::AbstractHomotopy, x, t) = HomotopyNullCache()

# Homotopy API

"""
    evaluate!(u, H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)

Evaluate the homotopy `H` at `(x, t)` and store the result in `u`.
"""
evaluate!(u, H::AbstractHomotopy, args...) = error(MethodError(evaluate!, tuple(u, H, args...)))

"""
    evaluate(H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)

Evaluate the homotopy `H` at `(x, t)`.
"""
evaluate(H::AbstractHomotopy, x, t) = evaluate(H, x, t, cache(H, x, t))
function evaluate(H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)
    u = fill(zero(eltype(x)), size(H)[1])
    evaluate!(u, H, x, t, cache)
    u
end

"""
    dt!(u, H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)

Evaluate the homotopy `H` at `(x, t)` and store the result in `u`.
"""
function dt! end

"""
    dt(H::AbstractHomotopy, x::AbstractVector, cache::AbstractHomotopyCache)

Evaluate the homotopy `H` at `(x, t)`.
"""
dt(H::AbstractHomotopy, x, t) = dt(H, x, t, cache(H, x, t))
function dt(H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)
    u = fill(zero(eltype(x)), size(H)[1])
    dt!(u, H, x, t, cache)
    u
end

"""
    jacobian!(u, H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)

Evaluate the Jacobian of the homotopy `H` at `(x, t)` and store the result in `u`.
"""
jacobian!(u, H::AbstractHomotopy, args...) = error(MethodError(jacobian!, tuple(u, H, args...)))

"""
    jacobian(H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)

Evaluate the Jacobian of the homotopy `H` at `(x, t)`.
"""
jacobian(H::AbstractHomotopy, x, t) = jacobian(H, x, t, cache(H, x, t))
function jacobian(H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)
    U = fill(zero(eltype(x)), size(H))
    jacobian!(U, H, x, t, cache)
    U
end

# Optional
"""
    evaluate_and_jacobian!(u, U, F, x, t, cache::AbstractHomotopyCache)

Evaluate the homotopy `H` and its Jacobian at `(x, t)` and store the results in `u` (evalution)
and `U` (Jacobian).
"""
function evaluate_and_jacobian!(u, U, H::AbstractHomotopy, x, t, c=cache(H, x, t))
    evaluate!(u, H, x, t, c)
    jacobian!(U, H, x, t, c)
    nothing
end

"""
    evaluate_and_jacobian(H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)

Evaluate the homotopy `H` and its Jacobian at `(x, t)`.
"""
function evaluate_and_jacobian(H::AbstractHomotopy, x, t, c=cache(H, x, t))
    u = evaluate(H, x, t, c)
    U = jacobian(H, x, t, c)
    u, U
end

"""
    jacobian_and_dt!(U, u, H, x, t, cache::AbstractHomotopyCache)

Evaluate the homotopy `H` and its derivative w.r.t. `t` at `(x, t)` and store the results in `u` (evalution)
and `v` (∂t).
"""
function jacobian_and_dt!(U, u, H::AbstractHomotopy, x, t, c=cache(H, x, t))
    jacobian!(U, H, x, t, c)
    dt!(u, H, x, t, c)
    nothing
end

"""
    jacobian_and_dt(H::AbstractHomotopy, x, t, cache::AbstractHomotopyCache)

Evaluate the homotopy `H` and its derivative w.r.t. `t` at `(x, t)`.
"""
function jacobian_and_dt(H::AbstractHomotopy, x, t, c=cache(H, x, t))
    U = jacobian(H, x, t, c)
    u = dt(H, x, t, c)
    U, u
end

"""
    nvariables(H::AbstractHomotopy)

Returns the number of variables of the homotopy `H`.
"""
nvariables(H::AbstractHomotopy) = last(size(H))

Base.size(F::AbstractHomotopy, i::Integer) = size(F)[i]
Base.length(F::AbstractHomotopy) = size(F, 1)

"""
    basehomotopy(H::AbstractHomotopy)

Returns the 'proper' homotopy describing the problem. Any wrapper homotopy
recursively calls `wrappedhomotopy` with the wrapped homotopy as argument.
"""
basehomotopy(H::AbstractHomotopy) = H

"""
    Base.size(H::AbstractHomotopy)

Returns a tuple `(m, n)` indicating that `H` is a homotopy of `m` polynomials `m` in `n` variables.
"""
Base.size(::AbstractHomotopy) = error("Obligatory to define `Base.size($H)`")

include("homotopies/homotopy_witch_cache.jl")
include("homotopies/straight_line.jl")
include("homotopies/parameter_homotopy.jl")
include("homotopies/fixed_point.jl")
include("homotopies/patched_homotopy.jl")
include("homotopies/log_homotopy.jl")
include("homotopies/coefficient_homotopy.jl")
include("homotopies/toric_homotopy.jl")
include("homotopies/binomial_homotopy.jl")

function homotopy_interface_test(H::AbstractHomotopy, x=rand(Complex{Float64}, size(H, 2)) )
    m, n = size(H)
    t = rand()
    u = zeros(Complex{Float64}, m)
    U = zeros(Complex{Float64}, m, n)

    homotopy_cache = cache(H, x, t)

    evaluate!(u, H, x, t, homotopy_cache)
    @test evaluate(H, x, t, homotopy_cache) ≈ u atol=1e-14
    @test evaluate(H, x, t) ≈ u atol=1e-14

    dt!(u, H, x, t, homotopy_cache)
    @test dt(H, x, t, homotopy_cache) ≈ u atol=1e-14
    @test dt(H, x, t) ≈ u atol=1e-14

    jacobian!(U, H, x, t, homotopy_cache)
    @test jacobian(H, x, t, homotopy_cache) ≈ U atol=1e-14
    @test jacobian(H, x, t) ≈ U atol=1e-14

    evaluate_and_jacobian!(u, U, H, x, t, homotopy_cache)
    (v, V) = evaluate_and_jacobian(H, x, t, homotopy_cache)
    @test v ≈ u
    @test V ≈ U
    @test evaluate(H, x, t) ≈ u atol=1e-14
    @test jacobian(H, x, t) ≈ U atol=1e-14

    jacobian_and_dt!(U, u, H, x, t, homotopy_cache)
    (V, v) = jacobian_and_dt(H, x, t, homotopy_cache)
    @test V ≈ U atol=1e-14
    @test v ≈ u atol=1e-14
    @test jacobian(H, x, t) ≈ U atol=1e-14
    @test dt(H, x, t) ≈ u atol=1e-14
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI.FMIImport.FMICore: fmi2Real

myFMU = fmiLoad("SpringPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])
comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknwon fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

@test fmiSetupExperiment(fmuStruct) == 0
@test fmiEnterInitializationMode(fmuStruct) == 0
@test fmiExitInitializationMode(fmuStruct) == 0

targetValues = [[0.0, -10.0], [1.0, 0.0]]
dir_ders_buffer = zeros(fmi2Real, 2)
sample_ders_buffer = zeros(fmi2Real, 2, 1)
for i in 1:fmiGetNumberOfStates(myFMU)

    if fmiProvidesDirectionalDerivative(myFMU)
        # multi derivatives calls
        sample_ders = fmiSampleDirectionalDerivative(fmuStruct, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]])
        fmiSampleDirectionalDerivative!(fmuStruct, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]], sample_ders_buffer)

        @test sum(abs.(sample_ders[:,1] - targetValues[i])) < 1e-3
        @test sum(abs.(sample_ders_buffer[:,1] - targetValues[i])) < 1e-3

        dir_ders = fmiGetDirectionalDerivative(fmuStruct, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]])
        @test fmiGetDirectionalDerivative!(fmuStruct, myFMU.modelDescription.derivativeValueReferences, [myFMU.modelDescription.stateValueReferences[i]], dir_ders_buffer) == 0

        @test sum(abs.(dir_ders - targetValues[i])) < 1e-3
        @test sum(abs.(dir_ders_buffer - targetValues[i])) < 1e-3

        # single derivative call
        dir_der = fmiGetDirectionalDerivative(fmuStruct, myFMU.modelDescription.derivativeValueReferences[1], myFMU.modelDescription.stateValueReferences[1])
        @test dir_der == targetValues[1][1]

    else
        @warn "Skipping directional derivative testing, FMU from $(ENV["EXPORTINGTOOL"]) doesn't support directional derivatives."
    end

end

fmiUnload(myFMU)


text:
using Test
using ITensors
using LinearAlgebra

using ITensors.Ops: α, ∏, ∑, expand

function heisenberg(N)
  os = ∑{Op}()
  for j in 1:(N - 1)
    os += "Sz", j, "Sz", j + 1
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
  end
  return os
end

@testset "Ops" begin
  x1 = Op("X", 1)
  x2 = Op("X", 2)
  I1 = Op(I, 1)
  I2 = Op(I, 2)
  y1 = Op("Y", 1)
  y2 = Op("Y", 2)
  CX12 = Op("CX", 1, 2)
  Ry4 = Op("Ry", 4; θ=π / 3)

  @test 2y2 isa α{Op}
  @test coefficient(2y2) == 2
  @test y2 / 2 isa α{Op}
  @test coefficient(y2 / 2) ≈ 0.5
  @test -y2 isa α{Op}
  @test 1y2 + x1 isa ∑{<:α{Op}}
  @test 1y2 + x1 isa ∑{α{Op,Int}}
  @test x1 * y2 isa ∏{Op}
  @test 2x1 * y2 isa α{∏{Op}}
  @test x1 * y2 + CX12 isa ∑{∏{Op}}
  @test x1 * y2 + x1 * CX12 isa ∑{∏{Op}}
  @test x1 * y2 + 2CX12 isa ∑{<:α{∏{Op}}}
  @test x1 * y2 - CX12 isa ∑{<:α{∏{Op}}}
  @test 2x1 * y2 + 2CX12 isa ∑{<:α{∏{Op}}}
  @test 2x1 * y2 - 2CX12 isa ∑{<:α{∏{Op}}}
  @test (2x1 * y2 - 2CX12) / 3 isa ∑{<:α{∏{Op}}}

  o = (2x1 * y2 - 2CX12) / 3
  @test coefficient(o[1]) ≈ 2 / 3
  @test coefficient(o[2]) ≈ -2 / 3

  N = 4
  s = siteinds("Qubit", N)
  t1 = ITensor(x1, s)
  @test hassameinds(t1, (s[1]', dag(s[1])))
  @test t1[1, 1] == 0
  @test t1[1, 2] == 1

  @test ITensor(2.3x1, s) ≈ 2.3 * t1
  @test ITensor(x1 + x1, s) ≈ 2t1
  @test ITensor(x1 + 2.3x1, s) ≈ 3.3t1

  @test ITensor(Op(I, 2), s) ≈ ITensor([1 0; 0 1], s[2]', dag(s[2]))

  c = x1 * y2 * CX12
  cdag = c'
  @test c[1]' == cdag[3]
  @test c[2]' == cdag[2]
  @test c[3]' == cdag[1]

  x = randn(2, 2)
  tx = ITensor(Op(x, 3), s)
  @test tx[s[3]' => 1, s[3] => 2] == x[1, 2]

  @test ITensor(x1 * x1, s) ≈ ITensor(Op([1 0; 0 1], 1), s)
  @test ITensor(x1 * x1 * x1, s) ≈ ITensor(Op([0 1; 1 0], 1), s)
  @test ITensor(2x1 * x1, s) ≈ ITensor(Op([2 0; 0 2], 1), s)
  @test ITensor(x1 * y1, s) ≈ ITensor(Op([im 0; 0 -im], 1), s)
  @test ITensor(y1 * x1, s) ≈ ITensor(Op([-im 0; 0 im], 1), s)
  @test ITensor(2x1 * x1 + y1, s) ≈
    ITensor(2 * [1 0; 0 1] + [0 -im; im 0], s[1]', dag(s[1]))
  @test ITensor(2y1 * x2 + x1, s) ≈
    2 * ITensor(y1, s) * ITensor(x2, s) + ITensor(x1, s) * ITensor(I2, s)

  @test y1'' == y1

  @test ITensor(y1', s) ≈ ITensor(Op([0 -im; im 0], 1), s)

  @test ITensor(exp(x1), s) ≈ ITensor(Op(exp([0 1; 1 0]), 1), s)
  @test ITensor(exp(2x1 * x1), s) ≈ ITensor(exp(2 * [1 0; 0 1]), s[1]', dag(s[1]))
  @test ITensor(exp(2x1 * x1 + y1), s) ≈
    ITensor(exp(2 * [1 0; 0 1] + [0 -im; im 0]), s[1]', dag(s[1]))

  @test ITensor(I1, s) ≈ ITensor([1 0; 0 1], s[1]', dag(s[1]))

  @test exp(Op("X", 1)) * Op("Y", 2) isa ∏{Any}
  @test ITensor(exp(Op("X", 1)) * Op("Y", 1), s) ≈
    product(exp(ITensor(Op("X", 1), s)), ITensor(Op("Y", 1), s))
  @test 2exp(Op("X", 1)) * Op("Y", 2) isa α{∏{Any}}

  H = ∑{<:α{∏{Op}}}() - Op("X", 1)
  @test H isa ∑
  @test H isa ∑{<:α}
  @test H isa ∑{<:α{<:∏}}
  @test H isa ∑{<:α{∏{Op}}}
  @test H isa ∑{α{∏{Op},T}} where {T}
  @test H isa ∑{α{∏{Op},Int}}
  @test length(H) == 1
  @test coefficient(H[1]) == -1

  H = ∑{Op}() - Op("X", 1)
  @test H isa ∑
  @test H isa ∑{<:α}
  @test H isa ∑{<:α{Op}}
  @test H isa ∑{α{Op,T}} where {T}
  @test H isa ∑{α{Op,Int}}
  @test length(H) == 1
  @test coefficient(H[1]) == -1

  # MPO conversion
  H = ∑{Op}()
  H -= 2.3, "X", 1, "X", 2
  H += 1.2, "Z", 1
  H += 1.3, "Z", 2, (θ=π / 3,)
  @test H isa ∑{α{∏{Op},Float64}}
  @test length(H) == 3
  @test coefficient(H[1]) == -2.3
  @test length(H[1]) == 2
  @test Ops.sites(H[1]) == [1, 2]
  @test coefficient(H[2]) == 1.2
  @test length(H[2]) == 1
  @test Ops.sites(H[2]) == [1]
  @test coefficient(H[3]) == 1.3
  @test length(H[3]) == 1
  @test Ops.sites(H[3]) == [2]
  @test Ops.params(H[3]) == (θ=π / 3,)

  @test ∑{Op}(("X", 1)) isa ∑{Op}
  @test ∑{Op}((2.3, "X", 1)) isa ∑{α{Op,Float64}}
  @test ∑{Op}("X", 1) isa ∑{Op}
  @test ∑{Op}(2, "X", 1) isa ∑{α{Op,Int}}
  @test ∑{Op}([Op("X", 1), 2Op("Y", 1)]) isa ∑
  @test ∑{Op}([Op("X", 1), 2Op("Y", 1)]) isa ∑{<:α}
  @test ∑{Op}([Op("X", 1), 2Op("Y", 1)]) isa ∑{<:α{Op}}
  @test ∑{Op}([Op("X", 1), 2Op("Y", 1)]) isa ∑{α{Op,Int}}

  @testset "Expand expression, 2 products" begin
    expr = (Op("X", 1) + Op("Y", 2)) * (Op("Z", 1) + Op("W", 2))
    expr_expanded =
      Op("X", 1) * Op("Z", 1) +
      Op("Y", 2) * Op("Z", 1) +
      Op("X", 1) * Op("W", 2) +
      Op("Y", 2) * Op("W", 2)
    @test expand(expr) == expr_expanded
  end

  @testset "Expand expression, 3 products" begin
    expr = (Op("X", 1) + Op("Y", 2)) * (Op("Z", 1) + Op("W", 2)) * (Op("A", 1) + Op("B", 2))
    expr_expanded =
      Op("X", 1) * Op("Z", 1) * Op("A", 1) +
      Op("Y", 2) * Op("Z", 1) * Op("A", 1) +
      Op("X", 1) * Op("W", 2) * Op("A", 1) +
      Op("Y", 2) * Op("W", 2) * Op("A", 1) +
      Op("X", 1) * Op("Z", 1) * Op("B", 2) +
      Op("Y", 2) * Op("Z", 1) * Op("B", 2) +
      Op("X", 1) * Op("W", 2) * Op("B", 2) +
      Op("Y", 2) * Op("W", 2) * Op("B", 2)
    @test expand(expr) == expr_expanded
  end

  H = heisenberg(4)
  @test length(H) == 9
  @test H^2 == H * H
  @test length(H^2) == 2
  @test length(expand(H^2)) == 81

  @testset "Conversion to Sum of ITensors" begin
    H = Sum{Op}() + ("X", 1) + ("Y", 2)
    @test H == Sum{Op}([("X", 1), ("Y", 2)])
    s = siteinds("Qubit", 2)
    Hₜ = Sum{ITensor}(H, s)
    @test Hₜ isa Sum{ITensor}
    @test Hₜ[1] ≈ ITensor(Op("X", 1), s)
    @test Hₜ[2] ≈ ITensor(Op("Y", 2), s)
  end

  @testset "Conversion to ∏ of ITensors" begin
    C = ∏{Op}() * ("X", 1) * ("Y", 2)
    @test C == ∏{Op}([("X", 1), ("Y", 2)])
    s = siteinds("Qubit", 2)
    Cₜ = ∏{ITensor}(C, s)
    @test Cₜ isa ∏{ITensor}
    @test Cₜ[1] ≈ ITensor(Op("X", 1), s)
    @test Cₜ[2] ≈ ITensor(Op("Y", 2), s)
  end
end


text:
Base.@kwdef struct MOEADDE <: AbstractOptimizer
    # number of the subproblems
    N::Int = 100
    # number of the weight vectors in the neighborhoor of each weight vector
    T::Int = ceil(N/10)
    crossover = BINX(0.5)
    mutation = PLM()
    metrics::ConvergenceMetrics = ConvergenceMetric[GD(), GD(true)]
end

population_size(method::MOEADDE) = method.N
default_options(method::MOEADDE) = (iterations = 1000,)
summary(m::MOEADDE) =
    "MOEA/D-DE[P=$(population_size(m)),]"
show(io::IO, m::MOEADDE) = print(io, summary(m))

mutable struct MOEADDEState <: AbstractOptimizerState
    iteration
    start_time
    stop_time
    metrics::ConvergenceMetrics # collection of convergence metrics

    population                  # population
    m                           # number of objectives
    W                           # weight vectors
    B                           # neighbours of each solution
    z                           # the best value found so far
    EP                          # nondominated solutions
end
pfront(s::MOEADDEState) = s.EP
value(s::MOEADDEState) = objectives(pfront(s))
minimizer(s::MOEADDEState) = objectives(pfront(s))
copy(s::MOEADDEState) = MOEADDEState(s.iteration, s.start_time, s.stop_time, copy(s.metrics),
                                 copy(s.population),
                                 copy(s.m),
                                 copy(s.W),
                                 copy(s.B),
                                 copy(s.z),
                                 copy(s.EP))

function initial_state(method::MOEADDEState, objective, population, options)
    value!(objective, population)

    m = length(objectives(first(population)))

    # Generate weight vectors
    W = uniform_point(population_size(method), m)

    # Set the neighbours of each weight vector
    B = pairwise(Euclidean(), W, W, dims=1)
    B = sortperm(B)
    B = B[:,2:method.T]

    # Set reference point
    z = Matrix(minimum(objectives(population), dims=1)')

    EP = Individual[]
    return MOEADDEState(0, 0, 0, copy(method.metrics), population, m, W, B, z, EP)
end

function update_state!(
    method::MOEADDE,
    state,
    objective,
    constraints,
    options
)
    for i = 1:method.N
        x = state.population[i]

        # Reproduction
        xi = rand(options.rng, state.B[i,:], 2)
        xs = state.population[xi]
        y,_ = method.crossover(xs[1], xs[2], rng=options.rng)

        # Improvement
        y = method.mutation(y)

        apply!(constraints, variables(y))
        value!(objective, y)

        # Update reference point
        for j = 1:state.m
            if state.z[j] > objectives(y)[j]
                state.z[j] = objectives(y)[j]
            end
        end

        # Update of neighboring solutions
        for j = state.B[i,:]
            x = state.population[j]
            g_x, _ = findmax(state.W[j,:] .* abs.(objectives(x) - state.z))
            g_y, _ = findmax(state.W[j,:] .* abs.(objectives(y) - state.z))
            if g_y <= g_x
                state.population[j] = y
            end
        end

        # Update of EP
        deleteat!(state.EP, findall((p) -> dominate(y, p) == 1, state.EP))
        if isempty(findall((p) -> dominate(p, y) == 1, state.EP))
            pushfirst!(state.EP, y)
        end
    end

    return false
end


text:
using ITensors, Test

@testset "Physics Sites" begin
  N = 10

  @testset "Qubit sites" begin
    s = siteind("Qubit")
    @test hastags(s, "Qubit,Site")
    @test dim(s) == 2

    s = siteinds("Qubit", N)
    @test val(s[1], "0") == 1
    @test val(s[1], "1") == 2
    @test_throws ArgumentError val(s[1], "Fake")

    s = siteind("Qubit"; conserve_parity=true)
    @test hastags(s, "Qubit,Site")
    @test dim(s) == 2
    @test nblocks(s) == 2
    @test qn(s, 1) == QN("Parity", 0, 2)
    @test qn(s, 2) == QN("Parity", 1, 2)

    s = siteind("Qubit"; conserve_number=true, conserve_parity=true)
    @test hastags(s, "Qubit,Site")
    @test dim(s) == 2
    @test nblocks(s) == 2
    @test qn(s, 1) == QN(("Parity", 0, 2), ("Number", 0))
    @test qn(s, 2) == QN(("Parity", 1, 2), ("Number", 1))

    s = siteinds("Qubit", N)

    Z = op("Z", s, 5)
    @test hasinds(Z, s[5]', s[5])

    @test_throws ArgumentError op(s, "Fake", 2)
    @test Array(op("Id", s, 3), s[3]', s[3]) ≈ [1.0 0.0; 0.0 1.0]
    @test Array(op("√NOT", s, 3), s[3]', s[3]) ≈
          [(1 + im)/2 (1 - im)/2; (1 - im)/2 (1 + im)/2]
    @test Array(op("H", s, 3), s[3]', s[3]) ≈ [1/sqrt(2) 1/sqrt(2); 1/sqrt(2) -1/sqrt(2)]
    @test Array(op("Phase", s, 3), s[3]', s[3]) ≈ [1 0; 0 im]
    @test Array(op("P", s, 3), s[3]', s[3]) ≈ [1 0; 0 im]
    @test Array(op("S", s, 3), s[3]', s[3]) ≈ [1 0; 0 im]
    @test Array(op("π/8", s, 3), s[3]', s[3]) ≈ [1 0; 0 (1 + im)/sqrt(2)]
    @test Array(op("T", s, 3), s[3]', s[3]) ≈ [1 0; 0 (1 + im)/sqrt(2)]
    θ = randn()
    @test Array(op("Rx", s, 3; θ=θ), s[3]', s[3]) ≈
          [cos(θ / 2) -im*sin(θ / 2); -im*sin(θ / 2) cos(θ / 2)]

    # Test obtaining S=1/2 operators using Qubit tag
    @test Array(op("X", s, 3), s[3]', s[3]) ≈ [0.0 1.0; 1.0 0.0]
  end

  @testset "Spin Half sites" begin
    s = siteind("S=1/2")
    @test hastags(s, "S=1/2,Site")
    @test dim(s) == 2

    s = siteind("S=1/2"; conserve_qns=true)
    @test hastags(s, "S=1/2,Site")
    @test dim(s) == 2
    @test nblocks(s) == 2
    @test qn(s, 1) == QN("Sz", +1)
    @test qn(s, 2) == QN("Sz", -1)

    s = siteind("S=1/2"; conserve_szparity=true)
    @test hastags(s, "S=1/2,Site")
    @test dim(s) == 2
    @test nblocks(s) == 2
    @test qn(s, 1) == QN("SzParity", 1, 2)
    @test qn(s, 2) == QN("SzParity", 0, 2)

    s = siteind("S=1/2"; conserve_sz=true, conserve_szparity=true)
    @test hastags(s, "S=1/2,Site")
    @test dim(s) == 2
    @test nblocks(s) == 2
    @test qn(s, 1) == QN(("SzParity", 1, 2), ("Sz", +1))
    @test qn(s, 2) == QN(("SzParity", 0, 2), ("Sz", -1))

    s = siteinds("S=1/2", N)
    @test val(s[1], "Up") == 1
    @test val(s[1], "Dn") == 2
    @test_throws ArgumentError val(s[1], "Fake")

    Sz5 = op("Sz", s, 5)
    @test hasinds(Sz5, s[5]', s[5])

    @test_throws ArgumentError op(s, "Fake", 2)
    @test Array(op("Id", s, 3), s[3]', s[3]) ≈ [1.0 0.0; 0.0 1.0]
    @test Array(op("F", s, 3), s[3]', s[3]) ≈ [1.0 0.0; 0.0 1.0]
    @test Array(op("S+", s, 3), s[3]', s[3]) ≈ [0.0 1.0; 0.0 0.0]
    @test Array(op("S⁺", s, 3), s[3]', s[3]) ≈ [0.0 1.0; 0.0 0.0]
    @test Array(op("S-", s, 4), s[4]', s[4]) ≈ [0.0 0.0; 1.0 0.0]
    @test Array(op("S⁻", s, 4), s[4]', s[4]) ≈ [0.0 0.0; 1.0 0.0]
    @test Array(op("Sx", s, 2), s[2]', s[2]) ≈ [0.0 0.5; 0.5 0.0]
    @test Array(op("Sˣ", s, 2), s[2]', s[2]) ≈ [0.0 0.5; 0.5 0.0]
    @test Array(op("iSy", s, 2), s[2]', s[2]) ≈ [0.0 0.5; -0.5 0.0]
    @test Array(op("iSʸ", s, 2), s[2]', s[2]) ≈ [0.0 0.5; -0.5 0.0]
    @test Array(op("Sy", s, 2), s[2]', s[2]) ≈ [0.0 -0.5im; 0.5im 0.0]
    @test Array(op("Sʸ", s, 2), s[2]', s[2]) ≈ [0.0 -0.5im; 0.5im 0.0]
    @test Array(op("Sz", s, 2), s[2]', s[2]) ≈ [0.5 0.0; 0.0 -0.5]
    @test Array(op("Sᶻ", s, 2), s[2]', s[2]) ≈ [0.5 0.0; 0.0 -0.5]
    @test Array(op("ProjUp", s, 2), s[2]', s[2]) ≈ [1.0 0.0; 0.0 0.0]
    @test Array(op("projUp", s, 2), s[2]', s[2]) ≈ [1.0 0.0; 0.0 0.0]
    @test Array(op("ProjDn", s, 2), s[2]', s[2]) ≈ [0.0 0.0; 0.0 1.0]
    @test Array(op("projDn", s, 2), s[2]', s[2]) ≈ [0.0 0.0; 0.0 1.0]

    # Test obtaining Qubit operators using S=1/2 tag:
    @test Array(op("√NOT", s, 3), s[3]', s[3]) ≈
          [(1 + im)/2 (1 - im)/2; (1 - im)/2 (1 + im)/2]
  end

  @testset "Spin One sites" begin
    s = siteinds("S=1", N)

    @test val(s[1], "Up") == 1
    @test val(s[1], "0") == 2
    @test val(s[1], "Dn") == 3
    @test_throws ArgumentError val(s[1], "Fake")

    Sz5 = op("Sz", s, 5)
    @test hasinds(Sz5, s[5]', s[5])

    @test_throws ArgumentError op(s, "Fake", 2)
    @test Array(op("Id", s, 3), s[3]', s[3]) ≈ [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    @test Array(op("S+", s, 3), s[3]', s[3]) ≈ [0 √2 0; 0 0 √2; 0 0 0]
    @test Array(op("S⁺", s, 3), s[3]', s[3]) ≈ [0 √2 0; 0 0 √2; 0 0 0]
    @test Array(op("S-", s, 3), s[3]', s[3]) ≈ [0 0 0; √2 0 0; 0.0 √2 0]
    @test Array(op("S⁻", s, 3), s[3]', s[3]) ≈ [0 0 0; √2 0 0; 0.0 √2 0]
    @test Array(op("Sx", s, 3), s[3]', s[3]) ≈ [0 1/√2 0; 1/√2 0 1/√2; 0 1/√2 0]
    @test Array(op("Sˣ", s, 3), s[3]', s[3]) ≈ [0 1/√2 0; 1/√2 0 1/√2; 0 1/√2 0]
    @test Array(op("iSy", s, 3), s[3]', s[3]) ≈ [0 1/√2 0; -1/√2 0 1/√2; 0 -1/√2 0]
    @test Array(op("iSʸ", s, 3), s[3]', s[3]) ≈ [0 1/√2 0; -1/√2 0 1/√2; 0 -1/√2 0]
    @test Array(op("Sy", s, 3), s[3]', s[3]) ≈ (1 / (√2im)) * [0 +1 0; -1 0 +1; 0 -1 0]
    @test Array(op("Sʸ", s, 3), s[3]', s[3]) ≈ (1 / (√2im)) * [0 +1 0; -1 0 +1; 0 -1 0]
    #@test Array(op("Sʸ", s, 3), s[3]', s[3]) ≈ [0 +1/√2im 0; +1/√2im 0 -1/√2im; 0 +1/√2im 0]
    @test Array(op("Sz", s, 2), s[2]', s[2]) ≈ [1.0 0 0; 0 0 0; 0 0 -1.0]
    @test Array(op("Sᶻ", s, 2), s[2]', s[2]) ≈ [1.0 0 0; 0 0 0; 0 0 -1.0]
    @test Array(op("Sz2", s, 2), s[2]', s[2]) ≈ [1.0 0 0; 0 0 0; 0 0 +1.0]
    @test Array(op("Sx2", s, 2), s[2]', s[2]) ≈ [0.5 0 0.5; 0 1.0 0; 0.5 0 0.5]
    @test Array(op("Sy2", s, 2), s[2]', s[2]) ≈ [0.5 0 -0.5; 0 1.0 0; -0.5 0 0.5]
  end

  @testset "Fermion sites" begin
    s = siteind("Fermion")

    @test val(s, "0") == 1
    @test val(s, "1") == 2
    @test_throws ArgumentError val(s, "Fake")

    N = op(s, "N")
    @test hasinds(N, s', s)

    @test_throws ArgumentError op(s, "Fake")
    N = Array(op(s, "N"), s', s)
    @test N ≈ [0.0 0; 0 1]
    C = Array(op(s, "C"), s', s)
    @test C ≈ [0.0 1; 0 0]
    Cdag = Array(op(s, "Cdag"), s', s)
    @test Cdag ≈ [0.0 0; 1 0]
    F = Array(op(s, "F"), s', s)
    @test F ≈ [1.0 0; 0 -1]

    @test has_fermion_string("C", s)
    @test has_fermion_string("Cdag", s)
    @test has_fermion_string("C*F", s)
    @test has_fermion_string("F*Cdag*F", s)
    @test !has_fermion_string("N", s)
    @test !has_fermion_string("N*F", s)

    s = siteind("Fermion"; conserve_nf=true)
    @test qn(s, 1) == QN("Nf", 0, -1)
    @test qn(s, 2) == QN("Nf", 1, -1)
    s = siteind("Fermion"; conserve_nfparity=true)
    @test qn(s, 1) == QN("NfParity", 0, -2)
    @test qn(s, 2) == QN("NfParity", 1, -2)
    s = siteind("Fermion"; conserve_qns=false)
    @test dim(s) == 2

    s = siteind("Fermion"; conserve_nf=true, conserve_sz=true)
    @test qn(s, 1) == QN(("Nf", 0, -1), ("Sz", 0))
    @test qn(s, 2) == QN(("Nf", 1, -1), ("Sz", 1))
    s = siteind("Fermion"; conserve_nfparity=true, conserve_sz=true)
    @test qn(s, 1) == QN(("NfParity", 0, -2), ("Sz", 0))
    @test qn(s, 2) == QN(("NfParity", 1, -2), ("Sz", 1))
    s = siteind("Fermion"; conserve_nf=true, conserve_sz="Up")
    @test qn(s, 1) == QN(("Nf", 0, -1), ("Sz", 0))
    @test qn(s, 2) == QN(("Nf", 1, -1), ("Sz", 1))
    s = siteind("Fermion"; conserve_nfparity=true, conserve_sz="Up")
    @test qn(s, 1) == QN(("NfParity", 0, -2), ("Sz", 0))
    @test qn(s, 2) == QN(("NfParity", 1, -2), ("Sz", 1))
    s = siteind("Fermion"; conserve_nf=true, conserve_sz="Dn")
    @test qn(s, 1) == QN(("Nf", 0, -1), ("Sz", 0))
    @test qn(s, 2) == QN(("Nf", 1, -1), ("Sz", -1))
    s = siteind("Fermion"; conserve_nfparity=true, conserve_sz="Dn")
    @test qn(s, 1) == QN(("NfParity", 0, -2), ("Sz", 0))
    @test qn(s, 2) == QN(("NfParity", 1, -2), ("Sz", -1))
  end

  @testset "Electron sites" begin
    s = siteind("Electron")

    @test val(s, "0") == 1
    @test val(s, "Up") == 2
    @test val(s, "Dn") == 3
    @test val(s, "UpDn") == 4
    @test_throws ArgumentError val(s, "Fake")

    Nup = op(s, "Nup")
    @test hasinds(Nup, s', s)

    @test_throws ArgumentError op(s, "Fake")
    Nup = Array(op(s, "Nup"), s', s)
    @test Nup ≈ [0.0 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]
    Ndn = Array(op(s, "Ndn"), s', s)
    @test Ndn ≈ [0.0 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1]
    Ntot = Array(op(s, "Ntot"), s', s)
    @test Ntot ≈ [0.0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 2]
    Cup = Array(op(s, "Cup"), s', s)
    @test Cup ≈ [0.0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0]
    Cdagup = Array(op(s, "Cdagup"), s', s)
    @test Cdagup ≈ [0.0 0 0 0; 1 0 0 0; 0 0 0 0; 0 0 1 0]
    Cdn = Array(op(s, "Cdn"), s', s)
    @test Cdn ≈ [0.0 0 1 0; 0 0 0 -1; 0 0 0 0; 0 0 0 0]
    Cdagdn = Array(op(s, "Cdagdn"), s', s)
    @test Cdagdn ≈ [0.0 0 0 0; 0 0 0 0; 1 0 0 0; 0 -1 0 0]
    F = Array(op(s, "F"), s', s)
    @test F ≈ [1.0 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1]
    Fup = Array(op(s, "Fup"), s', s)
    @test Fup ≈ [1.0 0 0 0; 0 -1 0 0; 0 0 1 0; 0 0 0 -1]
    Fdn3 = Array(op(s, "Fdn"), s', s)
    @test Fdn3 ≈ [1.0 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 -1]
    Sz3 = Array(op(s, "Sz"), s', s)
    @test Sz3 ≈ [0.0 0 0 0; 0 0.5 0 0; 0 0 -0.5 0; 0 0 0 0]
    Sx3 = Array(op(s, "Sx"), s', s)
    @test Sx3 ≈ [0.0 0 0 0; 0 0 0.5 0; 0 0.5 0 0; 0 0 0 0]
    Sp3 = Array(op(s, "S+"), s', s)
    @test Sp3 ≈ [0.0 0 0 0; 0 0 1 0; 0 0 0 0; 0 0 0 0]
    Sm3 = Array(op(s, "S-"), s', s)
    @test Sm3 ≈ [0.0 0 0 0; 0 0 0 0; 0 1 0 0; 0 0 0 0]

    @test has_fermion_string("Cup", s)
    @test has_fermion_string("Cup*F", s)
    @test has_fermion_string("Cdagup", s)
    @test has_fermion_string("F*Cdagup", s)
    @test has_fermion_string("Cdn", s)
    @test has_fermion_string("Cdn*F", s)
    @test has_fermion_string("Cdagdn", s)
    @test !has_fermion_string("N", s)
    @test !has_fermion_string("F*N", s)

    s = siteind("Electron"; conserve_nf=true)
    @test qn(s, 1) == QN("Nf", 0, -1)
    @test qn(s, 2) == QN("Nf", 1, -1)
    @test qn(s, 3) == QN("Nf", 2, -1)
    s = siteind("Electron"; conserve_sz=true)
    @test qn(s, 1) == QN(("Sz", 0), ("NfParity", 0, -2))
    @test qn(s, 2) == QN(("Sz", +1), ("NfParity", 1, -2))
    @test qn(s, 3) == QN(("Sz", -1), ("NfParity", 1, -2))
    @test qn(s, 4) == QN(("Sz", 0), ("NfParity", 0, -2))
    s = siteind("Electron"; conserve_nfparity=true)
    @test qn(s, 1) == QN("NfParity", 0, -2)
    @test qn(s, 2) == QN("NfParity", 1, -2)
    @test qn(s, 3) == QN("NfParity", 0, -2)
    s = siteind("Electron"; conserve_qns=false)
    @test dim(s) == 4
  end

  @testset "tJ sites" begin
    s = siteind("tJ")

    @test val(s, "0") == 1
    @test val(s, "Up") == 2
    @test val(s, "Dn") == 3
    @test_throws ArgumentError val(s, "Fake")

    @test_throws ArgumentError op(s, "Fake")
    Nup = op(s, "Nup")
    @test Nup[2, 2] ≈ 1.0
    Ndn = op(s, "Ndn")
    @test Ndn[3, 3] ≈ 1.0
    Ntot = op(s, "Ntot")
    @test Ntot[2, 2] ≈ 1.0
    @test Ntot[3, 3] ≈ 1.0
    Id = Array(op(s, "Id"), s', s)
    @test Id ≈ [1.0 0 0; 0 1 0; 0 0 1]
    Cup = Array(op(s, "Cup"), s', s)
    @test Cup ≈ [0.0 1 0; 0 0 0; 0 0 0]
    Cdup = Array(op(s, "Cdagup"), s', s)
    @test Cdup ≈ [0 0 0; 1.0 0 0; 0 0 0]
    Cdn = Array(op(s, "Cdn"), s', s)
    @test Cdn ≈ [0.0 0.0 1; 0 0 0; 0 0 0]
    Cddn = Array(op(s, "Cdagdn"), s', s)
    @test Cddn ≈ [0 0 0; 0.0 0 0; 1 0 0]
    FP = Array(op(s, "F"), s', s)
    @test FP ≈ [1.0 0.0 0; 0 -1.0 0; 0 0 -1.0]
    Fup = Array(op(s, "Fup"), s', s)
    @test Fup ≈ [1.0 0.0 0; 0 -1.0 0; 0 0 1.0]
    Fdn = Array(op(s, "Fdn"), s', s)
    @test Fdn ≈ [1.0 0.0 0; 0 1.0 0; 0 0 -1.0]
    Sz = Array(op(s, "Sz"), s', s)
    @test Sz ≈ [0.0 0.0 0; 0 0.5 0; 0 0 -0.5]
    Sx = Array(op(s, "Sx"), s', s)
    @test Sx ≈ [0.0 0.0 0; 0 0 0.5; 0 0.5 0]
    Sp = Array(op(s, "Splus"), s', s)
    @test Sp ≈ [0.0 0.0 0; 0 0 1.0; 0 0 0]
    Sm = Array(op(s, "Sminus"), s', s)
    @test Sm ≈ [0.0 0.0 0; 0 0 0; 0 1.0 0]

    @test has_fermion_string("Cup", s)
    @test has_fermion_string("Cdagup", s)
    @test has_fermion_string("Cdn", s)
    @test has_fermion_string("Cdagdn", s)
    @test !has_fermion_string("N", s)
  end
end

nothing


text:
export Periodic , Bloch, DirichletNeumann , NeumannDirichlet , NeumannNeumann , DirichletDirichlet

abstract type BCType end
struct Periodic <: BCType end
struct Bloch <: BCType end
struct DirichletNeumann <: BCType end
struct NeumannDirichlet <: BCType end
struct NeumannNeumann <: BCType end
struct DirichletDirichlet <: BCType end


text:
for i in 10:-1:0
    println(i)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Feather, DataFrames, StatsBase, BoskoDB
import DataFrames.DataFrame

using JuliaDB
addprocs()
@everywhere using JuliaDB, Dagger, BoskoDB

const lc_feather_path = "C:/Users/dzj/Documents/lc.feather"
const lc_csv_path = "C:/Users/dzj/Downloads/Lending Club/csv/"

max_stars_repo_path: src/julia/lending club/0_setup.jl

max_issues_repo_path: src/julia/lending club/0_setup.jl

max_forks_repo_path: src/julia/lending club/0_setup.jl


text:
### Julia OpenStreetMapXPlot Package ###
### MIT License                 ###

### Standard map display "layers." ###

const LAYER_STANDARD = Dict(
    1 => OpenStreetMapXPlot.Style("0xd97486", 4), # @motorway-fill
    2 => OpenStreetMapXPlot.Style("0xebb36a", 3), # @primary-fill
    3 => OpenStreetMapXPlot.Style("0xe8b995", 3), # @secondary-fill
    4 => OpenStreetMapXPlot.Style("0xecf27c", 3), # @tertiary-fill
    5 => OpenStreetMapXPlot.Style("0xb9b9b6", 2), # gray
    6 => OpenStreetMapXPlot.Style("0xc7c7c4", 2), # gray
    7 => OpenStreetMapXPlot.Style("0xd9d9d5", 2), # Light gray
    8 => OpenStreetMapXPlot.Style("0xd9d9d5", 1)) # Light gray

const LAYER_CYCLE = Dict(
    1 => OpenStreetMapXPlot.Style("0x006600", 3), # Green
    2 => OpenStreetMapXPlot.Style("0x5C85FF", 3), # Blue
    3 => OpenStreetMapXPlot.Style("0x5C85FF", 2), # Blue
    4 => OpenStreetMapXPlot.Style("0x999999", 2)) # Dark gray

const LAYER_PED = Dict(
    1 => OpenStreetMapXPlot.Style("0x999999", 3), # Dark gray
    2 => OpenStreetMapXPlot.Style("0x999999", 3), # Dark gray
    3 => OpenStreetMapXPlot.Style("0x999999", 2), # Dark gray
    4 => OpenStreetMapXPlot.Style("0xE0E0E0", 2)) # Light gray

const LAYER_FEATURES = Dict(
    1 => OpenStreetMapXPlot.Style("0x9966FF", 1.5, "."),  # Lavender
    2 => OpenStreetMapXPlot.Style("0xFF0000", 1.5, "."),  # Red
    3 => OpenStreetMapXPlot.Style("0x000000", 1.5, "."),  # Black
    4 => OpenStreetMapXPlot.Style("0xFF66FF", 1.5, "."),  # Pink
    5 => OpenStreetMapXPlot.Style("0x996633", 1.5, "."),  # Brown
    6 => OpenStreetMapXPlot.Style("0xFF9900", 2.0, "."),  # Orange
    7 => OpenStreetMapXPlot.Style("0xCC00CC", 1.5, "."),  # Brown
	8 => OpenStreetMapXPlot.Style("0xFFFF00", 1.5, "."),  # Yellow
    9 => OpenStreetMapXPlot.Style("0xF4CCCC", 1.5, "."),  # Vanilla Ice
    10 => OpenStreetMapXPlot.Style("0x351C75", 1.5, "."), # Persian Indigo
    11 => OpenStreetMapXPlot.Style("0x00FF00", 1.5, "."), # Lime
    12 => OpenStreetMapXPlot.Style("0x00FFFF", 1.5, "."), # Aqua
    13 => OpenStreetMapXPlot.Style("0x005353", 2.0, "."), # Sherpa Blue
    14 => OpenStreetMapXPlot.Style("0xBDAD7D", 1.5, "."), # Ecru
    15 => OpenStreetMapXPlot.Style("0xFF00FF", 1.5, "."), # Fuchsia
    16 => OpenStreetMapXPlot.Style("0xB9D1D6", 1.5, "."), # Light Blue
    17 => OpenStreetMapXPlot.Style("0x7A00CC", 1.5, "."), # Violet
    18 => OpenStreetMapXPlot.Style("0x004225", 1.5, "."), # British Racing Green
    19 => OpenStreetMapXPlot.Style("0x7E8386", 2.0, "."), # Silver
    20 => OpenStreetMapXPlot.Style("0xB87333", 1.5, "."), # Copper
    21 => OpenStreetMapXPlot.Style("0x800020", 1.5, "."), # Burgundy
    22 => OpenStreetMapXPlot.Style("0x5B718D", 1.5, "."), # Ultramarine
    23 => OpenStreetMapXPlot.Style("0x636F22", 1.5, "."), # Fiji Green
    24 => OpenStreetMapXPlot.Style("0xCAF4DF", 1.5, "."), # Mint
    25 => OpenStreetMapXPlot.Style("0x231F66", 2.0, "."), # Midnight Blue
    26 => OpenStreetMapXPlot.Style("0xE6DFE7", 1.5, ".")) # Selago

const LAYER_BUILDINGS = Dict(
    1 => OpenStreetMapXPlot.Style("0xE1E1EB", 1, "-"), # Lighter gray
    2 => OpenStreetMapXPlot.Style("0xB8DBFF", 1, "-"), # Light blue
    3 => OpenStreetMapXPlot.Style("0xB5B5CE", 1, "-"), # Light gray
    4 => OpenStreetMapXPlot.Style("0xFFFF99", 1, "-"), # Pale yellow
    5 => OpenStreetMapXPlot.Style("0x006600", 1, "-")) # Green


text:
@reexport module SAA

using JuMP
using MathOptInterface
using StochasticPrograms
using ProgressMeter
using Parameters

import StochasticPrograms: load_model!, optimizer_name, optimal_instance

const MOI = MathOptInterface
const MOIU = MOI.Utilities

include("solver.jl")
include("MOI_wrapper.jl")

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SLIC(k, m; tol=1e-4, maxiter=10, vars=nothing)

A method for partitioning spatial data into approximately `k`
clusters using Simple Linear Iterative Clustering (SLIC).
The method produces clusters of samples that are spatially
connected based on a distance `dₛ` and that, at the same
time, are similar in terms of `vars` with distance `dᵥ`.
The tradeoff is controlled with a hyperparameter parameter
`m` in an additive model `dₜ = √(dᵥ² + m²(dₛ/s)²)`.

## Parameters

* `k`       - Approximate number of clusters
* `m`       - Hyperparameter of SLIC model
* `tol`     - Tolerance of k-means algorithm (default to `1e-4`)
* `maxiter` - Maximum number of iterations (default to `10`)
* `vars`    - Variables (or features) to consider (default to all)

## References

* Achanta et al. 2011. [SLIC superpixels compared to state-of-the-art
  superpixel methods](https://ieeexplore.ieee.org/document/6205760)
"""
struct SLIC <: PartitionMethod
  k::Int
  m::Float64
  tol::Float64
  maxiter::Int
  vars::Union{Vector{Symbol},Nothing}
end

SLIC(k::Int, m::Real; tol=1e-4, maxiter=10, vars=nothing) =
  SLIC(k, m, tol, maxiter, vars)

function partition(data, method::SLIC)
  # variables used for clustering
  dvars = Tables.schema(values(data)).names
  vars = isnothing(method.vars) ? dvars : method.vars

  @assert vars ⊆ dvars "SLIC features not found in geospatial data"

  # view subset of variables
  ctor = constructor(typeof(data))
  dom  = domain(data)
  tab  = TableOperations.select(values(data), vars...)
  Ω    = ctor(dom, Dict(paramdim(dom) => tab))

  # SLIC hyperparameter
  m = method.m

  # initial spacing of clusters
  s = slic_spacing(Ω, method)

  # initialize cluster centers
  c = slic_initialization(Ω, s)

  # ball neighborhood search
  searcher = BallSearch(Ω, NormBall(s))

  # pre-allocate memory for label and distance
  l = fill(0, nelements(Ω))
  d = fill(Inf, nelements(Ω))

  # performance parameters
  tol     = method.tol
  maxiter = method.maxiter

  # k-means algorithm
  err, iter = Inf, 0
  while err > tol && iter < maxiter
    o = copy(c)

    slic_assignment!(Ω, searcher, m, s, c, l, d)
    slic_update!(Ω, c, l)

    err = norm(c - o) / norm(o)
    iter += 1
  end

  subsets = [findall(isequal(k), l) for k in 1:length(c)]

  Partition(data, subsets)
end

function slic_spacing(Ω, method)
  V = measure(boundingbox(Ω))
  d = embeddim(Ω)
  k = method.k
  (V/k) ^ (1/d)
end

function slic_initialization(Ω, s)
  # efficient neighbor search
  searcher = KNearestSearch(Ω, 1)

  # bounding box properties
  bbox = boundingbox(Ω)
  lo, up = coordinates.(extrema(bbox))

  # cluster centers
  clusters = Vector{Int}()
  neighbor = Vector{Int}(undef, 1)
  ranges = [(l+s/2):s:u for (l, u) in zip(lo, up)]
  for x in Iterators.product(ranges...)
    search!(neighbor, Point(x), searcher)
    push!(clusters, neighbor[1])
  end

  unique(clusters)
end

function slic_assignment!(Ω, searcher, m, s, c, l, d)
  for (k, cₖ) in enumerate(c)
    pₖ = centroid(Ω, cₖ)
    inds = search(pₖ, searcher)

    # distance between points
    X  = (coordinates(centroid(Ω, ind)) for ind in inds)
    xₖ = [coordinates(pₖ)]
    dₛ = pairwise(Euclidean(), X, xₖ)

    # distance between variables
    𝒮ᵢ = view(Ω, inds)
    𝒮ₖ = view(Ω, [cₖ])
    V  = Tables.matrix(values(𝒮ᵢ))
    vₖ = Tables.matrix(values(𝒮ₖ))
    dᵥ = pairwise(Euclidean(), V, vₖ, dims=1)

    # total distance
    dₜ = @. √(dᵥ^2 + m^2 * (dₛ/s)^2)

    @inbounds for (i, ind) in enumerate(inds)
      if dₜ[i] < d[ind]
        d[ind] = dₜ[i]
        l[ind] = k
      end
    end
  end
end

function slic_update!(Ω, c, l)
  for k in 1:length(c)
    inds = findall(isequal(k), l)
    X  = (coordinates(centroid(Ω, ind)) for ind in inds)
    μ  = [mean(X)]
    dₛ = pairwise(Euclidean(), X, μ)
    @inbounds c[k] = inds[argmin(vec(dₛ))]
  end
end


text:
@inline compute_prime(x::UInt8, i) = 30 * (i - 1) + ps[trailing_zeros(x) + 0x01]

"""
An iterator for finding small (<= 1_000_000) primes.

    for p in SmallSieve(1_000_000)
        println(p)
    end

Uses n ÷ 30 bytes of memory. Skips 2, 3, and 5; starts at 7.
"""
struct SmallSieve
    xs::Vector{UInt8}

    function SmallSieve(n::Integer)
        # Unrolled loop without segments
        n_bytes = cld(n, 30)
        xs = fill(0xFF, n_bytes)

        # Ensure `1` is not a prime number
        @inbounds xs[1] &= wheel_mask(1)

        # And ensure numbers > n are not prime since we are not considering them
        @inbounds for i = 8 : -1 : 1
            n >= 30 * (n_bytes - 1) + ps[i] && break
            xs[n_bytes] &= wheel_mask(ps[i])
        end

        hi = isqrt(n)

        @inbounds for i = eachindex(xs)
            x = xs[i]
            while x != 0x00
                # The next prime number
                p = compute_prime(x, i)

                # Are we done yet?
                p > hi && @goto done

                # Otherwise cross off multiples of p starting at p²
                p²        = p * p
                byte_idx  = p² ÷ 30 + 1
                wheel     = to_idx(p % 30)
                wheel_idx = 8 * (wheel - 1) + wheel
                increment = i - 1
                @sieve_loop :unroll # Just unroll -- no segmented business here
                x &= x - 0x01
            end
        end

        @label done

        return new(xs)
    end
end

# Yes, this is an O(n) computation, but it should be much faster than computing
# the prime numbers from the bitmask anways
Base.length(s::SmallSieve) = vec_count_ones(s.xs)
Base.eltype(s::SmallSieve) = Int

@inline function Base.iterate(s::SmallSieve)
    @inbounds for i = eachindex(s.xs)
        x = s.xs[i]
        x !== 0x00 && return compute_prime(x, i), (x & (x - 0x01), i)
    end

    return nothing
end

@inline function Base.iterate(s::SmallSieve, state::Tuple{UInt8, Int})
    x, i = state
    @inbounds while true
        x !== 0x00 && return compute_prime(x, i), (x & (x - 0x01), i)
        i === length(s.xs) && return nothing
        x = s.xs[i += 1]
    end
end

function small_primes(n)
    # This runs 20% faster than collect(SmallSieve(n)) :(
    xs = SmallSieve(n).xs
    primes = Vector{Int}(undef, vec_count_ones(xs))
    j = 0
    @inbounds for i = eachindex(xs)
        x = xs[i]
        while x != 0x00
            primes[j += 1] = compute_prime(x, i)
            x &= x - 0x01
        end
    end

    primes
end


text:
function normout!(A::SparseMatrixCSC{Float64,Int64})
    d = sum(A,2) # sum over rows
    # use some internal julia magic
    for i=1:length(A.nzval)
        A.nzval[i] = A.nzval[i] / d[A.rowval[i]]
    end
    return A
end

export normout!

##

""" Filter out the entries that are not in the topk list for each column. """
function filter_topk!(A::SparseMatrixCSC, k::Integer)

  k < 0 && throw(ArgumentError("The value of k = $k must be non-negative"))
  rows = rowvals(A)
  vals = nonzeros(A)
  m, n = size(A)
  sortarr = Array{eltype(A)}(m)
  permarr = Array{Int}(m)

  for i = 1:n
    currange = nzrange(A,i)
    if length(currange) <= k
      continue
    end
    first = currange[1]
    last = currange[end]
    len = last-first+1

    resize!(sortarr, len)
    resize!(permarr, len)
    copy!(sortarr, 1, vals, first, len)
    sortperm!(permarr, sortarr, order=Base.Order.Reverse)
    for j in nzrange(A, i) # clear out the results
      vals[j] = 0
    end
    # at this point, sortarr is still untouched, so
    # sortarr[j] -> currange[j]
    for (j, ji) in enumerate(permarr)
      if j > k # last one! we are done
        break
      end
      #sortarr[ji] is an entry we want to return!
      vals[currange[ji]] = sortarr[ji]
    end
  end
  dropzeros!(A)
end
filter_topk(A,k) = dropzeros!(filter_topk!(copy(A),k))

export filter_topk, filter_topk!


text:
module TestSyncBarriers

using Test
using SyncBarriers
using SyncBarriers: Internal

function test_barrier()
    barrier_factories = [
        SyncBarriers.CentralizedBarrier,
        SyncBarriers.DisseminationBarrier,
        SyncBarriers.TreeBarrier{2},
        SyncBarriers.TreeBarrier{3},
        SyncBarriers.FlatTreeBarrier{2},
        SyncBarriers.FlatTreeBarrier{3},
        SyncBarriers.StaticTreeBarrier{4,2},
        SyncBarriers.StaticTreeBarrier{5,3},
        # ...
    ]
    ntasks_list = unique([1, 2, Threads.nthreads(), 2 * Threads.nthreads()])
    @testset for f in barrier_factories, ntasks in ntasks_list
        test_barrier(f, ntasks)
    end
end

function test_barrier(factory, ntasks)
    @debug "`test_barrier($factory, $ntasks)`"
    states = use_barrier(factory, ntasks)
    ncycles = size(states, 1)
    desired = repeat((1:ncycles) .* ntasks, 1, ntasks)
    @test states[:, 1] == (1:ncycles) .* ntasks
    @test all(states[:, 1] .== states)
    @test states == desired
end

function use_barrier(factory, ntasks; kwargs...)
    barrier = factory(ntasks)
    if factory isa Type
        @test barrier isa factory
    end
    return use_barrier(barrier; kwargs...)
end

function use_barrier(barrier::SyncBarriers.Barrier; ncycles = 1000)
    ntasks = length(barrier)
    states = zeros(Int, ncycles, ntasks)
    return use_barrier!(states, barrier)
end

function use_barrier!(states::Matrix, barrier::SyncBarriers.Barrier, spin = nothing)
    ncycles, ntasks = size(states)
    @assert ntasks == length(barrier)
    value = Threads.Atomic{Int}(0)
    @sync for i in 1:ntasks
        Threads.@spawn try
            for k in 1:ncycles
                Threads.atomic_add!(value, 1)
                SyncBarriers.cycle!(barrier[i], spin)
                states[k, i] = value[]
                SyncBarriers.cycle!(barrier[i], spin)
            end
        catch err
            @error(
                "`use_barrier` failed",
                exception = (err, catch_backtrace()),
                i,
                current_task()
            )
            # TODO: close(barrier)
            rethrow()
        end
    end
    return states
end

function test_default_barriers()
    @test SyncBarriers.Barrier(2) isa SyncBarriers.CentralizedBarrier
    @test SyncBarriers.fuzzy_barrier(2) isa SyncBarriers.CentralizedBarrier
    @test SyncBarriers.fuzzy_reduce_barrier(+, Int, 2) isa SyncBarriers.FlatTreeBarrier{2,Int}
    @test SyncBarriers.reduce_barrier(+, Int, 2) isa SyncBarriers.FlatTreeBarrier{2,Int}

    @test SyncBarriers.Barrier(32) isa SyncBarriers.Barrier
    @test SyncBarriers.fuzzy_barrier(32) isa SyncBarriers.Barrier
    @test SyncBarriers.fuzzy_reduce_barrier(+, Int, 32) isa SyncBarriers.Barrier
    @test SyncBarriers.reduce_barrier(+, Int, 32) isa SyncBarriers.Barrier
end

function test_default_barriers_internal()
    @testset for n in 1:32, nthreads in 1:32
        @test Internal._Barrier(n, nthreads) isa SyncBarriers.Barrier
        @test Internal._fuzzy_barrier(n, nthreads) isa SyncBarriers.Barrier
        @test Internal._fuzzy_reduce_barrier(+, Int, n, nthreads) isa SyncBarriers.Barrier
        @test Internal._reduce_barrier(+, Int, n, nthreads) isa SyncBarriers.Barrier
    end
end

end  # module


text:
using Plots
using MOSLab

NpolyDoping = 5e17 # Gate Npoly doping concentration in cm^-3
N_a = 5e17 # Bulk doping concentration in cm^-3
E_a = 0.044 # The ground state energy of the acceptor in eV
t_ox = 4e-7 # oxide thickness in cm
t_si = 1e-6
MOSf(T) = MOSStructure(NPoly(NpolyDoping),SiO2(),SemiconductorData(T,BoltzmanDist(),PSilicon(N_a,E_a)),t_ox,t_si) ## Calculate Parameters of a MOS Structure the given parameters using Boltzman Distribution at temperature T
g = 1
d = 2
s = 0
M1 = CircuitComponent("M1",ACMModel(MOSf(300.0),1000.0e-4,1.0e-4))
netT = Netlist()
addComponent(netT,M1,Dict("g"=>1,"d"=>d,"s"=>s))
addComponent(netT,CircuitComponent("Vg",VoltageSource(0.3)),Dict("p"=>g,"n"=>s))
addComponent(netT,CircuitComponent("Vd",VoltageSource(1.0)),Dict("p"=>d,"n"=>s))
ckt = Circuit(netT)
res = dc_op(ckt)
for i in 2:length(netT.nodes) # exclude ground
    println("$(netT.nodes[i].name): $(res[i-1])")
end


text:
#Definitions of the different reference frames used and the relations between them.
#= Frames included by default are;
    - International celastial reference Frame (ICRF)
    - CBCI (CelestialBody Centred inertial/J2000 frame)
    - CBCF (Celestial body centred and fixed)
    - Vertical (vertical/local horizon frame)
    - Trajectory (Trajectory Frame)
    - Aerodynamic (Aerodynamic frame )
    - Body (Body frame)
Simulations will be performed in the ICRF, because this is the pseudo-inertial
reference frame that is considered are as inertial.

+ Preferably you want structs for all frames, but how then to refer them in a
tree-like structure?  Make the frames abstract struct will do the trick!
Supertypes will then determine which transformations needs to be explicetely written!
=#

#Define Supertype of the reference frames
abstract type AbstractFrame end

function NewFrame(name::String)
    FrameName = Symbol(name)
    @eval begin
        struct $FrameName <: AbstractFrame end
    end
end

#Global frames that do not require the position of the vehicle
NewFrame("ICRF")
NewFrame("ECI")
NewFrame("ECEF")

#Frames that requre the state of the vehicle
NewFrame("Vertical")
NewFrame("Trajectory")
NewFrame("Aerodynamic")
#NewFrame("Body") -> should be the vehicle type

#create transformation functions for pre-defined frames
struct TransformationQuaternion
    f1::Type{<:AbstractFrame}
    f2::Type{<:AbstractFrame}
    rel::Function
end

#=
#create transformation quaternion
function q_r2i(theta::Float64)
    return angleaxis2quat(theta, [0. 0. 1.])
end
Q_R2I = TransformationQuaternion(ECEF, ECI, q_r2i)

function q_i2r(theta::Float64)
    return angleaxis2quat(-theta, [0. 0. 1.])
end
Q_I2R = TransformationQuaternion(ECI, ECEF, q_i2r)

function q_r2v()

end

function q_aa2i(y::SimState)

end
Q_AA2I = TransformationQuaternion(Aerodynamic, ECI, q_aa2i)
=#


text:
module ChaosGame

    src = [
    ]

    for s = src
        include("$(s).jl")
    end

end # module


text:
## TEST of HBAI inequality - comment this out
using ReadStat
using StatFiles
using DataFrames

hbai_2003_4 = DataFrame( load( "/mnt/data/hbai/UKDA-5828-stata11/stata11/hbai0304_g4.dta" ))
iter = IteratorInterfaceExtensions.getiterator( hbai_2003_4 )
fieldnames( eltype( iter ))

povline = hbai_2003_4[1,:MDOEAHC]*0.6 # The median of S_OE_AHC see if we can replicate this
growth = 0.023

pov_2003_4 = makepoverty( hbai_2003_4, povline, growth, :G_NEWPP, :S_OE_AHC )

print( "Poverty 2003/4\n")
print( pov_2003_4 )
print( "\n\n")

ineq_2003_4 = makeinequality( hbai_2003_4, :G_NEWPP, :S_OE_AHC )

print( "Inequality 2003/4\n")
print( ineq_2003_4 )
print( "\n\n")

print( "deciles\n")
deciles = binify( hbai_2003_4, 10, :G_NEWPP,  :S_OE_AHC )
for i in 1:10
    print( "$i = " );print( deciles[i,:]);print("\n")
end

print( "vigintiles\n")
vigintiles = binify( hbai_2003_4, 20, :G_NEWPP,  :S_OE_AHC )
for i in 1:20
    print( "$i = " );print( vigintiles, [i,:]);print("\n")
end

print( "percentiles\n")
percentiles = binify( hbai_2003_4, 100, :G_NEWPP,  :S_OE_AHC )
for i in 1:100
    print( "$i = " );print( percentiles[i,:]);print("\n")
end


text:
export FdParams,
       get_helmholtz_LU,
       Source,
       get_wavefield_FDFD,
       Recordings,
       get_recordings!,
       amplitude_spectra,
       MonochromaticRecordings,
       sample_spt2rec!,
       inject_rec2spt!,
       get_residue,
       velocity_gradient

include("helmholtz.jl")
include("spectra.jl")
include("recordings.jl")
include("gradient.jl")


text:
module TestMetacommunity
using Test
using LinearAlgebra

using Diversity
using Diversity.API
using DataFrames
using Missings

three = [0.3
         0.3
         0.4]
oc_count = Onecommunity()
ab3 = [1 3 0
       2 0 4]
abnorm = ab3 / sum(ab3)
sc = Subcommunities(size(ab3, 2))
@testset "Communities" begin
    oc_2 = Onecommunity("All of it")
    @test countsubcommunities(Onecommunity("Hello")) == countsubcommunities(oc_2)
    @test getsubcommunitynames(oc_2) == ["All of it"]
    @test getsubcommunitynames(oc_count) == ["1"]
    @test getsubcommunitynames(sc) == map(x -> "$x", 1:countsubcommunities(sc))
    @test getsubcommunitynames(Subcommunities(["a", "b"])) == ["a", "b"]
end

sim = [1.0 0.0 0.0
       1.0 1.0 0.0
       1.0 1.0 1.0]
ms = GeneralTypes(sim)
@testset "GeneralTypes" begin
    @test calcsimilarity(ms, 1.0) == sim
    @test_throws DomainError GeneralTypes(-sim)
    @test_throws MethodError GeneralTypes(ab3)
    @test_throws DimensionMismatch GeneralTypes(convert(Matrix{Float64}, ab3))
end

tax = Taxonomy(DataFrame(Species=["This", "That"]), Dict(:Species=>1.0))
@testset "Taxonomy" begin
    @test counttypes(tax) == 2
    @test gettypenames(tax) == ["This", "That"]
    @test_throws ErrorException calcsimilarity(tax, 1.0)
    @test_throws ErrorException _calcordinariness(tax, abnorm, 1.0)
    @test_throws ErrorException Taxonomy(DataFrame(Species=["This", "That"]),
                                         Dict(:Species=>1.0, :Genus=>0.5))
    @test_throws ErrorException Taxonomy(DataFrame(Genus=["This", "That"]),
                                         Dict(:Genus=>1.0))
    @test gettypenames(Taxonomy(DataFrame(Genus=["This", "That"]),
                                Dict(:Genus=>1.0), :Genus)) == ["This", "That"]
end

@testset "Type names" begin
    @test gettypenames(Species(3)) == map(x -> "$x", 1:3)
    @test gettypenames(Species(["My species"])) == ["My species"]
    @test gettypenames(GeneralTypes(Diagonal([1.0]), ["My species"])) == ["My species"]
    @test gettypenames(UniqueTypes(["One", "Two"])) == ["One", "Two"]
end

meta = Metacommunity(three, ms, oc_count)
sp = Species(size(ab3, 1))
abf = ab3 ./ sum(ab3)
meta2 = Metacommunity(abf, sp, sc)
g2 = GeneralTypes(Matrix(1.0I, 2, 2))
@testset "Metacommunity" begin
    @test gettypes(meta) == ms
    @test getpartition(meta) == oc_count
    @test ismissing(meta.ordinariness)
    @test getordinariness!(meta) ≈ [0.3, 0.6, 1.0]
    @test !ismissing(meta.ordinariness)
    @test getabundance(Metacommunity(ab3, g2, sc)) ≈
        getabundance(Metacommunity(abf, g2, sc))
    @test_nowarn Metacommunity([0.5, 0.5], Matrix(1.0I, 2, 2))
    @test_throws ErrorException Metacommunity(abf, ms, sc)
    @test_throws DimensionMismatch Metacommunity([1, 2, 3]/6, meta2)
    @test_throws DimensionMismatch getabundance(Metacommunity([0.5, 0.5], meta2))
    @test getabundance(Metacommunity(abf, Matrix(1.0I, 2, 2))) ≈
        getabundance(Metacommunity(abf, meta2))
    #@test_throws ErrorException Metacommunity(-abf, g2, sc)
    @test calcsimilarity(gettypes(meta2), _getscale(meta2)) ≈ Matrix(1.0I, size(ab3, 1), size(ab3, 1))
    @test Diversity.floattypes(meta) ⊆ mapreduce(Diversity.floattypes, ∩,
                                                 [getabundance(meta),
                                                  getpartition(meta),
                                                  gettypes(meta)])
end

end