text:
"""
    stresstest(device_or_devices)

Run a GPU stress test (matrix multiplication) on one or multiple GPU devices, as specified by the positional argument.
If no argument is provided (only) the currently active GPU will be used.

**Keyword arguments:**

Choose one of the following (or none):
* `duration`: stress test will take about the given time in seconds. (StressTestBatched)
* `enforced_duration`: stress test will take almost precisely the given time in seconds. (StressTestEnforced)
* `approx_duration`: stress test will hopefully take approximately the given time in seconds. No promises made! (StressTestFixedIter)
* `niter`: stress test will run the given number of matrix-multiplications, however long that will take. (StressTestFixedIter)
* `mem`: number (`<:Real`) between 0 and 1, indicating the fraction of the available GPU memory that should be used, or a `<:UnitPrefixedBytes` indicating an absolute memory limit. (StressTestStoreResults)

General settings:
* `dtype` (default: `Float32`): element type of the matrices
* `monitoring` (default: `false`): enable automatic monitoring, in which case a [`MonitoringResults`](@ref) object is returned.
* `size` (default: `2048`): matrices of size `(size, size)` will be used
* `verbose` (default: `true`): toggle printing of information
* `parallel` (default: `true`): If `true`, will (try to) run each GPU test on a different Julia thread. Make sure to have enough Julia threads.
* `threads` (default: `nothing`): If `parallel == true`, this argument may be used to specify the Julia threads to use.
* `clearmem` (default: `false`): If `true`, we call [`clear_all_gpus_memory`](@ref) after the stress test.
* `io` (default: `stdout`): set the stream where the results should be printed.

When `duration` is specifiec (i.e. [`StressTestEnforced`](@ref)) there is also:
* `batch_duration` (default: `ceil(Int, duration/10)`): desired duration of one batch of matmuls.
"""
function stresstest(
    devices;
    mem=nothing,
    dtype=Float32,
    verbose=true,
    size=_default_matsize_stresstest(),
    duration=nothing,
    enforced_duration=nothing,
    approx_duration=nothing,
    niter=nothing,
    clearmem=false,
    monitoring=false,
    batch_duration=nothing,
    io::IO=stdout,
    kwargs...,
)
    logger = ConsoleLogger(io)

    Base.with_logger(logger) do
        if eltype(devices) != CuDevice
            throw(ArgumentError("Unknown devices iterator / collection?!"))
        end
        if all(isnothing, [enforced_duration, duration, approx_duration, niter, mem])
            duration = 60 # default
        end
        if (monitoring || ismonitoring()) && Threads.nthreads() < length(devices) + 1
            ngpus = length(devices)
            error(
                "To test $ngpus GPUs while monitoring requires $(ngpus + 1) Julia threads. Only $(Threads.nthreads()) available.",
            )
        end

        if !isnothing(duration)
            # StressTestBatched
            verbose && @info("Will run for about $(duration) seconds on each GPU!")
            ts = [
                StressTestBatched(dev; duration, dtype, verbose, size, batch_duration) for
                dev in devices
            ]
        elseif !isnothing(niter)
            # StressTestFixedIter
            verbose && @info("Will run $(niter) iterations on each GPU.")
            ts = [StressTestFixedIter(dev; niter, dtype, verbose, size) for dev in devices]
        elseif !isnothing(approx_duration)
            # StressTestFixedIter
            verbose && @info(
                "Will try to run for $(duration) seconds on each GPU, but no promises!"
            )
            ts = [
                StressTestFixedIter(dev; approx_duration, dtype, verbose, size) for
                dev in devices
            ]
        elseif !isnothing(enforced_duration)
            # StressTestEnforced
            verbose && @info(
                "Will run for almost precisely $(enforced_duration) seconds on each GPU."
            )
            ts = [
                StressTestEnforced(dev; enforced_duration, dtype, verbose, size) for
                dev in devices
            ]
        elseif !isnothing(mem)
            # StressTestStoreResults (gpu-burn-like)
            verbose && @info("Will run a `StressTestStoreResults`.")
            ts = [StressTestStoreResults(dev; mem, dtype, verbose, size) for dev in devices]
        end
        monitoring && monitoring_start(; devices=devices, verbose)
        Δt = @elapsed _run_stresstests(ts; verbose, kwargs...)
        if clearmem
            verbose && @info("Clearing GPU memory.")
            clear_all_gpus_memory(devices)
        end
        verbose && @info("Took $(round(Δt; digits=2)) seconds to run the tests.")
        if monitoring
            results = monitoring_stop(; verbose)
            return results
        else
            return nothing
        end
    end
end
stresstest(device::CuDevice=CUDA.device(); kwargs...) = stresstest([device]; kwargs...)

function _default_threads(ntests)
    nthreads = Threads.nthreads()
    monitoring = ismonitoring()
    if nthreads >= ntests + monitoring + 1
        # if there are enough "worker" threads
        # use them and skip the main thread (1).
        return collect(2:(ntests + 1))
    elseif nthreads >= ntests + monitoring
        return collect(1:ntests)
    else
        # not enough threads
        @warn("There are too few Julia threads to run all tests in parallel!")
        return [mod1(i, nthreads - monitoring) for i in 1:ntests] # round-robin distribution
    end
end

function _run_stresstests(
    tests; parallel=true, threads=_default_threads(length(tests)), verbose=true
)
    if parallel == true
        if length(threads) != length(tests)
            throw(ArgumentError("length(threads) != length(tests)"))
        end
        @sync for (i, test) in enumerate(tests)
            @tspawnat threads[i] test(; verbose=verbose)
        end
    else
        for test in tests
            test(; verbose=verbose)
        end
    end
    return nothing
end

# function stresstest_on_worker(devices; verbose=true, worker=first(workers()), duration=nothing, kwargs...)
#     verbose && @info("Using worker $worker.")
#     # pid = withenv("JULIA_NUM_THREADS"=>length(devices)) do
#     #     pid, = addprocs(1, exeflags = "--project=$(Base.active_project())")
#     #     return pid
#     # end

#     t = @spawnat worker stresstest(devices; verbose=false, duration, kwargs...)
#     monitor_temperature(duration+5; liveplot=true)
#     return fetch(t)
# end


text:
# FASTQ Read
# ==========
#
# Design discussions in https://github.com/BioJulia/FASTX.jl/pull/35

struct FASTQRead{A <: BioSequences.Alphabet}
    identifier::String
    description::String
    sequence::BioSequences.LongSequence{A}
    quality::Vector{UInt8} # in raw PHRED scores, not offset
end

"""
    FASTQ.Read(record::FASTQ.Record, offset::Integer=33)
    FASTQ.Read(record::FASTQ.Record, encoding_name::Symbol)

Create a FASTQ read object from a FASTQ.Record.
The FASTQ.Read has fields:
    identifier::String
    description::String
    sequence::BioSequences.LongSequence
    quality::Vector{UInt8} # in raw PHRED scores, not offset
"""
function Read(record::FASTQ.Record, offset::Integer=33)
    FASTQRead(
        identifier(record),
        description(record),
        sequence(record),
        quality(record, offset)
    )
end

function Read(record::FASTQ.Record, encoding_name::Symbol)
    FASTQRead(
        identifier(record),
        description(record),
        sequence(record),
        quality(record, encoding_name)
    )
end

"""
    sequence(read::FASTQ.FASTQRead)
Get the sequence of a FASTQ read. Same as read.sequence
"""
function sequence(read::FASTQ.FASTQRead)
    read.sequence
end

"""
    quality(read::FASTQ.FASTQRead)
Get the quality of a FASTQ read (Vector{UInt8}). Same as read.quality
"""
function quality(read::FASTQ.FASTQRead)
    read.quality
end

"""
    length(read::FASTQ.FASTQRead)
Get the length of a FASTQ read.
"""
function Base.length(read::FASTQ.FASTQRead)
    length(read.sequence)
end

"""
     Base.getindex(record::Record, i::UnitRange{Int})
Subset a FASTQRead using string syntax. Eg. read[3:7] or read[3]
"""
function Base.getindex(read::FASTQ.FASTQRead, i::UnitRange{<:Integer})
    return FASTQRead(
        join([read.identifier,string(i[1],"..",i[end])],"_"),
        join([read.description,string(i[1],"..",i[end])]," "),
        read.sequence[i],
        read.quality[i]
    )
end

function Base.getindex(read::FASTQRead, i::Integer)
    read[i:i]
end

function Base.show(io::IO, read::FASTQ.FASTQRead)
    print(io, summary(read), ':')
    println(io)
    println(io, "   identifier: ", read.identifier)
    println(io, "  description: ", read.description)
    println(io, "     sequence: ", String(read.sequence))
    print(io, "      quality: ", Int64.(quality(read)))
end

ext: jl

max_forks_repo_name: BioJulia/FASTX


text:
# Fusion trees:
#==============================================================================#
"""
    struct FusionTree{I, N, M, L, T}

Represents a splitting tree of simple objects of type `I<:Sector`, splitting a coupled
object to `N` uncoupled objects.

The `isdual` field indicates whether an isomorphism ``Z`` is present or not.

The field `uncoupled` contains the objects coming out of the splitting trees, before the
possible ``Z`` isomorphism.

`M=max(0, N-2)` is the number of inner lines.

`L=max(0, N-1)` is the number of vertices.

`T` is the type of the labels for the vertices. `T = Nothing` if
`FusionStyle(I) == UniqueFusion()` or `SimpleFusion()`.
"""
struct FusionTree{I<:Sector, N, M, L, T}
    uncoupled::NTuple{N, I}
    coupled::I
    isdual::NTuple{N, Bool}
    innerlines::NTuple{M, I}
    vertices::NTuple{L, T}
    function FusionTree{I, N, M, L, T}(
            uncoupled::NTuple{N, I},
            coupled::I,
            isdual::NTuple{N, Bool},
            innerlines::NTuple{M, I},
            vertices::NTuple{L, T}) where {I<:Sector, N, M, L, T}
        new{I, N, M, L, T}(uncoupled, coupled, isdual, innerlines, vertices)
    end
end

"""
    fusiontreetype(::Type{I}, N::Int) where {I<:Sector}

Return the correct fusiontree type `FusionTree{I<:Sector, N, M, L, T}` based on `I` and `N`.
"""
function fusiontreetype(::Type{I}, N::Int) where {I<:Sector}
    if N === 0
        FusionTree{I, 0, 0, 0, vertex_labeltype(I)}
    elseif N === 1
        FusionTree{I, 1, 0, 0, vertex_labeltype(I)}
    else
        FusionTree{I, N, N-2, N-1, vertex_labeltype(I)}
    end
end

function FusionTree{I}(uncoupled::NTuple{N, Any},
                       coupled,
                       isdual::NTuple{N, Bool},
                       innerlines,
                       vertices = ntuple(n->nothing, max(0, N-1))) where {I<:Sector, N}
    if FusionStyle(I) isa GenericFusion
        fusiontreetype(I, N)(map(s->convert(I, s), uncoupled),
            convert(I, coupled), isdual, map(s->convert(I, s), innerlines), vertices)
    else
        vertices′ = ntuple(n->nothing, max(0, N-1))
        if vertices == vertices′ || all(isone, vertices)
            fusiontreetype(I, N)(map(s->convert(I, s), uncoupled),
                convert(I, coupled), isdual, map(s->convert(I, s), innerlines), vertices′)
        else
            throw(ArgumentError("Incorrect fusion vertices"))
        end
    end
end

function FusionTree(uncoupled::NTuple{N, I},
                    coupled::I,
                    isdual::NTuple{N, Bool},
                    innerlines,
                    vertices = ntuple(n->nothing, max(0, N-1))) where {I<:Sector, N}
    if FusionStyle(I) isa GenericFusion
        fusiontreetype(I, N)(uncoupled, coupled, isdual, innerlines, vertices)
    else
        vertices′ = ntuple(n->nothing, max(0, N-1))
        if vertices == vertices′ || all(isone, vertices)
            fusiontreetype(I, N)(uncoupled, coupled, isdual, innerlines, vertices′)
        else
            throw(ArgumentError("Incorrect fusion vertices"))
        end
    end
end

"""
    _abelianinner(outer::NTuple{N, I}) where {I<:Sector, N}

Return the innerlines as `NTuple{N-3,I}` in the Abelian case.

Return `()` if ``N ≤ 3``.

The input `outer` corresponds to `(uncoupled..., dual(coupled))`.
"""
_abelianinner(outer::Tuple{}) = ()
_abelianinner(outer::Tuple{I}) where {I<:Sector} =
    outer[1] == one(I) ? () : throw(SectorMismatch())
_abelianinner(outer::Tuple{I, I}) where {I<:Sector} =
    outer[1] == dual(outer[2]) ? () : throw(SectorMismatch())
_abelianinner(outer::Tuple{I, I, I}) where {I<:Sector} =
    first(⊗(outer...)) == one(I) ? () : throw(SectorMismatch())
function _abelianinner(outer::NTuple{N, I}) where {I<:Sector, N}
    c = first(outer[1] ⊗ outer[2])
    return (c, _abelianinner((c, TupleLabXD.tail2(outer)...))...)
end

function FusionTree{I}(uncoupled::NTuple{N}, coupled = one(I),
                        isdual = ntuple(n->false, N)) where {I<:Sector, N}
    FusionStyle(I) isa UniqueFusion ||
        error("fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`")
    FusionTree{I}(map(s->convert(I, s), uncoupled), convert(I, coupled), isdual,
                    _abelianinner(map(s->convert(I, s), (uncoupled..., dual(coupled)))))
end

function FusionTree(uncoupled::NTuple{N, I}, coupled::I = one(I),
                        isdual = ntuple(n->false, N)) where {I<:Sector, N}
    FusionStyle(I) isa UniqueFusion ||
        error("fusion tree requires inner lines if `FusionStyle(I) <: MultipleFusion`")
    FusionTree{I}(uncoupled, coupled, isdual, _abelianinner((uncoupled..., dual(coupled))))
end

# Properties
"""
    sectortype(f::FusionTree) -> Sector
    sectortype(::Type{<:FusionTree{I}}) -> Sector

Return the sector type of the fusion tree.
"""
sectortype(::Type{<:FusionTree{I}}) where {I<:Sector} = I
sectortype(f::FusionTree) = sectortype(typeof(f))

"""
    FusionStyle(f::FusionTree) -> FusionStyle
    FusionStyle(::Type{<:FusionTree{I}}) -> FusionStyle

Return the fusion style of the fusion tree.
"""
FusionStyle(::Type{<:FusionTree{I}}) where {I<:Sector} = FusionStyle(I)
FusionStyle(f::FusionTree) = FusionStyle(typeof(f))

"""
    BraidingStyle(f::FusionTree) -> BraidingStyle
    BraidingStyle(::Type{<:FusionTree{I}}) -> BraidingStyle

Return the braiding style of the fusion tree.
"""
BraidingStyle(::Type{<:FusionTree{I}}) where {I<:Sector} = BraidingStyle(I)
BraidingStyle(f::FusionTree) = BraidingStyle(typeof(f))

"""
    length(f::FusionTree) -> Int
    length(::Type{<:FusionTree{<:Sector, N}}) -> Int

Extend `Base.length`. Return the number of the uncoupled objects in a fusion tree.
"""
Base.length(::Type{<:FusionTree{<:Sector, N}}) where {N} = N
Base.length(f::FusionTree) = length(typeof(f))

# Hashing, important for using fusion trees as key in a dictionary
function Base.hash(f::FusionTree{I}, h::UInt) where {I}
    h = hash(f.isdual, hash(f.coupled, hash(f.uncoupled, h)))
    if FusionStyle(I) isa MultipleFusion
        h = hash(f.innerlines, h)
    end
    if FusionStyle(I) isa GenericFusion
        h = hash(f.vertices, h)
    end
    return h
end

function Base.isequal(f1::FusionTree{I, N}, f2::FusionTree{I, N}) where {I<:Sector, N}
    f1.coupled == f2.coupled || return false
    @inbounds for i = 1:N
        f1.uncoupled[i] == f2.uncoupled[i] || return false
        f1.isdual[i] == f2.isdual[i] || return false
    end
    if FusionStyle(I) isa MultipleFusion
        @inbounds for i=1:N-2
            f1.innerlines[i] == f2.innerlines[i] || return false
        end
    end
    if FusionStyle(I) isa GenericFusion
        @inbounds for i=1:N-1
            f1.vertices[i] == f2.vertices[i] || return false
        end
    end
    return true
end
Base.isequal(f1::FusionTree, f2::FusionTree) = false

# converting to actual array
function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, 0}) where {I}
    X = convert(A, fusiontensor(one(I), one(I), one(I)))[1, 1, :]
    return X
end

function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, 1}) where {I}
    c = f.coupled
    if f.isdual[1]
        sqrtdc = sqrtdim(c)
        Zcbartranspose = sqrtdc * convert(A, fusiontensor(conj(c), c, one(c)))[:, :, 1, 1]
        X = conj!(Zcbartranspose) # we want Zcbar^†
    else
        X = convert(A, fusiontensor(c, one(c), c))[:, 1, :, 1, 1]
    end
    return X
end

function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I, 2}) where {I}
    a, b = f.uncoupled
    isduala, isdualb = f.isdual
    c = f.coupled
    da, db, dc = dim.((a, b, c))
    μ = (FusionStyle(I) isa GenericFusion) ? f.vertices[1] : 1
    C = convert(A, fusiontensor(a, b, c))[:, :, :, μ]
    X = C
    if isduala
        Xtemp = X
        X = similar(Xtemp)
        Za = convert(A, FusionTree((a,), a, (isduala,), ()))
        TO.contract!(1, Za, :N, Xtemp, :N, 0, X, (1,), (2,), (2, 3), (1,), (1,2,3))
    end
    if isdualb
        Xtemp = X
        X = similar(Xtemp)
        Zb = convert(A, FusionTree((b,), b, (isdualb,), ()))
        TO.contract!(1, Zb, :N, Xtemp, :N, 0, X, (1,), (2,), (1, 3), (2,), (2,1,3))
    end
    return X
end

function Base.convert(A::Type{<:AbstractArray}, f::FusionTree{I,N}) where {I,N}
    tailout = (f.innerlines[1], TupleLabXD.tail2(f.uncoupled)...)
    isdualout = (false, TupleLabXD.tail2(f.isdual)...)
    ftail = FusionTree(tailout, f.coupled, isdualout,
                        Base.tail(f.innerlines), Base.tail(f.vertices))
    Ctail = convert(A, ftail)
    f1 = FusionTree((f.uncoupled[1], f.uncoupled[2]), f.innerlines[1],
                    (f.isdual[1], f.isdual[2]), (), (f.vertices[1],))
    C1 = convert(A, f1)
    dtail = size(Ctail)
    d1 = size(C1)
    X = similar(C1, (d1[1], d1[2], Base.tail(dtail)...))
    trivialtuple = ntuple(identity, Val(N))
    TO.contract!(1, C1, :N, Ctail, :N, 0, X,
                    (1,2), (3,), Base.tail(trivialtuple), (1,), (trivialtuple..., N+1))
    return X
end

# Show methods
function Base.show(io::IO, t::FusionTree{I, N, M, K, Nothing}) where {I<:Sector, N, M, K}
    print(IOContext(io, :typeinfo => I), "FusionTree{", type_repr(I), "}(",
        t.uncoupled, ", ", t.coupled, ", ", t.isdual, ", ", t.innerlines, ")")
end
function Base.show(io::IO, t::FusionTree{I}) where {I<:Sector}
    print(IOContext(io, :typeinfo => I), "FusionTree{", type_repr(I), "}(",
        t.uncoupled, ", ", t.coupled, ", ", t.isdual, ",",
        t.innerlines, ", ", t.vertices, ")")
end

# Manipulate fusion trees
include("manipulations.jl")

# Fusion tree iterators
include("iterator.jl")


text:
using BeaData
using Base.Test
using Compat

try
    open(joinpath(homedir(),".beadatarc"), "r") do f
        key = readstring(f)
    end
    key = rstrip(key)
    include("test_with_key.jl")
catch err
end


text:
module SimpleMemoizeSample

export @simplememoize

macro simplememoize(ex)
    @assert Meta.isexpr(ex, :function)
    fname = ex.args[1].args[1]
    fname_escaped = esc(fname)
    fargs = esc.(ex.args[1].args[2:end])
    fbody = esc(ex.args[2])
    _cache = gensym(Symbol(fname, "_cache"))
    quote
        const $_cache = Dict{NTuple{$(length(fargs))}, Any}()
        function $fname_escaped($(fargs...))
            get!($_cache, ($(Expr(:tuple, fargs...)))) do
                $fbody
            end
        end
    end
end

end # module


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

""" Base class for defining points on a lattice. Lattice points are defined by basis vectors eᵢ and lattice indices indexᵢ:

point = Σᵢ indexᵢ*eᵢ

Example: define a 2D basis.

```julia
julia> a = LatticeBasis([1,2],[3,4])
LatticeBasis{2, Int64}(SVector{2, Int64}[[1, 2], [3, 4]])
```

Array indexing is used to generate a lattice point: a[1,2] = 1*[1,2] + 2*[3,4]

```julia
julia> a[1,2]
2-element SVector{2, Int64} with indices SOneTo(2):
  7
 10
```

Subtypes supporting the AbstractBasis interface should implement these functions:

Returns the neighbors in ring n surrounding centerpoint, excluding centerpoint
```
neighbors(::Type{B},centerpoint::Tuple{T,T},neighborhoodsize::Int) where{T<:Real,B<:AbstractBasis}
```
Returns the lattice basis vectors that define the lattice
```
basis(a::S) where{S<:AbstractBasis}
```
Returns the vertices of the unit polygon that tiles the plane for the basis
```
tilevertices(a::S) where{S<:AbstractBasis}
```
"""
abstract type AbstractBasis{N,T<:Real} end
export AbstractBasis

function Base.getindex(A::B1, indices::Vararg{Int, N}) where{N,T,B1<:AbstractBasis{N,T}}
    return basismatrix(A)*SVector{N,Int}(indices)
end

#     temp::SVector{N,T} = (basis(A)*SVector{N,Int}(indices))::SVector{N,T}
#     return temp
# end

Base.setindex!(A::AbstractBasis{N,T}, v, I::Vararg{Int, N}) where{T,N} = nothing #can't set lattice points. Might want to throw an exception instead.

struct LatticeBasis{N,T<:Real} <: AbstractBasis{N,T}
    basisvectors::SMatrix{N,N,T}

    """ Convenience constructor that lets you use tuple arguments to describe the basis instead of SVector

    Example:
    ```
    basis = LatticeBasis{2}((1.0,0.0),(0.0,1.0)) #basis for rectangular lattice
    ```
    This allocates 48 bytes for a 2D basis in Julia 1.6.2. It  shouldn't allocate anything but not performance critical.
    """
    function LatticeBasis(vectors::Vararg{NTuple{N,T},N}) where{T,N}
        temp = MMatrix{N,N,T}(undef)

        for (j,val) in pairs(vectors)
            for i in 1:N
                temp[i,j] = val[i]
            end
        end

        return new{N,T}(SMatrix{N,N,T}(temp))
    end

     LatticeBasis(vectors::SMatrix{N,N,T}) where{N,T<:Real} = new{N,T}(vectors)
end
export LatticeBasis

function basismatrix(a::LatticeBasis{N,T})::SMatrix{N,N,T,N*N} where{N,T}
    return a.basisvectors
end

"""Can access any point in a lattice so the range of indices is unlimited"""
function Base.size(a::LatticeBasis{N,T}) where{N,T}
    return ntuple((i)->Base.IsInfinite(),N)
end


text:
module StackOverflow

using HTTP
using CodecZlib
using JSON

include("questions.jl")
export
    getrecentquestionsfortag

include("privileges.jl")
export
    checkmyprivileges,
    getprivileges

include("utils.jl")
export
    convert_HTTP_Response_To_JSON

include("answers.jl")
export
    getanswers

include("badges.jl")
export
    getrecentbadgerecipients

include("siteinfo.jl")
export
    getsiteinfo

include("search.jl")
export
    searcherror

include("tags.jl")
export
    searchtag

include("users.jl")
export
    getusers

include("comments.jl")
export
    getcomments

include("suggestededits.jl")
export
    getsuggested_edits

include("errors.jl")
export
    geterrors

end # module


text:
###
#   Simple implmentation of the Metropolis Hasting algorithm in julia
#
#            Institute for Risk and Uncertainty, Uni of Liverpool
#
#                       Authors: Ander Gray, Adolphus Lye
#
#                       Email: Ander.Gray@liverpool.ac.uk,
#                              Adolphus.Lye@liverpool.ac.uk
#
#
#
#       W. K. Hastings (1970). Monte Carlo sampling methods using
#       Markov chains and their applications. Biometrika,
#       57(1), 97-109. doi:10.1093/biomet/57.1.97
#
###

###
#   Add description of function and inputs
###

function metropolis_hastings(
    Target::Function,
    Prop::Function,
    start::Vector{<:Real},
    Nsamples::Integer,
    burnin::Integer=50,
    thin::Integer=3;
    islogged::Bool=true)

    dims = length(start)                                    # Dimensions of input/ prior

    PropRnd = mu -> rand(Prop(mu))                          # Generates a sample from the proposal given mean mu

    chain = zeros(Nsamples * thin + burnin, dims)
    chain[1,:] = start
    accRate = 0

    islogged ? evalDen =  x -> exp(Target(x)) : evalDen = x -> Target(x)

    for i = 2:( Nsamples * thin + burnin)

        next = PropRnd(chain[i - 1, :])              # Draw candidate

        targDen = evalDen(next)                  # Target Density at next sample
        targPrevious = evalDen(chain[i - 1, :])     # Target Density at current sample

        propDen = pdf(Prop(chain[i - 1, :]), next)         # Proposal at next centred at current
        propPrevious = pdf(Prop(next), chain[i - 1, :])    # Propsoal at current centred at next

        α = targDen / targPrevious * propPrevious / propDen     # General formula for acceptance probability

        accepted = α >= rand()

        if accepted
            chain[i, :] = next
            accRate = accRate + 1
        else
            chain[i, :] = chain[i - 1, :]
        end

    end
    accRate = accRate / (Nsamples * thin + burnin)
    return chain[burnin + 1:thin:end,:], accRate
end

##
#   Symetric proposal and logged target
##
function metropolis_hastings_simple(
    Target::Function,
    PropRnd::Function,
    start::AbstractVector{<:Real},
    Nsamples::Integer,
    burnin::Integer=50,
    thin::Integer=3)

    dims = length(start)                                    # Dimensions of input/ prior
    chain = zeros(dims, Nsamples * thin + burnin)
    chain[:, 1] = start[:]
    accRate = 0

    for i = 2:( Nsamples * thin + burnin)

        next = PropRnd(chain[:, i - 1])                 # Draw candidate

        targDen = Target(next)[1]                    # Target Density at next sample
        targPrevious = Target(chain[:, i - 1])[1]       # Target Density at current sample

        α =  min(0, targDen - targPrevious)

        accepted = α >= log(rand())

        if accepted
            chain[:, i] = next
            accRate = accRate + 1
        else
            chain[:, i] = chain[:, i - 1]
        end

    end
    accRate = accRate / (Nsamples * thin + burnin)
    return chain[:, burnin + 1:thin:end]', accRate
end


text:
function _all_translation_vectors(dim::Int, max_component::Int)::Array{Int, 2}
    out = zeros(Int, dim, (max_component+1)^dim)
    @inbounds for i in axes(out, 2)
        k = i-1
        for l in dim:-1:1
            k, out[l, i] = divrem(k, max_component+1)
        end
    end
    return out
end

function _all_signed_translation_vectors(tv::Array{Int, 2})::Array{Int, 3}
    dim = size(tv, 1)
    stv = zeros(Int, dim, size(tv, 2), (2^dim))

    signs = 2*_all_translation_vectors(dim, 1) .- 1

    @inbounds for i in axes(signs, 2)
        @views stv[:, :, i] = signs[:, i] .* tv
    end
    return stv
end

# acts as a bias term before computing the norm
function _add_unit_cell_index(lattice::AbstractLattice, tv::Array{Int, 2})::Array{Int, 3}
    dim = size(tv, 1)
    ucsize = unitcellsize(lattice)
    utv = zeros(Int, dim+1, size(tv, 2), ucsize)

    for u in 1:ucsize
        utv[1:end-1, :, u] = tv
        utv[end, :, u] .= u
    end
    return utv
end

function _norm2(::Val{true}, lattice::AbstractLattice, v::AbstractVector)
    uv = unitcell_translation_vectors(lattice)
    u = v[1:end-1] + uv[:, Int(v[end])]
    return dot(u, metric(lattice), u)
end

_norm2(::Val{false}, lattice::AbstractLattice, v::AbstractVector) = dot(v, metric(lattice), v)

norm2(lattice::AbstractLattice, v::AbstractVector) = _norm2(hasunitcell(lattice), lattice, v)
norm(lattice::AbstractLattice, v::AbstractVector) = sqrt(norm2(lattice, v))

function _approx_uniq(x::AbstractVector, maximum_k::Int)
    uniq_ = filter(!iszero, sort!(unique(x)))
    uniq = Float64[]
    for i in axes(uniq_, 1)
        if length(uniq) == 0 || !(last(uniq) ≈ uniq_[i])
            push!(uniq, uniq_[i])
        end
        if length(uniq) >= maximum_k
            break
        end
    end
    return uniq
end

function _generate_all_translation_vectors!(::Val{true}, lattice::AbstractLattice, maximum_k::Int)
    N = ndims(lattice)

    tv = _all_translation_vectors(N, maximum_k)
    tv = _all_signed_translation_vectors(tv)
    tv = reshape(tv, N, :)
    tv = unique(tv; dims=2)
    tv = _add_unit_cell_index(lattice, tv)
    tv = reshape(tv, size(tv, 1), :)

    d = lattice.translation_vectors

    uv = unitcell_translation_vectors(lattice)

    for u in 1:unitcellsize(lattice)
        norms = [
            norm2(lattice, [v[1:end-1] - uv[:, u]; v[end]])
            for v in eachcol(tv)
        ]

        p = sortperm(norms)

        norms = norms[p]
        utv = tv[:, p]

        uniq = _approx_uniq(norms, maximum_k)

        for k in 1:maximum_k
            idx = findall(uniq[k] .≈ norms)
            d[(k, u)] = [Coordinate(v...) for v in eachcol(@views utv[:, idx])]
        end
    end

    return d
end

function _generate_all_translation_vectors!(::Val{false}, lattice::AbstractLattice, maximum_k::Int)
    N = ndims(lattice)

    tv = _all_translation_vectors(N, maximum_k)
    if !ismetricdiag(lattice)
        tv = _all_signed_translation_vectors(tv)
        tv = reshape(tv, N, :)
        tv = unique(tv; dims=2)
    end

    norms = [norm2(lattice, v) for v in eachcol(tv)]
    p = sortperm(norms)

    norms = norms[p]
    tv = tv[:, p]
    if ismetricdiag(lattice)
        stv = _all_signed_translation_vectors(tv)
    else
        stv = reshape(tv, size(tv)..., 1)
    end

    uniq = _approx_uniq(norms, maximum_k)

    d = lattice.translation_vectors
    for k in 1:maximum_k
        idx = findall(uniq[k] .≈ norms)
        d[k] = ([
            Coordinate(v...)
            for tv in eachslice(@view stv[:, idx, :]; dims=3)
            for v in eachcol(tv)
        ])
    end
    return d
end

precompute_translation_vectors!(lattice, maximum_k) = (
    maximum_k <= 0
    ? throw(ArgumentError("maximum_k must be positive!"))
    : _generate_all_translation_vectors!(hasunitcell(lattice), lattice, maximum_k)
)

function _translation_vectors(::Val{false}, lattice::AbstractLattice, k::Int, ::Nothing)
    if !(k in keys(lattice.translation_vectors))
        precompute_translation_vectors!(lattice, k)
    end
    return lattice.translation_vectors[k]
end

function _translation_vectors(::Val{true}, lattice::AbstractLattice, k::Int, u::Int)
    if !(k in first.(keys(lattice.translation_vectors)))
        precompute_translation_vectors!(lattice, k)
    end
    return lattice.translation_vectors[(k, u)]
end

function _translation_vectors(::Val{true}, lattice::AbstractLattice, k::Int, ::Nothing)
    if !(k in first.(keys(lattice.translation_vectors)))
        precompute_translation_vectors!(lattice, k)
    end
    d = Dict{Int, Vector{Coordinate{ncoordinates(lattice), Int}}}()
    for u in 1:unitcellsize(lattice)
        d[u] = lattice.translation_vectors[(k, u)]
    end
    return d
end

function translation_vectors(lattice::AbstractLattice, k::Int, u::Union{Int, Nothing}=nothing)
    return _translation_vectors(hasunitcell(lattice), lattice, k, u)
end

_neighbors(b::Val{false}, lattice::AbstractLattice, site::Coordinate, k::Int) = map(v -> site + v, _translation_vectors(b, lattice, k, nothing))
function _neighbors(b::Val{true}, lattice::AbstractLattice, site::Coordinate, k::Int)
    tv = _translation_vectors(b, lattice, k, site.coordinates[end])

    ns = Vector{Coordinate{ncoordinates(lattice), Int}}(undef, length(tv))

    for i in axes(ns, 1)
        ns[i] = Coordinate(
            map(+, site.coordinates[1:end-1], tv[i].coordinates[1:end-1])...,
            tv[i].coordinates[end]
        )
    end
    return ns
end

function neighbors(lattice::AbstractLattice, site::Coordinate{N, Int}, ::Val{K}) where {N,K}
    ns = _neighbors(hasunitcell(lattice), lattice, site, K)
    return filter(!isnothing, [apply_boundary_conditions(lattice, n) for n in ns])
end

neighbors(lattice::AbstractLattice, site::Coordinate, k::Int) = neighbors(lattice, site, Val{k}())
neighbors(lattice::AbstractLattice, site::Coordinate) = neighbors(lattice, site, Val{1}())


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -14.06559781669544 <= q <= 13.37751741699121)

                     add_NL_constraint(m, :(sigmoid(-0.7556338259975872 + -0.2850076752045205*sigmoid(0.7432962759699406 + 0.616365769922441*sigmoid(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*sigmoid(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*sigmoid(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.5819748667402531*sigmoid(0.6367931081099356 + 0.08862555157308405*sigmoid(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*sigmoid(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*sigmoid(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + 0.4871037750371179*sigmoid(0.12597161846431115 + 0.10770543162761159*sigmoid(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*sigmoid(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*sigmoid(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) + sigmoid(0.9186311390475046 + -0.8559712370251158*sigmoid(0.7432962759699406 + 0.616365769922441*sigmoid(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*sigmoid(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*sigmoid(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.2512416329037399*sigmoid(0.6367931081099356 + 0.08862555157308405*sigmoid(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*sigmoid(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*sigmoid(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + 0.7592758779705528*sigmoid(0.12597161846431115 + 0.10770543162761159*sigmoid(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*sigmoid(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*sigmoid(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) + sigmoid(-0.7151630404444922 + -0.15417874070690862*sigmoid(0.7432962759699406 + 0.616365769922441*sigmoid(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*sigmoid(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*sigmoid(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.016895941959414795*sigmoid(0.6367931081099356 + 0.08862555157308405*sigmoid(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*sigmoid(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*sigmoid(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.17047547651198514*sigmoid(0.12597161846431115 + 0.10770543162761159*sigmoid(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*sigmoid(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*sigmoid(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using DEBmicroTrait
using CSV, DataFrames, Statistics
using JLD
using HypothesisTests

########################################
# I/O
df_isolates     = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates2traits.csv", DataFrame, missingstring="N/A")
########################################

########################################
Genome_size     = convert(Array{Float64,1}, df_isolates.Genome_size)
V_cell          = DEBmicroTrait.genome_size_to_cell_volume(Genome_size)
Min_gen_time    = df_isolates.Min_gen_time
gmax            = log(2)./Min_gen_time
Gram_stain      = convert(Array{String,1}, df_isolates.gram_stain)########################################

########################################
dry_mass        = 0.47*DEBmicroTrait.cell_volume_to_dry_mass(V_cell, gmax, Gram_stain)
ρ_bulk          = 1.0 # g/cm^3
N_cells         = 1e3
Bio_0           = N_cells*1e6*ρ_bulk*dry_mass./12.011
########################################

########################################
# I/O
save("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates_batch_init.jld", "Bio0", Bio_0, "Md", dry_mass, "rhoB", ρ_bulk, "Ncells", N_cells)
########################################

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
function solve(problem::LeastSquaresProblem, x, method::Union{Adam, AdaMax, LineSearch, ParticleSwarm, TrustRegion}, options::LeastSquaresOptions)
    td = OptimizationProblem(problem.residuals, problem.bounds)
    res = solve(td, x, method, OptimizationOptions())
end


text:
# Benchmark on weighted sampling

using BenchmarkLite
using StatsBase

import StatsBase: direct_sample!, alias_sample!, xmultinom_sample!

### procedure definition

type WSampleProc{Alg} <: Proc end

abstract WithRep
abstract NoRep

type Direct <: WithRep end
tsample!(s::Direct, wv, x) = direct_sample!(1:length(wv), wv, x)

type Alias <: WithRep end
tsample!(s::Alias, wv, x) = alias_sample!(1:length(wv), wv, x)

type Xmultinom_S <: WithRep end
tsample!(s::Xmultinom_S, wv, x) = shuffle!(xmultinom_sample!(1:length(wv), wv, x))

type Xmultinom <: WithRep end
tsample!(s::Xmultinom, wv, x) = xmultinom_sample!(1:length(wv), wv, x)

type Direct_S <: WithRep end
tsample!(s::Direct_S, wv, x) = sort!(direct_sample!(1:length(wv), wv, x))

type Sample_WRep <: WithRep end
tsample!(s::Sample_WRep, wv, x) = sample!(1:length(wv), wv, x; ordered=false)

type Sample_WRep_Ord <: WithRep end
tsample!(s::Sample_WRep_Ord, wv, x) = sample!(1:length(wv), wv, x; ordered=true)

# config is in the form of (n, k)

Base.string{Alg}(p::WSampleProc{Alg}) = lowercase(string(Alg))

Base.length(p::WSampleProc, cfg::(Int, Int)) = cfg[2]
Base.isvalid{Alg<:WithRep}(p::WSampleProc{Alg}, cfg::(Int, Int)) = ((n, k) = cfg; n >= 1 && k >= 1)
Base.isvalid{Alg<:NoRep}(p::WSampleProc{Alg}, cfg::(Int, Int)) = ((n, k) = cfg; n >= k >= 1)

function Base.start(p::WSampleProc, cfg::(Int, Int))
    n, k = cfg
    x = Array(Int, k)
    w = weights(fill(1.0/n, n))
    return (w, x)
end

Base.run{Alg}(p::WSampleProc{Alg}, cfg::(Int, Int), s) = tsample!(Alg(), s[1], s[2])
Base.done(p::WSampleProc, cfg, s) = nothing

### benchmarking

const ns = 5 * (2 .^ [0:9])
const ks = 2 .^ [1:16]

## with replacement

const procs1 = Proc[ WSampleProc{Direct}(),
                     WSampleProc{Alias}(),
                     WSampleProc{Xmultinom_S}(),
                     WSampleProc{Sample_WRep}(),
                     WSampleProc{Xmultinom}(),
                     WSampleProc{Direct_S}(),
                     WSampleProc{Sample_WRep_Ord}() ]

const cfgs1 = vec([(n, k) for k in ks, n in ns])

rtable1 = run(procs1, cfgs1; duration=0.2)
println()

## show results

println("Sampling With Replacement")
println("===================================")
show(rtable1; unit=:mps, cfghead="(n, k)")
println()


text:
using Test
include("../../learner/RELIEF.jl")

data = [-1 0 1; 1 1 1; -1 -3 1; 1 -2 2; 9 -1 2; 1 2 2; 10 20 1]
labels = ['-', '+',  '-',   '-',  '-', '+' , '+' ]

@testset "Relief algoritmo" begin
    w = RELIEF(data, labels)
    @test length(w) == length(data[1,:])
    @test minimum(w) >= 0
    @test maximum(w) <= 1
end


text:
using PackageCompiler
create_sysimage([:Pluto, :PlutoUI,
                 :CodeTracking, :Combinatorics,
                 :DataFrames, :Distributions,
                 :FixedEffectModels, :ForwardDiff, :GLM,
                 :InteractiveUtils, :LineSearches,
                 :LinearMaps, :Markdown, :NLsolve,
                 :OffsetArrays, :Optim, :PrettyTables,
                 :Revise, :StatsPlots];
                #precompile_execution_file = "warmup.jl",
                precompile_statements_file = "combined_trace.jl",
                #replace_default = true,
                sysimage_path="/home/jovyan/sysimage.so",
                cpu_target = PackageCompiler.default_app_cpu_target())


text:
export ObjectAgent_MAG, calcMAGForce, jac!

const Gxy_params = [1.6602e-09, 0.0166];

mutable struct ObjectAgent_MAG{T,U} <: ObjectAgent{T,U}
    name::String;
    pos::Tuple{T, T};
    Fdes::Tuple{T,T};
    actList::Array{Tuple{U, U}, 1};
    act_used::Array{U, 1};
    Gxy::Array{T, 2};
    neighbors::Array{Tuple{ObjectAgent{T,U}, Array{U, 1}, Array{U, 1}},1};
    xk::Vector{T};
    zk::Vector{T};
    uk::Vector{T};
    rcvBuffer_x::Vector{T};
    rcvBuffer_u::Vector{T};
    rcvBuffer_N::Vector{T};
    J::Array{T, 2};
    fvk::Vector{T};
    decDir::Vector{T};
    Fdes_sc::T;
    D2::Array{T, 2};
    Q::Cholesky{T,Array{T,2}};
    R::Matrix{T};
    r::Vector{T};
    b0::Vector{T};
    b1::Array{T, 2};
    A_kkt::Array{T, 2};
    b_kkt::Array{T,2};

    function ObjectAgent_MAG(name::String,
        pos::Tuple{T, T},
        Fdes::Tuple{T, T},
        aa::ActuatorArray{T,U},
        actList::Array{Tuple{U, U}, 1},
        act_used::Array{U, 1},
        λ::T) where {T<:Real, U<:Unsigned}

        neighbrs = Array{Tuple{ObjectAgent{T,U}, Array{U, 1}, Array{U, 1}},1}[];

        # Compute Γ and Λ_a matrices matrix
        Gxy = genGxy(aa, pos, actList);

        # Normalize the optimization problem
        Fdes_sc = norm(Fdes);
        Fdes = Fdes./Fdes_sc;
        Gxy /= Fdes_sc;

        Gxy_red = Gxy[act_used, :];

        # number of actuators
        n = length(actList);

        Q = cholesky(Gxy_red*Gxy_red' + 1/λ*UniformScaling(1)); # G*G' + 1/λ*I
        b0 = -Gxy_red*[Fdes[1];Fdes[2]]; # b0 = -J*oa.Fdes
        b1 = Gxy_red*Gxy'; # b1 = J*oa.Gxy'

        Atilde = [Gxy_red;UniformScaling(1/sqrt(λ))];
        R = zeros(T, 2,2);
        LinSolvers.qr_r!(Atilde, R);
        r = zeros(T,2);

        # Matrices for the variant assuring the direction of the generated force
        c = 1e1;
        n_red = length(act_used);
        kktA =  vcat([(1/λ + c)*I zeros(n_red,1) -Gxy_red],
                     [zeros(1,n_red) 1 [Fdes[1] Fdes[2]]],
                     [-Gxy_red' [Fdes[1] Fdes[2]]' zeros(2,2)]);
        kktb = vcat(zeros(T, n_red, 1), [1.0; 0.0; 0.0]);

        new{T,U}(name, pos, Fdes, actList, act_used, Gxy, neighbrs,
        zeros(T, n), # xk
        zeros(T, n), # zk
        zeros(T, n), # uk
        zeros(T, n), # buffer xk
        zeros(T, n), # buffer uk
        zeros(T, n), # buffer N
        Gxy_red, # J
        zeros(T, 2), # fvk
        zeros(T, length(act_used)), # decDir
        Fdes_sc,
        Gxy_red*Gxy_red', # D2
        Q,
        R,
        r,
        b0,
        b1,
        kktA,
        kktb
        );
    end
end

function jac!(oa::ObjectAgent_MAG)
    return nothing
end

function Gxyf( x, y )
    val = -Gxy_params[1]*[x; y] ./ ((x.^2 + y.^2 + (Gxy_params[2]).^2).^3);
end

function genGxy(aa::ActuatorArray,
    pos::Tuple{T, T},
    al::Array{Tuple{U,U}, 1} = vec([(i,j) for j in oneunit(aa.ny):aa.ny, i in oneunit(aa.nx):aa.nx])) where {T<:Real, U<:Unsigned}

    Gxy = zeros(length(al), 2);
    for (k, a_pos) in enumerate(al)
            xa, ya = actuatorPosition(aa, a_pos[1:2]);
            Gxy[k, :] = Gxyf(pos[1] - xa, pos[2] - ya);
    end

    Gxy
end

function mag_fi( currents::Array{T, 2} ) where T<:Real
    -13.6.*abs.(currents).^3 + 11.1.*currents.^2;
end

function mag_fi_inv( fi2::Array{T, 2} ) where T<:Real
    # fi2 must be in the interval (0,1.1)
    a = -13.6;
    b = 11.1;
    d = Complex.(-fi2);

    D0 = b^2;
    D1 = 2*b^3 .+ 27*a^2*d;
    C = ((D1 .- sqrt.(D1.^2 .- 4*D0^3)) ./2 ).^(1/3);

    zeta = (-1/2 + 1/2*sqrt(3)*im);
    i = -1/(3*a) *(b .+ zeta .* C .+ D0 ./ (zeta .* C));
    ires = real.(i);
end

function calcMAGForce(aa::ActuatorArray{T,U}, pos::Tuple{T,T}, currents::Array{T,2}) where {T<:Real, U<:Unsigned}
    Gxy = genGxy(aa, pos);

    # currents = mag_fi(currents);
    currents[isnan.(currents)] .= 0;

    return tuple((Gxy'*vec(currents'))...);
end

function updatex_dirFixed!(oa::ObjectAgent_MAG{T,U}, λ::T) where {T<:Real, U<:Unsigned}
    oa.xk .= oa.zk .- oa.uk;

    oa.b_kkt[end-1:end] .= oa.Gxy' * oa.xk;

    sol = oa.A_kkt\oa.b_kkt;
    oa.xk[oa.act_used] .+= sol[1:end-3];

    nothing;
end

function updatex2!(oa::ObjectAgent_MAG{T,U}, λ::T) where {T<:Real, U<:Unsigned}
    oa.xk .= oa.zk .- oa.uk;
    oa.xk[oa.act_used] .-= oa.Q\(oa.b1*oa.xk + oa.b0);

    nothing;
end

function updatex!(oa::ObjectAgent_MAG{T,U}, λ::T) where {T<:Real, U<:Unsigned}
    oa.xk .= oa.zk .- oa.uk;

    mul!(oa.r, oa.Gxy', oa.xk); # b0 = fi = Gxy_red'*oa.xk
    oa.r[1] -= oa.Fdes[1];
    oa.r[2] -= oa.Fdes[2];

    LinSolvers.forward_substitution!(oa.R, oa.r)
    LinSolvers.backward_substitution!(oa.R, oa.r)

    mul!(oa.decDir, oa.J, oa.r);

    for i in 1:length(oa.act_used)
        oa.xk[oa.act_used[i]] -= oa.decDir[i];
    end

    nothing;
end

function saturate_control_action(oa::ObjectAgent_MAG{T,U}, x::T) where {T<:Real, U<:Unsigned}
    """
    This method is used to saturate the control aciton.
    """
    min(max(x, 0), 1)
end

function costFun!(oa::ObjectAgent_MAG)
    oa.fvk .= oa.Gxy' * oa.xk .- oa.Fdes;

    nothing
end

function costFun(oa::ObjectAgent_MAG{T, U}, xk::Vector{T}) where {T<:Real,U<:Unsigned}
    return  sum((oa.Gxy' * oa.xk .- oa.Fdes).^2);
end


text:
#DEMO1: Sample Multi-Dimensional (Vectorized) calculation:
#-------------------------------------------------------------------------------

using MDDatasets
using InspectDR
using Colors

#==Constants
===============================================================================#
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)
black = RGB24(0, 0, 0)
grey = RGB24(0.5, 0.5, 0.5)

dfltcolorlist = [red, green, black, blue, grey] #Different color for each
dfltline = line(color=blue, width=3)
redline = line(color=red, width=3)

#==Helper functions
===============================================================================#
#Add DataF1 signal to plot:
function _add(plot, sig::DataF1, args...; line=dfltline, kwargs...)
	wfrm = add(plot, sig.x, sig.y, args...; kwargs...)
	wfrm.line = line
	return wfrm
end

#Add DataRS multi-dimensional signal to plot (1-D sweep only):
#TODO: Use EasyPlot module, once registered
function _add(plot, sig::DataRS, args...; id="", line=dfltline, clist=dfltcolorlist, kwargs...)
	for (data, sweepv, color) in zip(sig.elem, sig.sweep.v, clist)
		wfrm = add(plot, data.x, data.y, args..., id="$id[$sweepv]"; kwargs...)
			wfrm.line = deepcopy(line)
			wfrm.line.color = color
	end
	return
end

#Display plot using "do" syntax:
function displayplot(fn, striplist; title="Plotting variable.y vs variable.x")
	plot = InspectDR.transientplot(striplist, title=title)
	fn(plot)
	display(InspectDR.GtkDisplay(), plot)
end

#==DEMO: Sample Multi-Dimensional (Vectorized) Calculation:
===============================================================================#

#Create (x,y) container pair, and call it "x":
x = DataF1(0:.1:20)
#NOTE: Both x & y coordinates of "x" object initialized as y = x = [supplied range]

#“Extract” maximum x-value from data:
xmax = maximum(x)

#Construct noramlized ramp dataset:
unity_ramp = x/xmax

#Observe x & unity_ramp
#-------------------------------------------------------------------------------
displayplot([:lin, :lin]) do plot
	_add(plot, x, id="x", strip = 1)
	_add(plot, unity_ramp, id="unity_ramp", strip = 2)
	InspectDR.write_png("samplemdcalc_1.png", plot)
end

#Compute cos(kx) & ksinkx = cos'(kx):
coskx = cos((2.5pi/10)*x)
ksinkx = deriv(coskx)

#Compute ramps with different slopes using unity_ramp (previously computed):
#(NOTE: Inner-most sweep, we need to specify leaf element type (DataF1 here))
ramp = fill(DataRS{DataF1}, PSweep("slope", [-1, -0.5, 0, 0.5, 1])) do slope
	return unity_ramp * slope
end
#NOTE: Above expression constructs a multi-dimensional DataRS structure,
#      and fills it with (x,y) values for each of the desired parameter
#      values (the slope).

#Observe coskx, ksinkx & ramp
#-------------------------------------------------------------------------------
displayplot([:lin, :lin]) do plot
	_add(plot, coskx, id="coskx", strip = 1)
	_add(plot, ksinkx, id="ksinkx", line=redline, strip = 1)
	_add(plot, ramp, id="ramp", strip = 2)
	plot.strips[1].yext_full = InspectDR.PExtents1D(min=-1.2, max=1.2)
	plot.strips[2].yext_full = InspectDR.PExtents1D(min=-1.5, max=1.5)
	InspectDR.write_png("samplemdcalc_2.png", plot)
end

#Merge 2 datasets with different # of sweeps (coskx & ramp):
r_cos = coskx+ramp

#Observe r_cos
#-------------------------------------------------------------------------------
displayplot([:lin], title = "r_cos = coskx+ramp") do plot
	_add(plot, r_cos, id="r_cos")
	plot.strips[1].yext_full = InspectDR.PExtents1D(min=-2, max=2)
	InspectDR.write_png("samplemdcalc_3.png", plot)
end

#Shift all ramped sin(x) waveforms to be centered at their mid-points:
midval = value(ramp, x=xmax/2)
#midval = (minimum(ramp) + maximum(ramp)) / 2 #Also possible
c_cos = r_cos - midval #Shift by midval (different for each swept slope of "ramp")

#Observe c_cos
#-------------------------------------------------------------------------------
displayplot([:lin], title = "c_cos = r_cos - value(ramp, x=xmax/2)") do plot
	_add(plot, c_cos, id="c_cos")
	plot.strips[1].yext_full = InspectDR.PExtents1D(min=-1.5, max=1.5)
	InspectDR.write_png("samplemdcalc_4.png", plot)
end

#Last line


text:
module ReverseDiffX

export NotTracked

using LinearAlgebra
using ..ReverseDiff
using StaticArrays
using Distributions
using PDMats
using ..ForwardDiff

using Base.Broadcast: BroadcastStyle, ArrayStyle, Broadcasted, broadcasted
using ..ForwardDiff: Dual
using ..ReverseDiff: SpecialInstruction, value, value!, deriv, track, record!,
                     tape, unseed!, @grad, TrackedReal, TrackedVector,
                     TrackedMatrix, TrackedArray
using ..DistributionsAD: DistributionsAD

import SpecialFunctions, NaNMath
import ..DistributionsAD: turing_chol, symm_turing_chol, _mv_categorical_logpdf, adapt_randn,
    simplex_logpdf
import Base.Broadcast: materialize
import StatsFuns: logsumexp

const TrackedVecOrMat{V,D} = Union{TrackedVector{V,D},TrackedMatrix{V,D}}
const RDBroadcasted{F, T} = Broadcasted{<:Any, <:Any, F, T}

import Random

import Distributions: logpdf,
                      _logpdf,
                      loglikelihood,
                      Gamma,
                      MvNormal,
                      MvLogNormal,
                      Dirichlet,
                      Wishart,
                      InverseWishart,
                      PoissonBinomial,
                      isprobvec

using ..DistributionsAD: TuringPoissonBinomial,
                         TuringMvNormal,
                         TuringMvLogNormal,
                         TuringWishart,
                         TuringInverseWishart,
                         TuringDirichlet,
                         TuringScalMvNormal,
                         TuringDiagMvNormal,
                         TuringDenseMvNormal,
                         VectorOfMultivariate,
                         FillVectorOfMultivariate

include("reversediffx.jl")

adapt_randn(rng::Random.AbstractRNG, x::TrackedArray, dims...) = adapt_randn(rng, value(x), dims...)

function PoissonBinomial(p::TrackedArray{<:Real}; check_args=true)
    return TuringPoissonBinomial(p; check_args = check_args)
end

Gamma(α::TrackedReal, θ::Real; check_args=true) = pgamma(α, θ, check_args = check_args)
Gamma(α::Real, θ::TrackedReal; check_args=true) = pgamma(α, θ, check_args = check_args)
Gamma(α::TrackedReal, θ::TrackedReal; check_args=true) = pgamma(α, θ, check_args = check_args)
pgamma(α, θ; check_args=true) = Gamma(promote(α, θ)..., check_args = check_args)
Gamma(α::T; check_args=true) where {T <: TrackedReal} = Gamma(α, one(T), check_args = check_args)
function Gamma(α::T, θ::T; check_args=true) where {T <: TrackedReal}
    check_args && Distributions.@check_args(Gamma, α > zero(α) && θ > zero(θ))
    return Gamma{T}(α, θ)
end

# Work around to stop TrackedReal of Inf and -Inf from producing NaN in the derivative
function Base.minimum(d::LocationScale{T}) where {T <: TrackedReal}
    if isfinite(minimum(d.ρ))
        return d.μ + d.σ * minimum(d.ρ)
    else
        return convert(T, ReverseDiff.@skip(minimum)(d.ρ))
    end
end
function Base.maximum(d::LocationScale{T}) where {T <: TrackedReal}
    if isfinite(minimum(d.ρ))
        return d.μ + d.σ * maximum(d.ρ)
    else
        return convert(T, ReverseDiff.@skip(maximum)(d.ρ))
    end
end

## MvNormal

for (f, T) in (
    (:_logpdf, :TrackedVector),
    (:logpdf, :TrackedMatrix),
    (:loglikelihood, :TrackedMatrix),
)
    @eval begin
        function $f(d::MvNormal{<:Real,<:PDMats.ScalMat}, x::$T{<:Real})
            return $f(TuringScalMvNormal(d.μ, sqrt(d.Σ.value)), x)
        end
        function $f(d::MvNormal{<:Real,<:PDMats.PDiagMat}, x::$T{<:Real})
            return $f(TuringDiagMvNormal(d.μ, sqrt.(d.Σ.diag)), x)
        end
        function $f(d::MvNormal{<:Real,<:PDMats.PDMat}, x::$T{<:Real})
            return $f(TuringDenseMvNormal(d.μ, d.Σ.chol), x)
        end

        function $f(d::MvLogNormal{<:Real,<:PDMats.ScalMat}, x::$T{<:Real})
            return $f(
                TuringMvLogNormal(TuringScalMvNormal(d.normal.μ, sqrt(d.normal.Σ.value))),
                x,
            )
        end
        function $f(d::MvLogNormal{<:Real,<:PDMats.PDiagMat}, x::$T{<:Real})
            return $f(
                TuringMvLogNormal(TuringDiagMvNormal(d.normal.μ, sqrt.(d.normal.Σ.diag))),
                x,
            )
        end
        function $f(d::MvLogNormal{<:Real,<:PDMats.PDMat}, x::$T{<:Real})
            return $f(
                TuringMvLogNormal(TuringDenseMvNormal(d.normal.μ, d.normal.Σ.chol)),
                x,
            )
        end
    end
end

# zero mean, dense covariance
MvNormal(A::TrackedMatrix) = TuringMvNormal(A)

# zero mean, diagonal covariance
MvNormal(σ::TrackedVector) = TuringMvNormal(σ)

# dense mean, dense covariance
MvNormal(m::AbstractVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::TrackedVector{<:Real}, A::Matrix{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::TrackedVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvNormal(m, A)

# dense mean, diagonal covariance
function MvNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvNormal(m, D)
end
function MvNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvNormal(m, D)
end
function MvNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{T, <:AbstractVector{T}} where {T<:Real},
)
    return TuringMvNormal(m, D)
end

# dense mean, diagonal covariance
MvNormal(m::TrackedVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::AbstractVector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::Vector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::AbstractVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvNormal(m, σ)

# dense mean, constant variance
MvNormal(m::TrackedVector{<:Real}, σ::TrackedReal) = TuringMvNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::Real) = TuringMvNormal(m, σ)
MvNormal(m::AbstractVector{<:Real}, σ::TrackedReal) = TuringMvNormal(m, σ)

# dense mean, constant variance
function MvNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvNormal(m, A)
end
function MvNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvNormal(m, A)
end
function MvNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:Real})
    return TuringMvNormal(m, A)
end

# zero mean,, constant variance
MvNormal(d::Int, σ::TrackedReal) = TuringMvNormal(d, σ)

# zero mean, dense covariance
MvLogNormal(A::TrackedMatrix) = TuringMvLogNormal(TuringMvNormal(A))

# zero mean, diagonal covariance
MvLogNormal(σ::TrackedVector) = TuringMvLogNormal(TuringMvNormal(σ))

# dense mean, dense covariance
MvLogNormal(m::TrackedVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::TrackedVector{<:Real}, A::Matrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::AbstractVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))

# dense mean, diagonal covariance
function MvLogNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end
function MvLogNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end
function MvLogNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{T, <:AbstractVector{T}} where {T<:Real},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end

# dense mean, diagonal covariance
MvLogNormal(m::TrackedVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::TrackedVector{<:Real}, σ::AbstractVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::TrackedVector{<:Real}, σ::Vector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::AbstractVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))

# dense mean, constant variance
function MvLogNormal(m::TrackedVector{<:Real}, σ::TrackedReal)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end
function MvLogNormal(m::TrackedVector{<:Real}, σ::Real)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end
function MvLogNormal(m::AbstractVector{<:Real}, σ::TrackedReal)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end

# dense mean, constant variance
function MvLogNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end
function MvLogNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end
function MvLogNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:Real})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end

# zero mean,, constant variance
MvLogNormal(d::Int, σ::TrackedReal) = TuringMvLogNormal(TuringMvNormal(d, σ))

# Dirichlet

Dirichlet(alpha::AbstractVector{<:TrackedReal}) = TuringDirichlet(alpha)
Dirichlet(d::Integer, alpha::TrackedReal) = TuringDirichlet(d, alpha)

function _logpdf(d::Dirichlet, x::AbstractVector{<:TrackedReal})
    return _logpdf(TuringDirichlet(d), x)
end
function logpdf(d::Dirichlet, x::AbstractMatrix{<:TrackedReal})
    return logpdf(TuringDirichlet(d), x)
end
function loglikelihood(d::Dirichlet, x::AbstractMatrix{<:TrackedReal})
    return loglikelihood(TuringDirichlet(d), x)
end

for func_header in [
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::Real, x::AbstractVector)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractVector)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::Real, x::AbstractVector{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::TrackedReal, x::AbstractVector)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractVector{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::Real, x::AbstractVector{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::TrackedReal, x::AbstractVector{<:TrackedReal})),

    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::Real, x::AbstractMatrix)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractMatrix)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::Real, x::AbstractMatrix{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::TrackedReal, x::AbstractMatrix)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractMatrix{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::Real, x::AbstractMatrix{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::TrackedReal, x::AbstractMatrix{<:TrackedReal})),
]
    @eval $func_header = track(simplex_logpdf, alpha, lmnB, x)
end
@grad function simplex_logpdf(alpha, lmnB, x::AbstractVector)
    _alpha = value(alpha)
    _lmnB = value(lmnB)
    _x = value(x)
    simplex_logpdf(_alpha, _lmnB, _x), Δ -> begin
        (Δ .* log.(_x), -Δ, Δ .* (_alpha .- 1) ./ _x)
    end
end
@grad function simplex_logpdf(alpha, lmnB, x::AbstractMatrix)
    _alpha = value(alpha)
    _lmnB = value(lmnB)
    _x = value(x)
    simplex_logpdf(_alpha, _lmnB, _x), Δ -> begin
        (log.(_x) * Δ, -sum(Δ), ((_alpha .- 1) ./ _x) * Diagonal(Δ))
    end
end

Distributions.Wishart(df::TrackedReal, S::Matrix{<:Real}) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::AbstractMatrix{<:Real}) = TuringWishart(df, S)
Distributions.Wishart(df::Real, S::AbstractMatrix{<:TrackedReal}) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::AbstractMatrix{<:TrackedReal}) = TuringWishart(df, S)
Distributions.Wishart(df::Real, S::TrackedMatrix) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::TrackedMatrix) = TuringWishart(df, S)
Distributions.Wishart(df::Real, S::AbstractPDMat{<:TrackedReal}) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::AbstractPDMat{<:TrackedReal}) = TuringWishart(df, S)

Distributions.InverseWishart(df::TrackedReal, S::Matrix{<:Real}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::AbstractMatrix{<:Real}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::Real, S::AbstractMatrix{<:TrackedReal}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::AbstractMatrix{<:TrackedReal}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::Real, S::TrackedMatrix) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::TrackedMatrix) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::Real, S::AbstractPDMat{<:TrackedReal}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::AbstractPDMat{<:TrackedReal}) = TuringInverseWishart(df, S)

function _logpdf(d::Wishart, X::TrackedMatrix)
    return _logpdf(TuringWishart(d), X)
end
function logpdf(d::Wishart, X::AbstractArray{<:TrackedMatrix})
    return logpdf(TuringWishart(d), X)
end
function loglikelihood(d::Wishart, X::AbstractArray{<:TrackedMatrix})
    return loglikelihood(TuringWishart(d), X)
end

function _logpdf(d::InverseWishart, X::TrackedMatrix)
    return _logpdf(TuringInverseWishart(d), X)
end
function logpdf(d::InverseWishart, X::AbstractArray{<:TrackedMatrix})
    return logpdf(TuringInverseWishart(d), X)
end
function loglikelihood(d::InverseWishart, X::AbstractArray{<:TrackedMatrix})
    return loglikelihood(TuringInverseWishart(d), X)
end

# isprobvec

function isprobvec(p::TrackedArray{<:Real})
    pdata = value(p)
    all(x -> x ≥ zero(x), pdata) && isapprox(sum(pdata), one(eltype(pdata)), atol = 1e-6)
end
function isprobvec(p::SubArray{<:TrackedReal, 1, <:TrackedArray{<:Real}})
    pdata = value(p)
    all(x -> x ≥ zero(x), pdata) && isapprox(sum(pdata), one(eltype(pdata)), atol = 1e-6)
end

_mv_categorical_logpdf(ps::TrackedMatrix, x) = track(_mv_categorical_logpdf, ps, x)
@grad function _mv_categorical_logpdf(ps, x)
    ps_data = value(ps)
    probs = view(ps_data, x, :)
    ps_grad = zero(ps_data)
    sum(log, probs), Δ -> begin
        ps_grad .= 0
        ps_grad[x,:] .= Δ ./ probs
        return (ps_grad, nothing)
    end
end

end


text:
module Mobile3D

using Modia

# Modia equation-based models
include("$(Modia.modelsPath)/AllModels.jl")

const depthMax = 5
const enableVisualization = true
const animationFile = string("Mobile", depthMax, ".json")

const rodLength       = 0.1
const rodDiameter     = 0.05*rodLength
const rodMaterial     = "DryWood"
const rodColor        = VisualMaterial(color="Brown", transparency=0.5)
const rodTranslation1 = [0.0,  rodLength/2, 0.0]
const rodTranslation2 = [0.0, -rodLength/2, 0.0]

const sphereDiameter  = 0.3*rodLength
const sphereMaterial  = "DryWood"
const sphereColor     = VisualMaterial(color="Blue")

const barWidth       = 0.3*rodLength
const barThickness   = 0.3*rodLength
const barMaterial    = "DryWood"
const barColor       = VisualMaterial(color="DeepSkyBlue3", transparency=0.5)
const leaveBarLength = 0.8*rodLength

const revDiameter = rodDiameter
const revLength   = 1.2*barThickness
const revColor    = VisualMaterial(color="Red")

const damping = 0.006

Rod = Model(
    frame0 = Object3D(feature = Solid(shape = Cylinder(axis=2, length=rodLength, diameter=rodDiameter),
                                    solidMaterial  = rodMaterial,
                                    visualMaterial = rodColor)),
    frame1 = Object3D(parent = :frame0, translation = rodTranslation1),
    frame2 = Object3D(parent = :frame0, translation = rodTranslation2)
)
Bar = Model(
    L = missing,
#    frame0 = Object3D(feature = Solid(shape = Beam(axis = 1, length = :L, width=barWidth, thickness=barThickness),
    frame0 = Object3D(feature = Solid(shape = Box(lengthX=:(L + barWidth), lengthY=barWidth, lengthZ=barThickness),
                                      solidMaterial  = barMaterial,
                                      visualMaterial = barColor)),
    frame1 = Object3D(parent = :frame0, translation = :[-L/2, 0.0, 0.0]),
    frame2 = Object3D(parent = :frame0, translation = :[ L/2, 0.0, 0.0])
)
RevoluteWithDamping(;obj1, obj2, axis=3, phi_start=0.0) = Model(
    rev     = RevoluteWithFlange(obj1=obj1, obj2=obj2, axis=axis, phi=Var(init=phi_start)),
    cyl     = Object3D(parent  = obj1,
                    feature = Visual(shape = Cylinder(axis=axis, diameter=revDiameter, length=revLength),
                                        visualMaterial = revColor)),
    damper  = Damper | Map(d=damping),
    fixed   = Fixed,
    connect = :[(damper.flange_b, rev.flange),
                (damper.flange_a, fixed.flange)]
)
function barLength(depth)
    if depth<=2
        return leaveBarLength
    end
    return 2*barLength(depth-1)
end
function createMobile(depth)
    if depth == 1
        Model(
            rod    = Rod,
            sphere = Object3D(parent = :(rod.frame0), translation = rodTranslation2,
                              feature = Solid(shape = Sphere(diameter=sphereDiameter),
                                              solidMaterial  = sphereMaterial,
                                              visualMaterial = sphereColor))
        )
    else
        Model(
            rod  = Rod,
            bar  = Bar | Map(L=barLength(depth)),
            sub1 = createMobile(depth-1),
            sub2 = createMobile(depth-1),
            rev0 = RevoluteWithDamping(obj1=:(rod.frame2), obj2=:(bar.frame0)),
            rev1 = RevoluteWithDamping(obj1=:(bar.frame1), obj2=:(sub1.rod.frame1)),
            rev2 = RevoluteWithDamping(obj1=:(bar.frame2), obj2=:(sub2.rod.frame1))
        )
    end
end
Mobile = Model(
    world      = Object3D(feature=Scene(gravityField=UniformGravityField(g=9.81, n=[0, -1, 0]), enableContactDetection=false, nominalLength=0.2*barLength(depthMax), enableVisualization=enableVisualization, animationFile=animationFile)),
    worldFrame = Object3D(parent = :world, feature = Visual(shape = CoordinateSystem(length=rodLength))),
    top        = createMobile(depthMax),
    rev0       = RevoluteWithDamping(obj1 = :world, obj2 = :(top.rod.frame1), phi_start=0.2)
)

mobile = @instantiateModel(buildModia3D(Mobile), unitless=true, log=false, logDetails=false, logModel=false, logStateSelection=false, logCode=false, logExecution=false, logTiming=false)

stopTime = 5.0
tolerance = 1e-5
requiredFinalStates = [0.011686313703318246, 0.01045668408948015, 0.0019009731577780986, -0.01786066760079631, 0.0038231810505064317, -0.02016203637563375, 0.0126793101749162, 0.01691192217141304, 0.0023271830236842326, -0.017388317974414997, 0.004235202194998209, -0.019458842081769283, 0.02369063170860951, 0.15558539617677464, -0.004928234187014311, -0.24357952222975443, -0.0008598315100012872, -0.25778191880090334, 0.011970016614214521, 0.009770151879848948, 0.0019715145725989034, -0.018428234746089393, 0.003927329453913458, -0.020712258229027464, 0.0130019833269834, 0.01646408641685571, 0.002416123119856374, -0.017925382312598832, 0.004356974507478808, -0.01996955076328113, 0.024253081005718227, 0.15740983132693404, -0.00508053867112354, -0.24798117927445137, -0.0008446712535811513, -0.26259974654104284, -0.04547737419362646, 0.262258261366544, 0.06092051704639535, -0.2992009041652458, 0.062163705985459186, -0.30201328017123014, 0.011774841775337119, 0.008484715801774874, 0.0017813540404277002, -0.018359678091449557, 0.003722833439382662, -0.020559247090548637, 0.012811196185191788, 0.014445725473740444, 0.002212787855986135, -0.018108079678312596, 0.00414068450155657, -0.02006172527918578, 0.024682227696403755, 0.15164227948198014, -0.0064617212778372585, -0.2427675425787477, -0.002449157908043746, -0.2581724674241067, 0.012050676208440915, 0.00785647627890249, 0.0018526658835679355, -0.01890009720306999, 0.003825074776328946, -0.02107627115366446, 0.013126442724564754, 0.014047087346111506, 0.0023025518844593472, -0.018620796297829156, 0.004260695947219754, -0.020542119547756133, 0.025232373917168024, 0.15350921040211213, -0.006601205580033179, -0.24702565773703974, -0.0024230921741892745, -0.2628525384776904, -0.0440032853139086, 0.269398657305908, 0.05921435657256203, -0.3091722470155659, 0.0604662087479324, -0.31187411328295117, 0.07674790919904127, 0.17815099357106337, -0.06512247859139587, -0.1963202500586452, -0.06616717655644311, -0.20502267390766496, 0.04770632153276008, -0.1850427647991364]
simulate!(mobile, stopTime=stopTime, tolerance=tolerance, log=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(mobile, "rev0.rev.phi")

end


text:
for x=10:10:20, y=1:2
    # y rotates than x
    println(x+y)

end

# List comprehension
[println(i) for i in collect(1:5)]
# x rotates then y
[println(x+y) for x in [10, 20], y in [1,2]]


text:
module TestBangBang
using Test

@testset "$file" for file in sort([file for file in readdir(@__DIR__) if match(
    r"^test_.*\.jl$",
    file,
) !== nothing])

    if file == "test_doctest.jl"
        if lowercase(get(ENV, "JULIA_PKGEVAL", "false")) == "true"
            @info "Skipping doctests on PkgEval."
            continue
        elseif !(v"1.6" <= VERSION < v"1.7-")
            @info "Skipping doctests on Julia $VERSION."
            continue
        end
    end

    if file == "test_zygote.jl"
        if lowercase(get(ENV, "JULIA_PKGEVAL", "false")) == "true"
            @info "Skipping tests with Zygote on PkgEval."
            continue
        elseif VERSION < v"1.2-"
            @info "Skip $file for Julia $VERSION to avoid segfault"
            # For example:
            # https://travis-ci.com/JuliaFolds/BangBang.jl/builds/134791701
            # https://travis-ci.com/JuliaFolds/BangBang.jl/builds/134795161
            continue
        elseif VERSION >= v"1.8-"
            @info "Skip $file for Julia $VERSION"
            # For example:
            continue
        end
    end

    include(file)
end

end  # module


text:
export getNode

function getNode(configDict, fg, requestDict)::Dict{String, Any}
    error("Please use getVariable or getFactor, this function is deprecated")
    # TODO: Build a cleaner contract to return this value.
    vert = RoME.getVert(fg, Symbol(requestDict["payload"]))
    return JSON.parse(JSON.json(vert))
end


text:
using Documenter, DocumenterMarkdown
using PhyloNetworks
using QuartetNetworkGoodnessFit
DocMeta.setdocmeta!(QuartetNetworkGoodnessFit, :DocTestSetup, :(using QuartetNetworkGoodnessFit); recursive=true)

makedocs(
    sitename = "QuartetNetwork GoF.jl",
    authors = "Cécile Ané and Ruoyi Cai",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"), # easier local build
    modules = [QuartetNetworkGoodnessFit],
    pages = [
        "home" => "index.md",
        "manual" => [
            "goodness of fit" => "man/gof.md",
        ],
        "library" => [
            "public" => "lib/public.md",
            "internal" => "lib/internal.md",
        ]
    ],
)

deploydocs(
    repo = "github.com/cecileane/QuartetNetworkGoodnessFit.jl.git",
    push_preview = true,
)


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libpng16"], :libpng),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/SimonDanisch/LibpngBuilder/releases/download/v1.6.31"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/libpng.v1.0.0.aarch64-linux-gnu.tar.gz", "066248a17adcbd9533a3abb61aeb92d96f1341e1b588ee938ca93aece66afeee"),
    Linux(:aarch64, :musl) => ("$bin_prefix/libpng.v1.0.0.aarch64-linux-musl.tar.gz", "f1f85dcb141b508301db19e9a11ef43a999b39b1a4ced29e64c8f52700eab43b"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/libpng.v1.0.0.arm-linux-gnueabihf.tar.gz", "28be4910b6ef152d604b432b2aef9617774e306b2744c312c7499196a660950c"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/libpng.v1.0.0.arm-linux-musleabihf.tar.gz", "0809cbce8c12a3f9632718a827ed8820225374553eb3eda92a3e8c56b033b77b"),
    Linux(:i686, :glibc) => ("$bin_prefix/libpng.v1.0.0.i686-linux-gnu.tar.gz", "c6937fcc23a15c3eadfbf13a1fa7aa8046e5e82192d0e79378d1d57e25c2a7c0"),
    Linux(:i686, :musl) => ("$bin_prefix/libpng.v1.0.0.i686-linux-musl.tar.gz", "50a7c43df57ef32d0ab61f5059f3bccf423e18d4d9e86d1a084174bbe2b62ffb"),
    Windows(:i686) => ("$bin_prefix/libpng.v1.0.0.i686-w64-mingw32.tar.gz", "7adfe9b2e2f91798d43bc98f79645b69cb611fd5221fa803ce8ba530ccf9bcfa"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/libpng.v1.0.0.powerpc64le-linux-gnu.tar.gz", "90430fa0a2858db0ffc87981e14c5864171389e6af50be9cf7c9ed5a85713e5c"),
    MacOS(:x86_64) => ("$bin_prefix/libpng.v1.0.0.x86_64-apple-darwin14.tar.gz", "53a530ca9251761cd3e0123fee3ff470a473097e7f112038e4c5b7ca7aa4eae9"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/libpng.v1.0.0.x86_64-linux-gnu.tar.gz", "eee0ffe168e0ba11129a544f666bc9840036a0761234f36e15421c9004c5be55"),
    Linux(:x86_64, :musl) => ("$bin_prefix/libpng.v1.0.0.x86_64-linux-musl.tar.gz", "458cb397e183772a87ca7bc7fbb309cdc082a74a70460ce3b90027f1565d5f40"),
    FreeBSD(:x86_64) => ("$bin_prefix/libpng.v1.0.0.x86_64-unknown-freebsd11.1.tar.gz", "afe7b58dc391ff1f2e982c5a734d33fc942ba611349d79f917d5c0647cc1148a"),
    Windows(:x86_64) => ("$bin_prefix/libpng.v1.0.0.x86_64-w64-mingw32.tar.gz", "ff95e2f2f7d5c410c5374e7e86ab443bb58ff132bf418ab0db99c3385cb0bee2"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps_png.jl"), products)