text:
"""
    padding(n::Integer)

Determines the total number of bytes needed to store `n` bytes with padding.
Note that the Arrow standard requires buffers to be aligned to 8-byte boundaries.
"""
padding(n::Integer) = ((n + ALIGNMENT - 1) ÷ ALIGNMENT) * ALIGNMENT

paddinglength(n::Integer) = padding(n) - n

function writezeros(io::IO, n::Integer)
    s = 0
    for i ∈ 1:n
        s += Base.write(io, 0x00)
    end
    s
end

# efficient writing of arrays
writearray(io, col) = writearray(io, maybemissing(eltype(col)), col)

function writearray(io::IO, ::Type{T}, col) where {T}
    if col isa Vector{T}
        n = Base.write(io, col)
    elseif isbitstype(T) && (col isa Vector{Union{T, Missing}} || col isa SentinelVector{T, T, Missing, Vector{T}})
        # need to write the non-selector bytes of isbits Union Arrays
        n = Base.unsafe_write(io, pointer(col), sizeof(T) * length(col))
    elseif col isa ChainedVector
        n = 0
        for A in col.arrays
            n += writearray(io, T, A)
        end
    else
        n = 0
        for x in col
            n += Base.write(io, coalesce(x, ArrowTypes.default(T)))
        end
    end
    return n
end

"""
    getbit

This deliberately elides bounds checking.
"""
getbit(v::UInt8, n::Integer) = Bool((v & 0x02^(n - 1)) >> (n - 1))

"""
    setbit

This also deliberately elides bounds checking.
"""
function setbit(v::UInt8, b::Bool, n::Integer)
    if b
        v | 0x02^(n - 1)
    else
        v & (0xff ⊻ 0x02^(n - 1))
    end
end

"""
    bitpackedbytes(n[, pad=true])

Determines the number of bytes used by `n` bits, optionally with padding.
"""
function bitpackedbytes(n::Integer)
    ℓ = cld(n, 8)
    return ℓ + paddinglength(ℓ)
end

# count # of missing elements in an iterable
nullcount(col) = count(ismissing, col)

# like startswith/endswith for strings, but on byte buffers
function _startswith(a::AbstractVector{UInt8}, pos::Integer, b::AbstractVector{UInt8})
    for i = 1:length(b)
        @inbounds check = a[pos + i - 1] == b[i]
        check || return false
    end
    return true
end

function _endswith(a::AbstractVector{UInt8}, endpos::Integer, b::AbstractVector{UInt8})
    aoff = endpos - length(b) + 1
    for i = 1:length(b)
        @inbounds check = a[aoff] == b[i]
        check || return false
        aoff += 1
    end
    return true
end

# read a single element from a byte vector
# copied from read(::IOBuffer, T) in Base
function readbuffer(t::AbstractVector{UInt8}, pos::Integer, ::Type{T}) where {T}
    GC.@preserve t begin
        ptr::Ptr{T} = pointer(t, pos)
        x = unsafe_load(ptr)
    end
end

# # argh me mateys, don't nobody else go pirating this method
# # this here be me own booty!
# if !applicable(iterate, missing)
# Base.iterate(::Missing, st=1) = st === nothing ? nothing : (missing, nothing)
# end

getK(::Type{Pair{K, V}}) where {K, V} = K
getV(::Type{Pair{K, V}}) where {K, V} = V
getN(::Type{NTuple{N, T}}) where {N, T} = N
getT(::Type{NTuple{N, T}}) where {N, T} = T
getnames(::Type{NamedTuple{names, T}}) where {names, T} = names
getT(::Type{NamedTuple{names, T}}) where {names, T} = T
getN(::Type{NamedTuple{names, T}}) where {names, T} = length(names)

encodingtype(n) = n < div(typemax(Int8), 2) ? Int8 : n < div(typemax(Int16), 2) ? Int16 : n < div(typemax(Int32), 2) ? Int32 : Int64

struct Converter{T, A} <: AbstractVector{T}
    data::A
end

converter(::Type{T}, x::A) where {T, A} = Converter{eltype(A) >: Missing ? Union{T, Missing} : T, A}(x)
converter(::Type{T}, x::ChainedVector{A}) where {T, A} = ChainedVector([converter(T, x) for x in x.arrays])

Base.IndexStyle(::Type{<:Converter}) = Base.IndexLinear()
Base.size(x::Converter) = (length(x.data),)
Base.eltype(x::Converter{T, A}) where {T, A} = T
Base.getindex(x::Converter{T}, i::Int) where {T} = ArrowTypes.arrowconvert(T, getindex(x.data, i))

maybemissing(::Type{T}) where {T} = T === Missing ? Missing : Base.nonmissingtype(T)

macro miss_or(x, ex)
    esc(:($x === missing ? missing : $(ex)))
end

function getfooter(filebytes)
    len = readbuffer(filebytes, length(filebytes) - 9, Int32)
    FlatBuffers.getrootas(Meta.Footer, filebytes[end-(9 + len):end-10], 0)
end

function getrb(filebytes)
    f = getfooter(filebytes)
    rb = f.recordBatches[1]
    return filebytes[rb.offset+1:(rb.offset+1+rb.metaDataLength)]
    # FlatBuffers.getrootas(Meta.Message, filebytes, rb.offset)
end

function readmessage(filebytes, off=9)
    @assert readbuffer(filebytes, off, UInt32) === 0xFFFFFFFF
    len = readbuffer(filebytes, off + 4, Int32)

    FlatBuffers.getrootas(Meta.Message, filebytes, off + 8)
end

# an AbstractVector version of Iterators.flatten
# code based on SentinelArrays.ChainedVector
struct ToList{T, A, I} <: AbstractVector{T}
    data::Vector{A} # A is AbstractVector or AbstractString
    inds::Vector{I}
end

stringtype(T) = false
stringtype(::Type{T}) where {T <: AbstractString} = true

function ToList(input; largelists::Bool=false)
    AT = eltype(input)
    ST = Base.nonmissingtype(AT)
    T = stringtype(ST) ? UInt8 : eltype(ST)
    len = stringtype(ST) ? sizeof : length
    data = AT[]
    I = largelists ? Int64 : Int32
    inds = I[0]
    sizehint!(data, length(input))
    sizehint!(inds, length(input))
    totalsize = I(0)
    for x in input
        if x === missing
            push!(data, missing)
        else
            push!(data, x)
            totalsize += len(x)
            if I === Int32 && totalsize > 2147483647
                I = Int64
                inds = convert(Vector{Int64}, inds)
            end
        end
        push!(inds, totalsize)
    end
    return ToList{T, AT, I}(data, inds)
end

Base.IndexStyle(::Type{<:ToList}) = Base.IndexLinear()
Base.size(x::ToList) = (length(x.inds) == 0 ? 0 : x.inds[end],)

function Base.pointer(A::ToList{UInt8}, i::Integer)
    chunk = searchsortedfirst(A.inds, i)
    return pointer(A.data[chunk - 1])
end

@inline function index(A::ToList, i::Integer)
    chunk = searchsortedfirst(A.inds, i)
    return chunk - 1, i - (@inbounds A.inds[chunk - 1])
end

Base.@propagate_inbounds getchunk(data, chunk) = @inbounds data[chunk]
Base.@propagate_inbounds getchunk(data::Vector{A}, chunk) where {A <: Union{Missing, AbstractString}} = @inbounds codeunits(data[chunk])

Base.@propagate_inbounds function Base.getindex(A::ToList, i::Integer)
    @boundscheck checkbounds(A, i)
    chunk, ix = index(A, i)
    @inbounds x = getchunk(A.data, chunk)[ix]
    return x
end

Base.@propagate_inbounds function Base.setindex!(A::ToList, v, i::Integer)
    @boundscheck checkbounds(A, i)
    chunk, ix = index(A, i)
    @inbounds getchunk(A.data, chunk)[ix] = v
    return v
end

# efficient iteration
@inline function Base.iterate(A::ToList)
    length(A) == 0 && return nothing
    i = 1
    chunk = 2
    chunk_i = 1
    chunk_len = A.inds[chunk]
    while i > chunk_len
        chunk += 1
        chunk_len = A.inds[chunk]
    end
    x = getchunk(A.data, chunk - 1)[1]
    # find next valid index
    i += 1
    if i > chunk_len
        while true
            chunk += 1
            chunk > length(A.inds) && break
            chunk_len = A.inds[chunk]
            i <= chunk_len && break
        end
    else
        chunk_i += 1
    end
    return x, (i, chunk, chunk_i, chunk_len, length(A))
end

@inline function Base.iterate(A::ToList, (i, chunk, chunk_i, chunk_len, len))
    i > len && return nothing
    @inbounds x = getchunk(A.data, chunk - 1)[chunk_i]
    i += 1
    if i > chunk_len
        chunk_i = 1
        while true
            chunk += 1
            chunk > length(A.inds) && break
            @inbounds chunk_len = A.inds[chunk]
            i <= chunk_len && break
        end
    else
        chunk_i += 1
    end
    return x, (i, chunk, chunk_i, chunk_len, len)
end

#
struct ToFixedSizeList{T, N, A} <: AbstractVector{T}
    data::A # A is AbstractVector or AbstractString
end

function ToFixedSizeList(input)
    NT = Base.nonmissingtype(eltype(input)) # NTuple{N, T}
    return ToFixedSizeList{getT(NT), getN(NT), typeof(input)}(input)
end

Base.IndexStyle(::Type{<:ToFixedSizeList}) = Base.IndexLinear()
Base.size(x::ToFixedSizeList{T, N}) where {T, N} = (N * length(x.data),)

Base.@propagate_inbounds function Base.getindex(A::ToFixedSizeList{T, N}, i::Integer) where {T, N}
    @boundscheck checkbounds(A, i)
    a, b = fldmod1(i, N)
    @inbounds x = A.data[a]
    return @inbounds x === missing ? ArrowTypes.default(T) : x[b]
end

# efficient iteration
@inline function Base.iterate(A::ToFixedSizeList{T, N}, (i, chunk, chunk_i, len)=(1, 1, 1, length(A))) where {T, N}
    i > len && return nothing
    @inbounds y = A.data[chunk]
    @inbounds x = y === missing ? ArrowTypes.default(T) : y[chunk_i]
    if chunk_i == N
        chunk += 1
        chunk_i = 1
    else
        chunk_i += 1
    end
    return x, (i + 1, chunk, chunk_i, len)
end

struct ToStruct{T, i, A} <: AbstractVector{T}
    data::A # eltype is NamedTuple
end

ToStruct(x::A, j::Integer) where {A} = ToStruct{fieldtype(Base.nonmissingtype(eltype(A)), j), j, A}(x)

Base.IndexStyle(::Type{<:ToStruct}) = Base.IndexLinear()
Base.size(x::ToStruct) = (length(x.data),)

Base.@propagate_inbounds function Base.getindex(A::ToStruct{T, j}, i::Integer) where {T, j}
    @boundscheck checkbounds(A, i)
    @inbounds x = A.data[i]
    return @miss_or(x, @inbounds getfield(x, j))
end

# Union arrays
# need a custom representation of Union types since arrow unions
# are ordered, and possibly indirected via separate typeIds array
# here, T is Meta.UnionMode.Dense or Meta.UnionMode.Sparse,
# typeIds is a NTuple{N, Int32}, and U is a Tuple{...} of the
# unioned types
struct UnionT{T, typeIds, U}
end

unionmode(::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U} = T
typeids(::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U} = typeIds
Base.eltype(::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U} = U

ArrowTypes.ArrowType(::Type{<:UnionT}) = ArrowTypes.UnionType()

# convenience wrappers for signaling that an array shoudld be written
# as with dense/sparse union arrow buffers
struct DenseUnionVector{T, U} <: AbstractVector{UnionT{Meta.UnionMode.Dense, nothing, U}}
    itr::T
end

DenseUnionVector(x::T) where {T} = DenseUnionVector{T, Tuple{eachunion(eltype(x))...}}(x)
Base.IndexStyle(::Type{<:DenseUnionVector}) = Base.IndexLinear()
Base.size(x::DenseUnionVector) = (length(x.itr),)
Base.eltype(x::DenseUnionVector{T, U}) where {T, U} = UnionT{Meta.UnionMode.Dense, nothing, U}
Base.iterate(x::DenseUnionVector, st...) = iterate(x.itr, st...)
Base.getindex(x::DenseUnionVector, i::Int) = getindex(x.itr, i)

struct SparseUnionVector{T, U} <: AbstractVector{UnionT{Meta.UnionMode.Sparse, nothing, U}}
    itr::T
end

SparseUnionVector(x::T) where {T} = SparseUnionVector{T, Tuple{eachunion(eltype(x))...}}(x)
Base.IndexStyle(::Type{<:SparseUnionVector}) = Base.IndexLinear()
Base.size(x::SparseUnionVector) = (length(x.itr),)
Base.eltype(x::SparseUnionVector{T, U}) where {T, U} = UnionT{Meta.UnionMode.Sparse, nothing, U}
Base.iterate(x::SparseUnionVector, st...) = iterate(x.itr, st...)
Base.getindex(x::SparseUnionVector, i::Int) = getindex(x.itr, i)

# iterate a Julia Union{...} type, producing an array of unioned types
function eachunion(U::Union, elems=nothing)
    if elems === nothing
        return eachunion(U.b, Type[U.a])
    else
        push!(elems, U.a)
        return eachunion(U.b, elems)
    end
end

function eachunion(T, elems)
    push!(elems, T)
    return elems
end

# produce typeIds, offsets, data tuple for DenseUnion
isatypeid(x::T, ::Type{types}) where {T, types} = isatypeid(x, fieldtype(types, 1), types, 1)
isatypeid(x::T, ::Type{S}, ::Type{types}, i) where {T, S, types} = x isa S ? i : isatypeid(x, fieldtype(types, i + 1), types, i + 1)

function todense(::Type{UnionT{T, typeIds, U}}, x) where {T, typeIds, U}
    typeids = typeIds === nothing ? (0:(fieldcount(U) - 1)) : typeIds
    len = length(x)
    types = Vector{UInt8}(undef, len)
    offsets = Vector{Int32}(undef, len)
    data = Tuple(Vector{i == 1 ? Union{Missing, fieldtype(U, i)} : fieldtype(U, i)}(undef, 0) for i = 1:fieldcount(U))
    for (i, y) in enumerate(x)
        typeid = y === missing ? 0x00 : UInt8(typeids[isatypeid(y, U)])
        @inbounds types[i] = typeid
        @inbounds offsets[i] = length(data[typeid + 1])
        push!(data[typeid + 1], y)
    end
    return types, offsets, data
end

# sparse union child array producer
# for sparse unions, we split the parent array into
# N children arrays, each having the same length as the parent
# but with one child array per unioned type; each child
# should include the elements from parent of its type
# and other elements can be missing/default
function sparsetypeids(::Type{UnionT{T, typeIds, U}}, x) where {T, typeIds, U}
    typeids = typeIds === nothing ? (0:(fieldcount(U) - 1)) : typeIds
    len = length(x)
    types = Vector{UInt8}(undef, len)
    for (i, y) in enumerate(x)
        typeid = y === missing ? 0x00 : UInt8(typeids[isatypeid(y, U)])
        @inbounds types[i] = typeid
    end
    return types
end

struct ToSparseUnion{T, A} <: AbstractVector{T}
    data::A
end

ToSparseUnion(::Type{T}, data::A) where {T, A} = ToSparseUnion{T, A}(data)

Base.IndexStyle(::Type{<:ToSparseUnion}) = Base.IndexLinear()
Base.size(x::ToSparseUnion) = (length(x.data),)

Base.@propagate_inbounds function Base.getindex(A::ToSparseUnion{T}, i::Integer) where {T}
    @boundscheck checkbounds(A, i)
    @inbounds x = A.data[i]
    return @inbounds x isa T ? x : ArrowTypes.default(T)
end

# ToMap
struct ToMap{T, A, I} <: AbstractVector{T}
    data::Vector{A} # A is a Vector{KeyValue{K, V}}
    inds::Vector{I}
end

function ToMap(input; largelists::Bool=false)
    AT = eltype(input)
    ST = Base.nonmissingtype(AT)
    PT = eltype(ST)
    K, V = getK(PT), getV(PT)
    T = KeyValue{K, V}
    VT = AT !== ST ? Union{Missing, Vector{T}} : Vector{T}
    data = VT[]
    I = largelists ? Int64 : Int32
    inds = I[0]
    sizehint!(data, length(input))
    sizehint!(inds, length(input))
    totalsize = I(0)
    for x in input
        if x === missing
            push!(data, missing)
        else
            push!(data, [KeyValue{K, V}(k, v) for (k, v) in pairs(x)])
            totalsize += length(x)
            if I === Int32 && totalsize > 2147483647
                I = Int64
                inds = convert(Vector{Int64}, inds)
            end
        end
        push!(inds, totalsize)
    end
    return ToMap{T, VT, I}(data, inds)
end

Base.IndexStyle(::Type{<:ToMap}) = Base.IndexLinear()
Base.size(x::ToMap) = (length(x.inds) == 0 ? 0 : x.inds[end],)

@inline function index(A::ToMap, i::Integer)
    chunk = searchsortedfirst(A.inds, i)
    return chunk - 1, i - (@inbounds A.inds[chunk - 1])
end

Base.@propagate_inbounds function Base.getindex(A::ToMap, i::Integer)
    @boundscheck checkbounds(A, i)
    chunk, ix = index(A, i)
    @inbounds x = A.data[chunk][ix]
    return x
end

Base.@propagate_inbounds function Base.setindex!(A::ToMap, v, i::Integer)
    @boundscheck checkbounds(A, i)
    chunk, ix = index(A, i)
    @inbounds A.data[chunk][ix] = v
    return v
end

# efficient iteration
@inline function Base.iterate(A::ToMap)
    length(A) == 0 && return nothing
    i = 1
    chunk = 1
    chunk_i = 1
    chunk_len = A.inds[chunk + 1]
    while i > chunk_len
        chunk += 1
        @inbounds chunk_len = A.inds[chunk + 1]
    end
    x = A.data[chunk][1]
    # find next valid index
    i += 1
    if i > chunk_len
        while true
            chunk += 1
            chunk > length(A.inds) && break
            @inbounds chunk_len = A.inds[chunk + 1]
            i <= chunk_len && break
        end
    else
        chunk_i += 1
    end
    return x, (i, chunk, chunk_i, chunk_len, length(A))
end

@inline function Base.iterate(A::ToMap, (i, chunk, chunk_i, chunk_len, len))
    i > len && return nothing
    @inbounds x = A.data[chunk][chunk_i]
    i += 1
    if i > chunk_len
        chunk_i = 1
        while true
            chunk += 1
            chunk > length(A.inds) && break
            @inbounds chunk_len = A.inds[chunk + 1]
            i <= chunk_len && break
        end
    else
        chunk_i += 1
    end
    return x, (i, chunk, chunk_i, chunk_len, len)
end

max_stars_repo_name: StefanKarpinski/arrow

max_issues_repo_name: StefanKarpinski/arrow

max_forks_repo_name: StefanKarpinski/arrow


text:
#
# Count the number of hydrogen bonds between two selections
#

# If the coordinates of a PDB file read into the atoms vector is given
# Bonds are needed, because we need to know to which atoms the hydrogens
# are bound to
#

dotp(x,y) = x[1]*y[1]+x[2]*y[2]+x[3]*y[3]
vnorm(x) = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
cosang(x,y) = dotp(x,y) / (vnorm(x)*vnorm(y))

function nhbonds(atoms,bonds,sel1,sel2; d = 3., angle = 20. )

  dhbond2 = d^2
  coshbond = cos(pi*angle/180.)

  n1 = length(sel1)
  n2 = length(sel2)

  nhbonds = 0
  for isel in 1:n1
    i = sel1[isel]
    for jsel in 1:n2
      j = sel2[jsel]
      # If one of the atoms is not electronegative, skip
      if ! ( atoms[i].element in [ "N", "O", "F" ] && atoms[j].element in [ "N", "O", "F" ] )
        continue
      end
      # If the distance is not compatible, skip
      d2 = (atoms[i].coor[1] - atoms[j].coor[1])^2 +
           (atoms[i].coor[2] - atoms[j].coor[2])^2 +
           (atoms[i].coor[3] - atoms[j].coor[3])^2
      if d2 > dhbond2
        continue
      end
      # Check if a hydrogen is bound to atom i
      foundhbond = false
      ibonds = filter( bond -> ( bond.i == i || bond.j == i ), bonds )
      for bond in ibonds
        if atoms[bond.i].element == "H"
          ib = bond.i
        elseif atoms[bond.j].element == "H"
          ib = bond.j
        else
          continue
        end
        v1H = [ atoms[ib].coor[1] - atoms[i].coor[1] ,
                atoms[ib].coor[2] - atoms[i].coor[2] ,
                atoms[ib].coor[3] - atoms[i].coor[3] ]
        vH2 = [ atoms[j].coor[1] - atoms[ib].coor[1] ,
                atoms[j].coor[2] - atoms[ib].coor[2] ,
                atoms[j].coor[3] - atoms[ib].coor[3] ]
        if cosang(v1H,vH2) > coshbond
          nhbonds = nhbonds + 1
          foundhbond = true
          break
        end
      end
      # If the hbond is already found between these to heavy atoms, skip next
      if foundhbond
        continue
      end
      # Check if a hydrogen is bound to atom j
      jbonds = filter( bond -> ( bond.i == j || bond.j == j ), bonds )
      for bond in jbonds
        if atoms[bond.i].element == "H"
          ib = bond.i
        elseif atoms[bond.j].element == "H"
          ib = bond.j
        else
          continue
        end
        v1H = [ atoms[ib].coor[1] - atoms[j].coor[1] ,
                atoms[ib].coor[2] - atoms[j].coor[2] ,
                atoms[ib].coor[3] - atoms[j].coor[3] ]
        vH2 = [ atoms[i].coor[1] - atoms[ib].coor[1] ,
                atoms[i].coor[2] - atoms[ib].coor[2] ,
                atoms[i].coor[3] - atoms[ib].coor[3] ]
        if cosang(v1H,vH2) > coshbond
          nhbonds = nhbonds + 1
          break
        end
      end
    end
  end

  return nhbonds
end

ext: jl


text:
function unwrap(in::AFArray, w1, w2, s1, s2, p1, p2, iscolumn::Bool)
  out = af_array[0]
  af_unwrap(out, in, w1, w2, s1, s2, p1, p2, iscolumn)
  AFArray(out[1])
end

function wrap(in::AFArray, o1, o2, w1, w2, s1, s2, p1, p2, iscolumn::Bool)
  out = af_array[0]
  af_wrap(out, in, o1, o2, w1, w2, s1, s2, p1, p2, iscolumn)
  AFArray(out[1])
end


text:
function LSBcostfun(stego::T,pidx) where {T<:AbstractImageDistortion}
    c0=stego.distortion
    c1=min(trypmone(stego,pidx)...)
    (stego[pidx]&0x1 == 0x0) ? (c0,c1) : (c1,c0)
end

function LSBcostfun2(stego::T,pidx) where {T<:AbstractImageDistortion}
    d = LSBcostfun(stego,pidx)
    [d[1],d[2]]
end

function LSBcostfun3(stego::T,pidx) where {T<:AbstractImageDistortion}
    d0 = stego.distortion
    dp,dm = trypmone(stego,pidx)
    m = mod(stego[pidx],3)
    if m == 0
        !((mod(stego[pidx],3)==0) && (mod(stego[pidx]+1,3)==1) && (mod(stego[pidx]-1,3)==2)) && println(stego[pidx])
        return[d0,dp,dm]
    elseif m == 1
        !((mod(stego[pidx]-1,3)==0) && (mod(stego[pidx],3)==1) && (mod(stego[pidx]+1,3)==2)) && println(stego[pidx])
        return[dm,d0,dp]
    else
        !((mod(stego[pidx]+1,3)==0) && (mod(stego[pidx]-1,3)==1) && (mod(stego[pidx],3)==2)) && println(stego[pidx])
        return[dp,dm,d0]
    end
end

LSBcostfun(stego::T,pidx,q) where {T<:AbstractImageDistortion} = (q==2) ? LSBcostfun2(stego,pidx) : LSBcostfun3(stego,pidx);

function matchlsb!(img::T,mbits,embpath::Array{Int,1}) where {T<:AbstractImageDistortion}
    for i in 1:length(mbits)
        matchlsb!(img,mbits[i],embpath[i])
    end
    return(img)
end

function matchlsb2!(img::T,mbit,pidx::Int) where {T<:AbstractImageDistortion}
    original=img[pidx]
    # println("matchlsb mbit = $mbit")
    if img[pidx]&0x1!=mbit
        if original==255
            img[pidx]=original-1
            return
        end
        if original==0
            img[pidx]=original+1
            return
        end
        (dp,dm)=trypmone(img,pidx)
        if abs(dp-dm)<1e-8
            img[pidx]+=rand([-1,+1])
        elseif dp>dm
            img[pidx]=original-1
        else
            img[pidx]=original+1
        end
    end
    @assert(mod(img[pidx],2) == mbit)
end

function matchlsb3!(img::T,mbit,pidx::Int) where {T<:AbstractImageDistortion}
    original=img[pidx]
    m = mod(original,3)
    if m != mbit
        if mod(original+1,3) == mbit
            img[pidx]=original+1
        elseif mod(original-1,3) == mbit
            img[pidx]=original-1
        else
            error("cannot embed")
        end
    end
    @assert(mod(img[pidx],3) == mbit)
end

matchlsb!(img::T,mbit,pidx::Int) where {T<:AbstractImageDistortion} = matchlsb2!(img,mbit,pidx)
matchlsb!(img::T,mbit,pidx::Int,q) where {T<:AbstractImageDistortion} = (q==2) ? matchlsb2!(img,mbit,pidx) : matchlsb3!(img,mbit,pidx);

readlsb(img::T,embpath) where {T<:AbstractImageDistortion} = readlsbq(img,embpath,2)

function readlsbq(img::T,embpath,q) where {T<:AbstractImageDistortion}
    y=zeros(UInt8,length(embpath))
    for (i,ii) in enumerate(embpath)
        y[i]=Int(mod(img.image[ii],q))
    end
    return(y)
end


text:
#!/usr/local/bin/julia

using RDatasets: dataset
using DecisionTree

iris = dataset("datasets", "iris")

head(iris)

features = Array.(iris[:, 1:4])
labels = Array(iris[:, 5])

# Single Tree
model = build_tree(labels, features)
print_tree(model)

"""
Wrapper for producing a Random Forest model
"""
function RandomForest(x, y, ntrees=10)
  features = Array.(x)
  labels = Array(y)

  nfeatures = size(features)[2] |> sqrt |> floor |> Int
  model = build_forest(labels, features, nfeatures, ntrees)
  return(model)
end

model = RandomForest(iris[:, 1:4], iris[:, 5], 50)

# Predict on new data
newdata = [5.9,3.0,5.1,1.9]
apply_forest(model, newdata)

# Estimate predicted probabilities for each class
apply_forest_proba(model, newdata, ["setosa", "versicolor", "virginica"])

# Cross-validation for random forests
accuracy = nfoldCV_forest(labels, features, 2, 10, 3, 0.5)

max_stars_repo_name: ellisvalentiner/julia-talk

max_issues_repo_name: ellisvalentiner/julia-talk

max_forks_repo_name: ellisvalentiner/julia-talk


text:
function mean( img::AbstractArray{T,2} ) where {T<:Real}

	ret = 0.0
	N = length( img );

	@simd for li in 1:N
		@inbounds ret += img[li]
	end

	return ret/N
end

function variance( img::AbstractArray{T,2} ) where {T<:Real}

	mean = mean( img )

	ret = 0.0
	N = length( img )

	@simd for li in 1:N
		@inbounds ret += ( img[li] - mean )*( img[li] - mean )
	end

	return ret/N
end

function variance( img::AbstractArray{T,2}, mean ) where {T<:Real}

	ret = 0.0
	N = length( img )

	@simd for li in 1:N
		@inbounds ret += ( img[li] - mean )*( img[li] - mean )
	end

	return ret/N
end

function std( img::AbstractArray{T,2} ) where {T<:Real}

	return sqrt( variance( img ) )
end

function std( img::AbstractArray{T,2}, mean ) where {T<:Real}

	return sqrt( variance( img, mean ) )
end

function cov( img1::Array{T1,2}, img2::Array{T2,2}) where {T1<:Real,T2<:Real}

	length( img1 ) !== length( img2 ) && ( error("Arrays have different lengths") )

	m1 = mean( img1 )
	s1 = std( img1, m1 )

	m2 = mean( img2 )
	s2 = std( img2, m2 )

	ret = 0.0
	N = length( img1 )

	@simd for li in 1:N
		@inbounds ret += ( img1[li] - m1 )*( img2[li] - m2 )
	end

	return ret/N
end

function correlation( img1::Array{T1,2}, img2::Array{T2,2} ) where {T1<:Real,T2<:Real}

	length( img1 ) !== length( img2 ) && ( error("Arrays have different lengths") )

	m1 = mean( img1 )
	s1 = std( img1, m1 )

	m2 = mean( img2 )
	s2 = std( img2, m2 )

	ret = 0.0
	N = length( img1 )

	@simd for li in 1:N
		@inbounds ret += ( img1[li] - m1 )*( img2[li] - m2 )
	end

	return ret/(N*s1*s2)
end

function crossCorrelationZNCC( f::Array{T,2}, g::Array{T,2} ) where {T<:Real}

		lN = div( length( f ), 16 );

		# statstics for ZNCC
		#=
		meanf  = mean( f );
		meanf2 = meanf * meanf;

		meang  = mean( g );
		meang2 = meang * meang;
		=#

		# integral Arrays to compute stdf at each possible translation
		sumf  = integralArray( f );
		sumf2 = integralArray( f, fun=(x)->(x*x) );
		sumg  = integralArray( g );
		sumg2 = integralArray( g, fun=(x)->(x*x) );

	    # correlation size = N + M - 1
	    csize = size(f) .+ size(g) .- 1;
	    h, w  = csize;
		num   = zeros( Float32, csize ); # we compute the numerator with FFT for each cross-correlation element
		corr  = zeros( Float32, csize );

		#=
		padf  = zeros( Complex{T}, csize );
		padf[ 1:size(f,1), 1:size(f,2) ] .= complex.( f .- meanf, 0.0 );

		padg  = zeros( Complex{T}, csize );
		padg[ 1:size(g,1), 1:size(g,2) ] .= complex.( g .- meang, 0.0 );

		plan  =  plan_fft!( padf );
		iplan = plan_ifft!( padf );

		# cross-correlation ( f - meanf ) and ( g - meang ) with FFT
		plan * padf;
		plan * padg;
		@inbounds @simd for e in 1:length(padf)
			padf[e] = conj( padf[e] ) * padg[e];
		end

		iplan * padf;
		@inbounds @simd for e in 1:length(padf)
			corr[e] = real( padf[e] )
		end

		shifts = div.( csize, 2 ).- 1
	    Base.circshift!( num, corr, shifts )
		=#

		# computing the denominator
		# Computing denominator through summed-area tables.
	    @inbounds for col in 1:csize[2]

			minCol_g = 1 + max( 0, col - size(f,2) )
			maxCol_g = min( col, size(g,2) )

			#
			minCol_f = max( 1, 1 + size(f,2) - col )
	        maxCol_f = size(f,2) - max( 0, col - size(g,2) )

			lenCol_g = length( minCol_g:maxCol_g );

	        for row in 1:csize[1]

				minRow_g = 1 + max( 0, row - size(f,1) );
            	maxRow_g = min( row, size(g,1) );

				#
				minRow_f = max( 1, 1 + size(f,1) - row )
	            maxRow_f = size(f,1) - max( 0, row - size(g,1) )

				lenRow_g = length( minRow_g:maxRow_g )
				N = lenRow_g*lenCol_g

				#
				F, F2 = integralAreaZNCC( sumf, sumf2, row, col, size(f), size(g) )
				G, G2 = integralAreaZNCC( sumg, sumg2, row, col, size(g), size(g) )

				#
				meanf  = F/N;
				meanf2 = meanf * meanf;
				stdf   = sqrt( F2 + N*meanf2 - 2*F*meanf );
				meang  = G/N;
				meang2 = meang * meang;
				stdg   = sqrt( G2 + N*meang2 - 2*G*meang );

				#
				subg = view( g, minRow_g:maxRow_g, minCol_g:maxCol_g )
				subf = view( f, minRow_f:maxRow_f, minCol_f:maxCol_f )

				#
				zncc = 0.0
				for e in 1:length(subg)
					zncc += ( subf[e] - meanf )*( subg[e] - meang )
				end

				#
				corr[ row, col ] = zncc/(stdg*stdf)

				#=
				if N < lN
					corr[ row, col ] = 0.0
				else
		            F, F2 = integralAreaZNCC( sumf, sumf2, row, col, size(f), size(g) )
					G, G2 = integralAreaZNCC( sumg, sumg2, row, col, size(g), size(g) )

		            corr[ row, col ] = num[ row, col ]/( sqrt( F2 + N*meanf2  - 2*F*meanf )*sqrt( G2 + N*meang2  - 2*G*meang ) );
				end
				=#
	        end
	    end

		return corr;
end


text:
"""
Module: SymbolicDiff (Symbolic Operation for Arithmetic)
"""

"""
AbstractSymbolic

An abstract structure for symbolic expr.

AbstractSymbolic has the filed
- params::Set{Symbol} A set of Symbol included in expr
"""

abstract type AbstractSymbolic{Tv} end

"""
AbstractNumberSymbolic
AbstractVectorSymbolic
AbstractMatrixSymbolic
AbstractGeneralSymbolic

They represent the types of zero values when the derivative becomes zero
"""

abstract type AbstractNumberSymbolic{Tv} <: AbstractSymbolic{Tv} end
abstract type AbstractVectorSymbolic{Tv} <: AbstractSymbolic{Tv} end
abstract type AbstractMatrixSymbolic{Tv} <: AbstractSymbolic{Tv} end
abstract type AbstractGeneralSymbolic{Tv} <: AbstractSymbolic{Tv} end

"""
SymbolicValue

A constant value
"""

mutable struct SymbolicValue{Tv} <: AbstractNumberSymbolic{Tv}
    params::Set{Symbol}
    val::Tv
end

SymbolicValue(value::Tv) where Tv = SymbolicValue(Set{Symbol}([]), value)

"""
SymbolicVariable

A variable to be derivatived
"""

mutable struct SymbolicVariable{Tv} <: AbstractNumberSymbolic{Tv}
    params::Set{Symbol}
    var::Symbol
end

"""
SymbolicExpr

An expr
"""

mutable struct SymbolicExpression{Tv} <: AbstractNumberSymbolic{Tv}
    params::Set{Symbol}
    op::Symbol
    args::Vector{<:AbstractSymbolic{Tv}}
    others
end

"""
convert
"""

function Base.convert(::Type{<:AbstractNumberSymbolic{T}}, x::SymbolicVariable{S}) where {T<:Number,S<:Number}
    symbolic(x.var, T)
end

function Base.convert(::Type{<:AbstractNumberSymbolic{T}}, x::SymbolicValue{S}) where {T<:Number,S<:Number}
    symbolic(T(x.val), T)
end

function Base.convert(::Type{<:AbstractNumberSymbolic{T}}, x::S) where {T<:Number,S<:Number}
    symbolic(T(x), T)
end

function Base.convert(::Type{<:AbstractNumberSymbolic{T}}, x::SymbolicExpression{S}) where {T<:Number,S<:Number}
    if T != S
        SymbolicExpression{T}(x.params, x.op, [convert(AbstractNumberSymbolic{T}, u) for u = x.args], [])
    else
        x
    end
end

"""
promotion
"""

function Base.promote_rule(::Type{<:AbstractNumberSymbolic{T}}, ::Type{S}) where {T<:Number,S<:Number}
    AbstractNumberSymbolic{promote_type(T,S)}
end

function Base.promote_rule(::Type{<:AbstractNumberSymbolic{T}}, ::Type{<:AbstractNumberSymbolic{S}}) where {T<:Number,S<:Number}
    AbstractNumberSymbolic{promote_type(T,S)}
end

"""
iszero
"""

function Base.iszero(x::AbstractNumberSymbolic{Tv}) where Tv
    false
end

function Base.iszero(x::SymbolicValue{Tv}) where Tv
    Base.iszero(x.val)
end

function Base.zero(::Type{AbstractNumberSymbolic{Tv}}) where Tv
    SymbolicValue(Tv(0))
end

"""
IO for structures
"""

function Base.show(io::IO, x::AbstractNumberSymbolic)
    expr = _toexpr(x)
    Base.show(io, expr)
end

function _toexpr(x::SymbolicExpression)
    args = [_toexpr(e) for e = x.args]
    Expr(:call, x.op, args...)
end

function _toexpr(x::SymbolicVariable)
    x.var
end

function _toexpr(x::SymbolicValue)
    x.val
end

"""
symbolicexpr(expr)

Build a SymbolicExpression
"""

const operations = [:+, :-, :*, :/, :^, :exp, :sqrt, :log]

function symbolic(expr::Expr, ::Type{Tv} = Float64) where Tv
    if Meta.isexpr(expr, :call) && expr.args[1] in operations
        args = [symbolic(x, Tv) for x = expr.args[2:end]]
        eval(Expr(:call, expr.args[1], args...))
    elseif Meta.isexpr(expr, :vect)
        vec = [symbolic(x, Tv) for x = expr.args]
        symbolic(vec, Tv)
    elseif Meta.isexpr(expr, :vcat)
        elem = [Expr(:row, [symbolic(y, Tv) for y = x.args]...) for x = expr.args]
        mat = eval(Expr(:vcat, elem...))
        symbolic(mat, Tv)
    else
        nothing
    end
end

function symbolic(expr::Symbol, ::Type{Tv} = Float64) where Tv
    SymbolicVariable{Tv}(Set([expr]), expr)
end

function symbolic(expr::Tv, ::Type{Tx}) where {Tv,Tx}
    SymbolicValue(expr)
end

function symbolic(expr::Tv) where Tv
    SymbolicValue(expr)
end

function symbolic(expr::Nothing, ::Type{Tv}) where Tv
    nothing
end

function symbolic(expr::Nothing) where Tv
    nothing
end


text:
"""
`Basics.read()`  ... reads in data from different files and sources.

  + `("CSF list: Grasp92", "cslFilename")`
    ... reads in the CSF list as given by the (existing) .csl file cslFilename; an basis::Basis is returned with
        basis.isDefined = true and with a proper sequence of orbitals, but with basis.orbitals = Orbital[].

  + `("orbital list: Grasp92", "orbFilename")`
    ... read in the orbitals from a (formatted) .rwf file orbFilename; a list of orbitals::Array{OrbitalsR,1}
        is returned with all subfields specified.
"""
function Basics.read(sa::String, filename::String)
    if      sa == "CSF list: Grasp92"                          wa = readCslFileGrasp92(filename)
    elseif  sa == "orbital list: Grasp92"                      wa = readOrbitalFileGrasp92(filename)
    elseif  sa == "energies & mixing coefficients: Grasp92"    wa = readMixingFileGrasp92(filename)
    else    error("Unsupported keystring = $sa ")
    end

    return( wa )
end

"""
Basics.readCslFileGrasp92(filename::String)`
    ... reads in the CSF list from a Grasp92 .csl file; a basis::Basis is returned.
"""
function Basics.readCslFileGrasp92(filename::String)
    coreSubshells = Subshell[];    peelSubshells =  Subshell[]

    f  = open(filename)
    sa = readline(f);   sa[1:14] != "Core subshells"   &&   error("Not a Grasp92 .cls file.")
    sa = readline(f);   if  length(sa) > 3    go = true    else  go = false    end
    i  = -1
    while go
        i = i + 1;    sh = JAC.subshellGrasp( strip(sa[5i+1:5i+5]) );    push!(coreSubshells, sh)
        if  5i + 6 >  length(sa)    break    end
    end
    sa = readline(f);   println("sa = $sa")
    sa = readline(f);   if  length(sa) > 3    go = true    else  go = false    end
    i  = -1
    while go
        i = i + 1;    sh = JAC.subshellGrasp( strip(sa[5i+1:5i+5]) );    push!(peelSubshells, sh)
        if  5i + 6 >  length(sa)    break    end
    end
    # Prepare a list of all subshells to define the common basis
    subshells = deepcopy(coreSubshells);    append!(subshells, peelSubshells)

    # Now, read the CSF in turn until the EOF
    readline(f)
    NoCSF = 0;    csfs = CsfR[]
    Defaults.setDefaults("relativistic subshell list", subshells)

    while true
        sa = readline(f);      if length(sa) == 0   break    end
        sb = readline(f);      sc = readline(f)
        NoCSF = NoCSF + 1
        push!(csfs, JAC.ManyElectron.CsfRGrasp92(subshells, coreSubshells, sa, sb, sc) )
    end

    println("  ... $NoCSF CSF read in from Grasp92 file  $filename)")
    NoElectrons = sum( csfs[1].occupation )

    Basis(true, NoElectrons, subshells, csfs, coreSubshells, Dict{JAC.Subshell,JAC.Radial.Orbital}() )
end

"""
`Basics.readOrbitalFileGrasp92(filename::String, grid.JAC.Radial.Grid)`
    ... reads in the orbitals list from a Grasp92 .rwf or .out file; a dictionary
        orbitals::Dict{JAC.Subshell,JAC.Radial.Orbital} is returned.
"""
function Basics.readOrbitalFileGrasp92(filename::String, grid::JAC.Radial.Grid)
    # using FortranFiles
    f = FortranFile(filename)

    orbitals = Dict{JAC.Subshell,JAC.Radial.Orbital}();    first = true

    if (String(Base.read(f, FString{6})) != "G92RWF")
	  close(f)
	  error("File \"", filename, "\" is not a proper Grasp92 RWF file")
    end

    Defaults.setDefaults("standard grid", grid);

    while true
	  try
		n, kappa, energy, mtp = Base.read(f, Int32, Int32, Float64, Int32)
		pz, P, Q = Base.read(f, Float64, (Float64, mtp), (Float64, mtp))
		ra = Base.read(f, (Float64, mtp))
            # Now place the data into the right fields of Orbital
            subshell = Subshell( Int64(n), Int64(kappa) )
            if  energy < 0   isBound = true    else    isBound = false   end
            useStandardGrid = true
            println("JAC.Basics.readOrbitalFileGrasp92-aa: warning ... Radial functions still need to be interpolated upon the given grid.")
            Px = P
            Qx = Q

            orbitals = Base.merge( orbitals, Dict( subshell => Orbital(subshell, isBound, useStandardGrid, energy, Px, Qx, JAC.Radial.Grid() ) ))
	  catch ex
		if     ex isa EOFError  break
		else   throw(ex)
		end
	  end
    end
    close(f)

    return( orbitals )
end

"""
`Basics.readMixFileRelci(filename::String)`
    ... reads in the mixing coefficients from a RELCI .mix file; a multiplet::Multiplet is returned.
"""
function Basics.readMixFileRelci(filename::String, basis::Basis)
    levels = Level[];    name =  "from Relci"

    f  = open(filename)
    sa = readline(f);   sa[1:14] != "G92MIX (format"   &&   error("Not a formatted relci.mix file.")
    sa = readline(f);   NoElectrons = Base.parse(Int64, strip(sa[1:6]) );    NoCsf = Base.parse(Int64, strip(sa[7:12]) )
                        NoSubshells = Base.parse(Int64, strip(sa[13:18]) )
    sa = readline(f);   NoLevels    = Base.parse(Int64, strip(sa[1:6]) )
    #
    sa = readline(f)
    sa = readline(f);   list2J = AngularJ64[AngularJ64(0)  for  i = 1:NoLevels];    listp = JAC.Parity[JAC.plus  for  i = 1:NoLevels]
    for  i = 1:NoLevels
        jj = Base.parse(Int64, strip(sa[(i-1)*6+1:(i-1)*6+5]) );   list2J[i] = AngularJ64(jj//2)
        pp = string( sa[(i-1)*6+6] );                              listp[i]  = JAC.Parity( pp )
    end
    #
    sa   = readline(f);   listEnergy = Float64[0.  for  i = 1:NoLevels]
    enav = Base.parse(Float64, strip(sa[1:26]) )
    for  i = 1:NoLevels
        en = Base.parse(Float64, strip(sa[26i+1:26i+26]) );   listEnergy[i] = en
    end
    #
    eigenvectors = zeros(NoCsf,NoLevels)
    for  j = 1:NoCsf
        sa   = readline(f)
        for  i = 1:NoLevels
           cm = Base.parse(Float64, strip(sa[16*(i-1)+1:16(i-1)+16]) );   eigenvectors[j,i] = cm
        end
    end
    #
    for  j = 1:NoLevels
        mc = Float64[0.  for  i = 1:NoCsf]
        for  i = 1:NoCsf    mc[i] = eigenvectors[i,j]    end
        Jx = list2J[j];   Mx = AngularM64( Jx.num//Jx.den )
        push!( levels, Level(Jx, Mx, listp[j], j, enav+listEnergy[j], 0., true, basis, mc) )
    end

    multiplet = JAC.Multiplet(name, levels)
    return( multiplet )
end

"""
`Basics.readMixingFileGrasp92(filename::String, basis::Basis)`
    ... reads in the mixing coefficients from a Grasp92 .mix file by using the given basis; A multiplet::Multiplet
        is returned.
"""
function Basics.readMixingFileGrasp92(filename::String, basis::Basis)
    name = "from Grasp"
    f = FortranFile(filename)

    if (String(Base.read(f, FString{6})) != "G92MIX")
        close(f)
	  error("File \"", filename, "\" does not seem to be a Grasp92 MIX file")
    end

    no_electrons, no_csf, no_subshells = read(f, Int32, Int32, Int32)
    no_asf = Base.read(f, Int32)
    level_numbers = Base.read(f, (Int32, no_asf))
    ijip = Base.read(f, (Int32, 2, no_asf))

    mjs = ijip[1,:] # 2j+1
    mparities = ijip[2,:] # 1=='+', -1=='-'

    mixfile_to_j(i::Int32) = J((i-1)//2)
    mixfile_to_p(i::Int32) = i==1 ?  JAC.plus : JAC.minus

    js = mixfile_to_j.(mjs)
    parities = mixfile_to_p.(mparities)

    avg_energy, eval = read(f, Float64, (Float64, no_asf))
    energies = avg_energy + eval
    eigenvectors = read(f, (Float64, no_csf, no_asf))

    ##x @show no_electrons
    ##x @show no_csf
    ##x @show no_subshells
    ##x @show no_asf
    ##x println("angular momenta = ", js)
    ##x @show parities
    ##x @show energies
    ##x println("first eigenvector = ", eigenvectors[:,1])

    close(f)
    multiplet = JAC.Multiplet(name, JAC.Level[] )
    return( multiplet )
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

find maximum volume hypercube with edges parallel to the axes inside a polyhedron
=#

struct MaxVolumeNative{T <: Real} <: ExampleInstanceNative{T}
    n::Int
    use_hypogeomean::Bool # use hypogeomean cone for geomean objective
    use_power::Bool # use power cones for geomean objective
    use_epipersquare::Bool # use epipersquare cones for geomean objective
end

function build(inst::MaxVolumeNative{T}) where {T <: Real}
    @assert +(inst.use_hypogeomean, inst.use_power, inst.use_epipersquare) == 1
    n = inst.n

    poly_hrep = Matrix{T}(I, n, n)
    poly_hrep .+= T.(randn(n, n)) / n
    c = vcat(-1, zeros(T, n))
    A = hcat(zeros(T, n), poly_hrep)
    b = ones(T, n)

    log_floor(n, i) = (n <= 2^i ? i : log_floor(n, i + 1))
    function log_floor(n::Integer)
        @assert n > zero(n)
        return log_floor(n, zero(n))
    end

    if inst.use_hypogeomean
        G = -Matrix{T}(I, n + 1, n + 1)
        h = zeros(T, n + 1)
        cones = Cones.Cone{T}[Cones.HypoGeoMean{T}(1 + n)]
    elseif inst.use_power
        @assert n > 3 # power cone formulation minimum
        cones = Cones.Cone{T}[]
        # n - 1 3-dim power cones needed, number of new variables is n - 2
        len_power = 3 * (n - 1)
        G_geo_orig = zeros(T, len_power, n)
        G_geo_newvars = zeros(T, len_power, n - 2)
        c = vcat(c, zeros(T, n - 2))
        A = hcat(A, zeros(T, n, n - 2))

        # first cone is a special case since two of the original variables
        # participate in it
        G_geo_orig[1, 1] = -1
        G_geo_orig[2, 2] = -1
        G_geo_newvars[3, 1] = -1
        push!(cones, Cones.GeneralizedPower{T}(fill(inv(T(2)), 2), 1))
        offset = 4
        # loop over new vars
        for i in 1:(n - 3)
            G_geo_newvars[offset + 2, i + 1] = -1
            G_geo_newvars[offset + 1, i] = -1
            G_geo_orig[offset, i + 2] = -1
            push!(cones, Cones.GeneralizedPower{T}([inv(T(i + 2)),
                T(i + 1) / T(i + 2)], 1))
            offset += 3
        end
        # last row also special because hypograph variable is involved
        G_geo_orig[offset, n] = -1
        G_geo_newvars[offset + 1, n - 2] = -1
        G = [
            vcat(zeros(T, len_power - 1), -one(T))  G_geo_orig  G_geo_newvars
            ]
        push!(cones, Cones.GeneralizedPower{T}([inv(T(n)), T(n - 1) / T(n)], 1))
        h = zeros(T, 3 * (n - 1))
    else
        @assert inst.use_epipersquare == true
        # number of variables inside geometric mean is n
        # number of layers of variables
        num_layers = log_floor(n)
        # number of new variables = 1 + 2 + ... + 2^(l - 1) = 2^l - 1
        num_new_vars = 2 ^ num_layers - 1

        c = vcat(c, zeros(T, num_new_vars))
        A = hcat(A, zeros(T, n, num_new_vars))
        rtfact = sqrt(T(2) ^ num_layers)
        # excludes original hypograph variable, padded later
        G_rsoc = zeros(T, 3 * num_new_vars, n + num_new_vars)
        cones = Cones.Cone{T}[]

        offset = offset_next = 0
        row = 1
        # loop over layers, layer 1 describes hypograph variable
        for i in 1:(num_layers - 1)
            num_lvars = 2 ^ (i - 1)
            offset_next = offset + num_lvars
            # loop over variables in each layer
            for j in 1:num_lvars
                G_rsoc[row, n + offset_next + 2j - 1] = -1
                G_rsoc[row + 1, n + offset_next + 2j] = -1
                G_rsoc[row + 2, n + offset + j] = -1
                push!(cones, Cones.EpiPerSquare{T}(3))
                row += 3
            end
            offset = offset_next
        end

        for j in 1:(2 ^ (num_layers - 1))
            # in the last layer, we use the original variables;
            # if beyond the number of variables in the actual
            # geometric mean, adding the buffer variable
            if 2j - 1 > n
                G_rsoc[row, n + 1] = -inv(rtfact)
            else
                G_rsoc[row, 2j - 1] = -1
            end
            if 2j > n
                G_rsoc[row + 1, n + 1] = -inv(rtfact)
            else
                G_rsoc[row + 1, 2j] = -1
            end
            G_rsoc[row + 2, n + offset + j] = -1
            push!(cones, Cones.EpiPerSquare{T}(3))
            row += 3
        end

        # account for original hypograph variable
        G = [
            zeros(T, 3 * num_new_vars)  G_rsoc;
            one(T)  zeros(T, 1, n)  -inv(rtfact)  zeros(T, 1, num_new_vars - 1);
            ]
        push!(cones, Cones.Nonnegative{T}(1))
        h = zeros(T, 3 * num_new_vars + 1)
    end

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Actors, Printf, Random, .Threads
import Actors: spawn

struct Player{S,T}
    name::S
    capa::T
end

struct Ball{T,S,L}
    diff::T
    name::S
    from::L
end

struct Serve{L}
    to::L
end

function (p::Player)(prn, b::Ball)
    if p.capa ≥ b.diff
        send(b.from, Ball(rand(), p.name, self()))
        send(prn, p.name*" serves "*b.name)
    else
        send(prn, p.name*" looses ball from "*b.name)
    end
end
function (p::Player)(prn, s::Serve)
    send(s.to, Ball(rand(), p.name, self()))
    send(prn, p.name*" serves ")
end

@threads for i in 1:nthreads()
    Random.seed!(2021+threadid())
end

prn = spawn(s->print(@sprintf("%s\n", s))) # a print server
ping = spawn(Player("Ping", 0.8), prn, thrd=3)
pong = spawn(Player("Pong", 0.75), prn, thrd=4)

send(ping, Serve(pong))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "SDL2"
version = v"2.0.12"

# Collection of sources required to build SDL2
sources = [
    ArchiveSource("http://www.libsdl.org/release/SDL2-$(version).tar.gz",
                  "349268f695c02efbc9b9148a70b85e58cefbbf704abd3e91be654db7f1e2c863"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SDL2-*/

FLAGS=()
if [[ "${target}" == *-linux-* ]] || [[ "${target}" == *-freebsd* ]]; then
    FLAGS+=(--with-x)
fi

export CPPFLAGS="-I${prefix}/include"
export LDFLAGS="-L${libdir}"

./configure --prefix=${prefix} --host=${target} \
    --enable-shared \
    --disable-static \
    "${FLAGS[@]}"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libSDL2", "SDL2"], :libsdl2)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Xorg_libX11_jll"),
    Dependency("Xorg_libXcursor_jll"),
    Dependency("Xorg_libXext_jll"),
    Dependency("Xorg_libXinerama_jll"),
    Dependency("Xorg_libXrandr_jll"),
    Dependency("Xorg_libXScrnSaver_jll"),
    Dependency("Libglvnd_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using UUIDs: uuid4
using Genie, Genie.Router, Genie.Requests, Genie.Renderer.Json

# Generate a globally unique address for this node
const NODE_ID = replace(string(uuid4()), "-"=>"")

# Instatiate the blockchain
snowchain = Blockchain()

route("/mine") do
	(:message => "We'll mine a new block") |> json
	last_block = snowchain[end]
	last_proof = last_block.proof
	last_hash = shash(last_block)
	proof = proof_of_work(last_proof, last_hash)

	new_transaction!(snowchain, "", NODE_ID, REWARD)

	# Add the block to the chain
	block = new_block(proof, last_hash)

	response = Dict(
		:message => "New block mined",
		:index => block.index,
		:transactions => block.transactions,
		:proof => block.proof,
		:previous_hash => last_hash
	)
	return json(response, status=200)
end

route("transactions/new", method=POST) do
	required = ["sender", "recipient", "amount"]
	@show rawpayload()
	fields = jsonpayload()
	if fields == nothing
		return json(:error => "Format request as JSON object with keys " * join(required, ", "), status=400)
	end
	@show fields

	if !(required ⊆ keys(fields))
		return json(:error => "Missing keys" => setdiff(required, keys(fields)), status=400)
	end

	# Create a new transaction
	index = new_transaction!(snowchain, fields["sender"], fields["recipient"], fields["amount"])

	return json(:message => "Transaction will be added to block $index", status=201)
end

route("/chain") do
	Dict(:chain => snowchain.chain, :length => length(snowchain)) |> json
end

# Start web server at http://127.0.0.1:5000
up(5000)

max_stars_repo_name: derikk/Snowcoin

max_issues_repo_name: derikk/Snowcoin

max_forks_repo_name: derikk/Snowcoin


text:
using MittagLefflerFunctions
import MittagLefflerFunctions.Quadrature.MLQuad
using PyPlot

α_vals = [ 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0 ]
β = 1.0
G(t) = (1-t)/(1+t)
t = range(-1, 1, length=400)
x = G.(t)
Eα = Vector{MLQuad}(undef, 10)

style = ("-", "--", ":", "-.", "-")

figure(1)
for j = 1:5
    α = α_vals[j]
    Eα[j] = MLQuad(α, β, 12, :hyperbola)
    plot(t, Eα[j].(-x), style[j], label="α = $α")
end
grid(true)
legend()
xlabel(L"$t$", fontsize=12)
savefig("figure12L.pdf")

figure(2)
s = range(0, 1, length=400)
t = 2 * (s.^2) .- 1
x = G.(t)
for j = 6:10
    α = α_vals[j]
    Eα[j] = MLQuad(α, β, 12, :hyperbola)
    plot(t[2:end], Eα[j].(-x[2:end]), style[j-5], label="α = $α")
end
grid(true)
legend()
xlabel(L"$t$", fontsize=12)
savefig("figure12R.pdf")


text:
abstract type AbstractVectorBounds end

function set_from_constants(
    b::AbstractVectorBounds,
    ::Type{<:MOI.EqualTo},
    index,
)
    return MOI.EqualTo(b.lower[index])
end

function set_from_constants(
    b::AbstractVectorBounds,
    S::Type{<:Union{MOI.GreaterThan,MOI.EqualTo}},
    index,
)
    # Lower and upper bounds are equal for `EqualTo`, we can take either of them.
    return S(b.lower[index])
end

function set_from_constants(
    b::AbstractVectorBounds,
    S::Type{<:MOI.LessThan},
    index,
)
    return S(b.upper[index])
end

function set_from_constants(
    b::AbstractVectorBounds,
    S::Type{<:Union{MOI.Interval,MOI.Semicontinuous,MOI.Semiinteger}},
    index,
)
    return S(b.lower[index], b.upper[index])
end

function set_from_constants(
    ::AbstractVectorBounds,
    S::Type{<:Union{MOI.Integer,MOI.ZeroOne}},
    index,
)
    return S()
end

"""
    SUPPORTED_VARIABLE_SCALAR_SETS{T}

The union of scalar sets for `VariableIndex` constraints supported by
`Utilities.Hyperrectangle` and `Utilities.VariablesContainer`.
"""
const SUPPORTED_VARIABLE_SCALAR_SETS{T} = Union{
    MOI.EqualTo{T},
    MOI.GreaterThan{T},
    MOI.LessThan{T},
    MOI.Interval{T},
    MOI.Integer,
    MOI.ZeroOne,
    MOI.Semicontinuous{T},
    MOI.Semiinteger{T},
}

# 0xcb = 0x0080 | 0x0040 | 0x0008 | 0x0002 | 0x0001
const _LOWER_BOUND_MASK = 0x00cb
# 0xcd = 0x0080 | 0x0040 | 0x0008 | 0x0004 | 0x0001
const _UPPER_BOUND_MASK = 0x00cd

const _DELETED_VARIABLE = 0x8000

_single_variable_flag(::Type{<:MOI.EqualTo}) = 0x0001
_single_variable_flag(::Type{<:MOI.GreaterThan}) = 0x0002
_single_variable_flag(::Type{<:MOI.LessThan}) = 0x0004
_single_variable_flag(::Type{<:MOI.Interval}) = 0x0008
_single_variable_flag(::Type{MOI.Integer}) = 0x0010
_single_variable_flag(::Type{MOI.ZeroOne}) = 0x0020
_single_variable_flag(::Type{<:MOI.Semicontinuous}) = 0x0040
_single_variable_flag(::Type{<:MOI.Semiinteger}) = 0x0080

function _flag_to_set_type(flag::UInt16, ::Type{T}) where {T}
    if flag == 0x0001
        return MOI.EqualTo{T}
    elseif flag == 0x0002
        return MOI.GreaterThan{T}
    elseif flag == 0x0004
        return MOI.LessThan{T}
    elseif flag == 0x0008
        return MOI.Interval{T}
    elseif flag == 0x0010
        return MOI.Integer
    elseif flag == 0x0020
        return MOI.ZeroOne
    elseif flag == 0x0040
        return MOI.Semicontinuous{T}
    else
        @assert flag == 0x0080
        return MOI.Semiinteger{T}
    end
end

# Julia doesn't infer `S1` correctly, so we use a function barrier to improve
# inference.
function _throw_if_lower_bound_set_inner(variable, S2, mask, T)
    S1 = _flag_to_set_type(mask, T)
    throw(MOI.LowerBoundAlreadySet{S1,S2}(variable))
    return
end

function _throw_if_lower_bound_set(variable, S2, mask, T)
    lower_mask = mask & _LOWER_BOUND_MASK
    if iszero(lower_mask)
        return  # No lower bound set.
    elseif iszero(_single_variable_flag(S2) & _LOWER_BOUND_MASK)
        return  # S2 isn't related to the lower bound.
    end
    return _throw_if_lower_bound_set_inner(variable, S2, lower_mask, T)
end

# Julia doesn't infer `S1` correctly, so we use a function barrier to improve
# inference.
function _throw_if_upper_bound_set_inner(variable, S2, mask, T)
    S1 = _flag_to_set_type(mask, T)
    throw(MOI.UpperBoundAlreadySet{S1,S2}(variable))
    return
end

function _throw_if_upper_bound_set(variable, S2, mask, T)
    upper_mask = mask & _UPPER_BOUND_MASK
    if iszero(upper_mask)
        return  # No upper bound set.
    elseif iszero(_single_variable_flag(S2) & _UPPER_BOUND_MASK)
        return  # S2 isn't related to the upper bound.
    end
    return _throw_if_upper_bound_set_inner(variable, S2, upper_mask, T)
end

function _lower_bound(
    set::Union{MOI.GreaterThan,MOI.Interval,MOI.Semicontinuous,MOI.Semiinteger},
)
    return set.lower
end

_lower_bound(set::MOI.EqualTo) = set.value

function _upper_bound(
    set::Union{MOI.LessThan,MOI.Interval,MOI.Semicontinuous,MOI.Semiinteger},
)
    return set.upper
end

_upper_bound(set::MOI.EqualTo) = set.value

# Use `-Inf` and `Inf` for `AbstractFloat` subtypes.
_no_lower_bound(::Type{T}) where {T} = zero(T)
_no_lower_bound(::Type{T}) where {T<:AbstractFloat} = typemin(T)
_no_upper_bound(::Type{T}) where {T} = zero(T)
_no_upper_bound(::Type{T}) where {T<:AbstractFloat} = typemax(T)

###
### VariablesContainer
###
### For use in MOI.Utilities.Model
###

"""
    struct VariablesContainer{T} <: AbstractVectorBounds
        set_mask::Vector{UInt16}
        lower::Vector{T}
        upper::Vector{T}
    end

A struct for storing variables and VariableIndex-related constraints. Used in
`MOI.Utilities.Model` by default.
"""
mutable struct VariablesContainer{T} <: AbstractVectorBounds
    set_mask::Vector{UInt16}
    lower::Vector{T}
    upper::Vector{T}
end

function VariablesContainer{T}() where {T}
    return VariablesContainer{T}(UInt16[], T[], T[])
end

function Base.:(==)(a::VariablesContainer, b::VariablesContainer)
    return a.set_mask == b.set_mask && a.lower == b.lower && a.upper == b.upper
end

function MOI.empty!(b::VariablesContainer)
    empty!(b.set_mask)
    empty!(b.lower)
    empty!(b.upper)
    return b
end

function MOI.is_empty(b::VariablesContainer)
    if length(b.set_mask) == 0
        return true
    end
    return all(isequal(_DELETED_VARIABLE), b.set_mask)
end

function Base.resize!(b::VariablesContainer, n)
    resize!(b.set_mask, n)
    resize!(b.lower, n)
    resize!(b.upper, n)
    return
end

function MOI.add_variable(b::VariablesContainer{T}) where {T}
    push!(b.set_mask, 0x0000)
    push!(b.lower, _no_lower_bound(T))
    push!(b.upper, _no_upper_bound(T))
    x = MOI.VariableIndex(length(b.set_mask))
    return x
end

function MOI.get(b::VariablesContainer, ::MOI.ListOfVariableIndices)
    return MOI.VariableIndex[
        MOI.VariableIndex(i) for
        i in 1:length(b.set_mask) if b.set_mask[i] != _DELETED_VARIABLE
    ]
end

function MOI.is_valid(b::VariablesContainer, x::MOI.VariableIndex)
    mask = get(b.set_mask, x.value, _DELETED_VARIABLE)
    return mask != _DELETED_VARIABLE
end

function MOI.get(b::VariablesContainer, ::MOI.NumberOfVariables)::Int64
    if length(b.set_mask) == 0
        return 0
    end
    return sum(x != _DELETED_VARIABLE for x in b.set_mask)
end

function MOI.supports_constraint(
    ::VariablesContainer{T},
    ::Type{MOI.VariableIndex},
    ::Type{<:SUPPORTED_VARIABLE_SCALAR_SETS{T}},
) where {T}
    return true
end

function MOI.add_constraint(
    b::VariablesContainer{T},
    f::MOI.VariableIndex,
    set::S,
) where {T,S<:SUPPORTED_VARIABLE_SCALAR_SETS}
    flag = _single_variable_flag(S)
    mask = b.set_mask[f.value]
    _throw_if_lower_bound_set(f, S, mask, T)
    _throw_if_upper_bound_set(f, S, mask, T)
    if !iszero(flag & _LOWER_BOUND_MASK)
        b.lower[f.value] = _lower_bound(set)
    end
    if !iszero(flag & _UPPER_BOUND_MASK)
        b.upper[f.value] = _upper_bound(set)
    end
    b.set_mask[f.value] = mask | flag
    return MOI.ConstraintIndex{MOI.VariableIndex,S}(f.value)
end

function MOI.delete(
    b::VariablesContainer{T},
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {T,S}
    MOI.throw_if_not_valid(b, ci)
    flag = _single_variable_flag(S)
    b.set_mask[ci.value] &= ~flag
    if !iszero(flag & _LOWER_BOUND_MASK)
        b.lower[ci.value] = _no_lower_bound(T)
    end
    if !iszero(flag & _UPPER_BOUND_MASK)
        b.upper[ci.value] = _no_upper_bound(T)
    end
    return
end

function MOI.delete(b::VariablesContainer, x::MOI.VariableIndex)
    MOI.throw_if_not_valid(b, x)
    b.set_mask[x.value] = _DELETED_VARIABLE
    return
end

function MOI.is_valid(
    b::VariablesContainer,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S}
    if !(1 <= ci.value <= length(b.set_mask))
        return false
    end
    return !iszero(b.set_mask[ci.value] & _single_variable_flag(S))
end

function MOI.get(
    model::VariablesContainer,
    ::MOI.ConstraintFunction,
    ci::CI{MOI.VariableIndex},
)
    MOI.throw_if_not_valid(model, ci)
    return MOI.VariableIndex(ci.value)
end

function MOI.set(
    ::VariablesContainer,
    ::MOI.ConstraintFunction,
    ::MOI.ConstraintIndex{MOI.VariableIndex},
    ::MOI.VariableIndex,
)
    return throw(MOI.SettingVariableIndexNotAllowed())
end

function MOI.get(
    model::VariablesContainer,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S}
    MOI.throw_if_not_valid(model, ci)
    return set_from_constants(model, S, ci.value)
end

function MOI.set(
    model::VariablesContainer{T},
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
    set::S,
) where {T,S<:SUPPORTED_VARIABLE_SCALAR_SETS{T}}
    MOI.throw_if_not_valid(model, ci)
    flag = _single_variable_flag(S)
    if !iszero(flag & _LOWER_BOUND_MASK)
        model.lower[ci.value] = _lower_bound(set)
    end
    if !iszero(flag & _UPPER_BOUND_MASK)
        model.upper[ci.value] = _upper_bound(set)
    end
    return
end

function MOI.get(
    b::VariablesContainer,
    ::MOI.NumberOfConstraints{MOI.VariableIndex,S},
)::Int64 where {S}
    flag = _single_variable_flag(S)
    return count(mask -> !iszero(flag & mask), b.set_mask)
end

function _add_constraint_type(
    list,
    b::VariablesContainer,
    S::Type{<:MOI.AbstractScalarSet},
)
    flag = _single_variable_flag(S)::UInt16
    if any(mask -> !iszero(flag & mask), b.set_mask)
        push!(list, (MOI.VariableIndex, S))
    end
    return
end

function MOI.get(
    b::VariablesContainer{T},
    ::MOI.ListOfConstraintTypesPresent,
) where {T}
    list = Tuple{Type,Type}[]
    _add_constraint_type(list, b, MOI.EqualTo{T})
    _add_constraint_type(list, b, MOI.GreaterThan{T})
    _add_constraint_type(list, b, MOI.LessThan{T})
    _add_constraint_type(list, b, MOI.Interval{T})
    _add_constraint_type(list, b, MOI.Semicontinuous{T})
    _add_constraint_type(list, b, MOI.Semiinteger{T})
    _add_constraint_type(list, b, MOI.Integer)
    _add_constraint_type(list, b, MOI.ZeroOne)
    return list
end

function MOI.get(
    b::VariablesContainer,
    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,S},
) where {S}
    list = MOI.ConstraintIndex{MOI.VariableIndex,S}[]
    flag = _single_variable_flag(S)
    for (index, mask) in enumerate(b.set_mask)
        if !iszero(mask & flag)
            push!(list, MOI.ConstraintIndex{MOI.VariableIndex,S}(index))
        end
    end
    return list
end

# These `MOI.ModelLike` fallbacks have to be redefined here as
# `VariablesContainer` is not a subtype of `MOI.ModelLike`

function MOI.throw_if_not_valid(b::VariablesContainer, index)
    if !MOI.is_valid(b, index)
        throw(MOI.InvalidIndex(index))
    end
end

function MOI.set(
    ::VariablesContainer,
    ::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex,
    func,
)
    return throw(MOI.FunctionTypeMismatch{MOI.func_type(ci),typeof(func)}())
end

function MOI.set(
    ::VariablesContainer,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex},
    set,
)
    return throw(MOI.SetTypeMismatch{MOI.set_type(ci),typeof(set)}())
end

function MOI.supports_constraint(
    ::VariablesContainer,
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return false
end

function MOI.add_constraint(
    ::VariablesContainer,
    func::MOI.AbstractFunction,
    set::MOI.AbstractSet,
)
    return throw(MOI.UnsupportedConstraint{typeof(func),typeof(set)}())
end

###
### Hyperrectangle
###

"""
    struct Hyperrectangle{T} <: AbstractVectorBounds
        lower::Vector{T}
        upper::Vector{T}
    end

A struct for the .constants field in MatrixOfConstraints.
"""
struct Hyperrectangle{T} <: AbstractVectorBounds
    lower::Vector{T}
    upper::Vector{T}
end

Hyperrectangle{T}() where {T} = Hyperrectangle{T}(T[], T[])

function Base.:(==)(a::Hyperrectangle, b::Hyperrectangle)
    return a.lower == b.lower && a.upper == b.upper
end

function Base.empty!(b::Hyperrectangle)
    empty!(b.lower)
    empty!(b.upper)
    return b
end

function Base.resize!(b::Hyperrectangle, n)
    resize!(b.lower, n)
    resize!(b.upper, n)
    return
end

function load_constants(
    b::Hyperrectangle{T},
    offset,
    set::SUPPORTED_VARIABLE_SCALAR_SETS{T},
) where {T}
    flag = _single_variable_flag(typeof(set))
    if iszero(flag & _LOWER_BOUND_MASK)
        b.lower[offset+1] = _no_lower_bound(T)
    else
        b.lower[offset+1] = _lower_bound(set)
    end
    if iszero(flag & _UPPER_BOUND_MASK)
        b.upper[offset+1] = _no_upper_bound(T)
    else
        b.upper[offset+1] = _upper_bound(set)
    end
    return
end

function_constants(::Hyperrectangle{T}, row) where {T} = zero(T)


text:
# Usually identical to stan_inference in DiffEqBayes. Included here for testing.

struct StanModel{R,C}
  return_code::R
  chain_results::C
end

function generate_differential_equation(f)
  theta_ex = MacroTools.postwalk(f.pfuncs[1]) do x
    i = findfirst(f.params,x)
    i != 0 ? Symbol("theta[$i]") : x
  end
  differential_equation = ""
  for i in 1:length(theta_ex.args)-1
    differential_equation = string(differential_equation,theta_ex.args[i], ";\n")
  end
  return differential_equation
end

function generate_priors(f,priors)
  priors_string = ""
  params = f.params
  if priors==nothing
    for i in 1:length(params)
      priors_string = string(priors_string,"theta[$i] ~ normal(0, 1)", " ; ")
    end
  else
    for i in 1:length(params)
      μ = priors[i].μ
      σ = priors[i].σ
      priors_string = string(priors_string,"theta[$i] ~ normal($μ, $σ)", " ; ")
    end
  end
  priors_string
end

function stan_inference(prob::DEProblem,t,data,priors = nothing;alg=:rk45,
                            num_samples=1000, num_warmup=1000, reltol=1e-3,
                            abstol=1e-6, maxiter=Int(1e5),kwargs...)
  length_of_y = string(length(prob.u0))
  f = prob.f
  length_of_parameter = string(length(f.params))
  if alg ==:rk45
    algorithm = "integrate_ode_rk45"
  elseif alg == :bdf
    algorithm = "integrate_ode_bdf"
  else
    error("The choices for alg are :rk45 or :bdf")
  end

  differential_equation = generate_differential_equation(f)
  priors_string = generate_priors(f,priors)

  parameter_estimation_model = "
  functions {
    real[] sho(real t,real[] internal_var___u,real[] theta,real[] x_r,int[] x_i) {
      real internal_var___du[$length_of_y];
      $differential_equation
      return internal_var___du;
      }
    }
  data {
    real u0[$length_of_y];
    int<lower=1> T;
    real internal_var___u[T,$length_of_y];
    real t0;
    real ts[T];
  }
  transformed data {
    real x_r[0];
    int x_i[0];
  }
  parameters {
    vector<lower=0>[$length_of_y] sigma;
    real theta[$length_of_parameter];
  }
  model{
    real u_hat[T,$length_of_y];
    sigma ~ inv_gamma(2, 3);
    $priors_string
    u_hat = $algorithm(sho, u0, t0, ts, theta, x_r, x_i, $reltol, $abstol, $maxiter);
    for (t in 1:T){
      internal_var___u[t] ~ normal(u_hat[t], sigma);
      }
  }
  "

  stanmodel = Stanmodel(num_samples=num_samples, num_warmup=num_warmup, name="parameter_estimation_model", model=parameter_estimation_model);
  parameter_estimation_data = Dict("u0"=>prob.u0, "T" => size(t)[1], "internal_var___u" => data', "t0" => prob.tspan[1], "ts" => t)
  return_code, chain_results = stan(stanmodel, [parameter_estimation_data]; CmdStanDir=CMDSTAN_HOME)
  return StanModel(return_code,chain_results)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule VideoCaptureWrap_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("VideoCaptureWrap")
JLLWrappers.@generate_main_file("VideoCaptureWrap", UUID("b69a920f-0f27-5c26-a11f-67bef8892cc7"))
end  # module VideoCaptureWrap_jll


text:
using Fermi.Orbitals
using TensorOperations

function boys_localize(orbs::AbstractRestrictedOrbitals)
    nothing
end

function ER_localize(orbs::AbstractRestrictedOrbitals, ints, mix_range)

    eri = ints["ERI"]

    # Maximaze the self-repulsion (pp|pp)
    newC = deepcopy(orbs.C)

    N = length(mix_range)

    println("Number of mixing pairs: $N")
    # Create a list of pairs to be rotated
    Npairs = N*(N-1) ÷ 2
    Δvals = zeros(Npairs)
    θvals = zeros(Npairs)
    ij_vals = Array{Tuple{Int64, Int64}}(undef, Npairs)

    ## Fill ij_vals
    k = 1
    for i in mix_range
        for j in mix_range
            if i > j
                ij_vals[k] = (i,j)
                k += 1
            end
        end
    end

    #@views Co = newC[:, o]
    #self = zeros(No)
    #J = zeros(No, No)
    #K = zeros(No, No)
    #L = zeros(No, No)
    #for p = o
    #    @views C1 = newC[:,p]
    #    @tensoropt A = C1[μ]*C1[ν]*C1[λ]*C1[σ]*eri[μ, ν, λ, σ]
    #    self[p] = A
    #    for q = o
    #        @views C2 = newC[:,q]
    #        @tensoropt B = C1[μ]*C1[ν]*C2[λ]*C2[σ]*eri[μ, ν, λ, σ]
    #        @tensoropt C = C1[μ]*C2[ν]*C1[λ]*C2[σ]*eri[μ, ν, λ, σ]
    #        @tensoropt D = C1[μ]*C1[ν]*C1[λ]*C2[σ]*eri[μ, ν, λ, σ]

    #        J[p,q] = B
    #        K[p,q] = C
    #        L[p,q] = D
    #    end
    #end

    println("Self repulsion: $(selfrepulsion(newC, mix_range, ints))")

    ite = 1
    for idx in eachindex(ij_vals)
        (i,j) = ij_vals[idx]

        #iiii = self[i]
        #jjjj = self[j]

        #iiij = L[i,j]
        #jjij = L[j,i]

        #iijj = J[i,j]
        #ijij = K[i,j]

        #B = iiij - jjij
        #A = ijij + 0.5*iijj - 0.25*(iiii + jjjj)
        #θ = 0.25 * atan(-B / A)

        #Δvals[idx] = A*(1 - cos(4θ)) + B*sin(4θ)
        #θvals[idx] = θ

        θvals[idx], Δvals[idx] = ER_rotation_angle(newC, i,j, ints["ERI"])
    end

    while true
        if ite > 30
            println("Ops did not converge")
            break
        end

        println("Iter $ite")

        RMS = √(sum((Δvals) .^2) / length(Δvals))

        println("ΔI RMS: $RMS")

        if RMS < 1e-5
            break
        end

        _, idx = findmax(Δvals)

        i,j = ij_vals[idx]
        θ = θvals[idx]
        Δ = Δvals[idx]

        if Δ < 1e-8
            println("Converged")
            break
        end

        println("Rotating $i and $j by $θ. Δ expected: $Δ")
        rotate_orbitals!(newC, θ, i,j)
        θvals[idx], Δvals[idx] = ER_rotation_angle(newC, i,j, ints["ERI"])

        #@assert abs(Θvals[idx]) < 1e-10
        #@assert abs(Δvals[idx]) < 1e-10

        println("Updating integrals")

        for k in mix_range
            # Update i vals
            if k > i
                idx, = findall(x->x==(k,i), ij_vals)  # PROVISORIO
                θvals[idx], Δvals[idx] = ER_rotation_angle(newC, k,i, ints["ERI"])
            elseif i > k
                idx, = findall(x->x==(i,k), ij_vals)  # PROVISORIO
                θvals[idx], Δvals[idx] = ER_rotation_angle(newC, i,k, ints["ERI"])
            end

            # Update j vals
            if k > j
                idx, = findall(x->x==(k,j), ij_vals)  # PROVISORIO
                θvals[idx], Δvals[idx] = ER_rotation_angle(newC, k,j, ints["ERI"])
            elseif j > k
                idx, = findall(x->x==(j,k), ij_vals)  # PROVISORIO
                θvals[idx], Δvals[idx] = ER_rotation_angle(newC, j,k, ints["ERI"])
            end
        end

        println("Itotal after rotation: $(selfrepulsion(newC, mix_range, ints))")
        ite += 1
    end

    energy = sd_energy(orbs, orbs.molecule.Nα, ints)
    neworbs = Fermi.Orbitals.GeneralRestrictedOrbitals(orbs.molecule, orbs.basis, energy, newC)
    return neworbs
end

function selfrepulsion(C, r, ints)
    sum([get_pppp(C, ints["ERI"], i) for i = r])
end

function ER_rotation_angle(C, i,j, eri)

    @views Ci = C[:,i]
    @views Cj = C[:,j]

    iiii = @tensoropt Ci[μ]*Ci[ν]*Ci[ρ]*Ci[σ]*eri[μ,ν,ρ,σ]
    jjjj = @tensoropt Cj[μ]*Cj[ν]*Cj[ρ]*Cj[σ]*eri[μ,ν,ρ,σ]

    iiij = @tensoropt Ci[μ]*Ci[ν]*Ci[λ]*Cj[σ]*eri[μ,ν,λ,σ]
    jjij = @tensoropt Cj[μ]*Cj[ν]*Ci[λ]*Cj[σ]*eri[μ,ν,λ,σ]

    iijj = @tensoropt Ci[μ]*Ci[ν]*Cj[λ]*Cj[σ]*eri[μ,ν,λ,σ]
    ijij = @tensoropt Ci[μ]*Cj[ν]*Ci[λ]*Cj[σ]*eri[μ,ν,λ,σ]

    B = iiij - jjij
    A = ijij + 0.5*iijj - 0.25*(iiii + jjjj)
    θ = 0.25 * atan(-B / A)

    Δ = A*(1 - cos(4θ)) + B*sin(4θ)

    return θ, Δ
end

function rotate_orbitals!(C, θ, i, j)
    ui = C[:,i]
    uj = C[:,j]

    C[:,i] =  ui .* cos(θ) + uj .* sin(θ)
    C[:,j] = -ui .* sin(θ) + uj .* cos(θ)
end

# Get mo (pp|pp) integral
function get_pppp(C, aoeri, p)
    @views Cp = C[:,p]
    @tensoropt I = Cp[μ]*Cp[ν]*Cp[ρ]*Cp[σ]*aoeri[μ,ν,ρ,σ]
    return I
end

# Get mo (pp|pq) integral
function get_pppq(C, aoeri, p, q)
    @views Cp = C[:,p]
    @views Cq = C[:,q]
    @tensoropt I = Cp[μ]*Cp[ν]*Cp[ρ]*Cq[σ]*aoeri[μ,ν,ρ,σ]
    return I
end

# Get mo (pp|qq) integral
function get_ppqq(C, aoeri, p, q)
    @views Cp = C[:,p]
    @views Cq = C[:,q]
    @tensoropt I = Cp[μ]*Cp[ν]*Cq[ρ]*Cq[σ]*aoeri[μ,ν,ρ,σ]
    return I
end

# Get mo (pq|pq) integral
function get_pqpq(C, aoeri, p, q)
    @views Cp = C[:,p]
    @views Cq = C[:,q]
    @tensoropt I = Cp[μ]*Cq[ν]*Cp[ρ]*Cq[σ]*aoeri[μ,ν,ρ,σ]
    return I
end

function sd_energy(orbs, ndocc, ints)

    C = orbs.C
    Co = C[:,1:ndocc]
    D = Co * Co'

    H = ints["T"] + ints["V"]
    F = similar(H)

    Fermi.HartreeFock.build_fock!(F, H, D, ints)

    E = Fermi.HartreeFock.RHFEnergy(D, H, F)

    return E + orbs.molecule.Vnuc
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
# Parameters

beta = .96
y = [1.0, 2.0]
b0 = 0.0
P = [0.8 0.2;
     0.4 0.6]

cp = ConsumptionProblem(beta, y, b0, P)
Q = beta*P
N_simul = 150

c_bar, b1, b2 = consumption_complete(cp)
debt_complete = [b1, b2]

println("P = ", P)
println("Q= ", Q, "\n")
println("Govt expenditures in peace and war =", y)
println("Constant tax collections = ", c_bar)
println("Govt assets in two states = ", debt_complete)

msg = """
Now let's check the government's budget constraint in peace and war.
Our assumptions imply that the government always purchases 0 units of the
Arrow peace security.
"""
println(msg)

AS1 = Q[1,2] * b2
println("Spending on Arrow war security in peace = ", AS1)
AS2 = Q[2,2]*b2
println("Spending on Arrow war security in war = ", AS2)

println("\n")
println("Government tax collections plus asset levels in peace and war")
TB1=c_bar+b1
println("T+b in peace = ",TB1 )
TB2 = c_bar + b2
println("T+b in war = ", TB2)

println("\n")
println("Total government spending in peace and war")
G1= y[1] + AS1
G2 = y[2] + AS2
println("total govt spending in peace = ", G1)
println("total govt spending in war = ", G2)

println("\n")
println("Let's see ex post and ex ante returns on Arrow securities")

Pi= 1./Q#reciprocal(Q)
exret= Pi
println("Ex post returns to purchase of Arrow securities = $exret")
exant = Pi.*P
println("Ex ante returns to purchase of Arrow securities = $exant")


text:
module FundamentalSolutions
#
# This module calculates:
#
# Φ(x,y,E), the fundamental solution of the Helmholtz equation in a linearly stratified medium:
#
# -(Δ + E + x_2)Φ(x,y,E) = δ(x-y)
#
# also known as the gravity Helmholtz equation.
#
# This is a Julia translation of the C-library described in:
#
# A. H. Barnett, B. J. Nelson, and J. M. Mahoney, High-order boundary integral equation solution of high frequency wave scattering from obstacles in an unbounded linearly stratified medium, J. Comp. Phys., 297:407–426, 2015.
#
# ℜ(x,y,E), the Riemann function of the Helmholtz equation in a linearly stratified medium. The bivariate series is described in:
#
# R. M. Slevinsky and S. Olver, A fast and well-conditioned spectral method for singular integral equations, arXiv:1507.00596, 2015.
#

export lhelmfs, lhelm_riemann

include("constants.jl")
include("integrand.jl")
include("gammaforbidden.jl")
include("contour.jl")
include("locate_minimum.jl")
include("find_endpoints.jl")
include("quad_nodes.jl")
include("lhfs.jl")

include("lhelmfs.jl")
include("lhelm_riemann.jl")

end #module

using .FundamentalSolutions

export lhelmfs, lhelm_riemann


text:
"""
    maxcurv(mags, ΔM)

Apply the maximum curvature method to compute the magnitude of completeness.

# Arguments
* `mags`: array storing the magnitude values of an incomplete catalog.
* `ΔM`: bin width to be used to group the magnitude values.`
"""
function maxcurv(mags::AbstractVector{T},ΔM::S) where {T,S<:Real}

	bincenters, binevents = histmag(mags, ΔM)

	Mc = maxcurv(bincenters, binevents)

	return Mc
end

"""
    maxcurv(bincenters, binevents)

Apply the maximum curvature method to compute the magnitude of completeness.

# Arguments
* `bincenters`: array storing the central values of the magnitude bins.
* `binevents`: array storing the number of events in the magnitude bins.

`bincenters` and `binevents` can be obtained with  [`histmag(mag, ΔM)`](@ref).
"""
function maxcurv(bincenters::AbstractVector{T},binevents::AbstractVector{S}) where {T,S<:Real}

    if length(bincenters) != length(binevents)
        error("The two input arrays must have the same number of elements")
    end

    Mc::Float64 = bincenters[ argmax(binevents) ]

    return Mc
end


text:
# Tests for Centre and LazyCentre types, and associated methods

@testset "Centre" begin
    @testset "Constructors" begin
        @test_throws DomainError Centre(-1)

        a = Centre([1,1,2,3,6])
        b = Centre(a)
        c = Centre{Integer}(4)
        d = Centre(4)
        e = Centre(LazyCentre(d))

        @test a == b == c == d == e

        b.data[2] = 4
        @test a != b

        a = Centre([1.0, 1.0, 2.0, 3.0, 6.0])
        b = Centre(a)
        c = Centre{Float64}(4)
        d = Centre(4) # different data type
        @test a == b == c == d
    end

    @testset "Accessors" begin
        a = Centre(4)
        @test values(a) == [1, 1, 2, 3, 6]
    end

    @testset "Array functions" begin
        a = Centre(8)
        b = Centre{Float64}(8)
        @test axes(a) == axes(b) == (PT.ZeroRange(8),)

        @testset "Data access $i" for i ∈ 0:8
            @test a[i] == b[i]
        end

        @test_throws BoundsError a[-1]
        @test_throws BoundsError a[10]

        @test firstindex(a) == firstindex(b) == 0
        @test lastindex(a) == lastindex(b) == 8
    end

    @testset "Checks" begin
        a = Centre(4)
        b = Centre([1, 1, 2, 3, 6])
        c = Centre([1, 5, 15, 35, 80])

        @test PT.isvalid(a)
        @test !PT.isvalid(c)

        @test toarray(a) == toarray(b)
    end
end

@testset "LazyCentre" begin
    @testset "Constructors" begin
        a = LazyCentre(Dict(0 => 1, 1 => 1))
        b = LazyCentre(a)
        c = LazyCentre{Integer}()
        d = LazyCentre()
        e = LazyCentre(Centre(1))
        @test a == b == c == d == e

        b.data[1] = 2
        @test a != b

        a = LazyCentre(Dict(0 => 1.0, 1 => 1.0))
        b = LazyCentre(a)
        c = LazyCentre{Float64}()
        d = LazyCentre() # different data type
        e = LazyCentre(Centre{Float64}(1))
        @test a == b == c == d == e
    end

    @testset "Array functions" begin
        a = LazyCentre()
        b = LazyCentre{Float64}()

        @testset "Data access $i" for i ∈ 0:8
            @test a[i] == b[i]
        end

        @test_throws BoundsError a[-1]

        @test firstindex(a) == firstindex(b) == 0
    end

    @testset "Checks" begin
        a = LazyCentre(Centre(4))
        b = LazyCentre(
                Dict(0 => 1,
                    1 => 1,
                    2 => 2,
                    3 => 3,
                    4 => 6))
        c = LazyCentre(Dict(0 => 0, 1 => 0))

        @test PT.isvalid(a)
        @test !PT.isvalid(c)

        @test toarray(a) == toarray(b)
    end
end


text:
"Rotates vector counter-clockwise with respecto to the x-axis."
rotx(θ::Real) = [1		0		0
				 0	cos(θ)	-sin(θ);
				 0	sin(θ)	cos(θ)]
"Rotates vector counter-clockwise with respecto to the y-axis."
roty(θ::Real) = [cos(θ) 0	sin(θ);
				 0		1		0;
				-sin(θ)	0 	cos(θ)]
"Rotates vector counter-clockwise with respecto to the z-axis."
rotz(θ::Real) = [cos(θ) -sin(θ)	0;
				 sin(θ) cos(θ)	0;
				 0 		0 		1]
#####################
## Gradient OBJECT ##
#####################
# @everywhere begin
"""
    Grad(A,T)

Gradient object.

# Arguments
- `A::Real` := Gradient amplitude [T].
- `T::Real` := Gradient duration [s].

# Examples
```@example
Grad(1,2)
```
"""
mutable struct Grad
	A #Amplitud [T], if this is a scalar it makes a trapezoid
	T     #Duration of flat-top [s]
	rise::Real  #Duration of rise [s]
	fall::Real  #Duration of fall [s]
	delay::Real #Duration of fall [s]
    function Grad(A,T,rise,fall,delay)
		all(T .< 0) || rise < 0 || fall < 0 || delay < 0 ? error("Gradient timings must be positive.") : new(A, T, rise, fall, delay)
    end
	function Grad(A,T,rise,delay)
		all(T .< 0) < 0 || rise < 0 || delay < 0 ? error("Gradient timings must be positive.") : new(A, T, rise, rise, delay)
    end
	function Grad(A,T,rise)
		all(T .< 0) < 0 || rise < 0 ? error("Gradient timings must be positive.") : new(A, T, rise, rise, 0)
    end
	function Grad(A,T)
		all(T .< 0) < 0 ? error("Gradient timings must be positive.") : new(A, T, 0, 0, 0)
    end
end
#Gradient operations
*(x::Grad,α::Real) = Grad(α*x.A,x.T,x.rise,x.fall,x.delay)
*(α::Real,x::Grad) = Grad(α*x.A,x.T,x.rise,x.fall,x.delay)
*(A::Matrix{Float64},GR::Matrix{Grad}) = begin
	N, M = size(GR)
	[sum(A[i,1:N] .* GR[:,j]) for i=1:N, j=1:M]
end
zero(::Grad) = Grad(0,0)
/(x::Grad,α::Real) = Grad(x.A/α,x.T,x.rise,x.fall,x.delay)
+(x::Grad,y::Grad) = Grad(x.A+y.A,x.T,x.rise,x.fall,x.delay)
+(x::Array{Grad,1}, y::Array{Grad,1}) = [x[i]+y[i] for i=1:length(x)]
-(x::Grad) = -1*x
-(x::Grad,y::Grad) = Grad(x.A-y.A,x.T,x.rise,x.fall,x.delay)
#Gradient functions
vcat(x::Array{Grad,1},y::Array{Grad,1}) = [i==1 ? x[j] : y[j] for i=1:2,j=1:length(x)]
vcat(x::Array{Grad,1},y::Array{Grad,1},z::Array{Grad,1}) = [i==1 ? x[j] : i==2 ? y[j] : z[j] for i=1:3,j=1:length(x)]
getproperty(x::Vector{Grad}, f::Symbol) = getproperty.(x,f)
getproperty(x::Matrix{Grad}, f::Symbol) = begin
	if f == :x
		x[1,:]
	elseif f == :y && size(x,1) >= 2
		x[2,:]
	elseif f == :z && size(x,1) >= 3
		x[3,:]
	elseif f == :T || f == :rise || f == :fall || f == :delay
		getproperty.(x,f)
	elseif f == :dur
		T, ζ1, ζ2, delay = x.T, x.rise, x.fall, x.delay
		ΔT = [sum(t) for t=T] .+ ζ1 .+ ζ2 .+ delay
		maximum(ΔT,dims=1)[:]
	else
		getproperty.(x,f)
	end
end

#TIMINGS
"Duration `T` [s] of Grad datatype."
dur(x::Grad) = x.delay + x.rise + sum(x.T) + x.fall
"Duration `T` [s] of Vector{Grad}."
dur(x::Vector{Grad}) = maximum(dur.(x), dims=1)[:]

"""
	Grad_fun(f,T,N)

Generates **arbitrary gradient waveform** defined by function `f` in the interval t∈[0,`T`].
It uses `N` square gradients uniformly spaced in the interval.

# Arguments
 - `f::Function` := Gradient waveform.
 - `T::Real`     := Gradient duration.
 - `N::Integer=1`:= Number of sample points.

# Examples
```julia-repl
julia> Grad_fun(x-> sin(π*x),1,4)
1×4 Array{Grad,2}:
 Grad(0.0, 0.333333)  Grad(0.866025, 0.333333)  Grad(0.866025, 0.333333)  Grad(1.22465e-16, 0.333333)
```
"""
Grad(f::Function,T::Real,N::Int64=300) = begin
	t = range(0,T,N)
	G = f.(t')
	Grad(G,T)
end
#aux
Base.show(io::IO,x::Grad) = begin
	r(x) = round.(x,digits=4)
	compact = get(io, :compact, false)
	if !compact
		wave = length(x.A) == 1 ? r(x.A*1e3) : "∿"
		if x.rise == x.fall == 0.
			print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"Grad($(wave) mT, $(r(x.T*1e3)) ms)")
		else
			print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"Grad($(wave) mT, $(r(x.T*1e3)) ms, ↑$(r(x.rise*1e3)) ms, ↓$(r(x.fall*1e3)) ms)")
		end
	else
		wave = length(x.A) == 1 ? "⊓" : "∿"
		print(io, (sum(abs.(x.A)) > 0 ? wave : "⇿")*"($(r((x.delay+x.rise+x.fall+x.T)*1e3)) ms)")
	end
end


text:
# # Perfect convective adjustment calibration with Ensemble Kalman Inversion
#
# This example calibrates a convective adjustment model in the "perfect model context".
# In this context, synthetic observations are generated by a convective adjustment model
# with "true" parameters. The true parameters are then "rediscovered" by calibrating the model
# to match the synthetic observations.
#
# We use the discrepancy between observed and modeled buoyancy ``b`` to calibrate
# the convective adjustment model.
# The calibration problem is solved by Ensemble Kalman Inversion.
# For more information about Ensemble Kalman Inversion, see the
# [EnsembleKalmanProcesses.jl documentation](https://clima.github.io/EnsembleKalmanProcesses.jl/stable/ensemble_kalman_inversion/).

# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add ParameterEstimocean, Oceananigans, Distributions, CairoMakie"
# ```

using ParameterEstimocean, LinearAlgebra, CairoMakie

# We reuse some code from a previous example to generate observations,
examples_path = joinpath(pathof(ParameterEstimocean), "..", "..", "examples")
include(joinpath(examples_path, "intro_to_inverse_problems.jl"))

data_path = generate_synthetic_observations()
observations = SyntheticObservations(data_path, field_names=:b, transformation=ZScore())

# and an ensemble simulation,

ensemble_simulation, closure★ = build_ensemble_simulation(observations; Nensemble=50)

# The handy utility function `build_ensemble_simulation` also tells us the optimal
# parameters that were used when generating the synthetic observations:

@show θ★ = (convective_κz = closure★.convective_κz, background_κz = closure★.background_κz)

# # The `InverseProblem`
#
# To build an inverse problem we first define free parameters.
# Here we calibrate `convective_κz` and `background_κz`, using
# log-normal priors to prevent the parameters from becoming negative:

priors = (convective_κz = lognormal(mean=0.3, std=0.5),
          background_κz = lognormal(mean=2.5e-4, std=2.5e-5))

free_parameters = FreeParameters(priors)

# The `InverseProblem` is then constructed from `observations`, `ensemble_simulation`, and
# `free_parameters`,

calibration = InverseProblem(observations, ensemble_simulation, free_parameters)

# For more information about the above steps, see [Intro to observations](@ref)
# and [Intro to `InverseProblem`](@ref).

# # Ensemble Kalman Inversion
#
# Next, we construct an `EnsembleKalmanInversion` (EKI) object,
#
# The calibration is done here using Ensemble Kalman Inversion. For more information about the
# algorithm refer to
# [EnsembleKalmanProcesses.jl documentation](https://clima.github.io/EnsembleKalmanProcesses.jl/stable/ensemble_kalman_inversion/).

eki = EnsembleKalmanInversion(calibration; pseudo_stepping = ConstantConvergence(0.5))

# and perform few iterations to see if we can converge to the true parameter values.

iterate!(eki; iterations = 10)

# Last, we visualize the outputs of EKI calibration.

θ̅(iteration) = [eki.iteration_summaries[iteration].ensemble_mean...]
varθ(iteration) = eki.iteration_summaries[iteration].ensemble_var

weight_distances = [norm(θ̅(iter) - [θ★[1], θ★[2]]) for iter in 0:eki.iteration]
output_distances = [norm(forward_map(calibration, θ̅(iter))[:, 1] - y) for iter in 0:eki.iteration]
ensemble_variances = [varθ(iter) for iter in 0:eki.iteration]

f = Figure()

lines(f[1, 1], 0:eki.iteration, weight_distances, color = :red, linewidth = 2,
      axis = (title = "Parameter distance",
              xlabel = "Iteration",
              ylabel = "|θ̅ₙ - θ★|"))

lines(f[1, 2], 0:eki.iteration, output_distances, color = :blue, linewidth = 2,
      axis = (title = "Output distance",
              xlabel = "Iteration",
              ylabel = "|G(θ̅ₙ) - y|"))

ax3 = Axis(f[2, 1:2],
           title = "Parameter convergence",
           xlabel = "Iteration",
           ylabel = "Ensemble variance",
           yscale = log10)

for (i, pname) in enumerate(free_parameters.names)
    ev = getindex.(ensemble_variances, i)
    lines!(ax3, 0:eki.iteration, ev / ev[1], label = String(pname), linewidth = 2)
end

axislegend(ax3, valign = :top, halign = :right)

save("summary_convective_adjustment_eki.svg", f); nothing #hide

# ![](summary_convective_adjustment_eki.svg)

# And also we plot the the distributions of the various model ensembles for few EKI iterations to see
# if and how well they converge to the true diffusivity values.

fig = Figure()

axtop = Axis(fig[1, 1])
axmain = Axis(fig[2, 1], xlabel = "convective_κz [m² s⁻¹]",
                       ylabel = "background_κz [m² s⁻¹]")

axright = Axis(fig[2, 2])
scatters = []
labels = String[]

for iteration in [0, 1, 2, 10]
    ## Make parameter matrix
    parameters = eki.iteration_summaries[iteration].parameters
    Nensemble = length(parameters)
    Nparameters = length(first(parameters))
    parameter_ensemble_matrix = [parameters[i][j] for i=1:Nensemble, j=1:Nparameters]

    label = iteration == 0 ? "Initial ensemble" : "Iteration $iteration"
    push!(labels, label)
    push!(scatters, scatter!(axmain, parameter_ensemble_matrix))
    density!(axtop, parameter_ensemble_matrix[:, 1])
    density!(axright, parameter_ensemble_matrix[:, 2], direction = :y)
end

vlines!(axmain, [θ★.convective_κz], color = :red)
vlines!(axtop, [θ★.convective_κz], color = :red)

hlines!(axmain, [θ★.background_κz], color = :red)
hlines!(axright, [θ★.background_κz], color = :red)

colsize!(fig.layout, 1, Fixed(300))
colsize!(fig.layout, 2, Fixed(200))
rowsize!(fig.layout, 1, Fixed(200))
rowsize!(fig.layout, 2, Fixed(300))

Legend(fig[1, 2], scatters, labels, valign = :bottom, halign = :left)

hidedecorations!(axtop, grid = false)
hidedecorations!(axright, grid = false)

xlims!(axmain, -0.25, 3.2)
xlims!(axtop, -0.25, 3.2)
ylims!(axmain, 5e-5, 35e-5)
ylims!(axright, 5e-5, 35e-5)

save("distributions_convective_adjustment_eki.svg", fig); nothing #hide

# ![](distributions_convective_adjustment_eki.svg)


text:
mutable struct IISData
    stat::Int
    rownumber::Int # number of rows participating in the IIS
    colnumber::Int # number of columns participating in the IIS
    miisrow::Vector{Cint} # index of the rows that participate
    miiscol::Vector{Cint} # index of the columns that participate
    constrainttype::Vector{UInt8} # sense of the rows that participate
    colbndtype::Vector{UInt8} # sense of the column bounds that participate
end

function getfirstiis(model::Model)
    # This function always calls iisfirst first, which computes a first IIS.
    # In other words, any call to this function is expensive.

    # First, compute the first IIS.
    # int XPRS_CC XPRSiisfirst(XPRSprob prob, int ifiis, int *status_code);
    status = Ref{Cint}()
    ret = @xprs_ccall(iisfirst, Cint, (Ptr{Nothing}, Cint, Ptr{Cint}),
        model.ptr_model, 1, status)
    if ret != 0
        throw(XpressError(model))
    end

    # If the computation did not succeed, stop right now.
    if status[] == 1
        # The problem is actually feasible.
        return IISData(status[], 0, 0, Vector{Cint}(), Vector{Cint}(), Vector{UInt8}(), Vector{UInt8}())
    elseif status[] == 2
        # There was a problem in the computation; this should never happen, as
        # in this case the function should return a nonzero code.
        throw(XpressError(model))
    end

    # Then, retrieve it.
    # int XPRS_CC XPRSgetiisdata(XPRSprob prob, int num, int *rownumber, int *colnumber, int miisrow[],
    # int miiscol[], char constrainttype[], char colbndtype[], double duals[], double rdcs[],
    # char isolationrows[], char isolationcols[]);

    # XPRESS' API works in two steps: first, retrieve the sizes of the arrays to
    # retrieve; then, the user is expected to allocate the needed memory,
    # before asking XPRESS to fill it. (Making one call would require having large
    # enough arrays, which would waste lots of memory.)
    rownumber = Ref{Cint}()
    colnumber = Ref{Cint}()
    ret = @xprs_ccall(
        getiisdata,
        Cint,
        (Ptr{Nothing}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt8}, Ptr{UInt8}),
        model.ptr_model, 1, rownumber, colnumber, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)
    if ret != 0
        throw(XpressError(model))
    end

    miisrow = Vector{Cint}(undef, rownumber[])
    miiscol = Vector{Cint}(undef, colnumber[])
    constrainttype = Vector{UInt8}(undef, colnumber[])
    colbndtype = Vector{UInt8}(undef, colnumber[])
    # duals ignored (dual multiplier for each IIS constraint)
    # rdcs ignored (reduced costs for each IIS variable)
    # isolationrows ignored (would also required running iisisolations)
    # isolationcols ignored (would also required running iisisolations)
    ret = @xprs_ccall(
        getiisdata,
        Cint,
        (Ptr{Nothing}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt8}, Ptr{UInt8}),
        model.ptr_model, 1, rownumber, colnumber, miisrow, miiscol, constrainttype, colbndtype, C_NULL, C_NULL, C_NULL, C_NULL)
    if ret != 0
        throw(XpressError(model))
    end

    return IISData(status[], rownumber[], colnumber[], miisrow, miiscol, constrainttype, colbndtype)
end


text:
#fin. gen. submodules of Z^n and F_p^n (and possibly more)
#
#import Base.show, Base.inv, Nemo.inv, Base.solve, Hecke.solve,
#       Hecke.lift, Hecke.rational_reconstruction, Hecke.elementary_divisors,
#       Hecke.rank, Hecke.det

export det_mc, id, isid, isupper_triangular, norm2, hadamard_bound2,
       hnf, hnf!, echelon_with_trafo

const p = next_prime(2^20)

add_verbose_scope(:HNF)

add_assert_scope(:HNF)
set_assert_level(:HNF, 0)

function show(io::IO, M::ModuleCtxNmod)
  println("Sparse module over $(M.R) of (current) rank $(rows(M.basis)) and $(rows(M.gens))")
end

function show(io::IO, M::ModuleCtx_fmpz)
  println("Sparse module over FlintZZ of (current) rank $(rows(M.bas_gens)) and further $(rows(M.rel_gens))")
  if isdefined(M, :basis_idx)
    println("current index: $(M.basis_idx)")
  end
  if isdefined(M, :essential_elementary_divisors)
    println("current structure: $(M.essential_elementary_divisors)")
  end
end

function add_gen!(M::ModuleCtxNmod, g::SRow{nmod})
  if M.basis.r == M.basis.c
    return false
  end
  h = reduce(M.basis, g)
  if !iszero(h)
    i = 1
    while i<= rows(M.basis) && M.basis.rows[i].pos[1] < h.pos[1]
      i += 1
    end
    @hassert :HNF 2  i > rows(M.basis) || M.basis[i].pos[1] > h.pos[1]
    insert!(M.basis.rows, i, h)
    M.basis.r += 1
    M.basis.nnz += length(h)
    M.basis.c = max(M.basis.c, h.pos[end])
    push!(M.gens, g)
    return true
  end
  return false
end

function add_gen!(M::ModuleCtx_fmpz, g::SRow{fmpz})
  gp = SRow(g, M.Mp.R)
  M.new = true
  if add_gen!(M.Mp, gp)
    push!(M.bas_gens, g)
    return true
  else
    push!(M.rel_gens, g)
  end
  return false
end

function check_index(M::ModuleCtx_fmpz)
  if rows(M.Mp.basis) < cols(M.Mp.basis)
    return fmpz(0)
  end

  if !M.new
    return M.basis_idx
  end

  if isdefined(M, :trafo)  #once we have trafo, we need to keep it!
    module_trafo_assure(M)
    return M.basis_idx
  end

  M.new = false

  if isdefined(M, :basis)
    C = copy(M.basis)
  else
    d = abs(det(M.bas_gens))
    C = M.max_indep
    C.c = M.bas_gens.c
    for ii = M.bas_gens
      h = reduce(C, ii, 2*d) #to avoid problems with diag being 1...1 d
      @hassert :HNF 2  !iszero(h)
      i = 1
      while i<= rows(C) && C.rows[i].pos[1] < h.pos[1]
        i += 1
      end
      @hassert :HNF 2  i > rows(C) || C[i].pos[1] > h.pos[1]
      insert!(C.rows, i, h)
      C.r += 1
      C.nnz += length(h)
      C.c = max(C.c, h.pos[end])
    end
    M.max_indep = copy(C)
  end

  for i=length(M.rel_reps_p)+1:length(M.rel_gens)
    reduce(C, M.rel_gens[i])
    push!(M.rel_reps_p, solve_ut(M.Mp.basis, SRow(M.rel_gens[i], M.Mp.R)))
  end

#=
  for l=1:5
    mis = find(i->C[i,i] != 1, 1:rows(C))
    if length(mis) == 0
      break
    end
#    println("mis: $mis")
    for i = mis
      if C[i,i] == 1
#        println("miracle for $i")
        continue
      end
      r = find(x->i in M.rel_reps_p[x].pos, 1:length(M.rel_reps_p))
#      println("found $(length(r)) rows")
      if length(r) == 0
        break
      end
      g = M.rel_gens[rand(r)]
      for j=1:min(5, div(length(r), 2))
        g += rand(-10:10)*M.rel_gens[rand(r)]
      end
      reduce(C, g)
      if C[i,i] == 1
#        println("bingo for i=$i")
      end
    end
  end
=#
  M.basis = C
  M.basis_idx = prod([C[i,i] for i=1:rows(C)])

  return M.basis_idx
end

function elementary_divisors(M::ModuleCtx_fmpz)
  if !isdefined(M, :basis)
    i = check_index(M)
    if i == 0
      return fmpz[]
    end
  end
  C = M.basis
  f = find(i -> C[i,i] != 1, 1:rows(C))
  if length(f) == 0
    M.essential_elementary_divisors = fmpz[]
    return M.essential_elementary_divisors
  end
  e = minimum(f)
  m = fmpz_mat(sub(C, e:rows(C), e:cols(C)))
  s = snf(m)
  M.essential_elementary_divisors = [s[i,i] for i=1:rows(s)]
  return M.essential_elementary_divisors
end

function missing_pivot(M::ModuleCtx_fmpz)
  C = M.Mp.basis
  return setdiff(IntSet(1:cols(C)), [x.pos[1] for x=C])
end

function non_trivial_pivot(M::ModuleCtx_fmpz)
  h = check_index(M)
  if h == 0
    return missing_pivot(M)
  end
  C = M.basis
  @hassert :HNF 2  C.r == C.c
  return setdiff(IntSet(1:cols(C)), find(i->C[i].values[1] == 1, 1:C.c))
end

function rank(M::ModuleCtx_fmpz)
  return M.bas_gens.r
end

function rank(M::ModuleCtxNmod)
  return M.basis.r
end

function module_trafo_assure(M::ModuleCtx_fmpz)
  if !M.new && isdefined(M, :trafo)
    return
  end
  z = vcat(M.bas_gens, M.rel_gens)
  h, t = hnf_kannan_bachem(z, Val{true})
  M.trafo = t
  M.basis = h
  M.basis_idx = det(h) # h is upp_triangular, hence det is trivial
  M.new = false
  nothing
end


text:
# This could work for  ZigZag1d as well
"""
    move_forward!(τ, t, x, θ, Z::Union{BouncyParticle, ZigZag})

Updates the position `x`, velocity `θ` and time `t` of the
process after a time step equal to `τ` according to the deterministic
dynamics of the Bouncy particle sampler (`BouncyParticle`) and `ZigZag`:
(x(τ), θ(τ)) = (x(0) + θ(0)*t, θ(0)).
`x`: current location, `θ`: current velocity, `t`: current time,
"""
function move_forward!(τ, t, x, θ, Z::Union{BouncyParticle, ZigZag})
    t += τ
    x .+= θ .* τ
    t, x, θ
end

# This one could work for Boomerang 1d as well
"""
    move_forward!(τ, t, x, θ, B::Boomerang)
Updates the position `x`, velocity `θ` and time `t` of the
process after a time step equal to `τ` according to the deterministic
dynamics of the `Boomerang` sampler which are the Hamiltonian
dynamics preserving the Gaussian measure:
: x_t = μ +(x_0 − μ)*cos(t) + v_0*sin(t), v_t = −(x_0 − μ)*sin(t) + v_0*cos(t)
`x`: current location, `θ`: current velocity, `t`: current time.
"""
function move_forward!(τ, t, x, θ, B::Union{Boomerang, FactBoomerang})
    s, c = sincos(τ)
    for i in eachindex(x)
        x[i], θ[i] = (x[i] - B.μ[i])*c + θ[i]*s + B.μ[i],
                    -(x[i] - B.μ[i])*s + θ[i]*c
    end
    t + τ, x, θ
end

"""
        reflect!(i, x, θ, F)

Reflection rule of sampler `F` at reflection time.
`i`: coordinate which flips sign, `x`: position, `θ`: velocity (position
not used for the `ZigZag` and `FactBoomerang`.)
"""
function reflect!(i, ∇ϕx::Number, x, θ, F::Union{ZigZag, FactBoomerang})
    θ[i] = -θ[i]
    θ
end
function reflect!(i, ∇ϕx, x, θ, F::Union{ZigZag, FactBoomerang})
    θ[i] = θ[i] - (2*dot(∇ϕx, θ[i])/normsq(∇ϕx))*∇ϕx
    θ
end

"""
    move_forward(τ, t, x, θ, ::ZigZag1d)

Updates the position `x`, velocity `θ` and time `t` of the
process after a time step equal to `τ` according to the deterministic
dynamics of the `ZigZag1d` sampler: (x(τ), θ(τ)) = (x(0) + θ(0)*t, θ(0)).
`x`: current location, `θ`: current velocity, `t`: current time,
"""
function move_forward(τ, t, x, θ, ::ZigZag1d)
    τ + t, x + θ*τ , θ
end

"""
    move_forward(τ, t, x, θ, B::Boomerang1d)
Updates the position `x`, velocity `θ` and time `t` of the
process after a time step equal to `τ` according to the deterministic
dynamics of the `Boomerang1d` sampler: x_t = μ +(x_0 − μ)*cos(t) + v_0*sin(t),
v_t = −(x_0 − μ)*sin(t) + v_0*cos(t)
`x`: current location, `θ`: current velocity, `t`: current time.
"""
function move_forward(τ, t, x, θ, B::Boomerang1d)
    s, c = sincos(τ)
    t + τ, (x - B.μ)*c + θ*s + B.μ, -(x - B.μ)*s + θ*c
end

"""
    reflect!(∇ϕx, θ, F::BouncyParticle, Boomerang)

Reflection rule of sampler `F` at reflection time.
x`: position, `θ`: velocity
"""
function reflect!(∇ϕx, x, θ, F::BouncyParticle)
    θ .-= (2*dot(∇ϕx, θ)/normsq(F.L\∇ϕx))*(F.L'\(F.L\∇ϕx))
    θ
end
function reflect!(∇ϕx, x, θ, F::Boomerang)
    θ .-= (2*dot(∇ϕx, θ)/normsq(F.L\∇ϕx))*(F.L'\(F.L\∇ϕx))
#    θ .-= (2*dot(∇ϕx, θ)/normsq(∇ϕx))*∇ϕx
    θ
end

waiting_time_ref(rng, F) = poisson_time(rng, F.λref)
waiting_time_ref(F) = poisson_time(F.λref)

function refresh!(rng, θ, F)
    ρ̄ = sqrt(1-F.ρ^2)
    @inbounds for i in eachindex(θ)
        θ[i] = F.ρ*θ[i] + ρ̄*randn(rng)
    end
    θ
end

function refresh!(rng, θ, F::BouncyParticle)
    ρ̄ = sqrt(1-F.ρ^2)
    θ .*= F.ρ
    u = ρ̄*(F.L'\randn(rng, length(θ)))
    θ .+= u
    θ
end
function refresh!(rng, θ, F::Boomerang)
    ρ̄ = sqrt(1-F.ρ^2)
    θ .*= F.ρ
    u = ρ̄*(F.L'\randn(rng, length(θ)))
    θ .+= u
    θ
end