text:
function test_1D()
    la1 = set_Lattice(1,[[1.0]])
    add_atoms!(la1,[0,0])
    t = 1.0
    add_hoppings!(la1,-t,1,1,[1])
    ham1 = hamiltonian_k(la1)
    kmin = [-π]
    kmax = [π]
    nk = 20
    vec_k,energies = calc_band(kmin,kmax,nk,la1,ham1)
    println(energies)
    return energies

    pls = plot(vec_k[1,:],energies[1,:],marker=:circle,label=["1D"])
    savefig("1Denergy.png")

end

function test_2Dsquare()
    println("Test1: ")
    println("2D Square lattice with nearest neighbor hoppings")
    la2 = set_Lattice(2,[[1,0],[0,1]])
    add_atoms!(la2,[0,0])

    show_neighbors(la2)

    t = 1.0
    add_hoppings!(la2,-t,1,1,[1,0])
    add_hoppings!(la2,-t,1,1,[0,1])
    ham2 = hamiltonian_k(la2)

    kmin = [-π,-π]
    kmax = [0.0,0.0]
    nk = 20
    vec_k,energies = calc_band(kmin,kmax,nk,la2,ham2)
    println("Energies on the line from (-π,π) to (0,0)")
    println(energies)
    return energies

    klines = set_Klines()
    kmin = [0,0]
    kmax = [π,0]

    add_Kpoints!(klines,kmin,kmax,"(0,0)","(pi,0)")

    kmin = [π,0]
    kmax = [π,π]
    add_Kpoints!(klines,kmin,kmax,"(pi,0)","(pi,pi)")

    kmin = [π,π]
    kmax = [0,0]
    add_Kpoints!(klines,kmin,kmax,"(pi,pi)","(pi,0)")

    pls = calc_band_plot(klines,la2)
    return pls

    #pls = plot(vec_k[:],energies[1,:],marker=:circle,label=["2D"])
    #savefig("2Denergy.png")

end

function test_pnictides()
    la = set_Lattice(2,[[1,0],[0,1]])
    add_atoms!(la,[0,0])
    add_atoms!(la,[0,0])
    t1 = -1.0
    t2 = 1.3
    t3 = -0.85
    t4 = t3
    μ = 1.45

    #dxz
    add_hoppings!(la,-t1,1,1,[1,0])
    add_hoppings!(la,-t2,1,1,[0,1])
    add_hoppings!(la,-t3,1,1,[1,1])
    add_hoppings!(la,-t3,1,1,[1,-1])

    #dyz
    add_hoppings!(la,-t2,2,2,[1,0])
    add_hoppings!(la,-t1,2,2,[0,1])
    add_hoppings!(la,-t3,2,2,[1,1])
    add_hoppings!(la,-t3,2,2,[1,-1])

    #between dxz and dyz
    add_hoppings!(la,-t4,1,2,[1,1])
    add_hoppings!(la,-t4,1,2,[-1,-1])
    add_hoppings!(la,t4,1,2,[1,-1])
    add_hoppings!(la,t4,1,2,[-1,1])

    set_μ!(la,μ) #set the chemical potential
#    add_diagonals!(la,[-μ,-μ])

    pls = plot_fermisurface_2D(la)
    #return pls

    klines = set_Klines()
    kmin = [0,0]
    kmax = [π,0]
    add_Kpoints!(klines,kmin,kmax,"(0,0)","(pi,0)")

    kmin = [π,0]
    kmax = [π,π]
    add_Kpoints!(klines,kmin,kmax,"(pi,0)","(pi,pi)")

    kmin = [π,π]
    kmax = [0,0]
    add_Kpoints!(klines,kmin,kmax,"(pi,pi)","(0,0)")

    pls = calc_band_plot(klines,la)
    #savefig("Fe.png")
    #return pls

    ham = hamiltonian_k(la)
    kx = 0.1
    ky = 0.2
    hamk = ham([kx,ky])
    #println(hamk)
    return hamk

end

function test_pnictides_5orbitals()
    la = set_Lattice(2,[[1,0],[0,1]])
    add_atoms!(la,[0,0])
    add_atoms!(la,[0,0])
    add_atoms!(la,[0,0])
    add_atoms!(la,[0,0])
    add_atoms!(la,[0,0])

    tmat = [
    -0.7    0 -0.4  0.2 -0.1
    -0.8    0    0    0    0
     0.8 -1.5    0    0 -0.3
       0  1.7    0    0 -0.1
    -3.0    0    0 -0.2    0
    -2.1  1.5    0    0    0
     1.3    0  0.2 -0.2    0
     1.7    0    0  0.2    0
    -2.5  1.4    0    0    0
    -2.1  3.3    0 -0.3  0.7
     1.7  0.2    0  0.2    0
     2.5    0    0  0.3    0
     1.6  1.2 -0.3 -0.3 -0.3
       0    0    0 -0.1    0
     3.1 -0.7 -0.2    0    0
    ]
    tmat = 0.1.*tmat
    imap = zeros(Int64,5,5)
    count = 0
    for μ=1:5
        for ν=μ:5
            count += 1
            imap[μ,ν] = count
        end
    end
    Is = [1,-1,-1,1,1,1,1,-1,-1,1,-1,-1,1,1,1]
    σds = [1,-1,1,1,-1,1,-1,-1,1,1,1,-1,1,-1,1]
    tmat_σy = tmat[:,:]
    tmat_σy[imap[1,2],:] = -tmat[imap[1,3],:] #(1,1)
    tmat_σy[imap[1,3],:] = -tmat[imap[1,2],:] #(1,2)
    tmat_σy[imap[1,4],:] = -tmat[imap[1,4],:] #(1,3)
    tmat_σy[imap[2,2],:] = tmat[imap[3,3],:] #(2,2)
    tmat_σy[imap[2,4],:] = tmat[imap[3,4],:] #()
    tmat_σy[imap[2,5],:] = -tmat[imap[3,5],:] #()
    tmat_σy[imap[3,3],:] = tmat[imap[2,2],:] #()
    tmat_σy[imap[3,4],:] = tmat[imap[2,4],:] #()
    tmat_σy[imap[3,5],:] = -tmat[imap[2,5],:] #()
    tmat_σy[imap[4,5],:] = -tmat[imap[4,5],:] #()

    hoppingmatrix = zeros(Float64,5,5,5,5)
    hops = [-2,-1,0,1,2]
    hopelements = [[1,0],[1,1],[2,0],[2,1],[2,2]]

    #println(tmat)

#    count = 0
    for μ = 1:5
        for ν=μ:5
            for ii=1:5
                ihop = hopelements[ii][1]
                jhop = hopelements[ii][2]
                #[a,b],[a,-b],[-a,-b],[-a,b],[b,a],[b,-a],[-b,a],[-b,-a]

                #[a,b]
                i = ihop +3
                j = jhop +3
                hoppingmatrix[μ,ν,i,j]=tmat[imap[μ,ν],ii]
                #[a,-b] = σy*[a,b] [1,1] -> [1,-1]
                if jhop != 0
                    i = ihop +3
                    j = -jhop +3
                    hoppingmatrix[μ,ν,i,j]=tmat_σy[imap[μ,ν],ii]
                end

                if μ != ν
                    #[-a,-b] = I*[a,b] [1,1] -> [-1,-1],[1,0]->[-1,0]
                    i = -ihop +3
                    j = -jhop +3
                    hoppingmatrix[μ,ν,i,j]=Is[imap[μ,ν]]*tmat[imap[μ,ν],ii]
                    #[-a,b] = I*[a,-b] = I*σy*[a,b]  #[2,0]->[-2,0]
                    if jhop != 0
                        i = -ihop +3
                        j = jhop +3
                        hoppingmatrix[μ,ν,i,j]=Is[imap[μ,ν]]*tmat_σy[imap[μ,ν],ii]
                    end
                end
                #[b,a],[b,-a],[-b,a],[-b,-a]
                if jhop != ihop
                    #[b,a] = σd*[a,b]
                    i = jhop +3
                    j = ihop +3
                    hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*tmat[imap[μ,ν],ii]
                    #[-b,a] = σd*σy*[a,b]
                    if jhop != 0
                        i = -jhop +3
                        j = ihop +3
                        hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*tmat_σy[imap[μ,ν],ii]
                    end

                    if μ != ν
                        #[-b,-a] = σd*[-a,-b] = σd*I*[a,b]
                        i = -jhop +3
                        j = -ihop +3
                        hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*Is[imap[μ,ν]]*tmat[imap[μ,ν],ii]
                        #[b,-a] = σd*[-a,b] = σd*I*[a,-b] = σd*I*σy*[a,b]  #[2,0]->[-2,0]
                        if jhop != 0
                            i = jhop +3
                            j = -ihop +3
                            hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*Is[imap[μ,ν]]*tmat_σy[imap[μ,ν],ii]
                        end
                    end
                end
            end

        end
    end

    for μ=1:5
        for ν=μ:5
            for i = 1:5
                ih = hops[i]
                for j = 1:5
                    jh = hops[j]
                    if hoppingmatrix[μ,ν,i,j] != 0.0
                        #println(ih,"\t",jh,"\t",hoppingmatrix[μ,ν,i,j] )
                        add_hoppings!(la,hoppingmatrix[μ,ν,i,j],μ,ν,[ih,jh])
                    end
                end
            end
        end
    end

    onsite = [10.75,10.96,10.96,11.12,10.62]
    set_onsite!(la,onsite)

    set_μ!(la,10.96) #set the chemical potential
#    add_diagonals!(la,[-μ,-μ])
#    ham = hamiltonian_k(la)
    #println(ham([0,0]))

    #pls = plot_fermisurface_2D(la)
    #return pls

    nk = 100
    klines = set_Klines()
    kmin = [0,0]
    kmax = [π,0]
    add_Kpoints!(klines,kmin,kmax,"(0,0)","(pi,0)",nk=nk)

    kmin = [π,0]
    kmax = [π,π]
    add_Kpoints!(klines,kmin,kmax,"(pi,0)","(pi,pi)",nk=nk)

    kmin = [π,π]
    kmax = [0,0]
    add_Kpoints!(klines,kmin,kmax,"(pi,pi)","(0,0)",nk=nk)

    pls = calc_band_plot(klines,la)
    #savefig("Fe5band.png")
    #return pls

    ham = hamiltonian_k(la)
    kx = 0.1
    ky = 0.2
    hamk = ham([kx,ky])
    #println(hamk)
    return hamk

end

function test_2DGraphene()
    a1 = [sqrt(3)/2,1/2]
    a2= [0,1]
    la2 = set_Lattice(2,[a1,a2])
    add_atoms!(la2,[1/3,1/3])
    add_atoms!(la2,[2/3,2/3])

    show_neighbors(la2)

    t = 1.0
    add_hoppings!(la2,-t,1,2,[1/3,1/3])
    add_hoppings!(la2,-t,1,2,[-2/3,1/3])
    add_hoppings!(la2,-t,1,2,[1/3,-2/3])

    klines = set_Klines()
    kmin = [-π]
    kmax = [π]
    add_Kpoints!(klines,kmin,kmax,"-pi","pi")

    pls = calc_band_plot_finite(klines,la2,1,periodic=false)
    #return pls

    pls = plot_lattice_2d(la2)
    #return pls

    klines = set_Klines()
    kmin = [0,0]
    kmax = [2π/sqrt(3),0]
    add_Kpoints!(klines,kmin,kmax,"G","K")

    kmin = [2π/sqrt(3),0]
    kmax = [2π/sqrt(3),2π/3]
    add_Kpoints!(klines,kmin,kmax,"K","M")

    kmin = [2π/sqrt(3),2π/3]
    kmax = [0,0]
    add_Kpoints!(klines,kmin,kmax,"M","G")

    println("Get DOS data...")
    hist = get_DOS(la2,100)
    println("Done.")
    #pls = plot_DOS(la2,100)
    #return pls

    pls2 = calc_band_plot(klines,la2)

    ham = hamiltonian_k(la2)
    kx = 0.1
    ky = 0.2
    hamk = ham([kx,ky])
    println(hamk)
    return hamk

end

function test_surface()
    Ax = 1
    Ay = 1
    m2x = 1
    m2y = m2x
    m0 = -2*m2x
    m(k) = m0 + 2m2x*(1-cos(k[1]))+2m2y*(1-cos(k[2]))
    Hk(k) = Ax*sin(k[1]).*σx +  Ay*sin(k[2]).*σy + m(k).*σz
    norb = 2 #The size of the matrix
    hamiltonian = surfaceHamiltonian(Hk,norb,numhop=3,L=32,kpara="kx",BC="OBC")

    nkx = 100
    kxs = range(-π,stop=π ,length=nkx)
    mat_e = zeros(Float64,nkx,32*2)
    for i=1:nkx
        kx = kxs[i]
        mat_h = hamiltonian(kx)
    #println(mat_h)

        e,v = eigen(Matrix(mat_h))
        #println(e)
        mat_e[i,:] = real.(e[:])
    end

    plot(kxs,mat_e,labels="")
    return mat_e

    #

    return
end


text:
using Twistronics

n_a1=90
n_a2=90
m=17          # in version 0.1.0, when m is even, it does not work: the way the I construct the lattice destroys bipartite nature of honeycomb (will fix in later versions)
r=1
r0=[0,0,0]   # inter-layer alignment
d=1          # inter-layer distance
lattice,mode="honeycomb","uniform_hop"
# m and r determine twist angle

# build the geometry, it contains all sites in the supercell
g=Twistronics.get_twisted_lattice_2(lattice,n_a1,n_a2,m,r,d,r0,mode)

#Twistronics.plot_R(g.sites) # plot the sites

t_xy=-1 # intra-layer hopping
t_z=0.3 # inter-layer hopping
lambda=10 # this parameter determines the decay of inter-layer hopping w.r.t distance

H=Twistronics.get_H_inter_twisted_bilayer(g,t_xy,t_z,lambda,mode)  # get the Hamiltonian of twisted bilayer
P=Twistronics.get_valley_operator_inter_twisted_bilayer(g,mode)    # get the valley operator of the twisted bilayer

#Twistronics.plot_H(H.intra,g.sites)              # Plot the intra-unitcell Hamiltonian to see if hoppings are correct

b1,b2=Twistronics.get_reciprocal_vector(g.inter_vector.x,g.inter_vector.y) # get reciprocal vector for the supercell
k_points=[i/100*(b1+2*b2) for i=1:100]  # plot bandstrcture along Gamma-K-K'-Gamma
#k_points=[i/100*(b1-b2) for i=1:100]
n_bands=40                              # compute 40 eigenvalues for each k point, using Arpack
momenta, energies, valley_polarization=Twistronics.get_band_twisted(H,P,k_points,n_bands) # get bandstructure with valley-polarization

# plot bandstructure
Twistronics.plt.figure(figsize=(6,6),dpi=80)
colors=[(1,0,0),(0,1,0),(0,0,1)]
cm=Twistronics.colormap.LinearSegmentedColormap.from_list("my_list", colors)
sc=Twistronics.plt.scatter(momenta, energies,c=valley_polarization,cmap=cm,vmin=-1,vmax=1)  # color indicates valley polarization
Twistronics.plt.colorbar(sc)
Twistronics.plt.axis([-1,100+1,-0.1,0.1])
Twistronics.plt.xlabel("k")
Twistronics.plt.ylabel("E")
Twistronics.plt.show()

# Note: for TBG, all low energy states are valley-degenerate

ext: jl


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: February 2021
# Copyright: Georgia Institute of Technology, 2020

module ReliabilityAwareImaging

using DrWatson
using Flux
using JOLI
using JLD2
using HDF5
using Random
using DataFrames
using LinearAlgebra
using Distributions
using Statistics
using InvertibleNetworks
using CUDA

import Base.getindex
import Distributions: logpdf, gradlogpdf

# Utilities
include("./utils/load_experiment.jl")
include("./utils/data_loader.jl")
include("./utils/savefig.jl")
include("./utils/logpdf.jl")

# Objective functions
include("./objectives/objectives.jl")

end


text:
# when string matching patterns, we always use "/" as the
# delimiter. good luck if your path system allows "/" in
# file names!
canonical_path(p::AbstractPath) = join(p.segments, "/")

struct NoValue end

"""
    FileTree(parent, name, children, value)

#### Fields:

- `parent::Union{FileTree, Nothing}` -- The parent node. `nothing` if it's the root node.
- `name::String` -- Name of the root.
- `children::Vector` -- children
- `value::Any` -- the value at the node, if no value is present, a `NoValue()` sentinal value.

    FileTree(tree::FileTree; parent, name, children, value)

Copy over all fields from `tree`, but use any fields provided as keyword arguments.

    FileTree(dirname::String)

Construct a `FileTree` to reflect directory from disk in the current working directory.
"""
struct FileTree
    parent::Union{FileTree, Nothing}
    name::String
    children::Vector
    value::Any
end

FileTree(parent, name, children) = FileTree(parent, name, children, NoValue())

function Base.isequal(f1::FileTree, f2::FileTree)
    !isequal(f1.name, f2.name) && return false

    c1 = issorted(f1.children, by=name) ? f1.children : sort(f1.children, by=name)
    c2 = issorted(f2.children, by=name) ? f2.children : sort(f2.children, by=name)

    !isequal(c1, c2) && return false
    isequal(f1.value, f2.value)
end

# convenience method to replace a few parameters
# and leave others unchanged
function FileTree(t::FileTree;
                  parent=parent(t),
                  name=t.name,
                  children=t.children,
                  value=t.value)
    FileTree(parent, name, children, value)
end

FileTree(dir) = FileTree(nothing, dir)

function FileTree(parent, dir)
    children = []
    parent′ = FileTree(parent, dir, children)

    ls = readdir(dir)
    cd(dir) do
        children′ = map(ls) do f
            if isdir(f)
                FileTree(parent′, f)
            else
                File(parent′, f)
            end
        end
        append!(children, children′)
    end

    parent′
end

"""
    children(node::Union{FileTree, File})

Get the immediate children of a `FileTree` node.
If node is `File` then returns `()`.
"""
children(d::FileTree) = d.children

"""
    parent(node::Union{FileTree, File})

Get the parent node. Returns `nothing` if there are no parents.
"""
parent(f::FileTree) = f.parent

"""
    name(node::Union{FileTree, File})

Get the file or directory name.
"""
name(f::FileTree) = f.name

Base.isempty(d::FileTree) = isempty(d.children)

Base.empty(d::FileTree) = FileTree(d; children=[])

"""
    rename(node::Union{FileTree, File}, newname)

Return a copy of node with `name` set to `newname`.
"""
function rename(x::FileTree, newname)
    setparent(FileTree(x, name=string(newname)))
end

"""
    setvalue(node::Union{FileTree, File}, val)

Return a copy of node with `val` set as the value.
"""
setvalue(x::FileTree, val) = FileTree(x; value=val)

"""
    setparent(node::Union{FileTree, File}, parent)

Return a copy of node with `parent` set as the parent.
"""
function setparent(x::FileTree, parent=parent(x))
    p = FileTree(x, parent=parent, children=copy(x.children))
    copy!(p.children, setparent.(x.children, (p,)))
    p
end

Base.show(io::IO, d::FileTree) = AbstractTrees.print_tree(io, d)

"""
    File(parent, name, value=NoValue())

#### Fields:

- `parent::Union{FileTree, Nothing}` -- The parent node. `nothing` if it's the root node.
- `name::String` -- Name of the root.
- `value::Any` -- the value at the node, if no value is present, a `NoValue()` sentinal value.
"""
struct File
    parent::Union{Nothing, FileTree}
    name::String
    value::Any
end

File(parent, name) = File(parent, name, NoValue())

function File(f::File; parent=parent(f), name=f.name, value=f.value)
    File(parent, name, value)
end

Base.show(io::IO, f::File) = print(io, string("File(", path(f), ")"))

function AbstractTrees.printnode(io::IO, f::Union{FileTree, File})
    print(io, name(f))
    if f isa FileTree
        print(io, "/")
    end
    if hasvalue(f)
        print(io," (", summary(f[]), ")")
    end
end

File(parent, name::String) = File(parent, name, NoValue())

function Base.isequal(f1::File, f2::File)
    f1.name == f2.name && isequal(f1.value, f2.value)
end

children(d::File) = ()

parent(f::File) = f.parent

name(f::File) = f.name

Base.isempty(d::File) = false

Base.empty(d::File) = d

setvalue(x::File, val) = File(x, value=val)

setparent(x::File, parent=parent(x)) = File(x, parent=parent)

Base.getindex(tree::FileTree, i::Int) = tree.children[i]

function Base.getindex(tree::FileTree, ix::Vector)
    FileTree(tree;
            children=vcat(map(i->(x=tree[i]; i isa Regex ? x.children :  x), ix)...))
end

Base.getindex(tree::FileTree, i::AbstractString) = tree[Path(i)]

function Base.getindex(tree::FileTree, i::AbstractPath)
    spath = i.segments

    if length(spath) > 1
        for s in spath
            tree = tree[s]
        end
        return tree
    end
    p = string(i)

    idx = findfirst(x->name(x)==p, children(tree))
    if idx === nothing
        error("No file matched getindex $i")
    end

    tree[idx]
end

function Base.getindex(tree::FileTree, subtree::FileTree; toplevel=true)
    if name(tree) != name(subtree)
        toplevel && error("Root name should match while getindexing with trees")
        return nothing
    else
        namesof_subtree = name.(children(subtree))
        cs = []
        for c in children(tree)
            idx = findfirst(==(name(c)), namesof_subtree)
            if !isnothing(idx)
                rest = c[subtree[idx]]
                if rest !== nothing
                    push!(cs, rest)
                end
            end
        end
        if isempty(cs)
            return nothing
        else
            return FileTree(tree; children=cs)
        end
    end
end

function Base.getindex(f::File, g::File; toplevel=true)
    if name(f) != name(g)
        if toplevel
            error("File not found")
        else
            return nothing
        end
    else
        return f
    end
end

rename(x::File, newname) = File(x, name=string(newname))

### Stuff agnostic to FileTree or File nature of "Node"s

const Node = Union{FileTree, File}

_maketree(node::String) = File(nothing, node, NoValue())
_maketree(node::NamedTuple) = File(nothing, node.name, node.value)
_maketree(node::Pair) = _maketree(node[1], node[2])
_maketree(node::Node) = node

function _maketree(node, children)
    cs = maketree.(children)
    if node isa Node
        return typeof(node)(node; children=cs)
    elseif node isa NamedTuple
        name = node.name
        value = node.value
    else
        name = node
        value = NoValue()
    end
    return FileTree(nothing, name, cs, value)
end

maketree(node) = setparent(_maketree(node))
maketree(node::Vector) = maketree("."=>node)

Base.basename(d::Node) = Path(d.name)

"""
    path(file::Union{File, FileTree)

Returns an [`AbstractPath`](https://rofinn.github.io/FilePathsBase.jl/stable/design/#Path-Types-1) object which is the path of the file from the
root node leading up to this file.
"""
function path(d::Node)
    parent(d) === nothing && return Path(d.name)
    path(parent(d)) / Path(d.name)
end

Base.dirname(d::Node) = dirname(path(d))

"""
    `node[]`

Get the value stored in the node. `NoValue()` is
returned if there is no value stored.
"""
Base.getindex(d::Node) = d.value

hasvalue(x::Node) = !(x[] isa NoValue)

## Tree walking

function prewalk(f, t::FileTree; collect_children=identity)
    x = f(t)
    if x isa FileTree
        cs = map(c->prewalk(f, c; collect_children=collect_children), t.children)
        FileTree(x; children=collect_children(cs))
    else
        return x
    end
end

prewalk(f, t::File; collect_children=identity) = f(t)

function postwalk(f, t::FileTree; collect_children=identity)
    cs = map(c->postwalk(f, c; collect_children=collect_children), t.children)
    f(FileTree(t; children=collect_children(cs)))
end

postwalk(f, t::File; collect_children=identity) = f(t)

"""
    map(f, tree::FileTree; walk=FileTrees.postwalk, dirs=true)

apply `f` to every node in the tree. To only visit File nodes, pass `dirs=false`.

walk can be either `FileTrees.postwalk` or `FileTrees.postwalk`.
"""
function Base.map(f, tree::FileTree; walk=postwalk, dirs=true)
    walk(tree, collect_children=identity) do n
        (dirs || n isa File) ? f(n) : n
    end
end

"""
    map(f, tree::FileTree; walk=FileTrees.postwalk, dirs=true)

apply `f` to every node in the tree. To only visit File nodes, pass `dirs=false`.

walk can be either `FileTrees.postwalk` or `FileTrees.postwalk`.

    filter(f, tree::FileTree; walk=FileTrees.postwalk, dirs=true)

remove every node `x` from `tree` where `f(x)` is `true`. `f(x)` must return a boolean value.
"""
function Base.filter(f, tree::FileTree; walk=prewalk, dirs=true)
    walk(tree, collect_children=cs->filter(!isnothing, cs)) do n
        (dirs || n isa File) ? (f(n) ? n : nothing) : n
    end
end


text:
srcdir = dirname(pwd())*"/src"
any(srcdir .== LOAD_PATH) || push!(LOAD_PATH, srcdir)
push!(LOAD_PATH, pwd())
using PyPlot, Statistics, Distributed, Revise, Tools3D, TransD_GP, MCMC_Driver
##
n, dz, extendfrac  = 50, 1.01, 1.1
z, znall = Tools3D.makezn(n=n, dz=dz, extendfrac=extendfrac)

opt, optdummy, x, y = Tools3D.makeopt(znall=znall, dxfrac = 0.01, dyfrac=0.01)
m = TransD_GP.init(opt)
mdummy = TransD_GP.init(optdummy, m)
# Tools3D.slicemodel(m, opt, slicesx=[], slicesy=[], slicesz=[16,19], dz=dz, extendfrac=extendfrac)
#
# ρ = Tools3D.makecubemodel(opt, dz=dz, extendfrac=extendfrac, ρanom=0.2)
# slicesx, slicesy, slicesz = [20], [20], [17]
# Tools3D.slicemodel(ρ,0,[0. 0.],[0.], opt, slicesx=slicesx, slicesy=slicesy, slicesz=slicesz, dz=dz, extendfrac=extendfrac)
# savefig("true_model.png",dpi=300)
#
# d, sd = Tools3D.get_training_data(ρ, opt, dz=dz, extendfrac=extendfrac, zbreak=2.0, fractrain=0.05)
# savefig("random_data.png")
# Tools3D.calc_simple_RMS(d, ρ, sd)
## timing
for i = 1:98 # 20 also works well if you switch off timebirth. Keep it on for an interesting model
    TransD_GP.birth!(m, opt, mdummy, optdummy)
end
timebirth = true
if timebirth
    ts = time(); ndo = 100
    for i = 1:ndo
        TransD_GP.birth!(m, opt, mdummy, optdummy)
        TransD_GP.death!(m, opt, mdummy, optdummy)
    end
    dt = time() - ts
    @info "avg time per move is $(dt/ndo/2)"
end


text:
function load_from_gfa(filename::String)

    graph = SequenceGraph(SequenceGraphNode[], SequenceGraphLink[])
    oldnames_to_ids = Dict{String, NodeID}()
    oldnames = Vector{String}()

    gfaf = open(filename, "r")
    if eof(gfaf)
        throw(ArgumentError("GFA file is empty"))
    end
    line = readline(gfaf)
    if line != "H\tVN:Z:1.0"
        @warn "The first line of the GFA does not correspond to GFA1"
    end
    faf = open(FASTA.Reader, string(splitext(filename), ".fasta"))
    if eof(faf)
        throw(ArgumentError("FASTA file is empty"))
    end

    # Load all sequences from fasta, if they're not canonical, flip them and
    # remember they are flipped.

    rcnodes = 0

    farec = FASTA.Record()
    while !eof(faf)
        read!(faf, farec)
        nodename = FASTA.identifier(farec)
        nodeseq = FASTA.sequence(BioSequence{DNAAlphabet{2}}, farec)
        if haskey(oldnames_to_ids, nodename)
            error("sequence ", nodename, " is already defined.")
        end
        oldnames_to_ids[nodename] = add_node!(graph, Node(seq, true));
        push!(oldnames, nodename)
        # Canonicalize the sequence
        if !iscanonical(nodeseq)
            reverse_complement!(nodeseq)
            oldnames_to_ids[nodename] = -oldnames_to_ids[nodename]
            rcnodes += 1
        end
    end
    @info string(length(nodes(graph)), " nodes loaded (", rcnodes, " canonised).")

    dist_egt0 = 0
    lcount = 0
    # Load and store all connections
    for line in eachline(gfaf)
        tokens = split(line)
        recordtype = tokens[1]

        if recordtype == "S"
            id = tokens[2]
            seq = tokens[3]
            len = tokens[4]
            if seq != '*'
                error("sequences should be in a seperate fasta file.")
            end
            # Check equal length seq and node length reported in gfa
            if haskey(oldnames_to_ids, id)
                if parse(Int, len[6:end]) != length(node(graph, oldnames_to_ids[id]))
                    error(string("Different length in node and fasta for sequence: ", id, " -> gfa: ",  parse(Int, len[6:end]),  ", fasta: ", length(node(graph, oldnames_to_ids[id]))))
                end
            end
        elseif recordtype == "L"
            source = tokens[2]
            sdir = tokens[3]
            dest = tokens[4]
            ddir = tokens[5]
            cigar = tokens[6]
            if haskey(oldnames_to_ids, source)
                oldnames_to_ids[source] = add_node!(graph, Node(BioSequences{DNAAlphabet{2}}()))
                @info string("added source ", source)
            end
            if haskey(oldnames_to_ids, dest)
                oldnames_to_ids[dest] = add_node!(graph, Node(BioSequences{DNAAlphabet{2}}()))
                @info string("added dest ", dest)
            end
            source_id = oldnames_to_ids[source]
            dest_id = oldnames_to_ids[dest]
            # Go from GFA's "Links as paths" to a normal "nodes have sinks (-/start/left) and sources (+/end/right)"
            if sdir == "+"
                source_id = -source_id
            end
            if ddir == "-"
                dest_id = -dest_id
            end
            dist = 0
            if length(cigar) > 1 && cigar[end] == 'M'
                # TODO: better checks for M
                dist = -parse(Int, cigar[1:end-1])
                if dist >= 0
                    dist_egt0 += 1
                end
            end
            add_link!(graph, source_id, dest_id, dist)
            lcount += 1
        end
    end
    close(gfaf)
    if dist_egt0 > (lcount / 2)
        @warn string("Warning: The loaded graph contains ", dist_egt0, " non-overlapping links out of ", lcount)
    end
    @info string(length(nodes(graph)) - 1, " nodes after connecting with ", lcount, "links.")
end


text:
struct AlwaysAssertionError <: Exception
    msg::String
end

struct NewPackage end
struct NewVersion end

abstract type AutoMergeException <: Exception end

struct AutoMergeAuthorNotAuthorized <: AutoMergeException
    msg::String
end

struct AutoMergeCronJobError <: AutoMergeException
    msg::String
end

struct AutoMergeGuidelinesNotMet <: AutoMergeException
    msg::String
end

struct AutoMergeNeitherNewPackageNorNewVersion <: AutoMergeException
    msg::String
end

struct AutoMergePullRequestNotOpen <: AutoMergeException
    msg::String
end

struct AutoMergeShaMismatch <: AutoMergeException
    msg::String
end

struct AutoMergeWrongBuildType <: AutoMergeException
    msg::String
end

struct GitHubAutoMergeData
    # Handle to the GitHub API. Used to query the PR and update
    # comments and status.
    api::GitHub.GitHubAPI

    # Whether the registry PR refers to a new package or a new version
    # of an existing package.
    registration_type::Union{NewPackage,NewVersion}

    # The GitHub pull request data.
    pr::GitHub.PullRequest

    # Name of the package being registered.
    pkg::String

    # Version of the package being registered.
    version::VersionNumber

    # Used for updating CI status.
    current_pr_head_commit_sha::String

    # The GitHub repo data for the registry.
    registry::GitHub.Repo

    # GitHub authorization data.
    auth::GitHub.Authorization

    # Type of authorization for automerge. This can be either:
    # :jll - special jll exceptions are allowed,
    # :normal - normal automerge rules.
    authorization::Symbol

    # Directory of a registry clone that includes the PR.
    registry_head::String

    # Directory of a registry clone that excludes the PR.
    registry_master::String

    # Whether to add a comment suggesting bumping package version to
    # 1.0 if appropriate.
    suggest_onepointzero::Bool

    # GitHub identity resulting from the use of an authentication token.
    whoami::String

    # List of dependent registries. Typically this would contain
    # "General" when running automerge for a private registry.
    registry_deps::Vector{String}

    # Location of the directory where the package code
    # will be downloaded into. Populated at construction time
    # via `mktempdir`.
    pkg_code_path::String

    # A list of public Julia registries (repository URLs) which will
    # be checked for UUID collisions in order to mitigate the
    # dependency confusion vulnerability. See the
    # `dependency_confusion.jl` file for details.
    public_registries::Vector{String}

    # whether only read-only actions should be taken
    read_only::Bool
end

# Constructor that requires all fields (except `pkg_code_path`) as named arguments.
function GitHubAutoMergeData(; kwargs...)
    pkg_code_path = mktempdir(; cleanup=true)
    kwargs = (; pkg_code_path=pkg_code_path, kwargs...)
    fields = fieldnames(GitHubAutoMergeData)
    always_assert(Set(keys(kwargs)) == Set(fields))
    always_assert(kwargs[:authorization] ∈ (:normal, :jll))
    return GitHubAutoMergeData(getindex.(Ref(kwargs), fields)...)
end

Base.@kwdef mutable struct Guideline
    # Short description of the guideline. Only used for logging.
    info::String

    # Documentation for the guideline
    docs::Union{String,Nothing} = info

    # Function that is run in order to determine whether a guideline
    # is met. Input is an instance of `GitHubAutoMergeData` and output
    # is passed status plus a user facing message explaining the
    # guideline result.
    check::Function

    # Saved result of the `check` function.
    passed::Bool = false

    # Saved output message from the `check` function.
    message::String = "Internal error. A check that was supposed to run never did: $(info)"
end

passed(guideline::Guideline) = guideline.passed
message(guideline::Guideline) = guideline.message
function check!(guideline::Guideline, data::GitHubAutoMergeData)
    return guideline.passed, guideline.message = guideline.check(data)
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

using Base.Test

let
debug = false

# Test givens rotations
for elty in (Float32, Float64, Complex64, Complex128)
    debug && println("elty is $elty")

    if elty <: Real
        A = convert(Matrix{elty}, randn(10,10))
    else
        A = convert(Matrix{elty}, complex.(randn(10,10),randn(10,10)))
    end
    for Atype in ("Array", "SubArray")
        if Atype == "Array"
            A = A
        else
            A = view(A, 1:10, 1:10)
        end
        Ac = copy(A)
        R = Base.LinAlg.Rotation(Base.LinAlg.Givens{elty}[])
        for j = 1:8
            for i = j+2:10
                G, _ = givens(A, j+1, i, j)
                A_mul_B!(G, A)
                A_mul_Bc!(A, G)
                A_mul_B!(G, R)

                @test A_mul_B!(G,eye(elty,10,10)) == [G[i,j] for i=1:10,j=1:10]

                # test transposes
                @test ctranspose(G)*G*eye(10) ≈ eye(elty, 10)
                @test ctranspose(R)*(R*eye(10)) ≈ eye(elty, 10)
                @test_throws ErrorException transpose(G)
                @test_throws ErrorException transpose(R)
            end
        end
        @test_throws ArgumentError givens(A, 3, 3, 2)
        @test_throws ArgumentError givens(one(elty),zero(elty),2,2)
        G, _ = givens(one(elty),zero(elty),11,12)
        @test_throws DimensionMismatch A_mul_B!(G, A)
        @test_throws DimensionMismatch A_mul_Bc!(A,G)
        @test abs.(A) ≈ abs.(hessfact(Ac)[:H])
        @test norm(R*eye(elty, 10)) ≈ one(elty)

        G, _ = givens(one(elty),zero(elty),9,10)
        @test ctranspose(G*eye(elty,10))*(G*eye(elty,10)) ≈ eye(elty, 10)
        K, _ = givens(zero(elty),one(elty),9,10)
        @test ctranspose(K*eye(elty,10))*(K*eye(elty,10)) ≈ eye(elty, 10)

        # test that Givens * work for vectors
        if Atype == "Array"
            x = A[:, 1]
        else
            x = view(A, 1:10, 1)
        end
        G, r = givens(x[2], x[4], 2, 4)
        @test (G*x)[2] ≈ r
        @test abs((G*x)[4]) < eps(real(elty))
        @inferred givens(x[2], x[4], 2, 4)

        G, r = givens(x, 2, 4)
        @test (G*x)[2] ≈ r
        @test abs((G*x)[4]) < eps(real(elty))
        @inferred givens(x, 2, 4)

        G, r = givens(x, 4, 2)
        @test (G*x)[4] ≈ r
        @test abs((G*x)[2]) < eps(real(elty))
    end
end
end #let


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Usrsctp_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Usrsctp")
JLLWrappers.@generate_main_file("Usrsctp", UUID("9f05be0d-5fa0-5d85-a018-f41a39aa9d88"))
end  # module Usrsctp_jll


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

# recipe for base case
@recipe function f(segment::Segment)
  seriestype --> :scatterpath
  seriescolor --> :black
  primary --> false

  Tuple.(coordinates.(vertices(segment)))
end

# plot facets in a recursion
@recipe function f(polytope::Polytope)
  for f in facets(polytope)
    @series begin
      primary --> false
      f
    end
  end
end


text:
module TestModule

using Test
using MLStyle

liftsym = Symbol("@lift")

function lift!(ex, lifted::Vector{Any})
    @switch ex begin
        @case Expr(:macrocall, &liftsym, _, arg)
        push!(lifted, arg)
        ex.args
        ex.args[1] = Symbol("@static")
        ex.args[3] = :(true ? true : true)
        return
        @case Expr(hd, args...)
        for arg in args
            lift!(arg, lifted)
        end
        return
        @case _
        return
    end
end

macro testcase(name, ex)
    lifted = []
    lift!(ex, lifted)
    m = gensym(name)
    __module__.eval(:(module $m
    using MLStyle
    using Test
    MODULE = $MODULE
    $(Symbol("@test_macro_throws")) = $(getfield(TestModule, Symbol("@test_macro_throws")))
    $(lifted...)
    @testset $name $ex
    end))
end

macro test_macro_throws(errortype, m)
    :(@test_throws $errortype try
        @eval $m
    catch err
        while err isa LoadError
            err = err.error
        end
        throw(err)
    end)
end

MODULE = TestModule

@use GADT

include("when.jl")
include("untyped_lam.jl")
include("active_patterns.jl")
include("uncomp.jl")
include("lambda.jl")
include("as_record.jl")
include("adt.jl")

include("exception.jl")
include("expr_template.jl")
include("gallery/simple.jl")
include("dot_expression.jl")

include("modules/cond.jl")
include("modules/ast.jl")

include("pervasive.jl")
include("match.jl")
include("pattern.jl")
include("typelevel.jl")
include("nothing.jl")

include("MQuery/test.jl")

include("issues/87.jl")

end


text:
using TOML
using UUIDs
using GitHub

Base.@kwdef struct VersionTokens
    major::String = "major"
    minor::String = "minor"
    patch::String = "patch"
end

const VERSION_TOKENS = VersionTokens()
Base.show(io::IO, vt::VersionTokens) = print(io, "(", vt.major, ", ", vt.minor, ", ", vt.patch, ")")
Base.in(version::String, tokens::VersionTokens) = (version == tokens.major) ||
    (version == tokens.minor) || (version == tokens.patch)

function is_version_number(version)
    occursin(r"[0-9]+.[0-9]+.[0-9]+", version) ||
        occursin(r"v[0-9]+.[0-9]+.[0-9]+", version)
end

struct ProjectToRelease
    name::String # package name
    path::String # path to the git repo
    subdir::String # subdir to the package
    toml::String # Project.toml
    uuid::UUID # project UUID
    version::Maybe{VersionNumber} # current version
    branch::String # branch to release

    # user inputs
    note::String
    version_spec::String
    release_version::VersionNumber
end

function ProjectToRelease(version_spec::String, path_to_project::String=pwd();
        branch::String=default_branch(path_to_project), note::String="")

    toml = Base.current_project(path_to_project)
    toml === nothing && cmd_error("cannot find (Julia)Project.toml in $path_to_project")
    path_to_project = dirname(toml)
    path_to_repo = git_get_toplevel_path(path_to_project)
    pkg = Pkg.Types.read_project(toml)
    subdir = relpath(path_to_project, path_to_repo)

    return ProjectToRelease(
        pkg.name, path_to_repo, subdir, abspath(toml),
        pkg.uuid, pkg.version, branch,
        note, version_spec,
        release_version(pkg.name, pkg.version, version_spec)
    )
end

Base.show(io::IO, ::MIME"text/plain", x::ProjectToRelease) = GarishPrint.pprint_struct(io, x)

# NOTE: About Registration Process
#
# 1. checkout to the release branch
# 2.

"""
release a package.

# Intro

This command can release packages and auto bump versions for you
semantically. One will not need to remember how to type `JuliaRegistrator`
correctly but just use `ion release patch` under the package folder. Or
via the REPL interface as `Ion.release("patch")`

# Arguments

- `version_spec`: version number you want to release. Can be a specific version, "current"
    or either of $(VERSION_TOKENS)
- `path`: path to the project you want to release.

# Options

- `-r,--registry <registry name>`: registry you want to register the package.
    If the package has not been registered, ion will try to register
    the package in the General registry. Or the user needs to specify
    the registry to register using this option.
- `-b, --branch <branch name>`: branch you want to register.
- `--note <release note>`: optional, release note you would like to specify.
"""
@cast function release(version_spec::String, path::String=pwd();
        registry::String="",
        branch::String=current_branch(path),
        note::String="", debug::Bool=false,
    )
    project = ProjectToRelease(version_spec, path; branch, note)
    checkout(project.path; project.branch, quiet=true) do
        check_if_repo_dirty(project)
        check_release_info(project)
        sync_with_remote(project)
        try_register(registry, project)
    end
    return
end

function check_if_repo_dirty(project::ProjectToRelease)
    isdirty(project.path) &&
        cmd_error("package repository is dirty, please commit or stash changes.")
    return
end

function check_release_info(project::ProjectToRelease)
    if should_update_version(project)
        update_version!(project)
        commit_version_update(project)
    elseif prompt() do
            print("do you want to release current version ")
            printstyled(project.version; color=:light_cyan)
            print(" ?")
        end
    else
        cmd_exit()
    end
    return
end

function sync_with_remote(project::ProjectToRelease)
    @info "push local changes to remote"
    cd(project.path) do
        redirect_stdio(;stdout=devnull, stderr=devnull) do
            run(`$(git()) push origin $(project.branch)`)
        end
    end
    return
end

function try_register(registry::String, project::ProjectToRelease)
    try
        register(registry, project)
    catch e
        if should_update_version(project)
            run(`$(git()) revert --no-edit --no-commit HEAD`)
            run(`$(git()) commit -m "revert version bump due to an error occured in IonCLI"`)
            sync_with_remote(project)
            printstyled(" ✔  "; color=:light_green)
            println("revert version bump commit:")
        end
    end
    return
end

# ion release patch lib/EaRydCore
# ion release patch
function release_version(name::String, version::VersionNumber, version_spec::String)
    version_spec == "current" && return version
    Pkg.Registry.update() # make sure we always get latest version
    latest_version = find_max_version(name)

    if is_version_number(version_spec)
        version_number = VersionNumber(version_spec)
    elseif version_spec in VERSION_TOKENS
        version_number = bump_version(latest_version, version_spec)
    else
        cmd_error("invalid version spec: $(version_spec)")
    end

    if latest_version === nothing
        println("package not found in local registries")
    end

    print(" "^7, "version to release: ")
    printstyled(version_number; color=:light_cyan)
    println()

    if version == latest_version
        # current is same as latest version
        print(" "^7, "current/latest version: ")
        printstyled(latest_version; color=:light_cyan)
        println()
    elseif version > latest_version
        # current is newer than latest version
        # this may because we had a failed manual
        # version bump before
        println("your current version is already larger than latest release")

        print(" "^10, "current version: ")
        printstyled(version; color=:light_cyan)
        println()

        print(" "^10, "latest version: ")
        printstyled(latest_version; color=:light_cyan)
        println()

        if is_version_continuously_greater(latest_version, version)
            # very likely the current version is the correct one
            #
            # we must know explicitly what user intended to
            # even we have the specific version spec since
            # this is very likely to be rejected in General
            if prompt("do you want to release current version instead?"; require=true)
                version_number = version
            end
        else
            cmd_error("cannot release discontinous versions, please check your version in Project.toml")
        end
    else
        # current is older than latest version
        # this may because one wants to release
        # a patch/minor version of previous major
        # version, e.g
        # old versions: 1.1.0, 1.2.0
        # wants to release: 1.1.1
        println("your current version is smaller than latest release")

        print(" "^10, "current version: ")
        printstyled(version; color=:light_cyan)
        println()

        print(" "^10, "latest version: ")
        printstyled(latest_version; color=:light_cyan)
        println()

        prompt("are you sure to release this version?"; require=true) || cmd_exit()
    end
    return version_number
end

function should_update_version(project::ProjectToRelease)
    project.version_spec == "current" && return false
    project.release_version == project.version && return false
    return true
end

function update_version!(project::ProjectToRelease)
    prompt("do you want to update Project.toml?") || cmd_exit()
    write_version(project)
    print(" ")
    printstyled("✔"; color=:light_green)
    print("  Project.toml has been updated to ")
    printstyled(project.release_version; color=:light_cyan)
    println()
    return project
end

function is_version_continuously_greater(latest::VersionNumber, release::VersionNumber)
    release > latest || return false
    # patch release
    latest.major == release.major && latest.minor == release.minor &&
        latest.patch+1 == release.patch && return true

    # minor release
    latest.major == release.major && latest.minor+1 == release.minor &&
        latest.patch == release.patch && return true

    # major release
    latest.major+1 == release.major && latest.minor == release.minor &&
        latest.patch == release.patch && return true
    return false
end

function bump_version(::Nothing, token::String)
    return bump_version(v"0.0.0", token)
end

function bump_version(version::VersionNumber, token::String)
    if token == VERSION_TOKENS.major
        return VersionNumber(version.major+1, 0, 0)
    elseif token == VERSION_TOKENS.minor
        return VersionNumber(version.major, version.minor+1, 0)
    elseif token == VERSION_TOKENS.patch
        return VersionNumber(version.major, version.minor, version.patch+1)
    else
        cmd_error("invalid version spec $token")
    end
end

function write_version(project::ProjectToRelease)
    d = TOML.parsefile(project.toml)
    d["version"] = string(project.release_version)
    write_project_toml(project.toml, d)
    return
end

function commit_version_update(project::ProjectToRelease)
    @info "commit version updates" project
    cd(project.path) do
        redirect_stdout(devnull) do
            run(`$(git()) add $(project.toml)`)
            run(`$(git()) commit -m"bump version to $(project.release_version)"`)
        end
    end
end

"""
    PRN{name}

Package Registry Name
"""
struct PRN{name} end

"""
    PRN(name::String)

Create a `PRN` (Pacakge Registry Name) object.
"""
PRN(name::String) = PRN{Symbol(name)}()

macro PRN_str(name::String)
    return PRN{Symbol(name)}
end

Base.show(io::IO, ::PRN{registry}) where {registry} = print(io, "Pacakge Registry ", string(registry))

function register(registry::String, project::ProjectToRelease)
    registry = isempty(registry) ? determinte_register(project) : registry
    return register(PRN(registry), project)
end

function determinte_register(project::ProjectToRelease)
    info = find_package(project.name)
    if info === nothing
        if prompt(
                "package is not registered in local registries, " *
                "do you want to register it in the General registry?"
            )
            return "General"
        else
            cmd_exit()
        end
    else
        return info.reg.name
    end
end

function register(registry::PRN, project::ProjectToRelease)
    cmd_error("register workflow is not defined for $registry")
end

function register(::PRN"General", project::ProjectToRelease)
    @info "register package in General"
    auth_done = Base.Event()
    summon_done = Base.Event()
    interrupted_or_done = Base.Event()
    print_task = useless_animation(auth_done, summon_done, interrupted_or_done)

    @info "authenticating GitHub account"
    github_token = read_github_auth()
    auth = GitHub.authenticate(github_token)
    notify(auth_done)

    HEAD = get_head_sha256(project.path; project.branch)
    comment_json = Dict{String, Any}(
        "body" => julia_registrator_comment(project),
    )

    repo = github_repo(project.path)
    if repo === nothing
        cmd_error(
            "not a GitHub repository, registering " *
            "non-github packages is not supported, " *
            "please use JuliaHub instead"
        )
    end

    comment = GitHub.create_comment(repo, HEAD, :commit; params=comment_json, auth=auth)
    notify(summon_done)
    notify(interrupted_or_done)
    wait(print_task)

    print("\e[1G ")
    printstyled("✔"; color=:light_green)
    println("  JuliaRegistrator has been summoned, check it in the following URL:")
    printstyled("  ", comment.html_url; color=:cyan)
    println()
    return comment
end

function julia_registrator_comment(project::ProjectToRelease)
    lines = [
        "Released via [Ion](https://github.com/Roger-luo/Ion.jl)",
        "",
    ]
    cmd = ["@JuliaRegistrator", "register", "branch=$(project.branch)"]
    if project.subdir != "."
        push!(cmd, "subdir=$(project.subdir)")
    end
    push!(lines, join(cmd, " "))

    if !isempty(project.note)
        push!(lines, "")
        title = project.version_spec == "patch" ? "Patch notes:" : "Release notes:"
        push!(lines, title)
        push!(lines, project.note)
    end
    push!(lines, "")
    return join(lines, '\n')
end

function useless_animation(auth::Base.Event, summon::Base.Event, interrupted_or_done::Base.Event)
    anim_chars = ["◐","◓","◑","◒"]
    ansi_enablecursor = "\e[?25h"
    ansi_disablecursor = "\e[?25l"
    t = Timer(0; interval=1/20)
    print_lock = ReentrantLock()
    printloop_should_exit = interrupted_or_done.set
    return @async begin
        try
            count = 1
            while !printloop_should_exit
                lock(print_lock) do
                    print(ansi_disablecursor)
                    print("\e[1G  ")
                    printstyled(anim_chars[mod1(count, 4)]; color=:cyan)
                    print("    ")
                    if !auth.set
                        print("authenticating...")
                    elseif !summon.set
                        print("summoning JuliaRegistrator...")
                    else
                        printloop_should_exit = true
                    end
                end
                printloop_should_exit = interrupted_or_done.set
                count += 1
                wait(t)
            end
        catch e
            notify(interrupted_or_done)
            lock(print_lock) do
                println("\e[1G  ", RED_FG("❌"), "  fail to register $(project.pkg.name), error msg:")
                println(e.msg)
            end
            rethrow(e)
        finally
            print(ansi_enablecursor)
        end
    end
end


text:
import Base: ==

function ==(r1::RVRP.Range, r2::RVRP.Range)
    return (
        r1.lb == r2.lb
        && r1.ub == r2.ub
    )
end

function ==(f1::RVRP.Flexibility, f2::RVRP.Flexibility)
    return (
        f1.flexibility_level == f2.flexibility_level
        && f1.fixed_price == f2.fixed_price
        && f1.unit_price == f2.unit_price
    )
end

function ==(fr1::RVRP.FlexibleRange, fr2::RVRP.FlexibleRange)
    return (
        fr1.soft_range == fr2.soft_range
        && fr1.hard_range == fr2.hard_range
        && fr1.lb_flex == fr2.lb_flex
        && fr1.ub_flex == fr2.ub_flex
    )
end

function ==(l1::RVRP.Location, l2::RVRP.Location)
    return (
        l1.id == l2.id
        && l1.index == l2.index
        && l1.latitude == l2.latitude
        && l1.longitude == l2.longitude
        && l1.opening_time_windows == l2.opening_time_windows
        && l1.energy_fixed_cost == l2.energy_fixed_cost
        && l1.energy_unit_cost == l2.energy_unit_cost
        && l1.energy_recharging_speeds == l2.energy_recharging_speeds
    )
end

function ==(lg1::RVRP.LocationGroup, lg2::RVRP.LocationGroup)
    return (
        lg1.id == lg2.id
        && lg1.location_ids == lg2.location_ids
    )
end

function ==(pc1::RVRP.ProductCompatibilityClass, pc2::RVRP.ProductCompatibilityClass)
    return (
        pc1.id == pc2.id
        && pc1.conflict_compatib_class_ids == pc2.conflict_compatib_class_ids
        && pc1.prohibited_predecessor_compatib_class_ids == pc2.prohibited_predecessor_compatib_class_ids
    )
end

function ==(pcc1::RVRP.ProductCompatibilityClass, pcc2::RVRP.ProductCompatibilityClass)
    return (
        pcc1.id == pcc2.id
        && pcc1.pickup_availabitilies_at_location_ids == pcc2.pickup_availabitilies_at_location_ids
        && pcc1.delivery_capacities_at_location_ids == pcc2.delivery_capacities_at_location_ids
    )
end

function ==(ps1::RVRP.ProductSpecificationClass, ps2::RVRP.ProductSpecificationClass)
    return (
        sp1.id == sp2.id
        && ps1.capacity_consumptions == ps2.capacity_consumptions
        && ps1.property_requirements == ps2.property_requirements
    )
end

function ==(req1::RVRP.Request, req2::RVRP.Request)
    return (
        req1.id == req2.id
        && req1.product_compatibility_class_id == req2.product_compatibility_class_id
        && req1.product_sharing_class_id == req2.product_sharing_class_id
        && req1.product_specification_class_id == req2.product_specification_class_id
        && req1.split_fulfillment == req2.split_fulfillment
        && req1.request_flexibility == req2.request_flexibility
        && req1.precedence_status == req2.precedence_status
        && req1.product_quantity_range == req2.product_quantity_range
        && req1.pickup_location_group_id == req2.pickup_location_group_id
        && req1.delivery_location_group_id == req2.delivery_location_group_id
        && req1.pickup_service_time == req2.pickup_service_time
        && req1.delivery_service_time == req2.delivery_service_time
        && req1.max_duration == req2.max_duration
        && req1.duration_unit_cost == req2.duration_unit_cost
        && req1.pickup_time_windows == req2.pickup_time_windows
        && req1.delivery_time_windows == req2.delivery_time_windows
    )
end

function ==(vc1::RVRP.VehicleCategory, vc2::RVRP.VehicleCategory)
    return (
        vc1.id == vc2.id
        && vc1.capacities == vc2.capacities
        && vc1.properties == vc2.properties
        && vc1.loading_option == vc2.loading_option
        && vc1.energy_interval_lengths == vc2.energy_interval_lengths
    )
end

function ==(vc1::RVRP.CostPeriod, vc2::RVRP.CostPeriod)
    return (
        && vc1.work_period == vc2.work_period
        && vc1.travel_distance_unit_cost == vc2.travel_distance_unit_cost
        && vc1.travel_time_unit_cost == vc2.travel_time_unit_cost
        && vc1.service_time_unit_cost == vc2.service_time_unit_cost
        && vc1.waiting_time_unit_cost == vc2.waiting_time_unit_cost
        && vc1.fixed_cost == vc2.fixed_cost
    )
end

function ==(hvs1::RVRP.HomogeneousVehicleSet, hvs2::RVRP.HomogeneousVehicleSet)
    return (
        hvs1.id == hvs2.id
        && hvs1.vehicle_category_id == hvs2.vehicle_category_id
        && hvs1.cost_periods == hvs2.cost_periods
        && hvs1.work_periods == hvs2.work_periods
        && hvs1.departure_location_group_id == hvs2.departure_location_group_id
        && hvs1.arrival_location_group_id == hvs2.arrival_location_group_id
        && hvs1.initial_energy_charge == hvs2.initial_energy_charge
        && hvs1.max_working_time == hvs2.max_working_time
        && hvs1.max_travel_distance == hvs2.max_travel_distance
        && hvs1.allow_shipment_over_multiple_work_periods == hvs2.allow_shipment_on_multiple_work_periods
        && hvs1.nb_of_vehicles_range == hvs2.nb_of_vehicles_range
    )
end

function ==(data1::RVRP.RvrpInstance, data2::RVRP.RvrpInstance)
    return (
        data1.id == data2.id
        && data1.coordinate_mode == data2.coordinate_mode
        && data1.distance_mode == data2.distance_mode
        && data1.travel_matrix_periods == data2.travel_matrix_periods
        && data1.period_to_matrix_id == data2.period_to_matrix_id
        && data1.travel_time_matrices == data2.travel_time_matrices
        && data1.travel_distance_matrices == data2.travel_distance_matrices
        && data1.energy_consumption_matrices == data2.energy_consumption_matrices
        && data1.locations == data2.locations
        && data1.location_groups == data2.location_groups
        && data1.product_compatibility_classes == data2.product_compatibility_classes
        && data1.product_sharing_classes == data2.product_sharing_classes
        && data1.product_specification_classes == data2.product_specification_classes
        && data1.requests == data2.requests
        && data1.vehicle_categories == data2.vehicle_categories
        && data1.vehicle_sets == data2.vehicle_sets
    )
end


text:
#== What I want from an End to end function:
3. Save the experimental parameters (as a serializable object?)
4. For each component, save the time and a sparse rep. of the result
5. An easy to plot visualization, separate from these functions
==#

"""
    safety_based_refinement

Performs k-step refinement over uncertain states according to safety verification.
"""
function safety_based_refinement(experiment_params::ExperimentParameters, results_path::String, verification_mat, refinement_steps::Int;
                                 safety_threshold=0.95, horizon=1, reuse_regions_flag=false, reuse_transition_mats_flag=false)

    working_dir = results_path
    refinement_result_dirs = []

    for i=1:refinement_steps
        results_dir = "$results_path/refinement$i"
        states_to_update = findall(x->x<safety_threshold, verification_mat[1:end,3]) ∩ findall(x->x>safety_threshold, verification_mat[1:end,4])
        imdp = refine_imdp(states_to_update, experiment_params, working_dir=working_dir, results_dir=results_dir, reuse_regions_flag=reuse_regions_flag, reuse_transition_mats_flag=reuse_transition_mats_flag)
        verification_mat = Globally(imdp, "safe", horizon, "$results_dir/imdp-r$i.txt")
        save_legacy_mats(verification_mat, results_dir, horizon)
        working_dir = results_dir
        push!(refinement_result_dirs, results_dir)
    end

    return refinement_result_dirs
end

"""
    safety_based_synthesis

Performs synthesis on a system with a focus on remaining in a safe set.
"""
function safety_based_synthesis(experiment_params_array, system_paths::Array{String,1}, results_path::String,
    system_tag::String, label_fcn, refinement_steps::Int;
    minimum_threshold=0.80, horizon=-1, reuse_regions_flag=false, reuse_transition_mats_flag=false)

    exp_dir = results_path

    synth_dir = "$results_path/safety"
    synthesis_mat, results_path, _ = perform_synthesis_from_result_dirs_safety(system_paths, synth_dir, system_tag, label_fcn);
    save_legacy_mats(synthesis_mat, synth_dir, horizon)

    refinement_result_dirs = []

    for i=1:refinement_steps
        results_dir = "$synth_dir/refinement$i"
        states_to_update = findall(x->x<minimum_threshold, synthesis_mat[1:end,3]) ∩ findall(x->x>minimum_threshold, synthesis_mat[1:end,4])

        new_imdp_dirs = Array{String,1}()
        # For each IMDP in the system
        for (j, component_dir) in enumerate(system_paths)
            imdp_results_dir = "$results_dir/mode$j"
            refine_imdp(states_to_update, experiment_params_array[j], working_dir=component_dir, results_dir=imdp_results_dir, reuse_regions_flag=reuse_regions_flag, reuse_transition_mats_flag=reuse_transition_mats_flag)
            push!(new_imdp_dirs, imdp_results_dir)
        end

        cp("$results_dir/mode1/regions.bson", "$results_dir/regions.bson", force=true)

        synthesis_mat, _, _ = perform_synthesis_from_result_dirs_safety(new_imdp_dirs, results_dir, system_tag, label_fcn; horizon=horizon);

        save_legacy_mats(synthesis_mat, results_dir, horizon)

        system_paths = new_imdp_dirs
        push!(refinement_result_dirs, results_dir)
    end

    return refinement_result_dirs
end

"""
Performs k-step refinement over uncertain states according to minimum spec satisfaction from synthesis.
"""
function spec_based_refinement(experiment_params_array, system_paths::Array{String,1}, results_path::String,
                               system_tag::String, label_fcn, refinement_steps::Int;
                               minimum_threshold=0.80, horizon=-1, reuse_regions_flag=false, reuse_transition_mats_flag=false,
                               ltl_formula=nothing)

    refinement_result_dirs = []
    synth_dir = "$results_path/$system_tag/reach-specification"
    !isdir(synth_dir) && mkpath(synth_dir)

    res_mat, res_dir, pimdp = perform_synthesis_from_result_dirs(system_paths, synth_dir, "", experiment_params_array[1].specification_file, label_fcn; add_opt=true, horizon=horizon, ltl_formula=ltl_formula)

    save_legacy_mats(res_mat, synth_dir, horizon)

    for i=1:refinement_steps
        results_dir = "$synth_dir/refinement$i"

        num_dfa_states = Int(length(pimdp.states)/(length(keys(pimdp.imdp_label_dict))))
        states_to_update = unique(Int.(pimdp_col_to_imdp_state.(findall(x->x<minimum_threshold, res_mat[1:end-1,3]) ∩ findall(x->x>minimum_threshold, res_mat[1:end-1,4]), num_dfa_states)))

        new_imdp_dirs = Array{String,1}()
        # For each IMDP in the system
        for (j, component_dir) in enumerate(system_paths)
            imdp_results_dir = "$results_dir/mode$j"
            refine_imdp(states_to_update, experiment_params_array[j], working_dir=component_dir, results_dir=imdp_results_dir, reuse_regions_flag=reuse_regions_flag, reuse_transition_mats_flag=reuse_transition_mats_flag)
            push!(new_imdp_dirs, imdp_results_dir)
        end

        res_mat, res_dir, pimdp = perform_synthesis_from_result_dirs(new_imdp_dirs, results_dir, "", experiment_params_array[1].specification_file, label_fcn; add_opt=true, horizon=horizon, ltl_formula=ltl_formula)
        cp("$results_dir/switched-system/regions.bson", "$results_dir/regions.bson", force=true)
        system_paths = new_imdp_dirs
        push!(refinement_result_dirs, results_dir)
    end

    return refinement_result_dirs
end

function end_to_end_transition_bounds(params; reuse_gps_flag=false, reuse_regions_flag=false)

	@info "Initializing experiment and saving the parameters..."
	exp_dir = create_experiment_directory(params)
	params.experiment_directory = exp_dir
    logfile = initialize_log(exp_dir)
	save_experiment_params(params)
	total_runtime = 0.
    timing_info = Dict()

	regions_filename = "$exp_dir/regions.bson"
	gps_filename = generate_gp_filename(params)

	if reuse_regions_flag && (isfile(regions_filename) && isfile(gps_filename))
		gp_set, gp_info_dict = load_gp_info(gps_filename)
		timing_info["data_generation_time_s"] = -1.
		timing_info["gp_regression_time_s"] = -1.
		region_data = load_region_data("$exp_dir/regions.bson")
		timing_info["data_generation_time_s"] = -1.
		timing_info["gp_regression_time_s"] = -1.
		timing_info["region_bound_time_s"] = -1.
		@goto reuse_regions
	elseif reuse_gps_flag && isfile(gps_filename)
		gp_set, gp_info_dict = load_gp_info(gps_filename)
		timing_info["data_generation_time_s"] = -1.
		timing_info["gp_regression_time_s"] = -1.
		@goto reuse_gps
	end

    @info "Generating the training data..."

    gen_time = @elapsed begin
        x_train, y_train = generate_training_data(params)
    end
    total_runtime += gen_time
    timing_info["data_generation_time_s"] = gen_time
    @info "Data generation time: " gen_time
	save_training_data(x_train, y_train, params.experiment_directory)

    @info "Generating the GP regressions..."
    gp_time = @elapsed begin
        gp_set, gp_info_dict = generate_estimates(params, x_train, y_train)
    end
    total_runtime += gp_time
    timing_info["gp_regression_time_s"] = gp_time
    @info "GP generation time: " gp_time
    save_gp_info(params, gp_set, gp_info_dict)

@label reuse_gps
    @info "Generating the region info..."
    region_time = @elapsed begin
        region_data = generate_region_images(params, gp_info_dict)
    end
    total_runtime += region_time
    timing_info["region_bound_time_s"] = region_time
    @info "Region generation time: " region_time
    save_region_data(exp_dir, region_data)

@label reuse_regions
    @info "Generating the transition bounds..."
    bound_time = @elapsed begin
        result_mats = generate_transition_bounds(params, gp_info_dict, region_data)
    end
    total_runtime += bound_time
    timing_info["transition_bound_time_s"] = bound_time
    @info "Bound generation time: " bound_time
    save_transition_matrices(params, result_mats)

    @info "Total runtime: " total_runtime
    timing_info["total_runtime_s"] = total_runtime
    save_time_info(exp_dir, timing_info)

    flush(logfile)
    close(logfile)

    return timing_info, result_mats
end

function end_to_end_transition_bounds_local_gps(params; single_mode_verification=false, reuse_regions_flag=false)
    exp_dir = create_experiment_directory(params)
	params.experiment_directory = exp_dir
    logfile = initialize_log(params)
    @info "Generating the training data..."
    total_runtime = 0.
    timing_info = Dict()
    gen_time = @elapsed begin
        x_train, y_train = generate_training_data(params)
    end
    total_runtime += gen_time
    timing_info["data_generation_time_s"] = gen_time
    @info "Data generation time: " gen_time
    @info "Generating the regressions and region data..."
    region_time = @elapsed begin
        region_data = generate_region_images(params, x_train, y_train)
    end
    total_runtime += region_time
    timing_info["region_bound_time_s"] = region_time
    @info "Region generation time: " region_time
    save_region_data(params, region_data)
    @info "Generating the transition bounds..."
    bound_time = @elapsed begin
        result_mats = generate_transition_bounds(params, region_data)
    end
    total_runtime += bound_time
    timing_info["transition_bound_time_s"] = bound_time
    @info "Bound generation time: " bound_time
    save_transition_matrices(params, result_mats)

    safety_result_mat = nothing
    if single_mode_verification
        @info "Performing safety verification on single mode for 1 step..."
        verification_time = @elapsed begin
            imdp = create_simple_imdp(result_mats["minPr"], result_mats["maxPr"])
            horizon=1
            verification_result_mat = Globally(imdp, "safe", horizon, "$exp_dir/imdp.txt")
            save_legacy_mats(verification_result_mat, exp_dir, horizon)
        end
        timing_info["single_verification_time_s"] = verification_time
        @info "Verification time: " verification_time
    end

    @info "Total runtime: " total_runtime
    timing_info["total_runtime_s"] = total_runtime
    save_time_info(params, timing_info)

    flush(logfile)
    close(logfile)

    return timing_info, safety_result_mat, result_mats
end

function perform_synthesis_from_result_dirs_safety(res_dirs, exp_dir, system_tag, label_fcn;
	plot_graphs=true, rerun_flag=false, modes=nothing, add_opt=true,
	region_specific="", num_sims=10, sim_points=nothing, horizon=-1)

    @info "Creating the IMDP..."
	imdp = create_imdp_from_result_dirs(res_dirs, "$exp_dir")
    @info("$exp_dir")
	regions_file = @sprintf("%s/regions%s.bson", res_dirs[1], region_specific)
	create_imdp_labels(label_fcn, imdp, regions_file)
	plot_graphs ? create_dot_graph(imdp, "$exp_dir/imdp.gv") : nothing

	# res_mat = run_pimdp_synthesis(pimdp, pimdp_filename, add_opt=add_opt)
    imdp_filename = "$exp_dir/imdp.txt"
    res_mat = CautiousEngine.Globally(imdp, "safe", horizon, imdp_filename, synthesis_flag=true)

    if add_opt
        @info "Synthesizing a controller... (maximize optimistic)"
        res_mat_opt = run_imdp_synthesis(imdp_filename, horizon, mode2="optimistic", save_mats=false)
        for j in 1:length(res_mat[:,1])
            if res_mat[j, 3] == res_mat_opt[j,3] && res_mat_opt[j,4] > res_mat[j,4]
                res_mat[j, 2] = res_mat_opt[j,2]
                res_mat[j, 4] = res_mat_opt[j,4]
            end
        end
    end

    dst_dir = "$exp_dir/safety"
	isdir(dst_dir) ? nothing : mkpath(dst_dir)
	# Copy region data
	sys_dir = res_dirs[1]
	cp("$sys_dir/regions.bson", "$dst_dir/regions.bson", force=true)

	# TODO Plot gamma values somewhere else, i.e. not here
	# if add_opt
	#     plot_gamma_value(dst_dir, res_mat, res_mat_opt, num_dfa_states=length(dfa.states))
	# end
	save_legacy_mats(res_mat, dst_dir, horizon)

	# plot_synthesis_results(dst_dir, res_mat, imdp, dfa, pimdp)
    return res_mat, dst_dir, -1
end

function perform_synthesis_from_result_dirs(res_dirs, exp_dir, system_tag, spec_file, label_fcn;
	plot_graphs=false, rerun_flag=false, modes=nothing, add_opt=false,
	region_specific="", num_sims=10, sim_points=nothing, horizon=-1, ltl_formula=nothing)

	@info "Creating the IMDP..."
	imdp = create_imdp_from_result_dirs(res_dirs, "$exp_dir/$system_tag")
	regions_file = @sprintf("%s/regions%s.bson", res_dirs[1], region_specific)
	create_imdp_labels(label_fcn, imdp, regions_file)
	plot_graphs ? create_dot_graph(imdp, "$exp_dir/$system_tag/imdp.gv") : nothing

	@info "Creating the DFA..."
    if !isnothing(ltl_formula)
        @info "Constructing dfa from ltl formula: $ltl_formula"
        dfa = construct_DFA_from_LTL(ltl_formula)
        spec_tag = string(ltl_formula)
    else
        include("$spec_file")
        spec_tag = basename(spec_file)[1:end-3]
        dfa = DFA(dfa_states, dfa_props, dfa_transitions, dfa_accepting_state, dfa_sink_state, dfa_initial_state)
    end
	dst_dir = "$exp_dir/$system_tag/$spec_tag"
	isdir(dst_dir) ? nothing : mkpath(dst_dir)
	plot_graphs ? create_dot_graph(dfa, "$dst_dir/specification.gv") : nothing

	@info "Creating pimdp..."
	pimdp = construct_DFA_IMDP_product(dfa, imdp)
	plot_graphs ? create_dot_graph(pimdp, "$dst_dir/pimdp.gv") : nothing

	pimdp_filename = "$dst_dir/pimdp.txt"
	res_mat = run_pimdp_synthesis(pimdp, pimdp_filename, add_opt=add_opt, horizon=horizon)

	# Copy region data
	sys_dir = res_dirs[1]
	cp("$sys_dir/regions.bson", "$dst_dir/regions.bson", force=true)

	save_legacy_mats(res_mat, dst_dir, horizon)

	plot_synthesis_results(dst_dir, res_mat, imdp, dfa, pimdp)

	if !isnothing(modes)
		@info "Running a simulation using the resulting controller..."
		trajectories = []
		mt = MersenneTwister(111)

		if !isnothing(sim_points)
			num_sims=length(sim_points)
		end

		for i=1:num_sims
			if !isnothing(sim_points)
				x0 = sim_points[i]
			else
				x0 = 4*rand(mt,2,1)[:]-[2.; 2.]
			end
			extents, _ = deserialize_region_data(regions_file)
			simulate_system(x0, modes, 100, imdp, pimdp, dfa, extents, res_mat[:,2])
			push!(trajectories, copy(pimdp.trajectory))
		end
		plot_synthesis_results(dst_dir, res_mat, imdp, dfa, pimdp, trajectories=trajectories, filename="sim.png")
	end

	return res_mat, dst_dir, pimdp
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct EntityType <: AbstractEntityType
    label::Label
    is_root::Bool
end

function create(x::RemoteConcept{EntityType})
    result = execute(x.transaction, EntityTypeRequestBuilder.create_req(x.concept.label))
    return Entity(result.type_res.entity_type_create_res.entity)
end


text:
function atan2(y,x)
    if x>0
        return atan(y/x)
    elseif (x<0 && y>=0)
        return atan(y/x)+pi
    elseif (x<0 && y<0)
        return atan(y/x)-pi
    elseif (x==0 && y>0)
        return pi/2
    elseif (x==0 && y<0)
        return -pi/2
    else
        return NaN
    end
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-SimpleHel

max_forks_repo_name: umbertosaetti/J-SimpleHel


text:
export ASB07

struct ASB07 <: ContinuousPost
    options::TwoLayerOptions

    function ASB07(𝑂::Options)
        𝑂new = validate_and_wrap_options(𝑂, options_ASB07())
        return new(𝑂new)
    end
end

# convenience constructor from pairs of symbols
ASB07(𝑂::Pair{Symbol,<:Any}...) = ASB07(Options(Dict{Symbol,Any}(𝑂)))

# default options (they are added in the function validate_and_wrap_options)
ASB07() = ASB07(Options())

include("options.jl")
include("init.jl")
include("post.jl")
include("reach.jl")


text:
# # Laguerre-Gauss profile

struct LaguerreGaussLaserConstantCache{IC,W,K,T,Z,E,N}
    inv_c::IC
    ω::W
    k::K
    T₀::T
    z_R::Z
    E₀::E
    Nₚₘ::N
end

function LaguerreGaussLaserConstantCache(;c, λ, w₀, a₀, mₑ, q, p, m)
    ω = 2π * c / λ
    k = 2π / λ
    T = 2π / ω
    z_R = w₀^2 * k / 2
    E₀ = a₀ * mₑ * c * ω / abs(q)
    mₐ = abs(m)
    Nₚₘ = √(pochhammer(p+1, mₐ))

    return LaguerreGaussLaserConstantCache(inv(c), ω, k, T, z_R, E₀, Nₚₘ)
end

@auto_hash_equals mutable struct LaguerreGaussLaserCache{L,S,CE,EE}
    x::L
    y::L
    σ::S
    wz::L
    Ex::CE
    Ey::CE
    Ez::CE
    Eg::CE
    NEgexp::EE
    rwz::S
end

function LaguerreGaussLaserCache(λ, E)
    ThreadLocal(LaguerreGaussLaserCache(
        zero(λ),                    # x
        zero(λ),                    # y
        zero(λ/λ),                  # σ
        zero(λ),                    # wz
        zero(E*im),                 # Ex
        zero(E*im),                 # Ey
        zero(E*im),                 # Ez
        zero(E*im),                 # Eg
        zero(E*im),                 # NEgexp
        zero(λ/λ),                  # rwz
    ))
end

function Base.fill!(cache::LaguerreGaussLaserCache, x::AbstractVector)
    cache.x = x[1]
    cache.y = x[2]

    return nothing
end

@doc """
    struct LaguerreGaussLaser <: AbstractLaser

The `LaguerreGaussLaser` is defined by the `units` of the laser
(a positional argument that can be `:SI`, `:atomic` or `:SI_unitful` and `:atomic_unitful`)
and the following parameters
- `λ` is the laser wavelangth
- `a₀` is the normalized vector potential (defined as ``a_0=\\frac{eA}{m_e c^2}``)
- `ϕ₀` is the initial phase with the default value 0.0
- `w₀` is the beam radius at the Rayleigh range or [beam waist](https://en.wikipedia.org/wiki/Gaussian_beam#Beam_waist)
- `ξx` and `ξy` give the polarization and have the default value `1.0 + 0im` and `0.0 + 0im`
- `orientation` specifies how the laser is oriented with respect to the
default coordinate system (default `(:x, :z)`). See [`LaserGeometry`](@ref) for more details.
- `propagation_dir` is the propagation direction of the wave in the
intrinsic coordinate system (default `:z`). See [`LaserGeometry`](@ref) for more details.
- `profile` is the temporal profile of the pulse and the default one is constant (infinite pules duration)
- `p` is the radial index of the mode, ``p ∈ ℤ, p ≥ 0``, with the default value 1
- `m` is the azimuthal index of the mode, ``m ∈ ℤ``, with the default value 0

During the initialization of the `LaguerreGaussLaser` type, some useful derived values are
also computed
- `ω` is the angular frequency and is given by ``2π c / λ``
- `k` is the wavenumber and is given by ``2π / λ``
- `z_R` is the [Rayleigh range](https://en.wikipedia.org/wiki/Rayleigh_length) and is given by ``k w_0^2 / 2``
- `T₀` is the laser period and is given by ``2π / ω``
- `E₀` is the amplitude of the electric field and is given by ``a_0\\frac{m_q c \\omega}{q}``
- `Nₚₘ` is a normalization factor given by ``\\sqrt{(p+1)_{|m|}}``, with ``(x)_n`` the [Pochhammer symbol](https://mathworld.wolfram.com/PochhammerSymbol.html)
"""
LaguerreGaussLaser

struct LaguerreGaussLaser{_P,_M,_MA,
                          C0,Q,M,Eps,Mu,U,
                          IC,W,K,T,Z,E,F,
                          L,S,CE,EE,
                          D,R,
                          C,
                          P,
                          I} <: AbstractLaser
    constants::FundamentalConstants{C0,Q,M,Eps,Mu,U}
    derived::LaguerreGaussLaserConstantCache{IC,W,K,T,Z,E,F}
    cache::ThreadLocal{LaguerreGaussLaserCache{L,S,CE,EE}}
    geometry::LaserGeometry{D,R}
    polarization::LaserPolarization{C}
    profile::P
    # laser parameters
    λ::L
    a₀::F
    ϕ₀::F
    w₀::L
    p::I
    m::I
end

# This hack stores the m and p values in the type domain
function LaguerreGaussLaser(
    constants::FundamentalConstants{C0,Q,M,Eps,Mu,U},
    derived::LaguerreGaussLaserConstantCache{IC,W,K,T,Z,E,F},
    cache::ThreadLocal{LaguerreGaussLaserCache{L,S,CE,EE}},
    geometry::LaserGeometry{D,R},
    polarization::LaserPolarization{C},
    profile::P,
    λ::L,
    a₀::F,
    ϕ₀::F,
    w₀::L,
    p::I,
    m::I
) where {C0,Q,M,Eps,Mu,U,
         IC,W,K,T,Z,E,F,
         L,S,CE,EE,
         D,R,
         C,
         P,
         I}
    LaguerreGaussLaser{Val{p},Val{m},Val{abs(m)},
                       C0,Q,M,Eps,Mu,U,
                       IC,W,K,T,Z,E,F,
                       L,S,CE,EE,
                       D,R,
                       C,
                       P,
                       I}(
        constants,
        derived,
        cache,
        geometry,
        polarization,
        profile,
        λ,
        a₀,
        ϕ₀,
        w₀,
        p,
        m
    )
end

function LaguerreGaussLaser(units;
        λ,
        a₀,
        ϕ₀ = 0.0,
        w₀,
        p = 1,
        m = 0,
        ξx = 1.0+0im,
        ξy = 0,
        orientation = (:x, :z),
        profile = ConstantProfile()
    )

    ξx, ξy = promote(ξx, ξy)
    a₀, ϕ₀ = promote(a₀, ϕ₀)
    λ, w₀ = promote(λ, w₀)
    @assert hypot(ξx, ξy) ≈ 1 "Invalid ξx and ξy"
    constants = FundamentalConstants(units)
    @unpack mₑ, c, q = constants

    derived = LaguerreGaussLaserConstantCache(; c, λ, w₀, a₀, mₑ, q, p, m)
    E₀ = derived.E₀

    cache = LaguerreGaussLaserCache(λ, E₀)

    geometry = LaserGeometry(orientation)

    polarization = LaserPolarization(ξx, ξy)

    return LaguerreGaussLaser(
        constants,
        derived,
        cache,
        geometry,
        polarization,
        profile,
        λ,
        a₀,
        ϕ₀,
        w₀,
        p,
        m
    )
end

"""
    convert_laser(::Type{GaussLaser}, laser::LaguerreGaussLaser)

Convert a `LaguerreGaussLaser` to a `GaussLaser` with similar parameters.
"""
function convert_laser(::Type{GaussLaser}, laser::LaguerreGaussLaser)
    @unpack constants, geometry, profile, λ, a₀, w₀ = laser

    derived = GaussLaserConstantCache(;
        constants.c,
        λ,
        w₀,
        a₀,
        constants.mₑ,
        constants.q)
    # E₀ = derived.E₀

    # Don't create a cache since it's not used and it allocates memory
    cache = nothing # GaussLaserCache(λ, E₀)

    polarization = LaserPolarization(1, 0)

    GaussLaser(
        constants,
        derived,
        cache,
        geometry,
        polarization,
        profile,
        λ,
        a₀,
        zero(a₀),
        w₀
    )
end

get_p(::Type{<:LaguerreGaussLaser{Val{P},Val{M}}}) where {P,M} = P
get_m(::Type{<:LaguerreGaussLaser{Val{P},Val{M}}}) where {P,M} = M
get_mₐ(::Type{<:LaguerreGaussLaser{Val{P},Val{M},Val{MA}}}) where {P,M,MA} = MA

function required_coords(::LaguerreGaussLaser, r)
    CylindricalFromCartesian()(r)
end

function Ex(laser::LaguerreGaussLaser, coords)
    @unpack Nₚₘ, z_R = immutable_cache(laser)
    mₐ = get_mₐ(typeof(laser))
    ξx = polarization(laser, :ξx)
    @unpack ϕ₀, p, m, cache = laser
    @unpack r, θ, z = coords

    gauss_laser = convert_laser(GaussLaser, laser)
    Eg = Ex(gauss_laser, coords)
    # We avoid computing wz in the GaussLaser call so that we can skip cache creation
    wz = w(z, laser)
    σ = (r/wz)^2
    rwz = ustrip(NoUnits, r*√2/wz)
    @pack! cache[] = Eg, wz, σ, rwz

    ξx * Eg * Nₚₘ * rwz^mₐ * _₁F₁(-p, mₐ+1, 2σ) * exp(im*((2p+mₐ)*atan(z, z_R) - m*θ + ϕ₀))
end

function Ez(laser::T, coords) where T <: LaguerreGaussLaser
    @unpack Nₚₘ, k, z_R = immutable_cache(laser)
    @unpack wz, rwz, σ, Eg, Ex, Ey, x, y = mutable_cache(laser)
    @unpack ξx, ξy = polarization(laser)
    @unpack ϕ₀, p, m = laser
    @unpack r, θ, z = coords

    P = get_p(T)
    M = get_m(T)
    mₐ = get_mₐ(T)

    sgn = sign(M)
    𝟘 = zero(typeof(Ex))/oneunit(typeof(x))
    NEgexp = Nₚₘ*Eg*exp(im*((2p+mₐ)*atan(z, z_R)-m*θ+ϕ₀))
    update_cache!(laser, :NEgexp, NEgexp)
    # Compute (√2/wz)^mₐ * r^(mₐ-1) with the static value of |m|, so that it can be inferred in the unitful case.
    rwz_term = Base.literal_pow(^, (√2/wz), Val{mₐ}()) * Base.literal_pow(^, r, Val{mₐ-1}())

    -im/k * (
       (iszero(M) ? 𝟘 : mₐ * (ξx - im*sgn*ξy) * rwz_term * _₁F₁(-p, mₐ+1, 2σ) * NEgexp * exp(im*sgn*θ))
     - 2/(wz^2) * (1 + im*z/z_R) * (x*Ex + y*Ey)
     - (iszero(P) ? 𝟘 : 4p/((mₐ+1) * wz^2) * (x*ξx + y*ξy) * rwz^mₐ * _₁F₁(-p+1, mₐ+2, 2σ) * NEgexp)
    )
end

function Bz(laser::T, coords) where T <: LaguerreGaussLaser
    @unpack ω, z_R = immutable_cache(laser)
    @unpack wz, rwz, σ, Ex, Ey, NEgexp, x, y = mutable_cache(laser)
    @unpack ξx, ξy = polarization(laser)
    @unpack p, m = laser
    @unpack r, θ, z = coords

    P = get_p(T)
    M = get_m(T)
    mₐ = get_mₐ(T)

    sgn = sign(M)
    𝟘 = zero(typeof(Ex))/oneunit(typeof(x))

    -im/ω * (
       - (iszero(M) ? 𝟘 : mₐ * (ξy + im*sgn*ξx) * (√2/wz)^mₐ*r^(mₐ-1) * _₁F₁(-p, mₐ+1, 2σ) * NEgexp * exp(im*sgn*θ))
       + 2/(wz^2) * (1 + im*z/z_R) * (x*Ey - y*Ex)
       + (iszero(P) ? 𝟘 : (4p)/((mₐ+1) * wz^2) * (x*ξy - y*ξx) * rwz^mₐ * _₁F₁(-p+1, mₐ+2, 2σ) * NEgexp)
    )
end


text:
module OutputWriters

export
    FieldSlicer,
    JLD2OutputWriter, NetCDFOutputWriter,
    Checkpointer, restore_from_checkpoint,
    WindowedTimeAverage,
    TimeInterval, IterationInterval, WallTimeInterval, AveragedTimeInterval

using CUDA

using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Fields
using Oceananigans.Models
using Oceananigans.Utils: TimeInterval, IterationInterval, WallTimeInterval

using Oceananigans: AbstractOutputWriter
using Oceananigans.Fields: OffsetArray

import Oceananigans: write_output!

Base.open(ow::AbstractOutputWriter) = nothing
Base.close(ow::AbstractOutputWriter) = nothing

include("output_writer_utils.jl")
include("fetch_output.jl")
include("windowed_time_average.jl")
include("time_average_outputs.jl")
include("jld2_output_writer.jl")
include("netcdf_output_writer.jl")
include("checkpointer.jl")

end


text:
using Revise
using Plots
using DynamicalSystems
using DifferentialEquations
using Basins

μ=0.47
ds = Systems.rikitake(μ = μ, α = 1.0)
integ=integrator(ds)
xg=range(-6.,6.,length=200)
yg=range(-6.,6.,length=200)
pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9)
@time bsn = Basins.basins_map2D(xg, yg, pmap)
plot(xg,yg,bsn.basin',seriestype=:heatmap)