text:
using PenaltyFunctions, LaTeXStrings, Plots
pyplot()

_mkpath(fname) = joinpath(dirname(@__FILE__()), fname)
srand(123)

# ---------------------------------------------------------------
# README Plots

fnt = Plots.Font("ubuntu",8,:hcenter,:vcenter,0.0,
                 RGB{FixedPointNumbers.Normed{UInt8, 8}}(0.0,0.0,0.0))
default(bg_outside = RGBA(1,1,1,0),
        grid = false,
        legendfont = fnt,
        palette = :darktest,
        aspect_ratio = 1,
        size = (400, 400))

p = plot(L1Penalty(); ylim = (0, 8))
for pen in [L2Penalty(),
            ElasticNetPenalty(.5),
            MCPPenalty(1.), SCADPenalty(3.7, 0.4),
            LogPenalty(1.)]
    plot!(p, pen)
end
plot!(p, xlab = L"\theta_1", ylab = L"g(\theta_1)", xlim = (-4, 4))

p[1][1][:linestyle] = :solid
p[1][2][:linestyle] = :solid
p[1][3][:linestyle] = :dash
p[1][4][:linestyle] = :dashdot
p[1][5][:linestyle] = :dash
p[1][6][:linestyle] = :dot

Plots.svg(_mkpath("univariate"))

function cont_xy(penalty, stepx = 0.005, stepy = 0.0001)
    x = -1:stepx:1
    yp = 0:stepy:1
    ref = round(value(penalty, [1., 0]), 4)
    xs = collect(i for j in yp, i in x if round(value(penalty, [i,j]), 4) ≈ ref)
    ys = collect(j for j in yp, i in x if round(value(penalty, [i,j]), 4) ≈ ref)
    vcat(xs, reverse(xs)), vcat(ys, -ys)
end

default(size = (400, 400))
p = plot(cont_xy(L1Penalty())...)
plot!(cont_xy(L2Penalty())...)
plot!(cont_xy(ElasticNetPenalty(0.5))...)
plot!(cont_xy(MCPPenalty(1.))...)
plot!(cont_xy(SCADPenalty(3.7, 0.4))...)
plot!(cont_xy(LogPenalty(1.))...)
plot!(p, xlab = L"\theta_1", ylab = L"\theta_2",
         legend = false,
         xlim = (-1.2, 1.2), ylim = (-1.2, 1.2))

p[1][1][:linestyle] = :solid
p[1][2][:linestyle] = :solid
p[1][3][:linestyle] = :dash
p[1][4][:linestyle] = :dashdot
p[1][5][:linestyle] = :dash
p[1][6][:linestyle] = :dot

Plots.svg(_mkpath("bivariate"))

max_stars_repo_name: JuliaTagBot/FileStorage

max_issues_repo_name: JuliaTagBot/FileStorage

max_forks_repo_name: JuliaTagBot/FileStorage


text:
#---------------------------------------------------------------------#
#This function computes the RK Coefficients.
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
include("create_rhs.jl")

function ti_SSP!(q0,u,Dhat,Fhat,intma,periodicity,time,ntime,dt,stages,DFloat)

    #Initialize RK coefficients
    a0=zeros(DFloat,stages)
    a1=zeros(DFloat,stages)
    beta=zeros(DFloat,stages)
    if stages == 1 #RK1-SSP
        a0[1]=1; a1[1]=0; beta[1]=1
    elseif stages == 2 #RK2-SSP
        a0[1]=1; a1[1]=0; beta[1]=1; #SSP
        a0[2]=1.0/2.0; a1[2]=1.0/2.0; beta[2]=1.0/2.0
    elseif stages == 3 #RK3-SSP
        a0[1]=1; a1[1]=0; beta[1]=1;
        a0[2]=3.0/4.0; a1[2]=1.0/4.0; beta[2]=1.0/4.0
        a0[3]=1.0/3.0; a1[3]=2.0/3.0; beta[3]=2.0/3.0
    end

    #Initialize arrays
    q1=copy(q0); qp=copy(q0)
    Npoin=length(q0)

    #Time Integration
    for itime=1:ntime
        time=time + dt
        #RK Stages
        for s=1:stages
            #Create RHS Matrix
#            R=-Dhat*qp
            R=create_rhs(qp,u,Dhat,Fhat,intma,DFloat)

            #Solve System
            @inbounds for I=1:Npoin
                qp[I]=a0[s]*q0[I] + a1[s]*q1[I] + dt*beta[s]*R[I]
            end #s
            if periodicity[Npoin] == periodicity[1]
                qp[Npoin]=qp[1] #periodicity
            end
            #Update
            q1 .= qp
        end #ik
#        println(" itime = ",itime," time = ",time," qmax = ",maximum(qp)," qmin = ",minimum(qp))

        #Update Q
        q0 .= qp
    end #itime

    return (q0)
end

ext: jl


text:
export Read, @tag_str, calc_distance, calc_ref_pos, calc_read_pos, map_to_read, map_to_ref, phred_to_prob, prob_to_phred

const seqcode = b"=ACMGRSVTWYHKDBN"

macro tag_str(x)
    reinterpret(u16, x.data)[1]
end

#===
NOTE: about pos of indels:
      relpos of insertion: first base of the insertion
      refpos of insertion: the base before the insertion
      relpos of deletion:  the base before the deletion
      refpos of deletion:  fitst base of the deletion
all positions are 1-based
===#

type Read
    refID::i32
    pos::i32
    mapq::Byte
    flag::u16
    next_refID::i32
    next_pos::i32
    tlen::i32
    qname::String
    cigar::Vector{u32}
    seq::Bytes
    qual::Bytes
    tags::Dict{u16, Any}
    muts::Vector{Mut}
    mate::Read

    function Read(f::IO)
        block_size = f >> i32
        refID      = f >> i32
        pos        = f >> i32
        l_qname    = f >> Byte
        mapq       = f >> Byte
        n_cigar_op = f >>> u16 >> u16
        flag       = f >> u16
        l_seq      = f >> i32
        next_refID = f >> i32
        next_pos   = f >> i32
        tlen       = f >> i32
        qname      = f >> l_qname |> del_end! |> String

        cigar = read(f, u32, n_cigar_op)

        seq = Bytes(l_seq)
        for i in 1:l_seq÷2
            c = f >> Byte
            seq[2i-1] = seqcode[c>>4+1]
            seq[2i] = seqcode[c&0x0f+1]
        end
        if isodd(l_seq)
            seq[l_seq] = seqcode[f>>Byte>>4+1]
        end

        qual = f >> l_seq
        tags = f >> (block_size - 32 - l_qname - 4*n_cigar_op - (l_seq+1)÷2 - l_seq) |> parse_tags

        muts = if n_cigar_op != 0
            if haskey(tags, tag"MD")
                try
                    reconstruct_mut_by_md(cigar, tags[tag"MD"], seq)
                catch
                    Mut[]
                end
            else
                println(STDERR, "no `MD` found, try with --reference")
                Mut[]
            end
        else
            Mut[]
        end

        new(refID, pos+1, mapq, flag, next_refID, next_pos+1, tlen, qname, cigar, seq, qual, tags, muts)
    end
end

function hastag(r::Read, tag::u16)
    tag in r.tags
end

function getindex(r::Read, tag::u16)
    get(r.tags, tag, nothing)
end

function setindex!(r::Read, value, tag)
    r.tags[reinterpret(u16, string(tag).data)[1]] = value
end

function parse_tags(x::Bytes)
    f    = IOBuffer(x)
    tags = Dict{u16, Any}()

    while !eof(f)
        tag = f >> u16
        c   = f >> Byte
        value = c == Byte('Z') ? readuntil(f, '\0') |> del_end! :
                c == Byte('i') ? f >> Int32 :
                c == Byte('I') ? f >> UInt32 :
                c == Byte('c') ? f >> Int8 :
                c == Byte('C') ? f >> UInt8 :
                c == Byte('s') ? f >> Int16 :
                c == Byte('S') ? f >> UInt16 :
                c == Byte('f') ? f >> Float32 :
                c == Byte('A') ? f >> Byte |> Char :
                c == Byte('H') ? error("TODO") :
                c == Byte('B') ? error("TODO") :
                error("unknown tag type $(Char(c))")
        tags[tag] = value
    end

    tags
end

function show(io::IO, r::Read)
    io << r.qname << '\n'

    @printf(io, "ChrID: %-2d  Pos(1-based): %-9d  MapQ(0-60): %-d\n", r.refID,      r.pos,      r.mapq)
    @printf(io, "RNext: %-2d  PNext       : %-9d  TempLength: %-d\n", r.next_refID, r.next_pos, r.tlen)

    io << "Cigar: "
    isempty(r.cigar) ? io << '*' : for i in r.cigar
        io << (i >> 4) << cigarcode[i&0x0f+1]
    end
    @printf(io, "  Flag: %d (", r.flag)
    showflag(io, r.flag)
    io << ")\n"

    io << r.seq << '\n'
    io << (r.qual[1] == 0xff ? '*' : map(x->x+0x21, r.qual)) << '\n'

    for (k,v) in r.tags
        write(io, k)
        io << ':' << tagtype(v) << ':' << v << "  "
    end

    io << '\n'
    for i in r.muts
        io << i << ' '
    end

    io << '\n'
end

function showflag(io::IO, flag::u16)
    is_first = true
    interpunct() = is_first ? (is_first=false; "") : " · "
    flag & 0x0001 != 0 && io << interpunct() << "pair_seq"
    flag & 0x0002 != 0 && io << interpunct() << "aligned"
    flag & 0x0004 != 0 && io << interpunct() << "unmapped"
    flag & 0x0008 != 0 && io << interpunct() << "mate_unmapped"
    flag & 0x0010 != 0 && io << interpunct() << "reverse"
    flag & 0x0020 != 0 && io << interpunct() << "mate_reverse"
    flag & 0x0040 != 0 && io << interpunct() << "r1"
    flag & 0x0080 != 0 && io << interpunct() << "r2"
    flag & 0x0100 != 0 && io << interpunct() << "secondary"
    flag & 0x0200 != 0 && io << interpunct() << "not_pass_filter"
    flag & 0x0400 != 0 && io << interpunct() << "duplicate"
    flag & 0x0800 != 0 && io << interpunct() << "supplementary"
    io
end

"distance between the first and last base in ref"
function calc_distance(r::Read)
    reduce(0, r.cigar) do len, cigar
        ifelse(cigar&0b1101 == 0, len + cigar>>4, len)
    end
end

# return (ref_pos, cigar_op)
function calc_ref_pos(r::Read, relpos::Integer)
    refpos = r.pos - 1
    for cigar in r.cigar
        λ"""
        switch((& cigar 0x0f)
               *0 .((min (>> cigar 4) relpos)
                    ?((== relpos .)
                      return('((+ refpos .) 0x00))
                      >(=(relpos (- relpos .))
                        =(refpos (+ refpos .)))))
               *1 .((min (>> cigar 4) relpos)
                    ?((== relpos .)
                      return('(refpos 0x01))
                      =(relpos (- relpos .))))
               *2 =(refpos (+ refpos (>> cigar 4)))
               *4 .((min (>> cigar 4) relpos)
                    ?((== relpos .)
                      return('(refpos 0x04))
                      =(relpos (- relpos .))))
               *5 >()
               (error "TODO: cigar op not supported"))
        """
    end
    i32(0), 0xff
end

function calc_read_pos(r::Read, refpos::Integer)
    relpos, refpos = 0, refpos - r.pos + 1
    for cigar in r.cigar
        λ"""
        switch((& cigar 0x0f)
               *0 .((min (>> cigar 4) refpos)
                    ?((== refpos .)
                      return('((+ relpos .) 0x00))
                      >(=(relpos (+ relpos .))
                        =(refpos (- refpos .)))))
               *1 =(relpos (+ relpos (>> cigar 4)))
               *2 .((min (>> cigar 4) refpos)
                    ?((== refpos .)
                      return('(relpos 0x02))
                      =(refpos (- refpos .))))
               *4 =(relpos (+ relpos (>> cigar 4)))
               *5 >()
               (error "TODO: cigar op not supported"))
        """
    end
    i32(0), 0xff
end

map_to_read(x::SNP, r::Read)       = SNP(car(calc_read_pos(r, x.pos)), x.ref, x.alt)
map_to_read(x::Insertion, r::Read) = Insertion(car(calc_read_pos(r, x.pos))+1, x.bases)
map_to_read(x::Deletion, r::Read)  = Deletion(car(calc_read_pos(r, x.pos)), x.bases)
map_to_ref(x::SNP, r::Read)        = SNP(car(calc_ref_pos(r, x.pos)), x.ref, x.alt)
map_to_ref(x::Insertion, r::Read)  = Insertion(car(calc_ref_pos(r, x.pos)), x.bases)
map_to_ref(x::Deletion, r::Read)   = Deletion(car(calc_ref_pos(r, x.pos))+1, x.bases)

phred_to_prob(x) = 10 ^ (-i64(x) / 10)
prob_to_phred(x) = rount(Byte, -10 * log(10, x))

function del_end!(s::Bytes)
    ccall(:jl_array_del_end, Void, (Any, UInt), s, 1)
    return s
end

function del_end!(s::String)
    ccall(:jl_array_del_end, Void, (Any, UInt), s.data, 1)
    return s
end


text:
using PyPlot

close("all")

## Uncomment for production
PyPlot.rc("font", family="serif")
PyPlot.rc("xtick", labelsize="x-small")
PyPlot.rc("ytick", labelsize="x-small")
PyPlot.rc("figure", figsize=(4,3))
PyPlot.rc("text", usetex=false)

nband = p.nband
Ks = zeros(Int,3,nband)
Ntot = p.N1*p.N2*p.N3
Ns = [p.N1,p.N2,p.N3]

ind = 1
for l1=-L1:L1
    for l2=-L2:L2
        for l3=-L3:L3
            Ks[:,ind] = [l1, l2, l3]
            ind += 1
        end
    end
end

Wfour = zeros(Complex128,nwannier,(2L1+1)*N1,(2L2+1)*N2,(2L3+1)*N3)
H = zeros(Complex128,nwannier,nwannier,N1,N2,N3)

jmax, kmax = ind2sub((N2,N3),indmax(loc))

ijk_to_kind(i,j,k) = k + j*N3 + i*N2*N3 + 1
for i=0:N1-1
    for j=0:N2-1
        for k=0:N3-1
            ijk = [i;j;k]
            kpt = ijk./Ns
            kind = ijk_to_kind(i,j,k)
            eig_k = λ(kpt,Ks)
            perm = (sortperm(eig_k))
            H[:,:,i+1,j+1,k+1] = A[:,:,i+1,j+1,k+1]'*diagm(sort(eig_k))*A[:,:,i+1,j+1,k+1]
            @assert A[:,:,i+1,j+1,k+1]'A[:,:,i+1,j+1,k+1] ≈ eye(nwannier)
            for iwann = 1:nwannier
                for iband = 1:nband
                    Wfour[iwann,
                          i+1+(Ks[1,perm[iband]]+L1)*N1,
                          j+1+(Ks[2,perm[iband]]+L2)*N2,
                          k+1+(Ks[3,perm[iband]]+L3)*N3] = A[iband,iwann,i+1,j+1,k+1]
                    # mark visually the max of loc
                    # if j==jmax && k==kmax
                    #     Wfour[iwann,
                    #           i+1+(Ks[1,perm[iband]]+L1)*N1,
                    #           j+1+(Ks[2,perm[iband]]+L2)*N2,
                    #           k+1+(Ks[3,perm[iband]]+L3)*N3] = 0
                    # end
                end
            end
        end
    end
end

Wreal = zeros(Complex128, nwannier, (2L1+1)*N1,(2L2+1)*N2,(2L3+1)*N3)
# TODO this is a huge hack because I'm too lazy to do it properly. It makes WF be real when they should be
if do_shift == false && corner == [0.,0.,0.] && (N1 == N2 == 1)
    new_Wfour = circshift(Wfour, (0,0,0,N3/2))
    for iwann = 1:nwannier
        Wreal[iwann,:,:,:] = ifftshift(ifft(fftshift(new_Wfour[iwann,:,:,:])))
    end
else
    for iwann = 1:nwannier
        Wreal[iwann,:,:,:] = ifftshift(ifft(fftshift(Wfour[iwann,:,:,:])))
    end
end

# See https://julialang.org/blog/2016/02/iteration for the trick used here for a completely generic (and fast!) multidimensional algorithm
function four_interp(arr, fac=1)
    R = CartesianRange(size(arr))
    Ns = [size(arr)...]
    arr_fine = zeros(Complex128,(2fac-1)*Ns...)
    shift = CartesianIndex((fac-1)*Ns...)
    # function barrier for type inference
    function fillshift!(dest,src,R,shift)
        for I in R
            dest[shift+I] = src[I]
        end
    end
    fillshift!(arr_fine,fftshift(fft(arr)), R, shift)
    return (2fac-1)^length(size(arr))*ifft(ifftshift(arr_fine))
end

function four_interp_ham(H, fac=1)
    return mapslices(arr -> four_interp(arr,fac),H,3:length(size(H)))
end

# 1D
if dim == 1
    kspace = (0:N3-1)/N3 + k_red_to_real([0.,0.,0.])[3]
    ξspace = linspace(-L3+kspace[1],L3+kspace[end], (2L3+1)*N3) #TODO not actually linspace, but not that important
    rspace = linspace(-pi*N3,pi*N3, (2L3+1)*N3)
    # plot the gauge
    figure()
    plot(kspace, sort(p.eigs,2))
    ylim(0,2)
    xlim(kspace[1],kspace[end])
    xlabel("k")
    ylabel("ɛk")
    tight_layout()
    savefig("bandplot.pdf")

    figure()
    markers = ("s","x","o","d")
    markers = ("-x","--s","o","d")
    for n=2:2
    # for n=1:1
        gca()[:set_color_cycle](nothing)
        # for m=1:p.nband
        for m=1:3
            plot(kspace, real.(A[m,n,1,1,:]),markers[n], label="A$m$n",markersize=4, markeredgewidth=2,linewidth=2)
            plot(kspace, imag.(A[m,n,1,1,:]),markers[n], label="A$m$n",markersize=4, markeredgewidth=2,linewidth=2)
            # plot(imag.(A[m,n,1,1,:]),markers[n], label="$m $n",markersize=8, markeredgewidth=2,linewidth=2)
        end
    end
    xlim(kspace[1],kspace[end])
    legend()
    xlabel("k")
    ylabel("Re(Amn)")
    tight_layout()
    savefig("gauge.pdf")

    N = N3
    Wfour = squeeze(Wfour,(2,3))
    Wreal = squeeze(Wreal,(2,3))
    H = squeeze(H,(3,4))

    # normalize here rather than try to figure out correct normalizations
    for n=1:p.nwannier
        Wfour[n,:] /= sqrt(sum(abs2,Wfour[n,:]) * (ξspace[2]-ξspace[1]))
        Wreal[n,:] /= sqrt(sum(abs2,Wreal[n,:]) * (rspace[2]-rspace[1]))
    end

    figure()
    plot(ξspace,real(Wfour)[1,:],"-")
    # plot(ξspace,real(Wfour)[2,:],"--")
    plot(ξspace,imag(Wfour)[1,:],"--")
    # plot(ξspace,imag(Wfour)[2,:],"-.")
    xlim((-3,3))
    legend(vcat(["Re(w$i)" for i=1:1],["Im(w$i)" for i=1:1]))
    xlabel("ξ")
    ylabel("w(ξ)")
    tight_layout()
    savefig("wfour.pdf")
    # plot(ξspace,imag(Wfour))

    figure()
    semilogy(rspace,abs.(Wreal)[1,:])
    # semilogy(rspace,1./rspace.^2,"k")
    # legend(["w1","1/r²"])
    xlim(rspace[1]/2,rspace[end]/2)
    ylim(1e-10,1)
    xlabel("r")
    ylabel("|w(r)|")
    tight_layout()
    savefig("wreal_semilog.pdf")
    # plot(rspace,imag(Wreal))

    figure()
    print(size(Wreal))
    plot(rspace,real(Wreal)[1,:],"-")
    plot(rspace,real(Wreal)[2,:],"--")
    # plot(rspace,real(Wreal)[3,:],"--")
    # plot(rspace,imag(Wreal)')
    legend(["w1","w2"])
    xlabel("r")
    ylabel("w(r)")
    tight_layout()
    savefig("wreal.pdf")
    # plot(rspace,imag(Wreal))

    figure()
    fac = 10
    Nint = (2fac-1)*N
    x = (0:N-1)/N
    xint = (0:Nint-1)/Nint

    if do_shift
        # This is a bit awkward because Fourier interpolation works on a (0:N-1)/N grid
        x += 1/2/N
        xint += 1/2/N # this is really N and not Nint

        ## Here's a test showing why Nint: Fourier interpolation on a staggered grid (at least as defined here) actually extrapolates a bit also
        # N = 10
        # x = ((0:N-1)/N+1/2/N)*2pi
        # f = sin.(x)
        # plot(x,f)
        # Nint = 5N
        # xint = ((0:Nint-1)/Nint+1/2/Nint)*2pi
        # xint = ((0:Nint-1)/Nint+1/2/N)*2pi
        # plot(xint,four_interp(f,3))

    end

    x += corner[3]
    xint += corner[3]

    println(x)
    println(xint)

    Hint = four_interp_ham(H,fac)

    # not in reduced BZ
    λ_real(k,Ks) = squeeze(sum(abs2, (Ks .+ k), 1),1)

    plot(xint,[sort(λ_real([0,0,xint[i]],Ks))[1] for i =1:Nint], "-k", label="Exact")
    plot(xint,[sort(λ_real([0,0,xint[i]],Ks))[2] for i =1:Nint], "-k")
    # plot(xint,[sort(λ_real([0,0,xint[i]],Ks))[3] for i =1:Nint], "-k")
    # ylim(-.1,1.1)
    # plot(xint,[sort(λ([0,0,xint[i]],Ks))[3] for i =1:Nint], "-k")
    plot(x,[sort(real(eig(H[:,:,i])[1]))[m] for i=1:N, m=1], "x", label="Samples")
    plot(x,[sort(real(eig(H[:,:,i])[1]))[m] for i=1:N, m=2], "x")
    # plot(x,[sort(real(eig(H[:,:,i])[1]))[m] for i=1:N, m=3], "o")
    gca()[:set_color_cycle](nothing)
    plot(xint,[sort(real(eig(Hint[:,:,i])[1]))[m] for i=1:Nint, m=1], "-", label="Interpolation")
    plot(xint,[sort(real(eig(Hint[:,:,i])[1]))[m] for i=1:Nint, m=2], "-")
    # plot(xint,[sort(real(eig(Hint[:,:,i])[1]))[m] for i=1:Nint, m=3], "-")
    legend()
    xlabel("k")
    ylabel("ɛ")
    xlim([corner[3],corner[3]+1]) #clip to actual BZ, not more. TODO wrap stuff above corner[3]+1 around
    tight_layout()
    savefig("interp.pdf")

    maxerr = maximum(abs.([sort(real(eig(Hint[:,:,i])[1]))[m] for i=1:Nint, m=1] - [sort(λ_real([0,0,xint[i]],Ks))[1] for i =1:Nint]))
    println("max err $maxerr")
end

if dim == 2
    npt = 100
    kx = linspace(0,1,npt)+k_red_to_real([0.,0.,0.])[2]
    ky = linspace(0,1,npt)+k_red_to_real([0.,0.,0.])[3]
    eigs = zeros(npt,npt,nband)
    for i=1:npt
        for j=1:npt
            eigs[i,j,:] = sort(λ([0,kx[i],ky[j]],Ks))
        end
    end
    # plot_surface(kspace*kspace',kspace*kspace',eigs[:,:,1])

    if( false )
      for n=1:min(p.nband-1,6)
        # plot_surface(kspace,kspace,eigs[:,:,n],cmap=ColorMap("coolwarm"),linewidth=0,rstride=1,cstride=1,antialiased=false)
        # plot_surface(kspace,kspace,eigs[:,:,n+1] - eigs[:,:,n],cmap=ColorMap("coolwarm"),linewidth=0,rstride=1,cstride=1,antialiased=false)
        figure()
        pcolormesh(kx,ky,log10.(eigs[:,:,n+1] - eigs[:,:,n] + 1e-3))
        title("$n -> $(n+1)")
        # pcolormesh(log10.(eigs[:,:,n]))
        colorbar()
      end
    end

    kx = linspace(0,1,N2)+k_red_to_real([0.,0.,0.])[2]
    ky = linspace(0,1,N3)+k_red_to_real([0.,0.,0.])[3]
    figure()
    pcolormesh(kx,ky,loc[1,:,:],cmap=ColorMap("gray_r"))
    xlabel("kx")
    ylabel("ky")
    colorbar()
    savefig("normgrad_$N3.pdf")
    @assert N2 == N3
    N = N2

    Wfour = squeeze(Wfour,(2))
    Wreal = squeeze(Wreal,(2))
    H = squeeze(H,(3))

    ξ = (0:N3-1)/N3 + k_red_to_real([0.,0.,0.])[3]
    ξspacex = linspace(-L2-kx[1],L2+kx[end], (2L2+1)*N3) #TODO not actually linspace, but not that important
    ξspacey = linspace(-L3-ky[1],L3+ky[end], (2L3+1)*N3) #TODO not actually linspace, but not that important

    for iwann=1:p.nwannier
        # plot the WF in Fourier space
        # figure()
        # pcolormesh(ξspacex,ξspacey,real(Wfour[iwann,:,:]))
        # figure()
        # pcolormesh(ξspacex,ξspacey,imag(Wfour[iwann,:,:]))
        figure()
        pl = pcolormesh(ξspacex,ξspacey,abs.(Wfour[iwann,:,:]),cmap=ColorMap("gray_r"),linewidth=0,rasterized=true) #force rasterization to avoid "kilt" effect when saving
        xlim(-2,2)
        ylim(-2,2)
        xlabel(L"\xi_x")
        ylabel(L"\xi_y")
        colorbar()
        tight_layout()
        savefig("2D_abs_$iwann.pdf",dpi=450)

        # # plot the gradient in Fourier space
        # arr = real(Wfour[iwann,:,:])
        # grad = zeros(size(arr))
        # for i=1:size(arr,1)-1
        #     for j=1:size(arr,2)-1
        #         grad[i,j] = sqrt((arr[i+1,j]-arr[i,j])^2 + (arr[i,j+1]-arr[i,j])^2)
        #     end
        # end
        # grad /= (ξspacex[2]-ξspacex[1])
        # figure()
        # pcolormesh(ξspacex,ξspacey, grad,cmap=ColorMap("gray_r"),linewidth=0,rasterized=true)
        # xlim(-2,2)
        # ylim(-2,2)
        # xlabel(L"\xi_x")
        # ylabel(L"\xi_y")
        # colorbar()
        # tight_layout()
        # savefig("2D_grad_$iwann.pdf",dpi=450)
        # println("Max gradient $iwann: $(maximum(grad))")

        # plot a slice
        figure()
        data = abs.(Wreal[iwann,div(end,2),:])
        rspace = linspace(-pi*N3,pi*N3, (2L3+1)*N3)
        semilogy(rspace,data)
        xlim(rspace[1]/2,rspace[end]/2)
        semilogy(rspace,1./rspace.^2/100,"k")
        legend(["w$iwann","1/r²"])
        xlabel("r")
        ylabel("|w(r)|")
        tight_layout()
        savefig("wreal_slice_$iwann.pdf")

        figure()
        rspace = linspace(-pi*N3,pi*N3, (2L3+1)*N3)
        beg = 160
        finish = 200
        pcolormesh(rspace[beg:finish],rspace[beg:finish], real(Wreal[iwann,beg:finish,beg:finish]),linewidth=0,rasterized=true)
        xlim(-10,10)
        ylim(-10,10)
        xlabel("x")
        ylabel("y")
        colorbar()
        tight_layout()
        savefig("wreal_2D_$iwann.pdf",dpi=450)
    end

    fac = 2
    Nint = (2fac-1)*N
    x = (0:N-1)/N
    xint = (0:Nint-1)/Nint

    Hint = four_interp_ham(H,fac)

    @assert do_shift == false #not implemented yet

    STOP

    figure()

    #exact
    plot(xint,[sort(λ([0,xint[i],0],Ks))[m] for i =1:Nint,m=1:1],"-k",label="Exact")
    plot(xint,[sort(λ([0,xint[i],0],Ks))[m] for i =1:Nint,m=1:nwannier],"-k")
    plot(xint+1,[sort(λ([0,0,xint[i]],Ks))[m] for i =1:Nint,m=1:nwannier],"-k")
    plot(xint+2,[sort(λ([0,xint[i],xint[i]],Ks))[m] for i =1:Nint,m=1:nwannier],"-k")

    #original
    plot(x,[sort(real(eig(H[:,:,i,1])[1]))[m] for i=1:N, m=1:1], "o",label="Samples")
    gca()[:set_color_cycle](nothing)
    plot(x,[sort(real(eig(H[:,:,i,1])[1]))[m] for i=1:N, m=1:nwannier], "o")
    gca()[:set_color_cycle](nothing)
    plot(x+1,[sort(real(eig(H[:,:,1,i])[1]))[m] for i=1:N, m=1:nwannier],"o")
    gca()[:set_color_cycle](nothing)
    plot(x+2,[sort(real(eig(H[:,:,i,i])[1]))[m] for i=1:N, m=1:nwannier], "o")
    gca()[:set_color_cycle](nothing)

    #interpolated
    plot(xint,[sort(real(eig(Hint[:,:,i,1])[1]))[m] for i=1:Nint, m=1:1], "-",label="Interpolation")
    gca()[:set_color_cycle](nothing)
    plot(xint,[sort(real(eig(Hint[:,:,i,1])[1]))[m] for i=1:Nint, m=1:nwannier], "-")
    gca()[:set_color_cycle](nothing)
    plot(xint+1,[sort(real(eig(Hint[:,:,1,i])[1]))[m] for i=1:Nint, m=1:nwannier],"-")
    gca()[:set_color_cycle](nothing)
    plot(xint+2,[sort(real(eig(Hint[:,:,i,i])[1]))[m] for i=1:Nint, m=1:nwannier], "-")

    legend()

    xlabel("k")
    ylabel("ɛk")
    savefig("interp_2D.pdf")

end

if dim==3
    # for iwann = 1:1
    #     figure()
    #     pl = pcolormesh(abs.(Wfour[iwann,div(end,4),:,:]),linewidth=0,rasterized=true) #force rasterization to avoid "kilt" effect when saving
    #     figure()
    #     pl = pcolormesh(abs.(Wfour[iwann,div(end,2),:,:]),linewidth=0,rasterized=true) #force rasterization to avoid "kilt" effect when saving
    #     figure()
    #     pl = pcolormesh(abs.(Wfour[iwann,div(3*end,4),:,:]),linewidth=0,rasterized=true) #force rasterization to avoid "kilt" effect when saving
    # end

    # xlim(-2,2)
    # ylim(-2,2)
    # xlabel("ξx")
    # ylabel("ξy")
    # savefig("2D_abs_$iwann.pdf",dpi=450)

    # data = abs((Wfour[1,:,:,:]))
    # using GLVisualize, GLWindow
    # window = glscreen()
    # timesignal = bounce(linspace(Float32(minimum(data)), Float32(maximum(data)),360))
    # vol = visualize(data, :iso, isovalue=timesignal)
    # _view(vol, window)
    # renderloop(window)

    iwann = 1
    arr = real(Wfour[iwann,:,:,:])
    grad = zeros(size(arr))
    for i=1:size(arr,1)-1
        for j=1:size(arr,2)-1
            for k=1:size(arr,3)-1
                grad[i,j,k] = sqrt((arr[i+1,j,k]-arr[i,j,k])^2 + (arr[i,j+1,k]-arr[i,j,k])^2 + (arr[i,j,k+1]-arr[i,j,k])^2)
            end
        end
    end
    data = grad
    using GLVisualize, GLWindow
    window = glscreen()
    timesignal = bounce(linspace(Float32(minimum(data)), Float32(maximum(data)),360))
    vol = visualize(data, :iso, isovalue=timesignal)
    _view(vol, window)
    renderloop(window)
end

max_stars_repo_name: VVitale/wannier

max_issues_repo_name: VVitale/wannier

max_forks_repo_name: VVitale/wannier


text:
"""
    @include_config(path::AbstractString)

Prefix `path` with the JustSayIt application config path and then call `include(path)`. If `path` is an absolut path, then call `include` with `path` unmodified.

!!! note "NOTE: JustSayIt application config"
    The content of the JustSayIt application config folder is not evaluated within JustSayIt. The folder's single purpose is to provide an easily accessible storage for scripts to start JustSayIt and/or for custom command functions: `@include_config` permits to conveniently `include` files from this folder (for details about the Julia built-in `include` type `?include`).
    Your JustSayIt application config path on this system is: `$CONFIG_PREFIX`
"""
macro include_config(args...) checkargs_include(args...); include_config(__module__, args...); end

function include_config(caller::Module, path)
    path = esc(path)
    return quote
        if !isabspath($path)
            Base.include($caller, joinpath(CONFIG_PREFIX, $path))
        else
            Base.include($caller, $path)
        end
    end
end

checkargs_include(args...) = if (length(args) != 1) @ArgumentError("wrong number of arguments.") end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _eci2geo{T<:AbstractFloat}(x::T, y::T, z::T, jd::T)
    Re    = planets["earth"].eqradius*1e-3
    theta = atan2(y, x) # Azimuth.
    gst   = ct2lst(0.0, jd)
    sid_angle = gst*pi/12.0 # Sidereal angle.
    long  = cirrange(rad2deg(theta - sid_angle)) # Longitude.
    r     = hypot(x, y)
    lat   = atan2(z, r) # Latitude.
    alt   = r/cos(lat) - Re # Altitude.
    lat   = rad2deg(lat)
    return lat, long, alt
end

"""
    eci2geo(x, y, z, jd) -> latitude, longitude, altitude

### Purpose ###

Convert Earth-centered inertial coordinates to geographic spherical coordinates.

### Explanation ###

Converts from ECI (Earth-Centered Inertial) (x, y, z) rectangular coordinates to
geographic spherical coordinates (latitude, longitude, altitude).  Julian day is
also needed as input.

ECI coordinates are in km from Earth center at the supplied time (True of Date).
Geographic coordinates assume the Earth is a perfect sphere, with radius equal
to its equatorial radius.

### Arguments ###

* `x`: ECI x coordinate at `jd`, in kilometers.
* `y`: ECI y coordinate at `jd`, in kilometers.
* `z`: ECI z coordinate at `jd`, in kilometers.
* `jd`: Julian days.

The three coordinates can be passed as a 3-tuple `(x, y, z)`.  In addition, `x`,
`y`, `z`, and `jd` can be given as arrays of the same length.

### Output ###

The 3-tuple of geographical coordinate (latitude, longitude, altitude).

* latitude: latitude, in degrees.
* longitude: longitude, in degrees.
* altitude: altitude, in kilometers.

If ECI coordinates are given as arrays, a 3-tuple of arrays of the same length
is returned.

### Example ###

Obtain the geographic direction of the vernal point on 2015-06-30T14:03:12.857,
in geographic coordinates, at altitude 600 km.  Note: equatorial radii of Solar
System planets in meters are stored into `AstroLib.planets` dictionary.

``` julia
x = AstroLib.planets["earth"].eqradius*1e-3 + 600;
lat, long, alt = eci2geo(x, 0, 0, jdcnv("2015-06-30T14:03:12.857"))
# => (0.0,230.87301833205856,600.0)
```

These coordinates can be further transformed into geodetic coordinates using
`geo2geodetic` or into geomagnetic coordinates using `geo2mag`.

### Notes ###

`geo2eci` converts geographic spherical coordinates to Earth-centered inertial
coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
eci2geo(x::Real, y::Real, z::Real, jd::Real) =
    _eci2geo(promote(float(x), float(y), float(z), float(jd))...)

eci2geo(xyz::Tuple{Real, Real, Real}, jd::Real) =
    eci2geo(xyz..., jd)

function eci2geo{X<:Real, Y<:Real, Z<:Real, JD<:Real}(x::AbstractArray{X},
                                                      y::AbstractArray{Y},
                                                      z::AbstractArray{Z},
                                                      jd::AbstractArray{JD})
    @assert length(x) == length(y) == length(z) == length(jd)
    typex = typeof(float(one(X)))
    lat  = similar(x, typex)
    long = similar(x, typex)
    alt  = similar(x, typex)
    for i in eachindex(x)
        lat[i], long[i], alt[i] =
            eci2geo(x[i], y[i], z[i], jd[i])
    end
    return lat, long, alt
end


text:
module KNeighbors
  using ...Statistics
  using Distances

  function KNearestNeighbors(trainingSet, testSet, training_y,K=1)
    # m = Statistics.meanColumn(trainingSet)
    # s = Statistics.stdColumn(trainingSet)
    # println("m and s calculatd!")
    # trainingSet=DataManipulation.DataTransform.zScale(trainingSet,m,s)
    # testSet=DataManipulation.DataTransform.zScale(testSet,m,s)
    # println("Scaled!")

    predicted_y=zeros(Float64, size(testSet)[1],1)
    dot_prod = zeros(Float64,size(testSet)[1], size(trainingSet)[1])
    for i=1:size(testSet)[1]
      testSetRow=testSet[i,:]
      for j=1:size(trainingSet)[1]
        dot_prod[i,j]=euclidean(testSetRow,trainingSet[j,:])
      end
    end
    #println("we got our matrix!")
    for i=1:size(dot_prod)[1]
      #find the K smallest distances and get their labels. Take a majority. predicted_y[i] = majority
      indices = Statistics.getMinIndices(collect(dot_prod[i,:]),K)
      # println(typeof(indices));
      # println(indices)
      # println(typeof(training_y))
      # println(training_y)
      predicted_y[i]=Statistics.getMajority(training_y[indices];)

    end
    return predicted_y
    #return fn(test_y, predicted_y)

  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
# Here, we illustrate a "julian" implementation that leverages
# some of the advantages of GLAbstraction over basic ModernGL usage

using ModernGL, ColorTypes, GeometryTypes, GLAbstraction, GLFW

const GLA = GLAbstraction

# Create the window. This sets all the hints and makes the context current.
window = GLFW.Window(name="Drawing polygons 5", resolution=(800,600))

# We assign the created window as the "current" context in GLAbstraction to which all GL objects are "bound", this is to avoid using GL objects in the wrong context, but actually currently no real checks are made except initially that at least there is a context initialized.
# Think of this as a way of bookkeeping.
GLFW.MakeContextCurrent(window)
GLA.set_context!(window)
GLA.FrameBuffer((1024,1024),(RGBA{Float32}, RGBA{Float32}, GLA.Depth{Float32}), true)
# The vertex shader---note the `vert` in front of """
vertex_shader = GLA.vert"""
#version 150

in vec2 position;
in vec3 color;

out vec3 Color;

void main()
{
    Color = color;
    gl_Position = vec4(position, 0.0, 1.0);
}
"""

# The fragment shader
fragment_shader = GLA.frag"""
# version 150

in vec3 Color;

out vec4 outColor;

void main()
{
    outColor = vec4(Color, 1.0);
}
"""

# First we combine these two shaders into the program that will be used to render
prog = GLA.Program(vertex_shader, fragment_shader)

# Some geometry to render
# First a triangle
triangle_positions = Point{2,Float32}[(-0.5, -0.5),
                                      (0.5, -0.5),
                                      (0.0, 0.5)]
triangle_colors = Vec3f0[(1, 0, 0),                     # top-left
                         (0, 1, 0),                     # top-right
                         (0, 0, 1)]                     # bottom-right

# This one we will render in a simple way, i.e. not using an index buffer
# This geometry now has to be uploaded into OpenGL Buffers and attached to the correct points in a VertexArray, to be used in the above program.
# We first generate the Buffers and their BufferAttachmentInfos, for each of the buffers, from the corresponding names and data
# The GEOMETRY_DIVISOR distinguishes that this is geometry data and not possible uniform data used in instanced vertexarrays
buffers = GLA.generate_buffers(prog, GLA.GEOMETRY_DIVISOR, position = triangle_positions, color = triangle_colors)
triangle_vao = GLA.VertexArray(buffers)
glClearColor(0,0,0,0)
while !GLFW.WindowShouldClose(window)
    glClear(GL_COLOR_BUFFER_BIT)
    GLA.bind(prog)
    GLA.bind(triangle_vao)
    GLA.draw(triangle_vao)
    GLA.unbind(triangle_vao) #optional in this case
    GLA.unbind(prog) #optional in this case
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end
GLFW.SetWindowShouldClose(window, false)

# Now we define another geometry that we will render, a rectangle, this one with an index buffer
# The positions of the vertices in our rectangle
rectangle_positions = Point{2,Float32}[(-0.5,  0.5),     # top-left
                                       ( 0.5,  0.5),     # top-right
                                       ( 0.5, -0.5),     # bottom-right
                                       (-0.5, -0.5)]     # bottom-left

# The colors assigned to each vertex
rectangle_colors = Vec3f0[(1, 0, 0),                     # top-left
                          (0, 1, 0),                     # top-right
                          (0, 0, 1),                     # bottom-right
                          (1, 1, 1)]                     # bottom-left

# Specify how vertices are arranged into faces
# Face{N,T} type specifies a face with N vertices, with index type
# T (you should choose UInt32), and index-offset O. If you're
# specifying faces in terms of julia's 1-based indexing, you should set
# O=0. (If you instead number the vertices starting with 0, set
# O=-1.)
elements = Face{3,UInt32}[(0,1,2),          # the first triangle
                          (2,3,0)]          # the second triangle

buffers = GLA.generate_buffers(prog, GLA.GEOMETRY_DIVISOR, position = rectangle_positions, color = rectangle_colors)

# Now we create a VertexArray from these buffers and use the elements as the indices
rectangle_vao = GLA.VertexArray(buffers, elements)

# Draw the rectangle until we receive a close event
glClearColor(0,0,0,0)
while !GLFW.WindowShouldClose(window)
    glClear(GL_COLOR_BUFFER_BIT)
    GLA.bind(prog)
    GLA.bind(rectangle_vao)
    GLA.draw(rectangle_vao)
    GLA.unbind(rectangle_vao) #optional in this case
    GLA.unbind(prog) #optional in this case
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end
GLFW.SetWindowShouldClose(window, false)

GLFW.DestroyWindow(window)  # needed if you're running this from the REPL


text:
using POMDPs
using POMDPModels
using POMDPModelTools
using StaticArrays
using JLD2, FileIO
using Random
using Logging
using Distributions
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using CMSSPs

inhor_file = ARGS[1]
outhor_file = ARGS[2]
params_fn = ARGS[3]
trials = 5

rng = MersenneTwister(9810)

modal_horizon_policy = load_modal_horizon_policy_localapprox(inhor_file, outhor_file)

# Create parameters and then CMSSP
params = toy2d_parse_params(params_fn)
cmssp = create_toy2d_cmssp(params)

# Set up the modal mdp
# Need to do it for just one mode
modal_mdp = Toy2DModalMDPType(1, params, cmssp.control_actions, 0.8, params.horizon_limit)

for i = 1:trials

    tot_cost = 0.0

    start_horizon = convert(Int64,rand(params.horizon_limit/2:params.horizon_limit))
    target_state = sample_toy2d(rng)
    curr_state = Toy2DContState(0.0,0.0)
    tp_dist = TPDistribution([start_horizon],[1.0])

    @show target_state, start_horizon

    for t = 0:1:start_horizon-1
        a = get_best_intramodal_action(modal_horizon_policy, t, tp_dist, curr_state, target_state)
        if a == nothing
            println("ABORT!")
            break
        end
        (new_state, reward) = generate_sr(modal_mdp, curr_state, a.action, rng)
        curr_state = new_state
        tot_cost += -1.0*reward
    end

    @show CMSSPs.CMSSPDomains.norm(CMSSPs.CMSSPDomains.get_relative_state(modal_mdp,curr_state,target_state))

    if isterminal(modal_mdp, CMSSPs.CMSSPDomains.get_relative_state(modal_mdp,curr_state,target_state))
        @show "Success!"
    else
        @show "Failure!"
    end
end


text:
function build_arg(kind, type, name, description)
    return Dict{String, Any}(
        "name" => "MyArg",
        "defaultValue" =>  nothing,
        "description" =>  nothing,
        "type" => build_type(kind, name; description, type)
    )
end

function build_type(kind, name; description="", type=nothing, ofType=nothing)
    type_dict =  Dict{String, Any}(
        "kind" => kind,
        "name" => name,
        "description" => description,
    )
    !isnothing(type) && push!(type_dict, "type" => type)
    !isnothing(ofType) && push!(type_dict, "ofType" => ofType)
    return type_dict
end


text:
import Pkg
using Pkg
# TEST
# Pkg.add("CSV")
# Pkg.add("DataFrames")
# Pkg.add("Query")
# Pkg.add("Plots")

using CSV
using DataFrames
using Dates
using Query
using Statistics
using Plots

CSV_fileName = readdir("/home/prabu/Research/mintsData/001e06305a6c/2019/07/20/")       # taking all the file names from the folder
num_files = length(CSV_fileName)                                                        # number of files

dfArray = []            # Define an array of DataFrames (all CSV files)
sensorName = []         # Define an array for abstract sensor names from CSV file.

############ read each CSV file and take average
for i in 1:num_files-1
    df = CSV.read("/home/prabu/Research/mintsData/001e06305a6c/2019/07/20/"*CSV_fileName[i], DataFrame)         # Reading ith CSV data file
    df.dateTime =  SubString.(string.(df.dateTime), 1, 19)                                                      # Removing milisecond part of dateTime
    df.dateTime = DateTime.(df.dateTime,"yyyy-mm-dd HH:MM:SS")                                                  # converting to dateTime

    df.dateTime = map((x) -> round(x, Dates.Second(30)), df.dateTime)                                           # rounding dateTime column for 30s
    gdf = groupby(df, :dateTime)                                                                                # making groups by same dateTime
    cgdf = combine(gdf, valuecols(gdf) .=> mean)                                                                # taking average of same group and combine
    push!(dfArray,cgdf)                                                                                # put the DataFrame in the array

    sName =  SubString.(string.(CSV_fileName[i]), 20, 25)                                                       # taking sensor name from CSV file
    push!(sensorName,sName)                                                                                # put sensor name in the array
end

########### change the column name of each DataFrames
for j in 1:length(sensorName)
    colFullName = names(dfArray[j])                                             # taking all the column name of one DataFrame

    for x in 2:length(colFullName)
        num_let = length(colFullName[x])                                        # taking number of characters in one column name
        colName = SubString.(string.(colFullName[x]), 1, num_let-5)             # removing last five characters of each column name
        rename!(dfArray[j],colFullName[x] => sensorName[j]*"_"*colName)         # adding sensor name to the column name
    end

end

########### combine all the DataFrames
new_df = dfArray[1]                                                                             # first DataFrame
for x in 1:length(dfArray)-1
    global new_df = outerjoin(new_df, dfArray[x+1], on = :dateTime, makeunique=true)            # combine all the DataFrames together
end

final_df = sort!(new_df)                                # arranging raws of total DataFrame in ascending order by dateTime

println(final_df)
CSV.write("/home/prabu/dfcom3.csv", final_df)

#plot(final_df.dateTime, final_df.AS7262_temperature)


text:
"""
Runs the MCT step of MergeSplitTours. Assumes that cost_fn(i, j)
returns the surrogate estimate for travel time between depot_sites[i] and depot_sites[j]
"""
function min_connecting_tour(n_depots::Int64, n_sites::Int64,
                             depot_sites::Vector{LOC},
                             cost_fn::F) where {LOC, F <: Function}

    # Construct the cost vector (no inf terms!)
    # cost_fn takes as argument u_idx and v_idx (1...n_depots....n_depots + n_sites)
    cost_vector = Vector{Float64}(undef, 0)
    depot_site_vector_idx = Dict{Tuple{Int64,Int64},Int64}()
    vector_idx_to_depot_site = Dict{Int64,Tuple{Int64,Int64}}()
    out_nbrs = Dict{Int64,Set{Int64}}()
    in_nbrs = Dict{Int64,Set{Int64}}()

    for i = 1:n_depots+n_sites
        for j = 1:n_depots+n_sites

            # Exclude self-edges and site-site edges
            if i != j && (i <= n_depots || j <= n_depots)
                edge_cost = cost_fn(i, j)

                # Also exclude edges that do not have an estimate, i.e. any dispatch constraints
                if edge_cost < Inf
                    # First insert cost and idx
                    push!(cost_vector, edge_cost)
                    depot_site_vector_idx[(i, j)] = length(cost_vector)
                    vector_idx_to_depot_site[length(cost_vector)] = (i, j)

                    if ~(haskey(out_nbrs, i))
                        out_nbrs[i] = Set{Int64}()
                    end
                    push!(out_nbrs[i], j)

                    if ~(haskey(in_nbrs, j))
                        in_nbrs[j] = Set{Int64}()
                    end
                    push!(in_nbrs[j], i)
                end
            end
        end
    end

    # Since x vector will be subset of (l+k) x (l+k), need to maintain a map
    # from (l,k) to index and index to (l,k)
    model = Model(with_optimizer(GLPK.Optimizer))
    n_idxs = length(cost_vector)

    @variable(model, x[1:n_idxs]) # True for all edges
    @constraint(model, x .>= 0)
    @objective(model, Min, cost_vector' * x)

    # Add {0,1} constraint on package-depot / depot-package edges
    depot_site_mask = spzeros(n_depots*n_sites*2, n_idxs)
    for i = 1:n_depots # Depots
        for j = n_depots+1:n_depots+n_sites # Sites
            # Depot -> site
            idx = (i - 1)*n_sites + (j - n_depots)

            # @show (i, j, idx, idx + n_depots*n_sites)

            if haskey(depot_site_vector_idx, (i, j))
                depot_site_mask[idx, depot_site_vector_idx[(i, j)]] = 1.0
            end

            # Site -> Depot
            if haskey(depot_site_vector_idx, (j, i))
                depot_site_mask[idx + n_depots*n_sites, depot_site_vector_idx[(j, i)]] = 1.0
            end
        end
    end

    @constraint(model, depot_site_mask * x .<= ones(n_depots*n_sites*2)) # Depot<->site edges at most once

    # Now add constraints for out-edges and in-edges for sites
    site_out_edge_mask = spzeros(n_sites, n_idxs)
    site_in_edge_mask = spzeros(n_sites, n_idxs)
    for i = n_depots+1:n_depots+n_sites

        if isempty(out_nbrs[i])
            @warn "Out nbrs for site $i is empty!"
            readline()
        end
        for onbr in out_nbrs[i] # Should only be depots
            idx = depot_site_vector_idx[(i, onbr)]
            site_out_edge_mask[i - n_depots, idx] = 1.
        end

        if isempty(in_nbrs[i])
            @warn "In nbrs for site $i is empty!"
            readline()
        end
        for inbr in in_nbrs[i]
            idx = depot_site_vector_idx[(inbr, i)]
            site_in_edge_mask[i - n_depots, idx] = 1.
        end
    end

    @constraint(model, site_out_edge_mask * x .== ones(n_sites))
    @constraint(model, site_in_edge_mask * x .== ones(n_sites))

    # Finally, add constraints for in-flow and out-flow from depots
    depot_out_edge_mask = zeros(n_depots, n_idxs)
    depot_in_edge_mask = zeros(n_depots, n_idxs)
    for i = 1 : n_depots

        for onbr in out_nbrs[i]
            idx = depot_site_vector_idx[(i, onbr)]
            depot_out_edge_mask[i, idx] = 1.
        end

        for inbr in in_nbrs[i]
            idx = depot_site_vector_idx[(inbr, i)]
            depot_in_edge_mask[i, idx] = 1
        end

        # @constraint(model, depot_out_edge_mask[i, :]' * x - depot_in_edge_mask[i, :]' *x == 0)
    end

    @constraint(model, depot_out_edge_mask * x - depot_in_edge_mask * x .== zeros(n_depots))

    # Now we optimize!
    JuMP.optimize!(model)

    x_edges = JuMP.value.(x)

    # Integrality gap 0
    x_edges = convert(Vector{Int64}, x_edges)

    edges = [vector_idx_to_depot_site[i] for (i, val) in enumerate(x_edges) if val > 0] # not just = 1 as may have > 1

    return (edges, x_edges, vector_idx_to_depot_site, depot_site_vector_idx, cost_vector)

end # function

"""
Get the strongly connected components for depots from the result of MCT.
"""
function get_connected_depots(n_depots::Int64, n_sites::Int64, edges::Vector{Tuple{Int64,Int64}})

    adj_mat = zeros(Int64, n_depots+n_sites, n_depots+n_sites)

    for t in edges
        adj_mat[t[1], t[2]] = 1
    end

    g = LightGraphs.SimpleDiGraph(adj_mat)

    cc = LightGraphs.strongly_connected_components(g)

    # cc is an array of arrays
    depot_components = Vector{Vector{Int64}}(undef, 0)
    for c in cc
        if length(c) > 1
            this_comp_depots = Vector{Int64}(undef, 0)
            for idx in c
                if idx <= n_depots
                    push!(this_comp_depots, idx)
                end
            end

            push!(depot_components, this_comp_depots)
        end
    end

    return depot_components
end

"""
Cluster the depots based on the logic in MergeSplitTours. Edits x_edges::Vector{Int64} in place.
"""
function add_merged_depot_edges!(x_edges::Vector{Int64}, depot_comps::Vector{Vector{Int64}},
                                 depot_site_vector_idx::Dict{Tuple{Int64,Int64},Int64}, n_depots::Int64,
                                 n_sites::Int64, cost_vector::Vector{Float64})

    merged_depot_comps = depot_comps

    # Continue until one component
    while length(merged_depot_comps) > 1

        min_to_merge = (0, 0)
        min_depots_in_merged_comps = (0, 0)
        cmin = Inf

        for comp1_idx = 1:length(merged_depot_comps)-1
            for comp2_idx = comp1_idx+1 : length(merged_depot_comps)

                cmin_for_dep_pair = Inf
                depots_to_merge = (0, 0)

                for dep1 in merged_depot_comps[comp1_idx]
                    for dep2 in merged_depot_comps[comp2_idx]

                        dep_cost = cost_vector[depot_site_vector_idx[(dep1, dep2)]] +
                                    cost_vector[depot_site_vector_idx[(dep2, dep1)]]
                        if dep_cost < cmin_for_dep_pair
                            cmin_for_dep_pair = dep_cost
                            depots_to_merge = (dep1, dep2)
                        end
                    end
                end

                if cmin_for_dep_pair < cmin
                    cmin = cmin_for_dep_pair
                    min_to_merge = (comp1_idx, comp2_idx)
                    min_depots_in_merged_comps = depots_to_merge
                end
            end
        end

        # Add edges for merged deps
        dep1to2 = depot_site_vector_idx[(min_depots_in_merged_comps[1], min_depots_in_merged_comps[2])]
        @assert x_edges[dep1to2] == 0.0 "$(min_to_merge) \n $(merged_depot_comps)"
        x_edges[dep1to2] = 1.0
        dep2to1 = depot_site_vector_idx[(min_depots_in_merged_comps[2], min_depots_in_merged_comps[1])]
        @assert x_edges[dep2to1] == 0.0
        x_edges[dep2to1] = 1.0

        # Now create new merged depot comps
        new_comps = length(merged_depot_comps) - 1
        new_merged_depot_comps = Vector{Vector{Int64}}(undef, 0)

        push!(new_merged_depot_comps, vcat(merged_depot_comps[min_to_merge[1]], merged_depot_comps[min_to_merge[2]]))

        for (i, deps) in enumerate(merged_depot_comps)

            if ~(i in min_to_merge)
                push!(new_merged_depot_comps, deps)
            end
        end

        @assert length(new_merged_depot_comps) == new_comps

        merged_depot_comps = deepcopy(new_merged_depot_comps)
    end

    return merged_depot_comps
end

"""
Compute the Eulerian tour which will be distributed over m agents.
"""
function get_multiedge_eulerian_tour(x_edges::Vector{Int64}, vector_idx_to_depot_site::Dict{Int64,Tuple{Int64,Int64}},
                                     n_depots::Int64, n_sites::Int64)

    edge_count = zeros(Int64, n_depots+n_sites)
    adj_list = [Int64[] for i = 1:n_depots+n_sites]

    for (i, val) in enumerate(x_edges)

        # Val is the number of edges
        if val > 0

            idx_pred, idx_succ = vector_idx_to_depot_site[i]

            # Add to the number of out_edges
            edge_count[idx_pred] += val

            # Append the neighbor val number of times
            for _ = 1:val
                push!(adj_list[idx_pred], idx_succ)
            end
        end
    end

    # Now we have the correct edge_count and adj_list
    curr_path = Stack{Int64}()
    circuit = Vector{Int64}(undef, 0)

    # Start from the first vtx with degree > 1 (should just be 1)
    first_idx = 0
    for i = 1:n_depots
        if edge_count[i] > 0
            first_idx = i
            break
        end
    end
    # Push now because will also end at it
    push!(curr_path, first_idx)
    curr_v = first_idx

    while ~(isempty(curr_path))

        if edge_count[curr_v] > 0

            push!(curr_path, curr_v)
            next_v = pop!(adj_list[curr_v])
            edge_count[curr_v] -= 1
            curr_v = next_v

        else

            push!(circuit, curr_v)
            curr_v = pop!(curr_path)

        end
    end

    return circuit
end

"""
Remove leading and trailing depots except the one just before the first package and the one just after the last package
"""
function trim_circuit!(circuit::Vector, n_depots::Int64)

    # Trim leading depots
    first_site_idx = findfirst(x -> x > n_depots, circuit)

    if first_site_idx == nothing
        empty!(circuit)
        return
    end

    for i = 1:first_site_idx-2
        popfirst!(circuit)
    end

    last_site_idx = findlast(x -> x > n_depots, circuit)
    for i = last_site_idx+2:length(circuit)
        pop!(circuit)
    end

end

"""
Distribute the tour over m agents. The final step of MergeSplitTours
"""
function cut_tours(circuit::Vector{Int64}, n_depots::Int64, n_agents::Int64,
                   depot_sites::Vector{LOC}, cost_fn::F) where {LOC, F <: Function}

    arc_costs = [cost_fn(i, j) for (i, j) in zip(circuit[1:end-1], circuit[2:end])]
    total_cost = sum(arc_costs)

    agent_tours = [Int64[] for _ = 1:n_agents]
    idx = 1 # The idx is over circuit

    for i = 1:n_agents

        this_agent_tour = [circuit[idx]]
        tour_cost = 0

        while tour_cost <= total_cost/n_agents && idx < length(circuit)

            tour_cost += arc_costs[idx]
            idx += 1
            push!(this_agent_tour, circuit[idx])

        end

        # Continue the tour as required
        if idx < length(circuit)

            # @assert idx <= length(circuit) - 2 "Last agent has no package!"

            if circuit[idx] > n_depots # Ends at a site - add next depot
                tour_cost += arc_costs[idx]
                idx += 1
                push!(this_agent_tour, circuit[idx])
            end
        end

        agent_tours[i] = this_agent_tour

        if idx == length(circuit)
            break
        end
    end

    return agent_tours
end

"""
Combine all of the steps of MergeSplitTours in one function. Also trim the agent sub-tours after to avoid redundancies in sub-tours.
"""
function task_allocation(n_depots::Int64, n_sites::Int64, n_agents, depot_sites::Vector{LOC}, cost_fn::F) where {LOC, F <: Function}

    edges, x_edges, vds, dsv, cv = min_connecting_tour(n_depots, n_sites, depot_sites, cost_fn)

    depot_comps = get_connected_depots(n_depots, n_sites, edges)

    add_merged_depot_edges!(x_edges, depot_comps, dsv, n_depots, n_sites, cv)

    circuit = get_multiedge_eulerian_tour(x_edges, vds, n_depots, n_sites)

    trim_circuit!(circuit, n_depots)

    agent_tours = cut_tours(circuit, n_depots, n_agents, depot_sites, cost_fn)

    # Trim the agent sub-tours also
    empty_tours = Int64[]
    for (i, atour) in enumerate(agent_tours)
        if ~(isempty(atour))
            trim_circuit!(agent_tours[i], n_depots)
        end
        if isempty(agent_tours[i])
            push!(empty_tours, i)
        end
    end

    # delete empty agent tours
    deleteat!(agent_tours, empty_tours)

    return agent_tours

end

"""
(This is only used by the MAPF layer) Get the set of AgentTask instances as a tuple of (origin, site, destination)
"""
function get_agent_task_set(agent_tours::Vector{Vector{Int64}}, n_depots::Int64,
                            n_sites::Int64)

    agent_tasks = AgentTask[]

    for (i, atour) in enumerate(agent_tours)

        if length(atour) < 3 # For anything other than dpd'
            continue
        end

        @assert atour[1] <= n_depots && atour[2] > n_depots && atour[3] <= n_depots

        # HAVE TO SUBTRACT
        push!(agent_tasks, (origin=atour[1], site=atour[2]-n_depots, dest=atour[3]))
    end

    return agent_tasks
end

"""
For the given agent tour, get the next task after the finished site
"""
function get_next_agent_task(agent_tour::Vector{Int64}, n_depots::Int64,
                            n_sites::Int64, finished_site::Int64)

    idx = 1

    while agent_tour[idx] != n_depots + finished_site # Need to add n_depots here
        idx += 1
    end

    # Go to next depot and ensure more is left
    # SINCE TRIMMED, must be another site
    if idx > length(agent_tour) - 3
        return nothing
    end

    idx = idx + 2

    while idx <= length(agent_tour)-1

        if agent_tour[idx-1] <= n_depots && agent_tour[idx] > n_depots &&
            agent_tour[idx+1] <= n_depots
            break
        end
        idx = idx+1
    end

    @assert idx < length(agent_tour) "Didn't find dpd' for site $(n_depots+finished_site) in tour $(agent_tour)"

    return (origin = agent_tour[idx-1], site = agent_tour[idx] - n_depots, dest = agent_tour[idx+1])
end


text:
using BusinessDays
using Test
import Dates

# Issue #18
Dict(Any[])

# Issue #30
# list holidays for all available calendars
for c in [:BRSettlement, :BrazilExchange, :USNYSE, :USGovernmentBond, :USSettlement, :CanadaTSX, :CanadaSettlement, :EuroZone, :UKSettlement, :AustraliaASX]
    listholidays(c, Dates.Date(1900,1,1), Dates.Date(2100,1,1))
end

# Types
bhc = BusinessDays.Brazil()
ushc = BusinessDays.USSettlement()
ukhc = BusinessDays.UKSettlement()
usnysehc = BusinessDays.USNYSE()
usgovbondhc = BusinessDays.USGovernmentBond()
targethc = BusinessDays.TARGET()
hc_composite_BR_USA = CompositeHolidayCalendar([BusinessDays.Brazil(), BusinessDays.USSettlement()])
all_calendars_vec = [bhc, ushc, ukhc, hc_composite_BR_USA, usnysehc, usgovbondhc]

# two different instances of the same HolidayCalendar subtype should be equal
@test bhc == BusinessDays.Brazil()
@test ushc == BusinessDays.USSettlement()
@test bhc != ushc
@test string(bhc) == "BusinessDays.BRSettlement"

# Check typing system
@test isa(bhc,HolidayCalendar)
@test isa(ushc,HolidayCalendar)
@test isa(ukhc,HolidayCalendar)

################
## Easter Dates
################

# easter dates from http://www.maa.mhn.de/StarDate/publ_holidays.html
@test_throws ErrorException BusinessDays.easter_date(Dates.Year(1100))
@test BusinessDays.easter_date(Dates.Year(1901)) == Dates.Date(1901, 04, 07)
@test BusinessDays.easter_date(Dates.Year(1902)) == Dates.Date(1902, 03, 30)
@test BusinessDays.easter_date(Dates.Year(1903)) == Dates.Date(1903, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1904)) == Dates.Date(1904, 04, 03)
@test BusinessDays.easter_date(Dates.Year(1905)) == Dates.Date(1905, 04, 23)
@test BusinessDays.easter_date(Dates.Year(1906)) == Dates.Date(1906, 04, 15)
@test BusinessDays.easter_date(Dates.Year(1907)) == Dates.Date(1907, 03, 31)
@test BusinessDays.easter_date(Dates.Year(1908)) == Dates.Date(1908, 04, 19)
@test BusinessDays.easter_date(Dates.Year(1909)) == Dates.Date(1909, 04, 11)
@test BusinessDays.easter_date(Dates.Year(1910)) == Dates.Date(1910, 03, 27)
@test BusinessDays.easter_date(Dates.Year(1911)) == Dates.Date(1911, 04, 16)
@test BusinessDays.easter_date(Dates.Year(1912)) == Dates.Date(1912, 04, 07)
@test BusinessDays.easter_date(Dates.Year(1913)) == Dates.Date(1913, 03, 23)
@test BusinessDays.easter_date(Dates.Year(1914)) == Dates.Date(1914, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1915)) == Dates.Date(1915, 04, 04)
@test BusinessDays.easter_date(Dates.Year(1916)) == Dates.Date(1916, 04, 23)
@test BusinessDays.easter_date(Dates.Year(1917)) == Dates.Date(1917, 04, 08)
@test BusinessDays.easter_date(Dates.Year(1918)) == Dates.Date(1918, 03, 31)
@test BusinessDays.easter_date(Dates.Year(1919)) == Dates.Date(1919, 04, 20)
@test BusinessDays.easter_date(Dates.Year(1920)) == Dates.Date(1920, 04, 04)
@test BusinessDays.easter_date(Dates.Year(1921)) == Dates.Date(1921, 03, 27)
@test BusinessDays.easter_date(Dates.Year(1922)) == Dates.Date(1922, 04, 16)
@test BusinessDays.easter_date(Dates.Year(1923)) == Dates.Date(1923, 04, 01)
@test BusinessDays.easter_date(Dates.Year(1924)) == Dates.Date(1924, 04, 20)
@test BusinessDays.easter_date(Dates.Year(1925)) == Dates.Date(1925, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1926)) == Dates.Date(1926, 04, 04)
@test BusinessDays.easter_date(Dates.Year(1927)) == Dates.Date(1927, 04, 17)
@test BusinessDays.easter_date(Dates.Year(1928)) == Dates.Date(1928, 04, 08)
@test BusinessDays.easter_date(Dates.Year(1929)) == Dates.Date(1929, 03, 31)
@test BusinessDays.easter_date(Dates.Year(1930)) == Dates.Date(1930, 04, 20)
@test BusinessDays.easter_date(Dates.Year(1931)) == Dates.Date(1931, 04, 05)
@test BusinessDays.easter_date(Dates.Year(1932)) == Dates.Date(1932, 03, 27)
@test BusinessDays.easter_date(Dates.Year(1933)) == Dates.Date(1933, 04, 16)
@test BusinessDays.easter_date(Dates.Year(1934)) == Dates.Date(1934, 04, 01)
@test BusinessDays.easter_date(Dates.Year(1935)) == Dates.Date(1935, 04, 21)
@test BusinessDays.easter_date(Dates.Year(1936)) == Dates.Date(1936, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1937)) == Dates.Date(1937, 03, 28)
@test BusinessDays.easter_date(Dates.Year(1938)) == Dates.Date(1938, 04, 17)
@test BusinessDays.easter_date(Dates.Year(1939)) == Dates.Date(1939, 04, 09)
@test BusinessDays.easter_date(Dates.Year(1940)) == Dates.Date(1940, 03, 24)
@test BusinessDays.easter_date(Dates.Year(1941)) == Dates.Date(1941, 04, 13)
@test BusinessDays.easter_date(Dates.Year(1942)) == Dates.Date(1942, 04, 05)
@test BusinessDays.easter_date(Dates.Year(1943)) == Dates.Date(1943, 04, 25)
@test BusinessDays.easter_date(Dates.Year(1944)) == Dates.Date(1944, 04, 09)
@test BusinessDays.easter_date(Dates.Year(1945)) == Dates.Date(1945, 04, 01)
@test BusinessDays.easter_date(Dates.Year(1946)) == Dates.Date(1946, 04, 21)
@test BusinessDays.easter_date(Dates.Year(1947)) == Dates.Date(1947, 04, 06)
@test BusinessDays.easter_date(Dates.Year(1948)) == Dates.Date(1948, 03, 28)
@test BusinessDays.easter_date(Dates.Year(1949)) == Dates.Date(1949, 04, 17)
@test BusinessDays.easter_date(Dates.Year(1950)) == Dates.Date(1950, 04, 09)
@test BusinessDays.easter_date(Dates.Year(1951)) == Dates.Date(1951, 03, 25)
@test BusinessDays.easter_date(Dates.Year(1952)) == Dates.Date(1952, 04, 13)
@test BusinessDays.easter_date(Dates.Year(1953)) == Dates.Date(1953, 04, 05)
@test BusinessDays.easter_date(Dates.Year(1954)) == Dates.Date(1954, 04, 18)
@test BusinessDays.easter_date(Dates.Year(1955)) == Dates.Date(1955, 04, 10)
@test BusinessDays.easter_date(Dates.Year(1956)) == Dates.Date(1956, 04, 01)
@test BusinessDays.easter_date(Dates.Year(1957)) == Dates.Date(1957, 04, 21)
@test BusinessDays.easter_date(Dates.Year(1958)) == Dates.Date(1958, 04, 06)
@test BusinessDays.easter_date(Dates.Year(1959)) == Dates.Date(1959, 03, 29)
@test BusinessDays.easter_date(Dates.Year(1960)) == Dates.Date(1960, 04, 17)
@test BusinessDays.easter_date(Dates.Year(1961)) == Dates.Date(1961, 04, 02)
@test BusinessDays.easter_date(Dates.Year(1962)) == Dates.Date(1962, 04, 22)
@test BusinessDays.easter_date(Dates.Year(1963)) == Dates.Date(1963, 04, 14)
@test BusinessDays.easter_date(Dates.Year(1964)) == Dates.Date(1964, 03, 29)
@test BusinessDays.easter_date(Dates.Year(1965)) == Dates.Date(1965, 04, 18)
@test BusinessDays.easter_date(Dates.Year(1966)) == Dates.Date(1966, 04, 10)
@test BusinessDays.easter_date(Dates.Year(1967)) == Dates.Date(1967, 03, 26)
@test BusinessDays.easter_date(Dates.Year(1968)) == Dates.Date(1968, 04, 14)
@test BusinessDays.easter_date(Dates.Year(1969)) == Dates.Date(1969, 04, 06)
@test BusinessDays.easter_date(Dates.Year(1970)) == Dates.Date(1970, 03, 29)
@test BusinessDays.easter_date(Dates.Year(1971)) == Dates.Date(1971, 04, 11)
@test BusinessDays.easter_date(Dates.Year(1972)) == Dates.Date(1972, 04, 02)
@test BusinessDays.easter_date(Dates.Year(1973)) == Dates.Date(1973, 04, 22)
@test BusinessDays.easter_date(Dates.Year(1974)) == Dates.Date(1974, 04, 14)
@test BusinessDays.easter_date(Dates.Year(1975)) == Dates.Date(1975, 03, 30)
@test BusinessDays.easter_date(Dates.Year(1976)) == Dates.Date(1976, 04, 18)
@test BusinessDays.easter_date(Dates.Year(1977)) == Dates.Date(1977, 04, 10)
@test BusinessDays.easter_date(Dates.Year(1978)) == Dates.Date(1978, 03, 26)
@test BusinessDays.easter_date(Dates.Year(1979)) == Dates.Date(1979, 04, 15)
@test BusinessDays.easter_date(Dates.Year(1980)) == Dates.Date(1980, 04, 06)
@test BusinessDays.easter_date(Dates.Year(1981)) == Dates.Date(1981, 04, 19)
@test BusinessDays.easter_date(Dates.Year(1982)) == Dates.Date(1982, 04, 11)
@test BusinessDays.easter_date(Dates.Year(1983)) == Dates.Date(1983, 04, 03)
@test BusinessDays.easter_date(Dates.Year(1984)) == Dates.Date(1984, 04, 22)
@test BusinessDays.easter_date(Dates.Year(1985)) == Dates.Date(1985, 04, 07)
@test BusinessDays.easter_date(Dates.Year(1986)) == Dates.Date(1986, 03, 30)
@test BusinessDays.easter_date(Dates.Year(1987)) == Dates.Date(1987, 04, 19)
@test BusinessDays.easter_date(Dates.Year(1988)) == Dates.Date(1988, 04, 03)
@test BusinessDays.easter_date(Dates.Year(1989)) == Dates.Date(1989, 03, 26)
@test BusinessDays.easter_date(Dates.Year(1990)) == Dates.Date(1990, 04, 15)
@test BusinessDays.easter_date(Dates.Year(1991)) == Dates.Date(1991, 03, 31)
@test BusinessDays.easter_date(Dates.Year(1992)) == Dates.Date(1992, 04, 19)
@test BusinessDays.easter_date(Dates.Year(1993)) == Dates.Date(1993, 04, 11)
@test BusinessDays.easter_date(Dates.Year(1994)) == Dates.Date(1994, 04, 03)
@test BusinessDays.easter_date(Dates.Year(1995)) == Dates.Date(1995, 04, 16)
@test BusinessDays.easter_date(Dates.Year(1996)) == Dates.Date(1996, 04, 07)
@test BusinessDays.easter_date(Dates.Year(1997)) == Dates.Date(1997, 03, 30)
@test BusinessDays.easter_date(Dates.Year(1998)) == Dates.Date(1998, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1999)) == Dates.Date(1999, 04, 04)
@test BusinessDays.easter_date(Dates.Year(2000)) == Dates.Date(2000, 04, 23)
@test BusinessDays.easter_date(Dates.Year(2001)) == Dates.Date(2001, 04, 15)
@test BusinessDays.easter_date(Dates.Year(2002)) == Dates.Date(2002, 03, 31)
@test BusinessDays.easter_date(Dates.Year(2003)) == Dates.Date(2003, 04, 20)
@test BusinessDays.easter_date(Dates.Year(2004)) == Dates.Date(2004, 04, 11)
@test BusinessDays.easter_date(Dates.Year(2005)) == Dates.Date(2005, 03, 27)
@test BusinessDays.easter_date(Dates.Year(2006)) == Dates.Date(2006, 04, 16)
@test BusinessDays.easter_date(Dates.Year(2007)) == Dates.Date(2007, 04, 08)
@test BusinessDays.easter_date(Dates.Year(2008)) == Dates.Date(2008, 03, 23)
@test BusinessDays.easter_date(Dates.Year(2009)) == Dates.Date(2009, 04, 12)
@test BusinessDays.easter_date(Dates.Year(2010)) == Dates.Date(2010, 04, 04)
@test BusinessDays.easter_date(Dates.Year(2011)) == Dates.Date(2011, 04, 24)
@test BusinessDays.easter_date(Dates.Year(2012)) == Dates.Date(2012, 04, 08)
@test BusinessDays.easter_date(Dates.Year(2013)) == Dates.Date(2013, 03, 31)
@test BusinessDays.easter_date(Dates.Year(2014)) == Dates.Date(2014, 04, 20)
@test BusinessDays.easter_date(Dates.Year(2015)) == Dates.Date(2015, 04, 05)
@test BusinessDays.easter_date(Dates.Year(2016)) == Dates.Date(2016, 03, 27)
@test BusinessDays.easter_date(Dates.Year(2017)) == Dates.Date(2017, 04, 16)
@test BusinessDays.easter_date(Dates.Year(2018)) == Dates.Date(2018, 04, 01)
@test BusinessDays.easter_date(Dates.Year(2019)) == Dates.Date(2019, 04, 21)
@test BusinessDays.easter_date(Dates.Year(2020)) == Dates.Date(2020, 04, 12)
@test BusinessDays.easter_date(Dates.Year(2021)) == Dates.Date(2021, 04, 04)
@test BusinessDays.easter_date(Dates.Year(2022)) == Dates.Date(2022, 04, 17)
@test BusinessDays.easter_date(Dates.Year(2023)) == Dates.Date(2023, 04, 09)
@test BusinessDays.easter_date(Dates.Year(2024)) == Dates.Date(2024, 03, 31)
@test BusinessDays.easter_date(Dates.Year(2025)) == Dates.Date(2025, 04, 20)
@test BusinessDays.easter_date(Dates.Year(2026)) == Dates.Date(2026, 04, 05)
@test BusinessDays.easter_date(Dates.Year(2027)) == Dates.Date(2027, 03, 28)
@test BusinessDays.easter_date(Dates.Year(2028)) == Dates.Date(2028, 04, 16)
@test BusinessDays.easter_date(Dates.Year(2029)) == Dates.Date(2029, 04, 01)
@test BusinessDays.easter_date(Dates.Year(2030)) == Dates.Date(2030, 04, 21)
@test BusinessDays.easter_date(Dates.Year(2031)) == Dates.Date(2031, 04, 13)
@test BusinessDays.easter_date(Dates.Year(2032)) == Dates.Date(2032, 03, 28)
@test BusinessDays.easter_date(Dates.Year(2033)) == Dates.Date(2033, 04, 17)
@test BusinessDays.easter_date(Dates.Year(2034)) == Dates.Date(2034, 04, 09)
@test BusinessDays.easter_date(Dates.Year(2035)) == Dates.Date(2035, 03, 25)
@test BusinessDays.easter_date(Dates.Year(2036)) == Dates.Date(2036, 04, 13)
@test BusinessDays.easter_date(Dates.Year(2037)) == Dates.Date(2037, 04, 05)
@test BusinessDays.easter_date(Dates.Year(2038)) == Dates.Date(2038, 04, 25)
@test BusinessDays.easter_date(Dates.Year(2039)) == Dates.Date(2039, 04, 10)
@test BusinessDays.easter_date(Dates.Year(2040)) == Dates.Date(2040, 04, 01)
@test BusinessDays.easter_date(Dates.Year(2041)) == Dates.Date(2041, 04, 21)
@test BusinessDays.easter_date(Dates.Year(2042)) == Dates.Date(2042, 04, 06)
@test BusinessDays.easter_date(Dates.Year(2043)) == Dates.Date(2043, 03, 29)
@test BusinessDays.easter_date(Dates.Year(2044)) == Dates.Date(2044, 04, 17)
@test BusinessDays.easter_date(Dates.Year(2045)) == Dates.Date(2045, 04, 09)
@test BusinessDays.easter_date(Dates.Year(2046)) == Dates.Date(2046, 03, 25)
@test BusinessDays.easter_date(Dates.Year(2047)) == Dates.Date(2047, 04, 14)
@test BusinessDays.easter_date(Dates.Year(2048)) == Dates.Date(2048, 04, 05)
@test BusinessDays.easter_date(Dates.Year(2049)) == Dates.Date(2049, 04, 18)
@test BusinessDays.easter_date(Dates.Year(2050)) == Dates.Date(2050, 04, 10)
@test BusinessDays.easter_date(Dates.Year(2051)) == Dates.Date(2051, 04, 02)
@test BusinessDays.easter_date(Dates.Year(2052)) == Dates.Date(2052, 04, 21)
@test BusinessDays.easter_date(Dates.Year(2053)) == Dates.Date(2053, 04, 06)
@test BusinessDays.easter_date(Dates.Year(2054)) == Dates.Date(2054, 03, 29)
@test BusinessDays.easter_date(Dates.Year(2055)) == Dates.Date(2055, 04, 18)
@test BusinessDays.easter_date(Dates.Year(2056)) == Dates.Date(2056, 04, 02)
@test BusinessDays.easter_date(Dates.Year(2057)) == Dates.Date(2057, 04, 22)
@test BusinessDays.easter_date(Dates.Year(2058)) == Dates.Date(2058, 04, 14)
@test BusinessDays.easter_date(Dates.Year(2059)) == Dates.Date(2059, 03, 30)
@test BusinessDays.easter_date(Dates.Year(2060)) == Dates.Date(2060, 04, 18)
@test BusinessDays.easter_date(Dates.Year(2061)) == Dates.Date(2061, 04, 10)
@test BusinessDays.easter_date(Dates.Year(2062)) == Dates.Date(2062, 03, 26)
@test BusinessDays.easter_date(Dates.Year(2063)) == Dates.Date(2063, 04, 15)
@test BusinessDays.easter_date(Dates.Year(2064)) == Dates.Date(2064, 04, 06)
@test BusinessDays.easter_date(Dates.Year(2065)) == Dates.Date(2065, 03, 29)
@test BusinessDays.easter_date(Dates.Year(2066)) == Dates.Date(2066, 04, 11)
@test BusinessDays.easter_date(Dates.Year(2067)) == Dates.Date(2067, 04, 03)
@test BusinessDays.easter_date(Dates.Year(2068)) == Dates.Date(2068, 04, 22)
@test BusinessDays.easter_date(Dates.Year(2069)) == Dates.Date(2069, 04, 14)
@test BusinessDays.easter_date(Dates.Year(2070)) == Dates.Date(2070, 03, 30)
@test BusinessDays.easter_date(Dates.Year(2071)) == Dates.Date(2071, 04, 19)
@test BusinessDays.easter_date(Dates.Year(2072)) == Dates.Date(2072, 04, 10)
@test BusinessDays.easter_date(Dates.Year(2073)) == Dates.Date(2073, 03, 26)
@test BusinessDays.easter_date(Dates.Year(2074)) == Dates.Date(2074, 04, 15)
@test BusinessDays.easter_date(Dates.Year(2075)) == Dates.Date(2075, 04, 07)
@test BusinessDays.easter_date(Dates.Year(2076)) == Dates.Date(2076, 04, 19)
@test BusinessDays.easter_date(Dates.Year(2077)) == Dates.Date(2077, 04, 11)
@test BusinessDays.easter_date(Dates.Year(2078)) == Dates.Date(2078, 04, 03)

include("easter-min-max.jl")

###############
# findweekday
###############
# weekday values:
# const Monday,Tuesday,Wednesday,Thursday,Friday,Saturday,Sunday = 1,2,3,4,5,6,7
# see query.jl on Dates module
# See also dayofweek(dt) function.
# function findweekday(weekday_target :: Int, yy :: Int, mm:: Int, occurrence :: Int, ascending :: Bool )
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 1, true) == Dates.Date(2015, 07, 06)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 2, true) == Dates.Date(2015, 07, 13)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 3, true) == Dates.Date(2015, 07, 20)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 4, true) == Dates.Date(2015, 07, 27)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 5, true) == Dates.Date(2015, 08, 03)

@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 1, false) == Dates.Date(2015, 07, 27)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 2, false) == Dates.Date(2015, 07, 20)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 3, false) == Dates.Date(2015, 07, 13)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 4, false) == Dates.Date(2015, 07, 06)
@test BusinessDays.findweekday(Dates.Monday, 2015, 07, 5, false) == Dates.Date(2015, 06, 29)

@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 1, true) == Dates.Date(2015, 07, 03)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 2, true) == Dates.Date(2015, 07, 10)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 3, true) == Dates.Date(2015, 07, 17)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 4, true) == Dates.Date(2015, 07, 24)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 5, true) == Dates.Date(2015, 07, 31)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 6, true) == Dates.Date(2015, 08, 07)

@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 1, false) == Dates.Date(2015, 07, 31)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 2, false) == Dates.Date(2015, 07, 24)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 3, false) == Dates.Date(2015, 07, 17)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 4, false) == Dates.Date(2015, 07, 10)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 5, false) == Dates.Date(2015, 07, 03)
@test BusinessDays.findweekday(Dates.Friday, 2015, 07, 6, false) == Dates.Date(2015, 06, 26)

@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 1, true) == Dates.Date(2015, 07, 01)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 2, true) == Dates.Date(2015, 07, 08)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 3, true) == Dates.Date(2015, 07, 15)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 4, true) == Dates.Date(2015, 07, 22)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 5, true) == Dates.Date(2015, 07, 29)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 6, true) == Dates.Date(2015, 08, 05)

@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 1, false) == Dates.Date(2015, 07, 29)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 2, false) == Dates.Date(2015, 07, 22)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 3, false) == Dates.Date(2015, 07, 15)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 4, false) == Dates.Date(2015, 07, 08)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 5, false) == Dates.Date(2015, 07, 01)
@test BusinessDays.findweekday(Dates.Wednesday, 2015, 07, 6, false) == Dates.Date(2015, 06, 24)

@test_throws AssertionError BusinessDays.findweekday(Dates.Wednesday, 2015, 07, -1, false)
@test_throws AssertionError BusinessDays.findweekday(Dates.Wednesday, 2015, 07, -1, true)

# Test this only on 64bit or higher systems
len = typemax(UInt32) + 1
if len > typemax(UInt32)
    d0 = Dates.Date(1950, 2, 1)
    d1 = d0 + Dates.Day(len)
    @test_throws AssertionError BusinessDays._create_bdays_cache_arrays(bhc, d0, d1)
end

# Create HolidayCalendar instances
hc_australiaasx = BusinessDays.AustraliaASX()
hc_australiaact = BusinessDays.Australia(:ACT)
hc_australiansw = BusinessDays.Australia(:NSW)
hc_australiant  = BusinessDays.Australia(:NT)
hc_australiaqld = BusinessDays.Australia(:QLD)
hc_australiasa  = BusinessDays.Australia(:SA)
hc_australiatas = BusinessDays.Australia(:TAS)
hc_australiawa  = BusinessDays.Australia(:WA)
hc_australiavic = BusinessDays.Australia(:VIC)
hc_brazil       = BusinessDays.Brazil()
hc_brazil_bmf   = BusinessDays.BrazilBMF()
hc_usa          = BusinessDays.USSettlement()
hc_uk           = BusinessDays.UKSettlement()
hc_usnyse       = BusinessDays.USNYSE()
hc_usgovbond    = BusinessDays.USGovernmentBond()
hc_canadatsx    = BusinessDays.CanadaTSX()
hc_canada       = BusinessDays.CanadaSettlement()

####################
# Calendar Tests
####################

usecache = false
include("calendar_tests.jl")

usecache = true
include("calendar_tests.jl")

####################
# Performance Tests
####################
include("perftests.jl")

BusinessDays.cleancache(hc_brazil)
BusinessDays.cleancache()

struct TestCalendar <: HolidayCalendar end
cc = TestCalendar()
@test_throws ErrorException isholiday(cc, Dates.Date(2015,1,1))
BusinessDays.isholiday(::TestCalendar, dt::Dates.Date) = dt == Dates.Date(2015,8,27)
@test isholiday(cc, Dates.Date(2015,8,26)) == false
@test isholiday(cc, Dates.Date(2015,8,27)) == true
BusinessDays.initcache(cc)
@test isholiday(cc, Dates.Date(2015,8,26)) == false
@test isholiday(cc, Dates.Date(2015,8,27)) == true
BusinessDays.cleancache(cc)
isholiday(:Brazil, Dates.Date(2016,2,1))
isholiday(:TestCalendar, Dates.Date(2016,2,1))
isholiday("TestCalendar", Dates.Date(2016,2,1))
isbday(:Brazil, Dates.Date(2016,2,1))
isbday(:TestCalendar, Dates.Date(2016,2,1))
isbday("TestCalendar", Dates.Date(2016,2,1))

sym_vec = [:Brazil, :UKSettlement]
BusinessDays.initcache(sym_vec)
BusinessDays.cleancache(sym_vec)

str_vec = ["Brazil", "UKSettlement", "Canada", "UnitedStates", "USNYSE", "USGovernmentBond", "CanadaTSX", "WeekendsOnly"]
BusinessDays.initcache(str_vec)
BusinessDays.cleancache(str_vec)

BusinessDays.initcache("UKSettlement", Dates.Date(2000,1,1), Dates.Date(2000,5,2))
BusinessDays.initcache("UKSettlement", Dates.Date(2000,1,1), Dates.Date(2000,5,2)) # repeating initcache should work
BusinessDays.initcache("UKSettlement", Dates.Date(2000,1,1), Dates.Date(2000,1,1)) # single date cache should work
BusinessDays.cleancache("UKSettlement")

# equality and hash function for Australia
@test BusinessDays.Australia(:ACT) == BusinessDays.Australia(:ACT)
@test BusinessDays.Australia(:ACT) != BusinessDays.Australia(:NSW)
@test hash(BusinessDays.Australia(:ACT)) == hash(BusinessDays.Australia(:ACT))
@test hash(BusinessDays.Australia(:ACT)) != hash(BusinessDays.Australia(:NSW))

BusinessDays.initcache(BusinessDays.Australia(:ACT))
@test haskey(BusinessDays.CACHE_DICT, BusinessDays.Australia(:ACT))
@test !haskey(BusinessDays.CACHE_DICT, BusinessDays.Australia(:NSW))
BusinessDays.cleancache(BusinessDays.Australia(:ACT))

include("customcalendar-example.jl")

#########################
# GenericHolidayCalendar
#########################

d0 = Dates.Date(1980,1,1)
d1 = Dates.Date(2050,1,1)
gen_brazil = GenericHolidayCalendar(listholidays(hc_brazil, d0, d1))
@test isbday(gen_brazil, Dates.Date(2014, 12, 31)) == true # wednesday
@test isbday(gen_brazil, Dates.Date(2015, 01, 01)) == false # new year
@test isbday(gen_brazil, Dates.Date(2015, 01, 02)) == true # friday

@test advancebdays(gen_brazil, Dates.Date(2015,9,1), [0, 1, 3, 4, 5]) == advancebdays(hc_brazil, Dates.Date(2015,9,1), [0, 1, 3, 4, 5])
@test advancebdays(gen_brazil, Dates.Date(2015,9,1), 0:5) == advancebdays(hc_brazil, Dates.Date(2015,9,1), 0:5)
@test listholidays(gen_brazil, Dates.Date(2016,1,1), Dates.Date(2016,5,30)) == listholidays(hc_brazil, Dates.Date(2016,1,1), Dates.Date(2016,5,30))

@test tobday(gen_brazil, Dates.Date(2013, 02, 08)) == tobday(hc_brazil, Dates.Date(2013, 02, 08))
@test tobday(gen_brazil, Dates.Date(2013, 02, 09)) == tobday(hc_brazil, Dates.Date(2013, 02, 09))
@test_throws AssertionError bdays(gen_brazil, Dates.Date(1900,2,1), Dates.Date(2000,2,1))
@test_throws AssertionError bdays(gen_brazil, Dates.Date(2000,2,1), Dates.Date(2100,2,1))

d0_test = Dates.Date(1980,1,2)
d1_test = Dates.Date(2049,12,28)
@test bdays(hc_brazil, d0_test, d1_test) == bdays(gen_brazil, d0_test, d1_test)
@test bdayscount(hc_brazil, d0_test, d1_test) == bdayscount(gen_brazil, d0_test, d1_test)
println("a million with GenericHolidayCalendar...")
@time for i in 1:1000000 bdays(gen_brazil, d0_test, d1_test) end
BusinessDays.cleancache(gen_brazil)

# Start all over, but without cache
gen_brazil = GenericHolidayCalendar(listholidays(hc_brazil, d0, d1), d0, d1, false)

@test isbday(gen_brazil, Dates.Date(2014, 12, 31)) == true # wednesday
@test isbday(gen_brazil, Dates.Date(2015, 01, 01)) == false # new year
@test isbday(gen_brazil, Dates.Date(2015, 01, 02)) == true # friday

@test advancebdays(gen_brazil, Dates.Date(2015,9,1), [0, 1, 3, 4, 5]) == advancebdays(hc_brazil, Dates.Date(2015,9,1), [0, 1, 3, 4, 5])
@test advancebdays(gen_brazil, Dates.Date(2015,9,1), 0:5) == advancebdays(hc_brazil, Dates.Date(2015,9,1), 0:5)
@test listholidays(gen_brazil, Dates.Date(2016,1,1), Dates.Date(2016,5,30)) == listholidays(hc_brazil, Dates.Date(2016,1,1), Dates.Date(2016,5,30))

@test tobday(gen_brazil, Dates.Date(2013, 02, 08)) == tobday(hc_brazil, Dates.Date(2013, 02, 08))
@test tobday(gen_brazil, Dates.Date(2013, 02, 09)) == tobday(hc_brazil, Dates.Date(2013, 02, 09))
@test_throws AssertionError bdays(gen_brazil, Dates.Date(1900,2,1), Dates.Date(2000,2,1))
@test_throws AssertionError bdays(gen_brazil, Dates.Date(2000,2,1), Dates.Date(2100,2,1))

@test BusinessDays.needs_cache_update(gen_brazil, d0, d1) == false
BusinessDays.initcache(gen_brazil)
@test BusinessDays.needs_cache_update(gen_brazil, d0, d1) == true
@test isbday(gen_brazil, Dates.Date(2014, 12, 31)) == true # wednesday
@test isbday(gen_brazil, Dates.Date(2015, 01, 01)) == false # new year
@test isbday(gen_brazil, Dates.Date(2015, 01, 02)) == true # friday

# does nothing, because cache is already there
BusinessDays.initcache(gen_brazil)

# A GenericHolidayCalendar is defined by its set of holidays, dtmin, dtmax
dtmin = Dates.Date(2018,1,15)
dtmax = Dates.Date(2018,1,19)
gen_1 = GenericHolidayCalendar(Set([Dates.Date(2018,1,16), Dates.Date(2018,1,18)]), dtmin, dtmax)
gen_2 = GenericHolidayCalendar(Set([Dates.Date(2018,1,16), Dates.Date(2018,1,18)]), dtmin, dtmax)
@test gen_1 == gen_2

# On a set, they are the same element
cal_set = Set([gen_1, gen_2])
@test length(cal_set) == 1

# On a dict, they represent the same key
cal_dict = Dict(gen_1 => "hey")
@test cal_dict[gen_2] == "hey"