text:
mutable struct Lattice{D,N}
    size::NTuple{D, Int} #linear extent of the lattice in number of unit cells
    length::Int #Number of sites N_sites
    unitcell::UnitCell{D}
    sitePositions::Vector{NTuple{D,Float64}}

    spins::Matrix{Float64} #3*N_sites matrix containing the spin configuration

    interactionSites::Vector{NTuple{N,Int}} #list of length N_sites, for every site contains all interacting sites
    interactionMatrices::Vector{NTuple{N,InteractionMatrix}} #list of length N_sites, for every site contains all interaction matrices
    interactionOnsite::Vector{InteractionMatrix} #list of length N_sites, for every site contains the local onsite interaction matrix
    interactionField::Vector{NTuple{3,Float64}} #list of length N_sites, for every site contains the local field
    Lattice(D,N) = new{D,N}()
end

function Lattice(uc::UnitCell{D}, L::NTuple{D,Int}) where D
    #parse interactions
    ##For every basis site b, generate list of sites which b interacts with and store the corresponding interaction sites and matrices.
    ##Interaction sites are specified by the target site's basis id, b_target, and the offset in units of primitive lattice vectors.
    ##If b has multiple interactions defined with the same target site, eliminate those duplicates by summing up the interaction matrices.
    interactionTargetSites = [ Vector{Tuple{Int,NTuple{D,Int},Matrix{Float64}}}(undef,0) for i in 1:length(uc.basis) ] #tuples of (b_target, offset, M)
    for x in uc.interactions
        b1, b2, offset, M = x
        b1 == b2 && offset == Tuple(zeros(D)) && error("Interaction cannot be local. Use setInteractionOnsite!() instead.")

        #locate existing coupling to target site and add interaction matrix
        for i in 1:length(interactionTargetSites[b1])
            if interactionTargetSites[b1][i][1] == b2 && interactionTargetSites[b1][i][2] == offset
                interactionTargetSites[b1][i] = (interactionTargetSites[b1][i][1], interactionTargetSites[b1][i][2], interactionTargetSites[b1][i][3] + M)
                @goto endb1
            end
        end
        #if coupling does not exist yet, push new entry
        push!(interactionTargetSites[b1], (b2, offset, M))
        @label endb1

        #locate existing coupling from target site and add interaction matrix
        for i in 1:length(interactionTargetSites[b2])
            if interactionTargetSites[b2][i][1] == b1 && interactionTargetSites[b2][i][2] == (x->-x).(offset)
                interactionTargetSites[b2][i] = (interactionTargetSites[b2][i][1], interactionTargetSites[b2][i][2], interactionTargetSites[b2][i][3] + transpose(M))
                @goto endb2
            end
        end
        #if coupling does not exist yet, push new entry
        push!(interactionTargetSites[b2], (b1, (x->-x).(offset), transpose(M)))
        @label endb2
    end
    Ninteractions = findmax([ length(interactionTargetSites[i]) for i in 1:length(uc.basis) ])[1]

    #create lattice struct
    lattice = Lattice(D,Ninteractions)
    lattice.size = L
    lattice.length = prod(L) * length(uc.basis)
    lattice.unitcell = uc

    #generate linear representation of lattice sites to assign integer site IDs
    ##Enumeration sequence is (a1, a2, ..., b) in row-major fashion
    sites = Vector{NTuple{D+1,Int}}(undef, lattice.length)
    function nextSite(site)
        next = collect(site)
        next[D+1] += 1
        if next[D+1] > length(uc.basis)
            next[D+1] = 1
            next[D] += 1
        end
        for d in reverse(1:D)
            if next[d] >= L[d]
                next[d] = 0
                d-1 > 0 && (next[d-1] += 1)
            end
        end
        return Tuple(next)
    end
    sites[1] = tuple(zeros(Int,D)..., 1)
    for i in 2:length(sites)
        sites[i] = nextSite(sites[i-1])
    end

    #init site positions
    lattice.sitePositions = Vector{NTuple{D,Float64}}(undef, length(sites))
    for i in 1:length(sites)
        site = sites[i]
        lattice.sitePositions[i] = .+([uc.primitive[j] .* site[j] for j in 1:D]...) .+ uc.basis[site[end]]
    end

    #init spins
    lattice.spins = Array{Float64,2}(undef, 3, length(sites))

    #write interactions to lattice
    lattice.interactionSites = repeat([ NTuple{Ninteractions,Int}(ones(Int,Ninteractions)) ], lattice.length)
    lattice.interactionMatrices = repeat([ NTuple{Ninteractions,InteractionMatrix}(repeat([InteractionMatrix()],Ninteractions)) ], lattice.length)
    lattice.interactionOnsite = repeat([InteractionMatrix()], lattice.length)
    lattice.interactionField = repeat([(0.0,0.0,0.0)], lattice.length)

    function applyPBC(n, L)
        while n < 0; n += L end
        while n >= L; n -= L end
        return n
    end
    function siteIndexFromParametrization(site)
       return findfirst(isequal(site), sites)
    end

    for i in 1:length(sites)
        site = sites[i]
        b = site[end]

        #onsite interaction
        lattice.interactionOnsite[i] = InteractionMatrix(uc.interactionsOnsite[b])

        #field interaction
        lattice.interactionField[i] = NTuple{3,Float64}(uc.interactionsField[b])

        #two-spin interactions
        interactionSites = repeat([i], Ninteractions)
        interactionMatrices = repeat([InteractionMatrix()], Ninteractions)
        for j in 1:Ninteractions
            if j <= length(interactionTargetSites[b])
                b2, offset, M = interactionTargetSites[b][j]

                primitiveTarget = [applyPBC(site[k] + offset[k], L[k]) for k in 1:D]
                targetSite = tuple(primitiveTarget..., b2)

                interactionSites[j] = siteIndexFromParametrization(targetSite)
                interactionMatrices[j] = InteractionMatrix(M)
            end
        end
        lattice.interactionSites[i] = NTuple{Ninteractions,Int}(interactionSites)
        lattice.interactionMatrices[i] = NTuple{Ninteractions,InteractionMatrix}(interactionMatrices)
    end

    #return lattice
    return lattice
end

function Base.size(lattice::Lattice{D,N}) where {D,N}
    return lattice.size
end

function Base.length(lattice::Lattice{D,N}) where {D,N}
    return lattice.length
end

function getSpin(lattice::Lattice{D,N}, site::Int) where {D,N}
    return (lattice.spins[1,site], lattice.spins[2,site], lattice.spins[3,site])
end

function setSpin!(lattice::Lattice{D,N}, site::Int, newState::Tuple{Float64,Float64,Float64}) where {D,N}
    lattice.spins[1,site] = newState[1]
    lattice.spins[2,site] = newState[2]
    lattice.spins[3,site] = newState[3]
end

function getSitePosition(lattice::Lattice{D,N}, site::Int)::NTuple{D,Float64} where {D,N}
    return lattice.sitePositions[site]
end

function getInteractionSites(lattice::Lattice{D,N}, site::Int)::NTuple{N,Int} where {D,N}
    return lattice.interactionSites[site]
end

function getInteractionMatrices(lattice::Lattice{D,N}, site::Int)::NTuple{N,InteractionMatrix} where {D,N}
    return lattice.interactionMatrices[site]
end

function getInteractionOnsite(lattice::Lattice{D,N}, site::Int)::InteractionMatrix where {D,N}
    return lattice.interactionOnsite[site]
end

function getInteractionField(lattice::Lattice{D,N}, site::Int)::NTuple{3,Float64} where {D,N}
    return lattice.interactionField[site]
end


text:
abstract type AbstractVectorSpace end
struct TrivialVectorSpace <:AbstractVectorSpace
    n :: Int
end
struct VectorSpace{S} <: AbstractVectorSpace
    sectors :: SortedDict{S, Int}
    isdual  :: Bool
    function VectorSpace{S}(s, id) where {S}
        allunique([c for (c,d) in s]) ||
            error("duplicate charges! $s")
        sectors = SortedDict{S, Int}(s)
        for (c, d) in sectors
            d > 0 || error("space dimension can only be positive!")
        end
        new{S}(sectors, id)
    end
    VectorSpace{S}(s) where{S} = VectorSpace{S}(s, false)
end
VectorSpace{S}(pairs::Pair{S,Int}...) where {S} = VectorSpace{S}(pairs, false)

@inline vtype(::VectorSpace{S}) where {S} = S

@inline isdual(V::VectorSpace) = V.isdual
@inline isdual(V1::T, V2::T) where {T<:VectorSpace} =
    V1.sectors == V2.sectors && V1.isdual != V2.isdual
@inline isdual(V1::TrivialVectorSpace, V2::TrivialVectorSpace) =
    isequal(V1, V2)

@inline isequal(V1::VectorSpace, V2::VectorSpace) =
    isequal(V1.sectors, V2.sectors) && V1.isdual == V2.isdual
@inline isequal(V1::TrivialVectorSpace, V2::TrivialVectorSpace) =
    V1.n == V2.n

==(V1::AbstractVectorSpace, V2::AbstractVectorSpace) = isequal(V1, V2)

@inline hascharge(V::VectorSpace{S}, charge::S) where{S}= haskey(V.sectors)
@inline findcharge(V::VectorSpace{S}, charge::S) where{S}= haskey(V.sectors)
@inline charges(V::VectorSpace) = [c for (c,d) in V.sectors]
@inline dims(V::VectorSpace) = [d for (c,d) in V.sectors]

"Find the dimension of the the given `charge`. Returns 0 if charges
doesn't exist. Return total dimension if charge is not specified."
@inline dim(V::VectorSpace{S}, charge::S) where{S}= get(V.sectors, charge, 0)
@inline dim(V::VectorSpace{S}, charge::Int) where{S}= get(V.sectors, convert(S, charge), 0)
@inline dim(V::VectorSpace) = sum(dims(V))
@inline dim(V::TrivialVectorSpace) = V.n

@inline Base.length(V::VectorSpace) = length(V.sectors)
@inline Base.iterate(V::VectorSpace) = iterate(V.sectors)
@inline Base.iterate(V::VectorSpace, i) = iterate(V.sectors, i)

"Return the dual of the vector space"
@inline dual(V::VectorSpace) = VectorSpace{vtype(V)}(V.sectors, !V.isdual)
@inline dual(V::AbstractVectorSpace) = V # is this correct!

function layout(V::VectorSpace)
    S = vtype(V)
    dict = SortedDict{S, UnitRange{Int}}()
    # if V.isdual
    #     p = dim(V)
    #     for (c, d) in V
    #         dict[c] = p-d+1:p
    #         p -= d
    #     end
    # else
        p = 0
        for (c,d) in V
            dict[c] = p+1:p+d
            p+=d
        end
    #end
    dict
end

const U1Space = VectorSpace{U1Charge}
U1Space(pairs::Pair{Int, Int}...) = U1Space(U1(c)=>d for (c,d) in pairs)

"Find sector intersections between two VSpace"
function intersect(V1::VectorSpace{S}, V2::VectorSpace{S}) where {S}
    V1.isdual != V2.isdual || error("Not the same dual-type")
    sectors = SortedDict{S, Int}()
    for (c, d) in V1.sectors
        if d == get(V2, c, 0)
            sectors[c] = d
        elseif d > 0
            error("same charge but different dims!")
        end
    end
    return VectorSpace(sectors)
end

"""
    fuse(Vs)

fuse any number of VSpaces into the tensor product VSpace. useful for
fuselegs fn. Note that fuse is a binary operation. That is the tensor
product is the bifunctor in the monoidal category of vector spaces. So
it should be by definition a recursive function. And the fusion tree
(order/style of recursion) has to specified. Here I assume a right to
left fusion, (V1 ⊗ (V2 ⊗ (V3 ⊗ ...))).

The order matters because the associator is not identity for
representation vector spaces as the catogory of tensor product of
representation spaces is not an strict monoidal cateogry. So one needs
the associator isomorphism to change the fusion tree.

"""

fuse(_dual::Bool, Vs::VectorSpace...) = fuse(_dual, Vs)
#fuse(Vs::VectorSpace...) = fuse(Vs)
#fuse(Vs::NTuple{VectorSpace{S}}) where{N,S} = fuse(false, Vs)
function fuse(_dual::Bool, Vs::NTuple{N, VectorSpace{S}}) where {N, S}
    if N < 2
        if _dual == Vs[1].isdual
            return Vs[1]
        else
            return dual(mapcharges(x->inv(x), Vs[1]))
        end
    elseif N > 2
        # lets for now generally always fuse from right to left!
        return fuse(_dual, Vs[1:N-2]..., fuse(_dual, Vs[N-1:N]...))
    end
    # N == 2 case
    sectors = SortedDict{S, Int}()
    for (c2, d2) in Vs[2]
        c2 = Vs[2].isdual ? inv(c2) : c2
        for (c1, d1) in Vs[1]
            c1 = Vs[1].isdual ? inv(c1) : c1
            c = _dual ? inv(c1+c2) : c1+c2
            if haskey(sectors, c)
                sectors[c] += d1*d2
            else
                sectors[c] = d1*d2
            end
        end
    end
    VectorSpace{S}(sectors, _dual)
end

"map the charges of the VSpace by some strictly ascending function `f`"
function mapcharges(f::Function, V::VectorSpace)
    mappedchrs = f.(charges(V))
    # issorted(mappedchrs) ||
    #     error("function not strictly ascending!")
    allunique(mappedchrs) ||
        error("mapping charges creates duplicates! $(l.chrs) -> $mappedchrs")
    VectorSpace{vtype(V)}(SortedDict(f(c)=>d for (c,d) in V.sectors), V.isdual)
end

isdummy(V::VectorSpace) = charges(V) == (zero(S),) && V[zero(S)] == 1
dummyleg() = VectorSpace{S}(0 => 1)

Base.show(io::IO, ::Type{U1Space}) = print(io, "U1Space")

function Base.show(io::IO, V::VectorSpace)
    show(io, typeof(V))
    V.isdual && print(io, "†")
    print(io, "(")
    seperator = ""
    comma = ", "
    io2 = IOContext(io, :typeinfo => typeof(V))
    for c in charges(V)
        if V.isdual
            print(io2, seperator, c, "'", "=>", dim(V, c))
        else
            print(io2, seperator, c, "=>", dim(V, c))
        end
        seperator = comma
    end
    print(io, ")")
    return nothing
end

ext: jl

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
@testset "Regression" begin
    @test has_fit_methods(LinearRegression)

    # Test for 10 random points
    for _ in 1:10
        X = rand(100, 5)
        y = rand(100)
        model = LinearRegression(X, y)
        fit!(model)
        β = X \ y
        for i in 1:size(X, 2)
            @test get_constrained_value(model, "β_$i") .≈ β[i] atol = 1e-5
        end
    end
    X = [
        1 1
        1 2
        1 3
        1 4
        1 5.0
    ]
    y = [
        6
        7
        8
        9
        10.0
    ]
    model = LinearRegression(X, y)
    fit!(model)
    forec = forecast(model, [1 6.0])
    forec.expected_value[1] == [11.0]
end


text:
using ITensors,
      Test

@testset "Tag Types" begin

  N = 10

  @testset "Star in operator strings" begin
    sites = siteinds("S=1/2",N)
    @test_throws ArgumentError op(sites, "Sp", 1)
    @test sites[1] isa Index
    Sz = op(sites,"Sz",2)
    SzSz = op(sites,"Sz * Sz",2)
    @test SzSz ≈ product(Sz, Sz)
    Sy = op(sites,"Sy",2)
    SySy = op(sites,"Sy * Sy",2)
    @test SySy ≈ product(Sy, Sy)

    sites = siteinds("S=1",N)
    @test_throws ArgumentError op(sites, "Sp", 1)
    Sz = op(sites,"Sz",2)
    SzSz = op(sites,"Sz * Sz",2)
    @test SzSz ≈ product(Sz, Sz)
    Sy = op(sites,"Sy",2)
    SySy = op(sites,"Sy * Sy",2)
    @test SySy ≈ product(Sy, Sy)
  end

  @testset "Custom TagType" begin

    # Use "_Custom_" tag even though this example
    # is for S=3/2, because we might define the
    # "S=3/2" TagType inside ITensors.jl later
    function ITensors.op(::TagType"_Custom_",
                         s::Index,
                         opname::AbstractString)
      Op = ITensor(s',dag(s))
      if opname == "Sz"
        Op[s'=>1,s=>1] = +3/2
        Op[s'=>2,s=>2] = +1/2
        Op[s'=>3,s=>3] = -1/2
        Op[s'=>4,s=>4] = -3/2
      end
      return Op
    end

    s = Index(4,"_Custom_")
    Sz = op(s,"Sz")
    @test Sz[s'=>1,s=>1] ≈ +3/2
    @test Sz[s'=>2,s=>2] ≈ +1/2
    @test Sz[s'=>3,s=>3] ≈ -1/2
    @test Sz[s'=>4,s=>4] ≈ -3/2
  end

end

nothing


text:
import Pkg; Pkg.add("InfrastructureSystems")

push!(LOAD_PATH, joinpath("..","src"))
push!(LOAD_PATH, joinpath("..","src","parse"))
push!(LOAD_PATH, joinpath("..","dev"))
push!(LOAD_PATH, joinpath("..","dev","buildstream"))

using DelimitedFiles
using Documenter

# Include SLiDE modules.
using SLiDE

# Now, generate the documentation.
makedocs(clean = true,
    modules = [SLiDE],
    format = Documenter.HTML(prettyurls = false),
    sitename="SLiDE",
    authors="Maxwell Brown, Caroline Hughes",
    pages=[
        "Home" => "index.md",
        "Data" => Any[
            "blueNOTE Data Set" => "lib/data.md",
            # "Data Stream" => "lib/datastream.md",
            "Build Stream" => "lib/buildstream.md",
            "Scaling" => "lib/scaling.md",
            "Parameters" => "lib/parameters.md"
        ],
        "Model" => "api/model.md",
        "Functions" => "SLiDE.md"
    ]
)

# pages = [
#         "Home" => "index.md",
#         "Manual" => Any[
#             "Guide" => "man/guide.md",
#             "man/examples.md",
#             "man/syntax.md",
#             "man/doctests.md",
#             "man/latex.md",
#             hide("man/hosting.md", [
#                 "man/hosting/walkthrough.md"
#             ]),
#             "man/other-formats.md",
#         ],

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
import CLBLAS
import OpenCL

const cl = OpenCL
const clblas = CLBLAS

for platform in cl.platforms()
    if platform[:name] == "Portable Computing Language"
        warn("Portable Computing Language platform not yet supported")
        continue
    end

    for device in cl.devices(platform)
        @printf("====================================================\n")
        @printf("Platform name:    %s\n",  platform[:name])
        @printf("Platform profile: %s\n",  platform[:profile])
        @printf("Platform vendor:  %s\n",  platform[:vendor])
        @printf("Platform version: %s\n",  platform[:version])
        @printf("----------------------------------------------------\n")
        @printf("Device name: %s\n", device[:name])
        @printf("Device type: %s\n", device[:device_type])
        @printf("Device mem: %i MB\n",           device[:global_mem_size] / 1024^2)
        @printf("Device max mem alloc: %i MB\n", device[:max_mem_alloc_size] / 1024^2)
        @printf("Device max clock freq: %i MHZ\n",  device[:max_clock_frequency])
        @printf("Device max compute units: %i\n",   device[:max_compute_units])
        @printf("Device max work group size: %i\n", device[:max_work_group_size])
        @printf("Device max work item size: %s\n",  device[:max_work_item_size])
        @printf("Device max work dimension size: %s\n",  device[:max_work_item_dims])
   end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@doc doc"""
    abstract type AbstractPosteriorDensity <: AbstractDensity end

Abstract super-type for posterior probability densities.
"""
abstract type AbstractPosteriorDensity <: AbstractDensity end
export AbstractPosteriorDensity

@doc doc"""
    getlikelihood(posterior::AbstractPosteriorDensity)::AbstractDensity

*BAT-internal, not part of stable public API.*

The likelihood density of `posterior`. The likelihood may or may not be
normalized.
"""
function getlikelihood end

@doc doc"""
    getprior(posterior::AbstractPosteriorDensity)::AbstractDensity

*BAT-internal, not part of stable public API.*

The prior density of `posterior`. The prior may or may not be normalized.
"""
function getprior end

@doc doc"""
    BAT.apply_bounds_and_eval_posterior_logval!(
        T::Type{<:Real},
        density::AbstractDensity,
        params::AbstractVector{<:Real}
    )

*BAT-internal, not part of stable public API.*

First apply bounds to the parameters, then compute and return posterior log
value. May modify `params`.

Guarantees that  :

* If parameters are still out of bounds after applying bounds,
  `density_logval` is not called for either prior or likelihood.
* If `density_logval` for prior returns `-Inf`, `density_logval` is not called
  for likelihood.

In both cases, `T(-Inf)` is returned for both prior and posterior.
"""
function apply_bounds_and_eval_posterior_logval! end

function apply_bounds_and_eval_posterior_logval!(
    T::Type{<:Real},
    posterior::AbstractPosteriorDensity,
    params::AbstractVector{<:Real}
)
    bounds = param_bounds(posterior)
    apply_bounds!(params, bounds)

    parshapes = params_shape(posterior)
    zero_prob_logval = convert(T, -Inf)

    prior_logval = if !isoob(params)
        convert(T, eval_density_logval(getprior(posterior), params, parshapes))
    else
        zero_prob_logval
    end

    likelihood_logval = if prior_logval > zero_prob_logval
        convert(T, eval_density_logval(getlikelihood(posterior), params, parshapes))
    else
        zero_prob_logval
    end

    convert(T, prior_logval + likelihood_logval)
end

@doc doc"""
    BAT.apply_bounds_and_eval_posterior_logval_strict!(
        posterior::AbstractPosteriorDensity,
        params::AbstractVector{<:Real}
    )

*BAT-internal, not part of stable public API.*

First apply bounds to the parameters, then compute and return posterior log
value. May modify `params`.

Guarantees that  :

* If parameters are still out of bounds after applying bounds,
  `density_logval` is not called for either prior or likelihood.
* If `density_logval` for prior returns `-Inf`, `density_logval` is not called
  for likelihood.

In both cases, an exception is thrown.
"""
function apply_bounds_and_eval_posterior_logval_strict!(
    posterior::AbstractPosteriorDensity,
    params::AbstractVector{<:Real}
)
    bounds = param_bounds(posterior)
    apply_bounds!(params, bounds)

    parshapes = params_shape(posterior)

    prior_logval = if !isoob(params)
        eval_density_logval(getprior(posterior), params, parshapes)
    else
        throw(ArgumentError("Parameter(s) out of bounds"))
    end

    likelihood_logval = if prior_logval > convert(typeof(prior_logval), -Inf)
        eval_density_logval(getlikelihood(posterior), params, parshapes)
    else
        throw(ErrorException("Prior density must not be zero."))
    end

    T = typeof(likelihood_logval)

    convert(T, convert(T, prior_logval) + likelihood_logval)
end

function density_logval(density::AbstractPosteriorDensity, params::AbstractVector{<:Real})
    parshapes = params_shape(density)
    eval_density_logval(getprior(density), params, parshapes) +
    eval_density_logval(getlikelihood(density), params, parshapes)
end

function density_logval(density::AbstractPosteriorDensity, params::Any)
    density_logval(getprior(density), params) +
    density_logval(getlikelihood(density), params)
end

function param_bounds(density::AbstractPosteriorDensity)
    li_bounds = param_bounds(getlikelihood(density))
    pr_bounds = param_bounds(getprior(density))
    if ismissing(li_bounds)
        pr_bounds
    else
        li_bounds ∩ pr_bounds
    end
end

function nparams(density::AbstractPosteriorDensity)
    li_np = nparams(getlikelihood(density))
    pr_np = nparams(getprior(density))
    ismissing(li_np) || li_np == pr_np || error("Likelihood and prior have different number of parameters")
    pr_np
end

@doc doc"""
    PosteriorDensity{
        Li<:AbstractDensity,
        Pr<:DistLikeDensity
    } <: AbstractPosteriorDensity

A representation of a PosteriorDensity, based a likelihood and prior, a
representation of the posterior density is cached internally. The densities
be accessed via

```julia
getlikelihood(posterior::PosteriorDensity)::Li
getprior(posterior::PosteriorDensity)::Pr
```

Constructors:

```julia
PosteriorDensity(likelihood::AbstractDensity, prior::AbstractDensity)
PosteriorDensity(likelihood::Any, prior::Any)
PosteriorDensity(log_likelihood::Function, prior::Any)
```
"""
struct PosteriorDensity{
    L<:AbstractDensity,
    P<:DistLikeDensity,
    B<:AbstractParamBounds,
    S<:AbstractValueShape
} <: AbstractPosteriorDensity
    likelihood::L
    prior::P
    parbounds::B
    parshapes::S
end

export PosteriorDensity

function PosteriorDensity(likelihood::Any, prior::Any)
    li = convert(AbstractDensity, likelihood)
    pr = convert(AbstractDensity, prior)

    parbounds = _posterior_parbounds(
        param_bounds(li),
        param_bounds(pr)
    )

    parshapes = _posterior_parshapes(
        params_shape(li),
        params_shape(pr)
    )

    PosteriorDensity(li, pr, parbounds, parshapes)
end

getlikelihood(posterior::PosteriorDensity) = posterior.likelihood

getprior(posterior::PosteriorDensity) = posterior.prior

param_bounds(posterior::PosteriorDensity) = posterior.parbounds

params_shape(posterior::PosteriorDensity) = posterior.parshapes

function _posterior_parshapes(li_ps::AbstractValueShape, pr_ps::AbstractValueShape)
    if li_ps == pr_ps
        li_ps
    else
        throw(ArgumentError("Variable shapes of likelihood and prior are incompatible"))
    end
end

function _posterior_parshapes(li_ps::AbstractValueShape, pr_ps::ArrayShape{T,1}) where T
    n = totalndof(li_ps)
    if n == totalndof(pr_ps)
        @assert size(pr_ps) == (n,)
        li_ps
    else
        throw(ArgumentError("Likelihood and prior have different number of free parameters"))
    end
end

_posterior_parshapes(li_ps::ArrayShape{T,1}, pr_ps::AbstractValueShape) where T =
    _posterior_parshapes(pr_ps, li_ps)

function _posterior_parshapes(li_ps::ArrayShape{T,1}, pr_ps::ArrayShape{U,1}) where {T,U}
    if size(li_ps) == size(pr_ps)
        li_ps
    else
        throw(ArgumentError("Likelihood and prior have different number of free parameters"))
    end
end

_posterior_parshapes(li_ps::Missing, pr_ps::AbstractValueShape) = pr_ps

_posterior_parshapes(li_ps::Union{AbstractValueShape,Missing}, pr_ps::Missing) =
    throw(ArgumentError("Parameter shapes of prior must not be missing"))

_posterior_parbounds(li_bounds::AbstractParamBounds, pr_bounds::AbstractParamBounds) =
     li_bounds ∩ pr_bounds

_posterior_parbounds(li_bounds::Missing, pr_bounds::AbstractParamBounds) = pr_bounds


text:
module GLMBasis

using LinearAlgebra

export make_stim_basis, make_spike_basis, spikefilter
# ============================================================================ #
@inline stretchx(x::Real) = log(x + 1e-20)
@inline stretchx_inv(x::Real) = exp(x) - 1e-20
@inline pimin(x::Real) = min(pi, x)
@inline pimax(x::Real) = max(-pi, x)
# ============================================================================ #
function orth(x::AbstractMatrix{<:Real})
    return try
        f = svd(x)
        k = rank(x)
        return f.U[:,1:k]
    catch err
        return fill(NaN, size(x))
    end
end
# ============================================================================ #
function get_indicies(s1::Integer, s2::Integer)
    pad = s2 - s1
    isodd(pad) && error("valid portion of result is not centered")
    n = floor(Int, pad/2)
    return n+1:s2-n
end
# ============================================================================ #
function spikefilter(idx::Vector{<:Integer}, h::AbstractArray{<:Real}, twin::Tuple{Int,Int})

    nrow, ncol = size(h)
    t1, t2 = twin
    t0 = t1 - nrow + 1
    len = (t2 - t1) + 1

    out = zeros(len, ncol)

    ks = searchsortedfirst(idx, t0)

    if ks < length(idx)

        @inbounds while (ks < length(idx)) && (idx[ks] < t2+1)
            isp = idx[ks]
            i1 = isp < t1 ? t1 : isp
            imx = isp + nrow
            imx = imx > t2 ? t2 : imx

            @inbounds for k = 1:ncol
                ycol = (k-1) * len - t1 + 1
                hcol = (k-1) * nrow + 1

                @inbounds for j = i1:(imx-1)
                    out[ycol + j] += h[hcol+j-isp]
                end

            end

            ks += 1
        end
    end

    return out
end
# ============================================================================ #
function make_basis(ncol::Integer, kpeak::AbstractVector{<:Real}, b::Real, dt::Real=0.0)
    yrange = stretchx.(kpeak .+ b)

    # it appears the original algorithm fails in some circumstances due to
    # round off error (I think...), so just use range() to guarantee the correct
    # size
    centers = range(yrange[1], yrange[2], length=ncol)
    spacing = step(centers)

    max_time = stretchx_inv(yrange[2] + 2 * spacing) - b

    if dt > 0.0
        t0 = 0.0:dt:max_time
    else
        t0 = 0:1:max_time
    end

    #this is actually faster and more efficient than the correspoding loop, not
    #sure how the compiler is pulling that off... well done Julia
    x = (repeat(stretchx.(t0 .+ b), 1, ncol) .- repeat(centers', length(t0), 1)) .* (pi / spacing / 2)

    return (cos.(pimax.(pimin.(x))) .+ 1.0) ./ 2.0, t0
end
# ============================================================================ #
function make_spike_basis(ncol::Integer, kpeak::AbstractVector{<:Real}, b::Real, dt::Real, refractory::Real=0, nt::Integer=0)

    if refractory >= dt
        ncol -= 1
    elseif refractory > 0
        @warn("Refractory period is smaller than bin size!")
    end

    basis, t0 = make_basis(ncol, kpeak, b, dt)

    if refractory >= dt
        k = findlast(x->x<=refractory, t0)
        basis[1:k, :] .= 0.0
        basis = hcat(zeros(size(basis, 1)), basis)
        basis[1:k, 1] .= 1
    end

    if nt > 0
        nrow = size(basis, 1)
        if nt > nrow
            basis = vcat(basis, zeros(nt-nrow, size(basis, 2)))
        elseif nt < nrow
            basis = basis[1:nt, :]
        end
    end

    return basis, orth(basis)
end
# ============================================================================ #
function make_stim_basis(neye::Integer, ncos::Integer, kpeak::AbstractVector{<:Integer}, b::Real, nt::Integer=0)

    basis, t0 = make_basis(ncos, kpeak, b)

    # reverse in time so that fine-scale is at the end
    basis2 = reverse(
        hcat(
            vcat(Matrix{eltype(basis)}(I, neye, neye), zeros(length(t0), neye)),
            vcat(zeros(neye, ncos), basis)
        ),
        dims=1
    )

    if nt > 0
        nrow = size(basis2, 1)
        if nt > nrow
            basis2 = vcat(zeros(nt-nrow, ncos+neye), basis2)
        elseif nt < nrow
            basis2 = basis2[end-nt+1:end, :]
        end
    end

    for k in 1:size(basis, 2)
        basis2[:,k] ./= sqrt(sum(abs2, basis[:,k]))
    end

    return basis2, orth(basis2)
end
# ============================================================================ #
end


text:
import Base: eye, rand, *, inv

# Type definitions

immutable OrthogonalGroup{T, N} <: Group
    α::T
    β::T
    γ::T
end

O3{T} = OrthogonalGroup{T, 3}
O3{T}(α::T, β::T, γ::T) = OrthogonalGroup{T, 3}(α, β, γ)

# Identity element

eye(::Type{O3}) = O3(0.0, 0.0, 0.0)

# Random group element

function rand(::Type{O3})
    α, β, γ = 2π*rand(), π*rand(), 2π*rand()
    O3(α, β, γ)
end

# Defining representation

@inbounds function defrep!{T}(M::Matrix{T}, g::O3{T})
    α, β, γ = g.α, g.β, g.γ
    M[1, 1] = cos(γ)⋅cos(α)⋅cos(β) - sin(γ)⋅sin(α)
    M[2, 1] = sin(α)⋅cos(γ)⋅cos(β) + sin(γ)⋅cos(α)
    M[3, 1] = -sin(β)⋅cos(γ)
    M[1, 2] = -sin(γ)⋅cos(α)⋅cos(β) - cos(γ)⋅sin(α)
    M[2, 2] = -sin(γ)⋅sin(α)⋅cos(β) + cos(γ)⋅cos(α)
    M[3, 2] = sin(β)⋅sin(γ)
    M[1, 3] = sin(β)⋅cos(α)
    M[2, 3] = sin(β)⋅sin(α)
    M[3, 3] = cos(β)

    return M
end

defrep{T}(g::O3{T}) = defrep!(Matrix{T}(3, 3), g)

# Hacky and suboptimal. Optimize this.
# Procedure from
#   http://mathworld.wolfram.com/EulerAngles.html
#   http://mathworld.wolfram.com/NonlinearLeastSquaresFitting.html

function nearest{T}(::Type{O3{T}}, M_defrep::Matrix{T}, tol = T(1e-8))
    b = vec(M_defrep)
    tol_sq = tol^2

    α, β, γ = 2π*rand(T), π*rand(T), 2π*rand(T)
    dα, dβ, dγ = T(Inf), T(Inf), T(Inf)

    while dα^2 + dβ^2 + dγ^2 > tol_sq
        A = [(-cos(γ)*sin(α)*cos(β) - sin(γ)*cos(α))    (-cos(γ)*cos(α)*sin(β))     (-sin(γ)*cos(α)*cos(β) - cos(γ)*sin(α))
             (cos(α)*cos(γ)*cos(β) - sin(γ)*sin(α))     (-sin(α)*cos(γ)*sin(β))     (-sin(α)*sin(γ)*cos(β) + cos(γ)*cos(α))
             (zero(T))                                  (-cos(β)*cos(γ))            (sin(β)*sin(γ))
             (sin(γ)*sin(α)*cos(β) - cos(γ)*cos(α))     (sin(γ)*cos(α)*sin(β))      (-cos(γ)*cos(α)*cos(β) + sin(γ)*sin(α))
             (-sin(γ)*cos(α)*cos(β) - cos(γ)*sin(α))    (sin(γ)*sin(α)*sin(β))      (-cos(γ)*sin(α)*cos(β) - sin(γ)*cos(α))
             (zero(T))                                  (cos(β)*sin(γ))             (sin(β)*cos(γ))
             (-sin(β)*sin(α))                           (cos(β)*cos(α))             (zero(T))
             (sin(β)*cos(α))                            (cos(β)*sin(α))             (zero(T))
             (zero(T))                                  (-sin(β))                   (zero(T))]

        dα, dβ, dγ = A \ b
        α = mod2pi(α + dα)
        β = mod2pi(2β + 2dβ) / 2
        γ = mod2pi(γ + dγ)
    end

    return O3(α, β, γ)
end

nearest{T}(::Type{O3}, M::Matrix{T}) = nearest(O3{T}, M)

# Multiplication
# Hacky.  Find an analytic expression for multiplying Euler-ZYZ angles
function *{T}(g₁::O3{T}, g₂::O3{T})
    nearest(O3, defrep(g₁) * defrep(g₂))
end

# Inversion
# Hacky. Find an analytic expression for inverting Euler-ZYZ angles
function inv{T}(g::O3{T})
    nearest(O3, inv(defrep(g)))
end


text:
using NQCDynamics
using Test
using OrdinaryDiffEq

@testset "select_reduction" begin
    @test Ensembles.select_reduction(:mean) isa Ensembles.MeanReduction
    @test Ensembles.select_reduction(:sum) isa Ensembles.SumReduction
    @test Ensembles.select_reduction(:append) isa Function
end

@testset "get_u_init" begin

    @testset "AbstractOutput, reduction=$reduction" for reduction in (:mean, :sum)
        reduction = Ensembles.select_reduction(reduction)
        tspan = (0.0, 1.0)
        u0 = 1.0
        saveat = 1
        output = Ensembles.OutputFinal()
        u_init = Ensembles.get_u_init(reduction, saveat, Dict(), tspan, u0, output)
        @test u_init ≈ zero(u0)
    end

    @testset "reduction=:append" begin
        reduction = Ensembles.select_reduction(:append)
        output = Ensembles.OutputFinal()
        u0 = 1
        saveat = 1
        u_init = Ensembles.get_u_init(reduction, saveat, Dict(), (0, 1), u0, output)
        @test u_init == []
    end

    @testset "TimeCorrelationFunction, reduction=$reduction" for reduction in (:mean, :sum)
        reduction = Ensembles.select_reduction(reduction)
        sim = Simulation{FSSH}(Atoms(1), DoubleWell())
        saveat = 0.1
        tspan = (0.0, 1.0)
        u0 = DynamicsVariables(sim, 1, 0, SingleState(1, Adiabatic()))
        output = TimeCorrelationFunctions.PopulationCorrelationFunction(sim, Diabatic())
        u_init = Ensembles.get_u_init(reduction, saveat, Dict(), tspan, u0, output)
        @test u_init ≈ [zeros(2, 2) for _ ∈ tspan[1]:0.1:tspan[2]]
    end

end

@testset "$reduction" for reduction in (Ensembles.MeanReduction(0), Ensembles.SumReduction())
    prob = ODEProblem((u,p,t) -> 1.01u, 0.5, (0.0,1.0))
    prob_func(prob, i, repeat) = remake(prob,u0=rand())
    output = Ensembles.OutputFinal()
    u_init = Ensembles.output_template(output, rand())
    ensemble = EnsembleProblem(prob, prob_func=prob_func, output_func=output, reduction=reduction, u_init=u_init)
    sol = solve(ensemble, Tsit5(), trajectories=1e3, batch_size=20)
    if reduction isa Ensembles.SumReduction
        @test sol.u / 1e3 ≈ 0.5 * exp(1.01) rtol=1e-1
    else
        @test sol.u ≈ 0.5 * exp(1.01) rtol=1e-1
    end
end


text:
using Latte
using FactCheck

facts("Testing simple RNN layer") do
    n_time_steps = 5
    net = Net(2; time_steps=n_time_steps)
    data, data_value = MemoryDataLayer(net, :data, (8,))
    rand!(data_value)

    rnn = RNNBlock(:rnn, net, data)
    init(net)

    context("Forward") do
        h_weights = get_buffer(net, :rnnstateweights)
        h_bias = get_buffer(net, :rnnstatebias)

        x_weights = get_buffer(net, :rnnxweights)
        x_bias = get_buffer(net, :rnnxbias)

        weights = get_buffer(net, :rnnweights)
        bias = get_buffer(net, :rnnbias)

        forward(net; phase=Latte.Test)
        h_prev = zeros(Float32, size(get_buffer(net, :rnntanhvalue))...)
        for t in 1:n_time_steps
            @fact get_buffer(net, :datavalue, t) --> data_value[:,:,t]

            x_expected = x_weights' * data_value[:, :, t] .+ reshape(x_bias, prod(size(x_bias)))
            @fact get_buffer(net, :rnnxvalue, t) --> roughly(x_expected)

            state_expected = h_weights' * h_prev .+ reshape(h_bias, prod(size(h_bias)))
            @fact get_buffer(net, :rnnstatevalue, t) --> roughly(state_expected)

            # h_expected = state_expected .+ x_expected
            h_expected = tanh(state_expected .+ x_expected)
            @fact get_buffer(net, :rnntanhvalue, t) --> roughly(h_expected)
            h_prev = h_expected

            output_expected = weights' * h_expected .+ reshape(bias, prod(size(bias)))
            @fact get_buffer(net, :rnnvalue, t) --> roughly(output_expected)
        end
    end
end

FactCheck.exitstatus()


text:
function calDIV!(;
    ASUM      :: AdvectionSpeedUpMatrix,
    u_U       :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny   )
    v_V       :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny+1 )
    div       :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny   )
    workspace :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny   )
)

    mul_autoflat!(div,       ASUM.T_DIVx_U, u_U)
    mul_autoflat!(workspace, ASUM.T_DIVy_V, v_V)
    @. div += workspace

end

function calDiffAdv_QUICK_SpeedUp!(
    ocn         :: Ocean;
    qs          :: AbstractArray{Float64, 3},
    dΔqdt       :: AbstractArray{Float64, 2},     # ( Nx, Ny )
    FLUX_CONV   :: AbstractArray{Float64, 3},
    FLUX_CONV_h :: AbstractArray{Float64, 3},
    FLUX_DEN_x  :: AbstractArray{Float64, 3},     # ( Nz_bone   ,  Nx+1, Ny   )
    FLUX_DEN_y  :: AbstractArray{Float64, 3},     # ( Nz_bone   ,  Nx  , Ny+1 )
    FLUX_DEN_z  :: AbstractArray{Float64, 3},     # ( Nz_bone+1 ,  Nx  , Ny   )
    Dh          :: Float64,
    Dv          :: Float64,
    Δt          :: Float64,
)

    wksp = ocn.wksp
    ASUM = ocn.ASUM

    tmp1 = getSpace!(wksp, :T)
    tmp2 = getSpace!(wksp, :T)

   #println("GRAD_CRUV")

    mul_autoflat!(ocn.GRAD_bnd_x, ASUM.U_∂x_T, qs)
    mul_autoflat!(ocn.GRAD_bnd_y, ASUM.V_∂y_T, qs)
    mul_autoflat!(ocn.GRAD_bnd_z, ASUM.W_∂z_T, qs)

    mul_autoflat!(ocn.CURV_x, ASUM.T_∂x_U, ocn.GRAD_bnd_x)
    mul_autoflat!(ocn.CURV_y, ASUM.T_∂y_V, ocn.GRAD_bnd_y)
    mul_autoflat!(ocn.CURV_z, ASUM.T_∂z_W, ocn.GRAD_bnd_z)

    calFluxDensity_abstract!(;
        X_T                 = qs,
        u_bnd               = ocn.u_bnd,
        GRAD_bnd            = ocn.GRAD_bnd_x,
        CURV_T              = ocn.CURV_x,
        FLUX_DEN_bnd        = FLUX_DEN_x,
        K                   = Dh,
        Δt                  = Δt,
        Δx_bnd              = ASUM.Δx_U,

        op_bnd_∂x_T         = ASUM.U_∂x_T,
        op_T_∂x_bnd         = ASUM.T_∂x_U,
        op_bnd_pos_dir_T    = ASUM.op.U_E_T,
        op_bnd_neg_dir_T    = ASUM.op.U_W_T,
        op_bnd_interp_T     = ASUM.U_interp_T,
        op_filter_bnd       = ASUM.U_fluxmask_U,

        pos_u_mask_bnd      = getSpace!(wksp, :U),
        CURV_on_pos_dir_bnd = getSpace!(wksp, :U),
        CURV_on_neg_dir_bnd = getSpace!(wksp, :U),
        CURV_r_bnd          = getSpace!(wksp, :U),
        uΔt_bnd             = getSpace!(wksp, :U),
        X_star_bnd          = getSpace!(wksp, :U),

    )

    calFluxDensity_abstract!(;
        X_T                 = qs,
        u_bnd               = ocn.v_bnd,
        GRAD_bnd            = ocn.GRAD_bnd_y,
        CURV_T              = ocn.CURV_y,
        FLUX_DEN_bnd        = FLUX_DEN_y,
        K                   = Dh,
        Δt                  = Δt,
        Δx_bnd              = ASUM.Δy_V,

        op_bnd_∂x_T         = ASUM.V_∂y_T,
        op_T_∂x_bnd         = ASUM.T_∂y_V,
        op_bnd_pos_dir_T    = ASUM.op.V_N_T,
        op_bnd_neg_dir_T    = ASUM.op.V_S_T,
        op_bnd_interp_T     = ASUM.V_interp_T,
        op_filter_bnd       = ASUM.V_fluxmask_V,

        pos_u_mask_bnd      = getSpace!(wksp, :V),
        CURV_on_pos_dir_bnd = getSpace!(wksp, :V),
        CURV_on_neg_dir_bnd = getSpace!(wksp, :V),
        CURV_r_bnd          = getSpace!(wksp, :V),
        uΔt_bnd             = getSpace!(wksp, :V),
        X_star_bnd          = getSpace!(wksp, :V),

    )

    calFluxDensity_abstract!(;
        X_T                 = qs,
        u_bnd               = ocn.w_bnd,
        GRAD_bnd            = ocn.GRAD_bnd_z,
        CURV_T              = ocn.CURV_z,
        FLUX_DEN_bnd        = FLUX_DEN_z,
        K                   = Dv,
        Δt                  = Δt,
        Δx_bnd              = ASUM.Δz_W,

        op_bnd_∂x_T         = ASUM.W_∂z_T,
        op_T_∂x_bnd         = ASUM.T_∂z_W,
        op_bnd_pos_dir_T    = ASUM.op.W_UP_T,
        op_bnd_neg_dir_T    = ASUM.op.W_DN_T,
        op_bnd_interp_T     = ASUM.W_interp_T,
        op_filter_bnd       = ASUM.W_fluxmask_W,

        pos_u_mask_bnd      = getSpace!(wksp, :W),
        CURV_on_pos_dir_bnd = getSpace!(wksp, :W),
        CURV_on_neg_dir_bnd = getSpace!(wksp, :W),
        CURV_r_bnd          = getSpace!(wksp, :W),
        uΔt_bnd             = getSpace!(wksp, :W),
        X_star_bnd          = getSpace!(wksp, :W),

    )

    mul_autoflat!(FLUX_CONV_h, ASUM.T_DIVx_U, FLUX_DEN_x)
    mul_autoflat!(tmp1,        ASUM.T_DIVy_V, FLUX_DEN_y)
    mul_autoflat!(tmp2,        ASUM.T_DIVz_W, FLUX_DEN_z)

    @. FLUX_CONV_h = -1.0 * (FLUX_CONV_h + tmp1)
    @. FLUX_CONV = FLUX_CONV_h - tmp2

   #println("CALMIXEDLAYER")
    calMixedLayer_dΔqdt!(
        Nx          = ocn.Nx,
        Ny          = ocn.Ny,
        Nz          = ocn.Nz,
        FLUX_CONV_h = FLUX_CONV_h,
        FLUX_DEN_z  = FLUX_DEN_z,
        dΔqdt       = dΔqdt,
        mask        = ocn.mask,
        FLDO        = ocn.FLDO,
        h_ML        = ocn.h_ML,
        hs          = ocn.hs,
        zs          = ocn.zs,
    )
end

function calFluxDensity_abstract!(;
    X_T                       :: AbstractArray{Float64, 3},
    u_bnd                     :: AbstractArray{Float64, 3},
    GRAD_bnd                  :: AbstractArray{Float64, 3},
    CURV_T                    :: AbstractArray{Float64, 3},
    FLUX_DEN_bnd              :: AbstractArray{Float64, 3},
    K                         :: Float64,
    Δt                        :: Float64,
    Δx_bnd                    :: AbstractArray{Float64, 3},

    op_bnd_∂x_T               :: AbstractArray{Float64, 2},
    op_T_∂x_bnd               :: AbstractArray{Float64, 2},
    op_bnd_pos_dir_T          :: AbstractArray{Float64, 2},
    op_bnd_neg_dir_T          :: AbstractArray{Float64, 2},
    op_bnd_interp_T           :: AbstractArray{Float64, 2},
    op_filter_bnd             :: AbstractArray{Float64, 2},

    pos_u_mask_bnd            :: AbstractArray{Float64, 3},
    CURV_on_pos_dir_bnd       :: AbstractArray{Float64, 3},
    CURV_on_neg_dir_bnd       :: AbstractArray{Float64, 3},
    CURV_r_bnd                :: AbstractArray{Float64, 3},
    uΔt_bnd                   :: AbstractArray{Float64, 3},
    X_star_bnd                :: AbstractArray{Float64, 3},
)

    mul_autoflat!(GRAD_bnd, op_bnd_∂x_T, X_T)
    mul_autoflat!(CURV_T, op_T_∂x_bnd, GRAD_bnd)

    @. pos_u_mask_bnd = u_bnd >= 0
    mul_autoflat!(CURV_on_pos_dir_bnd, op_bnd_neg_dir_T, CURV_T)
    mul_autoflat!(CURV_on_neg_dir_bnd, op_bnd_pos_dir_T, CURV_T)
    @. CURV_r_bnd = CURV_on_neg_dir_bnd * pos_u_mask_bnd + CURV_on_pos_dir_bnd * (1-pos_u_mask_bnd)

    @. uΔt_bnd    = u_bnd * Δt
    mul_autoflat!(X_star_bnd, op_bnd_interp_T, X_T)
    #verbose && println("interpolated: ", X_T[1,37,23])
    @. X_star_bnd = X_star_bnd - uΔt_bnd / 2.0 * GRAD_bnd + ( K * Δt / 2.0 - Δx_bnd^2.0/6.0 + uΔt_bnd^2.0 / 6.0 ) * CURV_r_bnd
    #verbose && println("X_star_bnd: ", X_T[1,37,23])

    # reuse
    empty = uΔt_bnd
    @. empty = u_bnd * X_star_bnd - K * ( GRAD_bnd - uΔt_bnd / 2.0 * CURV_r_bnd )

    #verbose && println("Flux: ", empty[1,37,23])

    mul_autoflat!(FLUX_DEN_bnd, op_filter_bnd, empty)

    #verbose && println("Flux after filter: ", FLUX_DEN_bnd[1,37,23])
    #verbose && println("speed bnd : ", u_bnd[1,37,23])

end

function calMixedLayer_dΔqdt!(;
    Nx          :: Integer,
    Ny          :: Integer,
    Nz          :: AbstractArray{Int64, 2},
    FLUX_CONV_h :: AbstractArray{Float64, 3},     # ( Nz_bone  ,  Nx, Ny )
    FLUX_DEN_z  :: AbstractArray{Float64, 3},     # ( Nz_bone+1,  Nx, Ny )
    dΔqdt       :: AbstractArray{Float64, 2},     # ( Nx, Ny )
    mask        :: AbstractArray{Float64, 2},     # ( Nx, Ny )
    FLDO        :: AbstractArray{Int64, 2},       # ( Nx, Ny )
    h_ML        :: AbstractArray{Float64, 2},     # ( Nx, Ny )
    hs          :: AbstractArray{Float64, 3},     # ( Nz_bone  ,  Nx, Ny )
    zs          :: AbstractArray{Float64, 3},     # ( Nz_bone+1,  Nx, Ny )
)

    for i=1:Nx, j=1:Ny

        if mask[i, j] == 0.0
            continue
        end

        _FLDO = FLDO[i, j]

        if _FLDO == -1
            continue
        end

        tmp = 0.0
        for k = 1:_FLDO-1
            tmp += FLUX_CONV_h[k, i, j] * hs[k, i, j]
        end
        tmp += (
              FLUX_CONV_h[_FLDO, i, j] * zs[_FLDO, i, j]
            + ( FLUX_DEN_z[_FLDO+1, i, j] * zs[_FLDO, i, j] - FLUX_DEN_z[_FLDO, i, j] * zs[_FLDO+1, i, j] ) / hs[_FLDO, i, j]
        )

        dΔqdt[i, j] = tmp / h_ML[i, j]
    end

end

function calDIV!(;
    gi       :: DisplacedPoleCoordinate.GridInfo,
    Nx       :: Integer,
    Ny       :: Integer,
    Nz       :: AbstractArray{Int64, 2},
    u_bnd    :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx+1, Ny   )
    v_bnd    :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny+1 )
    div      :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny   )
    mask3    :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny   )
)

#    local tmp = tmp_σ = 0.0
    for i=1:Nx, j=1:Ny

        for k=1:Nz[i, j]

            if mask3[k, i, j] == 0.0
                break
            end

            div[k, i, j] =  (
                u_bnd[k, i+1, j  ]  * gi.DY[i+1, j  ]
              - u_bnd[k, i,   j  ]  * gi.DY[i  , j  ]
              + v_bnd[k, i,   j+1]  * gi.DX[i  , j+1]
              - v_bnd[k, i,   j  ]  * gi.DX[i  , j  ]
            ) / gi.dσ[i, j]

        end

    end

end

function calVerVelBnd!(;
    gi       :: DisplacedPoleCoordinate.GridInfo,
    Nx       :: Integer,
    Ny       :: Integer,
    Nz       :: AbstractArray{Int64, 2},
    w_bnd    :: AbstractArray{Float64, 3},   # ( Nz_bone+1, Nx  , Ny   )
    hs       :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny   )
    div      :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny   )
    mask3    :: AbstractArray{Float64, 3},   # ( Nz_bone  , Nx  , Ny   )
)

#    local tmp = tmp_σ = 0.0
    for i=1:Nx, j=1:Ny

        w_bnd[1, i, j] = 0.0

        for k=1:Nz[i, j]

            if mask3[k, i, j] == 0.0
                break
            end

            w_bnd[k+1, i, j] = w_bnd[k, i, j] + div[k, i, j] * hs[k, i, j]
        end

#        tmp   += w_bnd[Nz[i, j]+1, i, j] * gi.dσ[i, j]
#        tmp_σ += gi.dσ[i, j]
    end

#   #println("tmp: ", tmp, "; tmp_σ: ", tmp_σ, "; Average w: ", tmp/tmp_σ)

end

function calHorVelBnd!(;
    Nx       :: Integer,
    Ny       :: Integer,
    Nz       :: AbstractArray{Int64, 2},
    weight_e :: AbstractArray{Float64, 2},   # (Nx+1, Ny)
    weight_n :: AbstractArray{Float64, 2},   # (Nx, Ny+1)
    u        :: AbstractArray{Float64, 3},   # (Nz_bone, Nx, Ny)
    v        :: AbstractArray{Float64, 3},   # (Nz_bone, Nx, Ny)
    u_bnd    :: AbstractArray{Float64, 3},   # (Nz_bone, Nx+1, Ny)
    v_bnd    :: AbstractArray{Float64, 3},   # (Nz_bone, Nx, Ny+1)
    mask3    :: AbstractArray{Float64, 3},   # (Nz_bone, Nx, Ny)
    noflux_x_mask3 :: AbstractArray{Float64, 3}, # (Nz_bone, Nx+1, Ny)
    noflux_y_mask3 :: AbstractArray{Float64, 3}, # (Nz_bone, Nx, Ny+1)
)

    # x
    for i=2:Nx, j=1:Ny
        for k=1:Nz[i, j]
            if noflux_x_mask3[k, i, j] == 0.0
                u_bnd[k, i, j] = 0.0
            else
                u_bnd[k, i, j] = u[k, i-1, j] * (1.0 - weight_e[i, j]) + u[k, i, j] * weight_e[i, j]
                #u_bnd[k, i, j] = (u[k, i-1, j] + u[k, i, j]) / 2.0
            end
        end
    end

    # x - periodic
    for j=1:Ny
        for k=1:Nz[1, j]
            if noflux_x_mask3[k, 1, j] == 0.0
                u_bnd[k, 1, j] = u_bnd[k, Nx+1, j] = 0.0
            else
                u_bnd[k, 1, j] = u_bnd[k, Nx+1, j] = u[k, Nx, j] * (1.0 - weight_e[1, j]) + u[k, 1, j] * weight_e[1, j]
                #u_bnd[k, 1, j] = u_bnd[k, Nx+1, j] = (u[k, Nx, j] + u[k, 1, j]) / 2.0
            end
        end
    end

    # y
    for i=1:Nx, j=2:Ny
        for k=1:Nz[i, j]
            if noflux_y_mask3[k, i, j] == 0.0
                v_bnd[k, i, j] = 0.0
            else
                v_bnd[k, i, j] = v[k, i, j-1] * (1.0 - weight_n[i, j]) + v[k, i, j] * weight_n[i, j]
                #v_bnd[k, i, j] = (v[k, i, j-1] + v[k, i, j]) / 2.0
            end
        end
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
"""
Double sha256 function
"""
function hash256(x::Vector{UInt8})
    return sha256(sha256(x))
end

"""
read_varint reads a variable integer from a stream
"""
 function read_varint(s::Base.GenericIOBuffer{Vector{UInt8}})
     i = UInt8[]
     readbytes!(s, i, 1)
     if i == [0xfd]
         # 0xfd means the next two bytes are the number
         readbytes!(s, i, 2)
         return reinterpret(Int16, i)[1]
     elseif i == [0xfe]
         # 0xfe means the next four bytes are the number
         readbytes!(s, i, 4)
         return reinterpret(Int32, i)[1]
     elseif i == [0xff]
         # 0xff means the next eight bytes are the number
         readbytes!(s, i, 8)
         return reinterpret(Int64, i)[1]
     else
         # anything else is just the integer
         return Int(i[1])
     end
 end

"""
Encodes an integer as a varint
"""
 function encode_varint(n::Integer)
    if n < 0xfd
        return [UInt8(n)]
    elseif n < 0x10000
        return prepend!(bytes(n, len=2, little_endian=true), [0xfd])
    elseif n < 0x100000000
        return prepend!(bytes(n, len=4, little_endian=true), [0xfe])
    elseif n < 0x10000000000000000
        return prepend!(bytes(n, len=8, little_endian=true), [0xff])
    else
        error("Integer, ", n, " is too large")
    end
 end

struct VarString <: AbstractString
    len::Integer
    str::String
    VarString(str::String) = new(length(str), str)
end

serialize(x::VarString) = append!(encode_varint(x.len), x.str)

io2varstring(io::IOBuffer) = VarString(String(read(io, read_varint(io))))

"""
    bytes2flags(bytes::Vector{UInt8}) -> Vector{Bool}

Returns an Vector{Bool} representing bits
"""
function bytes2flags(bytes::Vector{UInt8})
    result = Bool[]
    for byte in bytes
        for i in 0:7
            push!(result, (byte & (0x01 << i)) != 0)
        end
    end
    result
end

"""
    bytes2flags(bytes::Vector{UInt8}) -> Vector{Bool}

Returns an Vector{Bool} representing bits
"""
function flags2bytes(flags::Vector{Bool})
    if mod(length(flags), 8) != 0
        error("flags does not have a length that is divisible by 8")
    end
    result = fill(0x00, fld(length(flags), 8))
    i = 0
    for flag in flags
        byte_index, bit_index = fldmod(i, 8)
        if flag
            result[byte_index + 1] |= (1 << bit_index)
        end
        i += 1
    end
    result
end

# @deprecate read_varint(s::Base.GenericIOBuffer{Vector{UInt8}}) read(io::IOBuffer)::CompactSizeUInt
# @deprecate encode_varint(n::Integer) serialize(n::CompactSizeUInt)


text:
using Nextion
using Test

@testset "gauge" begin

    nexSerial = NexSerial("/dev/ttyUSB0")

    nexGauge = NexGauge(nexSerial, Name("z0"), cid=ComponentID(3))

    send(nexSerial, "page pg_gauge")

    nexGauge.backcolor = Colour.WHITE
    nexGauge.forecolor = Colour.RED

    #nexGauge.width = 3  # 0-5
    nexGauge.pointer.width = 3  # 0-5

    #for angle in 30:30:360
    #    sleep(1)
    #    nexGauge.value = angle
    #end
    for angle in 0:6:360
        sleep(0.1)
        nexGauge.value = (angle + 90) % 360
    end

    sleep(1)

    @test nexGauge.value == 90

    close(nexSerial)
    @test true

end


text:
using SwaggerMarkdown
using Test

const GET = "get"
const PUT = "put"
const POST = "post"

const ROOT = dirname(dirname(@__FILE__))
const TEST_ROOT = dirname(@__FILE__)
const SPECS_PATH = joinpath(TEST_ROOT, "specs")

include("utils.jl")
include("swagger_tests.jl")
include("parse_tests.jl")
include("build_tests.jl")
include("validation_tests.jl")


text:
mutable struct NDLCCBasis{T,M,P} <: Space{T}
    geo::M
    fns::Vector{Vector{Shape{T}}}
    pos::Vector{P}
end

NDLCCBasis(geo, fns) = NDLCCBasis(geo, fns, Vector{vertextype(geo)}(undef,length(fns)))

refspace(space::NDLCCBasis{T}) where {T} = NDLCCRefSpace{T}()
Base.similar(space::NDLCCBasis{T,M,P} where {T,M,P}, geo, fns, pos) = NDLCCBasis(geo, fns, pos)

"""
nedelecc3d(mesh, edges::Array{SArray{Tuple{2},Int64,1,2},1})

returns an object of the type NDLCCBasis.
"""

function nedelecc3d(mesh, edges)
    T = coordtype(mesh)
    P = vertextype(mesh)
    num_edges = numcells(edges)

    C = connectivity(edges, mesh, identity)
    rows = rowvals(C)
    vals = nonzeros(C)

    fns = Vector{Vector{Shape{T}}}(undef,num_edges)
    pos = Vector{P}(undef,num_edges)
    for (i,edge) in enumerate(cells(edges))

        fns[i] = Vector{Shape{T}}()
        pos[i] = cartesian(center(chart(edges,edge)))

        for k in nzrange(C,i)

            j = rows[k]
            s = vals[k]
            push!(fns[i], Shape{T}(j, abs(s), sign(s)))
        end
    end

    NDLCCBasis(mesh, fns, pos)
end

"""
nedelecc3d(mesh) returns an object of type NDLCCBasis
"""

function nedelecc3d(mesh)
    edges = skeleton(mesh,1)
    nedelecc3d(mesh, edges)
end

curl(space::NDLCCBasis, geo, fns) = NDLCDBasis(geo, fns, space.pos)

ttrace(X::NDLCCBasis, geo, fns) = RTBasis(geo, fns, deepcopy(X.pos))


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

function GitRevWalker(repo::GitRepo)
    w_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_revwalk_new, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}), w_ptr, repo.ptr)
    return GitRevWalker(repo, w_ptr[])
end

function Base.start(w::GitRevWalker)
    id_ptr = Ref(GitHash())
    err = ccall((:git_revwalk_next, :libgit2), Cint,
                (Ptr{GitHash}, Ptr{Void}), id_ptr, w.ptr)
    err != Int(Error.GIT_OK) && return (nothing, true)
    return (id_ptr[], false)
end

Base.done(w::GitRevWalker, state) = Bool(state[2])

function Base.next(w::GitRevWalker, state)
    id_ptr = Ref(GitHash())
    err = ccall((:git_revwalk_next, :libgit2), Cint,
                (Ptr{GitHash}, Ptr{Void}), id_ptr, w.ptr)
    err != Int(Error.GIT_OK) && return (state[1], (nothing, true))
    return (state[1], (id_ptr[], false))
end

Base.iteratorsize(::Type{GitRevWalker}) = Base.SizeUnknown()

function push_head!(w::GitRevWalker)
    @check ccall((:git_revwalk_push_head, :libgit2), Cint, (Ptr{Void},), w.ptr)
    return w
end

function Base.push!(w::GitRevWalker, cid::GitHash)
    @check ccall((:git_revwalk_push, :libgit2), Cint, (Ptr{Void}, Ptr{GitHash}), w.ptr, Ref(cid))
    return w
end

function Base.push!(w::GitRevWalker, range::AbstractString)
    @check ccall((:git_revwalk_push_range, :libgit2), Cint, (Ptr{Void}, Ptr{UInt8}), w.ptr, range)
    return w
end

function Base.sort!(w::GitRevWalker; by::Cint = Consts.SORT_NONE, rev::Bool=false)
    rev && (by |= Consts.SORT_REVERSE)
    ccall((:git_revwalk_sorting, :libgit2), Void, (Ptr{Void}, Cint), w.ptr, by)
    return w
end

repository(w::GitRevWalker) = w.owner

function Base.map(f::Function, walker::GitRevWalker;
                  oid::GitHash=GitHash(),
                  range::AbstractString="",
                  by::Cint = Consts.SORT_NONE,
                  rev::Bool=false,
                  count::Int=0)
    res = []
    sort!(walker, by=by, rev=rev)
    if !iszero(oid)
        push!(walker, oid)
    elseif !isempty(range)
        push!(walker, range)
    else
        push_head!(walker)
    end
    s = start(walker)

    c = 0
    repo = repository(walker)
    while !done(walker, s)
        val = f(s[1], repo)
        push!(res, val)
        val, s = next(walker, s)
        c +=1
        count == c && break
    end
    return res
end

function Base.count(f::Function, walker::GitRevWalker;
                  oid::GitHash=GitHash(),
                  by::Cint = Consts.SORT_NONE,
                  rev::Bool=false)
    c = 0
    sort!(walker, by=by, rev=rev)
    if !iszero(oid)
        push!(walker, oid)
    else
        push_head!(walker)
    end
    s = start(walker)

    repo = repository(walker)
    while !done(walker, s)
        val = f(s[1], repo)
        _, s = next(walker, s)
        c += (val == true)
    end
    return c
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
using NBInclude, Compat, Compat.Test

@testset "Scoping Tests" begin
    VERSION >= v"1.0.0" ? eval(:(@test_throws LoadError @nbinclude(joinpath("includes", "scoping.ipynb")))) : true
    @nbinclude(joinpath("includes", "scoping.ipynb"); softscope = true)
    @test a == 11
end

@test include(joinpath("includes", "test1.jl")) == 314159
@test f(5) == 6
@test normpath(myfile) == abspath("test.ipynb") * ":In[6]"
@test normpath(myfile2) == abspath("test.ipynb") * ":In[+7]"

x=[]; @nbinclude("test2.ipynb")
@test x == [1, 2, 3, 4, 5, 6]

x=[]; @nbinclude("test2.ipynb"; counters = [1, 4, 5])
@test x == [1, 4, 5]

x=[]; @nbinclude("test2.ipynb"; regex=r"#.*executeme")
@test x == [2, 4]

x=[]; @nbinclude("test2.ipynb"; counters = [1, 4, 5], regex=r"#.*executeme")
@test x == [4]

z = 0; @nbinclude("test2.ipynb"; anshook = x -> (global z += 1))
@test z == 6


text:
export Member

"""
A [`Guild`](@ref) member.
More details [here](https://discordapp.com/developers/docs/resources/guild#guild-member-object).
"""
struct Member
    user::Optional{User}
    nick::OptionalNullable{String}  # Not supposed to be nullable.
    roles::Vector{Snowflake}
    joined_at::DateTime
    premium_since::Optional{DateTime}
    deaf::Bool
    mute::Bool
end
@boilerplate Member :constructors :docs :lower :merge :mock


text:
"""
Initial query by controller receiving OF message
"""
function do_query!(msg::OFMessage,a::Agent,model)
    # If asset's network does not have any edge, there is no way to transport packets
    ignore = ne(a.ntw_graph) > 0 ? false : true
    path=[]

    if !ignore
        # src (from this sne) to dst
        query = (msg.dpid,msg.data.dst)
        query_time = model.ticks
        #existing precalc paths
        query_paths = a.paths
        #local graph to query
        query_graph = a.ntw_graph

        ####For benchmark#####
        sdir = model.data_dir
        record_benchmark!(sdir,model.run_label,a.id,query_time,query,query_graph,query_paths,model.benchmark)

        path = do_query(query_time,query,query_graph,query_paths)

        # log_info(model.ticks,a.id,"query: $(query) ----path found: $path ----> in precalc paths: $(query_paths)")

        # For no centralised control topo
        if model.ctrl_model != GraphModel(1) # ¬ centralised
            #check path validity as it is unknown how long ago the cache/graph has been updated
            if !isempty(path) && is_invalid_path(path,get_controlled_assets(a.id,model),model)
                path =  []
                clear_cache!(a,model)
            end

            #valid path not found?, need to send queries to nbs
            if isempty(path)
                #should always query nbs even if has already queried with no reply?
                sent_to = query_nbs!(a,msg,query_graph,query,[a.id],model)
                record_pending_query!(a,model.ticks,sent_to,query)
                push!(a.pending,(model.ticks,msg,false)) #false:not reprocess next time
            else
                clear_pending_query!(a,query)
            end
        end
        # update_agent_state!(a,msg,path)
        new_state = get_state(a)
        ## +1 query
        new_state.q_queries += 1.0
        ## +score path
        if !isempty(path) && msg.dpid != msg.data.dst
            push!(new_state.path_scores,(msg.dpid,msg.data.dst,path[2]))
        end
        set_state!(a,new_state)

    end
    return isempty(path) ? [] : last(path)
end

"""
    Query by neighbour control agent after receiving AGMessage
"""
function do_query!(msg::AGMessage,a::Agent,model)
    # log_info(model.ticks,a.id,"query msg is: $msg")

        # visited control ag
        trace = deepcopy(msg.body[:trace])
        push!(trace,a.id)
        query = msg.body[:query]

        if get_state(a).up
            # Aggregated information from other agents
            ## join graph received
            msg_ntw_g = create_subgraph(msg.body[:ntw_edgel],msg.body[:ntw_equiv],:eid)
            jg = join_subgraphs(a.ntw_graph,msg_ntw_g)

            # Update Knowledge Base
            # update local graph, however the problem it can grow too much
            # hence there is a periodic clear_cache to restore to original
            if model.clear_cache_graph_freq != 0
                a.ntw_graph = jg
            end

            query_time = model.ticks
            query_paths = a.paths
            query_graph = jg
            # path = Dict()

            #path for benchmark recording
            sdir = model.data_dir * model.run_label
            record_benchmark!(sdir,model.run_label,a.id,query_time,query,query_graph,query_paths,model.benchmark)

            #query

            path = do_query(query_time,query,query_graph,query_paths)

            if isempty(path)
                if !haskey(a.previous_queries,query)
                    sent_to = query_nbs!(a,msg,jg,query,trace,model)
                    record_pending_query!(a,model.ticks,sent_to,query)
                end
            else
                do_match!(path,msg,a,model)
                clear_pending_query!(a,query)
            end

            new_state = get_state(a)
            new_state.q_queries += 1.0
            set_state!(a,new_state)

        #TODO consider case where multiple paths are found
        else
            fw_query_nbs!(a,msg,trace,model)
        end

end

"""
Query local calculated paths and local graph
"""
function do_query(time::Int64,query::Tuple{Int64,Int64},lg::MetaGraph,paths::Dict{Tuple{Int64,Int64},Array{Tuple{Int64,Float64,Float64,Array{Int64}}}})
    path = []
    cp_paths = []
    lg_paths = []

    #query pre-calculated (cache) paths
    cp_paths = haskey(paths,query) ? paths[query] : []

    #query graph path regardless of cache, in case there is another
    # TODO: Do this only if cache path is too old
    path_state = query_paths(lg,query)

    #reshape paths found local graph
    for lg_path in path_state.paths
        default_confidence = 0.5
        push!(lg_paths,(time,default_confidence,last(path_state.dists),last(path_state.paths)))
    end

    #assumes query_paths is sorted by tick,score
    lg_ps = !isempty(lg_paths) ? first(lg_paths)[3] : -1.0
    cp_ps = !isempty(cp_paths) ? first(cp_paths)[3] : -1.0

    path = @match (lg_ps,cp_ps) begin
        #no path
        (-1.0,-1.0)  => []
        #path in local graph
        (a::Float64,-1.0), if a > 0 end => first(lg_paths)
        #path in cache
        (-1.0,b::Float64), if b > 0 end => first(cp_paths)
        #path in cache paths and local graph
        (_::Float64,_::Float64) => first(cp_paths)[3] < first(lg_paths)[3] ? first(cp_paths) : first(lg_paths) ##lower is better
    end
    return path
end

"""
Prepare query msg and send to neighbors
"""

function query_nbs!(a::Agent,msg::CTLMessage,query_graph::MetaGraph,query::Tuple{Int64,Int64},trace::Array{Int64},model::ABM)::Array{Int64}
    ntw_edgel = [ e for e in edges(query_graph) if src(e) <  dst(e) ]
    ntw_equiv = [(v,query_graph[v,:eid]) for v in vertices(query_graph)]
    # when length == 1, msg is the original OFMessage otherwise is AGMessage
    of_mid = length(trace) == 1 ? msg.id : msg.body[:of_mid]
    body = Dict(:query=>query,:trace=>trace,:ntw_edgel => ntw_edgel, :ntw_equiv=>ntw_equiv, :of_mid=>of_mid)
    msg_template = AGMessage(-1,model.ticks,a.id,-1,QUERY_PATH,body)
    return send_to_nbs!(msg_template,a,model)
end

"""
    Forward query to neighbors when control agent is down
"""

function fw_query_nbs!(a::Agent,msg::AGMessage,trace::Array{Int64},model::ABM)
    of_mid = msg.body[:of_mid]
    nbody = Dict(:query=>msg.body[:query],:trace=>trace,:ntw_edgel => msg.body[:ntw_edgel], :ntw_equiv=>msg.body[:ntw_equiv], :of_mid=>of_mid)
    msg_template = AGMessage(-1,model.ticks,a.id,-1,QUERY_PATH,nbody)
    send_to_nbs!(msg_template,a,model)
end

"""
records when control ag send a query to other control ags
"""

function record_pending_query!(a::Agent,time::Int64,sent_to::Array{Int64},query::Tuple{Int64,Int64})
    #track query sent, if not previously tracked
    if !haskey(a.previous_queries,query)
        a.previous_queries[query] = (time,sent_to)
    end
end

"""
clear when a solution is found/received to a query
"""

function clear_pending_query!(a::Agent,query::Tuple{Int64,Int64})
    #no need to clear pending?
    new_previous_queries = Dict()
    #path found, need to remove potential pending queries
    for k in keys(a.previous_queries)
        if k != query #exclude current query as path has been found
            new_previous_queries[k] = a.previous_queries[k]
        end
    end
    a.previous_queries = new_previous_queries
end

"""
clear queries requested by sne to controller that have expired, so the same query can be requested after some time model.ofmsg_reattempt
"""
function clear_pending_query!(a::SimNE,model::ABM)
    #no need to clear pending?
    new_requested_ctl = Dict()
    #path found, need to remove potential pending queries
    for k in keys(a.requested_ctl)
        if model.ticks - a.requested_ctl[k] < model.ofmsg_reattempt #
            new_requested_ctl[k] = a.requested_ctl[k]
        end
    end
    a.requested_ctl = new_requested_ctl
end

"""
When a solution is received for a query it marks the pending msg
so it can be reprocessed in the next tick
"""
function mark_reprocess_of_msg!(a::Agent,msg::AGMessage)
    #reprocess of msg right after, to do local query with new path found
    new_pending = []
    for p in a.pending
        if p[2].id == msg.body[:of_mid]
            push!(new_pending,(p[1],p[2],true))#reprocess next tick
        else
            push!(new_pending,p)
        end
    end
    a.pending = new_pending
end


text:
struct Loss
    target::TargetLoss
    regularization::Float32
end

gettotal(x::Loss) = gettotal(x.target) + x.regularization
gettarget(x::Loss) = x.target
getregularization(x::Loss) = x.regularization
getpolicy(x::Loss) = getpolicy(x.target)
getreward(x::Loss) = getreward(x.target)
getvalue(x::Loss) = getvalue(x.target)


text:
module GAIO

using LinearAlgebra
using SparseArrays
using StaticArrays
# using GLFW
# using ModernGL
using GeometryBasics
using LightGraphs
using ForwardDiff
using Arpack
using Base.Threads
using Base.Cartesian: @nexprs, @ntuple
using Base: @propagate_inbounds
using MuladdMacro
using HostCPUFeatures
using SIMD

using GLMakie
using WGLMakie

export Box

export AbstractBoxPartition, BoxPartition, TreePartition
export dimension, depth, key_to_box, point_to_key, tree_search

export BoxSet
export boxset_empty, subdivide, subdivide!

export BoxFun

export TransferOperator
export strongly_connected_components, matrix, eigs

export BoxMap, PointDiscretizedMap
export SampledBoxMap, AdaptiveBoxMap
export boxmap

export map_boxes, map_boxes_new

export rk4, rk4_flow_map

export relative_attractor, unstable_set!, chain_recurrent_set
export cover_roots, finite_time_lyapunov_exponents

export plot

# ENV["JULIA_DEBUG"] = all

include("box.jl")

abstract type AbstractBoxPartition{B <: Box} end

include("partition_regular.jl")
include("partition_tree.jl")
include("boxset.jl")
include("simd_helper.jl")
include("boxmap.jl")
include("boxfun.jl")
include("transfer_operator.jl")
include("algorithms.jl")
include("plot.jl")

end # module


text:
using Diversity.API
using EcoBase: AbstractAssemblage, AbstractThings, AbstractPlaces

"""
    gettypes(m::AbstractAssemblage)

Returns the AbstractTypes component of the metacommunity.
"""
function gettypes(m::AbstractAssemblage)
    return _gettypes(m)
end

"""
    getpartition(m::AbstractAssemblage)

Returns the AbstractPartition component of the metacommunity.

"""
function getpartition(m::AbstractAssemblage)
        return _getpartition(m)
end

"""
    counttypes(m::AbstractAssemblage[, raw::Bool = false])
    counttypes(t::AbstractTypes[, raw::Bool = false])

Returns number of types in an `AbstractTypes` object or the
`AbstractAssemblage` containing it. `raw` determines whether to
count the number of raw or processed types, which varies, for instance,
when the types are determined by a phylogeny.
"""
function counttypes end

function counttypes(m::AbstractAssemblage, raw::Bool = false)
    return _counttypes(_gettypes(m), raw)
end

function counttypes(t::AbstractThings, raw::Bool = false)
    return _counttypes(t, raw)
end

"""
    countsubcommunities(m::AbstractAssemblage)
    countsubcommunities(p::AbstractPartition)

Returns number of subcommunities in an `AbstractPartition` object or the
`AbstractAssemblage` containing it.
"""
function countsubcommunities end

function countsubcommunities(m::AbstractAssemblage)
    return _countsubcommunities(_getpartition(m))
end

function countsubcommunities(p::AbstractPartition)
    return _countsubcommunities(p)
end

"""
    gettypenames(m::AbstractAssemblage[, raw::Bool = false])
    gettypenames(t::AbstractTypes[, raw::Bool = false])

Returns the names of the types of the `AbstractTypes` object or the
`AbstractAssemblage` containing it. `raw` determines whether to
count the number of raw or processed types, which varies, for instance,
when the types are determined by a phylogeny.
"""
function gettypenames end

function gettypenames(m::AbstractAssemblage, raw::Bool = false)
    return _gettypenames(_gettypes(m), raw)
end

function gettypenames(t::AbstractTypes, raw::Bool = false)
    return _gettypenames(t, raw)
end

"""
    getdiversityname(m::AbstractAssemblage)
    getdiversityname(t::AbstractTypes)

Returns the name of the diversity type used.
"""
function getdiversityname end

getdiversityname(t::AbstractThings) = _getdiversityname(t)
getdiversityname(m::AbstractAssemblage) = _getdiversityname(_gettypes(m))

"""
    addedoutputcols(m::AbstractAssemblage)
    addedoutputcols(t::AbstractTypes)

Returns the name of any additional columns needed to disambiguate the
diversity type used.
"""
function addedoutputcols end

addedoutputcols(t::AbstractThings) = _addedoutputcols(t)
addedoutputcols(m::AbstractAssemblage) = addedoutputcols(_gettypes(m))

"""
    getaddedoutput(::AbstractTypes)

Returns the contents of any additional columns to be added to outputs.
"""
function getaddedoutput end

getaddedoutput(t::AbstractTypes) = _getaddedoutput(t)
getaddedoutput(m::AbstractAssemblage) = _getaddedoutput(_gettypes(m))

"""
    getsubcommunitynames(m::AbstractAssemblage)
    getsubcommunitynames(p::AbstractPartition)

Returns the names of the subcommunities in an `AbstractPartition` object or the
`AbstractAssemblage` containing it.
"""
function getsubcommunitynames end

function getsubcommunitynames(m::AbstractAssemblage)
    return _getsubcommunitynames(_getpartition(m))
end

function getsubcommunitynames(p::AbstractPartition)
    return _getsubcommunitynames(p)
end

"""
    getabundance(m::AbstractAssemblage, raw::Bool)

Returns the abundances array of the metacommunity.
"""
function getabundance(m::AbstractAssemblage, raw::Bool = false)
    return _getabundance(m, raw)
end

"""
    getmetaabundance(m::AbstractAssemblage)

Returns the metacommunity abundances of the metacommunity.
"""
function getmetaabundance(m::AbstractAssemblage, raw::Bool = false)
    return _getmetaabundance(m, raw)
end

"""
    getweight(m::AbstractAssemblage)

Returns the subcommunity weights of the metacommunity.
"""
function getweight(m::AbstractAssemblage)
    return _getweight(m)
end

"""
    getordinariness!(m::AbstractAssemblage)

Returns (and possibly calculates) the ordinariness array of the subcommunities.
"""
function getordinariness!(m::AbstractAssemblage)
    return _getordinariness!(m)
end

"""
    getmetaordinariness!(m::AbstractAssemblage)

Returns (and possibly calculates) the ordinariness of the
metacommunity as a whole.
"""
function getmetaordinariness!(m::AbstractAssemblage)
    return _getmetaordinariness!(m)
end

"""
    hassimilarity(t::AbstractAssemblage)
    hassimilarity(t::AbstractThings)

Is there similarity of some non-trivial type in the object?
"""
function hassimilarity end
hassimilarity(t::AbstractThings) = _hassimilarity(t)
hassimilarity(asm::AbstractAssemblage) = hassimilarity(gettypes(asm))

"""
    calcsimilarity(t::AbstractTypes, scale::Real)

Retrieves (and possibly calculates) a similarity matrix from t.
"""
function calcsimilarity(t::AbstractTypes, scale::Real)
    return _calcsimilarity(t, scale)
end

function createsummaryline(vec::AbstractVector{<:AbstractString})
    linefunc(vec) = mapreduce(x->x*", ", *, vec[1:(end-1)])*vec[end]
    length(vec) == 1 && return vec[1]
    length(vec) < 6 && return linefunc(vec)
    linefunc(vec[1:3])*"..."*linefunc(vec[(end-1):end])
end

import Base.show
function show(io::IO, mc::AbstractMetacommunity)
    sp = createsummaryline(gettypenames(mc))
    si = createsummaryline(getsubcommunitynames(mc))
    println(io, "$(typeof(mc)) with $(counttypes(mc)) species in $(countsubcommunities(mc)) sites measuring $(getdiversityname(gettypes(mc))) diversity.\n\nSpecies names:\n$(sp)\n\nSite names:\n$(si)")
end


text:
import Base.LinAlg.LU, Base.getindex, Base.setindex!

@doc """
`A_ldiv_B_md!(dest, F, src, dim)` solves a tridiagonal system along dimension `dim` of `src`,
storing the result in `dest`. Currently, `F` must be an LU-factorized tridiagonal matrix.
If desired, you may safely use the same array for both `src` and `dest`, so that this becomes an
in-place algorithm.
""" ->
function A_ldiv_B_md!(dest, F, src, dim::Integer)
    1 <= dim <= max(ndims(dest),ndims(src)) || throw(DimensionMismatch("The chosen dimension $dim is larger than $(ndims(src)) and $(ndims(dest))"))
    n = size(F, 1)
    n == size(src, dim) && n == size(dest, dim) || throw(DimensionMismatch("Sizes $n, $(size(src,dim)), and $(size(dest,dim)) do not match"))
    size(dest) == size(src) || throw(DimensionMismatch("Sizes $(size(dest)), $(size(src)) do not match"))
    check_matrix(F)
    R1 = CartesianRange(size(dest)[1:dim-1])
    R2 = CartesianRange(size(dest)[dim+1:end])
    _A_ldiv_B_md!(dest, F, src, R1, R2)
end
_A_ldiv_B_md(F, src, R1::CartesianRange, R2::CartesianRange) =
    _A_ldiv_B_md!(similar(src, promote_type(eltype(F), eltype(src))), F, src, R1, R2)

# Solving along the first dimension
function _A_ldiv_B_md!{T,CI<:CartesianIndex{0}}(dest, F::LU{T,Tridiagonal{T}}, src,  R1::CartesianRange{CI}, R2::CartesianRange)
    n = size(F, 1)
    dl = F.factors.dl
    d  = F.factors.d
    du = F.factors.du
    # Forward substitution
    @inbounds for I2 in R2
        dest[1, I2] = src[1, I2]
        for i = 2:n       # note: cannot use @simd here!
            dest[i, I2] = src[i, I2] - dl[i-1]*dest[i-1, I2]
        end
    end
    # Backward substitution
    dinv = 1./d
    @inbounds for I2 in R2
        dest[n, I2] /= d[n]
        for i = n-1:-1:1  # note: cannot use @simd here!
            dest[i, I2] = (dest[i, I2] - du[i]*dest[i+1, I2])*dinv[i]
        end
    end
    dest
end

# Solving along any other dimension
function _A_ldiv_B_md!{T}(dest, F::LU{T,Tridiagonal{T}}, src, R1::CartesianRange, R2::CartesianRange)
    n = size(F, 1)
    dl = F.factors.dl
    d  = F.factors.d
    du = F.factors.du
    # Forward substitution
    @inbounds for I2 in R2
        @simd for I1 in R1
            dest[I1, 1, I2] = src[I1, 1, I2]
        end
        for i = 2:n
            @simd for I1 in R1
                dest[I1, i, I2] = src[I1, i, I2] - dl[i-1]*dest[I1, i-1, I2]
            end
        end
    end
    # Backward substitution
    dinv = 1./d
    for I2 in R2
        @simd for I1 in R1
            dest[I1, n, I2] *= dinv[n]
        end
        for i = n-1:-1:1
            @simd for I1 in R1
                dest[I1, i, I2] = (dest[I1, i, I2] - du[i]*dest[I1, i+1, I2])*dinv[i]
            end
        end
    end
    dest
end

function check_matrix{T}(F::LU{T,Tridiagonal{T}})
    n = size(F,1)
    for i = 1:n
        F.ipiv[i] == i || error("For efficiency, pivoting is not supported")
    end
    nothing
end


text:
using ModelingToolkit, DifferentialEquations, Plots
@variables t
D = Differential(t)

function dynamics(; name, τ=0.01, a_max=300.0, V₀=300.0, x_T=5E3, y_T=0.0)
    states = @variables x(t) y(t) γ(t) a(t) a_cmd(t) r(t) σ(t) λ(t) V(t)
    ps = @parameters τ = τ a_max = a_max x_T = x_T y_T = y_T

    eqs = [
        D(x) ~ V * cos(γ)
        D(y) ~ V * sin(γ)
        D(γ) ~ a / V
        D(a) ~ (min(max(a_cmd, -a_max), a_max) - a) / τ
        r ~ sqrt((x_T - x)^2 + (y_T - y)^2)
        λ ~ atan(y_T - y, x_T - x)
        σ ~ γ - λ
        V ~ V₀
    ]

    return ODESystem(eqs, t, states, ps; name)
end

function BPNG(; name, σ_lim=deg2rad(60), m=10.0, δ=0.01, k=4.0, n=1, p=0.8, N=3.0, r₀=5E3, γ_f_d=-deg2rad(50.0))
    states = @variables r(t) γ(t) σ(t) λ(t) λ̇(t) V(t) e_γ_f(t) K_r(t) K_η(t) a_cmd(t)
    ps = @parameters σ_lim = σ_lim m = m δ = δ k = k n = n p = p N = N r₀ = r₀ γ_f_d = γ_f_d

    eqs = [
        λ̇ ~ -V * sin(σ) / r
        e_γ_f ~ (N * λ - γ) / (N - 1.0) - γ_f_d
        K_r ~ (N - 1.0 + δ) * (1.0 + r / r₀)^n
        K_η ~ 1.0 + k * (1.0 - abs(sin(σ) / sin(σ_lim))^m)
        a_cmd ~ N * V * λ̇ + V^2 * cos(σ) * K_r * K_η / r * abs(e_γ_f)^p * sign(e_γ_f)
    ]

    return ODESystem(eqs, t, states, ps; name)
end

function RPG(; name, σ_lim=deg2rad(60), m=2.0, n=3.0, γ_f_d=-deg2rad(50.0))
    states = @variables r(t) γ(t) σ(t) λ(t) λ′(t) λ″(t) V(t) a_cmd(t) c_m(t) c_n(t)
    ps = @parameters σ_lim = σ_lim m = m n = n γ_f_d = γ_f_d

    eqs = [
        c_m ~ (-n * (λ - γ_f_d) + tan(σ)) / (m - n) / r^m
        c_n ~ (m * (λ - γ_f_d) - tan(σ)) / (m - n) / r^n
        # λ′      ~ tan(σ)/r
        λ′ ~ n * c_n * r^(n - 1) + m * c_m * r^(m - 1)
        λ″ ~ n * (n - 1) * c_n * r^(n - 2) + m * (m - 1) * c_m * r^(m - 2)
        a_cmd ~ -V^2 * cos(σ) * ((r * λ″ + tan(σ) / r) / (sec(σ))^2 + λ′)
    ]

    return ODESystem(eqs, t, states, ps; name)
end

indexof(sym, syms) = findfirst(isequal(sym), syms)

function condition_stop(u, t, integrator)
    u_sys = states(sys_simplified)
    p_sys = parameters(sys_simplified)

    x = u[indexof(model.x, u_sys)]
    y = u[indexof(model.y, u_sys)]
    γ = u[indexof(model.γ, u_sys)]
    x_T = integrator.p[indexof(model.x_T, p_sys)]
    y_T = integrator.p[indexof(model.y_T, p_sys)]
    σ_lim = integrator.p[indexof(command.σ_lim, p_sys)]

    r = sqrt((x_T - x)^2 + (y_T - y)^2)
    λ = atan(y_T - y, x_T - x)
    σ = γ - λ

    return (r < 0.5) || (σ > σ_lim) #|| (r < 10 && ṙ >= 0)
end
affect!(integrator) = terminate!(integrator)
cb_stop = DiscreteCallback(condition_stop, affect!)
# @time sol = solve(prob, Tsit5(), callback = cb_stop)

##
case = 6       # [manual input required]
γ_f_d_deg = -170.0

if case <= 3
    if case == 1
        n_list = 1
        p_list = [0.6, 0.8, 1.0]
    elseif case == 2
        n_list = [0, 1, 2]
        p_list = 0.9
    elseif case == 3
        n_list = 0
        p_list = 0.6
    end

    @named command = BPNG()
    pv_list = [[command.n => n, command.p => p, command.γ_f_d => deg2rad(γ_f_d_deg)] for n in n_list for p in p_list]
else
    if case == 4
        n_list = [1.5, 2.0, 3.0]
        m_list = 1.0
    elseif case == 5
        n_list = 3.0
        m_list = [1.0, 1.5, 2.0]
    elseif case == 6
        n_list = 1.5
        m_list = 1.0
    end

    @named command = RPG()
    pv_list = [[command.n => n, command.m => m, command.γ_f_d => deg2rad(γ_f_d_deg)] for n in n_list for m in m_list]
end

@named model = dynamics()

eqs_c = [
    model.V ~ command.V
    model.r ~ command.r
    model.γ ~ command.γ
    model.σ ~ command.σ
    model.λ ~ command.λ
    model.a_cmd ~ command.a_cmd
]

sys = compose(ODESystem(eqs_c, name=:sys), model, command)
sys_simplified = structural_simplify(sys)

IC = [model.x => 0.0
    model.y => 0.0
    model.γ => deg2rad(55)
    model.a => 0.0]
t_span = (0.0, 30.0)
# pv = [command.n => 1
#       command.p => 1.0]
# prob = ODEProblem(sys_simplified, IC, t_span, pv; jac = true)

fs_val = 14
lw_val = 2.5
ls_list = [:solid :dash :dot :dashdot]

f_xy = plot()
f_r = plot()
f_σ = plot()
f_γ = plot()
f_e_γ_f = plot()
f_a = plot()
f_a_cmd = plot()

for pv in pv_list
    prob = ODEProblem(sys_simplified, IC, t_span, pv; jac=true)
    @time sol = solve(prob, solver=Tsit5(), callback=cb_stop, reltol=1e-8, abstol=1e-8)

    if case <= 3
        lgnd = "\$n = $(pv[1][2])\$, \$p = $(pv[2][2])\$"
        plot!(f_e_γ_f, sol, vars=rad2deg(command.e_γ_f), label=lgnd, xlabel="\$t\$ [s]", ylabel="\$e_{\\gamma_{f}}\$ [deg]", linewidth=lw_val, linestyle=ls_list[findfirst(isequal(pv), pv_list)], legendfontsize=fs_val, guidefontsize=fs_val, xtickfontsize=fs_val, ytickfontsize=fs_val)
    else
        lgnd = "\$n = $(pv[1][2])\$, \$m = $(pv[2][2])\$"
        plot!(f_e_γ_f, sol, vars=rad2deg(model.γ)-γ_f_d_deg, label=lgnd, xlabel="\$t\$ [s]", ylabel="\$\\gamma - \\gamma_{f}\$ [deg]", linewidth=lw_val, linestyle=ls_list[findfirst(isequal(pv), pv_list)], legendfontsize=fs_val, guidefontsize=fs_val, xtickfontsize=fs_val, ytickfontsize=fs_val)
    end

    plot!(f_xy, sol, vars=(model.x, model.y), aspect_ratio=:equal, label=lgnd, xlabel="\$x\$ [m]", ylabel="\$y\$ [m]", linewidth=lw_val, linestyle=ls_list[findfirst(isequal(pv), pv_list)], legendfontsize=fs_val, guidefontsize=fs_val, xtickfontsize=fs_val, ytickfontsize=fs_val)

    plot!(f_r, sol, vars=model.r, label=lgnd, xlabel="\$t\$ [s]", ylabel="\$r\$ [m]", linewidth=lw_val, linestyle=ls_list[findfirst(isequal(pv), pv_list)], legendfontsize=fs_val, guidefontsize=fs_val, xtickfontsize=fs_val, ytickfontsize=fs_val)

    plot!(f_σ, sol, vars=rad2deg(model.σ), label=lgnd, xlabel="\$t\$ [s]", ylabel="\$\\sigma\$ [deg]", linewidth=lw_val, linestyle=ls_list[findfirst(isequal(pv), pv_list)], legendfontsize=fs_val, guidefontsize=fs_val, xtickfontsize=fs_val, ytickfontsize=fs_val, ylims=(0, 60))

    plot!(f_γ, sol, vars=rad2deg(model.γ), label=lgnd, xlabel="\$t\$ [s]", ylabel="\$\\gamma\$ [deg]", linewidth=lw_val, linestyle=ls_list[findfirst(isequal(pv), pv_list)], legendfontsize=fs_val, guidefontsize=fs_val, xtickfontsize=fs_val, ytickfontsize=fs_val)

    plot!(f_a, sol, vars=model.a, label=lgnd, xlabel="\$t\$ [s]", ylabel="\$a\$ [m/s\$^{2}\$]", linewidth=lw_val, linestyle=ls_list[findfirst(isequal(pv), pv_list)], legendfontsize=fs_val, guidefontsize=fs_val, xtickfontsize=fs_val, ytickfontsize=fs_val)

    plot!(f_a_cmd, sol, vars=model.a_cmd, label=lgnd, xlabel="\$t\$ [s]", ylabel="\$a_{cmd}\$ [m/s\$^{2}\$]", linewidth=lw_val, linestyle=ls_list[findfirst(isequal(pv), pv_list)], legendfontsize=fs_val, guidefontsize=fs_val, xtickfontsize=fs_val, ytickfontsize=fs_val)

    println("\n sum(abs(a)) = $(sum(abs.(sol[model.a])))\n")
end

plot!(f_σ, [0.0, 20.0], [60, 60], linecolor=:red, linestyle=:dashdot, label="\$\\sigma_{\\lim}\$")

fig_dir = "Simulation/Julia_BPNG/Figures"
mkpath(fig_dir)

savefig(f_xy, joinpath(fig_dir, "Fig_xy_$(case).pdf"))
savefig(f_r, joinpath(fig_dir, "Fig_r_$(case).pdf"))
savefig(f_σ, joinpath(fig_dir, "Fig_sigma_$(case).pdf"))
savefig(f_γ, joinpath(fig_dir, "Fig_gamma_$(case).pdf"))
savefig(f_e_γ_f, joinpath(fig_dir, "Fig_e_gamma_f_$(case).pdf"))
savefig(f_a, joinpath(fig_dir, "Fig_a_$(case).pdf"))
savefig(f_a_cmd, joinpath(fig_dir, "Fig_a_cmd_$(case).pdf"))

display(f_xy)
display(f_r)
display(f_σ)
display(f_γ)
display(f_e_γ_f)
display(f_a)
display(f_a_cmd)