text:
################################################################################
#
#   TYPE DEFINITIONS OF OBJECT CLASSES IN JULIA
#
#   STRUCTURE OF THE FILE:
#
#   1) Definition of UNITCELL OBJECT
#       - Saving
#       - Loading
#
#   2) Definition of LATTICE OBJECT
#       - Saving
#       - Loading
#
#   3) CONVERSION LATTICE --> UNITCELL
#
#   4) OBTAINING GENERAL INFORMATION
#       - Testing a Unitcell object
#       - Print Information (printInfo(..))
#       - obtain the connection strengths of an object
#
#   5) OBTAINING CONNECTIVITY INFORMATION
#       - getConnectivityList(..)
#       - getConnectionList(..)
#
################################################################################

#-----------------------------------------------------------------------------------------------------------------------------
#
#   THE UNIT CELL OBJECT CLASS IN JULIA
#
#   FIELDS:
#   - lattice_vectors   lattice vectors of the underlying bravais lattice
#   - basis             basis sites of the unit cell
#   - connections       the connetions between sites of the unit cells
#   - filename          The filename that the unitcell should be stored into
#
#-----------------------------------------------------------------------------------------------------------------------------
"""
    mutable struct Unitcell

The type that contains information on a unitcell. Fields are

    lattice_vectors :: Array{Array{Float64, 1}, 1}
    basis           :: Array{Array{Float64, 1}, 1}
    connections     :: Array{Array{Any, 1}, 1}
    filename        :: String

Note that connections have the format

    [index_from::Int, index_to::Int, strength::Any, unitcell_offset::Tuple]

New objects can be created with

    Unitcell(filename::String)

to load from a given file with name `filename` or with

    Unitcell()

if the constructed object serves as a placeholder. In this case, the filename
is set to `UNITCELL_DUMMY_FILENAME`.

# Examples

```julia-repl
julia> Unitcell("myunitcellfile.jld")
```
"""
mutable struct Unitcell

    # FIELDS

    # basis vectors of the Bravais lattice
    lattice_vectors::Array{Array{Float64, 1}, 1}

    # basis sites within the unitcell
    basis::Array{Array{Float64, 1}, 1}

    # list of connections of the format
    # [index_from, index_to, strength, unitcell_offset]
    connections::Array{Array{Any, 1}, 1}

    # filename if saved
    filename::String

    # the general constructor
    function Unitcell(lattice_vectors::Array{Array{Float64, 1}, 1}, basis::Array{Array{Float64, 1}, 1}, connections::Array{Array{Any,1},1}, filename::String)
        return new(lattice_vectors, basis, connections, filename)
    end
    # the custom constructor when loading a unitcell from a file
    function Unitcell(filename::String)
        # define an empty unitcell
        uc = new(Array{Float64, 1}[],Array{Float64, 1}[],Array{Any, 1}[],filename)
        # load the uc
        loadUnitcell(uc)
        # return the uc
        return uc
    end
    # the dummy constructor (for lattices that do not need a UC)
    function Unitcell()
        # just initialize everything
        return new(Array{Float64, 1}[],Array{Float64, 1}[],Array{Any, 1}[],UNITCELL_DUMMY_FILENAME)
    end

end

#-----------------------------------------
#
#   METHODS FOR SAVING UNITCELLS
#
#-----------------------------------------

"""
    saveUnitcell(uc::Unitcell [, filename::String])

Save the `Unitcell` object `unitcell` to a file by using the module `JLD`.
A new filename can be given by passing an additional argument.

Note that the filename has to end with `.jld` or the ending will be added automatically.

Folders that lead to the specified file will be created on the way.

# Examples

```julia-repl
julia> saveUnitcell(unitcell)
"filename.jld"

julia> saveUnitcell(unitcell, "testfile.jld")
"testfile.jld"
```
"""
function saveUnitcell(uc::Unitcell, filename::String="NONE")
    # maybe set the new filename
    if filename != "NONE"
        uc.filename = filename
    end
    # make sure that the file ends with .jld
    if uc.filename[end-3:end] != ".jld"
        uc.filename = "$(uc.filename).jld"
    end
    # ensure default path is build
    if contains(uc.filename, "/")
		# get the containing folder
		folder = uc.filename[1:findlast(uc.filename, '/')]
		# build the path to that folder
		mkpath(folder)
	end
    # open the file and write the fields into it
    save(uc.filename,
        "basis",            uc.basis,
        "connections",      uc.connections,
        "lattice vectors",  uc.lattice_vectors,
        "filename",         uc.filename
    )
    # return the filename
    return uc.filename
end

#-----------------------------------------
#
#   METHODS FOR LOADING UNITCELLS
#
#-----------------------------------------

"""
    loadUnitcell(uc::Unitcell [, filename::String])
    loadUnitcell(filename::String)

Loads a `Unitcell` object `unitcell` from its file by using the module `JLD`.
A new filename can be given by passing an additional argument.
Moreover, if there is no `Unitcell` object which should be loaded but one wants to obtain a new
object, one can only pass the filename of the unitcell file.

# Examples

```julia-repl
julia> loadUnitcell(unitcell)
"filename.jld"

julia> loadUnitcell(unitcell, "testfile.jld")
"testfile.jld"

julia> loadUnitcell("testfile.jld")
LatticePhysics.Unitcell[...]
```
"""
function loadUnitcell(uc::Unitcell, filename::String="NONE")
    # maybe set the new filename
    if filename != "NONE"
        uc.filename = filename
    end
    # open the file and load all fields from it
    uc.basis            = load(uc.filename, "basis")
    uc.connections      = load(uc.filename, "connections")
    uc.lattice_vectors  = load(uc.filename, "lattice vectors")
    # return the filename
    return uc.filename
end
function loadUnitcell(filename::String)
    # return the filename
    return Unitcell(filename)
end

# make the type public accessible
export Unitcell
export loadUnitcell, saveUnitcell

#-----------------------------------------------------------------------------------------------------------------------------
#
#   THE LATTICE OBJECT CLASS IN JULIA
#
#   FIELDS:
#   - unitcell              The unit cell of the lattice (if present)
#   - unitcellRepetitions   Number of repetitions for each lattice direction (integer array with length of lattice vectors of UC)
#   - lattice_vectors       all periodic lattice vectors that span the complete lattice (empty list if finite lattice)
#   - positions             all sites of the lattice
#   - connections           the connetions between sites
#   - filename              The filename that the lattice should be stored into
#
#-----------------------------------------------------------------------------------------------------------------------------
"""
    mutable struct Lattice

The type that contains information on a lattice. Fields are

	unitcell            :: Unitcell
	unitcellRepetitions :: Array{Int64, 1}
    lattice_vectors     :: Array{Array{Float64, 1}, 1}
    positions           :: Array{Array{Float64, 1}, 1}
	positions_indices   :: Array{Int64, 1}
    connections         :: Array{Array{Any, 1}, 1}
    filename            :: String

Note that connections have the format

    [index_from::Int, index_to::Int, strength::Any, lattice_offset::Tuple]

where `lattice_offset` only describes the offset along the periodic directions

New objects can be created with

    Lattice(filename::String)

to load from a given file with name `filename` or with

    Lattice()

if the constructed object serves as a placeholder. In this case, the filename
is set to `LATTICE_DUMMY_FILENAME`.

For a complete list of all constructors please consult the documentation.

# Examples

```julia-repl
julia> Lattice("mylatticefile.jld")
```
"""
type Lattice

    # FIELDS

	# The unitcell of the lattice if applicable
    unitcell::Unitcell
	# how often the unitcell is repeated in the elementary directions
    unitcellRepetitions::Array{Int64, 1}

	# The lattice vectors that span the LATTICE in periodic dimensions
    lattice_vectors::Array{Array{Float64, 1}, 1}

	# Array with the positions of all sites within the lattice
    positions::Array{Array{Float64, 1}, 1}
	# array of ints that give the index of the site in the original unit cell
    positions_indices::Array{Int64, 1}

	# Array with connections between sites (same format as in Unitcell)
    connections::Array{Array{Any, 1}, 1}

	# Filename of the jld file
    filename::String

	# general Constructor if everything is given
    function Lattice(
				unitcell::Unitcell,
				unitcellRepetitions::Array{Int64, 1},
				lattice_vectors::Array{Array{Float64, 1}, 1},
				positions::Array{Array{Float64, 1}, 1},
				positions_indices::Array{Int64, 1},
				connections::Array{Array{Any, 1}, 1},
				filename::String
			)
        # just initialize everything and set the position indices to 1
        return new(unitcell, unitcellRepetitions, lattice_vectors, positions, positions_indices, connections, filename)
    end

	# Constructor if no position indices should be given
    function Lattice(
				unitcell::Unitcell,
				unitcellRepetitions::Array{Int64, 1},
				lattice_vectors::Array{Array{Float64, 1}, 1},
				positions::Array{Array{Float64, 1}, 1},
				connections::Array{Array{Any, 1}, 1},
				filename::String
			)
        # just initialize everything and set the position indices to 1
        return new(unitcell, unitcellRepetitions, lattice_vectors, positions,ones(size(positions,1)), connections, filename)
    end

    # Constructor if the lattice is not constructed from a unitcell
    function Lattice(
				lattice_vectors::Array{Array{Float64, 1}, 1},
				positions::Array{Array{Float64, 1}, 1},
				connections::Array{Array{Any, 1}, 1},
				filename::String
			)
        # just initialize everything
        return new(Unitcell(), [], lattice_vectors, positions,ones(size(positions,1)), connections, filename)
    end

    # Constructor when loading from a file
    function Lattice(filename::String)
        # define an empty unitcell
        lattice = new(Unitcell(),Array{Int64, 1}[], Array{Float64, 1}[],Array{Float64, 1}[],Array{Int64, 1}[],Array{Any, 1}[],filename)
        # load the uc
        loadLattice(lattice)
        # return the uc
        return lattice
    end

    # Dummy constructor
    function Lattice()
        # define an empty unitcell
        lattice = new(Unitcell(),Array{Int64, 1}[], Array{Float64, 1}[],Array{Float64, 1}[],Array{Int64, 1}[],Array{Any, 1}[],LATTICE_DUMMY_FILENAME)
        # return the uc
        return lattice
    end

end

#-----------------------------------------
#
#   METHODS FOR SAVING LATTICES
#
#-----------------------------------------

"""
    saveLattice(lattice::Lattice [, filename::String])

Save the `Lattice` object `lattice` to a file by using the module `JLD`.
A new filename can be given by passing an additional argument.

Note that the filename has to end with `.jld` or the ending will be added automatically.

Folders that lead to the specified file will be created on the way.

# Examples

```julia-repl
julia> saveLattice(lattice)
"filename.jld"

julia> saveLattice(lattice, "testfile.jld")
"testfile.jld"
```
"""
function saveLattice(lattice::Lattice, filename::String="NONE")
    # maybe set the new filename
    if filename != "NONE"
        lattice.filename = filename
    end
    # make sure that the file ends with .jld
    if lattice.filename[end-3:end] != ".jld"
        lattice.filename = "$(lattice.filename).jld"
    end
    # ensure default path is build
    if contains(lattice.filename, "/")
		# get the containing folder
		folder = lattice.filename[1:findlast(lattice.filename, '/')]
		# build the path to that folder
		mkpath(folder)
	end
    # open the file and write the fields into it
    save(lattice.filename,
        "unitcell",             lattice.unitcell,
        "unitcell repetitions", lattice.unitcellRepetitions,
        "positions",            lattice.positions,
        "positions indices",    lattice.positions_indices,
        "connections",          lattice.connections,
        "lattice vectors",      lattice.lattice_vectors,
        "filename",             lattice.filename
    )
    # return the filename
    return lattice.filename
end

#-----------------------------------------
#
#   METHODS FOR LOADING LATTICES
#
#-----------------------------------------

"""
    loadLattice(lattice::Lattice [, filename::String])
    loadLattice(filename::String)

Loads a `Lattice` object `lattice` from its file by using the module `JLD`.
A new filename can be given by passing an additional argument.
Moreover, if there is no `Lattice` object which should be loaded but one wants to obtain a new
object, one can only pass the filename of the unitcell file. (Note, that this can also be done
by using the constructor `Lattice(filename)`)

# Examples

```julia-repl
julia> loadLattice(lattice)
"filename.jld"

julia> loadLattice(lattice, "testfile.jld")
"testfile.jld"

julia> loadLattice("testfile.jld")
LatticePhysics.Lattice[...]
```
"""
function loadLattice(lattice::Lattice, filename::String="NONE")
    # maybe set the new filename
    if filename != "NONE"
        lattice.filename = filename
    end
    # open the file and load all fields from it
    lattice.unitcell            = load(lattice.filename, "unitcell")
    lattice.unitcellRepetitions = load(lattice.filename, "unitcell repetitions")
    lattice.positions           = load(lattice.filename, "positions")
    lattice.positions_indices   = load(lattice.filename, "positions indices")
    lattice.connections         = load(lattice.filename, "connections")
    lattice.lattice_vectors     = load(lattice.filename, "lattice vectors")
    # return the filename
    return lattice.filename
end
function loadLattice(filename::String)
    # just return the object that the constructor would return
	return Lattice(filename)
end

# EXPORT the relevant types and methods
export Lattice
export saveLattice, loadLattice

#-----------------------------------------------------------------------------
#
#   METHODS FOR CONVERTING LATTICE --> UNITCELL
#
#-----------------------------------------------------------------------------

"""
	toUnitcell(lattice::Lattice [, filename::String])

Converts a `Lattice` object `lattice` to a `Unitcell` object by keeping
connections and lattice vectors and setting all sites as basis sites.
Optionally, a further argument can specify a new filename for the newly
created unitcell.

# Examples

```julia-repl
julia> toUnitcell(lattice)
LatticePhysics.Unitcell[...]

julia> toUnitcell(lattice, "unitcell_test.jld")
LatticePhysics.Unitcell[...]
```
"""
function toUnitcell(lattice::Lattice, filename::String="AUTOMATIC")
    # maybe set the new filename
    if filename == "AUTOMATIC"
        filename = replace(lattice.filename, ".jld", "_to_unitcell.jld")
    end
	# return a newly built unitcell object
    return Unitcell(
        lattice.lattice_vectors,
        lattice.positions,
        lattice.connections,
		filename
            )
end

# Export the conversion method
export toUnitcell

#-----------------------------------------------------------------------------
#
#   METHODS FOR PRINTING OR OBTAINING GENERAL INFORMATION
#
#-----------------------------------------------------------------------------

"""
    testUnitcell(unitcell::Unitcell, N_lattice_vectors::Int64, N_space_dimension::Int64)

Performs a series of tests on a given `Unitcell` object so see if it is implemented correctly.
Checks include the following:
1) test if the number of `lattice_vectors` is as expected
2) test if the number of spatial dimensions is correct for every basis site
3) test if the filename does not contain `"FOLDER_UNITCELLS"` (as in this case a \$ would be missing somewhere)
4) test if the connectivity is okay (if every connection has a counterpart)

The return value is a `Bool` indicating if the test is passed or not. If the test fails, the
failure reason is printed as well.

# Examples

```julia-repl
julia> testUnitcell(unitcell, 2, 2)
...
```
"""
function testUnitcell(unitcell::Unitcell, N_lattice_vectors::Int64, N_space_dimension::Int64)
    # check if the unitcell has enough lattice vectors
    if size(unitcell.lattice_vectors,1) != N_lattice_vectors
        println("Uncorrect number of lattice vectors ($(size(unitcell.lattice_vectors,1)) instead of $(N_lattice_vectors))")
        return false
    end
    # check if the space dimension for all sites is okay
    for p in unitcell.basis
        if size(p,1) != N_space_dimension
            println("Site ($(p)) has uncorrect space dimension ($(size(p,1)) instead of $(N_space_dimension))")
            return false
        end
    end
    # check if the filename does not contain FOLDER_UNITCELLS (this would mean someone forgot a "$")
    if contains(unitcell.filename, "FOLDER_UNITCELLS")
        println("filename ($(unitcell.filename)) still contains \"FOLDER_UNITCELLS\"")
        return false
    end
    # check if connectivity okay
    for c1 in unitcell.connections
        # so far no counterpart found
        counterpart = false
        # search for a counterpart in all connections
        for c2 in unitcell.connections
            # skip identical
            if c1==c2
                continue
            end
            # skip if indices of sites are not correct
            if c1[1] != c2[2] || c1[2] != c2[1]
                continue
            end
            # skip if connection strength not equal
            if c1[3] != c2[3]
                continue
            end
            # skip if the unitcell wraps dont sum up
            if sum([abs(element) for element in collect(c1[4]) .+ collect(c2[4])]) > 1e-9
                continue
            end
            # if everything passed up to here, a counterpart has been found
            counterpart = true
            break
        end
        if counterpart == false
            println("connection $(c1) does not have counterpart, connectivity broken")
            return false
        end
    end
    # now, return true as all tests are passed
    return true
end
export testUnitcell

"""
    testLattice(lattice::Lattice, N_lattice_vectors::Int64, N_space_dimension::Int64)

Performs a series of tests on a given `Lattice` object so see if it is implemented correctly.
Checks include the following:
1) test if the number of `lattice_vectors` is as expected
2) test if the number of spatial dimensions is correct for every site in the lattice
3) test if the filename does not contain `"FOLDER_LATTICES"` (as in this case a \$ would be missing somewhere)
4) test if the connectivity is okay (if every connection has a counterpart)

The return value is a `Bool` indicating if the test is passed or not. If the test fails, the
failure reason is printed as well.

# Examples

```julia-repl
julia> testLattice(lattice, 2, 2)
...
```
"""
function testLattice(lattice::Lattice, N_lattice_vectors::Int64, N_space_dimension::Int64)
    # check if the unitcell has enough lattice vectors
    if size(lattice.lattice_vectors,1) != N_lattice_vectors
        println("Uncorrect number of lattice vectors ($(size(lattice.lattice_vectors,1)) instead of $(N_lattice_vectors))")
        return false
    end
    # check if the space dimension for all sites is okay
    for p in lattice.positions
        if size(p,1) != N_space_dimension
            println("Site ($(p)) has uncorrect space dimension ($(size(p,1)) instead of $(N_space_dimension))")
            return false
        end
    end
    # check if the filename does not contain FOLDER_UNITCELLS (this would mean someone forgot a "$")
    if contains(lattice.filename, "FOLDER_LATTICES")
        println("filename ($(lattice.filename)) still contains \"FOLDER_LATTICES\"")
        return false
    end
    # check if connectivity okay
    for c1 in lattice.connections
        # so far no counterpart found
        counterpart = false
        # search for a counterpart in all connections
        for c2 in lattice.connections
            # skip identical
            if c1==c2
                continue
            end
            # skip if indices of sites are not correct
            if c1[1] != c2[2] || c1[2] != c2[1]
                continue
            end
            # skip if connection strength not equal
            if c1[3] != c2[3]
                continue
            end
            # skip if the lattice wraps dont sum up
            if sum([abs(element) for element in collect(c1[4]) .+ collect(c2[4])]) > 1e-9
                continue
            end
            # if everything passed up to here, a counterpart has been found
            counterpart = true
            break
        end
        if counterpart == false
            println("connection $(c1) does not have counterpart, connectivity broken")
            return false
        end
    end
    # now, return true as all tests are passed
    return true
end
export testLattice

"""
    printInfo(unitcell::Unitcell [; detailed::Bool=false])
    printInfo(lattice::Lattice [; detailed::Bool=false])

Prints (detailed) information about a `Unitcell` object or a `Lattice` object.
Information Includes Bravais lattice vectors, number of sites and number of bonds.
When detailed output is enabled, furthermore all sites and bonds are printed.
Last but now least, statistics is printed about how many bonds per site are found.

# Examples

```julia-repl
julia> printInfo(unitcell)
...

julia> printInfo(unitcell, detailed=true)
...

```
"""
function printInfo(unitcell::Unitcell; detailed::Bool=false)

    # Header of the information
    println("Information on the unitcell stored in file \"$(unitcell.filename)\":")

    # General information on Bravais lattice vectors
    println(" - periodicity given by $(size(unitcell.lattice_vectors,1)) lattice vectors:")
    for l in unitcell.lattice_vectors
        println("     - $(l)")
    end

    # Information on basis sites
    if detailed
        println(" - $(size(unitcell.basis,1)) sites in unitcell of dimension $(length(unitcell.basis[1])):")
        for site in unitcell.basis
            println("     - $(site)")
        end
    else
        println(" - $(size(unitcell.basis,1)) sites in unitcell of dimension $(length(unitcell.basis[1]))")
    end

    # Information on bonds / connections
    if detailed
        println(" - $(size(unitcell.connections,1)) connections in the unitcell:")
        for c in unitcell.connections
            if typeof(c[3]) == String
                println("     - from $(c[1]) to $(c[2]) (with warp $(c[4])): \"$(c[3])\"")
            else
                println("     - from $(c[1]) to $(c[2]) (with warp $(c[4])): $(c[3])")
            end
        end
    else
        println(" - $(size(unitcell.connections,1)) connections in the unitcell")
    end

    # Information on bond / connection statistics
    println(" - $(size(unitcell.connections,1)/size(unitcell.basis,1)) connections per site")
    # check statistics of connections
    nc = zeros(Int64, size(unitcell.basis,1))
    for i in 1:length(nc)
        for c in unitcell.connections
            if Int(c[1]) == i
                nc[i] += 1
            end
        end
    end
    tc = 0
    c = 0
    print(" - statistics of connections per site:")
    while tc < size(unitcell.basis, 1)
        cc = 0
        for ncc in nc
            if ncc == c
                cc += 1
            end
        end
        if cc != 0
            print(" $(c)($(cc))")
        end
        c = c+1
        tc += cc
    end
    println("")
    broken = false
    for c1 in unitcell.connections
        counterpart = false
        for c2 in unitcell.connections
            if c1==c2
                continue
            end
            if c1[1] != c2[2] || c1[2] != c2[1]
                continue # site indices not correct
            end
            if c1[3] != c2[3]
                continue # connection strength not equal
            end
            if sum([abs(element) for element in collect(c1[4]) .+ collect(c2[4])]) > 1e-9
                continue # wrap not equal
            end
            counterpart = true
            break
        end
        if counterpart == false
            broken = true
        end
    end
    if broken
        println(" - connectivity of unitcell is broken (connections not vice versa)")
    else
        println(" - connectivity of unitcell is okay (including periodic BC)")
    end
end
function printInfo(lattice::Lattice; detailed=false)
	# Print the header informations
    println("Information on the lattice stored in file \"$(lattice.filename)\":")

	# Information on lattice vectors, print all of them
    println(" - periodicity given by $(size(lattice.lattice_vectors,1)) lattice vectors:")
    for l in lattice.lattice_vectors
        println("     - $(l)")
    end

	# Information on sites, depending on detailed print information
    if detailed
        println(" - $(size(lattice.positions,1)) sites in lattice of dimension $(length(lattice.positions[1])):")
        for site in lattice.positions
            println("     - $(site)")
        end
    else
        println(" - $(size(lattice.positions,1)) sites in lattice of dimension $(length(lattice.positions[1]))")
    end

	# Information on connectioins, depending on detailed print information
    if detailed
        println(" - $(size(lattice.connections,1)) connections in the lattice:")
        for c in lattice.connections
            if typeof(c[3]) == String
                println("     - from $(c[1]) to $(c[2]) (with warp $(c[4])): \"$(c[3])\"")
            else
                println("     - from $(c[1]) to $(c[2]) (with warp $(c[4])): $(c[3])")
            end
        end
    else
        println(" - $(size(lattice.connections,1)) connections in the lattice")
    end
	# Number of connections per site
    println(" - $(size(lattice.connections,1)/size(lattice.positions,1)) connections per site")
    # check statistics of connections
    nc = zeros(Int64, size(lattice.positions,1))
    for i in 1:length(nc)
        for c in lattice.connections
            if Int(c[1]) == i
                nc[i] += 1
            end
        end
    end
    tc = 0
    c = 0
	# print the statistics of the connections
    print(" - statistics of connections per site:")
    while tc < size(lattice.positions, 1)
        cc = 0
        for ncc in nc
            if ncc == c
                cc += 1
            end
        end
        if cc != 0
            print(" $(c)($(cc))")
        end
        c = c+1
        tc += cc
    end
    println("")
	# Test the connectivity of the lattice
    broken = false
    for c1 in lattice.connections
        counterpart = false
        for c2 in lattice.connections
            if c1==c2
                continue
            end
            if c1[1] != c2[2] || c1[2] != c2[1]
                continue # site indices not correct
            end
            if c1[3] != c2[3]
                continue # connection strength not equal
            end
            if sum([abs(element) for element in collect(c1[4]) .+ collect(c2[4])]) > 1e-9
                continue # wrap not equal
            end
            counterpart = true
            break
        end
        if counterpart == false
            broken = true
        end
    end
	# Print the connectivity of the lattice
    if broken
        println(" - connectivity of lattice is broken (connections not vice versa)")
    else
        println(" - connectivity of lattice is okay (including periodic BC)")
    end
end
export printInfo

"""
    getConnectionStrengthList(unitcell::Unitcell)
	getConnectionStrengthList(lattice::Lattice)

Obtains which connection strength values are taken in the definition of a given
`Lattice` object `lattice` or a `Unitcell` object `unitcell`.
The functions return a list which contains all connection strengths in the order
of appearence.

# Examples

```julia-repl
julia> connection_strength_list = getConnectionStrengthList(lattice)
["tx", "ty", "tz"]
```
"""
function getConnectionStrengthList(unitcell::Unitcell)
    # list of connetion strengths
    cs_list = []
    # iterate over all connections
    for c in unitcell.connections
        if !(c[3] in cs_list)
            push!(cs_list, c[3])
        end
    end
    # return the list
    return cs_list
end
function getConnectionStrengthList(lattice::Lattice)
    # list of connetion strengths
    cs_list = []
    # iterate over all connections
    for c in lattice.connections
        if !(c[3] in cs_list)
            push!(cs_list, c[3])
        end
    end
    # return the list
    return cs_list
end

# export the FUNCTIONS
export getConnectionStrengthList

#-----------------------------------------------------------------------------
#
#   METHODS FOR CONNECTIVITY INFORMATION
#
#-----------------------------------------------------------------------------

"""
	getConnectivityList(lattice::Lattice)
    getConnectivityList(unitcell::Unitcell)

Obtains connectivity information of a `Lattice` object `lattice` or a `Unitcell` object `unitcell`.
The functions return a list which contains a list of Tuples for every site in the
passed objects. Therefore, `connectivity_list[i][j]` can be used to access the connection `j`
from site `i` which will give a Tuple `(to_index, strength)`.

# Examples

```julia-repl
julia> connectivity_list = getConnectivityList(lattice)
Array{...}[...]

julia> connectivity_list[2][1]   # connection 1 outgoing from site 2
(12, "tx")
```
"""
function getConnectivityList(lattice::Lattice)
    # list of lists
    lists = Array[]
    for p in lattice.positions
        push!(lists, [])
    end
    for c in lattice.connections
        push!(lists[Int(c[1])], (c[2], c[3]) )
    end
    # return the lists
    return lists
end
function getConnectivityList(unitcell::Unitcell)
    # list of lists
    lists = Array[]
    for p in unitcell.basis
        push!(lists, [])
    end
    for c in unitcell.connections
        push!(lists[Int(c[1])], (c[2], c[3]) )
    end
    # return the lists
    return lists
end

# export the FUNCTIONS
export getConnectivityList

"""
	getConnectionList(lattice::Lattice)
    getConnectionList(unitcell::Unitcell)

Obtains connectivity information of a `Lattice` object `lattice` or a `Unitcell` object `unitcell`.
The functions return a list which contains a list of connection Arrays for every site in the
passed objects. Therefore, `connectivity_list[i][j]` can be used to access the connection `j`
from site `i` which will give an Array `[from_index, to_index, strength, wrap_tuple]`
like it is saved in the list of all connections of the respective object.

# Examples

```julia-repl
julia> connection_list = getConnectionList(lattice)
Array{...}[...]

julia> connection_list[2][1]   # connection 1 outgoing from site 2
[1, 12, "tx", (0,0)]
```
"""
function getConnectionList(lattice::Lattice)
    # list of lists
    lists = Array[]
    for p in lattice.positions
        push!(lists, Array[])
    end
    for c in lattice.connections
        push!(lists[Int(c[1])], c )
    end
    # return the lists
    return lists
end
function getConnectionList(unitcell::Unitcell)
    # list of lists
    lists = Array[]
    for p in unitcell.basis
        push!(lists, Array[])
    end
    for c in unitcell.connections
        push!(lists[Int(c[1])], c )
    end
    # return the lists
    return lists
end

# export the FUNCTIONS
export getConnectionList


text:
#-------------------------------------------------------------------------------
# STATE

# See docs/src/devnotes.md.
"""
A mutable state type of transcoding streams.

See Developer's notes for details.
"""
mutable struct State
    mode::Symbol        # current stream mode; in {:idle, :read, :write, :stop, :close, :panic}
    code::Symbol        # return code of the last method call; in {:ok, :end, :error}
    stop_on_end::Bool   # flag to go :stop on :end
    error::Error        # exception thrown while data processing
    buffer1::Buffer     # data buffers
    buffer2::Buffer     # "

    State(buffer1::Buffer, buffer2::Buffer) = new(:idle, :ok, false, Error(), buffer1, buffer2)
end

State(size::Integer) = State(Buffer(size), Buffer(size))

#-------------------------------------------------------------------------------
# ERROR

"""
Container of transcoding error.
An object of this type is used to notify the caller of an exception that
happened inside a transcoding method.  The `error` field is undefined at first
but will be filled when data processing failed. The error should be set by
calling the `setindex!` method (e.g. `error[] = ErrorException("error!")`).
"""
mutable struct Error
    error::Exception
    Error() = new()
end

# Test if an exception is set.
haserror(error::Error) = isdefined(error, :error)
Base.setindex!(error::Error, ex::Exception) = ((@assert !haserror(error) "an error is already set"); error.error = ex; error)
Base.getindex(error::Error) = (@assert haserror(error) "no error is set"); error.error

ext: jl


text:
        # intersect
        # -------------------------------------------
        # basic intersection
        l₁ = LineSegment2D(Point2D{T}(0,  1), Point2D{T}(2, -1))
        l₂ = LineSegment2D(Point2D{T}(0, -1), Point2D{T}(2,  1))
        hit, p₁ = intersect(l₁, l₂)
        @test hit
        @test p₁ ≈ Point2D{T}(1, 0)
        @test typeof(p₁) == Point2D{T}

        # vertex intersection
        l₂ = LineSegment2D(Point2D{T}(0, -1), Point2D{T}(2, -1))
        hit, p₁ = l₁ ∩ l₂
        @test hit
        @test p₁ ≈ Point2D{T}(2, -1)

        # vertical
        l₁ = LineSegment2D(Point2D{T}(0,  1), Point2D{T}(2,   1))
        l₂ = LineSegment2D(Point2D{T}(1, 10), Point2D{T}(1, -10))
        hit, p₁ = intersect(l₁, l₂)
        @test hit
        @test p₁ ≈ Point2D{T}(1, 1)

        # nearly vertical
        l₁ = LineSegment2D(Point2D{T}(-1, -100000), Point2D{T}(1,  100000))
        l₂ = LineSegment2D(Point2D{T}(-1,   10000), Point2D{T}(1,  -10000))
        hit, p₁ = l₁ ∩ l₂
        @test hit
        @test p₁ ≈ Point2D{T}(0, 0)

        # parallel
        l₁ = LineSegment2D(Point2D{T}(0, 1), Point2D{T}(1, 1))
        l₂ = LineSegment2D(Point2D{T}(0, 0), Point2D{T}(1, 0))
        hit, p₁ = intersect(l₁, l₂)
        @test !hit

        # collinear
        l₁ = LineSegment2D(Point2D{T}(0, 0), Point2D{T}(2, 0))
        l₂ = LineSegment2D(Point2D{T}(0, 0), Point2D{T}(1, 0))
        hit, p₁ = intersect(l₁, l₂)
        @test !hit

        # line intersects, not segment (invalid s)
        l₁ = LineSegment2D(Point2D{T}(0, 0), Point2D{T}(2, 0    ))
        l₂ = LineSegment2D(Point2D{T}(1, 2), Point2D{T}(1, 1//10))
        hit, p₁ = l₁ ∩ l₂
        @test !hit

        # line intersects, not segment (invalid r)
        l₂ = LineSegment2D(Point2D{F}(0, 0), Point2D{F}(2, 0    ))
        l₁ = LineSegment2D(Point2D{F}(1, 2), Point2D{F}(1, 1//10))
        hit, p₁ = intersect(l₁, l₂)
        @test !hit

        # intersect
        l₁ = LineSegment3D(Point3D{T}(0,  1, 0), Point3D{T}(2, -1, 0))
        l₂ = LineSegment3D(Point3D{T}(0, -1, 0), Point3D{T}(2,  1, 0))
        hit, p₁ = intersect(l₁, l₂)
        @test hit
        @test p₁ ≈ Point3D{T}(1, 0, 0)
        @test typeof(p₁) == Point3D{T}

        # intersect
        𝘅₁ = Point2D{T}(0, 0)
        𝘅₂ = Point2D{T}(2, 0)
        𝘅₃ = Point2D{T}(1, 1)
        𝘅₄ = Point2D{T}(1, 0)
        𝘅₅ = Point2D{T}(1, 2)

        # 1 intersection, straight
        q = QuadraticSegment2D(𝘅₁, 𝘅₂, Point2D{T}(1//2, 0))
        l = LineSegment2D(Point2D{T}(1,-1), Point2D{T}(1,1))
        npoints, (point1, point2) = intersect(l, q)
        @test npoints == 1
        @test point1 ≈ Point2D{T}(1, 0)

        # 1 intersection
        q = QuadraticSegment2D(𝘅₁, 𝘅₂, 𝘅₃)
        l = LineSegment2D(𝘅₄, 𝘅₅)
        npoints, (point1, point2) = intersect(l, q)
        @test npoints == 1
        @test point1 ≈ Point2D{T}(1, 1)

        # 2 intersections
        𝘅₄ = Point2D{T}(0, 3//4)
        𝘅₅ = Point2D{T}(2, 3//4)
        l = LineSegment2D(𝘅₄, 𝘅₅)
        npoints, (point1, point2) = l ∩ q
        @test npoints == 2
        @test point1 ≈ Point2D{T}(1//2, 3//4)
        @test point2 ≈ Point2D{T}(3//2, 3//4)

        # 0 intersections
        𝘅₄ = Point2D{T}(0, 3)
        𝘅₅ = Point2D{T}(2, 3)
        l = LineSegment2D(𝘅₄, 𝘅₅)
        npoints, (point1, point2) = intersect(l, q)
        @test npoints == 0

        # intersect
        hit, point = LineSegment(Point2D{F}(0, 1), Point2D{F}(1, 0)) ∩ plane
        @test hit
        @test point ≈ Point2D{F}(1//2, 1//2)

        # Line is in the plane
        hit, point = LineSegment(Point2D{F}(0, 0), Point2D{F}(1, 1)) ∩ plane
        @test !hit

        # Segment stops before plane
        hit, point = LineSegment(Point2D{F}(0, 2), Point2D{F}(1, 3//2)) ∩ plane
        @test !hit

        # Plane is before segment
        hit, point = LineSegment(Point2D{F}(1, 0), Point2D{F}(2, -1)) ∩ plane
        @test !hit

        # intersect
        hit, point = LineSegment(Point3D{F}(1, 2, 0), Point3D{F}(1, 2, 5)) ∩ plane
        @test hit
        @test point ≈ Point3D{F}(1,2,2)

        # Line is in the plane
        hit, point = LineSegment(Point3D{F}(0, 0, 2), Point3D{F}(1, 0, 2)) ∩ plane
        @test !hit

        # Segment stops before plane
        hit, point = LineSegment(Point3D{F}(1, 2, 0), Point3D{F}(1, 2, 1)) ∩ plane
        @test !hit

        # Plane is before segment
        hit, point = LineSegment(Point3D{F}(1, 2, 1), Point3D{F}(1, 2, 0)) ∩ plane
        @test !hit

        #isleft
        l = LineSegment(Point3D{F}(0, 0, 2), Point3D{F}(1, 0, 2))
        @test isleft(Point3D{F}(1,1,2), l, plane)
        @test isleft(Point3D{F}(1,0,2), l, plane)
        @test !isleft(Point3D{F}(1,-1,2), l, plane)

        # intersect
        # Horizontal
        hit, points = LineSegment2D(Point2D{T}(-1, 1), Point2D{T}(4, 1)) ∩ aab
        @test hit
        @test points[1] ≈ Point2D{T}(1, 1)
        @test points[2] ≈ Point2D{T}(3, 1)

        # Horizontal miss
        hit, points = LineSegment2D(Point2D{T}(-1, 5), Point2D{T}(4, 5)) ∩ aab
        @test !hit

        # Vertical
        hit, points = LineSegment2D(Point2D{T}(2, -1), Point2D{T}(2, 5)) ∩ aab
        @test hit
        @test points[1] ≈ Point2D{T}(2, 0)
        @test points[2] ≈ Point2D{T}(2, 2)

        # Vertical miss
        hit, points = LineSegment2D(Point2D{T}(5, -1), Point2D{T}(5, 5)) ∩ aab
        @test !hit

        # Angled
        hit, points = LineSegment2D(Point2D{T}(0, 0), Point2D{T}(5, 2)) ∩ aab
        @test hit
        @test points[1] ≈ Point2D{T}(1, 0.4)
        @test points[2] ≈ Point2D{T}(3, 1.2)

        # Angled miss
        hit, points = LineSegment2D(Point2D{T}(0, 5), Point2D{T}(5, 7)) ∩ aab
        @test !hit

        # intersect
        # Horizontal
        hit, points = LineSegment3D(Point3D{T}(-1, 1, 0), Point3D{T}(4, 1, 0)) ∩ aab
        @test hit
        @test points[1] ≈ Point3D{T}(1, 1, 0)
        @test points[2] ≈ Point3D{T}(3, 1, 0)

        # Horizontal miss
        hit, points = LineSegment3D(Point3D{T}(-1, 5, 0), Point3D{T}(4, 5, 0)) ∩ aab
        @test !hit

        # Vertical
        hit, points = LineSegment3D(Point3D{T}(2, -1, 0), Point3D{T}(2, 5, 0)) ∩ aab
        @test hit
        @test points[1] ≈ Point3D{T}(2, 0, 0)
        @test points[2] ≈ Point3D{T}(2, 2, 0)

        # Vertical miss
        hit, points = LineSegment3D(Point3D{T}(5, -1, 0), Point3D{T}(5, 5, 0)) ∩ aab
        @test !hit

        # Angled
        hit, points = LineSegment3D(Point3D{T}(0, 0, -1), Point3D{T}(5, 2, 1)) ∩ aab
        @test hit
        @test points[1] ≈ Point3D{T}(1, 0.4, -0.6)
        @test points[2] ≈ Point3D{T}(3, 1.2,  0.2)

        # Angled miss
        hit, points = LineSegment3D(Point3D{T}(0, 5, -10), Point3D{T}(5, 7, 20)) ∩ aab
        @test !hit

                # intersect
        # 3 intersections
        l = LineSegment2D(p₁, Point2D{T}(1, 1))
        hit, points = intersect(l, tri)
        @test hit
        @test points[1] ≈ p₁
        @test points[2] ≈ Point2D{T}(1//2, 1//2)

        # 2 intersections
        l = LineSegment2D(Point2D{T}(0, 1//2), Point2D{T}(1//2, 0))
        hit, points = intersect(l, tri)
        @test hit
        @test points[1] ≈ Point2D{T}(1//2, 0)
        @test points[2] ≈ Point2D{T}(0, 1//2)

        # 0 intersections
        l = LineSegment2D(Point2D{T}(-1, -1), Point2D{T}(2, -1))
        hit, points = intersect(l, tri)
        @test !hit

                # intersect
        # line is not coplanar with triangle
        p₄ = Point3D{T}(-1, 1//10, 1//10)
        p₅ = Point3D{T}( 1, 1//10, 1//10)
        l = LineSegment(p₄, p₅)
        hit, point = intersect(l, tri)
        @test hit
        @test point ≈ Point3D{T}(0, 1//10,  1//10)

        # line is coplanar with triangle
        p₄ = Point3D{T}(0, -1, 1//10)
        p₅ = Point3D{T}(0,  2, 1//10)
        l = LineSegment(p₄, p₅)
        hit, point = intersect(l, tri)
        @test !hit

        # no intersection non-coplanar
        p₄ = Point3D{T}(-1, 1//10, -1//10)
        p₅ = Point3D{T}( 1, 1//10, -1//10)
        l = LineSegment(p₄, p₅)
        hit, point = intersect(l, tri)
        @test !hit

        # no intersection coplanar
        p₄ = Point3D{T}(0, -1, 1)
        p₅ = Point3D{T}(0, -1, 0)
        l = LineSegment(p₄, p₅)
        hit, point = intersect(l, tri)
        @test !hit

        # intersects on boundary of triangle
        p₄ = Point3D{T}(-1, 0, 0)
        p₅ = Point3D{T}( 1, 0, 0)
        l = LineSegment(p₄, p₅)
        hit, point = intersect(l, tri)
        @test hit
        @test point ≈ Point3D{T}(0, 0, 0)

                # 4 intersections
        l = LineSegment2D(p₃, p₁)
        hit, points = intersect(l, quad)
        @test hit
        @test points[1] ≈ p₁
        @test points[2] ≈ p₃

        # 2 intersections
        l = LineSegment2D(Point2D{T}(0, 1//2), Point2D{T}(1, 1//2))
        hit, points = intersect(l, quad)
        @test hit
        @test points[1] ≈ Point2D{T}(1, 1//2)
        @test points[2] ≈ Point2D{T}(0, 1//2)

        # 0 intersections
        l = LineSegment2D(Point2D{T}(-1, -1), Point2D{T}(2, -1))
        hit, points = intersect(l, quad)
        @test !hit

                    # intersect
            # 0 intersection
            l = LineSegment2D(Point2D{T}(0, -1), Point2D{T}(4, -1))
            n, points = l ∩ tri6
            @test n == 0

            # 2 intersection
            l = LineSegment2D(Point2D{T}(0, 1), Point2D{T}(4, 1))
            n, points = l ∩ tri6
            @test n == 2
            @test points[1] ≈ Point2D{T}(3, 1)
            @test points[2] ≈ Point2D{T}(1, 1)

            # 4 intersection
            l = LineSegment2D(Point2D{T}(0, 1//10), Point2D{T}(4, 1//10))
            n, points = l ∩ tri6
            @test n == 4
            @test points[1] ≈ Point2D{T}(0.4254033307585166, 1//10)
            @test points[2] ≈ Point2D{T}(1.9745966692414834, 1//10)
            @test points[3] ≈ Point2D{T}(2.1900000000000000, 1//10)
            @test points[4] ≈ Point2D{T}(1//10,              1//10)

            # 6 intersection
            p₁ = Point2D{T}( 1, 0)
            p₂ = Point2D{T}( 0, 0)
            p₃ = Point2D{T}(-1, 0)
            p₄ = Point2D{T}( 1//2, -1//2)
            p₅ = Point2D{T}(-1//2, -1//2)
            p₆ = Point2D{T}(  0,    -2)
            tri6 = QuadraticTriangle(p₁, p₂, p₃, p₄, p₅, p₆)
            l = LineSegment2D(Point2D{T}(-2, -1//4), Point2D{T}(2, -1//4))
            n, points = l ∩ tri6
            @test n == 6
            @test points[1] ≈ Point2D{T}( 0.14644659, -1//4)
            @test points[2] ≈ Point2D{T}( 0.8535534,  -1//4)
            @test points[3] ≈ Point2D{T}(-0.8535534,  -1//4)
            @test points[4] ≈ Point2D{T}(-0.14644665, -1//4)
            @test points[5] ≈ Point2D{T}( 0.9354143,  -1//4)
            @test points[6] ≈ Point2D{T}(-0.9354143,  -1//4)

        # intersect
        # 0 intersections
        l = LineSegment2D(Point2D{T}(0, -1), Point2D{T}(4, -1))
        npoints, points = l ∩ quad8
        @test npoints == 0

        # 2 intersections
        l = LineSegment2D(Point2D{T}(-1, 1), Point2D{T}(4, 1))
        npoints, points = l ∩ quad8
        @test npoints == 2
        @test points[1] ≈ Point2D{T}(2.44444444444, 1)
        @test points[2] ≈ Point2D{T}(0, 1)

        # 4 intersections
        l = LineSegment2D(Point2D{T}(0, 1//10), Point2D{T}(4, 1//10))
        npoints, points = l ∩ quad8
        @test npoints == 4
        @test points[1] ≈ Point2D{T}(0.20557280900008415,       1//10)
        @test points[2] ≈ Point2D{T}(1.9944271909999158,        1//10)
        @test points[3] ≈ Point2D{T}(2.0644444444444447,        1//10)
        @test points[4] ≈ Point2D{T}(0,                         1//10)

        # 6 intersections
        p₁ = Point2D{T}( 1, 0)
        p₂ = Point2D{T}( 0, 0)
        p₃ = Point2D{T}(-1, 0)
        p₄ = Point2D{T}( 0, -2)
        p₅ = Point2D{T}( 1//2, -1//2)
        p₆ = Point2D{T}(-1//2, -1//2)
        p₇ = Point2D{T}(-1//2, -3//2)
        p₈ = Point2D{T}( 1//2, -3//2)
        quad8 = QuadraticQuadrilateral(p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈)
        l = LineSegment2D(Point2D{T}(-2, -1//4), Point2D{T}(2, -1//4))
        n, points = l ∩ quad8
        @test n == 6
        @test points[1] ≈ Point2D{T}(T[ 0.14644659, -1//4])
        @test points[2] ≈ Point2D{T}(T[ 0.8535534,  -1//4])
        @test points[3] ≈ Point2D{T}(T[-0.8535534,  -1//4])
        @test points[4] ≈ Point2D{T}(T[-0.14644665, -1//4])
        @test points[5] ≈ Point2D{T}(T[-0.9354143,  -1//4])
        @test points[6] ≈ Point2D{T}(T[ 0.9354143,  -1//4])

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
export findcontent
################################################################
# Content Extractor
################################################################
# Documents
################################################################
function findcontent(::Type{T}, name::String, doc::Document, argAmlType::Type{<:AbsNode}) where {T}
    findcontent(T, name, root(doc), argAmlType)
end

# if no type is provided consider it to be Vector{Union{String, Nothing}}
function findcontent(name::String, doc::Document, argAmlType::Type{<:AbsNode})
    findcontent(Vector{Union{String, Nothing}}, name, root(doc), argAmlType)
end
################################################################
# Nodes
################################################################
# Single extraction

"""

    findcontent(type, element::String, node, argAmlType)

Finds all the elements with the address of string in the node and converts to the passed type.

    findcontent(element::String, node, argAmlType)

If no type is provided consider it to be Vector{Union{String, Nothing}}

```julia
findcontent("instrument-name",node, AbsNormal)
findcontent(UInt8,"midi-channel",node, AbsNormal)
```
"""
function findcontent(::Type{T}, name::String, node::Node, argAmlType::Type{<:AbsNode}) where {T}
    elm = findfirst(name, node, argAmlType)
    if isnothing(elm) # return nothing if nothing is found
        return nothing
    else
        return nodeparse(T, elm)
    end
end

# Union with Nothing
@transform function findcontent(::Type{UN{T}}, name::String, node::Node, argAmlType::Type{allsubtypes(AbsDocOrNode)}) where {T}
    return findcontent(T, name, node, argAmlType)
end

# Nothing Alone
@transform function findcontent(::Type{Nothing}, name::String, node::Node, argAmlType::Type{allsubtypes(AbsDocOrNode)})
    return nothing
end

################################################################
# Vector extraction

function findcontent(::Type{Vector{T}},  name::String, node::Node, argAmlType::Type{<:Union{AbsNormal, AbsText}}) where{T}
    elms = findall(name, node, argAmlType)
    if isnothing(elms) # return nothing if nothing is found
        return nothing
    else
        return nodeparse(T, elms)
    end
end

# Fast attribute extractor
function findcontent(::Type{Vector{T}},  name::String, node::Node, argAmlType::Type{<:Union{AbsAttribute}}) where{T}
    if haskey(node, name)
        elms = node[name]
        return nodeparse(T, elms)
    else  # return nothing if nothing is found
        return nothing
    end
end

# Union with Nothing
@transform function findcontent(::Type{Vector{UN{T}}}, name::String,node::Node, argAmlType::Type{allsubtypes(AbsDocOrNode)}) where {T}
    return findcontent(Vector{T},name,node, argAmlType)
end

# Nothing Alone
@transform function findcontent(::Type{Vector{Nothing}}, name::String,node::Node, argAmlType::Type{allsubtypes(AbsDocOrNode)})
    return nothing
end

# vector of Any - consider it to be string
@transform function findcontent(::Type{Vector{Any}}, name::String,node::Node, argAmlType::Type{allsubtypes(AbsDocOrNode)})
    return findcontent(Vector{String},name,node, argAmlType)
end

# if no type is provided consider it to be Vector{Union{String, Nothing}}
@transform function findcontent(name::String, node::Node, argAmlType::Type{allsubtypes(AbsDocOrNode)})
    return findcontent(Vector{Union{String, Nothing}}, name, node, argAmlType)
end

################################################################
# Dict Extraction
function findcontent(::Type{AbstractDict}, name, node, argAmlType)
    throw(MethodError("Dicts are not supported for extraction/updating"))
end

# TODO: For Dict findcontent try to have some assumptions and fix it


text:
const HybridMatrixLike{T} = Union{
    HybridMatrix{<:Any, <:Any, T},
    Transpose{T, <:HybridMatrix{T}},
    Adjoint{T, <:HybridMatrix{T}},
    Symmetric{T, <:HybridMatrix{T}},
    Hermitian{T, <:HybridMatrix{T}},
    Diagonal{T, <:HybridMatrix{<:Any, T}}
}

const HybridVecOrMatLike{T} = Union{HybridVector{<:Any, T}, HybridMatrixLike{T}}

# Binary ops
# Between arrays
@inline +(a::HybridArray, b::HybridArray) = a .+ b
@inline +(a::AbstractArray, b::HybridArray) = a .+ b
@inline +(a::StaticArray, b::HybridArray) = a .+ b
@inline +(a::HybridArray, b::AbstractArray) = a .+ b
@inline +(a::HybridArray, b::StaticArray) = a .+ b

@inline -(a::HybridArray, b::HybridArray) = a .- b
@inline -(a::AbstractArray, b::HybridArray) = a .- b
@inline -(a::StaticArray, b::HybridArray) = a .- b
@inline -(a::HybridArray, b::AbstractArray) = a .- b
@inline -(a::HybridArray, b::StaticArray) = a .- b

# Scalar-array
@inline *(a::Number, b::HybridArray) = a .* b
@inline *(a::HybridArray, b::Number) = a .* b

@inline /(a::HybridArray, b::Number) = a ./ b
@inline \(a::Number, b::HybridArray) = a .\ b

@inline vcat(a::HybridVecOrMatLike) = a
@inline vcat(a::HybridVecOrMatLike, b::HybridVecOrMatLike) = _vcat(Size(a), Size(b), a, b)
@inline vcat(a::HybridVecOrMatLike, b::HybridVecOrMatLike, c::HybridVecOrMatLike...) = vcat(vcat(a,b), vcat(c...))

@generated function _vcat(::Size{Sa}, ::Size{Sb}, a::HybridVecOrMatLike, b::HybridVecOrMatLike) where {Sa, Sb}
    if Size(Sa)[2] != Size(Sb)[2]
        throw(DimensionMismatch("Tried to vcat arrays of size $Sa and $Sb"))
    end

    if a <: HybridVector && b <: HybridVector
        Snew = (Sa[1] + Sb[1],)
    else
        Snew = (Sa[1] + Sb[1], Size(Sa)[2])
    end

    return quote
        Base.@_inline_meta
        @inbounds return _h_similar_type($a, promote_type(eltype(a), eltype(b)), Size($Snew))(vcat(a.data, b.data))
    end
end

@inline hcat(a::HybridVecOrMatLike, b::HybridVecOrMatLike) = _hcat(Size(a), Size(b), a, b)
@inline hcat(a::HybridVecOrMatLike, b::HybridVecOrMatLike, c::HybridVecOrMatLike...) = hcat(hcat(a,b), c...)

# used in _vcat and _hcat
@inline +(::Dynamic, ::Dynamic) = Dynamic()

@generated function _hcat(::Size{Sa}, ::Size{Sb}, a::HybridVecOrMatLike, b::HybridVecOrMatLike) where {Sa, Sb}
    if Sa[1] != Sb[1]
        throw(DimensionMismatch("Tried to hcat arrays of size $Sa and $Sb"))
    end

    Snew = (Sa[1], Size(Sa)[2] + Size(Sb)[2])

    return quote
        Base.@_inline_meta
        return _h_similar_type($a, promote_type(eltype(a), eltype(b)), Size($Snew))(hcat(a.data, b.data))
    end
end


text:
using Test, DiffEqOperators
using OrdinaryDiffEq

@testset "Parabolic Heat Equation with Dirichlet BCs" begin
    x = collect(-pi : 2pi/511 : pi)
    u0 = @. -(x - 0.5).^2 + 1/12
    @. u_analytic(x)  = -(x .- 0.5).^2 + 1/12
    A = CenteredDifference(2,2,2π/511,512);
    bc = DirichletBC(u_analytic(-pi - 2pi/511),u_analytic(pi + 2pi/511))
    step = (u,p,t) ->A*bc*u
    heat_eqn = ODEProblem(step, u0, (0.,1.))
    soln = solve(heat_eqn,Tsit5())

    for t in 0:0.1:1
        @test soln(t)[1] ≈ u0[1] rtol=0.05
        @test soln(t)[end] ≈ u0[end] rtol=0.05
    end

    # UpwindDifference with equal no. of primay wind and offside points should behave like a CenteredDifference
    A2 = UpwindDifference(2,1,2π/511,512,1,offside=1);
    step = (u,p,t) ->A2*bc*u
    heat_eqn = ODEProblem(step, u0, (0.,1.))
    soln = solve(heat_eqn,Tsit5())

    for t in 0:0.1:1
        @test soln(t)[1] ≈ u0[1] rtol=0.05
        @test soln(t)[end] ≈ u0[end] rtol=0.05
    end
end

@testset "Parabolic Heat Equation with Neumann BCs" begin
    N = 128
    dx = 2π/(N-1)
    x = collect(-pi : dx : pi)
    u0 = @. -(x - 0.5)^2 + 1/12
    B = CenteredDifference(1,2,dx,N-2)
    deriv_start, deriv_end = (B*u0)[1], (B*u0)[end]

    A = CenteredDifference(2,2,dx,N)
    bc = NeumannBC((deriv_start,deriv_end),dx,1)

    step = (u,p,t) ->A*bc*u
    heat_eqn = ODEProblem(step, u0, (0.,1.))
    soln = solve(heat_eqn,Tsit5())

    first_order_coeffs_start = [-11/6, 3.0, -3/2, 1/3] * (1/dx)
    first_order_coeffs_end = -reverse([-11/6, 3.0, -3/2, 1/3] * (1/dx))

    for t in 0:0.1:1
        @test sum(first_order_coeffs_start .* soln(t)[1:4]) ≈ deriv_start atol=1e-1
        @test sum(first_order_coeffs_end .* soln(t)[end-3:end]) ≈ deriv_end atol=1e-1
    end

    # UpwindDifference with equal no. of primay wind and offside points should behave like a CenteredDifference
    A2 = UpwindDifference(2,1,dx,N,1,offside=1)
    B2 = UpwindDifference(1,2,dx,N-2,1,offside=1)
    deriv_start, deriv_end = (B2*u0)[1], (B2*u0)[end]
    bc = NeumannBC((deriv_start,deriv_end),dx,1)

    step = (u,p,t) ->A2*bc*u
    heat_eqn = ODEProblem(step, u0, (0.,1.))
    soln = solve(heat_eqn,Tsit5())

    for t in 0:0.1:1
        @test sum(first_order_coeffs_start .* soln(t)[1:4]) ≈ deriv_start atol=1e-1
        @test sum(first_order_coeffs_end .* soln(t)[end-3:end]) ≈ deriv_end atol=1e-1
    end
    # Testing for 2 offside points against Standard Vector input
    B3 = UpwindDifference(1,4,dx,N-2,-1,offside=2)
    deriv_start, deriv_end = (-1*B3*u0)[1], (-1*B3*u0)[end]
    bc = NeumannBC((deriv_start,deriv_end),dx,1)

    step = (u,p,t) ->A2*bc*u
    heat_eqn = ODEProblem(step, u0, (0.,1.))
    soln = solve(heat_eqn,Tsit5())

    for t in 0:0.1:1
        @test sum(first_order_coeffs_start .* soln(t)[1:4]) ≈ deriv_start atol=1e-1
        @test sum(first_order_coeffs_end .* soln(t)[end-3:end]) ≈ deriv_end atol=1e-1
    end
end

@testset "Parabolic Heat Equation with Robin BCs" begin
    N = 128
    dx = 2π/(N-1)
    x = collect(-pi : dx : pi)
    u0 = @. -(x - 0.5)^2 + 1/12
    B = CenteredDifference(1,2,dx,N-2);
    deriv_start, deriv_end = (B*u0)[1], (B*u0)[end]
    params = [1.0,0.5]

    left_RBC = params[1]*u0[1] - params[2]*deriv_start
    right_RBC = params[1]*u0[end] + params[2]*deriv_end

    A = CenteredDifference(2,2,dx,N);
    bc = RobinBC((params[1],-params[2],left_RBC), (params[1],params[2],right_RBC),dx,1);
    step1(u,p,t)=A*bc*u
    heat_eqn = ODEProblem(step1, u0, (0.,1.));
    soln = solve(heat_eqn,Rodas4(autodiff=false),reltol=1e-6);

    first_order_coeffs_start = [-11/6, 3.0, -3/2, 1/3] * (1/dx)
    first_order_coeffs_end = -reverse([-11/6, 3.0, -3/2, 1/3] * (1/dx))
    val = []

    for t in 0.2:0.1:1.0
        @test params[1]*soln(t)[1] - params[2]*sum(first_order_coeffs_start .* soln(t)[1:4]) ≈ left_RBC atol=4e-1
        # append!(val,params[1]*soln(t)[1] + -params[2]*sum(first_order_coeffs_start .* soln(t)[1:4]) - left_RBC)
        @test params[1]*soln(t)[end] + params[2]*sum(first_order_coeffs_end .* soln(t)[end-3:end]) ≈ right_RBC atol=4e-1
    end

    # UpwindDifference with equal no. of primay wind and offside points should behave like a CenteredDifference
    A2 = UpwindDifference(2,1,dx*ones(N+1),N,1,offside=1)
    B2 = UpwindDifference(1,2,dx*ones(N-1),N-2,1,offside=1)
    deriv_start, deriv_end = (B2*u0)[1], (B2*u0)[end]
    left_RBC = params[1]*u0[1] - params[2]*deriv_start
    right_RBC = params[1]*u0[end] + params[2]*deriv_end
    bc = RobinBC((params[1],-params[2],left_RBC), (params[1],params[2],right_RBC),dx,1);

    step2(u,p,t)=A2*bc*u
    heat_eqn = ODEProblem(step2, u0, (0.,1.));
    soln = solve(heat_eqn,Rodas4(autodiff=false),reltol=1e-6);

    for t in 0.2:0.1:1.0
        @test params[1]*soln(t)[1] - params[2]*sum(first_order_coeffs_start .* soln(t)[1:4]) ≈ left_RBC atol=4e-1
        @test params[1]*soln(t)[end] + params[2]*sum(first_order_coeffs_end .* soln(t)[end-3:end]) ≈ right_RBC atol=4e-1
    end
    # Testing for 2 offside points against Standard Vector input
    B3 = UpwindDifference(1,4,dx*ones(N-1),N-2,-1,offside=2)
    deriv_start, deriv_end = (-1*B3*u0)[1], (-1*B3*u0)[end]
    left_RBC = params[1]*u0[1] - params[2]*deriv_start
    right_RBC = params[1]*u0[end] + params[2]*deriv_end
    bc = RobinBC((params[1],-params[2],left_RBC), (params[1],params[2],right_RBC),dx,1);

    step3(u,p,t)=A2*bc*u
    heat_eqn = ODEProblem(step3, u0, (0.,1.));
    soln = solve(heat_eqn,Rodas4(autodiff=false),reltol=1e-6);

    for t in 0.2:0.1:1.0
        @test params[1]*soln(t)[1] - params[2]*sum(first_order_coeffs_start .* soln(t)[1:4]) ≈ left_RBC atol=4e-1
        @test params[1]*soln(t)[end] + params[2]*sum(first_order_coeffs_end .* soln(t)[end-3:end]) ≈ right_RBC atol=4e-1
    end
end


text:
# this file will test that both algorithms correctly infer parameters from a normal distribution

using Plots
using Statistics

function getnormal(params, constants, targetdata)
  simdata = rand(Normal(params...), 1000)
  ApproxBayes.ksdist(simdata, targetdata), 1
end

println("Test parameters of normal distribution are inferred correctly (mean within 5% of true value)")

Random.seed!(1)
p1 = 2.0
p2 = 0.4
targetdata = rand(Normal(p1, p2), 1000)

#setup ABC alogrithm specifications for Rejection algorithm
setup = ABCRejection(getnormal,
  2,
  0.1,
  Prior([Uniform(0.0, 20.0), Uniform(0.0, 2.0)]);
  maxiterations = 10^6,
  )
# run ABC inference
resrejection = runabc(setup, targetdata);

println("\t Check ABC rejection algorithm correctly infers parameters")

# test that mean value of posterior is within 10% of true value
@test isapprox(mean(resrejection.parameters, dims = 1)[1], p1, rtol = 0.05)
@test isapprox(mean(resrejection.parameters, dims = 1)[2], p2, rtol = 0.05)

println("\t Check no errors arising from plotting")
plot(resrejection)

println("\t Check no errors arising from printing results\n")
println("#########################################")
show(resrejection)
println("#########################################\n")

println("\t Check output is written to file")
writeoutput(resrejection)
@test isfile("Rejection-output.txt")
rm("Rejection-output.txt")

#now run SMC algorithm
setup = ABCSMC(getnormal,
  2,
  0.05,
  Prior([Uniform(0.0, 20.0), Uniform(0.0, 2.0)]),
  )
ressmc = runabc(setup, targetdata, verbose = false);

println("\t Check ABC SMC algorithm correctly infers parameters")
# test that mean value of posterior is within 5% of true value
@test isapprox(StatsBase.mean(ressmc.parameters, weights(ressmc.weights), 1)[1], p1, rtol = 0.05)
@test isapprox(StatsBase.mean(ressmc.parameters, weights(ressmc.weights), 1)[2], p2, rtol = 0.05)

#test weights sum to 1
@test isapprox(sum(ressmc.weights), 1.0, rtol = 0.0001)

println("\t Check no errors arising from plotting")
plot(ressmc)

println("\t Check no errors arising from printing results\n")
println("#########################################")
show(ressmc)
println("#########################################\n")

println("\t Check output is written to file")
writeoutput(ressmc)
@test isfile("SMC-output.txt")
rm("SMC-output.txt")

setup = ABCSMC(getnormal,
  2,
  0.05,
  Prior([Uniform(0.0, 20.0), Uniform(0.0, 2.0)]),
  kernel = ApproxBayes.gaussiankernel
  )
println("\t Check ABC SMC algorithm correctly infers parameters with gaussian perturbation kernel")
@test isapprox(StatsBase.mean(ressmc.parameters, weights(ressmc.weights), 1)[1], p1, rtol = 0.05)
@test isapprox(StatsBase.mean(ressmc.parameters, weights(ressmc.weights), 1)[2], p2, rtol = 0.05)

#test SMC is more efficient than rejection algorithm
println("\t Check ABC SMC is more efficient than ABC rejection")
setup = ABCSMC(getnormal,
  2,
  0.1,
  Prior([Uniform(0.0, 20.0), Uniform(0.0, 2.0)]),
  )
ressmc = runabc(setup, targetdata, verbose = false);
@test ressmc.accratio/resrejection.accratio > 1.0

#Check fallback to taking 100 particles with shortest distance works
setup = ABCRejection(getnormal,
  2,
  0.1,
  Prior([Uniform(0.0, 20.0), Uniform(0.0, 2.0)]);
  maxiterations = 10^3,
  )
# run ABC inference
resrejection = runabc(setup, targetdata);

@test resrejection.numsims == setup.maxiterations

#check results with gaussian kernel and non uniform prior
setup = ABCSMC(getnormal,
  2,
  0.05,
  Prior([Uniform(0.0, 20.0), Gamma(0.5, 0.5)]),
  kernel = gaussiankernel
  )
ressmc = runabc(setup, targetdata, verbose = false);
println(ressmc)

println("\t Check ABC SMC algorithm with gaussian perturbation kernel correctly infers parameters")
# test that mean value of posterior is within 5% of true value
@test isapprox(StatsBase.mean(ressmc.parameters, weights(ressmc.weights), 1)[1], p1, rtol = 0.05)
@test isapprox(StatsBase.mean(ressmc.parameters, weights(ressmc.weights), 1)[2], p2, rtol = 0.05)


text:
module NCTSSOS

using DynamicPolynomials
using MultivariatePolynomials
using JuMP
using MosekTools
using Graphs
using ChordalGraph
using MetaGraphs
using LinearAlgebra
using SparseArrays

export nctssos_first, nctssos_higher!, cs_nctssos_first, cs_nctssos_higher!

include("clique_merge.jl")
include("ncupop.jl")
include("nccpop.jl")
include("mixncpop.jl")

end


text:
module Day09
export day09

const ADJACENT = [CartesianIndex(1,0),CartesianIndex(-1,0),CartesianIndex(0,1),CartesianIndex(0,-1)]

parseLine(line) = parse.(Int, Iterators.take(line, length(line)))

function readInput(inputfile="inputs/input09.txt")
    return hcat([parseLine(line) for line in eachline(inputfile)]...)
end

function day09(heightmap = readInput("inputs/input09.txt"))
    low_points = _find_low_points(heightmap)
    return (day09Part1(heightmap, low_points), day09Part2!(heightmap, low_points))
end

function day09Part1(heightmap, low_points)
    # return sum(heightmap[low_points]) + length(low_points)
    # for some reason the broadcasted version is slightly faster
    # multi-index: 951.625 ns (2 allocations: 3.02 KiB)
    # broadcasted: 816.955 ns (3 allocations: 3.03 KiB)
    return sum(getindex.(Ref(heightmap), low_points)) + length(low_points)
end

function _find_low_points(heightmap)
    low_points = CartesianIndex{2}[]
    for I in CartesianIndices(heightmap)
        h = heightmap[I]
        islow = true
        for adj in ADJACENT
            if checkbounds(Bool, heightmap, I+adj) && heightmap[I+adj] < h
                islow = false
                break
            end
        end
        islow && push!(low_points, I)
        # adj_I = filter!(x->checkbounds(Bool, heightmap, x), Ref(I) .+ ADJACENT)
        # if all(heightmap[adj_I] .> h)
        #     push!(low_points, I)
        # end
    end
    return low_points
end

function day09Part2!(heightmap, low_points)
    basins = _flood!.(Ref(heightmap), low_points)
    return reduce(*, sort!(basins)[end-2:end])
end

# recursive destructive floodfill
function _flood!(heightmap, at)
    heightmap[at] = 9
    filled = 1
    for adj in ADJACENT
        new = at + adj
        if checkbounds(Bool, heightmap, new) && heightmap[new] < 9
            filled += _flood!(heightmap, new)
        end
    end
    return filled
end

## iterative desctructive floodfill with stack
function day09Part2!_iter(heightmap, low_points)
    # floodfill starting from each low point
    basins = Int[]
    for low_point in low_points
        l = 0
        todo = [low_point]
        while length(todo) > 0
            at = pop!(todo)
            heightmap[at] < 9 || continue # already seen
            heightmap[at] = 9# mark as seen
            l += 1
            for adj in ADJACENT
                new = at + adj
                if checkbounds(Bool, heightmap, new) && heightmap[new] < 9
                    push!(todo, new)
                end
            end
        end
        push!(basins, l)
    end
    return reduce(*, sort!(basins)[end-2:end])
end

## iterative non-destructive floodfill
function day09Part2(heightmap, low_points)
    # floodfill starting from each low point
    basins = Int[]
    for low_point in low_points
        todo = [low_point]
        seen = Set{CartesianIndex{2}}()
        while length(todo) > 0
            at = pop!(todo)
            for adj in ADJACENT
                new = at + adj
                if !in(new, seen) && checkbounds(Bool, heightmap, new) && heightmap[new] < 9
                    push!(todo, new)
                end
            end
            push!(seen, at) # mark as seen
        end
        push!(basins, length(seen))
    end
    return reduce(*, sort!(basins)[end-2:end])
end

end #module


text:
function jsrender(session::Session, md::Markdown.MD)
    md_div = DOM.div(md_html(session, md.content)...; class="markdown-body")
    return JSServe.jsrender(session, md_div)
end

function md_html(session::Session, content::Vector)
    return md_html.((session,), content)
end

function md_html(session::Session, header::Markdown.Header{l}) where l
    return [DOM.um("h$l", htmlinline(session, header.text)...)]
end

function md_html(session::Session, code::Markdown.Code)
    maybe_lang = !isempty(code.language) ? Any[:class=>"language-$(code.language)"] : []
    return [DOM.um("pre", DOM.um("code", code.code; maybe_lang...))]
end

function md_html(session::Session, md::Markdown.Paragraph)
    return [DOM.p(htmlinline(session, md.content)...)]
end

function md_html(session::Session, md::Markdown.BlockQuote)
    return [DOM.um("blockquote", md_html(session, md.content)...)]
end

function md_html(session::Session, f::Markdown.Footnote)
    return [DOM.div(DOM.p(f.id; class="footnote-title"), md_html(session, f.text)...;
                          class="footnote", id="footnote-$(f.id)")]
end

function md_html(session::Session, md::Markdown.Admonition)
    title = DOM.p(md.title; class="admonition-title")
    return [DOM.div(title; class="admonition $(md.category)")]
end

function md_html(session::Session, md::Markdown.List)
    maybe_attr = md.ordered > 1 ? Any[:start => string(md.ordered)] : []
    tag = Markdown.isordered(md) ? "ol" : "ul"
    return [DOM.um(
        tag, maybe_attr...,
        map(md.items) do item
            DOM.um("li", md_html(session, item)...)
        end...
    )]
end

function md_html(session::Session, md::Markdown.HorizontalRule)
    return [DOM.m_unesc("hr")]
end

function md_html(session::Session, md::Markdown.Table)
    content = map(enumerate(md.rows)) do (i, row)
        tr_content = map(enumerate(md.rows[i])) do (j, c)
            alignment = md.align[j]
            alignment = alignment === :l ? "left" : alignment === :r ? "right" : "center"
            return DOM.um(i == 1 ? "th" : "td", htmlinline(session, c)...; align=alignment)
        end
        return DOM.um("tr", tr_content...)
    end
    return [DOM.um("table", content...)]
end

function htmlinline(session::Session, content::Vector)
    return htmlinline.((session,), content)
end

function htmlinline(session::Session, code::Markdown.Code)
    return [DOM.um("code", code.code)]
end

function htmlinline(session::Session, md::Union{Symbol, AbstractString})
    return [md]
end

function htmlinline(session::Session, md::Markdown.Bold)
    return [DOM.um("strong", htmlinline(session, md.text)...)]
end

function htmlinline(session::Session, md::Markdown.Italic)
    return [DOM.um("em", htmlinline(session, md.text)...)]
end

function htmlinline(session::Session, md::Markdown.Image)
    return [DOM.um("img"; src=md.url, alt=md.alt)]
end

function htmlinline(session::Session, f::Markdown.Footnote)
    return [DOM.um("a", string("[", f.id, "]"); href="#footnote-$(f.id)", class="footnote")]
end

function htmlinline(session::Session, link::Markdown.Link)
    return [DOM.um("a", href = link.url, htmlinline(session, link.text)...)]
end

function htmlinline(session::Session, br::Markdown.LineBreak)
    return [DOM.um("br")]
end

htmlinline(session::Session, x) = [jsrender(session, x)]
md_html(session::Session, x) = [jsrender(session, x)]

function render_result(result)
    if Tables.istable(result)
        return Table(result)
    else
        return DOM.span(result)
    end
end

"""
    contextual_eval(parent, expr)
"Evals" expression without eval by only allowing getfield + getindex expressions

```
contextual_eval(context, :(lala.blalba)) == context.lala.blabla
contextual_eval(context, :(lala.blalba[1])) == context.lala.blabla[1]
contextual_eval(context, :(julia_func())) == error
```
"""
function contextual_eval(parent, expr)
    if expr isa Symbol
        return getfield(parent, expr)
    else
        if expr.head == :(.)
            return getfield(contextual_eval(parent, expr.args[1]), expr.args[2].value)
        elseif expr.head == :ref
            getindex(contextual_eval(parent, expr.args[1]), expr.args[2])
        else
            return contextual_eval(parent, expr.args[1])
        end
    end
end

"""
    replace_expressions(markdown, context)
Replaces all expressions inside `markdown` savely, by only supporting
getindex/getfield expression that will index into `context`
"""
function replace_expressions(markdown, context; eval_julia_code=false)
    if markdown isa Union{Expr, Symbol}
        result = contextual_eval(context, markdown)
        return render_result(result)
    end
    if hasproperty(markdown, :content)
        markdown.content .= replace_expressions.(markdown.content, (context,); eval_julia_code=eval_julia_code)
    elseif hasproperty(markdown, :text)
        markdown.text .= replace_expressions.(markdown.text, (context,); eval_julia_code=eval_julia_code)
    end
    return markdown
end

function parseall(str)
    pos = firstindex(str)
    exs = []
    while pos <= lastindex(str)
        ex, pos = Meta.parse(str, pos)
        push!(exs, ex)
    end
    if length(exs) == 0
        throw(Base.ParseError("end of input"))
    elseif length(exs) == 1
        return exs[1]
    else
        return Expr(:block, exs...)
    end
end

function replace_interpolation!(context, expr)
    return expr
end

function replace_interpolation!(context, expr::Expr)
    if expr.head == :$
        return contextual_eval(context, expr.args[1])
    else
        expr.args .= replace_interpolation!.((context,), expr.args)
        return expr
    end
end

function replace_expressions(markdown::Markdown.Code, context; eval_julia_code=false)
    if markdown.language == "julia" && eval_julia_code isa Module
        run = Button(">")
        result = Observable{Any}(DOM.span(""))
        on(run) do click
            expr = parseall(markdown.code)
            expr = replace_interpolation!(context, expr)
            evaled = eval_julia_code.eval(expr)
            result[] = render_result(evaled)
        end
        return md"""
        $(markdown)
        $(run)
        $(result)
        """
    else
        return markdown
    end
end

"""
    string_to_markdown(source::String, context; eval_julia_code=false)

Replaces all interpolation expressions inside `markdown` savely, by only supporting
getindex/getfield expression that will index into `context`.
You can eval Julia code blocks by setting `eval_julia_code` to a Module, into which
the code gets evaluated!
"""
function string_to_markdown(source::String, context; eval_julia_code=false)
    markdown = Markdown.parse(source)
    return replace_expressions(markdown, context; eval_julia_code=eval_julia_code)
end


text:
module CombOpt

    using DelimitedFiles

    export Node, Edge, Graph, adj
    export read_graph
    export dijkstra
    export bellman_ford

    include("graph.jl")
    include("forward_star.jl")
    include("utils.jl")
    include("dijkstra.jl")
    include("bellman_ford.jl")

end


text:
using DistributedObjectStores
using Test

@testset "DistributedObjectStores.jl" begin
    # Write your tests here.
end


text:
using ConfusionMatrix
using Test
using Flux: onehotbatch, softmax

@testset "ConfusionMatrix.jl" begin
    @testset "ConfMat construction" begin

# ConfMat([:truc, :bidule], [:truc, :truc, :truc, :bidule], [:truc :truc :bidule :bidule])
# ConfMat([:truc, :bidule, :machin], [:truc, :truc, :truc, :bidule], [:truc :truc :bidule :machin])
# ConfMat([0 1], [0, 1, 1, 0], [0, 1, 1, 0])

        labels = [:dog, :cat]
        y = rand(labels, 100)
        ŷ = rand(labels, 100)
        cm = ConfMat(labels, y, ŷ)
        @test size(cm.counts) == (2, 2)

        labels = [:dog, :cat, :frog]
        y = rand(labels, 100)
        ŷ = rand(labels, 100)
        cm = ConfMat(labels, y, ŷ)
        @test size(cm.counts) == (3, 3)

        labels = [1, 2, 3, 4, 5]
        y = rand(labels, 100)
        ŷ = rand(labels, 100)
        cm = ConfMat(labels, y, ŷ)
        @test size(cm.counts) == (5, 5)

        # Simulate the ouput of a Flux binary classifier
        labels = [:dog, :cat]
        cm = ConfMat(labels)
        nepochs = 10
        batchsize = 64
        for e ∈ 1:nepochs
            y = onehotbatch(rand(labels, batchsize), labels)
            ŷ = rand(1, batchsize)
            update!(cm, y, ŷ)
        end
        @test sum(cm.counts) == nepochs * batchsize

        # Simulate the ouput of a Flux multiclass classifier
        labels = [:dog, :cat, :frog, :bear]
        cm = ConfMat(labels)
        nepochs = 10
        batchsize = 64
        for e ∈ 1:nepochs
            y = onehotbatch(rand(labels, batchsize), labels)
            ŷ = rand(4, batchsize) |> softmax
            update!(cm, y, ŷ)
        end
        @test sum(cm.counts) == nepochs * batchsize

        io = IOBuffer()
        show(io, cm)
        @test startswith(String(take!(io)), "ConfMat{4,Symbol}")
    end
end


text:
# Interface that needs to be overwritten by backend
# Slightly difference behavior from buffer, since not all blas backends work directly with
# the gpu array buffer
function blas_module(A)
    error("$(typeof(A)) doesn't support BLAS operations")
end
function blasbuffer(A)
    error("$(typeof(A)) doesn't support BLAS operations")
end

for elty in (Float32, Float64, ComplexF32, ComplexF64)
    T = VERSION >= v"1.3.0-alpha.115" ? :(Union{($elty), Bool}) : elty
    @eval begin
        function BLAS.gemm!(
                transA::AbstractChar, transB::AbstractChar, alpha::$T,
                A::GPUVecOrMat{$elty}, B::GPUVecOrMat{$elty},
                beta::$T, C::GPUVecOrMat{$elty}
            )
            blasmod = blas_module(A)
            result = blasmod.gemm!(
                transA, transB, alpha,
                blasbuffer(A), blasbuffer(B), beta, blasbuffer(C)
            )
            C
        end
    end
end

for elty in (Float64, Float32)
    @eval begin
        function BLAS.scal!(
                n::Integer, DA::$elty,
                DX::GPUArray{$elty, N}, incx::Integer
            ) where N
            blasmod = blas_module(DX)
            blasmod.scal!(n, DA, blasbuffer(DX), incx)
            DX
        end
    end
end

LinearAlgebra.rmul!(s::Number, X::GPUArray) = rmul!(X, s)
function LinearAlgebra.rmul!(X::GPUArray{T}, s::Number) where T <: BLAS.BlasComplex
    R = typeof(real(zero(T)))
    N = 2*length(X)
    buff = unsafe_reinterpret(R, X, (N,))
    BLAS.scal!(N, R(s), buff, 1)
    X
end
function LinearAlgebra.rmul!(X::GPUArray{T}, s::Number) where T <: Union{Float32, Float64}
    BLAS.scal!(length(X), T(s), X, 1)
    X
end

for elty in (Float32, Float64, ComplexF32, ComplexF64)
    T = VERSION >= v"1.3.0-alpha.115" ? :(Union{($elty), Bool}) : elty
    @eval begin
        function BLAS.gemv!(trans::AbstractChar, alpha::$T, A::GPUVecOrMat{$elty}, X::GPUVector{$elty}, beta::$T, Y::GPUVector{$elty})
            m, n = size(A, 1), size(A, 2)
            if trans == 'N' && (length(X) != n || length(Y) != m)
                throw(DimensionMismatch("A has dimensions $(size(A)), X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'C' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("A' has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'T' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("A.' has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            end
            blasmod = blas_module(A)
            blasmod.gemv!(
                trans, alpha,
                blasbuffer(A), blasbuffer(X), beta, blasbuffer(Y)
            )
            Y
        end
    end
end

for elty in (Float32, Float64, ComplexF32, ComplexF64)
    @eval begin
        function BLAS.axpy!(
                alpha::Number, x::GPUArray{$elty}, y::GPUArray{$elty}
            )
            if length(x) != length(y)
                throw(DimensionMismatch("x has length $(length(x)), but y has length $(length(y))"))
            end
            blasmod = blas_module(x)
            blasmod.axpy!($elty(alpha), blasbuffer(vec(x)), blasbuffer(vec(y)))
            y
        end
    end
end

for elty in (Float32, Float64, ComplexF32, ComplexF64)
    @eval begin
        function BLAS.gbmv!(trans::AbstractChar, m::Integer, kl::Integer, ku::Integer, alpha::($elty), A::GPUMatrix{$elty}, X::GPUVector{$elty}, beta::($elty), Y::GPUVector{$elty})
            n = size(A, 2)
            if trans == 'N' && (length(X) != n || length(Y) != m)
                throw(DimensionMismatch("A has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'C' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("A' has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'T' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("A.' has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            end
            blasmod = blas_module(A)
            blasmod.gbmv!(
                trans, m, kl, ku, alpha,
                blasbuffer(A), blasbuffer(X), beta, blasbuffer(Y)
            )
            Y
        end
    end
end


text:
using PyCall
@pyimport pygtsam as gtsam

#file = open("/home/andrej/factor_graph.txt")
file = open("/home/andrej/anpl/code/debugging/ORB_SLAM2_SERIALIZED_GRAPH.txt")

global graph_str, val_str
graph_str = "";
for line in eachline(file)
    #println(line)
    graph_str *= line;
    graph_str *= "\n";
end
close(file)

file = open("/home/andrej/anpl/code/debugging/ORB_SLAM2_SERIALIZED_VALUES.txt")
val_str = "";
for line in eachline(file)
    #println(line)
    val_str *= line;
    val_str *= "\n";
end
close(file)

fg = gtsam.deserializeGraph(graph_str)
fg[:printf]()
#show(IOContext(STDOUT, limit=false), "text/plain", keys(fg))

# get first factor
f0 = fg[:__getitem__](0)

#relPose=gtsam.BetweenFactorPose2;
#relPose = gtsam.getNonlinearFactor(fg, 1)
#z = relPose[:measured]()
#z[:x](), z[:y](), z[:theta]()

initialEst = gtsam.deserializeValues(val_str)
initialEst[:printf]()

# symbol 'x0'
s = gtsam.Symbol( UInt8('x'),0)

robot_poses = gtsam.extractPose3(initialEst)
arr = collect(robot_poses)
# extract the first pose in array, not neccesarily 'x0'
s = arr[1][1]
pose3 = arr[1][2]
s[:print]("")
s[:chr]()
s[:index]()
p = robot_poses[s]
t = p[:translation]()[:vector]()
R = p[:rotation]()[:matrix]()

landmark_pos = gtsam.extractPoint3(initialEst)

# TODO No to_python (by-value) converter found for C++ type: gtsam::Value'
#k = gtsam.extractKeys(initialEst)
#initialEst[:at](k[1])

isam = gtsam.ISAM2()
#isam[:update](fg, initialEst)


text:
for T in [Float32, Float64, Int, Int8, UInt8, N0f8]

    a = rand(T, 5, 5)
    i = AstroImage(a)
    for stretch in [
        logstretch,
        powstretch,
        sqrtstretch,
        squarestretch,
        asinhstretch,
        sinhstretch,
        powerdiststretch
    ]
        # Easiest way to precompile everything we need is just to call these functions.
        # They have no side-effects.
        imview(a; stretch)

        # And precompile on an astroimage
        imview(i; stretch)
    end
    TI = typeof(i)
    precompile(parent, (TI,))
    precompile(header, (TI,))
    precompile(wcs, (TI,))
    precompile(getindex, (TI, Symbol))
    precompile(getindex, (TI, String))
    precompile(getindex, (TI, Int))
    precompile(getindex, (TI, Int, Int))
    precompile(getindex, (TI, Vector{Int}))
    precompile(getindex, (TI, Vector{Bool}))
    precompile(getindex, (TI, Matrix{Bool}))
    precompile(setindex!, (TI, Matrix{Bool}))
    precompile(world_to_pix, (TI, Vector{Float64}))
    precompile(pix_to_world, (TI, Vector{Float64}))
end

# From trace-compile:
precompile(Tuple{typeof(AstroImages.imview), Array{Float64, 2}})
precompile(Tuple{AstroImages.var"#imview##kw", NamedTuple{(:clims, :stretch, :cmap, :contrast, :bias), Tuple{AstroImages.Percent, typeof(Base.identity), Symbol, Int64, Float64}}, typeof(AstroImages.imview), AstroImages.AstroImage{Float64, 2, Tuple{DimensionalData.Dimensions.X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, DimensionalData.Dimensions.Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Array{Float64, 2}, Tuple{DimensionalData.Dimensions.X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, DimensionalData.Dimensions.Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}}})
precompile(Tuple{AstroImages.var"#imview_colorbar##kw", NamedTuple{(:clims, :stretch, :cmap, :contrast, :bias), Tuple{AstroImages.Percent, typeof(Base.identity), Symbol, Int64, Float64}}, typeof(AstroImages.imview_colorbar), AstroImages.AstroImage{Float64, 2, Tuple{DimensionalData.Dimensions.X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, DimensionalData.Dimensions.Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Array{Float64, 2}, Tuple{DimensionalData.Dimensions.X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, DimensionalData.Dimensions.Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}}})


text:
module PerfEvents

  using StructIO
  import Base: start, next, done, getindex, endof, length, show

  export EventIterator

  @struct immutable perf_file_section
    offset::UInt64
    size::UInt64
  end
  perf_file_section() = perf_file_section(0,0)

  @struct immutable perf_file_header
    magic::UInt64
    size::UInt64
    attr_size::UInt64
    attrs::perf_file_section
    data::perf_file_section
    event_types::perf_file_section
    features::NTuple{32, UInt8}
  end

  @struct immutable perf_trace_event_type
    event_id::UInt64
    name::NTuple{64,UInt8}
  end

  immutable PerfFileHandle{T}
    io::T
    header::perf_file_header
  end

  function readmeta(io::IO)
    PerfFileHandle(io, unpack(io,perf_file_header))
  end
  readmeta(path::String) = readmeta(open(path))

  # Record structures
  @struct immutable perf_event_header
    event_type::UInt32
    misc::UInt16
    size::UInt16
  end

  @struct immutable tid_record
    pid::UInt32
    tid::UInt32
  end
  const itrace_start = tid_record

  @struct immutable aux_record
    aux_offset::UInt64
    aux_size::UInt64
    flags::UInt64
  end

  @struct immutable mmap_record
    start::UInt64
    len::UInt64
    pgoff::UInt64
  end

  @struct immutable mmap2_extra_record
    maj::UInt32
    min::UInt32
    ino::UInt64
    ino_generation::UInt64
    prot::UInt32
    flags::UInt32
  end

  @struct immutable time_conv_record
    time_shift::UInt64
    time_mult::UInt64
    time_zero::UInt64
  end

  # Iterator for `attrs`. This us used both by the kernel and the file format
  @struct immutable perf_event_attr
    ev_type::UInt32
    size::UInt32
    config::UInt64
    sample_period_freq::UInt64
    sample_type::UInt64
    read_type::UInt64
    flags::UInt64
    wakeup::UInt32
    bp_type::UInt32
    config1::UInt64
    config2::UInt64
    branch_sample_type::UInt64
    sample_regs_user::UInt64
    sample_stack_user::UInt32
    clockid::Int32
    sample_regs_intr::UInt64
    aux_watermark::UInt32
    reserved_2::UInt32
  end

  @struct immutable perf_file_attr
    attr::perf_event_attr
    ids::perf_file_section
  end

  function show(io::IO, attr::perf_event_attr)
      if attr.ev_type == PERF_TYPE_HARDWARE
          print(io, PERF_HW[attr.config])
      else
          println(io, "Unkown event type"); return
      end
      print(io, ", sample_freq=", attr.sample_period_freq)
      print(io, ", size: ", attr.size)
      print(io, ", sample_type: ",
        join(map(x->PERF_SAMPLE_TYPE[x],
          sort(collect(filter(kind->(kind & attr.sample_type) != 0,
            keys(PERF_SAMPLE_TYPE))))),'|'))
  end

  immutable AttrIterator
    handle::PerfFileHandle
  end

  endof(it::AttrIterator) =
    div(it.handle.header.attrs.size, it.handle.header.attr_size)
  length(it::AttrIterator) = endof(it)

  function getindex(it::AttrIterator, n)
    @assert 1 <= n <= endof(it)
    @assert sizeof(perf_file_attr) <= it.handle.header.attr_size
    seek(it.handle.io,
      it.handle.header.attrs.offset +
      (n-1)*it.handle.header.attr_size)
    unpack(it.handle.io, perf_file_attr)
  end

  start(it::AttrIterator) = 1
  next(it::AttrIterator, s) = (it[s], s+1)
  done(it::AttrIterator, s) = s > endof(it)

  # Iterator for `event_types`
  immutable EventIterator
    handle::PerfFileHandle
  end

  endof(it::EventIterator) =
    div(it.handle.header.event_types.size,
      sizeof(perf_trace_event_type))
  length(it::EventIterator) = endof(it)

  function getindex(it::EventIterator, n)
    seek(it.handle.io,
      it.handle.header.event_types.offset +
      (n-1)*sizeof(perf_trace_event_type))
    unpack(it.handle.io, perf_trace_event_type)
  end

  start(it::EventIterator) = 1
  next(it::EventIterator, s) = (it[s], s+1)
  done(it::EventIterator, s) = s > endof(it)

  immutable Record
    event_type::UInt32
    misc::UInt16
    record::Any
  end
  Base.show(io::IO, rec::Record) =
    println(io, PERF_EVENT_TYPE[rec.event_type], ": ", rec.record)

  function read_record(io)
    header = unpack(io, perf_event_header)
    pos = position(io)
    function read_remaining_array()
      remaining_length = header.size -
       sizeof(perf_event_header) -
        (position(io) - pos)
      read(io, UInt8, remaining_length)
    end
    read_remaining_string() = String(read_remaining_array())
    rec = Record(header.event_type, header.misc,
      if header.event_type == PERF_RECORD_ITRACE_START
          unpack(io, itrace_start)
      elseif header.event_type == PERF_RECORD_AUX
          unpack(io, aux_record)
      elseif header.event_type == PERF_RECORD_COMM
          (unpack(io, itrace_start),
           read_remaining_string())
      elseif header.event_type == PERF_RECORD_MMAP
          x = unpack(io, itrace_start)
          y = unpack(io, mmap_record)
          z = read_remaining_string()
          idx = findfirst(z,'\0')
          z = idx == 0 ? z : z[1:idx-1]
          (x,y,z)
      elseif header.event_type == PERF_RECORD_MMAP2
          x = unpack(io, itrace_start)
          y = unpack(io, mmap_record)
          z = unpack(io, mmap2_extra_record)
          fname = read_remaining_string()
          idx = findfirst(fname,'\0')
          fname = idx == 0 ? fname : fname[1:idx-1]
          (x,y,z,fname)
      elseif header.event_type == PERF_RECORD_SAMPLE
          read_remaining_array()
      elseif header.event_type == PERF_RECORD_FORK
          (unpack(io, itrace_start),
           unpack(io, itrace_start),
           read_remaining_array())
      elseif header.event_type == PERF_RECORD_THROTTLE ||
        header.event_type == PERF_RECORD_UNTHROTTLE
          nothing
      elseif header.event_type == PERF_RECORD_EXIT
          (unpack(io, itrace_start),
           unpack(io, itrace_start),
           read(io, UInt64),
           read_remaining_array())
      elseif header.event_type == PERF_RECORD_TIME_CONV
          unpack(io, time_conv_record)
      elseif header.event_type == PERF_RECORD_FINISHED_ROUND
      else
          error(string("Unknown record type ",
            try PERF_EVENT_TYPE[header.event_type];
            catch; string("0x",hex(header.event_type)); end))
      end)
  end

  # Iterator for `data`
  immutable RecordIterator
    handle::PerfFileHandle
  end
  Base.iteratorsize(::Type{RecordIterator}) = Base.SizeUnknown()

  start(it::RecordIterator) = 0
  done(it::RecordIterator, s) = s >= it.handle.header.data.size
  function next(it::RecordIterator, s)
    seek(io, it.handle.header.data.offset + s)
    return (read_record(it.handle.io), s + header.size)
  end

  # Iterator for use in sorted_records
  immutable RecordChunkIterator
    handle::PerfFileHandle
    chunk::Vector{Tuple{UInt64, UInt64}}
  end
  Base.iteratorsize(::Type{RecordIterator}) = Base.HasLength()
  function Base.getindex(it::RecordChunkIterator, i)
    seek(it.handle.io, it.handle.header.data.offset + it.chunk[i][2])
    read_record(it.handle.io)
  end
  Base.start(it::RecordChunkIterator) = 1
  Base.done(it::RecordChunkIterator, i) = i > length(it.chunk)
  Base.next(it::RecordChunkIterator, i) = (it[i], i+1)

  function read_ids(h, ids::perf_file_section)
    seek(h.io, ids.offset)
    read(h.io, UInt64, div(ids.size,sizeof(UInt64)))
  end

  # Computes the offset to the timestamp from the start (for PERF_RECORD_SAMPLE)
  # or end (for everything else) of the record.
  function time_offsets_for_attr(attr)
    (sizeof(UInt64)*count_ones((attr.sample_type &
      (PERF_SAMPLE_IDENTIFIER | PERF_SAMPLE_IP | PERF_SAMPLE_TID))),
     sizeof(UInt64)*(1+count_ones((attr.sample_type &
      (PERF_SAMPLE_IDENTIFIER | PERF_SAMPLE_CPU | PERF_SAMPLE_STREAM_ID | PERF_SAMPLE_ID)))))
  end

  """
    Records in a perf.data file are not sorted in chronological order.
    This is because a single perf.data file may contain records from
    several CPUs and though the contents of each individual CPU buffer
    is ordered chronologicallly, the overall file is not.

    To be able to recover chronological order, perf_events has the ability
    to include a time stamp in all sample and non-sample events (the latter
    only if sample_id_all=1).

    One slight optimization is that perf marks the spot at which it is done
    going through all the CPU buffers and since we know that timestamps are
    monotonically increasing within one CPU buffer, we can use this to
    accelerate sorting.

    This function will sort the records by timestamp
  """
  function sorted_record_chunks(callback, h)
    header = h.header
    # Pair of the attribute and a vector of all corresponding ids
    attrs = Tuple{perf_event_attr,Vector{UInt64}}[]
    for file_attr in AttrIterator(h)
      push!(attrs,(file_attr.attr,read_ids(h, file_attr.ids)))
    end
    # If there are multiple attrs in the event stream, we'll need to disambiguate
    # them by id.
    id_to_attr_map = Dict{UInt64, perf_event_attr}()
    need_to_check_ids = length(attrs) > 1
    if need_to_check_ids
      for attr in attrs
        if (attr[1].sample_type & PERF_SAMPLE_IDENTIFIER) == 0
          error("sample_type needs to include PERF_SAMPLE_IDENTIFIER for sorting")
        end
        if !attr[1].sample_id_all
          error("sample_id_all needs to be set for sorting")
        end
        for id in attr[2]
          id_to_attr_map[id] = attr[1]
        end
      end
    end
    for attr in attrs
      if (attr[1].sample_type & PERF_SAMPLE_TIME) == 0
        error("sample_type needs to include PERF_SAMPLE_TIME for sorting")
      end
    end
    # Pair of a record's timestamp and its offset in the file. We keep two of
    # these buffers. At every FINISHED_ROUND barrier, we know that that no
    # future event will have a timestamp smaller than the largest we have
    # already seen. However, we do not know that no old events are straggling
    # and will be seen in the current round. Thus we keep two lists, with
    # timestamps split by the maximum timestamp of the last round.
    last_round_records = Tuple{UInt64,UInt64}[]
    records = Tuple{UInt64,UInt64}[]
    last_round_max_timestamp = 0
    current_round_max_timestamp = 0
    data_offset = 0
    data_size = header.data.size
    attr = attrs[1][1]
    while data_offset < data_size
      seek(h.io, header.data.offset + data_offset)
      record_header = unpack(h.io, perf_event_header)
      if record_header.event_type == PERF_RECORD_FINISHED_ROUND
        # Last round records is now complete. Sort them and off we go.
        callback(RecordChunkIterator(h, sort(last_round_records, by=x->x[1])))
        last_round_records = records
        records = Tuple{UInt64,UInt64}[]
        last_round_max_timestamp = current_round_max_timestamp
        data_offset += record_header.size
        continue
      end
      if need_to_check_ids
        # We verified above that we have PERF_SAMPLE_IDENTIFIER,
        # which is at a known location (at the beginning for PERF_RECORD_SAMPLE,
        # at the end for everything else)
        if record_header.event_type != PERF_RECORD_SAMPLE
          skip(h.io, record_header.size - sizeof(perf_event_header) - sizeof(UInt64))
          # Otherwise, we're already in the right spot
        end
        attr = id_to_attr_map[read(h.io, UInt64)]
      end
      sample_off, non_sample_off = time_offsets_for_attr(attr)
      @show (sample_off, non_sample_off)
      seek(h.io, header.data.offset + data_offset +
        ((record_header.event_type == PERF_RECORD_SAMPLE) ?
          sizeof(perf_event_header) + sample_off :
          record_header.size - non_sample_off))
      timestamp = read(h.io, UInt64)
      push!(timestamp <= last_round_max_timestamp ? last_round_records : records,
        (timestamp, data_offset))
      current_round_max_timestamp = max(timestamp, current_round_max_timestamp)
      data_offset += record_header.size
    end
    # We have all records, both the last and the current round should be ready to go
    callback(RecordChunkIterator(h, sort(last_round_records, by=x->x[1])))
    callback(RecordChunkIterator(h, sort(records, by=x->x[1])))
  end

  include("constants.jl")
  include("utils.jl")
  include("samples.jl")
  include("x86.jl")

end # module