text:
module MPIArrays

export MPIArray, localindices, getblock, getblock!, putblock!, allocate, forlocalpart, forlocalpart!, free, redistribute, redistribute!, sync, GlobalBlock, GhostedBlock, getglobal, globaltolocal, globalids

using MPI
import LinearAlgebra

"""
Store the distribution of the array indices over the different partitions.
This class forces a continuous, ordered distribution without overlap

    ContinuousPartitioning(partition_sizes...)

Construct a distribution using the number of elements per partition in each direction, e.g:

```
p = ContinuousPartitioning([2,5,3], [2,3])
```

will construct a distribution containing 6 partitions of 2, 5 and 3 rows and 2 and 3 columns.

"""
struct ContinuousPartitioning{N} <: AbstractArray{Int,N}
    ranks::LinearIndices{N,NTuple{N,Base.OneTo{Int}}}
    index_starts::NTuple{N,Vector{Int}}
    index_ends::NTuple{N,Vector{Int}}

    function ContinuousPartitioning(partition_sizes::Vararg{Any,N}) where {N}
        index_starts = Vector{Int}.(undef,length.(partition_sizes))
        index_ends = Vector{Int}.(undef,length.(partition_sizes))
        for (idxstart,idxend,nb_elems_dist) in zip(index_starts,index_ends,partition_sizes)
            currentstart = 1
            currentend = 0
            for i in eachindex(idxstart)
                currentend += nb_elems_dist[i]
                idxstart[i] = currentstart
                idxend[i] = currentend
                currentstart += nb_elems_dist[i]
            end
        end
        ranks = LinearIndices(length.(partition_sizes))
        return new{N}(ranks, index_starts, index_ends)
    end
end

Base.IndexStyle(::Type{ContinuousPartitioning{N}}) where {N} = IndexCartesian()
Base.size(p::ContinuousPartitioning) = length.(p.index_starts)
@inline function Base.getindex(p::ContinuousPartitioning{N}, I::Vararg{Int, N}) where {N}
    return UnitRange.(getindex.(p.index_starts,I), getindex.(p.index_ends,I))
end

function partition_sizes(p::ContinuousPartitioning)
    result = (p.index_ends .- p.index_starts)
    for v in result
        v .+= 1
    end
    return result
end

"""
  (private method)

Get the rank and local 0-based index
"""
function local_index(p::ContinuousPartitioning, I::NTuple{N,Int}) where {N}
    proc_indices = searchsortedfirst.(p.index_ends, I)
    lininds = LinearIndices(Base.Slice.(p[proc_indices...]))
    return (p.ranks[proc_indices...]-1, lininds[I...] - first(lininds))
end

# Evenly distribute nb_elems over parts partitions
function distribute(nb_elems, parts)
    local_len = nb_elems ÷ parts
    remainder = nb_elems % parts
    return [p <= remainder ? local_len+1 : local_len for p in 1:parts]
end

mutable struct MPIArray{T,N} <: AbstractArray{T,N}
    sizes::NTuple{N,Int}
    localarray::Array{T,N}
    partitioning::ContinuousPartitioning{N}
    comm::MPI.Comm
    win::MPI.Win
    myrank::Int

    function MPIArray{T}(comm::MPI.Comm, partition_sizes::Vararg{AbstractVector{<:Integer},N}) where {T,N}
        nb_procs = MPI.Comm_size(comm)
        rank = MPI.Comm_rank(comm)
        partitioning = ContinuousPartitioning(partition_sizes...)

        localarray = Array{T}(undef,length.(partitioning[rank+1]))
        win = MPI.Win_create(localarray, comm)
        sizes = sum.(partition_sizes)
        return new{T,N}(sizes, localarray, partitioning, comm, win, rank)
    end

    MPIArray{T}(sizes::Vararg{<:Integer,N}) where {T,N} = MPIArray{T}(MPI.COMM_WORLD, (MPI.Comm_size(MPI.COMM_WORLD), ones(Int,N-1)...), sizes...)
    MPIArray{T}(comm::MPI.Comm, partitions::NTuple{N,<:Integer}, sizes::Vararg{<:Integer,N}) where {T,N} = MPIArray{T}(comm, distribute.(sizes, partitions)...)

    function MPIArray(comm::MPI.Comm, localarray::Array{T,N}, nb_partitions::Vararg{<:Integer,N}) where {T,N}
        nb_procs = MPI.Comm_size(comm)
        rank = MPI.Comm_rank(comm)

        partition_size_array = reshape(MPI.Allgather(size(localarray), comm), Int.(nb_partitions)...)

        partition_sizes = ntuple(N) do dim
            idx = ntuple(i -> i == dim ? Colon() : 1,N)
            return getindex.(partition_size_array[idx...],dim)
        end

        win = MPI.Win_create(localarray, comm)
        result = new{T,N}(sum.(partition_sizes), localarray, ContinuousPartitioning(partition_sizes...), comm, win, rank)
        return result
    end

    MPIArray(localarray::Array{T,N}, nb_partitions::Vararg{<:Integer,N}) where {T,N} = MPIArray(MPI.COMM_WORLD, localarray, nb_partitions...)
end

Base.IndexStyle(::Type{MPIArray{T,N}}) where {T,N} = IndexCartesian()

Base.size(a::MPIArray) = a.sizes

# Individual element access. WARNING: this is slow
function Base.getindex(a::MPIArray{T,N}, I::Vararg{Int, N}) where {T,N}
    (target_rank, locind) = local_index(a.partitioning,I)

    result = Ref{T}()
    MPI.Win_lock(MPI.LOCK_SHARED, target_rank, 0, a.win)
    if target_rank == a.myrank
        result[] = a.localarray[locind+1]
    else
        MPI.Get(result, 1, target_rank, locind, a.win)
    end
    MPI.Win_unlock(target_rank, a.win)
    return result[]
end

# Individual element setting. WARNING: this is slow
function Base.setindex!(a::MPIArray{T,N}, v, I::Vararg{Int, N}) where {T,N}
    (target_rank, locind) = local_index(a.partitioning,I)

    MPI.Win_lock(MPI.LOCK_EXCLUSIVE, target_rank, 0, a.win)
    if target_rank == a.myrank
        a.localarray[locind+1] = v
    else
        result = Ref{T}(v)
        MPI.Put(result, 1, target_rank, locind, a.win)
    end
    MPI.Win_unlock(target_rank, a.win)
end

"""
    sync(a::MPIArray)

Collective call, making sure all operations modifying any part of the array are finished when it completes
"""
sync(a::MPIArray, ::Vararg{MPIArray, N}) where N = MPI.Barrier(a.comm)

function Base.similar(a::MPIArray, ::Type{T}, dims::NTuple{N,Int}) where {T,N}
    old_partition_sizes = partition_sizes(a.partitioning)
    old_dims = size(a)
    new_partition_sizes = Vector{Int}[]
    remaining_nb_partitons = prod(length.(old_partition_sizes))
    for i in eachindex(dims)
        if i <= length(old_dims)
            if dims[i] == old_dims[i]
                push!(new_partition_sizes, old_partition_sizes[i])
            else
                push!(new_partition_sizes, distribute(dims[i], length(old_partition_sizes[i])))
            end
        elseif remaining_nb_partitons != 1
            push!(new_partition_sizes, distribute(dims[i], remaining_nb_partitons))
        else
            push!(new_partition_sizes, [dims[i]])
        end
        @assert remaining_nb_partitons % length(last(new_partition_sizes)) == 0
        remaining_nb_partitons ÷= length(last(new_partition_sizes))
    end
    if remaining_nb_partitons > 1
        remaining_nb_partitons *= length(last(new_partition_sizes))
        new_partition_sizes[end] = distribute(dims[end], remaining_nb_partitons)
        remaining_nb_partitons ÷= length(last(new_partition_sizes))
    end
    @assert remaining_nb_partitons == 1
    return MPIArray{T}(a.comm, new_partition_sizes...)
end

function Base.filter(f,a::MPIArray)
    error("filter is only supported on 1D MPIArrays")
end

function Base.filter(f,a::MPIArray{T,1}) where T
    return MPIArray(forlocalpart(v -> filter(f,v), a), length(a.partitioning))
end

function Base.filter!(f,a::MPIArray)
    error("filter is only supported on 1D MPIArrays")
end

function copy_into!(dest::MPIArray{T,N}, src::MPIArray{T,N}) where {T,N}
    free(dest)
    dest.sizes = src.sizes
    dest.localarray = src.localarray
    dest.partitioning = src.partitioning
    dest.comm = src.comm
    dest.win = src.win
    dest.myrank = src.myrank
    return dest
end

Base.filter!(f,a::MPIArray{T,1}) where T = copy_into!(a, filter(f,a))

function redistribute(a::MPIArray{T,N}, partition_sizes::Vararg{Any,N}) where {T,N}
    rank = MPI.Comm_rank(a.comm)
    @assert prod(length.(partition_sizes)) == MPI.Comm_size(a.comm)
    partitioning = ContinuousPartitioning(partition_sizes...)
    localarray = getblock(a[partitioning[rank+1]...])
    return MPIArray(a.comm, localarray, length.(partition_sizes)...)
end

function redistribute(a::MPIArray{T,N}, nb_parts::Vararg{Int,N}) where {T,N}
    return redistribute(a, distribute.(size(a), nb_parts)...)
end

function redistribute(a::MPIArray)
    return redistribute(a, size(a.partitioning)...)
end

redistribute!(a::MPIArray{T,N}, partition_sizes::Vararg{Any,N})  where {T,N} = copy_into!(a, redistribute(a, partition_sizes...))
redistribute!(a::MPIArray) = redistribute!(a, size(a.partitioning)...)

function LinearAlgebra.mul!(y::MPIArray{T,1}, A::MPIArray{T,2}, b::MPIArray{T,1}) where {T}
    forlocalpart!(y) do ly
        fill!(ly,zero(T))
    end
    sync(y)
    (rowrng,colrng) = localindices(A)
    my_b = getblock(b[colrng])
    yblock = y[rowrng]
    my_y = allocate(yblock)
    forlocalpart(A) do my_A
        LinearAlgebra.mul!(my_y,my_A,my_b)
    end
    putblock!(my_y,yblock,+)
    sync(y)
    return y
end

"""
    localindices(a::MPIArray, rank::Integer)

Get the local index range (expressed in terms of global indices) of the given rank
"""
@inline localindices(a::MPIArray, rank::Integer=a.myrank) = a.partitioning[rank+1]

"""
    forlocalpart(f, a::MPIArray)

Execute the function f on the part of a owned by the current rank. It is assumed f does not modify the local part.
"""
function forlocalpart(f, a::MPIArray)
    MPI.Win_lock(MPI.LOCK_SHARED, a.myrank, 0, a.win)
    result = f(a.localarray)
    MPI.Win_unlock(a.myrank, a.win)
    return result
end

"""
    forlocalpart(f, a::MPIArray)

Execute the function f on the part of a owned by the current rank. The local part may be modified by f.
"""
function forlocalpart!(f, a::MPIArray)
    MPI.Win_lock(MPI.LOCK_EXCLUSIVE, a.myrank, 0, a.win)
    result = f(a.localarray)
    MPI.Win_unlock(a.myrank, a.win)
    return result
end

function linear_ranges(indexblock)
    cr = CartesianIndices(axes(indexblock)[2:end])
    result = Vector{UnitRange{Int}}(undef,length(cr))

    for (i,carti) in enumerate(cr)
        linrange = indexblock[:,carti]
        result[i] = linrange[1]:linrange[end]
    end
    return result
end

struct Block{T,N}
    array::MPIArray{T,N}
    ranges::NTuple{N,UnitRange{Int}}
    targetrankindices::CartesianIndices{N,NTuple{N,UnitRange{Int}}}

    function Block(a::MPIArray{T,N}, ranges::Vararg{AbstractRange{Int},N}) where {T,N}
        start_indices =  searchsortedfirst.(a.partitioning.index_ends, first.(ranges))
        end_indices = searchsortedfirst.(a.partitioning.index_ends, last.(ranges))
        targetrankindices = CartesianIndices(UnitRange.(start_indices, end_indices))

        return new{T,N}(a, UnitRange.(ranges), targetrankindices)
    end
end

Base.getindex(a::MPIArray{T,N}, I::Vararg{UnitRange{Int},N}) where {T,N} = Block(a,I...)
Base.getindex(a::MPIArray{T,N}, I::Vararg{Colon,N}) where {T,N} = Block(a,axes(a)...)

"""
Allocate a local array with the size to fit the block
"""
allocate(b::Block{T,N}) where {T,N} = Array{T}(undef,length.(b.ranges))

function blockloop(a::AbstractArray{T,N}, b::Block{T,N}, locktype::MPI.LockType, localfun, mpifun) where {T,N}
    for rankindex in b.targetrankindices
        r = b.array.partitioning.ranks[rankindex] - 1
        localinds = localindices(b.array,r)
        globalrange = b.ranges .∩ localinds
        a_range = map((g,r) -> g .- first(r) .+ 1, globalrange, b.ranges)
        b_range = map((g,r) -> g .- first(r) .+ 1, globalrange, localinds)
        MPI.Win_lock(locktype, r, 0, b.array.win)
        if r == b.array.myrank
            localfun(a, a_range, b.array.localarray, b_range)
        else
            alin = LinearIndices(a)
            blin = LinearIndices(length.(localinds))
            for (ai,bi) in zip(CartesianIndices(a_range[2:end]),CartesianIndices(b_range[2:end]))
                a_linear_range = alin[a_range[1][1],ai]:alin[a_range[1][end],ai]
                b_begin = blin[b_range[1][1],bi] - 1
                range_len = length(a_linear_range)
                @assert b_begin + range_len <= prod(length.(localinds))
                mpifun(pointer(a,first(a_linear_range)), range_len, r, b_begin, b.array.win)
            end
        end
        MPI.Win_unlock(r, b.array.win)
    end
end

"""
Get the (possibly remotely stored) elements of the block and store them in a
"""
function getblock!(a::AbstractArray{T,N}, b::Block{T,N}) where {T,N}
    localfun = function(local_a, a_range, local_b, b_range)
        local_a[a_range...] .= local_b[b_range...]
    end
    blockloop(a, b, MPI.LOCK_SHARED, localfun, MPI.Get)
end

"""
Allocate a matrix to store the data referred to by block and copy all elements from the global array to it
"""
function getblock(b::Block{T,N}) where {T,N}
    a = allocate(b)
    getblock!(a,b)
    return a
end

@inline _no_op(a,b) = b
_mpi_put_with_op(origin, count, rank, disp, win, op) = throw(ErrorException("Unsupported op $op"))
_mpi_put_with_op(origin, count, rank, disp, win, op::typeof(_no_op)) = MPI.Put(origin, count, rank, disp, win)
_mpi_put_with_op(origin, count, rank, disp, win, op::typeof(+)) = MPI.Accumulate(origin, count, rank, disp, MPI.SUM, win)

"""
Set the (possibly remotely stored) elements of the block and store them in a
"""
function putblock!(a::AbstractArray{T,N}, b::Block{T,N}, op::Function=_no_op) where {T,N}
    localfun = function(local_a, a_range, local_b, b_range)
        local_b[b_range...] .= op.(local_b[b_range...], local_a[a_range...])
    end
    mpifun = function(origin, count, target_rank, target_disp, win)
        _mpi_put_with_op(origin, count, target_rank, target_disp, win, op)
    end
    blockloop(a, b, MPI.LOCK_EXCLUSIVE, localfun, mpifun)
end

function free(a::MPIArray{T,N}) where {T,N}
    sync(a)
    MPI.free(a.win)
end

using CustomUnitRanges: filename_for_urange
include(filename_for_urange)

struct GlobalBlock{T,N} <: AbstractArray{T,N}
    array::Array{T,N}
    block::Block{T,N}
end

Base.IndexStyle(::Type{GlobalBlock{T,N}}) where {T,N} = IndexCartesian()
Base.axes(gb::GlobalBlock) = URange.(first.(gb.block.ranges), last.(gb.block.ranges))
@inline _convert_idx(rng, I) = I .- first.(rng) .+ 1
Base.getindex(gb::GlobalBlock{T,N}, I::Vararg{Int, N}) where {T,N} = gb.array[_convert_idx(gb.block.ranges, I)...]
Base.setindex!(gb::GlobalBlock{T,N}, value, I::Vararg{Int, N}) where {T,N} = (gb.array[_convert_idx(gb.block.ranges, I)...] = value)
Base.size(gb::GlobalBlock) = last.(gb.block.ranges), first.(gb.block.ranges)

"""
Read-only block, providing access to the local data and an arbitrary number of off-processer entries.
This is a vector because the represented region is not necessarily square
"""
struct GhostedBlock{T,N} <: AbstractArray{T,1}
    array::MPIArray{T,N}
    globaltolocal::Dict{Int, Int}
    localtoglobal::Vector{Int}
    ghostvalues::Vector{T}

    GhostedBlock(a::MPIArray{T,N}) where {T,N} = new{T,N}(a, Dict{Int,Int}(), Vector{Int}(), Vector{T}())
end

"""
Construct a sorted GhostedBlock that contains all the given gids
The passed list of gids may contain duplicates
"""
function GhostedBlock(a::MPIArray{T,1}, gids) where T
    gb = GhostedBlock(a)
    locinds = localindices(a)[1]
    nb_unique = 0
    for gid in gids
        if !(gid ∈ locinds || haskey(gb.globaltolocal, gid))
            nb_unique += 1
            gb.globaltolocal[gid] = nb_unique
        end
    end
    resize!(gb.localtoglobal, nb_unique)
    resize!(gb.ghostvalues, nb_unique)
    for (gid,i) in gb.globaltolocal
        gb.localtoglobal[i] = gid
    end
    sort!(gb)
    sync(gb)
    return gb
end

Base.IndexStyle(::Type{GhostedBlock{T,N}}) where {T,N} = IndexLinear()
Base.size(b::GhostedBlock) = (length(b.array.localarray) + length(b.localtoglobal),)
function Base.getindex(b::GhostedBlock, i)
    loclen = length(b.array.localarray)
    if i <= loclen
        return b.array.localarray[i]
    end
    return b.ghostvalues[i-loclen]
end

"""
Convert a global index to a local linear index into the union of the local and ghosted nodes
"""
globaltolocal(b::GhostedBlock{T,N}, gid::Integer) where {T,N} = globaltolocal(b, CartesianIndices(b.array)[gid])
function globaltolocal(b::GhostedBlock{T,N}, I::CartesianIndex{N}) where {T,N}
    locinds = localindices(b.array)
    if all(I.I .∈ locinds)
        return LinearIndices(b.array.localarray)[_convert_idx(locinds, I.I)...]
    end
    return length(b.array.localarray) + b.globaltolocal[LinearIndices(b.array)[I]]
end

function globaltolocal(b::GhostedBlock{T,1}, gid::Integer) where {T}
    locinds = localindices(b.array)[1]
    if gid ∈ locinds
        return gid - first(locinds) + 1
    end
    return length(b.array.localarray) + b.globaltolocal[gid]
end

"""
Linear list of all the global IDs referred by the GhostedBlock, starting with the local nodes and with the ghosts at the end
"""
function globalids(b::GhostedBlock)
    local_len = length(b.array.localarray)
    result = Vector{Int}(undef, local_len + length(b.localtoglobal))
    li = LinearIndices(b.array)
    for (i,I) in enumerate(CartesianIndices(localindices(b.array)))
        result[i] = li[I]
    end
    copyto!(result, local_len+1, b.localtoglobal, 1, length(b.localtoglobal))
    return result
end

"""
Get a value using a global index. Returns an error if I is not part of the local array or the ghosted nodes
"""
getglobal(b::GhostedBlock{T,N}, I::Vararg{Integer,N}) where {T,N} = b[globaltolocal(b, CartesianIndex(I...))]
getglobal(b::GhostedBlock{T,N}, I::CartesianIndex{N}) where {T,N} = getglobal(b, I.I...)
getglobal(b::GhostedBlock{T,1}, i::Integer) where {T} = b[globaltolocal(b, i)]

"""
Register the given index as a ghost, if it is not a local node or already registered
"""
function Base.push!(b::GhostedBlock{T,N}, I::Vararg{Int,N}) where{T,N}
    gid = LinearIndices(b.array)[I...]
    if all(I .∈ localindices(b.array)) || gid ∈ keys(b.globaltolocal)
        return b
    end

    push!(b.localtoglobal, gid)
    newidx = length(b.localtoglobal)
    resize!(b.ghostvalues, newidx)
    b.globaltolocal[gid] = newidx

    return b
end
Base.push!(b::GhostedBlock{T,N}, I::CartesianIndex{N}) where{T,N} = push!(b, I.I...)

"""
Sort the indices, so that they are contiguous per rank. Does not sort the data
"""
function Base.sort!(b::GhostedBlock)
    sort!(b.localtoglobal)
    for (i,gid) in enumerate(b.localtoglobal)
        @assert gid ∈ keys(b.globaltolocal)
        b.globaltolocal[gid] = i
    end
end

"""
Fetch the local data
"""
function sync(b::GhostedBlock)
    tocart = CartesianIndices(b.array)
    i = 1
    nb_ghosts = length(b.localtoglobal)
    if nb_ghosts == 0
        return
    end
    (newrank, lid) = local_index(b.array.partitioning, tocart[b.localtoglobal[i]].I)
    while i <= nb_ghosts
        lockedrank = newrank
        MPI.Win_lock(MPI.LOCK_SHARED, newrank, 0, b.array.win)
        while newrank == lockedrank && i <= nb_ghosts
            saved_i = i
            rngstart = lid
            rngend = lid
            i += 1
            while i <= nb_ghosts
                (newrank, lid) = local_index(b.array.partitioning, tocart[b.localtoglobal[i]].I)
                if (lid - rngend) == 1 && newrank == lockedrank
                    rngend = lid
                    i += 1
                else
                    break
                end
            end
            MPI.Get(pointer(b.ghostvalues,saved_i), rngend-rngstart+1, lockedrank, rngstart, b.array.win)
        end
        MPI.Win_unlock(lockedrank, b.array.win)
    end
end

end # module


text:
function henrici()
#
# M. Petkovic testing polynomials, page 123
#
    p = poly([-4.1;-3.8;-2.05;-1.85;1.95;2.15;3.9;4.05])
    z = [[-4.1;-3.8;-2.05;-1.85;1.95;2.15;3.9;4.05] ones(8)];

    p, PolyZeros(z)
end

ext: jl


text:
function convert_arguments(P::PlotFunc, d::KernelDensity.UnivariateKDE)
    ptype = plottype(P, Lines) # choose the more concrete one
    to_plotspec(ptype, convert_arguments(ptype, d.x, d.density))
end

function convert_arguments(::Type{<:Poly}, d::KernelDensity.UnivariateKDE)
    points = Vector{Point2f0}(undef, length(d.x) + 2)
    points[1] = Point2f0(d.x[1], 0)
    points[2:end-1] .= Point2f0.(d.x, d.density)
    points[end] = Point2f0(d.x[end], 0)
    (points,)
end

function convert_arguments(P::PlotFunc, d::KernelDensity.BivariateKDE)
    ptype = plottype(P, Heatmap)
    to_plotspec(ptype, convert_arguments(ptype, d.x, d.y, d.density))
end

"""
    density(values; npoints = 200, offset = 0.0, direction = :x)

Plot a kernel density estimate of `values`.
`npoints` controls the resolution of the estimate, the baseline can be
shifted with `offset` and the `direction` set to :x or :y.
`bandwidth` and `boundary` are determined automatically by default.

`color` is usually set to a single color, but can also be set to `:x` or
`:y` to color with a gradient. If you use `:y` when direction = `:x` (or vice versa),
note that only 2-element colormaps can work correctly.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Density) do scene
    Theme(
        color = theme(scene, :patchcolor),
        colormap = theme(scene, :colormap),
        colorrange = Makie.automatic,
        strokecolor = theme(scene, :patchstrokecolor),
        strokewidth = theme(scene, :patchstrokewidth),
        linestyle = nothing,
        strokearound = false,
        npoints = 200,
        offset = 0.0,
        direction = :x,
        boundary = automatic,
        bandwidth = automatic,
        cycle = [:color => :patchcolor],
        inspectable = theme(scene, :inspectable)
    )
end

function plot!(plot::Density{<:Tuple{<:AbstractVector}})
    x = plot[1]

    lowerupper = lift(x, plot.direction, plot.boundary, plot.offset,
        plot.npoints, plot.bandwidth) do x, dir, bound, offs, n, bw

        k = KernelDensity.kde(x;
            npoints = n,
            (bound === automatic ? NamedTuple() : (boundary = bound,))...,
            (bw === automatic ? NamedTuple() : (bandwidth = bw,))...,
        )

        if dir === :x
            lowerv = Point2f0.(k.x, offs)
            upperv = Point2f0.(k.x, offs .+ k.density)
        elseif dir === :y
            lowerv = Point2f0.(offs, k.x)
            upperv = Point2f0.(offs .+ k.density, k.x)
        else
            error("Invalid direction $dir, only :x or :y allowed")
        end
        (lowerv, upperv)
    end

    linepoints = lift(lowerupper, plot.strokearound) do lu, sa
        if sa
            ps = copy(lu[2])
            push!(ps, lu[1][end])
            push!(ps, lu[1][1])
            push!(ps, lu[1][2])
            ps
        else
            lu[2]
        end
    end

    lower = Node(Point2f0[])
    upper = Node(Point2f0[])

    on(lowerupper) do (l, u)
        lower.val = l
        upper[] = u
    end
    notify(lowerupper)

    colorobs = lift(Any, plot.color, lowerupper, plot.direction) do c, lu, dir
        if (dir == :x && c == :x) || (dir == :y && c == :y)
            dim = dir == :x ? 1 : 2
            [l[dim] for l in lu[1]]
        elseif (dir == :y && c == :x) || (dir == :x && c == :y)
            o = Float32(plot.offset[])
            dim = dir == :x ? 2 : 1
            vcat([l[dim] - o for l in lu[1]], [l[dim] - o for l in lu[2]])
        else
            c
        end
    end

    band!(plot, lower, upper, color = colorobs, colormap = plot.colormap,
        colorrange = plot.colorrange, inspectable = plot.inspectable)
    l = lines!(plot, linepoints, color = plot.strokecolor,
        linestyle = plot.linestyle, linewidth = plot.strokewidth,
        inspectable = plot.inspectable)
    plot
end


text:
isinteractive() && include("common.jl")

####
#### utilities
####

"""
    $SIGNATURES

Recursive comparison by fields; types and values should match by `==`.
"""
function ≂(x::T, y::T) where T
    fn = fieldnames(T)
    isempty(fn) ? x == y : all(getfield(x, f) ≂ getfield(y, f) for f in fn)
end

≂(x, y) = x == y

"""
Type for testing `≂`.
"""
struct Foo{T}
    a::T
    b::T
end

@testset "≂ comparisons" begin
    @test 1 ≂ 1
    @test 1 ≂ 1.0
    @test [1,2] ≂ [1,2]
    @test [1.0,2.0] ≂ [1,2]
    @test !(1 ≂ 2)
    @test Foo(1,2) ≂ Foo(1,2)
    @test !(Foo(1,2) ≂ Foo(1,3))
    @test !(Foo{Any}(1,2) ≂ Foo(1,2))
end

###
### random booleans
###

@testset "random booleans" begin
    for prob in (1:9) ./ 10
        logprob = log(prob)
        @test abs(mean(rand_bool_logprob(RNG, logprob) for _ in 1:10000) - prob) ≤ 0.01
    end

    # these operations don't call the RNG, this is checked
    RNG′ = copy(RNG)
    @test all(rand_bool_logprob(RNG, 0) for _ in 1:10000)
    @test all(rand_bool_logprob(RNG, 10) for _ in 1:10000)
    @test rand(RNG′) == rand(RNG)
end

###
### test turn statistics
###

@testset "low-level turn statistics" begin
    trajectory = TrajectoryNUTS(nothing, 0, 1, -1000, Val{:generalized})
    p = ones(3)                 # unit vector
    c = 0.1                     # a constant, just for consistency checking of combination
    # turn statistics constructed so that τ₁ + τ₂ won't be turning, τ₁ + τ₃ will be
    τ₁ = GeneralizedTurnStatistic(p, p .- c, p, p .- c, p)
    τ₂ = GeneralizedTurnStatistic(3 .* p, 3 .* p .+ c, 3 .* p, 3 .* p .+ c, 3 .* p)
    τ₃ = GeneralizedTurnStatistic(2 .* p, 2 .* p .+ c, 2 .* p, 2 .* p .+ c, -2 .* p)
    τ = combine_turn_statistics(trajectory, τ₁, τ₂)
    # test mechanics of combination
    @test τ ≂ GeneralizedTurnStatistic(τ₁.p₋, τ₁.p♯₋, τ₂.p₊, τ₂.p♯₊, τ₁.ρ .+ τ₂.ρ)
    # test non-turning
    @test !is_turning(trajectory, τ)
    # test turning
    @test is_turning(trajectory, combine_turn_statistics(trajectory, τ₁, τ₃))
end

@testset "low-level visited statistics" begin
    trajectory = TrajectoryNUTS(nothing, 0, 1, -1000, Val{:generalized})
    vs(p, is_initial = false) = leaf_acceptance_statistic(log(p), is_initial)
    x = vs(0.3)
    @test acceptance_rate(x) ≈ 0.3
    y = vs(0.6)
    @test acceptance_rate(y) ≈ 0.6
    x0 = vs(10, true)    # initial node, does not count
    z = reduce((x, y) -> combine_visited_statistics(trajectory, x, y),
               [x, x, y, x0])
    @test acceptance_rate(z) ≈ 0.4
end

# define a distribution which is divergent everywhere except at 0
struct AlwaysDivergentTest
    K::Int
end

function LogDensityProblems.capabilities(::Type{AlwaysDivergentTest})
    LogDensityProblems.LogDensityOrder{1}()
end
LogDensityProblems.dimension(d::AlwaysDivergentTest) = d.K
function LogDensityProblems.logdensity_and_gradient(d::AlwaysDivergentTest, x)
    if all(iszero.(x))
        0.0, zeros(length(x))
    else
        -Inf, nothing
    end
end

@testset "unconditional divergence" begin

    # test NUTS sampler where all movements are divergent
    K = 3
    ℓ = AlwaysDivergentTest(K)
    Q, tree_statistics = sample_tree(RNG, NUTS(), Hamiltonian(GaussianKineticEnergy(K), ℓ),
                                     evaluate_ℓ(ℓ, zeros(K)), 1.0)
    @test is_divergent(tree_statistics.termination)
    @test iszero(tree_statistics.acceptance_rate)
    @test iszero(tree_statistics.depth)
    @test tree_statistics.steps == 1
end

@testset "normal NUTS HMC transition mean and cov" begin
    # A test for sample_tree with a fixed ϵ and κ, which is perfectly adapted and should
    # provide excellent mixing
    for _ in 1:10
        K = rand(2:8)
        N = 10000
        μ = randn(K)
        Σ = rand_Σ(K)
        L = cholesky(Σ).L
        ℓ = multivariate_normal(μ, L)
        Q = evaluate_ℓ(ℓ, randn(K))
        H = Hamiltonian(GaussianKineticEnergy(Σ), ℓ)
        qs = Array{Float64}(undef, N, K)
        ϵ = 0.5
        algorithm = NUTS()
        for i in 1:N
            Q = first(sample_tree(RNG, algorithm, H, Q, ϵ))
            qs[i, :] = Q.q
        end
        m, C = mean_and_cov(qs, 1)
        @test vec(m) ≈ μ atol = 0.1 rtol = maximum(diag(C))*0.02 norm = x -> norm(x,1)
        @test cov(qs, dims = 1) ≈ L*L' atol = 0.1 rtol = 0.1
    end
end

###
### test proposals
###

# NOTE superseded by separating ω and ζ
# @testset "proposal" begin
#     trajectory = TrajectoryNUTS(nothing, 0, 1, -1000)
#     function test_sample(rng, prop1, prop2, bias, prob_prob2; atol = 0.02, N = 10000)
#         count = 0
#         for _ in 1:N
#             prop = combine_proposals(rng, trajectory, prop1, prop2,
#                                      # direction irrelevant for thie method
#                                      rand(Bool),
#                                      # test with given bias
#                                      bias)
#             if prop.z ≂ prop2.z
#                 count += 1
#             else
#                 @test prop.z ≂ prop1.z
#             end
#             @test prop.ω ≈ logaddexp(prop1.ω, prop2.ω)
#         end
#         @test count / N ≈ prob_prob2 atol = atol
#     end
#     prop1 = Proposal(1, log(1.0))
#     prop2 = Proposal(2, log(3.0))
#     prop3 = Proposal(3, log(1/3))
#     test_sample(RNG, prop1, prop2, true, 1; atol = 0, N = 100)
#     test_sample(RNG, prop1, prop2, false, 0.75)
#     test_sample(RNG, prop1, prop3, true, 1/3)
#     test_sample(RNG, prop1, prop3, false, 0.25)
# end


text:
using SKTB.FermiContour
using Test

# Compare the rate of convergence with the predicted rate.
# Not a proper unit test, but good enough for our purposes.

@testset "FermiContour" begin

E1 = 0
E2 = 1
β = 10.0
n = 2:2:18
x = range(E1,stop=E2,length=100)

# Compute empirical convergence rate
error = zeros(length(n))
for i = 1:length(n)
    w,z = fermicontour(E1,E2,β,n[i])
    fx = sum(real(w*fermidirac.(z,β) ./ (z  .- x)) for (w,z) in zip(w,z))
    error[i] = maximum(abs.(fermidirac.(x,β) .- fx))
end

# Compute theoretical convergence rate
m = E1^2 + π^2/β^2
M = E2^2 + π^2/β^2
k = (sqrt(M/m)-1)/(sqrt(M/m)+1)
K = SKTB.FermiContour.JacobiFunc.K(k^2)
iK = SKTB.FermiContour.JacobiFunc.iK(k^2)
predicted = 2 * exp.(-π*iK/(2K)*n)

println("n   error    predicted    err/pred")
display([n  error predicted error./predicted]); println()
@test maximum(abs.(error./predicted .- 1.0)) < 0.1

# TODO: implement analogous test for 0T contour!

end


text:
Pkg.activate("./docs")

using Documenter, EndowmentsLH, FilesLH

makedocs(
    modules = [EndowmentsLH],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "hendri54",
    sitename = "EndowmentsLH.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

pkgDir = rstrip(normpath(@__DIR__, ".."), '/');
@assert endswith(pkgDir, "EndowmentsLH")
deploy_docs(pkgDir; trialRun = false);

Pkg.activate(".")

# deploydocs(
#     repo = "github.com/hendri54/EndowmentsLH.jl.git",
#     push_preview = true
# )


text:
BigFloat(1.75)              # Accidentally OK, since 1.75 is exactly represented by Float64

BigFloat(2.1)               # First rounded to Float64, low precision

parse(BigFloat, "2.1")      # High precision

setprecision(512) do
    BigFloat(1) / BigFloat(3)
end

display(BigFloat(pi))
display(sin(BigFloat(pi)))

function halley_cuberoot(a, ε)
    x = 1    # Initial guess
    Δx = 1   # Need to initialize to get started
    while Δx ≥ ε
        xnew = x - 2(x^3 - a)*3x^2 / (2(3x^2)^2 - (x^3 - a)*(6x))
        Δx = abs(xnew - x)
        x = xnew
        println("Error = ", abs(x - cbrt(a)))
    end
    x
end

halley_cuberoot(BigFloat(2), 1e-60)


text:
@testset "DATA" begin
	# UCI datasets
    dataset = "abalone"
	data, _, _ = UCI.get_loda_data(dataset)
	cdata = UCI.normalize(hcat(data.normal, data.easy, data.medium))  # this is a control array
	na = size(data.easy,2) + size(data.medium,2)

	(tr_x, tr_y), (val_x, val_y), (tst_x, tst_y) = GenerativeAD.load_data(dataset)
	@test size(tr_x,1) == size(data.normal,1)
	@test size(tr_x,2) == floor(Int,size(data.normal,2)*0.6) == length(tr_y)
	@test sum(tr_y) == 0
	@test size(val_x,2) == floor(Int,size(data.normal,2)*0.2) + floor(Int,na*0.5) == length(val_y) > sum(val_y) > 0
	@test size(tst_x,2) == floor(Int,size(data.normal,2)*0.2) + floor(Int,na*0.5) == length(tst_y) > sum(tst_y) > 0

	# test contamination
	(tr_x, tr_y), (val_x, val_y), (tst_x, tst_y) = GenerativeAD.load_data("waveform-1", contamination=0.1)
	@test abs(sum(tr_y)/length(tr_y) - 0.1) < 0.001
	@test sum(val_y) == sum(tst_y)

	# img datasets
	function test_img_data(dataset::String, dims, contamination, method="leave-one-out")
		(tr_x, tr_y), (val_x, val_y), (tst_x, tst_y) = GenerativeAD.load_data(dataset, contamination=contamination, method=method)
		h,w,c,n = dims
		@test size(tr_x)[1:3] == size(val_x)[1:3] == size(tst_x)[1:3] == (h,w,c)
		@test n-3 <= size(tr_x,4) + size(val_x,4) + size(tst_x,4) <= n # some samples are left out due to rounding
		@test abs(sum(tr_y)/length(tr_y) - contamination) < 0.001
		if method == "leave-one-in"
			@test size(tr_x,4) < size(val_x,4)
		end
	end

	test_img_data("MNIST", (28, 28, 1, 70000), 0.0)
	test_img_data("FashionMNIST", (28, 28, 1, 70000), 0.0)
	test_img_data("MNIST", (28, 28, 1, 70000), 0.05)
	test_img_data("FashionMNIST", (28, 28, 1, 70000), 0.05)
	test_img_data("MNIST", (28, 28, 1, 70000), 0.2, "leave-one-in")
	test_img_data("FashionMNIST", (28, 28, 1, 70000), 0.2, "leave-one-in")

	# this takes a lot of time...
	if !complete
		@info "Skipping tests on CIFAR10 and SVHN2..."
	else
		@info "Testing CIFA10 and SVHN2 datasets"
		test_img_data("CIFAR10", (32, 32, 3, 60000), 0.0)
		test_img_data("SVHN2", (32, 32, 3, 99289), 0.0)
	end
end


text:
using FileIO, BedgraphFiles
using Bedgraph

using IteratorInterfaceExtensions
using TableTraits

using Test
using Logging

# old_logger = global_logger(ConsoleLogger(stdout, Logging.Debug))

module Bag

    using Bedgraph

    const chroms = ["chr19", "chr19", "chr19", "chr19", "chr19", "chr19", "chr19", "chr19", "chr19"]
    const firsts = [49302000, 49302300, 49302600, 49302900, 49303200, 49303500, 49303800, 49304100, 49304400]
    const lasts = [49302300, 49302600, 49302900, 49303200, 49303500, 49303800, 49304100, 49304400, 49304700]
    const values = [-1.0, -0.75, -0.50, -0.25, 0.0, 0.25, 0.50, 0.75, 1.00]

    const record = Bedgraph.Record("chr1", 1, 1, 0)
    const records = Bedgraph.Record.(Bag.chroms, Bag.firsts, Bag.lasts, Bag.values)

    const file = joinpath(@__DIR__, "data.bedgraph")
    const file_headerless = joinpath(@__DIR__, "data-headerless.bedgraph")

    const tmp_output_path = tempname() * ".bedgraph"

end # module Bag

using .Bag

@testset "BedgraphFiles" begin

    @test isfile(Bag.file) == true
    @test isfile(Bag.file_headerless) == true

    # Load tests.
    loader = load(Bag.file)
    @test IteratorInterfaceExtensions.isiterable(loader) == true
    @test TableTraits.isiterabletable(loader) == true

    loaded = Vector{Bedgraph.Record}(loader)
    @test Vector{Bedgraph.Record} == typeof(loaded)

    loader_from_headerless = load(Bag.file_headerless)
    @test isiterable(loader_from_headerless) == true
    @test TableTraits.isiterabletable(loader_from_headerless) == true

    loaded_from_headerless = Vector{Bedgraph.Record}(loader_from_headerless)
    @test Vector{Bedgraph.Record} == typeof(loaded_from_headerless)

    @test IteratorInterfaceExtensions.isiterable(Bag.records) == true
    @test TableTraits.isiterabletable(Bag.records) == true

    @test Bag.records == loaded
    @test Bag.records == loaded_from_headerless

    # Save and load from Vector{Bedgraph.Record}.
    save(Bag.tmp_output_path, Bag.records)

    @debug "direct load into Vector{Bedgraph.Record} - commencing"
    @test Bag.records == Vector{Bedgraph.Record}(load(Bag.tmp_output_path))
    @test Bag.records == collect(Bedgraph.Record, load(Bag.tmp_output_path))
    @debug "direct load into Vector{Bedgraph.Record} - complete"

    @test Bag.records == load(Bag.tmp_output_path) |> Vector{Bedgraph.Record}

    # Save using pipe.
    Bag.records |> save(Bag.tmp_output_path)
    @test Bag.records == Vector{Bedgraph.Record}(load(Bag.tmp_output_path))
    @test Bag.records == load(Bag.tmp_output_path) |> Vector{Bedgraph.Record}

    # Check return of data from save method.
    @test Bag.records == Bag.records |> save(Bag.tmp_output_path)

    # Check convert
    @test Bag.record == convert(Bedgraph.Record, (chrom = "chr1", first=1, last=1, value=0))
    @test Bag.record == convert(Bedgraph.Record, (chrom = "chr1", left=1, right=1, value=0))

    @testset "Integrations" begin
        include("integrations/test-Tables.jl")
        include("integrations/test-QueryOperators.jl")
    end # testset Integrations

    println()
    show(load(Bag.file))
    println()

end


text:
"""
    TwoSectionView{T<:Real} <: LightGraphs.SimpleGraphs.AbstractSimpleGraph{Int64}

Represents a 2-section view of a hypergraph `h`.
Note (1) this is a view - changes to the original hypergraph will be automatically reflected in the view.

Note (2) The view will only work correctly for hypergraphs not having overlapping hyperedges.
To check whether a graph has overlapping edges try `has_overlapping_hedges(h)` - for such graph
you need to fully materialize it rather than use a view.
This can be achieved via the `get_twosection_adjacency_mx(h)` method.

**Constructors**

TwoSectionView(::Hypergraph)

The 2-section view of a hypergraph is suitable for processing with the LightGraphs.jl package.
Several LightGraphs methods are provided for the compability.

"""
struct TwoSectionView{T<:Real} <: LightGraphs.SimpleGraphs.AbstractSimpleGraph{Int}
    h::Hypergraph{T}
    TwoSectionView(h::Hypergraph{T}) where { T <: Real } = begin
        has_overlapping_hedges(h) && error("A two section view can be created only for a graph with non overlapping edges")
        new{T}(h)
    end
end

"""
    has_overlapping_hedges(h::Hypergraph{T}) where { T <: Real }

Checks whether a hypergraph has hyperedges connecting the same pairs of vertices.
Such hypergraph cannot be presented as a two section view
"""

function has_overlapping_hedges(h::Hypergraph{T}) where { T <: Real }
    minimum(size(h)) == 0 && return false
    maximum(get_twosection_adjacency_mx(h;replace_weights=1)) > 1
end

"""
  Return the number of vertices in 2-section view `t` of a hypergraph.
"""
LightGraphs.nv(t::TwoSectionView) = length(t.h.v2he)

LightGraphs.vertices(t::TwoSectionView) = Base.OneTo(nv(t))

"""
  Return the number of edges in 2-section view `t` of a hypergraph.
"""
function LightGraphs.ne(t::TwoSectionView)
    s = 0
    for x in t.h.he2v
        s += length(x) * (length(x) - 1)
    end
    div(s, 2)
end

"""
    LightGraphs.all_neighbors(t::TwoSectionView, v::Integer)

Returns N(v) (the vertex v is not included in N(v))
"""
function LightGraphs.all_neighbors(t::TwoSectionView, v::Integer)
    neighbors = Set{Int}()
    for he in keys(t.h.v2he[v])
        union!(neighbors, keys(t.h.he2v[he]))
    end
    delete!(neighbors, v) #remove v from its neighborhood
    collect(neighbors) #returns the corresponding array
end

function LightGraphs.has_edge(t::TwoSectionView, s, d)
    s == d && return false
    !isempty(intersect(keys(t.h.v2he[s]), keys(t.h.v2he[d])))
end

LightGraphs.has_vertex(t::TwoSectionView, v::Integer) = 1 <= v <= LightGraphs.nv(t)

LightGraphs.outneighbors(t::TwoSectionView, v::Integer) =
    LightGraphs.all_neighbors(t::TwoSectionView, v)

LightGraphs.inneighbors(t::TwoSectionView, v::Integer) =
    LightGraphs.all_neighbors(t::TwoSectionView, v)

"""
    LightGraphs.SimpleGraph(t::TwoSectionView)

Creates a `LightGraphs.SimpleGraph` representation of a `TwoSectionView` t.

This creates a copy of the date. Note that the weights information is not stored
in the created `SimpleGraph`.
"""
function LightGraphs.SimpleGraph(t::TwoSectionView)
    g = SimpleGraph(nv(t))
    for v in LightGraphs.vertices(t)
        neighbors_v = LightGraphs.all_neighbors(t, v)
        for neighbor in neighbors_v
            add_edge!(g, v, neighbor)
        end
    end
    g
end

LightGraphs.is_directed(t::TwoSectionView{T}) where T = false

LightGraphs.is_directed(::Type{TwoSectionView{T}}) where T = false

Base.eltype(::TwoSectionView{T}) where T = Int

"""
    shortest_path(t::TwoSectionView,source::Int, target::Int)

Finds a single shortest path in a graph `b` between vertices
`source` and `target`.
Note that if several paths of the same length exist, only one
will be returned.

"""
function shortest_path(t::TwoSectionView, source::Int, target::Int)
    checkbounds(t.h.v2he, source)
    checkbounds(t.h.v2he, target)
    dj = dijkstra_shortest_paths(t, source)
    enumerate_paths(dj)[target]
end

"""
    LightGraphs.SimpleGraphs.fadj(t::TwoSectionView)

Generates an adjency list for this view of a hypergraph.
"""
function LightGraphs.SimpleGraphs.fadj(t::TwoSectionView)
    res = [Vector{Int}() for _ in 1:LightGraphs.nv(t)]
    for he in t.h.he2v
        vs = collect(keys(he))
        if length(vs) > 1
            for i in 1:length(vs)
                append!(res[vs[i]], vs[1:end .!= i])
            end
        end
    end
    sort!.(res)
end
LightGraphs.SimpleGraphs.fadj(t::TwoSectionView, v::Integer) = LightGraphs.all_neighbors(t,v)
LightGraphs.edges(t::TwoSectionView) = LightGraphs.SimpleGraphs.SimpleEdgeIter(t)

LightGraphs.edgetype(t::TwoSectionView{T}) where T = LightGraphs.SimpleGraphs.SimpleEdge{Int}

LightGraphs.zero(t::TwoSectionView{T}) where T = TwoSectionView(Hypergraph{T}(0,0))
LightGraphs.zero(::Type{TwoSectionView{T}}) where T = TwoSectionView(Hypergraph{T}(0,0))

"""
    get_twosection_adjacency_mx(h::Hypergraph{T,V,E}; count_self_loops::Bool=false,
                                replace_weights::Union{Nothing,Real}=nothing ) where {T<:Real, V, E}

Returns an adjacency matrix for a two section view of a hypergraph `h`.
"""
function get_twosection_adjacency_mx(h::Hypergraph{T,V,E}; count_self_loops::Bool=false,
                                     replace_weights::Union{Nothing,Real}=nothing ) where {T<:Real, V, E}
    mx = zeros(replace_weights==nothing ? T : typeof(replace_weights), nhv(h), nhv(h))
    for he in 1:nhe(h)
        for v1 in keys(h.he2v[he])
            for v2 in keys(h.he2v[he])
                v1 == v2 && !count_self_loops && continue
                mx[v1,v2] += replace_weights==nothing ? h.he2v[he][v1] : replace_weights
            end
        end
    end
    mx
end


text:
using OrdinaryDiffEq

function SinCosDiffEqToSolve!(ddu,du,u,p,t)
    ddu[1] = -u[1]
end

prob = SecondOrderODEProblem(SinCosDiffEqToSolve!,[2.0],[3.0],(0.0,1.0))
sol = solve(prob, Tsit5())
sol = solve(prob, Vern9())
sol = solve(prob, ROCK4())
sol = solve(prob, Rodas4())
sol = solve(prob, AutoVern9(Rodas4()))
sol = solve(prob, TRBDF2())


text:
## Use case for the paper submitted to JOSS
# We create two plots using Mangal.jl and EcologicalNetworks.jl
# First plot: association between the number of species and the number of links in different types of networks
# Second plot: association between richness, connectance, nestedness, and modularity in food webs

## upload packages
using EcologicalNetworks
using Mangal
using DataFrames
using Plots

## read data from mangal.io

# read the id and other metadata of all ecological networks archived on mangal.io
# id: network ID
# S: number of species in the network
# L: number of interactions in the network
# pred: number of interactions of predation in the network
# herb: number of interactions of herbivory in the network
# mutu: number of interactions of mutualism in the network
# para: number of interactions of parasitism in the network

number_of_networks = count(MangalNetwork)
count_per_page = 100
number_of_pages = convert(Int, ceil(number_of_networks/count_per_page))

mangal_networks = DataFrame(fill(Int64, 7),
                 [:id, :S, :L, :pred, :herb, :mutu, :para],
                 number_of_networks)

global cursor = 1
@progress "Paging networks" for page in 1:number_of_pages
    global cursor
    networks_in_page = Mangal.networks("count" => count_per_page, "page" => page-1)
    @progress "Counting items" for current_network in networks_in_page
        S = count(MangalNode, current_network)
        L = count(MangalInteraction, current_network)
        pred = count(MangalInteraction, current_network, "type" => "predation")
        herb = count(MangalInteraction, current_network, "type" => "herbivory")
        mutu = count(MangalInteraction, current_network, "type" => "mutualism")
        para = count(MangalInteraction, current_network, "type" => "parasitism")
        mangal_networks[cursor,:] .= (current_network.id, S, L, pred, herb, mutu, para)
        cursor = cursor + 1
    end
end

## classify all networks according to their most frequent type of interactions

# number of interactions in food webs (sum of predation and herbivory interactions)
mangal_networks.foodweb = mangal_networks.pred .+ mangal_networks.herb

# number of other types of interactions
mangal_networks.other = mangal_networks.L .- mangal_networks.foodweb .- mangal_networks.mutu .- mangal_networks.para

# remove networks with less than 5 links
mangal_networks2 = mangal_networks[mangal_networks.L .> 4, :]

# classify network according to their most frequent type of interactions
number_interactions_max_type = maximum.(eachrow(mangal_networks2[:,[:mutu, :para,:foodweb,:other]]))

foodwebs = mangal_networks2[mangal_networks2[:foodweb] .== number_interactions_max_type ,:]
parasitism_webs = mangal_networks2[mangal_networks2[:para] .== number_interactions_max_type ,:]
mutualism_webs = mangal_networks2[mangal_networks2[:mutu] .== number_interactions_max_type ,:]
other_webs = mangal_networks2[mangal_networks2[:other] .== number_interactions_max_type ,:]

## read food webs from mangal.io and convert them to unipartie networks using EcologicalNetworks.jl

# read all food webs archived on mangal.io
mangal_foodwebs = network.(foodwebs.id)

# convert MangalNetworks to UnipartiteNetworks when possible
unipartite_foodwebs=[]

for i in eachindex(mangal_foodwebs)
    try
        unipartite_foodweb = convert(UnipartiteNetwork, mangal_foodwebs[i])
        push!(unipartite_foodwebs, unipartite_foodweb)
    catch
        println("Cannot convert mangal food web $(i) to a unipartite network")
    end
end

## compute food-web measures (i.e. richness, connectance, nestedness and modularity)

foodweb_measures = DataFrame(fill(Float64, 4),
                 [:rich, :connect, :nested, :modul],
                 length(unipartite_foodwebs))

# compute richness and connectance
foodweb_measures.rich = richness.(unipartite_foodwebs)
foodweb_measures.connect = connectance.(unipartite_foodwebs)

# compute nestedness (spectral radius of the adjacency matrix)
foodweb_measures.nested = ρ.(unipartite_foodwebs)

# compute modularity (BRIM algorithm)
number_of_modules = repeat(3:15, outer=100)
modules = Array{Dict}(undef, length(number_of_modules))

for i in eachindex(unipartite_foodwebs)
    current_network = unipartite_foodwebs[i]
    for j in eachindex(number_of_modules)
        _, modules[j] = n_random_modules(number_of_modules[j])(current_network) |> x -> brim(x...)
    end
    partition_modularity = map(x -> Q(current_network,x), modules);
    foodweb_measures.modul[i] = maximum(partition_modularity)
end

## plot 1: association between the number of species and the number of links in different types of networks

# color palette
pal = (
    foodwebs=RGB(204/255,121/255,167/255),
    parasitism=RGB(230/255,159/255,0/255),
    mutualism=RGB(0/255,158/255,115/255),
    other=RGB(86/255,190/255,233/255),
    )

scatter(mutualism_webs.S, mutualism_webs.L,
    alpha=0.2, color=pal.mutualism, markersize=5, markershape=:rect,
    lab="mutualism",
    legend=:topleft, foreground_color_legend=nothing, background_color_legend=:white,
    dpi=1000, framestyle=:box)
scatter!(parasitism_webs.S, parasitism_webs.L,
    alpha=0.2, color=pal.parasitism,  markersize=5, markershape=:diamond,
    lab="parasitism")
scatter!(foodwebs.S, foodwebs.L,
    alpha=0.2, color=pal.foodwebs, markersize=5,
    lab="food webs")
scatter!(other_webs.S, other_webs.L,
    alpha=0.2, color=pal.other,  markersize=5, markershape=:utriangle,
    lab="other types")
xaxis!(:log, "Number of species")
yaxis!(:log, "Number of links")
savefig("paper/fig/LS.png")

## plot 2: association between richness, connectance, nestedness, and modularity in food webs

plotA = scatter(foodweb_measures.connect, foodweb_measures.nested,
    markersize=foodweb_measures.rich ./ 30,
    alpha=0.6, color=pal.foodwebs,
    lab="",
    dpi=1000, framestyle=:box)
xaxis!("Connectance")
yaxis!((0.4,0.9), "Nestedness")

plotB = scatter(foodweb_measures.connect, foodweb_measures.modul,
    markersize=foodweb_measures.rich ./ 30,
    alpha=0.6, color=pal.foodwebs,
    lab="",
    dpi=1000, framestyle=:box)
xaxis!("Connectance")
yaxis!("Modularity")

plotC = scatter(foodweb_measures.modul, foodweb_measures.nested,
    markersize=foodweb_measures.rich ./ 30,
    alpha=0.6, color=pal.foodwebs,
    lab="",
    dpi=1000, framestyle=:box)
xaxis!("Modularity")
yaxis!((0.4,0.9), "Nestedness")

plot(plotA, plotB, plotC, layout=(1,3),
    title=["($i)" for j in 1:1, i in 1:3], titleloc=:right, titlefont = font(12),
    size=(700,350), margin=5Plots.mm)
savefig("paper/fig/nestmod.png")


text:
abstract Position

type GridPosition
  x :: Float64
  y :: Float64
end
uGrid = GridPosition(0.0,0.5)
vGrid = GridPosition(0.5,0.0)
mGrid = GridPosition(0.0,0.0)
sGrid = GridPosition(0.5,0.5)

type GridMoving
  xshift :: Float64
  yshift :: Float64
  xband :: Float64
  yband :: Float64
end

function GridMoving( o :: GridPosition, t :: GridPosition, xband :: Float64, yband :: Float64)
  xshift = t.x - o.x
  yshift = t.y - o.y
  return GridMoving(xshift, yshift, xband, yband)
end

function intercept( ss :: Array{Float64,2}, g :: GridMoving; keep_when_notmoving :: Bool = false)
  signal = ss
  xband = g.xband
  yband = g.yband
  xshift = g.xshift
  yshift = g.yshift
  kx = collect(1:(size(signal,1)/xband)) / size(signal,1)
  ky = collect(1:(size(signal,2)/yband)) / size(signal,2)

  fresx = zeros(signal)
  if ((xshift == 0.0) & keep_when_notmoving)
    fresx = signal
  else
    for ix = kx
      psix = ix*pi*collect(1:size(signal,1))
      psix1 = ix*pi*(xshift + collect(1:size(signal,1)))
      fsinx  = mean(sin(psix) .* signal,1).* sin(psix1)/(0.5)
      fcosx  = mean(cos(psix) .* signal,1).* cos(psix1)/(0.5)
      fresx += fsinx + fcosx
    end
  end

  signal = fresx
  fresy = zeros(size(signal))
  if ((yshift == 0.0) & keep_when_notmoving)
    fresy = signal
  else
    for iy = ky
      psiy = iy*pi*collect(1:size(signal,2))
      psiy1 = iy*pi*(yshift + collect(1:size(signal,2)))
      psiy = psiy'
      psiy1 = psiy1'
      fsiny  = mean(sin(psiy) .* signal,2).* sin(psiy1)/(0.5)
      fcosy  = mean(cos(psiy) .* signal,2).* cos(psiy1)/(0.5)
      fresy += fsiny + fcosy
    end
  end

  return fresy
end

function intercept(ss:: Array{Float64,3}, g:: GridMoving; kws...)
  res = zeros(ss)
  for iz = 1:size(ss,3)
    res[:,:,iz] = intercept(ss[:,:,iz],g;kws...)
  end
  return res
end
export intercept
export GridMoving


text:
using Test
using PyCall
using BenchmarkTools
using Fermi.Wavefunction
using Fermi.MollerPlesset
using LinearAlgebra.BLAS
using Fermi
BLAS.set_num_threads(4)

psi4.core.set_num_threads(4)
psi4.core.be_quiet() #turn off output
# > setup
tol = 1E-14
psi4.set_options(Dict("D_CONVERGENCE" => 10,
					  "scf_type" => "pk"))
mol2 = psi4.geometry("""
                     pubchem:ethane
					 symmetry c1
					 """)
e2,wfn2 = psi4.energy("hf/cc-pvtz",mol=mol2,return_wfn=true)
JuWfn2 = Wfn(wfn2,Float64,false,false)
#mol3 = psi4.geometry("""
#					 1 2
#					 O
#					 H 1 1.1
#					 H 1 1.1 2 104.0
#					 symmetry c1
#					 """)
#psi4.set_options(Dict("reference" => "uhf"))
#e3,wfn3 = psi4.energy("hf/cc-pvdz",mol=mol3,return_wfn=true)
println("Starting benchmark!")
#println(@benchmark Wfn(wfn2,Float64,false,false))
println(@benchmark Wfn(wfn2,Float64,false,false))
#println(@benchmark Wfn(wfn3,Float64,false,true))
println(@benchmark do_rmp2(JuWfn2))
#println(@benchmark do_direct_rmp2(JuWfn2))
println(@benchmark do_df_rmp2(JuWfn2))
println(do_rmp2(JuWfn2))
println(do_df_rmp2(JuWfn2))
#println(@benchmark do_ump2(JuWfn3))
#println(@benchmark do_ump2(JuWfn4))


text:
using Gadfly, DataArrays, RDatasets

 plot(dataset("ggplot2", "diamonds"), x="Price", Geom.density)


text:
#########################################
# different types of domains and ranges #
#########################################

abstract type AbstractFiniteSet end
const FinSet = AbstractFiniteSet

abstract type AbstractSegment <: FinSet end
abstract type AbstractProduct <: FinSet end

####### INTERFACES ###########

iterator(finset :: FinSet) = error("not implemented")
factors(aprod :: AbstractProduct) = error("not implemented")

#######################
# iteration utilities #
#######################

first(s :: AbstractSegment) = first(iterator(s))
last(s :: AbstractSegment) = last(iterator(s))
iterate(s :: FinSet) = iterate(iterator(s))
iterate(s :: FinSet, state) = iterate(iterator(s), state)
length(s :: FinSet) = length(iterator(s))
ndims(s :: FinSet) = ndims(iterator(s))
size(s :: FinSet) = size(iterator(s))

proj(aprod :: AbstractProduct, i) = factors(aprod)[i]
getindex(aprod :: AbstractProduct, i) = proj(aprod, i)

#####################
# main files ########
#####################
include("finsets/segment.jl")
include("finsets/direct_prod.jl")

function Base.show(io :: IO, s :: AbstractSegment)
    print(io, "[", first(s), "..", last(s), "]")
end

function Base.show(io :: IO, dp :: DirectProduct)
    res = ""
    for x in factors(dp)
        res *= repr(x)
        res *= "×"
    end
    print(io, res[1 : end-1])
end

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
function sub_quote(W, @nospecialize(T), fast::Bool)
  vtyp = vtype(W, T)
  instrs = "%res = fneg $(fast_flags(fast)) $vtyp %0\nret $vtyp %res"
  quote
    $(Expr(:meta, :inline))
    Vec($LLVMCALL($instrs, _Vec{$W,$T}, Tuple{_Vec{$W,$T}}, data(v)))
  end
end

@generated Base.:(-)(v::Vec{W,T}) where {W, T <: Union{Float32,Float64}} = sub_quote(W, T, false)
@generated Base.FastMath.sub_fast(v::Vec{W,T}) where {W, T <: Union{Float32,Float64}} = sub_quote(W, T, true)

@inline Base.:(-)(v::Vec{<:Any,<:NativeTypes}) = FastMath.sub_fast(zero(v), v)
@inline Base.FastMath.sub_fast(v::Vec{<:Any,<:NativeTypes}) = FastMath.sub_fast(zero(v), v)

@inline Base.inv(v::AbstractSIMD{W,<:FloatingTypes}) where {W} = FastMath.div_fast(one(v), v)
@inline Base.inv(v::AbstractSIMD{W,<:IntegerTypesHW}) where {W} = inv(float(v))
@inline FastMath.inv_fast(v::AbstractSIMD{<:Any,<:IntegerTypesHW}) = FastMath.inv_fast(float(v))
@inline FastMath.inv_fast(v::AbstractSIMD{<:Any,<:FloatingTypes}) = FastMath.div_fast(one(v), v)

@inline Base.abs(v::AbstractSIMD{W,<:Unsigned}) where {W} = v
@inline Base.abs(v::AbstractSIMD{W,<:Signed}) where {W} = ifelse(v > 0, v, -v)

@inline Base.round(v::AbstractSIMD{W,<:Integer}) where {W} = v
@inline Base.round(v::AbstractSIMD{W,<:Integer}, ::RoundingMode) where {W} = v


text:
# interpolate point at (x, y) from 4 nearby pixel values
function interp_bilinear{T <: Union(Real, Integer)}(dat::Array{T, 2},
                                                    x::Float64, y::Float64,
                                                    x1, x2, y1, y2)
    q11 = dat[y1, x1]
    q12 = dat[y2, x1]
    q21 = dat[y1, x2]
    q22 = dat[y2, x2]
    if x1 != x2
        r1 = (x2 - x) / (x2 - x1) * q11 + (x - x1) / (x2 - x1) * q21
        r2 = (x2 - x) / (x2 - x1) * q12 + (x - x1) / (x2 - x1) * q22
    else
        # special case of x1 == x2, no interpolation needed
        r1 = q11
        r2 = q12
    end
    if y1 != y2
        r = (y - y1) / (y2 - y1) * r1 + (y2 - y) / (y2 - y1) * r2
    else
        # special case of y1 == y2, no interpolation needed
        r = r1
    end
    if typeof(r) <: Real
        r = round(r)
    end
    r = convert(T, r)
    return r
end

function imresize{T <: Union(Real, Integer)}(dat::Array{T, 2},
                                             new_size::(Int, Int))
    new_dat = similar(dat, new_size)
    h, w = size(dat)
    new_h, new_w = new_size
    for new_j=1:new_w, new_i=1:new_h
        # coordinates in original image
        x = new_j * w / new_w
        y = new_i * h / new_h
        # coordinates of 4 points to interpolate from
        x1, x2 = max(1, floor(x)), min(w, ceil(x))
        y1, y2 = max(1, floor(y)), min(h, ceil(y))
        new_dat[new_i, new_j] = interp_bilinear(dat, x, y, x1, x2, y1, y2)
    end
    return new_dat
end

function imresize{T <: Union(Real, Integer)}(dat::Array{T, 2},
                                             new_size...)
    return imresize(dat, new_size)
end

## if isinteractive()
##     using Images
##     using ImageView
##     im = imread("../../data/CK/faces_aligned/S999_003_00000054.png")
##     dat = convert(Array, im)
## end

using Images, Cairo

function imresize_julia!(resized, original)
    scale1 = (size(original,1)-1)/(size(resized,1)-0.999f0)
    scale2 = (size(original,2)-1)/(size(resized,2)-0.999f0)
    for jr = 0:size(resized,2)-1
        jo = scale2*jr
        ijo = itrunc(jo)
        fjo = jo - oftype(jo, ijo)
        @inbounds for ir = 0:size(resized,1)-1
            io = scale1*ir
            iio = itrunc(io)
            fio = io - oftype(io, iio)
            tmp = (1-fio)*((1-fjo)*original[iio+1,ijo+1] +
                           fjo*original[iio+1,ijo+2])
            + fio*((1-fjo)*original[iio+2,ijo+1] +
                   fjo*original[iio+2,ijo+2])
            resized[ir+1,jr+1] = convertsafely(eltype(resized), tmp)
        end
    end
    resized
end
imresize_julia(original, new_size) =
    imresize_julia!(similar(original,
                            new_size), original)

convertsafely{T<:FloatingPoint}(::Type{T}, val) = convert(T, val)
convertsafely{T<:Integer}(::Type{T}, val::Integer) = convert(T, val)
convertsafely{T<:Integer}(::Type{T}, val::FloatingPoint) =
    itrunc(T,
           val+oftype(val, 0.5))

function imresize_cairo(dat::Array{Uint32, 2}, new_size::(Int, Int))
    cs = CairoImageSurface(dat, 0)
    new_dat = zeros(Uint32, new_size)
    new_cs = CairoImageSurface(new_dat, 0)
    pat = CairoPattern(cs)
    pattern_set_filter(pat, Cairo.FILTER_BILINEAR)
    c = CairoContext(new_cs)
    h, w = size(dat)
    new_h, new_w = new_size
    scale(c, new_h / h, new_w / w)
    set_source(c, pat)
    paint(c)
    return new_cs.data
end

## img = rand(0x00:0xff, 774, 512)
## new_size = (3096, 2048)
## imresize_cairo(convert(Array{Uint32}, img), new_size)
## println("Cairo:")
## @time imresize_cairo(convert(Array{Uint32}, img), new_size)
## imresize_julia(img, new_size)
## println("Julia:")
## @time imresize_julia(img, new_size)

max_stars_repo_name: dfdx/masque

max_issues_repo_name: dfdx/masque

max_forks_repo_name: dfdx/masque


text:
using Gadfly
using Compose
using Distributions
import Cairo, Fontconfig
using DataFrames

# %%
g = Geometric(1/2)

# %%
x_min = [0 1/2 3/4 7/8 15/16]
x_max = [1/2 3/4 7/8 15/16 31/32]
y = [2 4 8 16 32]

label = ["H" "TH" "TTH" "..."]
text_x = [1/4 5/8 13/16 29/32]

tick_labels=Dict(zip(x_min, ["0" "1/2" "1/4" "1/8" "1/16"]))

theme = Theme(major_label_font="PragmataPro Mono",
              minor_label_font="PragmataPro Mono",
              key_title_font="PragmataPro Mono Bold",
              highlight_width=0.005w,
              bar_highlight=color("black"))
Gadfly.push_theme(theme)

p = plot(
    x_min=x_min,
    x_max=x_max,
    y=y,
    Geom.bar,
    Coord.cartesian(ymax=20),
    Guide.xticks(ticks=x_min),
    Guide.yticks(ticks=y),
    Scale.x_continuous(labels = x -> tick_labels[x]),
    Guide.xlabel("PROBABILITY"),
    Guide.ylabel("PAYOUT"),
    Guide.title("<b>Distribution of Payouts</b>"),
    Guide.annotation(compose(
        context(),
        text(text_x, ones(5), label, repeat([hcenter], 4), repeat([vcenter], 4)),
        # font("PragmataPro Mono"),
        fontsize(8pt),
        stroke("black"))))
# %%
draw(PDF("../../figures/distribution.pdf"), p)

max_stars_repo_name: akoen/petersburg

max_issues_repo_name: akoen/petersburg

max_forks_repo_name: akoen/petersburg


text:
using MD5
using Test

@testset "String Tests" begin
    # From the reference https://tools.ietf.org/html/rfc1321 (final page)
    @test bytes2hex(md5("")) == "d41d8cd98f00b204e9800998ecf8427e"
    @test bytes2hex(md5("a")) == "0cc175b9c0f1b6a831c399e269772661"
    @test bytes2hex(md5("abc")) == "900150983cd24fb0d6963f7d28e17f72"
    @test bytes2hex(md5("message digest")) == "f96b697d7cb7938d525a2f31aaf161d0"
    @test bytes2hex(md5("abcdefghijklmnopqrstuvwxyz")) == "c3fcd3d76192e4007dfb496cca67e13b"
    @test bytes2hex(md5("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789")) ==  "d174ab98d277d9f5a5611c2c9f419d9f"
    @test bytes2hex(md5("12345678901234567890123456789012345678901234567890123456789012345678901234567890")) == "57edf4a22be3c955ac49da2e2107b67a"

end

@testset "Different forms of input data consistent" begin
    for len in [0,1,10, 100,1000]
        bytes = rand(UInt8, len)
        str = String(copy(bytes)) # String will take ownership and empty it's array inputs
        stream = IOBuffer(bytes)
        @test md5(bytes) == md5(str)
        @test md5(stream) == md5(str)
    end
end

@testset "misc" begin
    @test MD5.digestlen(MD5.MD5_CTX) == length(md5(""))
    @test occursin("MD5", sprint(show, MD5.MD5_CTX()))
end

include("nettle.jl")
include("perf.jl")


text:
# This file was generated, do not modify it. # hide
@show macrorecall(goldI, I1)
@show macrorecall(goldI, I2)
@show macrorecall(goldI, I3)

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
using Printf, PyCall, PyPlot
using AbaqusReader, JLD, MAT

using AD4SM

mean(x) = sum(x)/length(x)
norm(x) = sqrt(sum(x.^2))
#
#
# sort_bnd sorts the nodes of the boundary to
# form a path without self-intersections
function sort_bnd(nodes, nid_bnd)
  closest(node0, nodes) = findmin([norm(node0-node) for node in nodes])[2]

  id_srtd   = [nid_bnd[1],]
  lastnode  = nid_bnd[1]
  bnd_nodes = nid_bnd[2:end]

  while !isempty(bnd_nodes)
    idd = closest(nodes[lastnode], nodes[bnd_nodes])
    lastnode = bnd_nodes[idd]
    push!(id_srtd, bnd_nodes[idd])
    deleteat!(bnd_nodes, idd)
  end
  return id_srtd
end

# get_vol(p0) = π/2*sum((p0[2,1:end-1]+p0[2,2:end]) .*
#                      (p0[1,2:end]-p0[1,1:end-1]) .*
#                      (p0[1,2:end]+p0[1,1:end-1]) )
#
# get_vol calculates the volume of the cavity
get_vol(p0) = π/2*sum(
  (p0[2:2:end-2]+p0[4:2:end]) .*
  (p0[3:2:end]-p0[1:2:end-2]) .*
  (p0[3:2:end]+p0[1:2:end-2]) )
#
# functions for plotting
#
function get_I1(elems, u0)
    F   = Elements.getinfo(elems,u0,info=:F)
    [ Materials.getInvariants(transpose(F)*F)[1] for F in F]
end
function get_I3(elems, u0)
    F   = Elements.getinfo(elems,u0,info=:F)
    [ Materials.getInvariants(transpose(F)*F)[3] for F in F]
end
patch = pyimport("matplotlib.patches")
coll  = pyimport("matplotlib.collections")
function plot_model(elems, nodes;
                    u = zeros(length(nodes[1]), length(nodes)),
                    Φ = [],
                    linewidth = 0.25,
                    facecolor = :c,
                    edgecolor = :b,
                    alpha     = 1,
                    cmap      = :hsv,
                    clim      = [],
                    dTol      = 1e-6,
                    cfig      = figure(),
                    ax        = cfig.add_subplot(1,1,1))

  nodes     = [node + u[:,ii] for (ii,node) in enumerate(nodes)]
  patchcoll = coll.PatchCollection([patch.Polygon(nodes[elem.nodes])
                                    for elem ∈ elems], cmap=cmap)
    if !isempty(Φ)
    patchcoll.set_array(Φ)
    cfig.colorbar.(patchcoll, ax=ax)

    if isempty(clim)
      clim = patchcoll.get_clim()
      if abs(clim[2]-clim[1]) < dTol
        clim  = sum(clim)/2*[0.9, 1.1]
      end
    end
    patchcoll.set_clim(clim)
  else
    patchcoll.set_color(facecolor)
  end

  patchcoll.set_edgecolor(edgecolor)
  patchcoll.set_alpha(alpha)
  patchcoll.set_linewidth(linewidth)

  ax.set_aspect("equal")
  ax.add_collection(patchcoll)
  ax.autoscale()
  (cfig, ax, patchcoll)
end
;

sMeshFile = "AxSymDomainj.inp"
mat       = Materials.MooneyRivlin(1e1, 1e1, 1e3)
Δu        = 18
sPosFix   = "NHd"
N         = 200
LF        = range(1/10N, 1, length=N)
dTol      = 1e-6
ballus    = true
;

sFileName = splitext(basename(sMeshFile))[1];
mymodel   = AbaqusReader.abaqus_read_mesh(sMeshFile)
nodes     = [mymodel["nodes"][ii] for ii in 1:mymodel["nodes"].count]
el_nodes  = [item[2]              for item in mymodel["elements"]]

# for node in nodes
#   node[1]-=10
# end

id_btm   = mymodel["node_sets"]["BTM"]
id_top   = mymodel["node_sets"]["TOP"]
id_bnd   = mymodel["node_sets"]["BND"]

elems     = []
for item ∈ el_nodes
  if length(item)==3
    push!(elems, Elements.ASTria(item, nodes[item], mat=mat))
  elseif length(item)==4
    push!(elems, Elements.ASQuad(item, nodes[item], mat=mat))
  end
end

@show (nNodes, nElems) = (length(nodes), length(elems)); flush(stdout)

# sort the internal boundary nodes
#
id_srtd   = sort_bnd(nodes, id_bnd);
push!(id_srtd, id_srtd[1])
pos0      = hcat(nodes[id_srtd]...)
idxu      = LinearIndices((2,nNodes))
@show V0  = get_vol(pos0)

pos_btm   = hcat(nodes[sort_bnd(nodes, id_btm)]...)
pos_top   = hcat(nodes[sort_bnd(nodes, id_top)]...)

plot_model(elems, nodes)
PyPlot.plot(pos0[1,:], pos0[2,:], color=:b)
PyPlot.plot(pos_btm[1,:], pos_btm[2,:], color=:r)
PyPlot.plot(pos_top[1,:], pos_top[2,:], color=:r)
PyPlot.title("Undeformed domain")
;

# displacement boundary conditions
u0     = 1e-5randn(2, nNodes)
idxu   = LinearIndices(u0)
bifree = trues(size(u0))

bifree[:,id_btm] .= false
bifree[2,id_top] .= false

u0[:,id_btm] .= 0
u0[2,id_top] .= 0

;

@time Elements.solvestep!(elems, u0, bifree,
  bprogress=true, becho=false, dTol=dTol)
;

u            = copy(u0)
u[:,id_btm] .= 0
u[2,id_top] .= -Δu

allus_e = Elements.solve(elems, u, ifree=bifree, LF=range(1/10N, 1, length=N),
                     becho=true, bechoi=false, ballus=ballus, bprogress=false,
                     dTolu=1e-6, dTole=1e-2, maxiter=15)
;

cfig = PyPlot.figure()

ax1  = cfig.add_subplot(1,1,1)
plot_model(elems, nodes, cfig=cfig, ax=ax1,
  u=allus_e[end][1], Φ=get_I3(elems, allus_e[end][1]).-1)
;

nid_bnd = length(id_srtd)
nDoFsu  = 2*nNodes
q       = 10
nGropus = ceil(Int, nid_bnd/q)

id_eqs_DoFs = [(ii-1)*q+1:ii*q+1 for ii in 1:nGropus-1]
append!(id_eqs_DoFs, [id_eqs_DoFs[end][end]:nid_bnd])

eqns = vcat(
  [begin
      id_bnd_ii = id_srtd[idx]
      idx_ii = idxu[:, id_bnd_ii][:]
      Elements.ConstEq(x->get_vol(pos0[:,idx][:]+x[1:end-1])-x[end],
                        vcat(idx_ii, nDoFsu+ii), adiff.D2)
      end for (ii,idx) in enumerate(id_eqs_DoFs)] ...,
  Elements.ConstEq(x->sum(x)-V0, nDoFsu+1:nDoFsu+nGropus, adiff.D2))
@show nEqns = length(eqns)
;

# displacement boundary conditions
idxu   = LinearIndices((2,nNodes))
u0     = 1e-4randn(2, nNodes+ceil(Int, nGropus/2))
bifree = trues(size(u0))
λ      = zeros(nEqns)

# make sure that the volume constraint is initially satisfied
for (ii, Eq) in enumerate(eqns[1:end-1])
  u0[nDoFsu+ii] = Eq.func(zeros(size(Eq.iDoFs)))
end

bifree[:,id_btm] .= false
bifree[2,id_top] .= false

u0[:,id_btm]     .= 0
u0[2,id_top]     .= 0
;

@time Elements.solvestep!(elems,
        u0, bifree, λ=λ, bprogress=true, becho=false,
        dTolu=1e-5, dTole=1e-3, eqns=eqns)
;

u           = copy(u0)
u[:,id_btm] .= 0
u[2,id_top] .= -Δu

allus_d = Elements.solve(elems, u, ifree=bifree, λ=λ, eqns=eqns, LF=range(0, 1, length=200),
                     becho=true, bechoi=false, ballus=ballus, bprogress=false,
                     dTolu=1e-5, dTole=5e-2, maxiter=11)
;

cfig = PyPlot.figure()

ax1  = cfig.add_subplot(1,1,1)
plot_model(elems, nodes, cfig=cfig, ax=ax1,
  u=allus_d[end][1], Φ=get_I1(elems, allus_d[end][1]).-3)
;

cfig = PyPlot.figure()

ax1  = cfig.add_subplot(2,1,1)
plot_model(elems, nodes, cfig=cfig, ax=ax1,
  u=allus_d[end][1], Φ=get_I1(elems, allus_d[end][1]).-3)
PyPlot.plot(pos0[1,:]+allus_d[end][1][1,id_srtd],
  pos0[2,:]+allus_d[end][1][2,id_srtd], color=:b)

ax2  = cfig.add_subplot(2,1,2)
plot_model(elems, nodes, cfig=cfig, ax=ax2,
  u=allus_e[end][1], Φ=get_I1(elems, allus_e[end][1]).-3)
PyPlot.plot(pos0[1,:]+allus_e[end][1][1,id_srtd],
  pos0[2,:]+allus_e[end][1][2,id_srtd], color=:b)

ax1.set_xlim([15, 85]); ax1.set_ylim([-26, 10])
ax1.set_title("with internal volume constraint, \$I_1-3\$")
ax1.set_xticks([]); ax1.set_yticks([])
ax2.set_xlim([15, 85]); ax2.set_ylim([-26, 10])
ax2.set_title("without internal volume constraint, \$I_1-3\$")
ax2.set_xticks([]); ax2.set_yticks([])

rf_tot_d = [sum(item[2][2, id_top])  for item in allus_d]
Δu_tot_d = [mean(item[1][2,id_top])  for item in allus_d]

rf_tot_e = [sum(item[2][2, id_top])  for item in allus_e]
Δu_tot_e = [mean(item[1][2,id_top])  for item in allus_e]

PyPlot.figure()
PyPlot.plot(Δu_tot_d/50, rf_tot_d/1e5, color=:r)
PyPlot.plot(Δu_tot_e/50, rf_tot_e/1e5, color=:b)
PyPlot.xlabel("\$\\Delta u_2/L_2\$")
PyPlot.ylabel("\$F_2\\times10^{-5}\$")
PyPlot.title("Total reaction force")
;

JLD.save(sFileName*".jld", Dict(
            "nodes"=>nodes, "elems"=>elems, "el_nodes"=>el_nodes,
            "allus_d"=>allus_d, "allus_e"=>allus_e,
            "LF"=>LF, "mat"=>mat,
            "id_btm"=>id_btm, "id_top"=>id_top, "id_bnd"=>id_bnd,
            "sMeshFile"=>sMeshFile))
@printf("results written to %s.jld\n", sFileName); flush(stdout)
;

MAT.matwrite(sFileName*".mat", Dict(
        "nodes"=>nodes, "elems"=>el_nodes, "nNodes"=>nNodes,
        "u0_d"=>allus_d[end][1], "u0_e"=>allus_e[end][1],
        "I1_d"=>get_I1(elems, allus_d[end][1]),
        "I1_e"=>get_I1(elems, allus_e[end][1]),
        "id_btm"=>id_btm, "id_top"=>id_top, "id_bnd"=>id_bnd,
        "u_tot_d"=>Δu_tot_d, "u_tot_e"=>Δu_tot_e,
        "rf_tot_d"=>rf_tot_d, "rf_tot_e"=>rf_tot_e))
@printf("results written to %s.mat\n", sFileName); flush(stdout)
;


text:
struct RankAdaptiveUnconventionalAlgorithm_Params{sType, lType, kType}
    S_rhs # rhs of S step (core projected rhs)
    L_rhs # rhs of L step (range projected rhs)
    K_rhs # rhs of K step (corange projected rhs)
    S_kwargs
    L_kwargs
    K_kwargs
    S_alg::sType
    L_alg::lType
    K_alg::kType
    tol
    r_max::Int
end

struct RankAdaptiveUnconventionalAlgorithm{sType, lType, kType} <: AbstractDLRAlgorithm
    alg_params::RankAdaptiveUnconventionalAlgorithm_Params{sType, lType, kType}
end
function RankAdaptiveUnconventionalAlgorithm(tol; S_rhs = nothing, L_rhs = nothing, K_rhs = nothing,
                                             S_kwargs = Dict(), L_kwargs = Dict(), K_kwargs = Dict(),
                                             S_alg = Tsit5(), L_alg = Tsit5(), K_alg = Tsit5(), rmax = 2^62)
    params = RankAdaptiveUnconventionalAlgorithm_Params(S_rhs, L_rhs, K_rhs, S_kwargs, L_kwargs, K_kwargs, S_alg, L_alg, K_alg, tol, rmax)
    return RankAdaptiveUnconventionalAlgorithm(params)
end

struct RankAdaptiveUnconventionalAlgorithm_Cache{uType,SIntegratorType,LIntegratorType,KIntegratorType,yType} <: AbstractDLRAlgorithm_Cache
    US::Matrix{uType}
    Uhat::Matrix{uType}
    VS::Matrix{uType}
    Vhat::Matrix{uType}
    M::Matrix{uType}
    N::Matrix{uType}
    SIntegrator::SIntegratorType
    SProblem
    LIntegrator::LIntegratorType
    LProblem
    KIntegrator::KIntegratorType
    KProblem
    r::Int # current rank
    tol
    r_max
    y
    ycurr::yType
    yprev::yType
    Δy::yType
end

function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorithm, u, dt; t0 = prob.tspan[1])
    # allocate memory for frequently accessed arrays
    tspan = (t0, t0+dt)
    r = rank(u)
    n, m = size(u)
    US = u.U*u.S
    VS = u.V*u.S'
    Uhat = zeros(n,2*r)
    Vhat = zeros(m,2*r)
    M = zeros(2*r,r)
    N = zeros(2*r,r)

    if isnothing(alg.alg_params.K_rhs)
        K_rhs = function (US, V, t)
                    return Matrix(prob.f(TwoFactorRepresentation(US,V),t)*V)
                end
    else
        K_rhs = alg.alg_params.K_rhs
    end
    KProblem = ODEProblem(K_rhs, US, tspan, u.V)
    KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...)
    if isnothing(alg.alg_params.L_rhs)
        L_rhs = function (VS, U, t)
                    return Matrix(prob.f(TwoFactorRepresentation(U,VS),t)'*U)
                end
    else
        L_rhs = alg.alg_params.L_rhs
    end

    LProblem = ODEProblem(L_rhs, VS, tspan, u.U)
    LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...)
    if isnothing(alg.alg_params.S_rhs)
        S_rhs = function (S, (U,V), t)
                    return Matrix(U'*prob.f(SVDLikeRepresentation(U,S,V),t)*V)
                end
    else
        S_rhs = alg.alg_params.S_rhs
    end
    SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, [Uhat, Vhat])
    SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...)

    return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N,
                                                     SIntegrator, SProblem, LIntegrator, LProblem,
                                                     KIntegrator, KProblem, r, alg.alg_params.tol, alg.alg_params.r_max,
                                                     nothing, nothing, nothing, nothing)
end

function init(prob::AbstractDLRProblem, alg::RankAdaptiveUnconventionalAlgorithm, dt)
    t0, tf = prob.tspan
    @assert tf > t0 "Integration in reverse time direction is not supported"
    u = deepcopy(prob.u0)
    # initialize solution
    sol = init_sol(dt, t0, tf, prob.u0)
    # initialize cache
    cache = alg_cache(prob, alg, u, dt; t0 = t0)
    sol.Y[1] = deepcopy(prob.u0) # add initial point to solution object
    return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0)
end

function alg_cache(prob::MatrixDataProblem, alg::RankAdaptiveUnconventionalAlgorithm, u, dt; t0 = prob.tspan[1])
    # creates caches for frequently used arrays by performing the first time step
    @unpack y = prob
    t0 = prob.tspan[1]
    n, r = size(u.U)
    m = size(u.V, 1)

    yprev = y isa AbstractArray ? deepcopy(y[1]) : y(t0)
    ycurr = deepcopy(yprev)
    Δy = similar(yprev)

    US = zeros(n,r)
    VS = zeros(m,r)
    Uhat = zeros(n,2*r)
    Vhat = zeros(m,2*r)
    M = zeros(2*r,r)
    N = zeros(2*r,r)
    KIntegrator = MatrixDataIntegrator(Δy, US, I, u.V, 1)
    LIntegrator = MatrixDataIntegrator(Δy', VS, I, u.U, 1)
    SIntegrator = MatrixDataIntegrator(Δy, zeros(2r,2r), Uhat, Vhat, 1)

    return  RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N,
                                                      SIntegrator, nothing, LIntegrator, nothing,
                                                      KIntegrator, nothing, r, alg.alg_params.tol, alg.alg_params.r_max,
                                                      y, ycurr, yprev, Δy)
end

function alg_recache(cache::RankAdaptiveUnconventionalAlgorithm_Cache, alg::RankAdaptiveUnconventionalAlgorithm, u, t)
    @unpack SProblem, KProblem, LProblem, tol, r_max, y, yprev, ycurr, Δy = cache
    r = rank(u)
    n, m = size(u)

    US = zeros(n,r)
    Uhat = zeros(n,2r)
    M = zeros(2r,r)

    VS = zeros(m,r)
    Vhat = zeros(m,2r)
    N = zeros(2r,r)

    # the following sequence can probably somehow be handled via dispatch in a simpler way
    if isnothing(KProblem)
        KIntegrator = MatrixDataIntegrator(Δy, US, I, u.V, 1)
    else
        KProblem = remake(KProblem, u0 = US, p = u.V, tspan = (t, t+1.0))
        KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...)
    end
    if isnothing(LProblem)
        LIntegrator = MatrixDataIntegrator(Δy', VS, I, u.U, 1)
    else
        LProblem = remake(LProblem, u0 = VS, p = u.U, tspan = (t, t+1.0))
        LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...)
    end

    if isnothing(SProblem)
        SIntegrator = MatrixDataIntegrator(Δy, zeros(2r,2r), Uhat, Vhat, 1)
    else
        SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = [Uhat, Vhat], tspan = (t, t+1.0))
        SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...)
    end

    return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, SIntegrator, SProblem, LIntegrator, LProblem,
                                                     KIntegrator, KProblem, r, tol, r_max, y, ycurr, yprev, Δy)
end

function step!(integrator::DLRIntegrator, alg::RankAdaptiveUnconventionalAlgorithm, dt)
    @unpack u, t, iter, cache, probType = integrator
    u_new, rank_adjusted = rankadaptive_unconventional_step!(u, cache, t, dt, probType)
    if rank_adjusted
        integrator.u = u_new
        integrator.cache = alg_recache(cache, alg, u_new, t+dt)
    end
    integrator.t += dt
    integrator.iter += 1
end

function rankadaptive_unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDataProblem})
    @unpack y, ycurr, yprev, Δy = cache
    update_data!(ycurr, y, t, dt)
    Δy .= ycurr - yprev
    yprev .= ycurr
    rankadaptive_unconventional_step!(u, cache, t, dt)
end

function rankadaptive_unconventional_step!(u, cache, t, dt, ::Type{<:MatrixDEProblem})
    rankadaptive_unconventional_step!(u, cache, t, dt)
end

function rankadaptive_unconventional_step!(u, cache, t, dt)
    @unpack r, r_max, tol, US, Uhat, VS, Vhat, M, N, KIntegrator, LIntegrator, SIntegrator = cache

    # K step
    mul!(US,u.U,u.S)
    set_u!(KIntegrator, US)
    step!(KIntegrator, dt, true)
    US .= KIntegrator.u
    Uhat[:, 1:r] .= US
    Uhat[:, r+1:end] .= u.U
    QRK = qr(Uhat)
    Uhat .= Matrix(QRK.Q)
    mul!(M,Uhat',u.U)

    # L step
    mul!(VS,u.V,u.S')
    set_u!(LIntegrator, VS)
    step!(LIntegrator, dt, true)
    VS .= LIntegrator.u
    Vhat[:, 1:r] .= VS
    Vhat[:, r+1:end] .= u.V
    QRL = qr(Vhat)
    Vhat .= Matrix(QRL.Q)
    mul!(N,Vhat',u.V)

    set_u!(SIntegrator, M*u.S*N')
    step!(SIntegrator, dt, true)

    U, S, V = svd(SIntegrator.u)
    r_new = min(r_max, LowRankArithmetic.truncate_to_tolerance(S, tol))
    if r_new == r
        mul!(u.U,Uhat,U[:,1:r_new])
        u.S .= Matrix(Diagonal(S[1:r_new]))
        mul!(u.V,Vhat,V[:,1:r_new])
        return nothing, false
    else
        println("rank adjusted: new rank = $r_new")
        return SVDLikeRepresentation(Uhat*U[:,1:r_new], Matrix(Diagonal(S[1:r_new])), Vhat*V[:,1:r_new]), true
    end
end


text:
function lb_volume(model::SAFTModel, z = SA[1.0])
    seg = model.params.segment.values
    σᵢᵢ = model.params.sigma.diagvalues
    val = π/6*N_A*sum(z[i]*seg[i]*σᵢᵢ[i]^3 for i in 1:length(z))
    return val
end

function x0_crit_pure(model::SAFTModel)
    lb_v = lb_volume(model)
    (2.0, log10(lb_v/0.3))
end

function T_scale(model::SAFTModel,z=SA[1.0])
    ϵ = model.params.epsilon.diagvalues
    return prod(ϵ[i]^z[i] for i in 1:length(z))^(1/sum(z))
end

function p_scale(model::SAFTModel,z=SA[1.0])
    ϵ = model.params.epsilon.diagvalues
    σᵢᵢ = model.params.sigma.diagvalues
    val =  sum(z[i]*σᵢᵢ[i]^3/ϵ[i] for i in 1:length(z))*N_A/R̄
    return 1/val
end

function  Δ(model::Union{SAFTModel,CPAModel}, V, T, z)
    κ = model.params.bondvol.values
    Δres = zero_assoc(κ,typeof(V+T+first(z)))
    for (idx,(i,j),(a,b)) in indices(Δres)
        Δres[idx] =@f(Δ,i,j,a,b)
    end
    return Δres
end

function  Δ(model::Union{SAFTModel,CPAModel}, V, T, z,data)
    κ = model.params.bondvol.values
    Δres = zero_assoc(κ,typeof(V+T+first(z)))
    for (idx,(i,j),(a,b)) in indices(Δres)
        Δres[idx] =@f(Δ,i,j,a,b,data)
    end
    return Δres
end
function issite(i,a,ij,ab)
    ia = (i,a)
    i1,i2 = ij
    a1,a2 = ab
    ia1 = (i1,a1)
    ia2 = (i2,a2)
    return (ia == ia1) | (ia == ia2)
end

function complement_index(i,ij)
    i1,i2 = ij
    ifelse(i1 == i,i2,i1)
end

function X!(output,input,pack_indices,delta,modelsites,ρ,z)

    _0 = zero(eltype(output))
    _1 = one(_0)
    Xinput = PackedVofV(pack_indices,input)
    Xoutput = PackedVofV(pack_indices,output)
    n = modelsites.n_sites
    for i ∈ 1:length(modelsites.components) #for i ∈ comps
        sitesᵢ = 1:length(n[i]) #sites are normalized, with independent indices for each component
        iszero(length(sitesᵢ)) && continue
        for a ∈ sitesᵢ #for a ∈ sites(comps(i))
            ∑X = _0
            for (idx,ij,ab) in indices(delta) #iterating for all sites
                if issite(i,a,ij,ab)
                    j = complement_index(i,ij)
                    b = complement_index(a,ab)
                    nj = n[j]
                    njb = nj[b]
                    ∑X += ρ*njb*z[j]*Xinput[j][b]*delta.values[idx]
                end
            end
            rhs = _1/(_1 +∑X)
            Xoutput[i][a] = rhs
        end
    end
    return output
end

function X(model::Union{SAFTModel,CPAModel}, V, T, z,data = nothing)
    _1 = one(V+T+first(z))
    X_ = PackedVectorsOfVectors.packed_ones(typeof(_1),length(@sites(i)) for i ∈ @comps)
    idxs = indices(X_)
    X0 = X_.v
    bv = model.params.bondvol.values
    nn = length(bv.values)
    if isone(nn)
        xia,xjb = X_exact1(model,V,T,z,data)
        i,j = only(bv.outer_indices)
        a,b = only(bv.inner_indices)
        #in the case that i = j, a = b, this does assignment twice
        #we do xia last because xjb is calculated from xia
        X_[j][b] = xjb
        X_[i][a] = xia
        return X_
    end
    ρ = N_A/V
    if data === nothing
        _Δ = @f(Δ)
    else
        _Δ = @f(Δ,data)
    end
    fX(out,in) = X!(out,in,idxs,_Δ,model.sites,ρ,z)

    options = model.assoc_options
    atol = options.atol
    rtol = options.rtol
    max_iters = options.max_iters
    α = options.dampingfactor
    Xsol = Solvers.fixpoint(fX,X0,Solvers.SSFixPoint(α),atol=atol,rtol = rtol,max_iters = max_iters)
    return PackedVofV(idxs,Xsol)
end

#exact calculation of site non-bonded fraction when there is only one site
function X_exact1(model,V,T,z,data=nothing)
    κ = model.params.bondvol.values
    i,j = κ.outer_indices[1]
    a,b = κ.inner_indices[1]

    if data === nothing
        _Δ = @f(Δ,i,j,a,b)
    else
        _Δ = @f(Δ,i,j,a,b,data)
    end
    ρ = N_A/V
    zi = z[i]
    zj = z[j]
    ni = model.sites.n_sites[i]
    na = ni[a]
    nj = model.sites.n_sites[j]
    nb = nj[b]
    ρ = N_A/V
    kia = na*zi*ρ*_Δ
    kjb = nb*zj*ρ*_Δ
    #kia*x*x + x(kjb-kia+1) - 1 = 0
    _a = kia
    _b = 1 -kia + kjb
    _c = -1
    xia = -2*_c/(_b + sqrt(_b*_b - 4*_a*_c))
    xjb = 1/(1+kia*xia)
    return (xia,xjb)
end

function a_assoc(model::Union{SAFTModel,CPAModel}, V, T, z,data=nothing)
    _0 = zero(V+T+first(z))
    nn = length(model.params.bondvol.values.values)
    iszero(nn) && return _0
    X_ = @f(X,data)

    return @f(_a_assoc,X_)
end

function _a_assoc(model::Union{SAFTModel,CPAModel}, V, T, z,X_)
    _0 = zero(first(X_.v))
    n = model.sites.n_sites
    res = _0
    resᵢₐ = _0
    for i ∈ @comps
        ni = n[i]
        iszero(length(ni)) && continue
        Xᵢ = X_[i]
        resᵢₐ = _0
        for (a,nᵢₐ) ∈ pairs(ni)
            Xᵢₐ = Xᵢ[a]
            nᵢₐ = ni[a]
            resᵢₐ +=  nᵢₐ* (log(Xᵢₐ) - Xᵢₐ/2 + 0.5)
        end
        res += resᵢₐ*z[i]
    end
    return res/sum(z)

end

function show_param(param)
    for i in fieldnames(typeof(param))
        display(getfield(param,i))
    end
end

#res =  ∑(z[i]*∑(n[i][a] * (log(X_[i][a]) - X_[i][a]/2 + 0.5) for a ∈ @sites(i)) for i ∈ @comps)/sum(z)

#=
on one site:
Xia = 1/(1+*nb*z[j]*rho*Δ*Xjb)
Xjb = 1/(1+*na*z[i]*rho*Δ*Xia)

kia = na*z[i]*rho*Δ
kjb = nb*z[j]*rho*Δ

Xia = 1/(1+kjb*Xjb)
Xjb = 1/(1+kia*Xia)

Xia = 1/(1+kjb*(1/(1+kia*Xia)))
Xia = 1/(1+kjb/(1+kia*Xia))
Xia = 1/((1+kia*Xia+kjb)/(1+kia*Xia))
Xia = (1+kia*Xia)/(1+kia*Xia+kjb)
Xia*(1+kia*Xia+kjb) = 1+kia*Xia #x = Xia
x*(1+kia*x+kjb) = 1+kia*x
x + kia*x*x + kjb*x - 1 - kia*x = 0
kia*x*x + x(kjb-kia+1) - 1 = 0
x = - (kjb-kia+1) +

x = 1/1+kiax
x(1+kx) - 1 = 0
kx2 +x - 1 = 0
end
=#


text:
export changeParam

# Support for handling complex datatypes in HDF5 files
function writeComplexArray(file, dataset, A::AbstractArray{Complex{T},D}) where {T,D}
  d_type_compound = HDF5.h5t_create(HDF5.H5T_COMPOUND,2*sizeof(T))
  HDF5.h5t_insert(d_type_compound, "r", 0 , HDF5.hdf5_type_id(T))
  HDF5.h5t_insert(d_type_compound, "i", sizeof(T) , HDF5.hdf5_type_id(T))

  shape = collect(reverse(size(A)))
  space = HDF5.h5s_create_simple(D, shape, shape)

  dset_compound = HDF5.h5d_create(file, dataset, d_type_compound, space,
                                  HDF5.H5P_DEFAULT,HDF5.H5P_DEFAULT,HDF5.H5P_DEFAULT)
  HDF5.h5s_close(space)

  HDF5.h5d_write(dset_compound, d_type_compound, HDF5.H5S_ALL, HDF5.H5S_ALL, HDF5.H5P_DEFAULT, A)

  HDF5.h5d_close(dset_compound)
  HDF5.h5t_close(d_type_compound)
end

function isComplexArray(file, dataset)
  if eltype(file[dataset]) <: Complex
    return true
  else
    return false
  end
end

function getComplexType(file, dataset)
  T = HDF5.get_jl_type(
            HDF5.Datatype(
              HDF5.h5t_get_member_type( datatype(file[dataset]).id, 0 )
          )
        )
    return Complex{T}
end

function readComplexArray(file::HDF5.File, dataset)
  T = getComplexType(file, dataset)
  A = copy(HDF5.readmmap(file[dataset],getComplexType(file,dataset)))
  return A
end

function readComplexArray(filename::AbstractString, dataset)
  h5open(filename, "r") do file
    return readComplexArray(file, dataset)
  end
end

function changeParam(filename::AbstractString, paramName::AbstractString, paramValue)
  #GC.gc()
  h5open(filename, "r+") do file
	  if haskey(file, paramName)
      delete_object(file, paramName)
    end
    write(file, paramName, paramValue)
  end
end

function makeAxisArray(array::Array{T,5}, pixelspacing, offset, dt) where T
  N = size(array)
  im = AxisArray(array, Axis{:color}(1:N[1]),
		 Axis{:x}(range(offset[1],step=pixelspacing[1],length=N[2])),
		 Axis{:y}(range(offset[2],step=pixelspacing[2],length=N[3])),
		 Axis{:z}(range(offset[3],step=pixelspacing[3],length=N[4])),
		 Axis{:time}(range(0*unit(dt),step=dt,length=N[5])))
  return im
end

# RawFile removed see commit https://github.com/MagneticParticleImaging/MPIFiles.jl/commit/c2b49833cc00127770fb2e1f5342883f213ff002

"Macro for making required fields return `missing` instead of a key error."
macro keyrequired(expr)
  return quote
    try
      $(esc(expr))
    catch e
      if isa(e, KeyError)
        missing
      else
        rethrow()
      end
    end
  end
end

"Macro for making optional fields return `nothing` instead of a key error."
macro keyoptional(expr)
  return quote
    try
      $(esc(expr))
    catch e
      if isa(e, KeyError)
        nothing
      else
        rethrow()
      end
    end
  end
end


text:
using BinDeps

@BinDeps.setup

#deps = [
	rsvg = library_dependency("xcb", aliases = ["libxcb"])
#]
@BinDeps.install [:xcb => :_jl_libxcb ]


text:
#=
datasets.jl:
- Julia version:
- Author: Jirik
- Date: 2019-03-03
=#

function random_image(shape, obj_min, obj_max, level)
    # print("Random image start")
    im = rand(shape...)

    for i=obj_min[1]:obj_max[1]
        for j=obj_min[2]:obj_max[2]
            for k=obj_min[3]:obj_max[3]
                im[i,j,k] = im[i,j,k] + level
            end
        end
    end
    return im
end

function generate_slope(data_size)
    data = zeros(Int8, data_size[1], data_size[2], data_size[3])
    for i=1:data_size[1]
        for j=1:data_size[2]
            for k=1:data_size[3]
                if ((j - k) < i) & ((j - k + 3) >= i)
                    data[i,j,k] = 1
                end
#                 if i > (j - k)
#                 if (k - j) < i
#                     data[i,j,k] = 1
#                 end
            end
        end
    end
    return data
end

function generate_truncated_sphere(r=10, data_size=nothing)
	if data_size == nothing
		data_size = [Int(r*2 + 3), Int(r*2 + 3), Int(r*2 + 3)]
	end
	center = data_size / 2
    data = zeros(Int8, data_size[1], data_size[2], data_size[3])
    for i=1:data_size[1]
        for j=1:data_size[2]
            for k=1:data_size[3]
				coords = [i, j, k]
                if ((j - k) < i) # & ((j - k + 3) >= i)
					dist = sum((center - coords).^2)^0.5
					if dist < r
                    	data[i,j,k] = 1
					end
                end
#                 if i > (j - k)
#                 if (k - j) < i
#                     data[i,j,k] = 1
#                 end
            end
        end
    end
    return data
end

function data234()
    @info "data234 init"
    data_size = [2,3,4]
    segmentation = zeros(Int8, 2, 3, 4)
    segmentation[1:2,2:3,3:4] .= 1
    return segmentation
end

function tetris_brick()
    @info "data tetris init"
    segmentation = zeros(Int8, 3, 4, 5)
    segmentation[2,2,2:4] .= 1
    segmentation[2,3,3] = 1
    return segmentation
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0
#     language: julia
#     name: julia-1.7
# ---

# %%
function f(x)
    (; a, b, c, d) = x
    @show a b c d
end

# %%
a = 1
b = 2.0
c = "three"
d = 0x04

nt = (; a, b, c, d)
@show nt;

# %%
f(nt);

# %%
function f(x)
    (; a, b, c, d) = x
    @show a b c d
end

# %%
Base.@kwdef struct Foo{A, B, C, D}
    a::A = 1
    b::B = 2.0
    c::C = "three"
    d::D = 0x04
end

foo = Foo()
@show foo;

# %%
f(foo);

# %%
function f(x)
    (; a, b, c, d) = x
    @show a b c d
end

# %%
using DataFrames

a = collect(1:5)
b = collect(1.0:5.0)
c = ["one", "two", "three", "four", "five"]
d = collect(0x01:0x05)

df = DataFrame(; a, b, c, d)

# %%
f(df);

# %%
@show nt;

# %%
typeof(nt)

# %%
dump(nt)

# %%
@show foo;

# %%
dump(foo)

# %%
g(x::NamedTuple{(:a, :b, :c, :d)}) = @show x
g(x::NamedTuple{(:a, :b, :c)}) = @show x

# %%
g(nt)

# %%
@which g(nt)

# %%
nt2 = (; a = 10, b = 20.0, c = "thirty")

# %%
g(nt2)

# %%
@which g(nt2)

# %%

max_stars_repo_path: 0021/unpack with named tuple style.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0021/unpack with named tuple style.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0021/unpack with named tuple style.jl

max_forks_repo_name: genkuroki/public


text:
@testset "lrsnash.jl" begin

    @testset "test 3 by 2 non-degenerate normal form game(Int)" begin
        g = NormalFormGame(Player([3 3; 2 5; 0 6]),
                           Player([3 2 3; 2 6 1]))
        NEs = [([1//1, 0//1, 0//1], [1//1, 0//1]),
               ([4//5, 1//5, 0//1], [2//3, 1//3]),
               ([0//1, 1//3, 2//3], [1//3, 2//3])]
        NEs_computed = @inferred(lrsnash(g))

        @test sort(NEs_computed) == sort(NEs)
    end

    @testset "test 3 by 2 non-degenerate normal form game(Rational)" begin
        g = NormalFormGame(Player([3//1 3//1; 2//1 5//1; 0//1 6//1]),
                           Player([3//1 2//1 3//1; 2//1 6//1 1//1]))
        NEs = [([1//1, 0//1, 0//1], [1//1, 0//1]),
               ([4//5, 1//5, 0//1], [2//3, 1//3]),
               ([0//1, 1//3, 2//3], [1//3, 2//3])]
        NEs_computed = @inferred(lrsnash(g))

        @test sort(NEs_computed) == sort(NEs)
    end

    @testset "test 3 by 2 degenerate normal form game(Int)" begin
        g = NormalFormGame(Player([1 -1; -1 1; 0 0]),
                           Player([1 0 0; 0 0 0]))
        NEs = [([1, 0, 0], [1, 0]),
               ([0, 1, 0], [0, 1]),
               ([0, 1, 0], [1//2, 1//2]),
               ([0, 0, 1], [1//2, 1//2])]
        NEs_computed = @inferred(lrsnash(g))

        @test sort(NEs_computed) == sort(NEs)
    end

end


text:
using Test
using SumOfSquares
using DynamicPolynomials

function quadratic_test(
    optimizer, config::MOIT.TestConfig,
    cone::SumOfSquares.PolyJuMP.PolynomialSet, bivariate::Bool)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @variable(model, α)

    @polyvar x
    if bivariate
        @polyvar y
        poly = x^2 + α*x*y + y^2
        cert_monos = [x, y]
        monos = [x^2, x*y, y^2]
    else
        poly = x^2 + α*x + 1
        cert_monos = [x, 1]
        monos = [x^2, x, 1]
    end
    cref = @constraint(model, poly in cone)

    # See https://github.com/JuliaOpt/MathOptInterface.jl/issues/676
    @objective(model, Max, α + 1)
    optimize!(model)

    @test certificate_monomials(cref) == cert_monos

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 3.0 atol=atol rtol=rtol

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(α) ≈ 2.0 atol=atol rtol=rtol

    @test_throws SumOfSquares.ValueNotSupported value(cref)
    p = gram_matrix(cref)
    @test getmat(p) ≈ ones(2, 2) atol=atol rtol=rtol
    @test p.x == cert_monos

    a = moment_value.(moments(dual(cref)))
    @test a[2] ≈ -1.0 atol=atol rtol=rtol
    @test a[1] + a[3] ≈ 2.0 atol=atol rtol=rtol

    @test dual_status(model) == MOI.FEASIBLE_POINT
    for μ in [dual(cref), moments(cref)]
        @test μ isa AbstractMeasure{Float64}
        @test length(moments(μ)) == 3
        @test a ≈ moment_value.(moments(μ)) atol=atol rtol=rtol
        @test monomial.(moments(μ)) == monos
    end

    ν = moment_matrix(cref)
    @test getmat(ν) ≈ [a[1] a[2]
                       a[2] a[3]] atol=atol rtol=rtol
    @test ν.x == cert_monos

    S = SumOfSquares.SOSPolynomialSet{
        SumOfSquares.FullSpace, Monomial{true}, MonomialVector{true}, SumOfSquares.Certificate.Remainder{typeof(cone), SumOfSquares.MonomialBasis, Tuple{}}
    }
    @test list_of_constraint_types(model) == [(Vector{AffExpr}, S)]
    test_delete_bridge(
        model, cref, 1,
        ((MOI.VectorOfVariables, MOI.Nonnegatives, 0),
         (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
         (MOI.VectorAffineFunction{Float64},
          SumOfSquares.PolyJuMP.ZeroPolynomialSet{
              SumOfSquares.FullSpace, SumOfSquares.MonomialBasis,
              Monomial{true}, MonomialVector{true}},
          0)))
end
sos_univariate_quadratic_test(optimizer, config)   = quadratic_test(optimizer, config, SOSCone(), false)
sd_tests["sos_univariate_quadratic"] = sos_univariate_quadratic_test
sos_bivariate_quadratic_test(optimizer, config)    = quadratic_test(optimizer, config, SOSCone(), true)
sd_tests["sos_bivariate_quadratic"] = sos_bivariate_quadratic_test
sdsos_univariate_quadratic_test(optimizer, config) = quadratic_test(optimizer, config, SDSOSCone(), false)
soc_tests["sdsos_univariate_quadratic"] = sdsos_univariate_quadratic_test
sdsos_bivariate_quadratic_test(optimizer, config)  = quadratic_test(optimizer, config, SDSOSCone(), true)
soc_tests["sdsos_bivariate_quadratic"] = sdsos_bivariate_quadratic_test
dsos_univariate_quadratic_test(optimizer, config)  = quadratic_test(optimizer, config, DSOSCone(), false)
linear_tests["dsos_univariate_quadratic"] = dsos_univariate_quadratic_test
dsos_bivariate_quadratic_test(optimizer, config)   = quadratic_test(optimizer, config, DSOSCone(), true)
linear_tests["dsos_bivariate_quadratic"] = dsos_bivariate_quadratic_test


text:
mutable struct Node
    camera::Union{Nothing,Int}
    children::Union{Nothing,Vector{Int}}
    skin::Union{Nothing,Int}
    matrix::Union{Nothing,Vector{Cfloat}}
    mesh::Union{Nothing,Int}
    rotation::Union{Nothing,Vector{Cfloat}}
    scale::Union{Nothing,Vector{Cfloat}}
    translation::Union{Nothing,Vector{Cfloat}}
    weights::Union{Nothing,Vector{Cfloat}}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Node()
        obj = new()
        obj.camera = nothing
        obj.children = nothing
        obj.skin = nothing
        obj.matrix = nothing
        obj.mesh = nothing
        obj.rotation = nothing
        obj.scale = nothing
        obj.translation = nothing
        obj.weights = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.getproperty(obj::Node, sym::Symbol)
    x = getfield(obj, sym)
    sym === :matrix && x === nothing && return Cfloat[1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1]
    sym === :rotation && x === nothing && return Cfloat[0,0,0,1]
    sym === :scale && x === nothing && return Cfloat[1,1,1]
    sym === :translation && x === nothing && return Cfloat[0,0,0]
    return x
end

function Base.setproperty!(obj::Node, sym::Symbol, x)
    if sym === :camera && x !== nothing
        x ≥ 0 || throw(ArgumentError("the index of the camera referenced by this node should be ≥ 0"))
    elseif sym === :children && x !== nothing
        all(c->c ≥ 0, x) || throw(ArgumentError("the indices of this node's children should be ≥ 0"))
        allunique(x) || throw(ArgumentError("children must be unique."))
    elseif sym === :skin && x !== nothing
        x ≥ 0 || throw(ArgumentError("the index of the skin referenced by this node should be ≥ 0"))
    elseif sym === :matrix && x !== nothing
        length(x) == 16 || throw(ArgumentError("the length of matrix should be exactly 16"))
    elseif sym === :mesh && x !== nothing
        x ≥ 0 || throw(ArgumentError("the index of the mesh in this node should be ≥ 0"))
    elseif sym === :rotation && x !== nothing
        length(x) == 4 || throw(ArgumentError("the length of rotation should be exactly 4"))
        all(-1 .≤ x .≤ 1) || throw(ArgumentError("each element in the unit quaternion rotation should be ≥ -1 and ≤ 1"))
    elseif sym === :scale && x !== nothing
        length(x) == 3 || throw(ArgumentError("the length of scale should be exactly 3"))
    elseif sym === :translation && x !== nothing
        length(x) == 3 || throw(ArgumentError("the length of translation should be exactly 3"))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Node}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Node}) = (:camera, :children, :skin, :matrix, :mesh, :rotation, :scale, :translation, :weights, :name, :extensions, :extras)


text:
# TODO: Make uppercase
# ids related with experiments performed by Rath2017
const EXPS = ["A", "B", "C", "D", "E", "F01","F02","F03"]
# experiments that reached steady states, ordered by ξ value
const STSTS = ["E", "D", "F01", "A", "B", "C"]
# measured mets
# Mets which concentration and exchange was measured
const MSD_METS = ["GLC", "LAC", "GLN", "NH4", "GAL","PYR","GLU","ALA","ASP"]
# ALL_METS
# Counting the feed medium mets and the measured mets
const ALL_METS = ["LAC", "NH4", "HCO3", "ASN", "MET",
    "GLC", "GLY", "TYR", "LEU", "PHE", "GLU", "ASP",
    "GLN", "CYS", "ILE", "PYR", "H", "LYS", "TRP", "SER",
    "HIS", "VAL", "ALA", "ARG", "PRO", "GAL", "THR"]
# exchange rates
const RATE_Q = ["q$met" for met in MSD_METS]
# feed concentrations
const FEED_C = ["cGLC", "cGAL", "cGLN"]
# medium concentration
const MED_S = ["s$met" for met in MSD_METS];
const GROWTH_IDER = "μ"