text:
"""
    CartesianIndexHalfInt(i, j, k...) -> I
    CartesianIndexHalfInt((i, j, k...)) -> I

Create a multidimensional index `I`, which can be used for
indexing a multidimensional `AbstractHalfIntegerArray`.
In particular, for an array `A`, the operation `A[I]` is
equivalent to `A[i,j,k...]`.  One can freely mix integer, half-integer and
`CartesianIndex` indices; for example, `A[Ipre, i, Ipost]` (where
`Ipre` and `Ipost` are `CartesianIndex` indices and `i` is an
integer or a half-integer)
can be a useful expression when writing algorithms that
work along a single dimension of an array of arbitrary dimensionality.

# Examples
```jldoctest
julia> h = HalfIntArray(reshape(1:4, 2, 2), -1//2:1//2, -1//2:1//2)
2×2 HalfIntArray(reshape(::UnitRange{Int64}, 2, 2), -1/2:1/2, -1/2:1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
 1  3
 2  4

julia> h[CartesianIndexHalfInt(-1//2, 1//2)]
3

julia> h[CartesianIndexHalfInt(-1//2), 1//2]
3
```
"""
struct CartesianIndexHalfInt{N} <: Base.AbstractCartesianIndex{N}
    parent :: CartesianIndex{N}
    offsets :: NTuple{N,HalfInt}
end

Base.parent(c::CartesianIndexHalfInt) = c.parent

CartesianIndexHalfInt(::Tuple{}) = CartesianIndexHalfInt{0}(())
function CartesianIndexHalfInt(index::NTuple{N,Int},
    offsets::NTuple{N,HalfInt} = Base.fill_to_length((),zero(HalfInt),Val(N))) where {N}

    CartesianIndexHalfInt{N}(index, offsets)
end
function CartesianIndexHalfInt{N}(index::NTuple{N,Int},
    offsets::NTuple{N,HalfInt} = Base.fill_to_length((),zero(HalfInt),Val(N))) where {N}

    CartesianIndexHalfInt{N}(CartesianIndex(index), offsets)
end
function CartesianIndexHalfInt(offsets::NTuple{N,HalfInt}) where {N}
    CartesianIndexHalfInt{N}(offsets)
end
function CartesianIndexHalfInt{N}(offsets::NTuple{N,HalfInt}) where {N}
    index = Base.fill_to_length((),0,Val(N))
    CartesianIndexHalfInt{N}(CartesianIndex(index), offsets)
end
# Specifying Ints populates the parent CaratesianIndex
CartesianIndexHalfInt() = CartesianIndexHalfInt{0}((),())
CartesianIndexHalfInt(index::Int...) = CartesianIndexHalfInt(index)
# Specifying HalfInts populates the offsets
CartesianIndexHalfInt(offsets::HalfInt...) = CartesianIndexHalfInt(offsets)
CartesianIndexHalfInt(offsets::Real...) = CartesianIndexHalfInt(HalfInt.(offsets))
# # Allow passing tuples smaller than N
CartesianIndexHalfInt{0}(index::Tuple{}) = CartesianIndexHalfInt{0}((),())
CartesianIndexHalfInt{N}(index::Tuple{}) where {N} = CartesianIndexHalfInt{N}(Base.fill_to_length(index, 1, Val(N)))
CartesianIndexHalfInt{N}(index::Tuple{Vararg{Int}}) where {N} = CartesianIndexHalfInt{N}(Base.fill_to_length(index, 1, Val(N)))
CartesianIndexHalfInt{N}(index::Integer...) where {N} = CartesianIndexHalfInt{N}(index)
CartesianIndexHalfInt{N}(offsets::Tuple{Vararg{HalfInt}}) where {N} = CartesianIndexHalfInt{N}(Base.fill_to_length(offsets, one(HalfInt), Val(N)))
CartesianIndexHalfInt{N}(offsets::Tuple{Vararg{Real}}) where {N} = CartesianIndexHalfInt{N}(HalfInt.(offsets))
CartesianIndexHalfInt{N}(offsets::HalfInt...) where {N} = CartesianIndexHalfInt{N}(offsets)
CartesianIndexHalfInt{N}(offsets::Real...) where {N} = CartesianIndexHalfInt{N}(HalfInt.(offsets))
CartesianIndexHalfInt{N}() where {N} = CartesianIndexHalfInt{N}(())
# Un-nest passed CartesianIndexes
CartesianIndexHalfInt(index::Union{Integer, HalfInteger, CartesianIndexHalfInt}...) = CartesianIndexHalfInt(index)
CartesianIndexHalfInt(index::Tuple{Vararg{Union{Integer, HalfInteger, CartesianIndexHalfInt}}}) = CartesianIndexHalfInt(flatten(index))

flatten(I::Tuple{}) = I
flatten(I::Tuple{Any}) = I
flatten(I::Tuple{HalfIntOrInt}) = (HalfInt(I[1]),)
flatten(I::Tuple{<:CartesianIndexHalfInt}) = Tuple(I[1])
@inline flatten(I) = _flatten(I...)
@inline _flatten() = ()
@inline _flatten(i, I...) = (HalfInt(i), _flatten(I...)...)
@inline _flatten(i::CartesianIndexHalfInt, I...) = (Tuple(i)..., _flatten(I...)...)

# length
Base.length(::CartesianIndexHalfInt{N}) where {N} = N

# indexing
Base.getindex(index::CartesianIndexHalfInt, i::Integer) = parent(index)[i] .+ index.offsets[i]
Base.eltype(::CartesianIndexHalfInt) = HalfInt

# access to index tuple
Base.Tuple(index::CartesianIndexHalfInt) = map(+,Tuple(parent(index)), index.offsets)

# equality
Base.:(==)(a::CartesianIndexHalfInt{N}, b::CartesianIndexHalfInt{N}) where N = Tuple(a) == Tuple(b)

# zeros and ones
Base.zero(::CartesianIndexHalfInt{N}) where {N} = zero(CartesianIndexHalfInt{N})
Base.zero(::Type{CartesianIndexHalfInt{N}}) where {N} = CartesianIndexHalfInt(ntuple(x -> 0, Val(N)))
Base.oneunit(::CartesianIndexHalfInt{N}) where {N} = oneunit(CartesianIndexHalfInt{N})
Base.oneunit(::Type{CartesianIndexHalfInt{N}}) where {N} = CartesianIndexHalfInt(ntuple(x -> 1, Val(N)))

# arithmetic, min/max
@inline Base.:(-)(index::CartesianIndexHalfInt{N}) where {N} =
    CartesianIndexHalfInt{N}(-parent(index), map(-,index.offsets))
@inline Base.:(+)(index1::CartesianIndexHalfInt{N}, index2::CartesianIndexHalfInt{N}) where {N} =
    CartesianIndexHalfInt{N}(parent(index1) + parent(index2), map(+, index1.offsets, index2.offsets))
@inline Base.:(-)(index1::CartesianIndexHalfInt{N}, index2::CartesianIndexHalfInt{N}) where {N} =
    CartesianIndexHalfInt{N}(parent(index1) - parent(index2), map(-, index1.offsets, index2.offsets))

for f in [:min, :max]
    @eval @inline function Base.$f(index1::CartesianIndexHalfInt{N}, index2::CartesianIndexHalfInt{N}) where {N}
        index = map($f, Tuple(parent(index1)), Tuple(parent(index2)))
        of = map($f, index1.offsets, index2.offsets)
        CartesianIndexHalfInt{N}(index, of)
    end
end

@inline Base.:(*)(a::Integer, index::CartesianIndexHalfInt{N}) where {N} = CartesianIndexHalfInt{N}(a*parent(index), map(x->a*x, index.offsets))
@inline Base.:(*)(index::CartesianIndexHalfInt, a::Integer) = *(a,index)

# comparison
@inline Base.isless(I1::CartesianIndexHalfInt{N}, I2::CartesianIndexHalfInt{N}) where {N} = _isless(0, Tuple(I1), Tuple(I2))
@inline function _isless(ret, I1::NTuple{N,HalfInt}, I2::NTuple{N,HalfInt}) where N
    newret = ifelse(ret==0, icmp(I1[N], I2[N]), ret)
    _isless(newret, Base.front(I1), Base.front(I2))
end
_isless(ret, ::Tuple{}, ::Tuple{}) = ifelse(ret==1, true, false)
icmp(a, b) = ifelse(isless(a,b), 1, ifelse(a==b, 0, -1))

# conversions
Base.convert(::Type{T}, index::CartesianIndexHalfInt{1}) where {T<:Number} = convert(T, index[1])
Base.convert(::Type{T}, index::CartesianIndexHalfInt) where {T<:Tuple} = convert(T, Tuple(index))

# hashing
const cartindexhash_seed = UInt == UInt64 ? 0xd60ca92f8284b8b0 : 0xf2ea7c2e
function Base.hash(ci::CartesianIndexHalfInt, h::UInt)
    h += cartindexhash_seed
    h = hash(ci.offsets, h)
    h = hash(parent(ci), h)
    h = hash(:CartesianIndexHalfInt, h)
    return h
end

Base.iterate(::CartesianIndexHalfInt) =
    error("iteration is deliberately unsupported for CartesianIndexHalfInt. Use `I` rather than `I...`, or use `Tuple(I)...`")

"""
    CartesianIndicesHalfInt((istart:istop, jstart:jstop, ...)) -> R

Define a region `R` spanning a multidimensional rectangular range
of integer indices. These are most commonly encountered in the
context of iteration, where `for I in R ... end` will return
[`CartesianIndexHalfInt`](@ref) indices `I` equivalent to the nested loops

    for j = jstart:jstop
        for i = istart:istop
            ...
        end
    end

Consequently these can be useful for writing algorithms that
work in arbitrary dimensions.

A `CartesianIndicesHalfInt` type is equivalent to a `CartesianIndices` type
for integer ranges. The difference is that a `CartesianIndicesHalfInt` allows the
ranges to be half-integer `AbstractUnitRange`s. They are therefore suitable for
indexing into `AbstractHalfIntegerArray`s.

    CartesianIndicesHalfInt(A::AbstractArray) -> R

As a convenience, constructing a `CartesianIndicesHalfInt` from an array makes a
range of its indices.

# Examples
```jldoctest
julia> h = HalfIntArray(reshape(1:4, 2, 2), -1//2:1//2, -1//2:1//2)
2×2 HalfIntArray(reshape(::UnitRange{Int64}, 2, 2), -1/2:1/2, -1/2:1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
 1  3
 2  4

julia> c = CartesianIndicesHalfInt(h)
2×2 CartesianIndicesHalfInt{2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}} with indices -1/2:1/2×-1/2:1/2:
 CartesianIndexHalfInt(-1/2, -1/2)  CartesianIndexHalfInt(-1/2, 1/2)
 CartesianIndexHalfInt(1/2, -1/2)   CartesianIndexHalfInt(1/2, 1/2)

julia> c[1/2,1/2]
CartesianIndexHalfInt(1/2, 1/2)
```
"""
struct CartesianIndicesHalfInt{N,R} <: AbstractHalfIntegerWrapper{CartesianIndexHalfInt{N},N}
    cartinds :: CartesianIndices{N,R}
    offsets :: NTuple{N,HalfInt}
end

Base.parent(iter::CartesianIndicesHalfInt) = CartesianIndicesHalfIntParent(iter.cartinds, iter.offsets)

struct CartesianIndicesHalfIntParent{N,R} <: AbstractArray{CartesianIndexHalfInt{N},N}
    cartinds :: CartesianIndices{N,R}
    offsets :: NTuple{N,HalfInt}
end

Base.parent(iter::CartesianIndicesHalfIntParent) = iter.cartinds

const CartesianIndicesHalfIntAny{N,R} = Union{CartesianIndicesHalfInt{N,R},CartesianIndicesHalfIntParent{N,R}}

_offsets(iter::Tuple{IdOffsetRange,Vararg{IdOffsetRange}}) = (iter[1].offset,_offsets(Base.tail(iter))...)
_offsets(iter::Tuple{IdOffsetRange}) = (iter[1].offset,)
function _offsets(iter::Tuple{AbstractUnitRange{HalfInt},Vararg{AbstractUnitRange{HalfInt}}})
    x = first(iter[1])
    (x - one(x), _offsets(Base.tail(iter))...)
end
_offsets(iter::Tuple{AbstractUnitRange{HalfInt}}) = (x = first(iter[1]); (x - one(x),))
_parents(iter::Tuple{IdOffsetRange,Vararg{IdOffsetRange}}) = (parent(iter[1]),_parents(Base.tail(iter))...)
_parents(iter::Tuple{IdOffsetRange}) = (parent(iter[1]),)
function _parents(iter::Tuple{AbstractUnitRange{HalfInt},Vararg{AbstractUnitRange{HalfInt}}})
    (Base.OneTo(length(iter[1])), _parents(Base.tail(iter))...)
end
_parents(iter::Tuple{AbstractUnitRange{HalfInt}}) = (Base.OneTo(length(iter[1])),)

CartesianIndicesHalfInt() = CartesianIndicesHalfInt(())
function CartesianIndicesHalfInt(::Tuple{})
    CartesianIndicesHalfInt(CartesianIndices(()), ())
end
function CartesianIndicesHalfInt(iter::NTuple{N,AbstractUnitRange{HalfInt}}) where {N}
    of = _offsets(iter)
    c = CartesianIndices(_parents(iter))
    CartesianIndicesHalfInt(c, of)
end
function CartesianIndicesHalfInt(iter::NTuple{N,AbstractUnitRange{Int}},
    of = ntuple(x->zero(HalfInt),Val(N))) where {N}
    CartesianIndicesHalfInt(CartesianIndices(iter), of)
end

# Mix of integer and half-integer axes, promote to half-integer
function CartesianIndicesHalfInt(iter::Tuple{Vararg{AbstractUnitRange}})
    iterIdOfR = map(IdOffsetRange, iter)
    CartesianIndicesHalfInt(iterIdOfR)
end

CartesianIndicesHalfInt(A::AbstractArray) = CartesianIndicesHalfInt(axes(A))

Base.IndexStyle(::Type{<:CartesianIndicesHalfIntAny{N,R}}) where {N,R} = IndexCartesian()

_convertIntOrHalfInt(t::NTuple{N,Int}) where {N} = t
_convertIntOrHalfInt(t::NTuple{N,Real}) where {N} = map(HalfInt, t)
_convertIntOrHalfInt(i::Int) = i
_convertIntOrHalfInt(i::Real) = HalfInt(i)

@propagate_inbounds function Base.getindex(iter::CartesianIndicesHalfInt{N}, I::Real...) where {N}
    @boundscheck checkbounds(iter, I...)
    J = trimtoN(I, Val(N))
    CartesianIndexHalfInt(_convertIntOrHalfInt(J))
end

@propagate_inbounds function Base.getindex(iter::CartesianIndicesHalfIntParent, I::Int...)
    @boundscheck checkbounds(iter, I...)
    CartesianIndexHalfInt(iter.cartinds[I...], iter.offsets)
end

# 1D arrays use Cartesian Indexing
for DT in [:Int, :Real]
    @eval @propagate_inbounds function Base.getindex(iter::CartesianIndicesHalfInt{1}, i::$DT)
        @boundscheck checkbounds(iter, i)
        CartesianIndexHalfInt(_convertIntOrHalfInt(i))
    end
end

# Linear Indexing
@propagate_inbounds function Base.getindex(A::CartesianIndicesHalfInt, i::Real)
    @boundscheck checkbounds(A, i)
    @boundscheck ensureInt(i)
    CartesianIndexHalfInt(A.cartinds[unsafeInt(i)], A.offsets)
end

Base.IteratorSize(::Type{<:CartesianIndicesHalfInt{N}}) where {N} = Base.HasShape{N}()

@inline function _iterate(iter, newstate)
    ind,_ = newstate
    cind = CartesianIndexHalfInt(ind, iter.offsets)
    cind, ind
end
@inline _iterate(iter, ::Nothing) = nothing

@inline function Base.iterate(iter::CartesianIndicesHalfIntAny)
    newstate = iterate(iter.cartinds)
    _iterate(iter, newstate)
end
@inline function Base.iterate(iter::CartesianIndicesHalfIntAny, state)
    newstate = iterate(iter.cartinds, state)
    _iterate(iter, newstate)
end

# 0-d cartesian ranges are special-cased to iterate once and only once
function Base.iterate(iter::CartesianIndicesHalfIntAny{0}, done=false)
    done ? nothing : (CartesianIndexHalfInt(), true)
end

@inline function Base.index_ndims(i1::CartesianIndexHalfInt, I...)
    (map(x->true, parent(i1).I)..., Base.index_ndims(I...)...)
end
@inline function Base.index_ndims(i1::AbstractArray{CartesianIndexHalfInt{N}}, I...) where N
    (ntuple(x->true, Val(N))..., Base.index_ndims(I...)...)
end

for f in [:size, :length]
    @eval @inline function Base.$f(iter::CartesianIndicesHalfIntAny)
        $f(iter.cartinds)
    end
end

Base.axes(c::CartesianIndicesHalfInt) = map(IdOffsetRange, axes(c.cartinds), c.offsets)
Base.axes(c::CartesianIndicesHalfInt, d) = d <= ndims(c) ? IdOffsetRange(axes(c.cartinds, d), c.offsets[d]) : IdOffsetRange(axes(c.cartinds, d))

for f in [:first, :last]
    @eval @inline function Base.$f(iter::CartesianIndicesHalfIntAny)
        CartesianIndexHalfInt($f(iter.cartinds), iter.offsets)
    end
end

@inline function Base.in(i::CartesianIndexHalfInt{N}, r::CartesianIndicesHalfInt{N}) where {N}
    _in(true, Tuple(i), Tuple(first(r)), Tuple(last(r)))
end
_in(b, ::Tuple{}, ::Tuple{}, ::Tuple{}) = b
@inline function _in(b, i, start, stop)
    _in(b & (start[1] <= i[1] <= stop[1]),
        Base.tail(i), Base.tail(start), Base.tail(stop))
end

function Base.hash(c::CartesianIndicesHalfInt, h::UInt)
    h = hash(c.offsets, h)
    h = hash(c.cartinds, h)
    h = hash(:CartesianIndicesHalfInt, h)
    return h
end

@inline Base.to_indices(A, I::Tuple{Vararg{Union{HalfInteger, CartesianIndexHalfInt}}}) = to_indices(A, axes(A), I)
@inline Base.to_indices(A::AbstractHalfIntegerArrayOrWrapper, inds::Tuple{AbstractRange{Int},Vararg{Any}}, I::Tuple{CartesianIndexHalfInt, Vararg{Any}}) =
    to_indices(A, inds, (Tuple(first(I))..., Base.tail(I)...))
@inline Base.to_indices(A::AbstractHalfIntegerArrayOrWrapper, inds, I::Tuple{CartesianIndexHalfInt, Vararg{Any}}) =
    to_indices(A, inds, (Tuple(first(I))..., Base.tail(I)...))

# In general assume that arrays have integer indices
@inline Base.to_indices(A::AbstractArray, inds::Tuple{AbstractRange{Int},Vararg{Any}}, I::Tuple{CartesianIndexHalfInt, Vararg{Any}}) =
    to_indices(A, inds, (map(Int, Tuple(first(I)))..., Base.tail(I)...))

# But for arrays of CartesianIndex, we just skip the appropriate number of inds
@inline function Base.to_indices(A, inds, I::Tuple{AbstractArray{CartesianIndexHalfInt{N}}, Vararg{Any}}) where N
    _, indstail = Base.IteratorsMD.split(inds, Val(N))
    (Base.to_index(A, first(I)), to_indices(A, indstail, Base.tail(I))...)
end

function Base.checkbounds(::Type{Bool}, A::AbstractArray, i::Union{CartesianIndexHalfInt, AbstractArray{<:CartesianIndexHalfInt}})
    Base.checkbounds_indices(Bool, axes(A), (i,))
end

@inline Base.checkbounds_indices(::Type{Bool}, ::Tuple{}, I::Tuple{CartesianIndexHalfInt,Vararg{Any}}) =
    Base.checkbounds_indices(Bool, (), (Tuple(first(I))..., Base.tail(I)...))
@inline Base.checkbounds_indices(::Type{Bool}, IA::Tuple{Any}, I::Tuple{CartesianIndexHalfInt,Vararg{Any}}) =
    Base.checkbounds_indices(Bool, IA, (Tuple(first(I))..., Base.tail(I)...))
@inline Base.checkbounds_indices(::Type{Bool}, IA::Tuple, I::Tuple{CartesianIndexHalfInt,Vararg{Any}}) =
    Base.checkbounds_indices(Bool, IA, (Tuple(first(I))..., Base.tail(I)...))

# Support indexing with an array of CartesianIndex{N}s
# Here we try to consume N of the indices (if there are that many available)
# The first two simply handle ambiguities
@inline function Base.checkbounds_indices(::Type{Bool}, ::Tuple{},
        I::Tuple{AbstractArray{CartesianIndexHalfInt{N}},Vararg{Any}}) where N
    checkindex(Bool, (), I[1]) & Base.checkbounds_indices(Bool, (), Base.tail(I))
end
@inline function Base.checkbounds_indices(::Type{Bool}, IA::Tuple{Any},
        I::Tuple{AbstractArray{CartesianIndexHalfInt{0}},Vararg{Any}})
    Base.checkbounds_indices(Bool, IA, Base.tail(I))
end
@inline function Base.checkbounds_indices(::Type{Bool}, IA::Tuple{Any},
        I::Tuple{AbstractArray{CartesianIndexHalfInt{N}},Vararg{Any}}) where N
    checkindex(Bool, IA, I[1]) & Base.checkbounds_indices(Bool, (), Base.tail(I))
end
@inline function Base.checkbounds_indices(::Type{Bool}, IA::Tuple,
        I::Tuple{AbstractArray{CartesianIndexHalfInt{N}},Vararg{Any}}) where N

    IA1, IArest = Base.IteratorsMD.split(IA, Val(N))
    checkindex(Bool, IA1, I[1]) & Base.checkbounds_indices(Bool, IArest, Base.tail(I))
end

function Base.checkindex(::Type{Bool}, inds::Tuple, I::AbstractArray{<:CartesianIndexHalfInt})
    b = true
    for i in I
        b &= Base.checkbounds_indices(Bool, inds, (i,))
    end
    b
end

# It might be possible to construct a LinearIndices range for integer indices
function Base.CartesianIndices(A::AbstractHalfIntegerArrayOrWrapper)
    CartesianIndices(IdentityUnitRange.(UnitRange{Int}.(axes(A))))
end

"""
    LinearIndicesHalfInt(A::AbstractArray)

Return a `LinearIndicesHalfInt` array with the same shape and `axes` as `A`,
holding the linear index of each entry in `A`. Indexing this array with
cartesian indices allows mapping them to linear indices.

For arrays with conventional indexing (indices start at 1), or any multidimensional
array, linear indices range from 1 to `length(A)`. However, for `AbstractVector`s
linear indices are `axes(A, 1)`, and therefore do not start at 1 for vectors with
unconventional indexing.

Calling this function is the "safe" way to write algorithms that
exploit linear indexing.

`LinearIndicesHalfInt` are equivalent to `LinearIndices` for arrays with integer axes,
however `LinearIndicesHalfInt` also support arrays with half-integer axes. They are
therefore the suitable for working with `AbstractHalfIntegerArray`s.

# Examples
```jldoctest
julia> h = HalfIntArray(reshape(1:4, 2, 2), -1//2:1//2, -1//2:1//2)
2×2 HalfIntArray(reshape(::UnitRange{Int64}, 2, 2), -1/2:1/2, -1/2:1/2) with eltype Int64 with indices -1/2:1/2×-1/2:1/2:
 1  3
 2  4

julia> linds = LinearIndicesHalfInt(h)
2×2 LinearIndicesHalfInt{2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}} with indices -1/2:1/2×-1/2:1/2:
 1  3
 2  4

julia> linds[1/2, 1/2]
4

julia> linds[CartesianIndexHalfInt(1/2, 1/2)]
4

julia> v = HalfIntArray([1,2,3], -1:1)
3-element HalfIntArray(::Array{Int64,1}, -1:1) with eltype Int64 with indices -1:1:
 1
 2
 3

julia> lindsv = LinearIndicesHalfInt(v)
3-element LinearIndicesHalfInt{1,Tuple{Base.OneTo{Int64}}} with indices -1:1:
 -1
  0
  1

julia> lindsv[0]
0
```
"""
struct LinearIndicesHalfInt{N,R<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractHalfIntegerWrapper{HalfInt,N}
    lininds :: LinearIndices{N,R}
    offsets :: NTuple{N,HalfInt}
end

Base.parent(l::LinearIndicesHalfInt) = LinearIndicesHalfIntParent(l.lininds, l.offsets)

struct LinearIndicesHalfIntParent{N,R<:NTuple{N,AbstractUnitRange{Int}}} <: AbstractArray{HalfInt,N}
    lininds :: LinearIndices{N,R}
    offsets :: NTuple{N,HalfInt}
end

Base.parent(l::LinearIndicesHalfIntParent) = l.lininds

function LinearIndicesHalfInt(iter::LinearIndices{N,R},
    offsets = first.(iter.indices) .- 1) where {N,R}

    LinearIndicesHalfInt{N,R}(iter, offsets)
end

function LinearIndicesHalfInt(inds::NTuple{N,AbstractUnitRange{Int}},
    offsets = Base.fill_to_length((),zero(HalfInt),Val(N))) where {N}

    LinearIndicesHalfInt(LinearIndices(inds), offsets)
end

function LinearIndicesHalfInt(inds::Tuple{Vararg{IdOffsetRange}})
    LinearIndicesHalfInt(parent.(inds), offset.(inds))
end

# Mix of integer and half-integer axes, promote to half-integer
function LinearIndicesHalfInt(iter::Tuple{Vararg{AbstractUnitRange}})
    iterIdOfR = map(IdOffsetRange, iter)
    LinearIndicesHalfInt(iterIdOfR)
end

LinearIndicesHalfInt(::Tuple{}) = LinearIndicesHalfInt{0,Tuple{}}(LinearIndices(()),())

LinearIndicesHalfInt(A::AbstractArray) = LinearIndicesHalfInt(axes(A))

# AbstractArray implementation
Base.IndexStyle(::Type{<:LinearIndicesHalfInt}) = IndexLinear()
Base.axes(iter::LinearIndicesHalfInt) = map(IdOffsetRange,axes(iter.lininds), iter.offsets)
Base.axes(iter::LinearIndicesHalfInt, d) = d <= ndims(iter) ? IdOffsetRange(axes(iter.lininds,d), iter.offsets[d]) : IdOffsetRange(axes(iter.lininds, d))
for f in [:length, :size]
    @eval function Base.$f(iter::Union{LinearIndicesHalfInt, LinearIndicesHalfIntParent})
        $f(iter.lininds)
    end
end

# Linear indexing
function Base.getindex(iter::LinearIndicesHalfInt, i::Real)
    @boundscheck checkbounds(iter, i)
    ensureInt(i)
    unsafeInt(i)
end

# Cartesian indexing
@propagate_inbounds function Base.getindex(iter::LinearIndicesHalfInt, I::Real...)
    @boundscheck checkbounds(iter, I...)
    J = to_parentindices(axes(iter),I)
    iter.lininds[J...]
end

function Base.getindex(iter::LinearIndicesHalfIntParent, i::Int)
    @boundscheck checkbounds(iter, i)
    i
end
@propagate_inbounds function Base.getindex(iter::LinearIndicesHalfIntParent, I::Int...)
    @boundscheck checkbounds(iter, I...)
    iter.lininds[I...]
end

@inline function Base.getindex(iter::LinearIndicesHalfInt{1}, i::Real)
    @boundscheck checkbounds(iter, i)
    i
end
@inline function Base.getindex(iter::LinearIndicesHalfInt{1}, i::Real, I::Real...)
    @boundscheck checkbounds(iter, i, I...)
    @inbounds iter[i]
end
@inline function Base.getindex(iter::LinearIndicesHalfIntParent{1}, i::Int)
    @boundscheck checkbounds(iter, i)
    @inbounds i + iter.offsets[1]
end
@inline function Base.getindex(iter::LinearIndicesHalfIntParent{1}, i::Int, I::Int...)
    @boundscheck checkbounds(iter, i, I...)
    @inbounds iter[i]
end

@propagate_inbounds function Base.getindex(iter::LinearIndicesHalfInt, i::AbstractRange{<:HalfIntegerOrInteger})
    @boundscheck checkbounds(iter, i)
    (first(iter):last(iter))[i]
end
@propagate_inbounds function Base.getindex(iter::LinearIndicesHalfInt{1}, i::AbstractRange{<:HalfIntegerOrInteger})
    @boundscheck checkbounds(iter, i)
    j = unsafeUnitRangeInt(i .- iter.offsets[1])
    IdOffsetRange(iter.lininds.indices[1][j], iter.offsets[1])
end
@propagate_inbounds function Base.getindex(iter::LinearIndicesHalfIntParent{1}, i::AbstractRange{<:HalfIntegerOrInteger})
    @boundscheck checkbounds(iter, i)
    iter.lininds.indices[1][i] .+ iter.offsets[1]
end

Base.iterate(iter::LinearIndicesHalfInt{1}, s...) = iterate(Base.axes1(iter), s...)
Base.iterate(iter::LinearIndicesHalfInt, i=1) = i > length(iter) ? nothing : (i, i+1)

# Needed since firstindex and lastindex are defined in terms of LinearIndicesHalfInt
Base.first(iter::LinearIndicesHalfInt) = 1
Base.first(iter::LinearIndicesHalfInt{1}) = first(iter.lininds.indices[1]) + iter.offsets[1]
Base.last(iter::LinearIndicesHalfInt) = length(iter)
Base.last(iter::LinearIndicesHalfInt{1}) = last(iter.lininds.indices[1]) + iter.offsets[1]

function Base.hash(c::LinearIndicesHalfInt, h::UInt)
    h = hash(c.offsets, h)
    h = hash(c.lininds, h)
    h = hash(:LinearIndicesHalfInt, h)
    return h
end

# It might be possible to construct a LinearIndices range for integer indices
function Base.LinearIndices(A::AbstractHalfIntegerArrayOrWrapper)
    LinearIndices(IdentityUnitRange.(UnitRange{Int}.(axes(A))))
end


text:
using Test
using KrylovKit, SparseArrays, ArnoldiMethod, LinearMaps, LinearAlgebra
using PseudoArcLengthContinuation

N = 110
A = I + sprand(N, N, 0.05)

T, vecs, vals, info = @time KrylovKit.schursolve(A, rand(N), 10, :LM, Arnoldi(krylovdim = 30, verbosity = 0))
vals = eigvals(Array(A))

# test some definitions
P = PrecPartialSchurKrylovKit(A, rand(N), 4, :LM)
P = PrecPartialSchurArnoldiMethod(A, 4, LM())
P = PrecPartialSchurArnoldiMethod(A, N, 4, LM())

# test some function
ldiv!(P, rand(N))
ldiv!(rand(N), P, rand(N))

# test that it deflates the eigenvalues
Jmap = LinearMap(x-> A * (P \ x), N, N ; ismutating = false)
decomp, history = ArnoldiMethod.partialschur(Jmap, nev = 110, tol=1e-9, which=LM());
@test Base.setdiff(round.(vals, digits = 5), round.(decomp.eigenvalues, digits = 5)) |> length in [4,5]

#
# using Plots
# scatter(real.(vals), imag.(vals), label = "")
# 	scatter!(real.(decomp.eigenvalues), imag.(decomp.eigenvalues), label = "", marker=:cross)

ext: jl


text:
using Test
using ADBCVUMPS
using ADBCVUMPS: energy, num_grad, diaglocal
using OMEinsum, Zygote, Random
using Optim, LineSearches
using LinearAlgebra: svd, norm

@testset "non-interacting" for Ni = [2], Nj = [2]
    Random.seed!(100)
    rd = randn()
    rd2 = randn(2,2,2,2,2)
    model = diaglocal(Ni,Nj,[1,-1.0])
    h = hamiltonian(model)
    as = (reshape([rand(3,3,3,3,2) for i = 1:Ni*Nj],(Ni,Nj)) for _ in 1:10)
    @test all(a -> -1 < energy(h, model, SquareBCIPEPS(a); χ=5, tol=1e-10, maxiter=10)/2 < 1, as)

    a = reshape([zeros(2,2,2,2,2) .+ 1e-12 * rd2 for i = 1:Ni*Nj],(Ni,Nj))
    for j=1:Nj, i=1:Ni
        a[i,j][1,1,1,1,2] = rd
    end
    @test energy(h, model, SquareBCIPEPS(a); χ=4, tol=1e-10, maxiter=10)/2 ≈ -1

    a = reshape([zeros(2,2,2,2,2) .+ 1e-12 * rd2 for i = 1:Ni*Nj],(Ni,Nj))
    for j=1:Nj, i=1:Ni
        a[i,j][1,1,1,1,1] = rd
    end
    @test energy(h, model, SquareBCIPEPS(a); χ=4, tol=1e-10, maxiter=10)/2 ≈ 1

    a = reshape([zeros(2,2,2,2,2) .+ 1e-12 * rd2 for i = 1:Ni*Nj],(Ni,Nj))
    for j=1:Nj, i=1:Ni
        a[i,j][1,1,1,1,1] = a[i,j][1,1,1,1,2] = rd
    end
    @test abs(energy(h,model,SquareBCIPEPS(a); χ=4, tol=1e-10, maxiter=10)) < 1e-9

    grad = let energy = x -> real(energy(h, model, SquareBCIPEPS(reshape([x for i = 1:Ni*Nj],(Ni, Nj))); χ=4, tol=1e-10, maxiter=10))
        res = optimize(energy,
            Δ -> Zygote.gradient(energy,Δ)[1], a[1,1], LBFGS(m=20), inplace = false)
    end
    @test grad != Nothing

    hdiag = [0.3,-0.43]
    model = diaglocal(Ni,Nj,hdiag)
    bcipeps, key = init_ipeps(model; D=2, χ=4, tol=1e-10, maxiter=20)
    res = optimiseipeps(bcipeps, key; f_tol = 1e-6)
    e = minimum(res)/2
    @test isapprox(e, minimum(hdiag), atol=1e-3)
end

@testset "gradient" for Ni = [2], Nj = [2]
    Random.seed!(0)
    model = TFIsing(Ni,Nj,1.0)
    h = hamiltonian(model)
    bcipeps, key = init_ipeps(model; D=2, χ=4, tol=1e-10, maxiter=20)
    gradzygote = first(Zygote.gradient(bcipeps) do x
        energy(h,model,x; χ=4, tol=1e-10, maxiter=20)
    end).bulk
    gradnum = num_grad(bcipeps.bulk, δ=1e-3) do x
        energy(h,model,SquareBCIPEPS(x); χ=4, tol=1e-10, maxiter=20)
    end
    @test isapprox(gradzygote, gradnum, atol=1e-3)

    # Random.seed!(3)
    # model = Heisenberg(Ni,Nj)
    # h = hamiltonian(model)
    # bcipeps, key = init_ipeps(model; D=2, χ=4, tol=1e-10, maxiter=20)
    # gradzygote = first(Zygote.gradient(bcipeps) do x
    #     energy(h,model,x; χ=4, tol=1e-10, maxiter=20)
    # end).bulk
    # gradnum = num_grad(bcipeps.bulk, δ=1e-3) do x
    #     energy(h,model, SquareBCIPEPS(x); χ=4, tol=1e-10, maxiter=20)
    # end
    # @test isapprox(gradzygote , gradnum, atol=1e-3)
end

@testset "TFIsing" for Ni = [2], Nj = [2]
    # comparison with results from https://github.com/wangleiphy/tensorgrad
    Random.seed!(3)
    model = TFIsing(Ni,Nj,1.0)
    bcipeps, key = init_ipeps(model; D=2, χ=5, tol=1e-10, maxiter=20)
    res = optimiseipeps(bcipeps, key; f_tol = 1e-6)
    e = minimum(res)
    @test isapprox(e, -2.12566, atol = 1e-2)

    # Random.seed!(3)
    # model = TFIsing(Ni,Nj,0.5)
    # bcipeps, key = init_ipeps(model; D=2, χ=5, tol=1e-10, maxiter=20)
    # res = optimiseipeps(bcipeps, key; f_tol = 1e-6)
    # e = minimum(res)
    # @test isapprox(e, -2.0312, atol = 1e-2)

    # Random.seed!(3)
    # model = TFIsing(Ni,Nj,2.0)
    # bcipeps, key = init_ipeps(model; D=2, χ=5, tol=1e-10, maxiter=20)
    # res = optimiseipeps(bcipeps, key; f_tol = 1e-6)
    # e = minimum(res)
    # @test isapprox(e, -2.5113, atol = 1e-2)
end

@testset "heisenberg" for Ni = [2], Nj = [2]
    # comparison with results from https://github.com/wangleiphy/tensorgrad
    Random.seed!(100)
    model = Heisenberg(Ni,Nj,1.0,1.0,1.0)
    bcipeps, key = init_ipeps(model; D=2, χ=5, tol=1e-10, maxiter=20)
    res = optimiseipeps(bcipeps, key; f_tol = 1e-6)
    e = minimum(res)
    @test isapprox(e, -0.66023, atol = 1e-4)

    # Random.seed!(100)
    # model = Heisenberg(Ni,Nj,1.0,2.0,2.0)
    # bcipeps, key = init_ipeps(model; D=2, χ=5, tol=1e-10, maxiter=20)
    # res = optimiseipeps(bcipeps, key; f_tol = 1e-6)
    # e = minimum(res)
    # @test isapprox(e, -1.190, atol = 1e-3)

    # Random.seed!(100)
    # model = Heisenberg(Ni,Nj,2.0,0.5,0.5)
    # bcipeps, key = init_ipeps(model; D=2, χ=5, tol=1e-10, maxiter=20)
    # res = optimiseipeps(bcipeps, key; f_tol = 1e-6)
    # e = minimum(res)
    # @test isapprox(e, -1.0208, atol = 1e-3)
end


text:
#!NOTE: Make Distribution.Distributions work in ModelWrappers and map a Bijector to each distribution.
############################################################################################
function construct_transform(d::Distributions.Distribution, val)
    transform = Bijectors.bijector(d)
    return transform, Bijectors.inverse(transform)
end

############################################################################################
function _check(
    _rng::Random.AbstractRNG,
    d::Distributions.Distribution,
    val::Union{R,Array{R},AbstractArray},
) where {R<:Real}
    _val = rand(_rng, d)
    return typeof(val) == typeof(_val) && size(val) == size(_val) ? true : false
end

############################################################################################
# Implicit constraint through density by specialization -> else just flatten
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::U,
    constraint::Distributions.Distribution,
    x::Union{R,Array{R}},
) where {
    T<:AbstractFloat,
    F<:FlattenTypes,
    U<:UnflattenTypes,
    R<:Real,
}
    return construct_flatten(T, flattentype, unflattentype, x)
end

############################################################################################
#=
!NOTE: For this constraint, we use a Bijector as Transformer, so we do not need to add a new functors
1. MyTransformer <: AbstractTransformer, MyInverseTransformer <: AbstractTransformer
2. define a function construct_transform(MyConstraint, val) -> MyTransformer, MyInverseTransformer
3. overload unconstrain and log_abs_det_jac on MyTransformer.
4. overload constrain on MyInverseTransformer.
=#

############################################################################################
#Export
export
    construct_flatten,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
type Collapsible
  header
  content
end

Collapsible(a::AbstractString, b::AbstractString = "") = Collapsible(HTML(a), HTML(b))
Collapsible(a) = Collapsible(a, HTML(""))

function register_collapsible(c)
  _currentresult_ == nothing && return
  id = uuid4()
  _currentresult_.data[id] = c
  return id
end

collapsibleclickjs(id) = """
    var content = this.parentNode.querySelector('.collapsible-content');
    if (content.classList.contains('lazy')) {
      lt.objs.notifos.working();
      $(jlcall(""" LightTable.collapsibleclick("$id") """))
    } else {
      \$(content).toggle(200);
      setTimeout(function () {
        lt.objs.editor.refresh(lt.objs.editor.pool.last_active());
      }, 210);
    }
  """

collapsibleclick(id) = jscall("""
  node = document.getElementById('$(id)');
  node.innerHTML = '$(jsescapestring(stringmime("text/html", _currentresult_.data[UUID(id)].content)))';
  lt.objs.langs.julia.process(node);
  node.classList.remove('lazy');
  lt.objs.notifos.done_working();
  \$(node).show(200)
  """)

function Base.writemime(io::IO, m::MIME"text/html", c::Collapsible)
  id = register_collapsible(c)

  println(io, """
""")

  println(io, """""")
  writemime(io, MIME"text/html"(), c.header)
  println(io, """""")
  if _currentresult_ != nothing
    println(io, """

""")
  else
    println(io, """
""")
    writemime(io, m, c.content)
    println(io, """
""")
  end

  println(io, """
""")
end

displayinline!(content::Collapsible, opts) =
  showresult(stringmime("text/html", content), opts, html=true)


text:
module escape_num_dimer
export escape_exit_num_dimer
using DifferentialEquations
using Roots
using Printf

# using MATLAB
# using ProgressMeter

include("leap_frog_definitions.jl")
max_hit_q1 = 10000
max_hit_p2 = 10000
barrier = 5

condition_max_hits(u, t, integrator) = u[3] > max_hit_q1 || u[6]>max_hit_p2 || maximum([abs(u[1]),abs(u[2]),abs(u[4]),abs(u[5])]) > barrier
affect_stop!(integrator) = terminate!(integrator)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end

function condition_hits_PSS_p2(u, t, integrator) # Event when event_f(u,t) == 0
    u[5]
end

function affect_update_iterator_q1!(integrator)
    integrator.u[3] = integrator.u[3] + 1
end
function affect_update_iterator_p2!(integrator)
    integrator.u[6] = integrator.u[6] + 1
end

callback_max_hits = DiscreteCallback(condition_max_hits, affect_stop!)
callback_hits_PSS_q1 = ContinuousCallback(condition_hits_PSS_q1, affect_update_iterator_q1!, nothing)
callback_hits_PSS_p2 = ContinuousCallback(condition_hits_PSS_p2, affect_update_iterator_p2!, nothing)

cb = CallbackSet(callback_hits_PSS_q1, callback_hits_PSS_p2, callback_max_hits)

function escape_exit_num_dimer(mesh_list, t_end, H)
    Q2 = mesh_list[1]; P2 = mesh_list[2];
    P1 = P1_find_dimer(Q2, P2, H)

    if isempty(P1)
        P1 = P1_find_dimer_second(Q2, P2, H)
    end

    if ~isempty(P1)
        # println(Q1)
        q0, p0 = [zeros(3) for i in 1:2]
        q0[1] = 0
        q0[2] = Q2
        q0[3] = 0
        p0[1] = P1[1]
        p0[2] = P2
        p0[3] = 0
       #constructor for ODE
        prob = HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., t_end));
        #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        sol = solve(prob, Vern9(),  reltol=1e-13,abstol=1e-15,maxiters = 1e10,  callback = cb, save_start = true, save_end = true, save_everystep = false)
        # uf=zeros(6)
        # uf[1]=sol[1,end]
        # uf[2]=sol[2,end]
        # uf[3]=sol[3,end]
        # uf[4]=sol[4,end]
        # uf[5]=sol[5,end]
        # uf[6]=sol[6,end]

        if sol.u[end][6] > max_hit_p2 - 1 ||  sol.u[end][3] > max_hit_q1 - 1
            return t_end
        else
            return sol.t[end]

        end
        # if dH<1e-5
        #     if is_it(uf,H,tol_dist)
        #         return sol.u[end][5]
        #     else
        #         return -1
        #     end
        # end
    else
        return NaN
    end
end

end


text:
## ExoplanetsSysSim/examples/setup_abc_christiansen.jl
## (c) 2016 Eric B. Ford & Danley C. Hsu

module EvalSysSimModel
  export setup, get_param_vector, get_ss_obs #, evaluate_model
  export gen_data, calc_summary_stats, calc_distance, is_valid
  using ExoplanetsSysSim
  include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","multiplanet_systems", "model.jl"))
  include(joinpath(pwd(), "param_file.jl"))

  sim_param_closure = SimParam()
  summary_stat_ref_closure =  CatalogSummaryStatistics()
  dist_weights = ones(Float64,0)

    function is_valid(param_vector::Vector{Float64})
      global sim_param_closure
      update_sim_param_from_vector!(param_vector,sim_param_closure)
      if haskey(sim_param_closure,"mean_log_planet_radius")
        const mu_log_r::Float64 = get_real(sim_param_closure,"mean_log_planet_radius")
        const sigma_log_r::Float64 = get_real(sim_param_closure,"sigma_log_planet_radius")
        const min_radius::Float64 = get_real(sim_param_closure,"min_radius")
        const max_radius::Float64 = get_real(sim_param_closure,"max_radius")
        if sigma_log_r <= 0. || !(log(min_radius)<=mu_log_r<=log(max_radius))
         return false
        end
      end
      if haskey(sim_param_closure,"mean_log_planet_period")
        const mu_log_P::Float64 = get_real(sim_param_closure,"mean_log_planet_period")
        const sigma_log_P::Float64 = get_real(sim_param_closure,"sigma_log_planet_period")
        const min_period::Float64 = get_real(sim_param_closure,"min_period")
        const max_period::Float64 = get_real(sim_param_closure,"max_period")
        if sigma_log_P<=0. || !(log(min_period)<=mu_log_P<=log(max_period))
         return false
        end
      end
      if haskey(sim_param_closure,"frac_zero_planet")
         const frac_zero_planet::Float64 = get_real(sim_param_closure,"frac_zero_planet")
         if !(0.0<=frac_zero_planet<=1.0)
            return false
         end
      end
      if haskey(sim_param_closure,"frac_one_planet")
         const frac_one_planet::Float64 = get_real(sim_param_closure,"frac_one_planet")
         if !(0.0<=frac_one_planet<=1.0)
            return false
         end
      end

      if haskey(sim_param_closure,"fracs_num_planets")
         fracs::Array{Float64,1} = get(sim_param_closure,"fracs_num_planets",ones(1))
         if any(.!(0.0.<=fracs.<=1))
            return false
         end
      end

      if haskey(sim_param_closure,"eta_pl")
         eta_pl = get_real(sim_param_closure,"eta_pl")
         if any(!(0.0.<=eta_pl))
            return false
         end
      end

#      const rate_tab::Array{Float64,2} = get_any(sim_param_closure, "obs_par", Array{Float64,2})
#      const lambda = sum_kbn(rate_tab)
    # println("Lambda = ", lambda)
#      if sigma_log_r <= 0. || sigma_log_P<=0. || !(log(min_radius)<=mu_log_r<=log(min_radius)) || !(log(min_period)<=mu_log_P<=log(min_period)) || !(0.0<=frac_zero_planet<=1.0) || !(0.0<=frac_one_planet<=1.0) # || lambda > 10. || any(x -> x < 0., rate_tab)
#         return false
#      end
      return true
    end

#   Force the fracts_num_planets vector to be a non-negative simplex.
#   If negative values happen often, then should change or adjust IS density to compensate
#   TODO FUTURE: I think Jessi has a better way to draw such vectors, but implementing that that might complicate the proposal step
    function normalize_categorical(param_vector::Vector{Float64})
      global sim_param_closure
      r::Range = get_range_for_sim_param("fracs_num_planets",sim_param_closure)
      param_vector[r] = abs(param_vector[r])
      param_vector[r] = param_vector[r]/sum(param_vector[r])
      update_sim_param_from_vector!(param_vector,sim_param_closure)

      if is_active(sim_param_closure,"sigma_incl")
        sig = get_real(sim_param_closure,"sigma_incl")
        if sig<0
           update_param(sim_param_closure,"sigma_incl",abs(sig))
        end
      end
      param_vector[1:end] = get_param_vector()
      return param_vector
    end

    function normalize_poisson_mixture(param_vector::Vector{Float64})
      global sim_param_closure
      update_sim_param_from_vector!(param_vector,sim_param_closure)
      frac_zero::Float64 = min(max(get_real(sim_param_closure,"frac_zero_planet"), 0.0),1.0)
      frac_one::Float64  = min(max(get_real(sim_param_closure,"frac_one_planet"), 0.0),1.0)
      frac_multi::Float64  = min(max(get_real(sim_param_closure,"frac_multi_planet"),0.0), 1.0)
      frac_sum = frac_zero+frac_one+frac_multi
      frac_zero /= frac_sum
      frac_one /= frac_sum
      frac_multi /= frac_sum
      update_param(sim_param_closure, "frac_zero_planet",frac_zero)
      update_param(sim_param_closure, "frac_one_planet",frac_one)
      update_param(sim_param_closure, "frac_multi_planet",frac_multi)

      if is_active(sim_param_closure,"sigma_incl")
        sig = get_real(sim_param_closure,"sigma_incl")
        if sig<0
           update_param(sim_param_closure,"sigma_incl",abs(sig))
        end
      end
      if is_active(sim_param_closure,"eta_pl")
        eta = get_real(sim_param_closure,"eta_pl")
        if eta<0
           update_param(sim_param_closure,"eta_pl",abs(eta))
        end
      end
      if is_active(sim_param_closure,"sigma_hk")
        sig = get_real(sim_param_closure,"sigma_hk")
        if sig<0
           update_param(sim_param_closure,"sigma_hk",abs(sig))
        end
      end
      param_vector[1:end] = get_param_vector()
      return param_vector
    end

#=
    function norm_log_simplex(param_vector::Vector{Float64})
      log_z_sum = logsumexp(abs(param_vector))
      param_vector = exp(abs(param_vector)-z_sum)
    end
=#

    function gen_data(param_vector::Vector{Float64})
      global sim_param_closure
      update_sim_param_from_vector!(param_vector,sim_param_closure)
      cat_phys = generate_kepler_physical_catalog(sim_param_closure)
      cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys, sim_param_closure)
      cat_obs = ExoplanetsSysSim.observe_kepler_targets_single_obs(cat_phys_cut, sim_param_closure)
      return cat_obs
    end

    function calc_summary_stats(cat::KeplerObsCatalog )
      global sim_param_closure
      sum_stat = calc_summary_stats_model(cat, sim_param_closure)
      return sum_stat
    end

    function calc_distance_no_weights(sum_stat_obs::CatalogSummaryStatistics,sum_stat_sim::CatalogSummaryStatistics, n::Integer = 0)
      global sim_param_closure
      dist = calc_distance_model(sum_stat_obs,sum_stat_sim, sim_param_closure)
      num_available = length(dist)
      num_to_use = n>0 ? min(n,num_available) : num_available
      dist_scalar = calc_scalar_distance(dist[1:num_to_use])
      println("param = ", get_param_vector()," dist = ",dist_scalar," : ", dist )
      return dist_scalar
    end
    function calc_distance(sum_stat_obs::CatalogSummaryStatistics,sum_stat_sim::CatalogSummaryStatistics, n::Integer = 0)
      global sim_param_closure, dist_weights
      dist = calc_distance_model(sum_stat_obs,sum_stat_sim, sim_param_closure)
      num_available = length(dist)
      num_to_use = n>0 ? min(n,num_available) : num_available
      @assert length(dist_weights)>=num_to_use
      idx_to_use = dist_weights[1:num_to_use] .> 0.0
      dist_scalar = calc_scalar_distance(dist[idx_to_use].*dist_weights[idx_to_use])
      println("param = ", get_param_vector()," dist = ",dist_scalar," : ", dist, "; weighted = ", dist[idx_to_use].*dist_weights[idx_to_use])
      return dist_scalar
    end

  function calc_distance_scales(param::Vector{Float64}; n::Integer = 20)
     global sim_param_closure
     println("# param =",param)
     cat_obs = gen_data(param)
     sum_stat_obs = calc_summary_stats(cat_obs)
     dist_test = calc_distance_model(sum_stat_obs,sum_stat_obs,sim_param_closure)
     dist_arr = Array{Float64}(length(dist_test),n)
     for i in 1:n
        cat_sim = gen_data(param)
        sum_stat_sim = calc_summary_stats(cat_sim)
        dist_arr[:,i] = calc_distance_model(sum_stat_obs,sum_stat_sim,sim_param_closure)
     end
     dist_medians = median(dist_arr,2)
     println("# Distance Means  =", mean(dist_arr,2))
     println("# Distance Medians=", median(dist_arr,2))
     println("# Distance Stddevs=", std(dist_arr,2))
     println("# 25th Percentile= ", map(i->quantile(dist_arr[i,:],0.25),1:size(dist_arr,1)))
     println("# 75th Percentile= ", map(i->quantile(dist_arr[i,:],0.75),1:size(dist_arr,1)))
     global dist_weights = zeros(Float64,length(dist_medians))
     for i in 1:length(dist_medians)
         if dist_medians[i]>0.0
            dist_weights[i] = 1.0/dist_medians[i]
         end
     end
     println("# Distance Weights=", dist_weights)
  end

  function setup()
    global sim_param_closure = setup_sim_param_model()

    sim_param_closure = set_test_param(sim_param_closure)
    add_param_fixed(sim_param_closure,"stellar_catalog","q1_q17_dr25_stellar.jld") # "q1_q12_christiansen.jld")

    ### Use simulated planet candidate catalog data
    add_param_fixed(sim_param_closure,"num_kepler_targets",150969)  # For "observed" catalog
    #cat_obs = simulated_read_kepler_observations(sim_param_closure)
    #inv_det_prob(cat_obs, sim_param_closure)
    #println("------------------------")
    #global summary_stat_ref_closure = calc_summary_stats_model(cat_obs,sim_param_closure, trueobs_cat = true)
    ###

    ### Use real planet candidate catalog data
    add_param_fixed(sim_param_closure,"koi_catalog","q1_q16_koi_cand.csv")
    df_star = setup_star_table_christiansen(sim_param_closure)
    df_koi,usable_koi = read_koi_catalog(sim_param_closure)
    cat_obs = setup_actual_planet_candidate_catalog(df_star, df_koi, usable_koi, sim_param_closure)
    global summary_stat_ref_closure = calc_summary_stats_model(cat_obs,sim_param_closure, trueobs_cat = true)
    ###

    num_targ_sim = floor(Int64,get_int(sim_param_closure,"num_kepler_targets")/10)  # TODO WARNING TEMPORARY: This is just for testing to speed things along...
    add_param_fixed(sim_param_closure,"num_targets_sim_pass_one",num_targ_sim)  # Set universal simulated catalog size

    calc_distance_scales(get_param_vector())
  end

  get_param_vector() = make_vector_of_sim_param(sim_param_closure)
  get_ss_obs() = summary_stat_ref_closure

  function set_simparam_ss(sim_param::ExoplanetsSysSim.SimParam, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics)
    global sim_param_closure = sim_param
    global summary_stat_ref_closure = ss_true
  end

  function profile_model_eval(n::Integer)
    global sim_param_closure
    setup()
    global cat_phys_try_closure = generate_kepler_physical_catalog(sim_param_closure)
    cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys, sim_param_closure)
    global cat_obs_try_closure = observe_kepler_targets_single_obs(cat_phys_try_closure,sim_param_closure)
    global summary_stat_ref_closure = calc_summary_stats(cat_obs_try_closure)

    param_base = make_vector_of_sim_param(sim_param_closure)

    evaluate_model_demo_scalar_param_scalar_ret(param::Float64) = sum(evaluate_model([param;param_base[2:end]]))

    io = open("profile_syssim.tmp","w")
    Profile.print(io) # ,cols=800)
    close(io)
    #rm(io)

    Profile.init(10000000,0.001)
    eta_min = 0.5*exp(get_real(sim_param_closure,"log_eta_pl"))
    eta_max = 2.0*exp(get_real(sim_param_closure,"log_eta_pl"))
    #gc(); if VERSION >= VersionNumber(0,4,0)     gc_enable(false);   else     gc_disable();   end
    tic()
    for (i,eta) in enumerate(linspace(eta_min,eta_max,n) )
      local dist
      if i==2
         Profile.clear()
         Profile.clear_malloc_data()
      end
      @profile dist = evaluate_model_demo_scalar_param_scalar_ret(log(eta))
      println("eta= ",eta," d= ", dist)
    end
    toc()
    #gc_enable()
    pwd()
    io = open("profile_syssim.out","w")
    Profile.print(io,cols=800)
    close(io)
  end

end  # module EvalSysSimModel

include(joinpath(Pkg.dir("ABC"),"src/composite.jl"))

module SysSimABC
  export setup_abc, run_abc
  import ABC
  import Distributions
  using CompositeDistributions
  import ExoplanetsSysSim
  import EvalSysSimModel
  include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","multiplanet_systems", "model.jl"))
  #include(joinpath(pwd(), "param_file.jl"))

  function setup_abc(num_dist::Integer = 0)
    EvalSysSimModel.setup()
    #println("# setup_abc: ",EvalSysSimModel.sim_param_closure)
    theta_true = EvalSysSimModel.get_param_vector()
    prior_dict = Dict("log_eta_pl" => Distributions.Uniform(0.1, log(10.)),"eta_pl" => Distributions.Uniform(0., 10.), "frac_zero_planet" => Distributions.Uniform(0.0, 1.0), "frac_one_planet" => Distributions.Uniform(0.0, 1.0), "frac_multi_planet" => Distributions.Uniform(0.0,1.0),
       "fracs_num_planets" => Distributions.Dirichlet(get_int(EvalSysSimModel.sim_param_closure,"max_planets_in_sys")+1,1.0),
       "mean_log_planet_radius" => Distributions.Uniform(log(0.5*ExoplanetsSysSim.earth_radius), log(4.*ExoplanetsSysSim.earth_radius)), "mean_log_planet_period" => Distributions.Uniform(log(1.), log(100.)),
       "sigma_log_planet_radius" => Distributions.Uniform(0., log(4.)), "sigma_log_planet_period" => Distributions.Uniform(0., log(4.)),
       "power_law_P" => Distributions.Uniform(-1.0, 1.0), "power_law_r" => Distributions.Uniform(-4.0, -1.0),
       "sigma_hk" => Distributions.Uniform(0.0,1.0),  "sigma_hk_one" => Distributions.Uniform(0.0,1.0),  "sigma_hk_multi" => Distributions.Uniform(0.0,1.0),
       "sigma_incl" => Distributions.Uniform(0.0,30.0)   # degrees; 0 = coplanar w/ generate_kepler_target_simple; ignored by generate_planetary_system_uncorrelated_incl
)
    active_param_keys_sorted = make_vector_of_active_param_keys(EvalSysSimModel.sim_param_closure)
    array_of_priors = Array{Distributions.ContinuousDistribution}(length(active_param_keys_sorted) )
    for i in 1:length(array_of_priors)
        array_of_priors[i]  = prior_dict[active_param_keys_sorted[i]];
    end
    param_prior = CompositeDist( array_of_priors )
    # param_prior = CompositeDist( Distributions.ContinuousDistribution[Distributions.Uniform(0., 2.) for x in 1:length(theta_true)] )
    in_parallel = nworkers() > 1 ? true : false

    calc_distance_ltd(sum_stat_obs::ExoplanetsSysSim.CatalogSummaryStatistics,sum_stat_sim::ExoplanetsSysSim.CatalogSummaryStatistics) = EvalSysSimModel.calc_distance(sum_stat_obs,sum_stat_sim,num_dist)

    global abc_plan = ABC.abc_pmc_plan_type(EvalSysSimModel.gen_data,EvalSysSimModel.calc_summary_stats,calc_distance_ltd, param_prior, make_proposal_dist=ABC.make_proposal_dist_gaussian_diag_covar, is_valid=EvalSysSimModel.is_valid,
         #normalize=EvalSysSimModel.normalize_poisson_mixture,
         normalize=EvalSysSimModel.normalize_categorical,
                                     num_part=20, num_max_attempt=200, num_max_times=60, epsilon_init=9.9e99, target_epsilon=1.0e-100, in_parallel=in_parallel, adaptive_quantiles = true, epsilon_reduction_factor=0.8, tau_factor=2.0);
  end

  function run_abc_largegen(pop::ABC.abc_population_type, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics, epshist_targ::Float64, npart::Integer = 1000, num_dist::Integer = 0)
    sim_param_closure = setup_sim_param_model()
    sim_param_closure = EvalSysSimModel.set_test_param(sim_param_closure)
    EvalSysSimModel.set_simparam_ss(sim_param_closure, ss_true)

    theta_true = EvalSysSimModel.get_param_vector()
    prior_dict = Dict("log_eta_pl" => Distributions.Uniform(0.1, log(10.)), "eta_pl" => Distributions.Uniform(0., 10.),"frac_zero_planet" => Distributions.Uniform(0.0, 1.0), "frac_one_planet" => Distributions.Uniform(0.0, 1.0), "frac_multi_planet" => Distributions.Uniform(0.0,1.0),
       "fracs_num_planets" => Distributions.Dirichlet(10,1.0),
       "mean_log_planet_radius" => Distributions.Uniform(log(0.5*ExoplanetsSysSim.earth_radius), log(4.*ExoplanetsSysSim.earth_radius)), "mean_log_planet_period" => Distributions.Uniform(log(1.), log(100.)),
       "sigma_log_planet_radius" => Distributions.Uniform(0., log(4.)), "sigma_log_planet_period" => Distributions.Uniform(0., log(4.)),
       "power_law_P" => Distributions.Uniform(-1.0, 1.0), "power_law_r" => Distributions.Uniform(-4.0, -1.0),
       "sigma_hk" => Distributions.Uniform(0.0,1.0),  "sigma_hk_one" => Distributions.Uniform(0.0,1.0),  "sigma_hk_multi" => Distributions.Uniform(0.0,1.0),
)
    active_param_keys_sorted = make_vector_of_active_param_keys(EvalSysSimModel.sim_param_closure)
    array_of_priors = Array{Distributions.ContinuousDistribution}(length(active_param_keys_sorted) )
    for i in 1:length(array_of_priors)
        array_of_priors[i]  = prior_dict[active_param_keys_sorted[i]];
    end
    param_prior = CompositeDist( array_of_priors )
    # param_prior = CompositeDist( Distributions.ContinuousDistribution[Distributions.Uniform(0., 2.) for x in 1:length(theta_true)] )
    in_parallel = nworkers() > 1 ? true : false

    calc_distance_ltd(sum_stat_obs::ExoplanetsSysSim.CatalogSummaryStatistics,sum_stat_sim::ExoplanetsSysSim.CatalogSummaryStatistics) = EvalSysSimModel.calc_distance(sum_stat_obs,sum_stat_sim,num_dist)

    global abc_plan = ABC.abc_pmc_plan_type(EvalSysSimModel.gen_data,EvalSysSimModel.calc_summary_stats,calc_distance_ltd, param_prior, is_valid=EvalSysSimModel.is_valid,
                                     num_part=npart, num_max_attempt=200, num_max_times=1, epsilon_init=9.9e99, target_epsilon=1.0e-100, in_parallel=in_parallel);

    println("# run_abc_largegen: ",EvalSysSimModel.sim_param_closure)
    sampler_plot = abc_plan.make_proposal_dist(pop, abc_plan.tau_factor)
    theta_plot = Array{Float64}(size(pop.theta, 1), npart)
    for i in 1:npart
      theta_plot[:,i], dist_plot, attempts_plot = ABC.generate_theta(abc_plan, sampler_plot, ss_true, epshist_targ)
    end
    return theta_plot
  end

  function run_abc(abc_plan::ABC.abc_pmc_plan_type)
    println("# run_abc: ",EvalSysSimModel.sim_param_closure)
    ss_true = EvalSysSimModel.get_ss_obs()
    #println("True catalog SS: ", ss_true)
    pop_out = ABC.run_abc(abc_plan,ss_true;verbose=true)
  end

  function run_abc(abc_plan::ABC.abc_pmc_plan_type, pop::ABC.abc_population_type)
    dist_threshold = maximum(pop.dist)
    EvalSysSimModel.add_param_fixed(EvalSysSimModel.sim_param_closure,"minimum ABC dist skip pass 2",dist_threshold)
    println("# run_abc: ",EvalSysSimModel.sim_param_closure)
    ss_true = EvalSysSimModel.get_ss_obs()
    pop_out = ABC.run_abc(abc_plan,ss_true,pop;verbose=true)
  end

end # module SysSimABC


text:
"""
    inject!(x::AbstractInteractor, y::AbstractPopulation)

Add an interactor to a population.
"""
function inject!(x::AbstractInteractor, y::AbstractPopulation)
  push!(y.census, x)
end

"""
    eject!(x::Int, y::AbstractPopulation)

Remove interactor with ID `x` from population `y`.
"""
function eject!(x::Int, y::AbstractPopulation)
  splice!(y.census, x)
end

"""
    eject!(x::AbstractPopulation)

Remove a randomly chosen interactor from a population.
"""
function eject!(x::AbstractPopulation)
  splice!(x.census, rand(1:length(x.census)))
end

"""
    rendezvous!(x::AbstractPopulation, y::Int, z::Int)

Conduct a rendezvous between interactors indexed by `y` and `z` in population `x`.
(The first mentioned acts first.)
"""
function rendezvous!(x::AbstractPopulation, y::Int, z::Int)
  act!(x.census[y], x.census[z])
  act!(x.census[z], x.census[y])
  x.nactions += 2
end

"""
    rendezvous!(x::AbstractPopulation)

Conduct a rendezvous between two interactors chosen at random in population `x`.
"""
function rendezvous!(x::AbstractPopulation)
  length(x.census) > 1 || return false
  int1 = 0
  int2 = 0
  while true
    int1 = rand(1:length(x.census))
    int2 = rand(1:length(x.census))
    int1 != int2 && break
  end
  rendezvous!(x, int1, int2)
end

"""
    characterize(x::AbstractPopulation)

Characterize the state of a population. Tries to form either a single number
or a vector that reflects the population's average behaviour, depending on
the types of replicators carried by the interactions in the population. In
practice, this is achieved by calling on `characterize` functions on the
interactors.
"""
function characterize(x::AbstractPopulation)
  length(x.census) > 0 || return false
  result = characterize(x.census[1])
  for person in x.census[2:end]
    result += characterize(person)
  end
  result/length(x.census)
end

"""
    characterize_by_location(x::AbstractLatticePopulation)

Characterize a lattice population cell by cell.
"""
function characterize_by_location(x::AbstractLatticePopulation)
  out = zeros(x.size, x.size)
  speakers_by_cell = zeros(x.size, x.size)

  for s in 1:length(x.census)
    speaker = x.census[s]
    sx = x.locations[s].x
    sy = x.locations[s].y
    out[sx, sy] += characterize(speaker)
    speakers_by_cell[sx, sy] += 1
  end

  out ./ speakers_by_cell
end


text:
function findPi(num_samples, trials_per_sample)
  samples = Array{Float64}(undef, num_samples)
  mean = 0.0

  for sample = 1:num_samples
    samples[sample] = integrateCircle(trials_per_sample)
  end

  mean = sum(samples) / size(samples)[1]

  std_dev = getStdDev(mean, samples)

  return mean, std_dev
end

function getStdDev(mean, samples)
  variance = 0

  for sample in samples
    variance += (sample - mean) ^ 2
  end

  variance /= size(samples)[1]
  return sqrt(variance)
end

function integrateCircle(points_placed)
  counter = 0

  for trial = 1:points_placed
    y_pos = rand()
    x_pos = rand()

    dist = y_pos ^ 2 + x_pos ^ 2
    if dist <= 1
      counter += 1
    end
  end

  return 4.0 * counter / points_placed
end

function monte_carlo()
  res = [0.0, 0.0]
  for i in 1:20
    res1, res2 = findPi(1000, 50000)
    res[1] += res1 / 20.0
    res[2] += res2 / 20.0
  end
  print(res)
end

max_stars_repo_name: pnadon/monte-carlo-comparison

max_issues_repo_name: pnadon/monte-carlo-comparison

max_forks_repo_name: pnadon/monte-carlo-comparison


text:
# Run examples in BoltzmannMachines.jl with plotting
import BoltzmannMachines
const BMs = BoltzmannMachines
include(joinpath(dirname(pathof(BoltzmannMachines)), "..", "test", "examples.jl"))

using Test
import Gadfly

using BoltzmannMachinesPlots

function test_scatterhidden()
   nsamples = 10
   x = BMs.barsandstripes(10,4)
   rbm = fitrbm(x, epochs = 1)
   dbm = fitdbm(x, epochs = 1)
   labels = rand(["1", "2", "3"], nsamples)
   @test scatterhidden(rbm, x,
         opacity = 0.5, labels = labels) isa Gadfly.Plot
   @test scatterhidden(dbm, x,
         labels = labels, opacity = 0.5) isa Gadfly.Plot
end
test_scatterhidden()

function test_plotevaluation_noribbon()
   monitor, rbm = BMs.monitored_fitrbm(BMs.barsandstripes(10, 9),
      monitoring = monitorloglikelihood!)
   @test plotevaluation(monitor; sdrange = 0.0) isa Gadfly.Plot
end
test_plotevaluation_noribbon()

function test_crossvalidation()
   nsamples = 30;
   nvariables = 4;
   x = barsandstripes(nsamples, nvariables);
   # Determine the optimal number of training epochs for a RBM
   monitor = crossvalidation(x,
         (x, datadict) ->
            begin
               monitors, dbm = BMs.monitored_fitdbm(x, nhiddens = [2, 2],
                     learningrate = 0.1, epochs = 3,
                     monitoring = monitorexactloglikelihood!)
               monitors[end]
            end);
   @test crossvalidationcurve(monitor) isa Gadfly.Plot
   nothing
end
test_crossvalidation()

@test BoltzmannMachinesPlots.plotcurvebundles(BMs.curvebundles(nvariables = 10, nbundles = 3,
      nperbundle = 4, noisesd = 0.03,
      addlabels = true)) isa Gadfly.Plot

@test BoltzmannMachinesPlots.plotcurvebundles(BMs.curvebundles(nvariables = 10, nbundles = 3,
      nperbundle = 4, noisesd = 0.03,
      addlabels = false)) isa Gadfly.Plot


text:
using ArgParse, StatsBase
include("../src/debruijnnets.jl");
include("../src/sampling.jl");
include("../src/motifs.jl");
include("../src/randomgraphs.jl");
include("../src/randomwalks.jl");

function empty_dict(counts)
    return Dict(m=>0 for m in keys(counts))
end

function divergence_by_nodes(G, k, true_counts, true_dist, walk_interval,
                              num_intervals,  all_kwalks, methods)

    nodes, A, A_bin, path_counts = filter_nodes(G, k)
    node_probabilities = compute_node_probs(nodes, path_counts)
    sampled_counts = Dict(key=>empty_dict(true_counts) for key in methods)
    divergences = Dict(key=>Vector{Float64}(undef, num_intervals) for key in keys(sampled_counts))
    for i in range(1, num_intervals)
        for met in methods
            if met == "all_walks"
                walks = uniform_walk_sample(G, k, walk_interval,
                                            -1, true, true)
	    elseif met == "uniform"
		walks = sample(all_kwalks, walk_interval)
            elseif met == "rw"
                walks = random_walks(G, k, walk_interval, nodes, node_probabilities, false, false)
            else
                walks = uniform_walk_sample(G, k, walk_interval,
                                            met, true, true)
            end

            curr_counts = count_motifs(walks, ones(length(walks)))
            for (motif, val) in curr_counts
                sampled_counts[met][motif] = val
            end

            sampled_dist = [Float64(sampled_counts[met][motif]) for motif in keys(true_counts)]
            sampled_dist .+= 10^-12
            sampled_dist ./= sum(sampled_dist)
            div = kldivergence(true_dist, sampled_dist)
            divergences[met][i] = div
        end
    end
    return divergences
end

function write_to_file(output, output_file)
    open(output_file, "w") do file
        for (k, kdict) in output
	    for (m, mtup) in kdict
		out_str = "$k|$(m)|$(join(mtup[1], ","))|$(join(mtup[2], ","))\n"
                write(file, out_str)
            end
        end
    end
    return nothing
end

function get_truth(G, adj, k)
    adjacency_dict = get_adj_dict(G, adj)
    walklist = Vector{Tuple}()
    for node in range(1, nv(G))
        nodewalks = all_walks(adjacency_dict, node, k)
        append!(walklist, nodewalks)
    end
    true_counts = count_motifs(walklist)
    true_dist = [true_counts[motif] for motif in keys(true_counts)]
    true_dist ./= sum(true_dist)
    return true_dist, true_counts
end

function get_all_truths(G, kvals)
    true_dist_dict = Dict()
    true_count_dict = Dict()
    all_kedge_walks = Dict()
    total_walks = Dict()
    adj = adjacency_matrix(G)
    for k in kvals
        true_dist_dict[k], true_count_dict[k] = get_truth(G, adj, k)
        total_walks[k] = sum(adj^k)
        all_kedge_walks[k] = get_all_walks(G, k)
    end
    return true_dist_dict, true_count_dict, total_walks, all_kedge_walks
end

function graph_and_motifs(N, er, p, pa, m, kvals)
   if er
        G = gnp_graph(N, p)
    elseif pa
        G = pa_graph(N, m)
   end
   return G, get_all_truths(G, kvals)
end

function run_gpr(N::Integer, er::Bool, p::Float64, pa::Bool, m::Integer,
                      kvals, methods, walk_interval::Integer, iterations::Integer, num_runs::Integer)
    output = Dict(k=>Dict(method=>Array{Float64}(undef, num_runs, iterations) for method in methods) for k in kvals)
    Threads.@threads for i in range(1, num_runs)
        # Create a new graph for each run
        G, (true_dist_dict, true_count_dict, total_walks, all_kedge_walks) = graph_and_motifs(N, er, p, pa, m, kvals)
        Threads.@threads for k in kvals
            divergences = divergence_by_nodes(G, k, true_count_dict[k], true_dist_dict[k], walk_interval, iterations, all_kedge_walks[k], methods)
            for method in keys(divergences)
                for j in range(1, length(divergences[method]))
                    output[k][method][i,j] = divergences[method][j]
                end
            end
        end
        println(i)
     end
   return output
end

function run_oneg(N::Integer, er::Bool, p::Float64, pa::Bool, m::Integer,
                      kvals, methods, walk_interval::Integer, iterations::Integer, num_runs::Integer)
    output = Dict(k=>Dict(method=>Array{Float64}(undef, num_runs, iterations) for method in methods) for k in kvals)
    # Create a single graph
    G, (true_dist_dict, true_count_dict, total_walks, all_kedge_walks) = graph_and_motifs(N, er, p, pa, m, kvals)
    # Sample from it many times
    Threads.@threads for i in range(1, num_runs)
        Threads.@threads for k in kvals
           divergences = divergence_by_nodes(G, k, true_count_dict[k], true_dist_dict[k], walk_interval, iterations, all_kedge_walks[k], methods)
            for method in keys(divergences)
                for j in range(1, length(divergences[method]))
                   output[k][method][i,j] = divergences[method][j]
                end
            end
        end
        println(i)
    end
    return output
end

function parse_commandline(args)
    s = ArgParseSettings()
    @add_arg_table s begin
        "-k"
	    help = "Order of the motifs."
	    arg_type = Int
        "-N", "--num_nodes"
	    help = "Number of nodes"
	    arg_type = Int
	"--pa"
	    action = :store_true
	    help = "Preferential attachment model."
	"-m"
	    arg_type = Int
	    default = 5
	    help = "Preferential attachment param.
                      Ignored for -e."
	"-e", "--er"
	    action = :store_true
	    help = "ER model."
	"-p"
	    arg_type = Float64
            default = 0.01
	    help = "ER model param. Ignored for -b."
	"-w", "--walk_interval"
	    help = "Number of walks per iteration."
	    arg_type = Int
	"-i", "--iterations"
	    help = "Number of iterations."
	    arg_type = Int
	"-r", "--runs"
	    help = "Number of runs."
	    arg_type = Int
	"--gpr"
	    help = "Compute a new graph every run."
	    action = :store_true
end
return parse_args(args, s)
end

function main()
	arguments = parse_commandline(ARGS)
	println(arguments)
	k = arguments["k"]
	N = arguments["num_nodes"]
	pa = arguments["pa"]
	m = arguments["m"]
	er = arguments["er"]
	p = arguments["p"]
	walk_interval = arguments["walk_interval"]
	iterations = arguments["iterations"]
	runs = arguments["runs"]
	gpr = arguments["gpr"]
	num_cpus = Threads.nthreads()
	kvals = [2, 3, 4]
	methods = ["uniform", "all_walks", "rw"]

	if !gpr
	    output = run_oneg(N, er, p, pa, m, kvals, methods, walk_interval, iterations, runs)
	else
	    output = run_gpr(N, er, p, pa, m, kvals, methods, walk_interval, iterations, runs)
	end

	mean_output = Dict(k=>Dict() for k in keys(output))
	for (k, kdict) in output
	    for (method, method_arr) in kdict
		    mean_output[k][method] = mean_and_std(method_arr, 1)
	    end
	end

	if er
	    model = "er"
	    param = p
	else
	    model = "pa"
	    param = m
	end

	output_filename = "../../debruijn-nets/results/motifs/$(model)_N-$(N)_param-$(param)_interval-$(walk_interval)_iters-$(iterations)_runs-$(runs)_kl"
	if gpr
	    output_filename = output_filename * "_gpr.csv"
	else
	    output_filename = output_filename * "_1g.csv"
	end

	write_to_file(mean_output, output_filename)
end

main()


text:
using BLG_SQUID
using Test

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


text:
"""
    choosereaction(reactionrates::AbstractArray{R}...)where R<:Real

Calculate the reaction time `τ` and randomly choose an reaction from reactions with weights `reactions`.
"""
function choosereaction(reactions::AbstractArray{<:Real}...)
    sum_rs = [sum(i) for i in reactions]
    sum_all = sum(sum_rs)
    τ = -log(rand()) / sum_all
    r_idx = sample(sum_rs) # reaction index
    idx = sample(reactions[r_idx]) # index in reaction
    return τ, r_idx, idx
end

"""
    Dynamics{R<:Real, I}

A struct records changes over time.
"""
struct Dynamics{R <: Real,T}
    history::Vector{Tuple{R,T}}
end

Dynamics(t::Real, n) = Dynamics([(t, n)])

function record!(d::Dynamics, t::Real, n)
    push!(d.history, (t, copy(n)))
    d
end

"""
    zeroatdiag!(m::AbstractMatrix{T}) where T

Make diagonal elements of given matrix `m` zero.

# Example
```jldoctest
julia> m = rand(2, 2)
2×2 Array{Float64,2}:
 0.301516  0.101581
 0.216389  0.672865

julia> zeroatdiag!(m)
2×2 Array{Float64,2}:
 0.0       0.101581
 0.216389  0.0
```
"""
function zeroatdiag!(m::AbstractMatrix{T})where T
    zero_t = zero(T)
    for i in 1:size(m, 1)
        m[i, i] = zero_t
    end
    m
end

"""
    apply(f, vars::Symbol...)

Apply the `f` to  all `vars` in place.

# Example
```jldoctest
julia> @macroexpand @apply copy a b c
quote
    a = copy(a)
    b = copy(b)
    c = copy(c)
end

julia> a = 1
1

julia> @apply x->x+1 a
2

julia> a
2
```
"""
macro apply(f, vars::Symbol...)
    ex = Expr(:block)
    for var in vars
        push!(ex.args, :($var = $f($var)))
    end
    esc(ex)
end

@inline function slicecol(m::AbstractMatrix, i::Integer, colsize::Integer)
    m[(i-1)*colsize+1:i*colsize]
end

@inline slicecol(m::AbstractMatrix, i::Integer) = slicecol(m, i, size(m, 1))

Base.getindex(a::Array, t::Tuple) = Core.arrayref(true, a, t...)
Base.setindex!(a::Array{T}, x::T, t::Tuple) where {T} = Core.arrayset(true, a, x, t...)
Base.setindex!(a::Array{T}, x, t::Tuple) where {T} = Core.arrayset(true, a, convert(T, x)::T, t...)


text:
# Autogenerated wrapper script for lib4ti2_jll for x86_64-w64-mingw32-cxx03
export circuits, exe4ti2gmp, exe4ti2int32, exe4ti2int64, genmodel, gensymm, graver, groebner, hilbert, lib4ti2, lib4ti2int32, lib4ti2int64, libzsolve, markov, minimize, normalform, output, ppi, qsolve, rays, walk, zbasis, zsolve

using GMP_jll
using GLPK_jll
JLLWrappers.@generate_wrapper_header("lib4ti2")
JLLWrappers.@declare_file_product(circuits)
JLLWrappers.@declare_executable_product(exe4ti2gmp)
JLLWrappers.@declare_executable_product(exe4ti2int32)
JLLWrappers.@declare_executable_product(exe4ti2int64)
JLLWrappers.@declare_executable_product(genmodel)
JLLWrappers.@declare_executable_product(gensymm)
JLLWrappers.@declare_file_product(graver)
JLLWrappers.@declare_file_product(groebner)
JLLWrappers.@declare_file_product(hilbert)
JLLWrappers.@declare_library_product(lib4ti2, "lib4ti2gmp-0.dll")
JLLWrappers.@declare_library_product(lib4ti2int32, "lib4ti2int32-0.dll")
JLLWrappers.@declare_library_product(lib4ti2int64, "lib4ti2int64-0.dll")
JLLWrappers.@declare_library_product(libzsolve, "libzsolve-0.dll")
JLLWrappers.@declare_file_product(markov)
JLLWrappers.@declare_file_product(minimize)
JLLWrappers.@declare_file_product(normalform)
JLLWrappers.@declare_executable_product(output)
JLLWrappers.@declare_executable_product(ppi)
JLLWrappers.@declare_file_product(qsolve)
JLLWrappers.@declare_file_product(rays)
JLLWrappers.@declare_file_product(walk)
JLLWrappers.@declare_file_product(zbasis)
JLLWrappers.@declare_executable_product(zsolve)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, GLPK_jll)
    JLLWrappers.@init_file_product(
        circuits,
        "bin\\circuits",
    )

    JLLWrappers.@init_executable_product(
        exe4ti2gmp,
        "bin\\4ti2gmp.exe",
    )

    JLLWrappers.@init_executable_product(
        exe4ti2int32,
        "bin\\4ti2int32.exe",
    )

    JLLWrappers.@init_executable_product(
        exe4ti2int64,
        "bin\\4ti2int64.exe",
    )

    JLLWrappers.@init_executable_product(
        genmodel,
        "bin\\genmodel.exe",
    )

    JLLWrappers.@init_executable_product(
        gensymm,
        "bin\\gensymm.exe",
    )

    JLLWrappers.@init_file_product(
        graver,
        "bin\\graver",
    )

    JLLWrappers.@init_file_product(
        groebner,
        "bin\\groebner",
    )

    JLLWrappers.@init_file_product(
        hilbert,
        "bin\\hilbert",
    )

    JLLWrappers.@init_library_product(
        lib4ti2,
        "bin\\lib4ti2gmp-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        lib4ti2int32,
        "bin\\lib4ti2int32-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        lib4ti2int64,
        "bin\\lib4ti2int64-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libzsolve,
        "bin\\libzsolve-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        markov,
        "bin\\markov",
    )

    JLLWrappers.@init_file_product(
        minimize,
        "bin\\minimize",
    )

    JLLWrappers.@init_file_product(
        normalform,
        "bin\\normalform",
    )

    JLLWrappers.@init_executable_product(
        output,
        "bin\\output.exe",
    )

    JLLWrappers.@init_executable_product(
        ppi,
        "bin\\ppi.exe",
    )

    JLLWrappers.@init_file_product(
        qsolve,
        "bin\\qsolve",
    )

    JLLWrappers.@init_file_product(
        rays,
        "bin\\rays",
    )

    JLLWrappers.@init_file_product(
        walk,
        "bin\\walk",
    )

    JLLWrappers.@init_file_product(
        zbasis,
        "bin\\zbasis",
    )

    JLLWrappers.@init_executable_product(
        zsolve,
        "bin\\zsolve.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Healpix
export std2fwhm, fwhm2std, gaussian_beam

@doc raw"""
    std2fwhm(stddev)

Convert a standard deviation into a Full Width Half Maximum (FWHM) value. The
measure unit for `stddev` can be arbitrary, as the transformation is linear.
See also `fwhm2std`.
"""
std2fwhm(stddev) = 2 * sqrt(2 * log(2)) * stddev

@doc raw"""
    std2fwhm(stddev)

Convert a standard deviation into a Full Width Half Maximum (FWHM) value. The
measure unit for `stddev` can be arbitrary, as the transformation is linear.
See also `std2fwhm`.
"""
fwhm2std(fwhm) = fwhm / (2 * sqrt(2 * log(2)))

@doc raw"""
    gaussian_beam!(beam_map::Healpix.Map{T,O}, angle_std_rad; normalization = 1.0) where {T,O <: Healpix.Order}

Save a map of a circual Gaussian beam in `beam_map`. The width of the beam is
provided through the parameter `angle_std_rad`, which is expressed in radians.
If you want to specify the FWHM instead, use `fwhm2std` like in the following
example:

````julia
import Healpix
nside = 1024
beam_map = Healpix.Map{Float64, RingOrder}(nside)

# Assume 2.5° of FWHM
gaussian_beam!(beam_map, 2.5 |> fwhm2std)
````

See also `gaussian_beam`.
"""
function gaussian_beam!(beam_map::Healpix.Map{T,O}, angle_std_rad; normalization = 1.0) where {T,O <: Healpix.Order}
    for pixidx in 1:beam_map.resolution.numOfPixels
        theta, _ = Healpix.pix2ang(beam_map, pixidx)
        beam_map[pixidx] = normalization * exp(-theta^2 / (2 * angle_std_rad))
    end
end

function gaussian_beam(nside::Integer, angle_std_rad; normalization = 1.0)
    beam_map = Healpix.Map{Float64,Healpix.RingOrder}(nside)
    gaussian_beam!(beam_map, angle_std_rad, normalization = normalization)

    beam_map
end

function gaussian_beam(res::Healpix.Resolution, angle_std_rad; normalization = 1.0)
    gaussian_beam(res.nside, angle_std_rad, normalization = normalization)
end

@doc raw"""
    gaussian_beam(res::Healpix.Resolution, angle_std_rad; normalization = 1.0) where {T,O <: Healpix.Order}
    gaussian_beam(nside::Integer, angle_std_rad; normalization = 1.0) where {T,O <: Healpix.Order}

Save a map of a circual Gaussian beam in `beam_map`. The width of the beam is
provided through the parameter `angle_std_rad`, which is expressed in radians.
If you want to specify the FWHM instead, use `fwhm2std` like in the following
example:

````julia
# Assume 2.5° of FWHM
beam_map = gaussian_beam(1024, 2.5 |> fwhm2std)
````

See also `gaussian_beam!`.
"""
gaussian_beam

struct BeamMoments
    S::Array{Float32,1}
    M::Array{Float32,2}
end

@doc raw"""

    beam_m(beam_map, l, m, n)

Calculate the moment `(l, m, n)` of the beam in the Healpix map `beam_map`. Beam
moments are defined in Appendix A of the paper "Planck 2013 results. LFI
calibration" (Planck collaboration, A&A, 2013),
https://dx.doi.org/10.1051/0004-6361/201321527.

This function does not take into account the factor `N` (normalization of the beam).
You can compute it using `beam_m(beam_map, 0, 0, 0)`.

"""
function beam_m(beam_map, l::Integer, m::Integer, n::Integer)
    result = 0.0
    for idx in 1:length(beam_map)
        θ, ϕ = Healpix.pix2ang(beam_map, idx)

        # We could have used `Healpix.ang2vec`, but we need the value
        # of sinθ in order to calculate dΩ = sinθ⋅dθ⋅dϕ
        sinθ, cosθ = sincos(θ)
        sinϕ, cosϕ = sincos(ϕ)

        x = sinθ * cosϕ
        y = sinθ * sinϕ
        z = cosθ

        result += x^l * y^m * z^n * sinθ * beam_map[idx]
    end
    result
end

function beam_moments(beam_map)
    N = beam_m(beam_map, 0, 0, 0)

    S = [
        beam_m(beam_map, 1, 0, 0) / N,
        beam_m(beam_map, 0, 1, 0) / N,
        beam_m(beam_map, 0, 0, 1) / N,
    ]

    M = [
        [(beam_m(beam_map, 2, 0, 0) / N)  (beam_m(beam_map, 1, 1, 0) / N)  (beam_m(beam_map, 1, 0, 1) / N)];
        [(beam_m(beam_map, 1, 1, 0) / N)  (beam_m(beam_map, 0, 2, 0) / N)  (beam_m(beam_map, 0, 1, 1) / N)];
        [(beam_m(beam_map, 1, 0, 1) / N)  (beam_m(beam_map, 0, 1, 1) / N)  (beam_m(beam_map, 0, 0, 2) / N)];
    ]

    BeamMoments(S, M)
end


text:
# functions for derivative matrices
function stencil1d(m::Int)

    ic = Int(ceil(m/2.0))
    M = zeros((m, m))
    for i=1:m
        for j=1:m
            M[i, j] = (i - ic)^(j-1) / factorial(j-1)
        end
    end
    return inv(M)
end

# NOTE: Minv is the output of running stencil1d()
function regular_derivative_1d(N, n, Minv)

    D = sparse(zeros(Float64, (N, N)))
    stencil = Minv[n+1,:]
    m = length(stencil)
    ic = Int(ceil(m/2.0))

    for i=1:m
        k = i - ic
        vec = stencil[i] * ones(N - abs(k))
        D += spdiagm(k => vec)
    end

    return D
end

function derivative_matrix(nx, ny, Mxinv, Myinv, grd::Grid)

    Dx = regular_derivative_1d(grd._Nx, nx, Mxinv)
    Dy = regular_derivative_1d(grd._Ny, ny, Myinv)

    D = kron(Dy, Dx)
    P = grd.Proj
    Pinv = transpose(grd.Proj)
    D2d = P * D * Pinv / ((grd.dx)^nx *(grd.dy)^ny)

    return kron(sparse(I, grd.Nz, grd.Nz), D2d)
end

function line_average(points::Matrix{Float64}, mx::Int64, my::Int64, MinvT::Matrix{Float64}, grd::Grid)

    N = size(points)[2]
    is = ones(Int64, N*mx*my)
    js = Int64[]
    dats = Float64[]

    for i=1:N
        row_dat, row_j = interpolation_row(points[:,i]..., mx, my, MinvT, grd)
        js = [js; row_j]
        dats = [dats; row_dat]
    end

    return sparse(is, js, dats, 1, grd._Nx*grd._Ny) * transpose(grd.Proj) / N
end