text:
# if index format is changed, than increase the format version number
const INDEX_FORMAT_VER = "v2"
const DNA2BIT = Dict('A'=>0, 'T'=>1, 'C'=>2, 'G'=>3)
const KMER = 16
const REF_SAMPLE_STEP = 4

"""
the position of a k-mer sequence is the position of its first base, normalized to the direction of strand +
this is why a kmer sequence and its reverse compelement has a position offset of KMER-1
for example
contig: AAAATTTTCCCCGGG[G]TCATGATTACCAACCAATACCGTGGGATGG
position of the first 16-mer AAAATTTTCCCCGGGG is: 1
position of reverse compelement of the first 16-mer, CCCCGGGGAAAATTTT is: 1+16-1 = 16, the position bracketed
"""
immutable Coord
    # of which contig, usually a chr, a gene, an exon or any sequence
    contig::Int16
    # offset in this contig, the sign of pos is the strand
    pos::Int32

    Coord(contig, pos) = new(Int16(contig), Int32(pos))
    Coord(contig, pos, strand) = new(Int16(contig), Int32(strand * pos))
end

# for debugging
function display_coords(coords::Array{Coord, 1})
    i = 0
    for coord in coords
        if valid(coord)
            print("(", coord.pos>0?"+":"", coord.contig, ":", coord.pos, ")\t")
        elseif coord.contig == -1
            print("(duplicate)\t")
        elseif coord.contig == -2
            print("(unknown)\t")
        end
        i+=1
        if i%10 == 0
            print("\n")
        end
    end
    print("\n")
end

# generate a coord represents a kmer key collision
function dup_coord()
    return Coord(-1, 0)
end

# generate a coord represents a unknown key
function unknown_coord()
    return Coord(-2, 0)
end

# distance of two coords
function distance(c1::Coord, c2::Coord)
    if c1.contig != c2.contig || sign(c1.pos) != sign(c2.pos)
        return Inf
    else
        return abs(c1.pos) - abs(c2.pos)
    end
end

# abs distance of two coords
function abs_distance(c1::Coord, c2::Coord)
    return abs(distance(abs(c1), abs(c2)))
end

-(c1::Coord, c2::Coord) = distance(c1::Coord, c2::Coord)
abs(c::Coord) = Coord(c.contig, abs(c.pos))

is_dup(coord::Coord) = (coord.contig == -1)
is_unknown(coord::Coord) = (coord.contig == -2)
valid(coord::Coord) = (coord.contig >= 0)

typealias KmerCoord Dict{Int64, Coord}
typealias KmerCoordList Dict{Int64, Array{Coord, 1}}

function kmer2key(seq::Sequence)
    kmer2key(seq.seq)
end

# convert a kmer seq to a int64 key
function kmer2key(str::ASCIIString)
    str = uppercase(str)
    if contains(str, "N")
        return -1
    end
    key = 0
    try
        for c in str
            key = key*4 + DNA2BIT[c]
        end
    catch e
        info(e)
        error("Invalid k-mer $str")
    end
    return key
end

# add and index a kmer of a panel
function add_to_panel_index(kmer_coord::KmerCoord, seq::Sequence, coord::Coord)
    key = kmer2key(seq)
    if haskey(kmer_coord, key)
        kmer_coord[key] = dup_coord()
        return false
    end
    kmer_coord[key] = coord
    return true
end

# add and index all kmers of a contig
function index_contig(kmer_coord::KmerCoord, contig_seq::Sequence, contig_number::Int)
    len = length(contig_seq)
    for i in 1:len-KMER+1
        seq = contig_seq[i:i+KMER-1]
        add_to_panel_index(kmer_coord, seq, Coord(contig_number, i))
        add_to_panel_index(kmer_coord, ~seq, Coord(contig_number, - (i+KMER-1)))
    end
end

function load_bed(bed_file::AbstractString)
    io = open(bed_file)
    bed_file = readall(io)
    lines = split(bed_file, '\n')
    contig_number = 0
    panel = Dict{Int16, Dict{}}()
    chr_contigs = Dict{ASCIIString, Array{Int}}()
    for line in lines
        line = rstrip(line, '\n')
        cols = split(line)
        if length(cols)<4
            continue
        end
        contig_number += 1
        chr = ASCIIString(cols[1])
        if !haskey(chr_contigs, chr)
            chr_contigs[chr]=Array{Int,1}()
        end
        push!(chr_contigs[chr], contig_number)
        from = parse(Int64, ASCIIString(cols[2]))
        to = parse(Int64, ASCIIString(cols[3]))
        contig_name = ASCIIString(cols[4])
        panel[contig_number] = Dict("chr"=>chr, "name"=>contig_name, "from"=>from, "to"=>to)
    end
    return panel, chr_contigs
end

function get_ref_files(ref_path::AbstractString)
    ref_files = []
    # if ref_path is a folder, then find all fa files in it
    if isdir(ref_path)
        files = readdir(ref_path)
        for file in files
            if endswith(file, ".fa")
                push!(ref_files, joinpath(ref_path, file))
            end
        end
    else
        push!(ref_files, ref_path)
    end
    return ref_files
end

function index_chr_bed(chrfa::FastaRead, panel_kmer_coord::KmerCoord, chr_contigs, panel, panel_seq)
    contig_numbers = chr_contigs[chrfa.name]
    chr_seq = chrfa.sequence
    for contig_number in contig_numbers
        from = panel[contig_number]["from"]
        to = panel[contig_number]["to"]
        contig_seq = chr_seq[from:to]
        index_contig(panel_kmer_coord, contig_seq, contig_number)
        panel_seq[contig_number] = contig_seq
    end
end

# ref_path is a folder contains fasta files by chromosomes
# like chr1.fa, chr2.fa ...
function make_panel_index(ref_path::AbstractString, bed_file::AbstractString)
    panel_kmer_coord = KmerCoord()
    panel_seq = Dict{Int16, Sequence}()
    panel, chr_contigs = load_bed(bed_file)
    ref_files = get_ref_files(ref_path)

    ref = []
    # for each reference file, get its chromosome and index them
    for ref_file in ref_files
        io = fasta_open(ref_file)
        pos = position(io)
        while (chrfa = fasta_read(io))!=false
            if haskey(chr_contigs, chrfa.name)
                index_chr_bed(chrfa, panel_kmer_coord, chr_contigs, panel, panel_seq)
                info("index_chr_bed " * chrfa.name)
            end
            push!(ref, Dict("file"=>ref_file, "position"=>pos, "length"=>length(chrfa)))
            pos = position(io)
        end
    end

    ref_kmer_coords = make_kmer_coord_list(ref, panel_kmer_coord)
    return Dict("panel"=>panel, "seq"=>panel_seq, "kmer_coord"=>panel_kmer_coord, "ref_kmer_coords"=>ref_kmer_coords)
end

function make_kmer_shared_array(panel_kmer_coord::KmerCoord)
    kmer_buffer_size = div(4^KMER, 8)
    shared_kmer = SharedArray(UInt8, kmer_buffer_size, init=true)
    for k in keys(panel_kmer_coord)
        offset = div(k, 8) + 1
        shift = k % 8
        shared_kmer[offset] |= (0b00000001 << shift)
    end
    return shared_kmer
end

# make an index with the reference data and a panel
# for each kmer of the panel, create an array, and store the coordinations of same kmer
function make_kmer_coord_list(ref, panel_kmer_coord::KmerCoord)
    ref_index = KmerCoordList()
    total = 0
    for k in keys(panel_kmer_coord)
        ref_index[k]=Array{Coord, 1}()
    end

    shared_kmer = make_kmer_shared_array(panel_kmer_coord)

    # create parallel tasks
    tasks = []
    for chrid in 1:length(ref)
        chrinfo = ref[chrid]
        task = Dict("chrid"=>chrid, "chrinfo"=>chrinfo, "shared_kmer"=>shared_kmer)
        push!(tasks, task)
    end
    if length(workers()) <= 1
        #addprocs()
    end

    # run parallel for indexing
    info("tasks:" * string(length(tasks)))
    info("start worker processes")
    results = pmap(make_kmer_coord_list_chr, tasks)
    #results = []
    #for task in tasks
    #    push!(results, make_kmer_coord_list_chr(task))
    #end

    # merge the result index
    i = 0
    for result in results
        i += 1
        info("merging " * string(i))
        if isa(result, RemoteException)
            info(result)
            continue
        end
        for (k, v) in result
            append!(ref_index[k], v)
        end
        gc()
    end

    info("merge done")

    # destroy worker processes
    worker_procs = workers()
    if length(worker_procs)>1 || (length(worker_procs) == 1 && worker_procs[1]!=1)
        rmprocs(worker_procs)
    end

    return ref_index
end

function inkmer(shared_kmer::SharedArray{UInt8, 1}, k::Int64)
    offset = div(k, 8) + 1
    shift = k % 8
    return (shared_kmer[offset] & (0b00000001 << shift))>0
end

# run in parallel
# make an index with a chromosome of the reference data and a panel
function make_kmer_coord_list_chr(task)
    chrid=task["chrid"]
    chrinfo=task["chrinfo"]
    io = fasta_open(chrinfo["file"])
    chrname = basename(chrinfo["file"])
    seek(io, chrinfo["position"])
    len = chrinfo["length"]
    info("indexing $chrid:" * chrinfo["file"])
    fa = fasta_read(io)
    chrseq = fa.sequence
    shared_kmer=task["shared_kmer"]
    ref_index = KmerCoordList()
    total = 0
    for i in 1:REF_SAMPLE_STEP:len-KMER+1
        if i%10000000 < REF_SAMPLE_STEP
            info("$chrid-$chrname:$i/$len:$total")
        end
        seq = chrseq[i:i+KMER-1]
        key = kmer2key(seq)
        if inkmer(shared_kmer, key)
            if !haskey(ref_index, key)
                ref_index[key]=Array{Coord, 1}()
            end
            push!(ref_index[key], Coord(Int16(chrid), i))
            total+=1
        end
        key = kmer2key(~seq)
        if inkmer(shared_kmer, key)
            if !haskey(ref_index, key)
                ref_index[key]=Array{Coord, 1}()
            end
            push!(ref_index[key], Coord(Int16(chrid), -(i+KMER-1)))
            total+=1
        end
    end
    gc()
    info("finished $chrid:" * chrinfo["file"])
    return ref_index
end

function get_cache_path(ref_path::AbstractString, bed_file::AbstractString)
    ref_name = basename(ref_path) * "." * string(filesize(ref_path))
    bed_name = basename(bed_file) * "." * string(filesize(bed_file))
    cache_name = bed_name * "_" * ref_name * "_" *  INDEX_FORMAT_VER * ".idx"
    cache_path = joinpath(dirname(bed_file), cache_name)
    return cache_path
end

function index_bed(ref_path::AbstractString, bed_file::AbstractString)
    cache_path = get_cache_path(ref_path, bed_file)
    loaded = false
    # load the index from a cache file
    if isfile(cache_path)
        try
            io = open(cache_path)
            index = deserialize(io)
            if !isa(index, Dict)
                error("invalid index file, index should be a serialized Dict")
            elseif !haskey(index, "panel") || !haskey(index, "seq") || !haskey(index, "kmer_coord") || !haskey(index, "ref_kmer_coords")
                error("invalid index file, wrong Dict keys")
            end

            loaded = true
            return index
        catch(e)
            warn(e)
            warn("Failed to load the pre-built index ($cache_path), maybe it is not completed when saving it. Attemping to delete it now!")
            try
                rm(cache_path)
            catch(e)
                warn(e)
                warn("Failed to delete $cache_path, please do it manually")
                warn("FusionDirect exited")
                exit(-1)
            end
        end
    end

    if !loaded
        info("Index doesn't exist, indexing now, it may take several minutes to a few hours")
        info("After the index is created, loading it will be very fast")
        index = make_panel_index(ref_path, bed_file)
        # save the index to a cache file
        try
            io = open(cache_path, "w")
            serialize(io, index)
        catch(e)
            warn("Unable to save index to $cache_path. Please check if the folder is writable for current user, if the index is not cached, it will be built every time you run FusionDirect, which is slow.")
        end
        return index
    end
end


text:
"""
    Calculators

This module exists to bridge the gap between the `Models` and the `Dynamics`.

Here we provide functions and types for evaluating and storing quantities obtained from the
`Models`.
In addition any further manipulation of those quantities, such as computing eigenvalues,
is included here.

This module is largely needed to facilitate seemless integration of both ring polymer and
classical dynamics, using the same models and functions for both.
Specific ring polymer types are provided that have the extra fields and methods needed
to evaluate the quantities for each bead.
"""
module Calculators

using LinearAlgebra: LinearAlgebra, Hermitian, I, Eigen, tr
using StaticArrays: SMatrix, SVector

using NQCModels: NQCModels, Model, nstates
using NQCModels.AdiabaticModels: AdiabaticModel
using NQCModels.DiabaticModels: DiabaticModel, DiabaticFrictionModel
using NQCModels.FrictionModels: AdiabaticFrictionModel

using NQCDynamics: RingPolymers, ndofs

"""
    AbstractCalculator{M<:Model}

Top-level type for all calculators.

Each concrete calculator contains the `Model` and the fields to store the quantities
obtained from the model.
"""
abstract type AbstractCalculator{T,M<:Model} end
abstract type AbstractAdiabaticCalculator{T,M<:AdiabaticModel} <: AbstractCalculator{T,M} end
abstract type AbstractDiabaticCalculator{T,M<:Union{DiabaticFrictionModel,DiabaticModel}} <: AbstractCalculator{T,M} end
abstract type AbstractStaticDiabaticCalculator{T,M} <: AbstractDiabaticCalculator{T,M} end
abstract type AbstractFrictionCalculator{T,M<:AdiabaticFrictionModel} <: AbstractCalculator{T,M} end

Base.eltype(::AbstractCalculator{T}) where {T} = T

mutable struct AdiabaticCalculator{T,M} <: AbstractAdiabaticCalculator{T,M}
    model::M
    potential::T
    derivative::Matrix{T}
    function AdiabaticCalculator{T}(model::M, atoms::Integer) where {T,M<:Model}
        new{T,M}(model, 0, zeros(ndofs(model), atoms))
    end
end

struct RingPolymerAdiabaticCalculator{T,M} <: AbstractAdiabaticCalculator{T,M}
    model::M
    potential::Vector{T}
    derivative::Array{T,3}
    function RingPolymerAdiabaticCalculator{T}(model::M, atoms::Integer, beads::Integer) where {T,M<:Model}
        new{T,M}(model, zeros(beads), zeros(ndofs(model), atoms, beads))
    end
end

mutable struct DiabaticCalculator{T,M,S,L} <: AbstractDiabaticCalculator{T,M}
    model::M
    potential::Hermitian{T,SMatrix{S,S,T,L}}
    derivative::Matrix{Hermitian{T,SMatrix{S,S,T,L}}}
    eigen::LinearAlgebra.Eigen{T,T,SMatrix{S,S,T,L},SVector{S,T}}
    adiabatic_derivative::Matrix{SMatrix{S,S,T,L}}
    nonadiabatic_coupling::Matrix{SMatrix{S,S,T,L}}
    tmp_position::Matrix{T}
    tmp_mat::Matrix{T}
    tmp_mat_complex1::Matrix{Complex{T}}
    tmp_mat_complex2::Matrix{Complex{T}}
    function DiabaticCalculator{T}(model::M, atoms::Integer) where {T,M<:Model}
        n = nstates(model)
        matrix_template = NQCModels.DiabaticModels.matrix_template(model, T)
        vector_template = NQCModels.DiabaticModels.vector_template(model, T)

        potential = Hermitian(matrix_template)
        derivative = [Hermitian(matrix_template) for _=1:ndofs(model), _=1:atoms]
        eigen = Eigen(vector_template, matrix_template + I)
        adiabatic_derivative = [matrix_template for _ in CartesianIndices(derivative)]
        nonadiabatic_coupling = [matrix_template for _ in CartesianIndices(derivative)]
        tmp_position = zeros(T, ndofs(model), atoms)
        tmp_mat = zeros(T, n, n)
        tmp_mat_complex1 = zeros(Complex{T}, n, n)
        tmp_mat_complex2 = zeros(Complex{T}, n, n)

        new{T,M,n,n^2}(model,
            potential, derivative, eigen,
            adiabatic_derivative, nonadiabatic_coupling,
            tmp_position, tmp_mat, tmp_mat_complex1, tmp_mat_complex2)
    end
end

mutable struct RingPolymerDiabaticCalculator{T,M,S,L} <: AbstractDiabaticCalculator{T,M}
    model::M
    potential::Vector{Hermitian{T,SMatrix{S,S,T,L}}}
    traceless_potential::Vector{Hermitian{T,SMatrix{S,S,T,L}}}
    V̄::Vector{T}
    derivative::Array{Hermitian{T,SMatrix{S,S,T,L}},3}
    traceless_derivative::Array{Hermitian{T,SMatrix{S,S,T,L}},3}
    D̄::Array{T,3}
    eigen::Vector{LinearAlgebra.Eigen{T,T,SMatrix{S,S,T,L},SVector{S,T}}}
    adiabatic_derivative::Array{SMatrix{S,S,T,L},3}
    traceless_adiabatic_derivative::Array{SMatrix{S,S,T,L},3}
    nonadiabatic_coupling::Array{SMatrix{S,S,T,L},3}

    centroid_potential::Hermitian{T,SMatrix{S,S,T,L}}
    centroid_derivative::Matrix{Hermitian{T,SMatrix{S,S,T,L}}}
    centroid_eigen::LinearAlgebra.Eigen{T,T,SMatrix{S,S,T,L},SVector{S,T}}
    centroid_adiabatic_derivative::Matrix{SMatrix{S,S,T,L}}
    centroid_nonadiabatic_coupling::Matrix{SMatrix{S,S,T,L}}
    tmp_centroid::Matrix{T}

    tmp_mat::Matrix{T}
    tmp_mat_complex1::Matrix{Complex{T}}
    tmp_mat_complex2::Matrix{Complex{T}}
    function RingPolymerDiabaticCalculator{T}(model::M, atoms::Integer, beads::Integer) where {T,M<:Model}
        n = nstates(model)
        matrix_template = NQCModels.DiabaticModels.matrix_template(model, T)
        vector_template = NQCModels.DiabaticModels.vector_template(model, T)

        potential = [Hermitian(matrix_template) for _=1:beads]
        traceless_potential = [Hermitian(matrix_template) for _=1:beads]
        V̄ = zeros(beads)

        derivative = [Hermitian(matrix_template) for _=1:ndofs(model), _=1:atoms, _=1:beads]
        traceless_derivative = [Hermitian(matrix_template) for _=1:ndofs(model), _=1:atoms, _=1:beads]
        D̄ = zeros(ndofs(model), atoms, beads)

        adiabatic_derivative = [matrix_template for _=1:ndofs(model), _=1:atoms, _=1:beads]
        traceless_adiabatic_derivative = [Hermitian(matrix_template) for _=1:ndofs(model), _=1:atoms, _=1:beads]

        eigen = [Eigen(vector_template, matrix_template + I) for _=1:beads]
        nonadiabatic_coupling = [matrix_template for _=1:ndofs(model), _=1:atoms, _=1:beads]

        centroid_potential = Hermitian(matrix_template)
        centroid_derivative = [Hermitian(matrix_template) for _=1:ndofs(model), _=1:atoms]
        centroid_eigen = Eigen(vector_template, matrix_template + I)
        centroid_adiabatic_derivative = [matrix_template for _ in CartesianIndices(centroid_derivative)]
        centroid_nonadiabatic_coupling = [matrix_template for _ in CartesianIndices(centroid_derivative)]
        tmp_centroid = zeros(T, ndofs(model), atoms)

        tmp_mat = zeros(T, n, n)
        tmp_mat_complex1 = zeros(Complex{T}, n, n)
        tmp_mat_complex2 = zeros(Complex{T}, n, n)
        new{T,M,n,n^2}(model,
            potential, traceless_potential, V̄,
            derivative, traceless_derivative, D̄,
            eigen, adiabatic_derivative, traceless_adiabatic_derivative, nonadiabatic_coupling,
            centroid_potential, centroid_derivative, centroid_eigen,
            centroid_adiabatic_derivative, centroid_nonadiabatic_coupling, tmp_centroid,
            tmp_mat, tmp_mat_complex1, tmp_mat_complex2)
    end
end

function Calculator(model::DiabaticModel, atoms::Integer, t::Type{T}) where {T}
    DiabaticCalculator{t}(model, atoms)
end
function Calculator(model::AdiabaticModel, atoms::Integer, t::Type{T}) where {T}
    AdiabaticCalculator{t}(model, atoms)
end
function Calculator(model::DiabaticModel, atoms::Integer, beads::Integer, t::Type{T}) where {T}
    RingPolymerDiabaticCalculator{t}(model, atoms, beads)
end
function Calculator(model::AdiabaticModel, atoms::Integer, beads::Integer, t::Type{T}) where {T}
    RingPolymerAdiabaticCalculator{t}(model, atoms, beads)
end

function evaluate_potential!(calc::AbstractCalculator, R)
    calc.potential = NQCModels.potential(calc.model, R)
end

function evaluate_potential!(calc::AbstractCalculator, R::AbstractArray{T,3}) where {T}
    @views for i in axes(R, 3)
        calc.potential[i] = NQCModels.potential(calc.model, R[:,:,i])
    end
end

function evaluate_V̄!(calc::RingPolymerDiabaticCalculator)
    for i in 1:length(calc.V̄)
        calc.V̄[i] = tr(calc.potential[i]) / nstates(calc.model)
    end
end

function evaluate_traceless_potential!(calc::RingPolymerDiabaticCalculator)
    n = nstates(calc.model)
    for I in 1:length(calc.traceless_potential)
        calc.traceless_potential[I] = Hermitian(SMatrix{n,n}(
            i != j ? calc.potential[I][j,i] : calc.potential[I][j,i] - calc.V̄[I] for j=1:n, i=1:n
        ))
    end
end

function evaluate_centroid_potential!(calc::AbstractCalculator, R::AbstractMatrix)
    calc.centroid_potential = NQCModels.potential(calc.model, R)
end

function evaluate_derivative!(calc::AbstractCalculator, R)
    NQCModels.derivative!(calc.model, calc.derivative, R)
end

function evaluate_derivative!(calc::AbstractCalculator, R::AbstractArray{T,3}) where {T}
    @views for i in axes(R, 3)
        NQCModels.derivative!(calc.model, calc.derivative[:,:,i], R[:,:,i])
    end
end

function evaluate_D̄!(calc::RingPolymerDiabaticCalculator)
    for I in eachindex(calc.D̄)
        calc.D̄[I] = tr(calc.derivative[I]) / nstates(calc.model)
    end
end

function evaluate_traceless_derivative!(calc::RingPolymerDiabaticCalculator)
    n = nstates(calc.model)
    for I in eachindex(calc.derivative)
        calc.traceless_derivative[I] = Hermitian(SMatrix{n,n}(
            i != j ? calc.derivative[I][j,i] : calc.derivative[I][j,i] - calc.D̄[I] for j=1:n, i=1:n
        ))
    end
end

function evaluate_traceless_adiabatic_derivative!(calc::RingPolymerDiabaticCalculator)
    n = nstates(calc.model)
    for I in eachindex(calc.adiabatic_derivative)
        calc.traceless_adiabatic_derivative[I] = SMatrix{n,n}(
            i != j ? calc.adiabatic_derivative[I][j,i] : calc.adiabatic_derivative[I][j,i] - calc.D̄[I] for j=1:n, i=1:n
        )
    end
end

function evaluate_centroid_derivative!(calc::AbstractCalculator, R::AbstractMatrix)
    NQCModels.derivative!(calc.model, calc.centroid_derivative, R)
end

function eigen!(calc::DiabaticCalculator)
    eig = LinearAlgebra.eigen(calc.potential)
    corrected_vectors = correct_phase(eig.vectors, calc.eigen.vectors)
    calc.eigen = Eigen(eig.values, corrected_vectors)
    return nothing
end

function centroid_eigen!(calc::RingPolymerDiabaticCalculator)
    eig = LinearAlgebra.eigen(calc.centroid_potential)
    corrected_vectors = correct_phase(eig.vectors, calc.centroid_eigen.vectors)
    calc.centroid_eigen = Eigen(eig.values, corrected_vectors)
    return nothing
end

function eigen!(calc::RingPolymerDiabaticCalculator)
    for i=1:length(calc.potential)
        eig = LinearAlgebra.eigen(calc.potential[i])
        corrected_vectors = correct_phase(eig.vectors, calc.eigen[i].vectors)
        calc.eigen[i] = Eigen(eig.values, corrected_vectors)
    end
    return nothing
end

function correct_phase(new_vectors::SMatrix, old_vectors::SMatrix)
    n = size(new_vectors, 1)
    vect = SVector{n}(sign(LinearAlgebra.dot(new_vectors[:,i], old_vectors[:,i])) for i=1:n)
    return new_vectors .* vect'
end

function transform_derivative!(calc::AbstractDiabaticCalculator)
    for I in eachindex(calc.derivative)
        calc.adiabatic_derivative[I] = calc.eigen.vectors' * calc.derivative[I] * calc.eigen.vectors
    end
end

function transform_centroid_derivative!(calc::RingPolymerDiabaticCalculator)
    for I in eachindex(calc.centroid_derivative)
        calc.centroid_adiabatic_derivative[I] = calc.centroid_eigen.vectors' * calc.centroid_derivative[I] * calc.centroid_eigen.vectors
    end
end

function transform_derivative!(calc::RingPolymerDiabaticCalculator)
    for i in axes(calc.derivative, 3) # Beads
        for j in axes(calc.derivative, 2) # Atoms
            for k in axes(calc.derivative, 1) # DoFs
                calc.adiabatic_derivative[k,j,i] = calc.eigen[i].vectors' * calc.derivative[k,j,i] * calc.eigen[i].vectors
            end
        end
    end
end

function evaluate_nonadiabatic_coupling!(calc::AbstractDiabaticCalculator)
    for I in eachindex(calc.adiabatic_derivative)
        calc.nonadiabatic_coupling[I] = evaluate_nonadiabatic_coupling(calc.adiabatic_derivative[I], calc.eigen.values)
    end
end

function evaluate_centroid_nonadiabatic_coupling!(calc::RingPolymerDiabaticCalculator)
    for I in eachindex(calc.centroid_adiabatic_derivative)
        calc.centroid_nonadiabatic_coupling[I] = evaluate_nonadiabatic_coupling(calc.centroid_adiabatic_derivative[I], calc.centroid_eigen.values)
    end
end

function evaluate_nonadiabatic_coupling!(calc::RingPolymerDiabaticCalculator)
    for i in axes(calc.nonadiabatic_coupling, 3) # Beads
        for I in CartesianIndices(size(calc.adiabatic_derivative)[1:2])
            calc.nonadiabatic_coupling[I,i] = evaluate_nonadiabatic_coupling(calc.adiabatic_derivative[I,i], calc.eigen[i].values)
        end
    end
end

"""
# References

- HammesSchifferTully_JChemPhys_101_4657_1994 Eq. (32)
- SubotnikBellonzi_AnnuRevPhyschem_67_387_2016, section 2.3
"""
function evaluate_nonadiabatic_coupling(adiabatic_derivative::SMatrix, eigenvalues::SVector)
    n = length(eigenvalues)
    SMatrix{n,n}(
        (i != j ? adiabatic_derivative[j,i] / (eigenvalues[i] - eigenvalues[j]) : 0
        for j=1:n, i=1:n))
end

"""
Evaluates all electronic properties for the current position `r`.

# Properties evaluated:
- Diabatic potential
- Diabatic derivative
- Eigenvalues and eigenvectors
- Adiabatic derivative
- Nonadiabatic coupling
"""
function update_electronics!(calculator::AbstractDiabaticCalculator, r::AbstractArray)
    evaluate_potential!(calculator, r)
    evaluate_derivative!(calculator, r)
    eigen!(calculator)
    transform_derivative!(calculator)
    evaluate_nonadiabatic_coupling!(calculator)
end

function update_electronics!(calculator::RingPolymerDiabaticCalculator, r::AbstractArray{T,3}) where {T}
    evaluate_potential!(calculator, r)
    evaluate_derivative!(calculator, r)
    eigen!(calculator)
    transform_derivative!(calculator)
    evaluate_nonadiabatic_coupling!(calculator)

    update_centroid_electronics!(calculator, r)
end

function update_centroid_electronics!(calculator::RingPolymerDiabaticCalculator, r::AbstractArray{T,3}) where {T}
    RingPolymers.get_centroid!(calculator.tmp_centroid, r)
    evaluate_centroid_potential!(calculator, calculator.tmp_centroid)
    evaluate_centroid_derivative!(calculator, calculator.tmp_centroid)
    centroid_eigen!(calculator)
    transform_centroid_derivative!(calculator)
    evaluate_centroid_nonadiabatic_coupling!(calculator)
end

include("large_diabatic.jl")
include("friction.jl")

position(calc) = calc.tmp_position

for field ∈ (:potential, :derivative)
    get_field = Symbol(:get, :_, field)
    evaluate_field = Symbol(:evaluate, :_, field, :!)

    eval(quote
        function $get_field(calc, r)
            if position(calc) ≈ r
                return calc.$field
            else
                copyto!(position(calc), r)
                return $evaluate_field(calc, r)
            end
        end
    end)
end

end # module

ext: jl


text:
struct Pos
    ts::DateTime
    event_id::Int
    visible::Bool
    ll::LatLon{Float64}
    ecef::ECEF{Float64}
end

function Pos(row, datum = wgs84)
    ts = DateTime(row.timestamp[1:19], "yyyy-mm-dd HH:MM:SS")
    lat = row[Symbol("location-lat")]
    lon = row[Symbol("location-long")]
    ismissing(lat) | ismissing(lon) && return nothing
    ll = LatLon{Float64}(lat, lon)
    ecef = ECEF(LLA(ll), datum)

    return Pos(ts, row[Symbol("event-id")], row.visible, ll, ecef)
end

struct AnimalInfo
    taxon_canonical_name::String
    # tag_id::Int
    name_id::String
end

struct Track
    info::AnimalInfo
    route::Vector{Pos}
end

function load_animal_data(filename, datum = wgs84)
    data = CSV.File(filename)
    data = sort(data, by = x -> (x[Symbol("individual-local-identifier")], x.timestamp))
    tracks = Track[]
    taxon = ""
    name_id = ""
    # tag_id = -1
    route = Pos[]
    for row in data
        taxon2 = row[Symbol("individual-taxon-canonical-name")]
        name_id2 = row[Symbol("individual-local-identifier")]
        # tag_id2 = row[Symbol("tag-local-identifier")]
        # if (taxon != taxon2) | (name_id2 != name_id) | (tag_id2 != tag_id)
        if (taxon != taxon2) | (name_id2 != name_id)
            taxon, name_id = taxon2, name_id2
            info = AnimalInfo(taxon, name_id)
            route = Pos[]
            push!(tracks, Track(info, route))
        end
        pos = Pos(row, datum)
        if !isnothing(pos)
            push!(route, Pos(row, datum))
        end
    end

    return tracks
end


text:
module Weave
import Highlights
using Mustache
using Requires

function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" Base.include(Main, joinpath(@__DIR__, "plots.jl"))
    @require Gadfly="c91e804a-d5a3-530f-b6f0-dfbca275c004" Base.include(Main, joinpath(@__DIR__, "gadfly.jl"))
end

"""
    list_out_formats()

List supported output formats
"""
function list_out_formats()
  for format = keys(formats)
      println(string(format,": ",  formats[format].description))
  end
end

"""
    tangle(source::AbstractString; kwargs...)

Tangle source code from input document to .jl file.

## Keyword options

- `informat::Union{Symbol,AbstractString} = :auto`: Input document format. `:auto` will set it automatically based on file extension. You can also specify either of `"script"`, `"markdown"`, `"notebook"`, or `"noweb"`
- `out_path::Union{Symbol,AbstractString} = :doc`: Path where the output is generated can be either of:
  * `:doc`: Path of the source document (default)
  * `:pwd`: Julia working directory
  * `"somepath"`: `String` of output directory e.g. `"~/outdir"`, or of filename e.g. `"~/outdir/outfile.tex"`
"""
function tangle(
    source::AbstractString;
    out_path::Union{Symbol,AbstractString} = :doc,
    informat::Union{Symbol,AbstractString} = :auto
)
    doc = read_doc(source, informat)
    doc.cwd = get_cwd(doc, out_path)

    outname = get_outname(out_path, doc, ext = "jl")

    open(outname, "w") do io
    for chunk in doc.chunks
      if typeof(chunk) == CodeChunk
          options = merge(doc.chunk_defaults, chunk.options)
          if options[:tangle]
            write(io, chunk.content*"\n")
          end
      end
    end
  end
  doc.cwd == pwd()  && (outname = basename(outname))
  @info("Writing to file $outname")
end

"""
    weave(source::AbstractString; kwargs...)

Weave an input document to output file.

## Keyword options

- `doctype::Union{Symbol,AbstractString} = :auto`: Output document format. `:auto` will set it automatically based on file extension. You can also manually specify it; see [`list_out_formats()`](@ref) for the supported formats
- `informat::Union{Symbol,AbstractString} = :auto`: Input document format. `:auto` will set it automatically based on file extension. You can also specify either of `"script"`, `"markdown"`, `"notebook"`, or `"noweb"`
- `out_path::Union{Symbol,AbstractString} = :doc`: Path where the output is generated can be either of:
  * `:doc`: Path of the source document (default)
  * `:pwd`: Julia working directory
  * `"somepath"`: `String` of output directory e.g. `"~/outdir"`, or of filename e.g. `"~/outdir/outfile.tex"`
- `args::Dict = Dict()`: Arguments to be passed to the weaved document; will be available as `WEAVE_ARGS` in the document
- `mod::Union{Module,Symbol} = :sandbox`: Module where Weave `eval`s code. Defaults to `:sandbox` to create new sandbox module. You also can also pass a `Module` e.g. `Main`
- `fig_path::AbstractString = "figures"`: Where figures will be generated, relative to `out_path`
- `fig_ext::Union{Nothing,AbstractString} = nothing`: Extension for saved figures e.g. `".pdf"`, `".png"`. Default setting depends on `doctype`
- `cache_path::AbstractString = "cache"`: Where of cached output will be saved
- `cache::Symbol = :off`: Controls caching of code:
  * `:off` means no caching (default)
  * `:all` caches everything
  * `:user` caches based on chunk options
  * `:refresh` runs all code chunks and save new cache
- `throw_errors::Bool = false`: If `false` errors are included in output document and the whole document is executed. If `true` errors are thrown when they occur
- `template::Union{Nothing,AbstractString,Mustache.MustacheTokens} = nothing`: Template (file path) or `Mustache.MustacheTokens`s for `md2html` or `md2tex` formats
- `highlight_theme::Union{Nothing,Type{<:Highlights.AbstractTheme}} = nothing`: Theme used for syntax highlighting (defaults to `Highlights.Themes.DefaultTheme`)
- `css::Union{Nothing,AbstractString} = nothing`: Path of a CSS file used for md2html format
- `pandoc_options::Vector{<:AbstractString} = String[]`: `String`s of options to pass to pandoc for `pandoc2html` and `pandoc2pdf` formats, e.g. `["--toc", "-N"]`
- `latex_cmd::AbstractString = "xelatex"`: The command used to make PDF file from .tex
- `latex_keep_unicode::Bool = false`: If `true`, do not convert unicode characters to their respective latex representation. This is especially useful if a font and tex-engine with support for unicode characters are used

!!! note
    Run Weave from terminal and try to avoid weaving from IJulia or ESS; they tend to mess with capturing output.
"""
function weave(
    source::AbstractString;
    doctype::Union{Symbol,AbstractString} = :auto,
    informat::Union{Symbol,AbstractString} = :auto,
    out_path::Union{Symbol,AbstractString} = :doc,
    args::Dict = Dict(),
    mod::Union{Module,Symbol} = :sandbox,
    fig_path::AbstractString = "figures",
    fig_ext::Union{Nothing,AbstractString} = nothing,
    cache_path::AbstractString = "cache",
    cache::Symbol = :off,
    throw_errors::Bool = false,
    template::Union{Nothing,AbstractString,Mustache.MustacheTokens} = nothing,
    highlight_theme::Union{Nothing,Type{<:Highlights.AbstractTheme}} = nothing,
    css::Union{Nothing,AbstractString} = nothing,
    pandoc_options::Vector{<:AbstractString} = String[],
    latex_cmd::AbstractString = "xelatex",
    latex_keep_unicode::Bool = false
)
    doc = read_doc(source, informat)
    doctype == :auto && (doctype = detect_doctype(doc.source))
    doc.doctype = doctype

    # Read args from document header, overrides command line args
    if haskey(doc.header, "options")
        (doctype, informat, out_path, args, mod, fig_path, fig_ext,
        cache_path, cache, throw_errors, template, highlight_theme, css,
        pandoc_options, latex_cmd) = header_args(doc, out_path, mod,
                                    fig_ext, fig_path,
                                    cache_path, cache, throw_errors,
                                    template, highlight_theme, css,
                                    pandoc_options, latex_cmd)
    end

    template != nothing && (doc.template = template)
    highlight_theme != nothing && (doc.highlight_theme = highlight_theme)
    #theme != nothing && (doc.theme = theme) #Reserved for themes
    css != nothing && (doc.css = css)

    doc = run(doc, doctype = doctype,
                mod = mod,
                out_path=out_path, args = args,
                fig_path = fig_path, fig_ext = fig_ext, cache_path = cache_path, cache=cache,
                throw_errors = throw_errors,latex_keep_unicode=latex_keep_unicode)
    formatted = format(doc)

    outname = get_outname(out_path, doc)

    open(outname, "w") do io
        write(io, formatted)
    end

    #Special for that need external programs
    if doc.doctype == "pandoc2html"
        mdname = outname
        outname = get_outname(out_path, doc, ext = "html")
        pandoc2html(formatted, doc, outname, pandoc_options)
        rm(mdname)
    elseif doc.doctype == "pandoc2pdf"
        mdname = outname
        outname = get_outname(out_path, doc, ext = "pdf")
        pandoc2pdf(formatted, doc, outname, pandoc_options)
        rm(mdname)
    elseif doc.doctype == "md2pdf"
        success = run_latex(doc, outname, latex_cmd)
        success || return
        outname = get_outname(out_path, doc, ext = "pdf")
    end

    doc.cwd == pwd() && (outname = basename(outname))
    @info("Report weaved to $outname")
    return abspath(outname)
end

weave(doc::AbstractString, doctype::Union{Symbol,AbstractString}) =
    weave(doc; doctype = doctype)

"""
    notebook(source::AbstractString; kwargs...)

Convert Weave document `source` to Jupyter Notebook and execute the code
using [`nbconvert`](https://nbconvert.readthedocs.io/en/latest/).
**Ignores** all chunk options.

## Keyword options

- `out_path::Union{Symbol,AbstractString} = :pwd`: Path where the output is generated can be either of:
  * `:doc`: Path of the source document
  * `:pwd`: Julia working directory (default)
  * `"somepath"`: `String` of output directory e.g. `"~/outdir"`, or of filename e.g. `"~/outdir/outfile.tex"`
- `timeout = -1`: nbconvert cell timeout in seconds. Defaults to `-1` (no timeout)
- `nbconvert_options::AbstractString = ""`: `String` of additional options to pass to nbconvert, such as `"--allow-errors"`
- `jupyter_path::AbstractString = "jupyter"`: Path/command for the Jupyter you want to use. Defaults to `"jupyter"`, which runs whatever is linked/alias to that

!!! warning
    The code is _**not**_ executed by Weave, but by [`nbconvert`](https://nbconvert.readthedocs.io/en/latest/).
    This means that the output doesn't necessarily always work properly; see [#116](https://github.com/mpastell/Weave.jl/issues/116).

!!! note
    In order to _just_ convert Weave document to Jupyter Notebook,
    use [`convert_doc`](@ref) instead.
"""
function notebook(
    source::AbstractString;
    out_path::Union{Symbol,AbstractString} = :pwd,
    timeout = -1,
    nbconvert_options::AbstractString = "",
    jupyter_path::AbstractString = "jupyter",
)
    doc = read_doc(source)
    converted = convert_doc(doc, NotebookOutput())
    doc.cwd = get_cwd(doc, out_path)
    outfile = get_outname(out_path, doc, ext = "ipynb")

    open(outfile, "w") do f
        write(f, converted)
    end

    @info "Running nbconvert"
    out = read(
        `$jupyter_path nbconvert --ExecutePreprocessor.timeout=$timeout --to notebook --execute $outfile  $nbconvert_options --output $outfile`,
        String,
    )
end

"""
    include_weave(source::AbstractString, informat::Union{Symbol,AbstractString} = :auto)
    include_weave(m::Module, source::AbstractString, informat::Union{Symbol,AbstractString} = :auto)

Include code from Weave document calling `include_string` on all code from doc.
Code is run in the path of the include document.
"""
function include_weave(
    m::Module,
    source::AbstractString,
    informat::Union{Symbol,AbstractString} = :auto,
)
    old_path = pwd()
    doc = read_doc(source, informat)
    cd(doc.path)
    try
        code = join([x.content for x in filter(x -> isa(x, Weave.CodeChunk), doc.chunks)], "\n")
        include_string(m, code)
    catch e
        throw(e)
    finally
        cd(old_path)
    end
end

include_weave(source, informat=:auto) = include_weave(Main, source, informat)

#Hooks to run before and after chunks, this is form IJulia,
#but note that Weave hooks take the chunk as input
const preexecute_hooks = Function[]
push_preexecute_hook(f::Function) = push!(preexecute_hooks, f)
pop_preexecute_hook(f::Function) = splice!(preexecute_hooks, findfirst(x -> x == f, preexecute_hooks))

const postexecute_hooks = Function[]
push_postexecute_hook(f::Function) = push!(postexecute_hooks, f)
pop_postexecute_hook(f::Function) = splice!(postexecute_hooks, findfirst(x -> x == f, postexecute_hooks))

include("chunks.jl")
include("config.jl")
include("WeaveMarkdown/markdown.jl")
include("display_methods.jl")
include("readers.jl")
include("run.jl")
include("cache.jl")
include("formatters.jl")
include("format.jl")
include("pandoc.jl")
include("writers.jl")

export weave, list_out_formats, tangle, convert_doc, notebook,
        set_chunk_defaults, get_chunk_defaults, restore_chunk_defaults,
        include_weave
end


text:
module LeastAction

using FindSteadyStates
using SteadyStateDiffEq
using LinearAlgebra
using ModelingToolkit
using ForwardDiff
using Optim
using NLsolve
using Parameters

include("utils.jl")
include("fixedpoint_gen.jl")
include("DiffusionMatrix.jl")
include("gradient_gen.jl")
include("action.jl")

export fixedpoint_gen
export gradient_gen
export action

end


text:
# Aggregation on the five region level

const ssp5 = Dict(
	# R5.2OECD = Includes the OECD 90 and EU member states and candidates.
	"OECD" => [
		"Albania", "Australia", "Austria", "Belgium", "Bosnia and Herzegovina", "Bulgaria", "Canada",
		"Croatia", "Cyprus", "Czech Republic", "Denmark", "Estonia", "Finland", "France", "Germany",
		"Greece", "Guam", "Hungary", "Iceland", "Ireland", "Italy", "Japan", "Latvia", "Lithuania",
		"Luxembourg", "Malta", "Montenegro", "Netherlands", "New Zealand", "Norway", "Poland",
		"Portugal", "Puerto Rico", "Romania", "Serbia", "Slovakia", "Slovenia", "Spain", "Sweden",
		"Switzerland", "Macedonia", "Turkey", "United Kingdom", "United States"
	],

	# R5.2REF = Countries from the Reforming Economies of Eastern Europe and the Former Soviet Union.
	"REF" => [
		"Armenia", "Azerbaijan", "Belarus", "Georgia", "Kazakhstan", "Kyrgyzstan", "Moldova",
		"Russia", "Tajikistan", "Turkmenistan", "Ukraine", "Uzbekistan"
	],

	# R5.2ASIA = The region includes most Asian countries with the exception of the Middle East, Japan
	# and Former Soviet Union states. China includes Hong Kong and Macao, excludes Taiwan.
	"ASIA" => [
		"Afghanistan", "Bangladesh", "Bhutan", "Brunei", "Cambodia", "China", "North Korea", "Fiji",
		"French Polynesia", "India", "Indonesia", "Laos", "Malaysia", "Maldives", "Micronesia",
		"Mongolia", "Myanmar", "Nepal", "New Caledonia", "Pakistan", "Papua New Guinea", "Philippines",
		"South Korea", "Samoa", "Singapore", "Solomon Islands", "Sri Lanka", "Taiwan", "Thailand",
		"Timor-Leste", "Vanuatu", "Vietnam"
	],

	# R5.2MAF = This region includes the countries of the Middle East and Africa.
	"MAF" => [
		"Algeria", "Angola", "Bahrain", "Benin", "Botswana", "Burkina Faso", "Burundi", "Cameroon",
		"Cape Verde", "Central African Republic", "Chad", "Comoros", "Republic of Congo", "Côte d'Ivoire",
		"Democratic Republic of the Congo", "Djibouti", "Egypt", "Equatorial Guinea", "Eritrea", "Ethiopia",
		"Gabon", "Gambia", "Ghana", "Guinea", "Guinea-Bissau", "Iran", "Iraq", "Israel", "Jordan", "Kenya",
		"Kuwait", "Lebanon", "Lesotho", "Liberia", "Libya", "Madagascar", "Malawi", "Mali", "Mauritania",
		"Mauritius", "Mayotte", "Morocco", "Mozambique", "Namibia", "Niger", "Nigeria", "Palestina", "Oman",
		"Qatar", "Rwanda", "Reunion", "Saudi Arabia", "Senegal", "Sierra Leone", "Somalia", "South Africa",
		"South Sudan", "Sudan", "Swaziland", "Syria", "Togo", "Tunisia", "Uganda", "United Arab Emirates",
		"Tanzania", "Western Sahara", "Yemen", "Zambia", "Zimbabwe"
	],

	# R5.2LAM = This region includes the countries of Latin America and the Caribbean.
	"LAM" => [
		"Argentina", "Aruba", "Bahamas", "Barbados", "Belize", "Bolivia", "Brazil", "Chile", "Colombia",
		"Costa Rica", "Cuba", "Dominican Republic", "Ecuador", "El Salvador", "French Guiana", "Grenada",
		"Guadeloupe", "Guatemala", "Guyana", "Haiti", "Honduras", "Jamaica", "Martinique", "Mexico",
		"Nicaragua", "Panama", "Paraguay", "Peru", "Suriname", "Trinidad and Tobago",
		"Virgin Islands, U.S.", "Uruguay", "Venezuela"
	]
)

# Aggregation on the 32 region level

const ssp32 = Dict(
	"ANUZ" => ["Australia", "New Zealand"],
	"BRA" => ["Brazil"],
	"CAN" => ["Canada"],
	# Central Asia
	"CAS" => ["Armenia", "Azerbaijan", "Georgia", "Kazakhstan", "Kyrgyzstan", "Tajikistan",
				"Turkmenistan", "Uzbekistan"],
	# China (Mainland, Hongkong, Macao; excl. Taiwan)
	"CHN" => ["China", "Hong Kong", "Macao"],
	# Eastern Europe (excl. former Soviet Union and EU member states)
	"EEU" => ["Albania", "Bosnia and Herzegovina", "Croatia", "Montenegro", "Serbia", "Macedonia"],
	# Eastern Europe, former Soviet Union (excl. Russia and EU members)
	"EEU-FSU" => ["Belarus", "Moldova", "Ukraine"],
	"EFTA" => ["Iceland", "Norway", "Switzerland"],
	# New EU member states that joined as of 2004 - high income.
	"EU12-H" => ["Cyprus", "Czech Republic", "Estonia", "Hungary", "Malta", "Poland", "Slovakia", "Slovenia"],
	# New EU member states that joined as of 2004 - medium income.
	"EU12-M" => ["Bulgaria", "Latvia", "Lithuania", "Romania"],
	# This region includes European Union member states that joined prior to 2004.
	"EU15" => ["Austria", "Belgium", "Denmark", "Finland", "France", "Germany", "Greece", "Ireland", "Italy",
				"Luxembourg", "Netherlands", "Portugal", "Spain", "Sweden", "United Kingdom"],
	"IDN" => ["Indonesia"],
	"IND" => ["India"],
	"JPN" => ["Japan"],
	"KOR" => ["South Korea"],
	# This region includes the countries of Latin America (excl. Brazil, Mexico) - low income.
	"LAM-L" => ["Belize", "Guatemala", "Haiti", "Honduras", "Nicaragua"],
	# This region includes the countries of Latin America (excl. Brazil, Mexico) - medium and high income.
	# [Netherlands Antilles not in GADM (new status), replaced with islands listed after Venezuela]
	"LAM-M" => ["Antigua and Barbuda", "Argentina", "Bahamas", "Barbados", "Bermuda", "Bolivia", "Chile",
				"Colombia", "Costa Rica", "Cuba", "Dominica", "Dominican Republic", "Ecuador", "El Salvador",
				"French Guiana", "Grenada", "Guadeloupe", "Guyana", "Jamaica", "Martinique",
				"Panama", "Paraguay", "Peru", "Saint Kitts and Nevis", "Saint Lucia",
				"Saint Vincent and the Grenadines", "Suriname", "Trinidad and Tobago", "Uruguay", "Venezuela",
				"Aruba", "Bonaire, Sint Eustatius and Saba", "Sint Maarten", "Curaçao"],
	# This region includes the countries of Middle East Asia - high income.
	"MEA-H" => ["Bahrain", "Israel", "Kuwait", "Oman", "Qatar", "Saudi Arabia", "United Arab Emirates"],
	# This region includes the countries of Middle East Asia - low and medium income.
	"MEA-M" => ["Iran", "Iraq", "Jordan", "Lebanon", "Palestina", "Syria", "Yemen"],
	"MEX" => ["Mexico"],
	# This region includes the countries of North Africa.
	"NAF" => ["Algeria", "Egypt", "Libya", "Morocco", "Tunisia", "Western Sahara"],
	# This region includes the countries of Other Asia - former Centrally Planned Asia.
	"OAS-CPA" => ["Cambodia", "Laos", "Mongolia", "Vietnam"],
	# This region includes the countries of Other Asia - low income.
	"OAS-L" => ["Bangladesh", "North Korea", "Fiji", "Micronesia", "Myanmar", "Nepal", "Papua New Guinea",
				"Philippines", "Samoa", "Solomon Islands", "Timor-Leste", "Tonga", "Vanuatu"],
	# This region includes the countries of Other Asia - medium and high income.
	"OAS-M" => ["Bhutan", "Brunei", "French Polynesia", "Guam", "Malaysia", "Maldives", "New Caledonia",
				"Singapore", "Sri Lanka", "Thailand"],
	"PAK" => ["Pakistan", "Afghanistan"],
	"RUS" => ["Russia"],
	"SAF" => ["South Africa"],
	# This region includes the countries of Subsahara Africa (excl. South Africa) - low income.
	"SSA-L" => ["Benin", "Burkina Faso", "Burundi", "Cameroon", "Cape Verde", "Central African Republic",
				"Chad", "Comoros", "Republic of Congo", "Côte d'Ivoire", "Democratic Republic of the Congo",
				"Djibouti", "Eritrea", "Ethiopia", "Gambia", "Ghana", "Guinea", "Guinea-Bissau", "Kenya",
				"Lesotho", "Liberia", "Madagascar", "Malawi", "Mali", "Mauritania", "Mozambique", "Niger",
				"Nigeria", "Rwanda", "São Tomé and Príncipe", "Senegal", "Sierra Leone", "Somalia",
				"South Sudan", "Sudan", "Swaziland", "Togo", "Uganda", "Tanzania", "Zambia", "Zimbabwe"],
	# This region includes the countries of Subsahara Africa (excl. South Africa) - medium and high income.
	"SSA-M" => ["Angola", "Botswana", "Equatorial Guinea", "Gabon", "Mauritius", "Mayotte", "Namibia",
				"Reunion", "Seychelles"],
	"TUR" => ["Turkey"],
	"TWN" => ["Taiwan"],
	"USA" => ["United States", "Puerto Rico", "Virgin Islands, U.S."]
)

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
baremodule UnaryOps
    using ..Types
end

UnaryUnion = Union{AbstractUnaryOp, libgb.GrB_UnaryOp}

function unarynames(name)
    simple = splitconstant(name)[2]
    containername = Symbol(simple, "_UNARY_T")
    exportedname = Symbol(simple)
    return containername, exportedname
end

function createunaryops()
    builtins = ["GrB_IDENTITY", "GrB_AINV", "GxB_LNOT", "GrB_MINV", "GxB_ONE", "GrB_ABS", "GrB_BNOT", "GxB_SQRT", "GxB_LOG", "GxB_EXP", "GxB_LOG2", "GxB_SIN", "GxB_COS", "GxB_TAN", "GxB_ACOS", "GxB_ASIN", "GxB_ATAN", "GxB_SINH", "GxB_COSH", "GxB_TANH", "GxB_ASINH", "GxB_ACOSH", "GxB_ATANH", "GxB_SIGNUM", "GxB_CEIL", "GxB_FLOOR", "GxB_ROUND", "GxB_TRUNC", "GxB_EXP2", "GxB_EXPM1", "GxB_LOG10", "GxB_LOG1P", "GxB_LGAMMA", "GxB_TGAMMA", "GxB_ERF", "GxB_ERFC", "GxB_FREXPE", "GxB_FREXPX", "GxB_CONJ", "GxB_CREAL", "GxB_CIMAG", "GxB_CARG", "GxB_ISINF", "GxB_ISNAN", "GxB_ISFINITE","GxB_POSITIONI", "GxB_POSITIONI1", "GxB_POSITIONJ", "GxB_POSITIONJ1"]
    for name ∈ builtins
        containername, exportedname = unarynames(name)
        structquote = quote
            struct $containername <: AbstractUnaryOp
                pointers::Dict{DataType, libgb.GrB_UnaryOp}
                name::String
                $containername() = new(Dict{DataType, libgb.GrB_UnaryOp}(), $name)
            end
        end
        @eval(Types, $structquote)
        constquote = quote
            const $exportedname = Types.$containername()
            export $exportedname
        end
        @eval(UnaryOps, $constquote)
    end
end

function load(unaryop::AbstractUnaryOp)
    booleans = ["GrB_IDENTITY", "GrB_AINV", "GrB_MINV", "GxB_LNOT", "GxB_ONE", "GrB_ABS"]
    integers = ["GrB_IDENTITY", "GrB_AINV", "GrB_MINV", "GxB_LNOT", "GxB_ONE", "GrB_ABS", "GrB_BNOT"]
    unsignedintegers = ["GrB_IDENTITY", "GrB_AINV", "GrB_MINV", "GxB_LNOT", "GxB_ONE", "GrB_ABS", "GrB_BNOT"]
    floats = ["GrB_IDENTITY", "GrB_AINV", "GrB_MINV", "GxB_LNOT", "GxB_ONE", "GrB_ABS", "GxB_SQRT", "GxB_LOG", "GxB_EXP", "GxB_LOG2", "GxB_SIN", "GxB_COS", "GxB_TAN", "GxB_ACOS", "GxB_ASIN", "GxB_ATAN", "GxB_SINH", "GxB_COSH", "GxB_TANH", "GxB_ASINH", "GxB_ACOSH", "GxB_ATANH", "GxB_SIGNUM", "GxB_CEIL", "GxB_FLOOR", "GxB_ROUND", "GxB_TRUNC", "GxB_EXP2", "GxB_EXPM1", "GxB_LOG10", "GxB_LOG1P", "GxB_LGAMMA", "GxB_TGAMMA", "GxB_ERF", "GxB_ERFC", "GxB_FREXPE", "GxB_FREXPX", "GxB_ISINF", "GxB_ISNAN", "GxB_ISFINITE"]
    positionals = ["GxB_POSITIONI", "GxB_POSITIONI1", "GxB_POSITIONJ", "GxB_POSITIONJ1"]
    name = unaryop.name
    if name ∈ booleans
        constname = name * "_BOOL"
        unaryop.pointers[Bool] = load_global(constname)
    end

    if name ∈ integers
        unaryop.pointers[Int8] = load_global(name * "_INT8")
        unaryop.pointers[Int16] = load_global(name * "_INT16")
        unaryop.pointers[Int32] = load_global(name * "_INT32")
        unaryop.pointers[Int64] = load_global(name * "_INT64")
    end

    if name ∈ unsignedintegers
        unaryop.pointers[UInt8] = load_global(name * "_UINT8")
        unaryop.pointers[UInt16] = load_global(name * "_UINT16")
        unaryop.pointers[UInt32] = load_global(name * "_UINT32")
        unaryop.pointers[UInt64] = load_global(name * "_UINT64")
    end

    if name ∈ floats
        unaryop.pointers[Float32] = load_global(name * "_FP32")
        unaryop.pointers[Float64] = load_global(name * "_FP64")
    end
end

ztype(op::UnaryUnion) = ptrtogbtype(libgb.GxB_UnaryOp_ztype(op))
xtype(op::UnaryUnion) = ptrtogbtype(libgb.GxB_UnaryOp_xtype(op))


text:
# D2Q17
# See lectures in LBM

# Picture is D2Q9, but with each vector extended to a 3 x 3 box

# Contains description:
# Multiple-relaxation time model for the correct thermodynamic equations
#
# https://www.researchgate.net/publication/23315587_Multiple-relaxation_time_model_for_the_correct_thermodynamic_equations
#
# https://www.researchgate.net/profile/Zhaoli_Guo2/publication/23315587_Multiple-relaxation_time_model_for_the_correct_thermodynamic_equations/links/0046352b40f45699f8000000.pdf
struct D2Q17{
    T <: Real,
    Abscissaes <: AbstractMatrix{<:Integer},
    Weights <: AbstractVector{<:Real},
} <: Quadrature
    abscissae::Abscissaes
    weights::Weights
    speed_of_sound_squared::T
end
function D2Q17()
    sq = sqrt(193)
    w_0 = (575 + 193sq) / 8100
    w_1 = (3355 - 91sq) / 18000
    w_2 = (655 + 17sq) / 27000
    w_3 = (685 - 49sq) / 54000
    w_4 = (1445 - 101sq) / 162000

    return D2Q17(
        [   #A #B          #C          #D          #D
            0 1 -1 0 0 1 -1 1 -1 2 -2 2 -2 3 -3 0 0
            0 0 0 1 -1 1 -1 -1 1 2 -2 -2 2 0 0 3 -3
        ],
        [
            w_0,
            w_1,
            w_1,
            w_1,
            w_1,
            w_2,
            w_2,
            w_2,
            w_2,
            w_3,
            w_3,
            w_3,
            w_3,
            w_4,
            w_4,
            w_4,
            w_4,
        ],
        (125 + 5 * sqrt(193)) / 72,
    )
end

order(q::D2Q17) = 7
Base.show(io::IO, q::D2Q17) = show(io, "D2Q17")
Base.string(q::D2Q17) = "D2Q17"


text:
# This file includes RODESystems

import DifferentialEquations: RandomEM, RODEProblem
import UUIDs: uuid4

"""
    @def_rode_system ex

where `ex` is the expression to define to define a new AbstractRODESystem component type. The usage is as follows:
```julia
@def_rode_system mutable struct MyRODESystem{T1,T2,T3,...,TN,OP,RH,RO,ST,IP,OP} <: AbstractRODESystem
    param1::T1 = param1_default                 # optional field
    param2::T2 = param2_default                 # optional field
    param3::T3 = param3_default                 # optional field
        ⋮
    paramN::TN = paramN_default                 # optional field
    righthandside::RH = righthandside_function  # mandatory field
    readout::RO = readout_function              # mandatory field
    state::ST = state_default                   # mandatory field
    input::IP = input_default                   # mandatory field
    output::OP = output_default                 # mandatory field
end
```
Here, `MyRODESystem` has `N` parameters. `MyRODESystem` is represented by the `righthandside` and `readout` function. `state`, `input` and `output` is the initial state, input port and output port of `MyRODESystem`.

!!! warning
    `righthandside` must have the signature
    ```julia
    function righthandside((dx, x, u, t, W, args...; kwargs...)
        dx .= .... # update dx
    end
    ```
    and `readout` must have the signature
    ```julia
    function readout(x, u, t)
        y = ...
        return y
    end
    ```

!!! warning
    New RODE system must be a subtype of `AbstractRODESystem` to function properly.

# Example
```julia
julia> @def_rode_system mutable struct MySystem{RH, RO, IP, OP} <: AbstractRODESystem
           A::Matrix{Float64} = [2. 0.; 0 -2]
           righthandside::RH = (dx, x, u, t, W) -> (dx .= A * x * W)
           readout::RO = (x, u, t) -> x
           state::Vector{Float64} = rand(2)
           input::IP = nothing
           output::OP = Outport(2)
       end

julia> ds = MySystem();
```
"""
macro def_rode_system(ex)
    checksyntax(ex, :AbstractRODESystem)
    appendcommonex!(ex)
    foreach(nex -> appendex!(ex, nex), [
        :( alg::$ALG_TYPE_SYMBOL = Causal.RandomEM() ),
        :( integrator::$INTEGRATOR_TYPE_SYMBOL = Causal.construct_integrator(
            Causal.RODEProblem, input, righthandside, state, t, modelargs, solverargs; alg=alg, modelkwargs=modelkwargs,
            solverkwargs = :dt in keys(solverkwargs) ? solverkwargs :
                (; zip((keys(solverkwargs)..., :dt), (values(solverkwargs)..., 0.01))...),
            numtaps=3) )
        ])
    quote
        Base.@kwdef $ex
    end |> esc
end

##### Define RODE sytem library

"""
    $TYPEDEF

Constructs a generic RODE system

# Fields

    $TYPEDFIELDS
"""
@def_rode_system mutable struct RODESystem{RH,
                                           RO,
                                           ST <: AbstractVector{<:Real},
                                           IP <: Union{<:Inport, <:Nothing},
                                           OP <: Union{<:Outport,<:Nothing}} <: AbstractRODESystem
    "Right-hand-side function"
    righthandside::RH
    "Readout function"
    readout::RO
    "State"
    state::ST
    "Input. Expected to be an `Inport` or `Nothing`"
    input::IP
    "Output port"
    output::OP
end

"""
    $TYPEDEF

Constructs a `MultiplicativeNoiseLinearSystem` with the dynamics
```math
\\begin{array}{l}
    \\dot{x} = A x W
\\end{array}
where `W` is the noise process.
```

# Fields

    $TYPEDFIELDS
"""
@def_rode_system mutable struct MultiplicativeNoiseLinearSystem{T1 <: AbstractMatrix{<:Real},
                                                                RH,
                                                                RO,
                                                                ST <: AbstractVector{<:Real},
                                                                IP <: Union{<:Inport, <:Nothing},
                                                                OP <: Union{<:Outport,<:Nothing}} <: AbstractRODESystem
    "A"
    A::T1 = [2. 0.; 0 -2]
    "Right-hand-side function"
    righthandside::RH = (dx, x, u, t, W) -> (dx .= A * x * W)
    "Readout function"
    readout::RO = (x, u, t) -> x
    "State"
    state::ST = rand(2)
    "Input. Expected to be an `Inport` or `Nothing`"
    input::IP = nothing
    "Output port"
    output::OP = Outport(2)
end

##### Pretty printing
show(io::IO, ds::RODESystem) = print(io, "RODESystem(righthandside:$(ds.righthandside), ",
    "readout:$(ds.readout), state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::MultiplicativeNoiseLinearSystem) = print(io,
    "MultiplicativeNoiseLinearSystem(A:$(ds.A), state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")


text:
    import AlgebraicMultigrid: ruge_stuben, smoothed_aggregation,
    poisson, aspreconditioner

import IterativeSolvers: cg

function test_cycles()

    A = poisson((50,50))
    b = A * ones(size(A,1))

    reltol = 1e-8

    for method in [ruge_stuben, smoothed_aggregation]
        ml = method(A)

        for cycle in [AlgebraicMultigrid.V(),AlgebraicMultigrid.W(),AlgebraicMultigrid.F()]
            x,convhist = solve(ml, b, cycle; reltol = reltol, log = true)

            @debug "number of iterations for $cycle using $method: $(length(convhist))"
            @test norm(b - A*x) < reltol * norm(b)
        end

        for cycle in [AlgebraicMultigrid.V(),AlgebraicMultigrid.W(),AlgebraicMultigrid.F()]
            p = aspreconditioner(ml,cycle)

            x,convhist = cg(A, b, Pl = p; reltol = reltol, log = true)
            @debug "CG: number of iterations for $cycle using $method: $(convhist.iters)"
            @test norm(b - A*x) <= reltol * norm(b)
        end
    end
end


text:
# Test core utility functions
function t_coreutils()
	Xv = [0,1,2,3]
	Xm=[0 0 1 1; 1 1 2 2]
	y = [1,1,2,2]
	y2= [1 1 2 2;0 0 1 1]

	f=x->x

	D = j4pr.datacell(Xv,y)
	Du = j4pr.datacell(Xv)
	W = j4pr.functioncell(f)

	Test.@test j4pr.getx(D) == Xv
	Test.@test !(j4pr.getx(D) === Xv)
	Test.@test j4pr.getx!(D) == Xv
	Test.@test j4pr.getx!(D) === Xv

	Test.@test j4pr.gety(D) == y
	Test.@test !(j4pr.gety(D) === y)
	Test.@test j4pr.gety!(D) == y
	Test.@test j4pr.gety!(D) === y

	Test.@test j4pr.getf(W) == f
	Test.@test j4pr.getf(W) === f
	Test.@test j4pr.getf!(W) == f
	Test.@test j4pr.getf!(W) === f

	Test.@test j4pr.nvars(Xv) == 1
	Test.@test j4pr.nvars(Xm) == 2
	Test.@test j4pr.nvars(D) == 1

	Test.@test j4pr.strip(D) == (Xv,y)
	Test.@test j4pr.strip(Du) == Xv
	Test.@test j4pr.size(nothing) == 0
	Test.@test j4pr.size(nothing,rand(Int)) == 0
	Test.@test j4pr.size(D) ==j4pr.size(Xv)
	Test.@test j4pr.size(D,1) == j4pr.size(Xv,1)
	Test.@test j4pr.ndims(D) == 1

	Test.@test j4pr.isnan(nothing) == false
	Test.@test j4pr.isnan(randstring()) == false
	Test.@test j4pr.isnan([randstring() for i in 1:3]) == falses(3)
	Test.@test j4pr.isnan(DataArrays.data([randstring() for i in 1:3])) == falses(3)
	Test.@test j4pr.isnan([1,NaN,2]) == j4pr.isnan.([1,NaN,2])
	Test.@test j4pr.isvoid(nothing) == true
	Test.@test j4pr.isvoid(Any[1]) == false

	Test.@test j4pr.classsizes(Du) == DataStructures.SortedDict("unlabeled"=>4)
	Test.@test j4pr.classsizes(D) == DataStructures.SortedDict(1=>2,2=>2)
	Test.@test j4pr.classsizes(j4pr.datacell(Xv,y2)) == [DataStructures.SortedDict(1=>2,2=>2),DataStructures.SortedDict(0=>2,1=>2)]
	Test.@test j4pr.nclasssizes(Du) == DataStructures.SortedDict("unlabeled"=>1)
	Test.@test j4pr.nclasssizes(D) == DataStructures.SortedDict(1=>0.5,2=>0.5)
	Test.@test j4pr.nclasssizes(j4pr.datacell(Xv,y2)) == [DataStructures.SortedDict(1=>0.5,2=>0.5),DataStructures.SortedDict(0=>0.5,1=>0.5)]
	Test.@test j4pr.countmapn([1,2]) == DataStructures.SortedDict(1=>0.5,2=>0.5)
	Test.@test j4pr.countmapn(nothing) == DataStructures.SortedDict()
	Test.@test j4pr.classnames(Du) == []
	Test.@test j4pr.classnames(D) == [1,2]
	Test.@test j4pr.classnames(j4pr.datacell(Xv,y2)) == [[1,2],[0,1]]
	Test.@test j4pr.nclass(Du) == 0
	Test.@test j4pr.nclass(D) == 2
	Test.@test j4pr.nclass(j4pr.datacell(Xv,y2)) == [2,2]

end


text:
module DatesBenchmarks
# Based on https://github.com/quinnj/Dates.jl/blob/master/perf/perf.jl

using BenchmarkTools
using Dates

const SUITE = BenchmarkGroup()

g = addgroup!(SUITE, "construction")

g["Date"] = @benchmarkable Dates.Date(2014,1,1)
g["DateTime"] = @benchmarkable Dates.DateTime(2014,1,1,0,0,0,0)

g = addgroup!(SUITE, "accessor")

_date = Date(2014,7,19)
_datetime = DateTime(2014,7,19,15,13,12,512)

g["year"] = @benchmarkable Dates.year($_date)
g["month"] = @benchmarkable Dates.month($_date)
g["day"] = @benchmarkable Dates.day($_date)
g["hour"] = @benchmarkable Dates.hour($_datetime)
g["minute"] = @benchmarkable Dates.minute($_datetime)
g["second"] = @benchmarkable Dates.second($_datetime)
g["millisecond"] = @benchmarkable Dates.millisecond($_datetime)

g = addgroup!(SUITE, "conversion")

g["DateTime -> Date"] = @benchmarkable Dates.Date($_datetime)
g["Date -> DateTime"] = @benchmarkable Dates.DateTime($_date)

g = addgroup!(SUITE, "string")

g["Date"] = @benchmarkable string($_date)
g["DateTime"] = @benchmarkable string($_datetime)

g = addgroup!(SUITE, "query")

g["isleapyear","Date"]            = @benchmarkable Dates.isleapyear($_date)
g["isleapyear","DateTime"]        = @benchmarkable Dates.isleapyear($_datetime)
g["firstdayofmonth","Date"]       = @benchmarkable Dates.firstdayofmonth($_date)
g["firstdayofmonth","DateTime"]   = @benchmarkable Dates.firstdayofmonth($_datetime)
g["lastdayofmonth","Date"]        = @benchmarkable Dates.lastdayofmonth($_date)
g["lastdayofmonth","DateTime"]    = @benchmarkable Dates.lastdayofmonth($_datetime)
g["dayofweek","Date"]             = @benchmarkable Dates.dayofweek($_date)
g["dayofweek","DateTime"]         = @benchmarkable Dates.dayofweek($_datetime)
g["dayofweekofmonth","Date"]      = @benchmarkable Dates.dayofweekofmonth($_date)
g["dayofweekofmonth","DateTime"]  = @benchmarkable Dates.dayofweekofmonth($_datetime)
g["daysofweekinmonth","Date"]     = @benchmarkable Dates.daysofweekinmonth($_date)
g["daysofweekinmonth","DateTime"] = @benchmarkable Dates.daysofweekinmonth($_datetime)
g["firstdayofweek","Date"]        = @benchmarkable Dates.firstdayofweek($_date)
g["firstdayofweek","DateTime"]    = @benchmarkable Dates.firstdayofweek($_datetime)
g["lastdayofweek","Date"]         = @benchmarkable Dates.lastdayofweek($_date)
g["lastdayofweek","DateTime"]     = @benchmarkable Dates.lastdayofweek($_datetime)
g["dayofyear","Date"]             = @benchmarkable Dates.dayofyear($_date)
g["dayofyear","DateTime"]         = @benchmarkable Dates.dayofyear($_datetime)

for b in values(g)
    b.params.time_tolerance = 0.25
end

g = addgroup!(SUITE, "arithmetic")

_years = Dates.Year(3)
_months = Dates.Month(7)
_days = Dates.Day(78)
_hours = Dates.Hour(13)
_minutes = Dates.Minute(22)
_seconds = Dates.Second(22)
_milliseconds = Dates.Millisecond(333)

g["Date","Year"] = @benchmarkable $_date + $_years
g["DateTime","Year"] = @benchmarkable $_datetime + $_years
g["Date","Month"] = @benchmarkable $_date + $_months
g["DateTime","Month"] = @benchmarkable $_datetime + $_months
g["Date","Day"] = @benchmarkable $_date + $_days
g["DateTime","Day"] = @benchmarkable $_datetime + $_days
g["DateTime","Hour"] = @benchmarkable $_datetime + $_hours
g["DateTime","Minute"] = @benchmarkable $_datetime + $_minutes
g["DateTime","Second"] = @benchmarkable $_datetime + $_seconds
g["DateTime","Millisecond"] = @benchmarkable $_datetime + $_milliseconds

g = addgroup!(SUITE, "parse")

_datetime_str = "2016-02-19T12:34:56.78"
_date_str = "2016-02-19"

g["DateTime"] = @benchmarkable DateTime($_datetime_str)
g["Date"] = @benchmarkable Date($_date_str)

_custom_date_str = "20160219"
_custom_date_fmt = "yyyymmdd"
_custom_datetime_str = "20160219 123456.78"
_custom_datetime_fmt = "yyyymmdd HHMMSS.sss"

# Note: Also benchmarks the time it takes to create the custom format
g["Date","DateFormat"] = @benchmarkable Date($_custom_date_str, $_custom_date_fmt)
g["DateTime","DateFormat"] = @benchmarkable DateTime($_custom_datetime_str, $_custom_datetime_fmt)

g["DateTime","ISODateTimeFormat"] = @benchmarkable DateTime($_datetime_str, $(Dates.ISODateTimeFormat))
g["Date","ISODateFormat"] = @benchmarkable Date($_date_str, $(Dates.ISODateFormat))
g["DateTime","RFC1123Format","Titlecase"] = @benchmarkable DateTime("Sat, 12 Nov 2016 07:45:36", $(Dates.RFC1123Format))
g["DateTime","RFC1123Format","Lowercase"] = @benchmarkable DateTime("sat, 12 Nov 2016 07:45:36", $(Dates.RFC1123Format))
g["DateTime","RFC1123Format","Mixedcase"] = @benchmarkable DateTime("sAt, 12 Nov 2016 07:45:36", $(Dates.RFC1123Format))

end


text:
using JSON
function load!(net::Network,file_str::String)::Network
    dict = Dict()

    open(file_str, "r") do file
        global dict
        dict = JSON.parse(file)
    end

    net.weights = dict["weights"]
    net.biases = dict["biases"]

    return net
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
using RoboGrid
using Test

include("test/grid.jl")
include("test/robot.jl")
include("test/episode.jl")
include("test/reward.jl")
include("test/fitness.jl")


text:
using Femshop2
using Documenter

DocMeta.setdocmeta!(Femshop2, :DocTestSetup, :(using Femshop2); recursive=true)

makedocs(;
    modules=[Femshop2],
    authors="Aadesh, Eric, Hari",
    repo="https://github.com/aa25desh/Femshop2.jl/blob/{commit}{path}#{line}",
    sitename="Femshop2.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://aa25desh.github.io/Femshop2.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/aa25desh/Femshop2.jl",
)


text:
# Breadth-first search / traversal

#################################################
#
#  Breadth-first visit
#
#################################################

type BreadthFirst <: AbstractGraphVisitAlgorithm
end

function breadth_first_visit_impl!(
    graph::AbstractGraph,   # the graph
    queue,                  # an (initialized) queue that stores the active vertices
    colormap::Vector{Int},          # an (initialized) color-map to indicate status of vertices
    visitor::AbstractGraphVisitor)  # the visitor

    while !isempty(queue)
        u = dequeue!(queue)
        open_vertex!(visitor, u)

        for v in out_neighbors(u, graph)
            vi = vertex_index(v, graph)
            v_color::Int = colormap[vi]
            # TODO: Incorporate edge colors to BFS
            examine_neighbor!(visitor, u, v, v_color, -1)

            if v_color == 0
                colormap[vi] = 1
                if !discover_vertex!(visitor, v)
                    return
                end
                enqueue!(queue, v)
            end
        end

        colormap[vertex_index(u, graph)] = 2
        close_vertex!(visitor, u)
    end
    nothing
end

function traverse_graph{V,E}(
    graph::AbstractGraph{V,E},
    alg::BreadthFirst,
    s::V,
    visitor::AbstractGraphVisitor;
    colormap = zeros(Int, num_vertices(graph)))

    @graph_requires graph adjacency_list vertex_map

    que = Queue(V)

    colormap[vertex_index(s, graph)] = 1
    if !discover_vertex!(visitor, s)
        return
    end
    enqueue!(que, s)

    breadth_first_visit_impl!(graph, que, colormap, visitor)
end

function traverse_graph{V,E}(
    graph::AbstractGraph{V,E},
    alg::BreadthFirst,
    sources::AbstractVector{V},
    visitor::AbstractGraphVisitor;
    colormap = zeros(Int, num_vertices(graph)))

    @graph_requires graph adjacency_list vertex_map

    que = Queue(V)

    for s in sources
        colormap[vertex_index(s, graph)] = 1
        if !discover_vertex!(visitor, s)
            return
        end
        enqueue!(que, s)
    end

    breadth_first_visit_impl!(graph, que, colormap, visitor)
end

#################################################
#
#  Useful applications
#
#################################################

# Get the map of the (geodesic) distances from vertices to source by BFS

immutable GDistanceVisitor{G<:AbstractGraph} <: AbstractGraphVisitor
    graph::G
    dists::Vector{Int}
end

GDistanceVisitor{G<:AbstractGraph}(g::G, dists) = GDistanceVisitor{G}(g, dists)

function examine_neighbor!(visitor::GDistanceVisitor, u, v, vcolor::Int, ecolor::Int)
    if vcolor == 0
        g = visitor.graph
        dists = visitor.dists
        dists[vertex_index(v, g)] = dists[vertex_index(u, g)] + 1
    end
end

function gdistances!{V,E,DMap}(graph::AbstractGraph{V,E}, s::V, dists::DMap)
    visitor = GDistanceVisitor(graph, dists)
    dists[vertex_index(s, graph)] = 0
    traverse_graph(graph, BreadthFirst(), s, visitor)
    dists
end

function gdistances!{V,E,DMap}(graph::AbstractGraph{V,E}, sources::AbstractVector{V}, dists::DMap)
    visitor = GDistanceVisitor(graph, dists)
    for s in sources
        dists[vertex_index(s, graph)] = 0
    end
    traverse_graph(graph, BreadthFirst(), sources, visitor)
    dists
end

function gdistances(graph::AbstractGraph, sources; defaultdist::Int=-1)
    dists = fill(defaultdist, num_vertices(graph))
    gdistances!(graph, sources, dists)
end


text:
# Run lint Playground.jl for Error or Critical level message
msgs = lintpkg("Playground", returnMsgs=true)
for m in msgs
    @test m.level < 2
end


text:
module GLOQ

using LinearAlgebra
using DifferentialEquations

export about
# Write your package code here.
include("about.jl")
# make up the operators
include("QuantumOperator.jl")
# utility
include("utility.jl")
# Lindblad solver
include("ForwardLindbladSolver.jl")
# Provide interface to do characterization with Ramsey, Echo, T1 experiment
include("ExperimentAPI.jl")
# Phhysical constants and units
include("Constants.jl")
export hello_world

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Author: N�bia do Carmo Ferreira
#* $Id: HeatExchangerDetailed.mso 197 2007-03-08 14:31:57Z bicca $
#*------------------------------------------------------------------
# This model is valid if 14m^2 < A < 1100 m^2
type HeatExchanger_NTU_cost
	HeatExchanger_NTU_cost()=begin
		new(
			Heatex_NTU(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["Stainless steel 316", "Stainless steel 304", "Stainless steel 347", "Nickel 200", "Monel 400", "Inconel 600", "Incoloy 825", "Titanium", "Hastelloy"],
				:Default=>"Stainless steel 316"
			)),
			fill(DanaReal()),
			currency(Dict{Symbol,Any}(
				:Brief=>"Capital Cost"
			)),
			currency(Dict{Symbol,Any}(
				:Brief=>"Basic Cost"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Cost Factor based on the type of the heat exchanger"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Cost Factor based on the project pressure"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Cost Factor based on the construction material"
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Average  Pressure"
			)),
			[
				:(Pmax = max( [(0.5*_base_1._base_1.InletHot.P+0.5*_base_1._base_1.OutletHot.P) , (0.5*_base_1._base_1.InletCold.P+0.5*_base_1._base_1.OutletCold.P)] )),
				:(Ce = Cb*Fd*Fp*Fm),
				:(Cb = "US\$"*exp(Cost(1,1) + Cost(1,2)*ln(A/"m^2") + Cost(1,3)*(ln(A/"m^2"))^2)),
				:(Fd = exp(Cost(2,1) + Cost(2,2)*ln(A/"m^2") + Cost(2,3)*ln(A/"m^2"))),
				:(Fp = Cost(5,1) + Cost(5,2)*ln(A/"m^2")),
				:(Fp = Cost(5,1) + Cost(5,2)*ln(A/"m^2")),
				:(Fp = Cost(6,1) + Cost(6,2)*ln(A/"m^2")),
				:(Fp = Cost(7,1) + Cost(7,2)*ln(A/"m^2")),
				:(Fp = Cost(7,1) + Cost(7,2)*ln(A/"m^2")),
				:(Fm = Cost(8,1) + Cost(8,2)*ln(A/"m^2")),
			],
			[
				"Average pressure","Capital Cost","Basic Cost","Cost Factor based on the type of the heat exchanger","Cost Factor based on the project pressure","","","","","Cost Factor based on the construction material",
			],
			[:Material,:Cost,],
			[:Ce,:Cb,:Fd,:Fp,:Fm,:Pmax,]
		)
	end
	_base_1::Heatex_NTU
	Material::DanaSwitcher
	Cost::Array{DanaReal}
	Ce::currency
	Cb::currency
	Fd::positive
	Fp::positive
	Fm::positive
	Pmax::pressure
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export HeatExchanger_NTU_cost
function setEquationFlow(in::HeatExchanger_NTU_cost)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	if Pmax <= 700 * "kPa"
		# verificar
		addEquation(5)
	else
		if 700 * "kPa" < Pmax && Pmax < 2100 * "kPa"
			addEquation(6)
		else
			if 2100 * "kPa" < Pmax && Pmax < 4200 * "kPa"
				addEquation(7)
			else
				if 4200 * "kPa" < Pmax && Pmax < 6200 * "kPa"
					addEquation(8)
				else
					addEquation(9)
					# verificar

				end
			end
		end
	end
	addEquation(10)
end
function atributes(in::HeatExchanger_NTU_cost,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/HeatExchanger_NTU"
	fields[:Brief]="Heat Exchanger Block - NTU Method"
	fields[:Info]="to be documented."
	drive!(fields,_)
	return fields
end
HeatExchanger_NTU_cost(_::Dict{Symbol,Any})=begin
	newModel=HeatExchanger_NTU_cost()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
import Base.mean
function mean(ds::Array{Dipole})

    mean_x = mean([d.x for d in ds])
    mean_y = mean([d.y for d in ds])
    mean_z = mean([d.z for d in ds])
    mean_s = mean([d.size for d in ds])

    Dipole(ds[1].coord_system, mean_x, mean_y, mean_z, 0, 0, 0, 0, 0, mean_s)
end

import Base.std
function std(ds::Array{Dipole})

    std_x = std([d.x / Meter for d in ds])
    std_y = std([d.y / Meter for d in ds])
    std_z = std([d.z / Meter for d in ds])
    std_s = std([d.size for d in ds])

    Dipole(ds[1].coord_system, std_x, std_y, std_z, 0, 0, 0, 0, 0, std_s)
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# for converting tntp files with nodes and road network to format usable for JEMSS
# tntp files from https://github.com/bstabler/TransportationNetworks

using JEMSS

include("graph_tools.jl")

function readTntpNodesFile(tntpNodesFilename::String; delim::Char = '\t')
	# expected node file format (if tab delimited):
	# node	x	y
	# 1	153.2983311	-28.00211673	;
	# 2	153.5293645	-28.19705754	;
	# etc.

	# read nodes file
	println("Reading nodes file")
	data = readDlmFile(tntpNodesFilename; delim = delim)
	cols = findall(.!isempty.(data[1,:])) # need to ignore any empty headers
	table = Table("tntpNodes", data[1,cols], data[2:end,cols])
	numNodes = size(table.data, 1)
	nodes = Vector{Node}(undef, numNodes)
	for i = 1:numNodes
		nodes[i] = Node()
		nodes[i].index = table.columns["node"][i]
		nodes[i].location.x = table.columns["x"][i]
		nodes[i].location.y = table.columns["y"][i]
		@assert(nodes[i].index == i)
	end

	return nodes
end

function readTntpNetworkFile(tntpNetworkFilename::String; delim::Char = '\t')
	# expected network file format (if tab delimited):
	# <NUMBER OF ZONES> 1068
	# <NUMBER OF NODES> 4807
	# <FIRST THRU NODE> 1069
	# <NUMBER OF LINKS> 11140
	# [...intermediate lines...]
	# ~	from	to	capacity (per lane)	length (km)	fftt (min)	alpha	beta	ff speed (km/h)	critical speed (km/h)	lanes
	#	 1	1371	900	0.300	0.327	0.282	4	55	42.9	2	;
	#	 2	2012	1600	0.260	0.173	0.282	4	90	70.2	2	;
	# etc.

	# read network file
	println("Reading network file")
	data = readDlmFile(tntpNetworkFilename; delim = delim)

	# first few lines have metadata
	f(r::Regex, s::AbstractString) = Meta.parse(match(r,s).captures[1])
	numZones = f(r"<NUMBER OF ZONES> (\d+)", data[1,1])
	numNodes = f(r"<NUMBER OF NODES> (\d+)", data[2,1])
	firstThruNode = f(r"<FIRST THRU NODE> (\d+)", data[3,1])
	@assert(firstThruNode == numZones + 1)
	numArcs = f(r"<NUMBER OF LINKS> (\d+)", data[4,1])

	# find and read table data
	row = first(findall(data[:,1] .== "~"))
	header = data[row,:]
	cols = findall(.!isempty.(header))[2:end] # need to ignore any empty headers, and "~"
	table = Table("tntpNetwork", data[row, cols], data[row+1:end, cols])
	lengthColString = first(filter(s -> occursin("length", s), table.header)) # to match "length", "length (km)", etc
	rowsFields = tableRowsFieldDicts(table, setdiff(table.header, ["from", "to"])) # will keep "length" field
	numArcs = size(table.data,1)
	arcs = Vector{Arc}(undef, numArcs)
	travelTimes = Vector{Float}(undef, numArcs)
	for i = 1:numArcs
		arcs[i] = Arc()
		arcs[i].index = i
		arcs[i].fromNodeIndex = table.columns["from"][i]
		arcs[i].toNodeIndex = table.columns["to"][i]
		arcs[i].distance = table.columns[lengthColString][i]
		arcs[i].fields = rowsFields[i]
	end

	return arcs, numZones
end

# Process nodes and arcs from tntp, processing performed:
# - remove all zones and their connected arcs
# - keep the largest strongly connected component
function convertTntpNetwork!(nodes::Vector{Node}, arcs::Vector{Arc}, numZones::Int)

	@assert(!graphContainsDuplicateArcs(arcs))

	# remove all zones and their connected arcs
	nodeFilter(node::Node) = node.index > numZones
	arcFilter(arc::Arc) = arc.fromNodeIndex > numZones && arc.toNodeIndex > numZones
	graphRemoveElts!(nodes, arcs; nodeFilter = nodeFilter, arcFilter = arcFilter)

	graphRemoveDisconnectedArcs!(nodes, arcs) # just in case, but should not be necessary
	graphKeepLargestComponent!(nodes, arcs)
end

function convertTntpNetworkFile(tntpNodesFilename::String, tntpNetworkFilename::String;
	delim::Char = '\t')
	nodes = readTntpNodesFile(tntpNodesFilename; delim = delim)
	(arcs, numZones) = readTntpNetworkFile(tntpNetworkFilename; delim = delim)
	convertTntpNetwork!(nodes, arcs, numZones)
	return nodes, arcs
end


text:
"""
    abstract type AbstractScalarDiffusivity <: AbstractTurbulenceClosure end

Abstract type for closures with *isotropic* diffusivities.
"""
abstract type AbstractScalarDiffusivity{TD, F} <: AbstractTurbulenceClosure{TD} end

"""
    struct ThreeDimensional end

Specifies a three-dimensionally-isotropic `ScalarDiffusivity`.
"""
struct ThreeDimensionalFormulation end

"""
    struct Horizontal end

Specifies a horizontally-isotropic, `VectorInvariant`, `ScalarDiffusivity`.
"""
struct HorizontalFormulation end

"""
    struct Vertical end

Specifies a `ScalarDiffusivity` acting only in the vertical direction.
"""
struct VerticalFormulation end

"""
    viscosity(closure, diffusivities)

Returns the scalar viscosity associated with `closure`.
"""
function viscosity end

"""
    diffusivity(closure, tracer_index, diffusivity_fields)

Returns the scalar diffusivity associated with `closure` and `tracer_index`.
"""
function diffusivity end

# Fallback locations
@inline viscosity_location(::AbstractScalarDiffusivity) = (Center(), Center(), Center())
@inline diffusivity_location(::AbstractScalarDiffusivity) = (Center(), Center(), Center())

# For tuples (note that kernel functions are "untupled", so these are for the user API)
viscosity(closure::Tuple, K) = sum(Tuple(viscosity(closure[n], K[n]) for n = 1:length(closure)))
diffusivity(closure::Tuple, K, id) = sum(Tuple(diffusivity(closure[n], K[n], id) for n = 1:length(closure)))

@inline formulation(::AbstractScalarDiffusivity{TD, F}) where {TD, F} = F()

Base.summary(::VerticalFormulation) = "VerticalFormulation"
Base.summary(::HorizontalFormulation) = "HorizontalFormulation"
Base.summary(::ThreeDimensionalFormulation) = "ThreeDimensionalFormulation"

#####
##### Coefficient extractors
#####

const ASD = AbstractScalarDiffusivity

@inline νᶜᶜᶜ(i, j, k, grid, clo::ASD, K, clock)     = νᶜᶜᶜ(i, j, k, grid, clock, viscosity_location(clo), viscosity(clo, K))
@inline νᶠᶠᶜ(i, j, k, grid, clo::ASD, K, clock)     = νᶠᶠᶜ(i, j, k, grid, clock, viscosity_location(clo), viscosity(clo, K))
@inline νᶠᶜᶠ(i, j, k, grid, clo::ASD, K, clock)     = νᶠᶜᶠ(i, j, k, grid, clock, viscosity_location(clo), viscosity(clo, K))
@inline νᶜᶠᶠ(i, j, k, grid, clo::ASD, K, clock)     = νᶜᶠᶠ(i, j, k, grid, clock, viscosity_location(clo), viscosity(clo, K))

@inline κᶠᶜᶜ(i, j, k, grid, clo::ASD, K, id, clock) = κᶠᶜᶜ(i, j, k, grid, clock, diffusivity_location(clo), diffusivity(clo, K, id))
@inline κᶜᶠᶜ(i, j, k, grid, clo::ASD, K, id, clock) = κᶜᶠᶜ(i, j, k, grid, clock, diffusivity_location(clo), diffusivity(clo, K, id))
@inline κᶜᶜᶠ(i, j, k, grid, clo::ASD, K, id, clock) = κᶜᶜᶠ(i, j, k, grid, clock, diffusivity_location(clo), diffusivity(clo, K, id))

@inline νzᶜᶜᶜ(i, j, k, grid, clo::ASD, K, clock)     = νᶜᶜᶜ(i, j, k, grid, clock, viscosity_location(clo), viscosity(clo, K))
@inline νzᶠᶠᶜ(i, j, k, grid, clo::ASD, K, clock)     = νᶠᶠᶜ(i, j, k, grid, clock, viscosity_location(clo), viscosity(clo, K))
@inline νzᶠᶜᶠ(i, j, k, grid, clo::ASD, K, clock)     = νᶠᶜᶠ(i, j, k, grid, clock, viscosity_location(clo), viscosity(clo, K))
@inline νzᶜᶠᶠ(i, j, k, grid, clo::ASD, K, clock)     = νᶜᶠᶠ(i, j, k, grid, clock, viscosity_location(clo), viscosity(clo, K))

@inline κzᶠᶜᶜ(i, j, k, grid, clo::ASD, K, id, clock) = κᶠᶜᶜ(i, j, k, grid, clock, diffusivity_location(clo), diffusivity(clo, K, id))
@inline κzᶜᶠᶜ(i, j, k, grid, clo::ASD, K, id, clock) = κᶜᶠᶜ(i, j, k, grid, clock, diffusivity_location(clo), diffusivity(clo, K, id))
@inline κzᶜᶜᶠ(i, j, k, grid, clo::ASD, K, id, clock) = κᶜᶜᶠ(i, j, k, grid, clock, diffusivity_location(clo), diffusivity(clo, K, id))

#####
##### Stress divergences
#####

const AID = AbstractScalarDiffusivity{<:Any, <:ThreeDimensionalFormulation}
const AHD = AbstractScalarDiffusivity{<:Any, <:HorizontalFormulation}
const AVD = AbstractScalarDiffusivity{<:Any, <:VerticalFormulation}

@inline viscous_flux_ux(i, j, k, grid, closure::AID, K, U, C, clock, b) = - 2 * ν_σᶜᶜᶜ(i, j, k, grid, closure, K, clock, Σ₁₁, U.u, U.v, U.w)
@inline viscous_flux_vx(i, j, k, grid, closure::AID, K, U, C, clock, b) = - 2 * ν_σᶠᶠᶜ(i, j, k, grid, closure, K, clock, Σ₂₁, U.u, U.v, U.w)
@inline viscous_flux_wx(i, j, k, grid, closure::AID, K, U, C, clock, b) = - 2 * ν_σᶠᶜᶠ(i, j, k, grid, closure, K, clock, Σ₃₁, U.u, U.v, U.w)
@inline viscous_flux_uy(i, j, k, grid, closure::AID, K, U, C, clock, b) = - 2 * ν_σᶠᶠᶜ(i, j, k, grid, closure, K, clock, Σ₁₂, U.u, U.v, U.w)
@inline viscous_flux_vy(i, j, k, grid, closure::AID, K, U, C, clock, b) = - 2 * ν_σᶜᶜᶜ(i, j, k, grid, closure, K, clock, Σ₂₂, U.u, U.v, U.w)
@inline viscous_flux_wy(i, j, k, grid, closure::AID, K, U, C, clock, b) = - 2 * ν_σᶜᶠᶠ(i, j, k, grid, closure, K, clock, Σ₃₂, U.u, U.v, U.w)

@inline viscous_flux_ux(i, j, k, grid, closure::AHD, K, U, C, clock, b) = - ν_δᶜᶜᶜ(i, j, k, grid, closure, K, clock, U.u, U.v)
@inline viscous_flux_vx(i, j, k, grid, closure::AHD, K, U, C, clock, b) = - ν_ζᶠᶠᶜ(i, j, k, grid, closure, K, clock, U.u, U.v)
@inline viscous_flux_uy(i, j, k, grid, closure::AHD, K, U, C, clock, b) = + ν_ζᶠᶠᶜ(i, j, k, grid, closure, K, clock, U.u, U.v)
@inline viscous_flux_vy(i, j, k, grid, closure::AHD, K, U, C, clock, b) = - ν_δᶜᶜᶜ(i, j, k, grid, closure, K, clock, U.u, U.v)
@inline viscous_flux_wx(i, j, k, grid, closure::AHD, K, U, C, clock, b) = - ν_σᶠᶜᶠ(i, j, k, grid, closure, K, clock, ∂xᶠᶜᶠ, U.w)
@inline viscous_flux_wy(i, j, k, grid, closure::AHD, K, U, C, clock, b) = - ν_σᶜᶠᶠ(i, j, k, grid, closure, K, clock, ∂yᶜᶠᶠ, U.w)

@inline viscous_flux_uz(i, j, k, grid, closure::AID, K, U, C, clock, b) = - 2 * ν_σᶠᶜᶠ(i, j, k, grid, closure, K, clock, Σ₁₃, U.u, U.v, U.w)
@inline viscous_flux_vz(i, j, k, grid, closure::AID, K, U, C, clock, b) = - 2 * ν_σᶜᶠᶠ(i, j, k, grid, closure, K, clock, Σ₂₃, U.u, U.v, U.w)
@inline viscous_flux_wz(i, j, k, grid, closure::AID, K, U, C, clock, b) = - 2 * ν_σᶜᶜᶜ(i, j, k, grid, closure, K, clock, Σ₃₃, U.u, U.v, U.w)

@inline viscous_flux_uz(i, j, k, grid, closure::AVD, K, U, C, clock, b) = - ν_σᶠᶜᶠ(i, j, k, grid, closure, K, clock, ∂zᶠᶜᶠ, U.u)
@inline viscous_flux_vz(i, j, k, grid, closure::AVD, K, U, C, clock, b) = - ν_σᶜᶠᶠ(i, j, k, grid, closure, K, clock, ∂zᶜᶠᶠ, U.v)
@inline viscous_flux_wz(i, j, k, grid, closure::AVD, K, U, C, clock, b) = - ν_σᶜᶜᶜ(i, j, k, grid, closure, K, clock, ∂zᶜᶜᶜ, U.w)

#####
##### Diffusive fluxes
#####

const AIDorAHD = Union{AID, AHD}
const AIDorAVD = Union{AID, AVD}

@inline diffusive_flux_x(i, j, k, grid, cl::AIDorAHD, K, ::Val{id}, U, C, clk, b) where id = - κᶠᶜᶜ(i, j, k, grid, cl, K, Val(id), clk) * ∂xᶠᶜᶜ(i, j, k, grid, C[id])
@inline diffusive_flux_y(i, j, k, grid, cl::AIDorAHD, K, ::Val{id}, U, C, clk, b) where id = - κᶜᶠᶜ(i, j, k, grid, cl, K, Val(id), clk) * ∂yᶜᶠᶜ(i, j, k, grid, C[id])
@inline diffusive_flux_z(i, j, k, grid, cl::AIDorAVD, K, ::Val{id}, U, C, clk, b) where id = - κᶜᶜᶠ(i, j, k, grid, cl, K, Val(id), clk) * ∂zᶜᶜᶠ(i, j, k, grid, C[id])

#####
##### Zero out not used fluxes
#####

for (dir, closure) in zip((:x, :y, :z), (:AVD, :AVD, :AHD))
    diffusive_flux = Symbol(:diffusive_flux_, dir)
    viscous_flux_u = Symbol(:viscous_flux_u, dir)
    viscous_flux_v = Symbol(:viscous_flux_v, dir)
    viscous_flux_w = Symbol(:viscous_flux_w, dir)
    @eval begin
        @inline $diffusive_flux(i, j, k, grid, closure::$closure, K, ::Val, args...) = zero(eltype(grid))
        @inline $viscous_flux_u(i, j, k, grid, closure::$closure, args...) = zero(eltype(grid))
        @inline $viscous_flux_v(i, j, k, grid, closure::$closure, args...) = zero(eltype(grid))
        @inline $viscous_flux_w(i, j, k, grid, closure::$closure, args...) = zero(eltype(grid))
    end
end

#####
##### Support for VerticallyImplicit
#####

const VITD = VerticallyImplicitTimeDiscretization

@inline ivd_viscous_flux_uz(i, j, k, grid, closure, K, U, C, clock, b) = - ν_σᶠᶜᶠ(i, j, k, grid, closure, K, clock, ∂xᶠᶜᶠ, U.w)
@inline ivd_viscous_flux_vz(i, j, k, grid, closure, K, U, C, clock, b) = - ν_σᶜᶠᶠ(i, j, k, grid, closure, K, clock, ∂yᶜᶠᶠ, U.w)

# General functions (eg for vertically periodic)
@inline viscous_flux_uz(i, j, k, grid,  ::VITD, closure::AIDorAVD, args...) = ivd_viscous_flux_uz(i, j, k, grid, closure, args...)
@inline viscous_flux_vz(i, j, k, grid,  ::VITD, closure::AIDorAVD, args...) = ivd_viscous_flux_vz(i, j, k, grid, closure, args...)
@inline viscous_flux_wz(i, j, k, grid,  ::VITD, closure::AIDorAVD, args...) = zero(eltype(grid))
@inline diffusive_flux_z(i, j, k, grid, ::VITD, closure::AIDorAVD, args...) = zero(eltype(grid))

# Vertically bounded grids
#
# For vertically bounded grids, we calculate _explicit_ fluxes on the boundaries,
# and elide the implicit vertical flux component everywhere else. This is consistent
# with the formulation of the tridiagonal solver, which requires explicit treatment
# of boundary contributions (eg boundary contributions must be moved to the right
# hand side of the tridiagonal system).

@inline function viscous_flux_uz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, args...)
    return ifelse((k == 1) | (k == grid.Nz+1),
                  viscous_flux_uz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  ivd_viscous_flux_uz(i, j, k, grid, closure, args...))
end

@inline function viscous_flux_vz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, args...)
    return ifelse((k == 1) | (k == grid.Nz+1),
                  viscous_flux_vz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  ivd_viscous_flux_vz(i, j, k, grid, closure, args...))
end

@inline function viscous_flux_wz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, args...)
    return ifelse((k == 1) | (k == grid.Nz+1),
                  viscous_flux_wz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  zero(eltype(grid)))
end

@inline function diffusive_flux_z(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AIDorAVD, args...)
    return ifelse((k == 1) | (k == grid.Nz+1),
                  diffusive_flux_z(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  zero(eltype(grid)))
end


text:
# strange error in simple network
# Claudia November 2015

include("../src/types.jl")
include("../src/functions.jl")

tree= "(A,((B,#H1),(C,(D)#H1)));" #fails
tree= "(A,((C,(D)#H1),(B,#H1)));" #works
net = readTopology(tree)
net = readTopologyUpdate(tree)
printEdges(net)
printNodes(net)


text:
function solve(
    instance::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
    ::DynamicProgramming,
) where {T, U}
    # A dynamic-programming stance on Prim's algorithm.

    # Base case: one edge, the solution is either to take it or not based
    # on the weight.
    if ne(instance.instance.graph) == 1
        e = collect(edges(instance.instance.graph))[1]
        v = instance.weights[e]
        w = reward(instance.instance, e)

        S = Dict{Int, Vector{Edge{T}}}(β => ifelse(w < β, Edge{T}[], [e]) for β in 0:(instance.min_budget))

        return BudgetedSpanningTreeDynamicProgrammingSolution(
            instance,
            S[instance.min_budget],
            S,
        )
    end

    # Get an indexed list of *all* edges, sorted by increasing value.
    sorted_edges = Tuple{Edge, Float64, Int}[]
    for e in edges(instance.instance.graph)
        push!(
            sorted_edges,
            (e, reward(instance.instance, e), instance.weights[e]),
        )
    end
    sort!(sorted_edges, by=(p) -> -p[2])

    # Fill the DP table, considering edges one by one (from the best to the
    # poorest), Prim-like.
    # Index:
    # - min budget, starting at 0
    # - considering edges up to that index (included), starting at 1
    V = Dict{Tuple{T, Int}, Float64}()
    S = Dict{Tuple{T, Int}, Vector{Edge{T}}}()
    LD = Dict{Tuple{T, Int}, LoopDetector}()

    for β in 0:(instance.min_budget)
        for v in 1:ne(instance.instance.graph)
            LD[β, v] = LoopDetector(nv(instance.instance.graph))
        end
    end

    first_edge, first_edge_value, first_edge_weight = sorted_edges[1]
    V[0, 1] = first_edge_value
    S[0, 1] = Edge{T}[first_edge]
    visit_edge(LD[0, 1], src(first_edge), dst(first_edge))

    # Initialisation for a zero min-budget (i.e. no constraint).
    for i in 2:ne(instance.instance.graph)
        # Can edges_sorted[i] be taken?
        edge, edge_value, edge_weight = sorted_edges[i]

        if edge_would_create_loop(LD[0, i - 1], src(edge), dst(edge)) # Don't take i.
            V[0, i] = V[0, i - 1]
            S[0, i] = S[0, i - 1]
            copy!(LD[0, i], LD[0, i - 1])
        else # Take i.
            V[0, i] = V[0, i - 1] + edge_value
            S[0, i] = copy(S[0, i - 1])
            push!(S[0, i], edge)

            copy!(LD[0, i], LD[0, i - 1])
            visit_edge(LD[0, i], src(edge), dst(edge))
        end
    end

    # Dynamic part.
    for β in 1:(instance.min_budget)
        # First edge (base case).
        if first_edge_weight >= β
            # The first edge is allowable on its own: take it.
            V[β, 1] = V[0, 1]
            S[β, 1] = S[0, 1]

            visit_edge(LD[β, 1], src(first_edge), dst(first_edge))
        else
            # The first edge is not allowable on its own, it does not have
            # a sufficient weight.
            V[β, 1] = -1
            S[β, 1] = Edge{T}[]
        end

        # Other edges (recursive cases).
        for i in 2:ne(instance.instance.graph)
            # Can edges_sorted[i] be taken?
            edge, edge_value, edge_weight = sorted_edges[i]

            remaining_budget = β - edge_weight
            if remaining_budget < 0
                remaining_budget = 0
            end

            if V[remaining_budget, i - 1] == -1 || length(S[remaining_budget, i - 1]) == 0
                # Are the previous solutions infeasible (recursive call)? If so,
                # don't consider using it.
                V[β, i] = V[β, i - 1]
                S[β, i] = S[β, i - 1]
            elseif edge_would_create_loop(LD[remaining_budget, i - 1], edge)
                # Don't take i with that previous solution, it would create a
                # loop.
                # Rather, try to build solutions without one of the ends of i,
                # to ensure there will be no loop when combining. Use the
                # remaining budget
                g_up_to_i = copy(instance.instance.graph)
                for j in i:ne(instance.instance.graph)
                    rem_edge!(g_up_to_i, sorted_edges[j][1])
                end

                g_no_src = copy(g_up_to_i)
                rem_vertex!(g_no_src, src(edge)) # Changes indices!
                @assert nv(g_no_src) < nv(g_up_to_i)
                @assert ne(g_no_src) < ne(g_up_to_i)
                map_index = (idx) -> ifelse(idx < src(edge), idx, idx - 1)
                map_edge = (e) -> Edge(map_index(src(e)), map_index(dst(e)))
                unmap_index = (idx) -> ifelse(idx < src(edge), idx, idx + 1)
                unmap_edge =
                    (e) -> Edge(unmap_index(src(e)), unmap_index(dst(e)))
                r_no_src = Dict(
                    filter(map(map_edge, collect(keys(instance.instance.rewards))) .=> values(instance.instance.rewards)) do r
                        r[1] ∈ edges(g_no_src)
                    end
                )
                w_no_src = Dict(
                    filter(map(map_edge, collect(keys(instance.weights))) .=> values(instance.weights)) do r
                        r[1] ∈ edges(g_no_src)
                    end
                )
                i_no_src = MinimumBudget(
                    SpanningTreeInstance(g_no_src, r_no_src),
                    w_no_src,
                    remaining_budget,
                )
                s_no_src = solve(i_no_src, DynamicProgramming())
                t_no_src = Edge{T}[unmap_edge(e) for e in s_no_src.variables]
                v_no_src =
                    _budgeted_spanning_tree_compute_value(instance, t_no_src)

                g_no_dst = copy(g_up_to_i)
                rem_vertex!(g_no_dst, dst(edge)) # Changes indices!
                @assert nv(g_no_dst) < nv(g_up_to_i)
                @assert ne(g_no_dst) < ne(g_up_to_i)
                map_index = (idx) -> ifelse(idx < dst(edge), idx, idx - 1)
                map_edge = (e) -> Edge(map_index(src(e)), map_index(dst(e)))
                unmap_index = (idx) -> ifelse(idx < dst(edge), idx, idx + 1)
                unmap_edge =
                    (e) -> Edge(unmap_index(src(e)), unmap_index(dst(e)))
                r_no_dst = Dict(
                    filter(map(map_edge, collect(keys(instance.instance.rewards))) .=> values(instance.instance.rewards)) do r
                        r[1] ∈ edges(g_no_dst)
                    end
                )
                w_no_dst = Dict(
                    filter(map(map_edge, collect(keys(instance.weights))) .=> values(instance.weights)) do r
                        r[1] ∈ edges(g_no_dst)
                    end
                )
                i_no_dst = MinimumBudget(
                    SpanningTreeInstance(g_no_dst, r_no_dst),
                    w_no_dst,
                    remaining_budget,
                )
                s_no_dst = solve(i_no_dst, DynamicProgramming())
                t_no_dst = Edge{T}[unmap_edge(e) for e in s_no_dst.variables]
                v_no_dst =
                    _budgeted_spanning_tree_compute_value(instance, t_no_dst)

                if length(t_no_src) == 0 && length(t_no_dst) == 0
                    V[β, i] = V[β, i - 1]
                    S[β, i] = S[β, i - 1]
                    copy!(LD[β, i], LD[β, i - 1])
                elseif length(t_no_src) == 0 ||
                       (length(t_no_src) > 0 && v_no_dst >= v_no_src)
                    # Includes the case where both subproblems have the same
                    # value. Pick to continue with t_no_src arbitrarily.
                    V[β, i] = v_no_dst + edge_value
                    S[β, i] = t_no_dst
                    push!(S[β, i], edge)

                    # No available data structure to copy, as this solution
                    # has not been built within the dynamic-programming
                    # mechanism, regenerate it from scratch.
                    for e in S[β, i]
                        visit_edge(LD[β, i], e)
                    end
                elseif length(t_no_dst) == 0 ||
                       (length(t_no_dst) > 0 && v_no_src > v_no_dst)
                    V[β, i] = v_no_src + edge_value
                    S[β, i] = t_no_src
                    push!(S[β, i], edge)

                    # No available data structure to copy, regenerate it from
                    # scratch.
                    for e in S[β, i]
                        visit_edge(LD[β, i], e)
                    end
                else
                    error("Assertion failed")
                end
            else
                # Take i, there is no chance of it adding a loop.
                V[β, i] = V[remaining_budget, i - 1] + edge_value
                S[β, i] = copy(S[remaining_budget, i - 1])
                push!(S[β, i], edge)

                copy!(LD[β, i], LD[remaining_budget, i - 1])
                visit_edge(LD[β, i], edge)
            end
        end
    end

    S_refined =
        Dict(β => S[β, dimension(instance)] for β in 0:(instance.min_budget))
    return BudgetedSpanningTreeDynamicProgrammingSolution(
        instance,
        S[instance.min_budget, dimension(instance)],
        S_refined,
    )
end