text:
export BitStr, @bit_str, @lbit_str, BitStr64, LongBitStr, bit_literal
export bcat, onehot, onehot_batch, buffer

"""
    BitStr{N,T} <: Integer

struct for bit string with fixed length `N`, the storage type is `T`.

    BitStr{N,T}(value)
    BitStr64{N}(value)
    LongBitStr{N}(value)

Returns a `BitStr`.

## Example

`BitStr` supports some basic arithmetic operations. It acts like an integer, but supports
some frequently used methods for binary basis.

```julia
julia> bit"0101" * 2
1010 ₍₂₎

julia> bcat(bit"101" for i in 1:10)
101101101101101101101101101101 (766958445)

julia> repeat(bit"101", 2)
101101 ₍₂₎

julia> bit"1101"[2]
0
```
"""
struct BitStr{N,T} <: Integer
    buf::T
end

const BitStr64{N} = BitStr{N,Int64}
const LongBitStr{N} = BitStr{N,BigInt}

BitStr{N,T}(val::BitStr{N,T}) where {N,T} = val
Base.zero(::Type{BitStr{N,T}}) where {N,T} = BitStr{N,T}(zero(T))
Base.zero(::BitStr{N,T}) where {N,T} = BitStr{N,T}(zero(T))

buffer(b::BitStr) = b.buf
Base.reinterpret(::Type{BitStr{N,T}}, x::Integer) where {N,T} = BitStr{N,T}(reinterpret(T,x))
Base.reinterpret(::Type{T}, x::BitStr) where {T} = reinterpret(T, buffer(x))
Base.convert(::Type{T}, b::BitStr) where T<:Integer = convert(T, buffer(b))
Base.convert(::Type{T}, b::Integer) where T<:BitStr = T(b)
Base.convert(::Type{T1}, b::BitStr{N2,T2}) where {T1<:BitStr,N2,T2} = convert(T1, buffer(b))
#Base.promote_rule(::Type{BitStr{N,T1}}, ::Type{BitStr{N,T2}}) where {N,T1,T2} = BitStr{N,promote_rule(T1,T2)}
for IT in [:BigInt, :Int128, :UInt128, :Int64,:UInt64, :Int32, :UInt32, :Int16, :UInt16, :Int8, :UInt8, :Bool]
    @eval Base.$IT(b::BitStr) = $IT(buffer(b))
end
for op in [:+, :-, :*, :÷, :|, :⊻, :&, :%, :mod, :mod1]
    @eval Base.$op(a::T, b::Integer) where T<:BitStr = T($op(buffer(a),b))
    @eval Base.$op(a::Integer, b::T) where T<:BitStr = T($op(a,buffer(b)))
    @eval Base.$op(a::BitStr{N,T}, b::BitStr{N,T}) where {N,T} = BitStr{N,T}($op(buffer(a), buffer(b)))
    @eval Base.$op(a::BitStr, b::BitStr) = error("type mismatch: $(typeof(a)), $(typeof(b))")
end
Base.:-(x::BitStr{N,T}) where {N,T} = BitStr{N,T}(-buffer(x))

for op in [:(>>), :(<<)]
    @eval Base.$op(a::BitStr{N,T}, b::Int) where {N,T} = BitStr{N,T}(Base.$op(buffer(a),b))
    #@eval Base.$op(a::T, b::T) where T<:BitStr = T(Base.$op(buffer(a),buffer(b)))
end

for op in [:<, :>, :(<=), :(>=)]
    @eval Base.$op(a::T, b::T) where T<:BitStr = Base.$op(buffer(a),buffer(b))
end

for op in [:(==)]
    @eval Base.$op(a::T, b::Number) where T<:BitStr = Base.$op(buffer(a),b)
    @eval Base.$op(a::Number, b::T) where T<:BitStr = Base.$op(a,buffer(b))
    @eval Base.$op(a::BitStr{N}, b::BitStr{N}) where N = Base.$op(buffer(a),buffer(b))
end
for op in [:count_ones, :count_zeros, :leading_ones, :leading_zeros]
    @eval Base.$op(a::BitStr) = Base.$op(buffer(a))
end

# Note: the transitivity of == is not satisfied here.
Base.:(==)(lhs::BitStr, rhs::BitStr) = false
Base.isapprox(a::BitStr, b::Number; kwargs...) = Base.isapprox(buffer(a),b; kwargs...)
Base.isapprox(a::Number, b::BitStr; kwargs...) = Base.isapprox(a,buffer(b); kwargs...)
Base.isapprox(lhs::BitStr, rhs::BitStr; kwargs...) = false
Base.isapprox(a::T, b::T; kwargs...) where T<:BitStr = Base.isapprox(buffer(a),buffer(b); kwargs...)

# Note: it is a bit confusing, with x::BitStr == y::Integer,
# they behave different when used for indexing.
Base.to_index(x::BitStr) = error("please do not use bit string for indexing, you may want to use `buffer(x)+1` for indexing to avoid ambiguity.")
Base.to_index(x::UnitRange{<:BitStr}) = error("please do not use bit string for indexing, you may want to use `buffer(x)+1` for indexing to avoid ambiguity.")

# use system interface
Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, i::BitStr) =
    checkindex(Bool, inds, Base.to_index(i))
Base.length(bits::BitStr{N}) where N = N
Base.lastindex(bits::BitStr) = length(bits)

Base.typemax(::Type{BitStr{N,T}}) where {N,T} = BitStr{N,T}(1<<N-1)
Base.typemin(::Type{BitStr{N,T}}) where {N,T} = BitStr{N,T}(0)
Base.typemax(::BitStr{N,T}) where {N,T} = BitStr{N,T}(1<<N-1)
Base.typemin(::BitStr{N,T}) where {N,T} = BitStr{N,T}(0)

"""
    @bit_str -> BitStr64

Construct a bit string. such as `bit"0000"`. The bit strings also supports string `bcat`. Just use
it like normal strings.

## Example

```jldoctest
julia> bit"10001"
10001 ₍₂₎

julia> bit"100_111_101"
100111101 ₍₂₎

julia> bcat(bit"1001", bit"11", bit"1110")
1001111110 ₍₂₎

julia> onehot(bit"1001")
16-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 1.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

```
"""
macro bit_str(str)
    return parse_bit(Int64, str)
end

"""
    @bit_str -> LongBitStr

Long bit string version of `@bit_str` macro.
"""
macro lbit_str(str)
    return parse_bit(BigInt, str)
end

function parse_bit(::Type{T}, str::String) where {T <: Integer}
    val = zero(T); k = 1
    for each in reverse(filter(x->x!='_', str))
        if each == '1'
            val += one(T) << (k - 1)
            k += 1
        elseif each == '0'
            k += 1
        elseif each == '_'
            continue
        else
            error("expect 0 or 1, got $each at $k-th bit")
        end
        (isbitstype(T) && k>bsizeof(T)) && error("string length is larger than $(bsizeof(T)), use @lbit_str instead")
    end
    return BitStr{k-1,T}(val)
end

sum_length(a::BitStr, bits::BitStr...) = length(a) + sum_length(bits...)
sum_length(a::BitStr) = length(a)

function bcat(bits::(BitStr{N,T} where N)...) where T
    total_bits = sum_length(bits...)
    val, len = zero(T), 0

    for k in length(bits):-1:1
        val += buffer(bits[k]) << len
        len += length(bits[k])
    end
    return BitStr{total_bits,T}(val)
end

# expand iterator to tuple
bcat(bits) = bcat(bits...)

Base.@propagate_inbounds function Base.getindex(bit::BitStr{N}, index::Int) where N
    @boundscheck 1 <= index <= N || throw(BoundsError(bit, index))
    return buffer(readbit(bit, index))
end

Base.@propagate_inbounds function Base.getindex(bit::BitStr{N}, itr::Union{AbstractVector, AbstractRange}) where N
    @boundscheck all(x->1<=x<=N, itr) || throw(BoundsError(bit, itr))
    return map(x->buffer(readbit(bit, x)), itr)
end

# TODO: support AbstractArray, should return its corresponding shape

Base.@propagate_inbounds function Base.getindex(bit::BitStr{N,T}, mask::Union{Vector{Bool}, BitArray}) where {N,T}
    @boundscheck N == length(mask) || error("length of bits and mask does not match.")

    out = T[]
    for k in eachindex(mask)
        if mask[k]
            push!(out, bit[k])
        end
    end
    return out
end

Base.eltype(::BitStr{N,T}) where {N,T} = T

function Base.iterate(bit::BitStr, state::Integer=1)
    if state > length(bit)
        return nothing
    else
        return bit[state], state + 1
    end
end
Base.IteratorSize(::BitStr) = Base.HasLength()

Base.repeat(s::BitStr, n::Integer) = bcat(s for i in 1:n)
Base.show(io::IO, bitstr::BitStr64{N}) where N = print(io, string(buffer(bitstr), base=2, pad=N), " ₍₂₎")
Base.show(io::IO, bitstr::LongBitStr{N}) where N = print(io, join(map(string, [bitstr[end:-1:1]...])), " ₍₂₎")

"""
    onehot([T=Float64], bit_str[, nbatch])

Returns an onehot vector in type `Vector{T}`, or a batch of onehot
vector in type `Matrix{T}`, where the `bit_str`-th element is one.
"""
onehot(::Type{T}, n::BitStr{N}) where {T,N} = onehot(T, N, buffer(n))
onehot(n::BitStr) = onehot(Float64, n)

onehot(::Type{T}, n::BitStr{N}, nbatch::Int) where {T,N} = onehot(T, N, buffer(n), nbatch)
onehot(n::BitStr, nbatch::Int) = onehot(Float64, n, nbatch)

# operations
"""
    breflect(bit_str[, masks])

Return left-right reflected bit string.
"""
breflect(b::BitStr{N}) where N = breflect(b; nbits=N)
breflect(b::BitStr{N,T}, masks::Vector{<:BitStr{N,T}}) where {N,T} = BitStr{N,T}(breflect(buffer(b), reinterpret(T,masks); nbits=N))

"""
    neg(b::BitStr) -> BitStr
"""
neg(b::BitStr{N}) where N = neg(b, N)

"""
    bfloat(b::BitStr) -> Float64
"""
bfloat(b::BitStr{N}) where N = bfloat(buffer(b); nbits=N)

"""
    bfloat_r(b::BitStr) -> Float64
"""
bfloat_r(b::BitStr{N}) where N = bfloat_r(buffer(b); nbits=N)

"""
    bint_r(b::BitStr) -> Integer
"""
bint_r(b::BitStr{N}) where N = buffer(breflect(b))

"""
    bint(b::BitStr) -> Integer
"""
bint(b::BitStr) = buffer(b)

"""
    bit_literal(xs...)

Create a [`BitStr`](@ref) by input bits `xs`.

# Example

```jldoctest
julia> bit_literal(1, 0, 1, 0, 1, 1)
110101 ₍₂₎
```
"""
bit_literal(xs...) = bit_literal(xs)
function bit_literal(xs::NTuple{N, T}) where {N, T<:Integer}
    val = T(0)
    for k in 1:N
        xs[k] == 0 || xs[k] == 1 || error("expect 0 or 1, got $(xs[k])")
        val += xs[k] << (k - 1)
    end
    return BitStr64{N}(val)
end

basis(b::BitStr) = typemin(b):typemax(b)
basis(::Type{BT}) where BT<:BitStr = typemin(BT):typemax(BT)


text:
using DataDeps
using Test
using ExpectationStubs

# HACK: todo, work out how ExpectationStubs should be changed to make this make sense
Base.open(stub::Stub, t::Any, ::AbstractString) = stub(t)

withenv("DATADEPS_ALWAYS_ACCEPT"=>"true") do
    @testset "automatic download" begin
        @stub dummyhash
        @expect(dummyhash(::Any) = [0x12, 0x34])

        @stub dummydown
        @expect(dummydown("http://www.example.com/eg.zip", ::String) = joinpath(@__DIR__, "eg.zip"))
            # HACK: this gives a directory, because we can't trivially mock out the cd

        register(DataDep( "Test1",
         "A dummy message",
         "http://www.example.com/eg.zip",
         (dummyhash, "1234"),#
         fetch_method=dummydown
        ))

        @test endswith(datadep"Test1", "Test1") || endswith(datadep"Test1", "Test1/") ||  endswith(datadep"Test1", "Test1\\")

        @test all_expectations_used(dummyhash)
        @test all_expectations_used(dummydown)

        rm(datadep"Test1"; recursive=true) # delete the directory
    end

    @testset "sanity check the macro's behaviour with variables" begin
        var = "foo/bar"
        macroexpand(@__MODULE__, :(@datadep_str var)) # this line would throw an error if the varibles were being handle wrong
        @test true
    end

end

ext: jl


text:
module SentinelMissings
struct SentinelMissing{T,SV} <: Number
    x::T
end
Base.promote_rule(::Type{SentinelMissing{SM,SV}},::Union{T,Missing}) where {SM, SV, T}=Union{promote_type(SM,T),Missing}
Base.promote_rule(::Type{SentinelMissing{SM,SV}},::Type{T}) where {SM,SV,T} = Union{promote_type(SM,T),Missing}
Base.promote_rule(::Type{SentinelMissing{SM,SV}},::Type{SentinelMissing{SM,SV}}) where {SM,SV} = Union{SM,Missing}
Base.convert(::Type{Union{T,Missing}},sm::SentinelMissing) where {T}  = ismissing(sm) ? missing : convert(T,sm[])
Base.convert(::Type{T},sm::SentinelMissing) where T<:Number = convert(T,sm[])
Base.getindex(x::SentinelMissing)=ismissing(x) ? missing : x.x
Base.convert(::Type{Any}, sm::SentinelMissing) = sm[]
Base.ismissing(x::SentinelMissing) = isequal(x.x,senval(x))
Base.nonmissingtype(::Type{<:SentinelMissing{T}}) where T = T
Base.isless(x::SentinelMissing,y::SentinelMissing) = isless(x[],y[])
function Base.convert(::Type{<:T}, x::Number) where T<:SentinelMissing
  sv = senval(T)
  et = eltype(T)
  SentinelMissing{et,sv}(convert(et,x))
end
Base.convert(::Type{<:T}, x::Missing) where T<:SentinelMissing = SentinelMissing{eltype(T),senval(T)}(senval(T))
for op in (:(+),:(-),:(*),:(/),:(^))
  eval(:(Base.$(op)(x1::T,x2::T) where T<:SentinelMissing{SM} where SM = $(op)(convert(Union{SM,Missing},x1),convert(Union{SM,Missing},x2))))
end
Base.zero(T::Type{<:SentinelMissing{SM}}) where SM = T(zero(SM))
Base.one(T::Type{<:SentinelMissing{SM}}) where SM  = T(one(SM))
Base.convert(::Type{<:T}, x::T) where T<:SentinelMissing = x
senval(::SentinelMissing{<:Any,SV}) where SV = SV
senval(::Type{<:SentinelMissing{<:Any,SV}}) where SV = SV
Base.eltype(::Type{<:SentinelMissing{T}}) where T = T
Base.show(io::IO,x::SentinelMissing)=print(io,x[])
Base.similar(x::AbstractArray{<:SentinelMissing{T,SV}}) where {T,SV} = zeros(Union{T,Missing},size(x))
Base.similar(x::AbstractArray{<:SentinelMissing{T,SV}},dims::Tuple{Vararg{Int64,N}}) where {T,SV,N} = zeros(Union{T,Missing},dims)

"""
    as_sentinel(x, v)

Reinterprets a Number Array or a Number `x` so that values in x that equal v will be treated as missing.
This is done by reinterpreting the array as a `SentinelMissing` without copying the data.
"""
as_sentinel(x::AbstractArray{T},v) where T<:Number = reinterpret(SentinelMissing{T,convert(T,v)},x)
as_sentinel(x::T,v) where T<:Number = SentinelMissing{T,convert(T,v)}(x)
export as_sentinel
end


text:
"""
    harmonic(n::Integer)

 Calculates harmonic numbers,
   e.g., see [http://mathworld.wolfram.com/HarmonicNumber.html](http://mathworld.wolfram.com/HarmonicNumber.html)

## Arguments
* ``n`` `::Integer`: non-negative index of the Harmonic number to calculate

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> harmonic(2)
1.5
```
"""
function harmonic(n::Integer)
    if n < 0
        throw(DomainError(n))
    elseif n==0
        return 0.0 # by convention (see Crandall, p.22)
    elseif n <= 10
        # perform exact sum for small n
        total = 0.0
        for k=1:n
            total +=  1.0 / k
        end
        return total
    else
        return γ + SpecialFunctions.digamma(n+1) # digamma(m) = ψ(m)
    end
end

"""
    harmonic(x::ComplexOrReal{Float64})

 Calculates harmonic numbers extended to non-integer arguments using the
 digamma form.

## Arguments
* ``x`` `::ComplexOrReal{Float64}`: index of the Harmonic number to calculate

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> harmonic(2.0)
1.5000000000000016
```
"""
function harmonic(x::ComplexOrReal{Float64})
    return γ + SpecialFunctions.digamma(x+1)
end
#        ``digamma(m,x) = (-1)^{m+1} m! \zeta(m+1,x), m>=1``, and integer

"""
    harmonic(n::Integer,r::Real)

 Calculates generalized harmonic numbers,
   e.g., see [http://mathworld.wolfram.com/HarmonicNumber.html](http://mathworld.wolfram.com/HarmonicNumber.html)

## Arguments
* ``n`` `::Integer`: non-negative index 1 of the Harmonic number to calculate
* ``r`` `::Real`: index 2 of the Harmonic number to calculate

It should be possible to extend this to complex r, but that requires more testing.

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> harmonic(2,1.5)
1.3535533905932737
```
"""
function harmonic(n::Integer, r::Real)
    if n < 0
        throw(DomainError(n))
    end
    if n == 0
        return 0.0
    end
    if r==1
        return harmonic(n)
    end
    total = 0.0
    for k=1:n
        total +=  1.0 / k^r
    end
    return total
end

"""
    harmonic(n::Integer,r::Integer)

 Calculates generalized harmonic numbers
   e.g., see [http://mathworld.wolfram.com/HarmonicNumber.html](http://mathworld.wolfram.com/HarmonicNumber.html)
 using a better approach which works when both inputs are integers
 [https://carma.newcastle.edu.au/resources/jon/Preprints/Papers/Published-InPress/Oscillatory%20(Tapas%20II)/Papers/coffey-zeta.pdf](https://carma.newcastle.edu.au/resources/jon/Preprints/Papers/Published-InPress/Oscillatory%20(Tapas%20II)/Papers/coffey-zeta.pdf), p.341

## Arguments
* ``n`` `::Integer`: non-negative index 1 of the Harmonic number to calculate
* ``r`` `::Integer`: index 2 of the Harmonic number to calculate

## Examples
```jldoctest; setup = :(using Polylogarithms)
julia> harmonic(2,1)
1.5000000000000002
```
"""
function harmonic(n::Integer, r::Integer)
    if r<1
        throw(DomainError(r))
    end
    return (-1)^(r-1) * ( SpecialFunctions.polygamma(r-1,n+1) - SpecialFunctions.polygamma(r-1,1) ) / SpecialFunctions.gamma(r)
end

# # better approach to calculation, but doesn't seem to work
# #  using
# function harmonic3(x::Real, r::Integer)
#     if (r < 1)
#         throw(DomainError(r))
#     end
#     return SpecialFunctions.zeta(r,1) - SpecialFunctions.zeta(r,x+1)
# end


text:
export ReferenceContigs, ReferenceContigs_hg38

type ReferenceContigs
    count::Int64
    names::Array{ASCIIString}
    sizes::Array{Int64}
    offsets::Array{Int64}

    function ReferenceContigs(count, names, sizes)
        new(count, names, sizes, [sum(sizes[1:i-1]) for i in 1:length(sizes)])
    end
end

ReferenceContigs_hg38 = ReferenceContigs(455, [
    "chr1", "chr10", "chr11", "chr11_KI270721v1_random", "chr12", "chr13", "chr14", "chr14_GL000009v2_random",
    "chr14_GL000225v1_random", "chr14_KI270722v1_random", "chr14_GL000194v1_random", "chr14_KI270723v1_random",
    "chr14_KI270724v1_random", "chr14_KI270725v1_random", "chr14_KI270726v1_random", "chr15", "chr15_KI270727v1_random",
    "chr16", "chr16_KI270728v1_random", "chr17", "chr17_GL000205v2_random", "chr17_KI270729v1_random",
    "chr17_KI270730v1_random", "chr18", "chr19", "chr1_KI270706v1_random", "chr1_KI270707v1_random",
    "chr1_KI270708v1_random", "chr1_KI270709v1_random", "chr1_KI270710v1_random", "chr1_KI270711v1_random",
    "chr1_KI270712v1_random", "chr1_KI270713v1_random", "chr1_KI270714v1_random", "chr2", "chr20", "chr21",
    "chr22", "chr22_KI270731v1_random", "chr22_KI270732v1_random", "chr22_KI270733v1_random", "chr22_KI270734v1_random", "chr22_KI270735v1_random", "chr22_KI270736v1_random", "chr22_KI270737v1_random", "chr22_KI270738v1_random",
    "chr22_KI270739v1_random", "chr2_KI270715v1_random", "chr2_KI270716v1_random", "chr3", "chr3_GL000221v1_random",
    "chr4", "chr4_GL000008v2_random", "chr5", "chr5_GL000208v1_random", "chr6", "chr7", "chr8", "chr9",
    "chr9_KI270717v1_random", "chr9_KI270718v1_random", "chr9_KI270719v1_random", "chr9_KI270720v1_random",
    "chr1_KI270762v1_alt", "chr1_KI270766v1_alt", "chr1_KI270760v1_alt", "chr1_KI270765v1_alt", "chr1_GL383518v1_alt",
    "chr1_GL383519v1_alt", "chr1_GL383520v2_alt", "chr1_KI270764v1_alt", "chr1_KI270763v1_alt", "chr1_KI270759v1_alt",
    "chr1_KI270761v1_alt", "chr2_KI270770v1_alt", "chr2_KI270773v1_alt", "chr2_KI270774v1_alt", "chr2_KI270769v1_alt",
    "chr2_GL383521v1_alt", "chr2_KI270772v1_alt", "chr2_KI270775v1_alt", "chr2_KI270771v1_alt", "chr2_KI270768v1_alt",
    "chr2_GL582966v2_alt", "chr2_GL383522v1_alt", "chr2_KI270776v1_alt", "chr2_KI270767v1_alt", "chr3_JH636055v2_alt",
    "chr3_KI270783v1_alt", "chr3_KI270780v1_alt", "chr3_GL383526v1_alt", "chr3_KI270777v1_alt", "chr3_KI270778v1_alt",
    "chr3_KI270781v1_alt", "chr3_KI270779v1_alt", "chr3_KI270782v1_alt", "chr3_KI270784v1_alt", "chr4_KI270790v1_alt",
    "chr4_GL383528v1_alt", "chr4_KI270787v1_alt", "chr4_GL000257v2_alt", "chr4_KI270788v1_alt", "chr4_GL383527v1_alt",
    "chr4_KI270785v1_alt", "chr4_KI270789v1_alt", "chr4_KI270786v1_alt", "chr5_KI270793v1_alt", "chr5_KI270792v1_alt",
    "chr5_KI270791v1_alt", "chr5_GL383532v1_alt", "chr5_GL949742v1_alt", "chr5_KI270794v1_alt", "chr5_GL339449v2_alt",
    "chr5_GL383530v1_alt", "chr5_KI270796v1_alt", "chr5_GL383531v1_alt", "chr5_KI270795v1_alt", "chr6_GL000250v2_alt",
    "chr6_KI270800v1_alt", "chr6_KI270799v1_alt", "chr6_GL383533v1_alt", "chr6_KI270801v1_alt", "chr6_KI270802v1_alt",
    "chr6_KB021644v2_alt", "chr6_KI270797v1_alt", "chr6_KI270798v1_alt", "chr7_KI270804v1_alt", "chr7_KI270809v1_alt",
    "chr7_KI270806v1_alt", "chr7_GL383534v2_alt", "chr7_KI270803v1_alt", "chr7_KI270808v1_alt", "chr7_KI270807v1_alt",
    "chr7_KI270805v1_alt", "chr8_KI270818v1_alt", "chr8_KI270812v1_alt", "chr8_KI270811v1_alt", "chr8_KI270821v1_alt",
    "chr8_KI270813v1_alt", "chr8_KI270822v1_alt", "chr8_KI270814v1_alt", "chr8_KI270810v1_alt", "chr8_KI270819v1_alt",
    "chr8_KI270820v1_alt", "chr8_KI270817v1_alt", "chr8_KI270816v1_alt", "chr8_KI270815v1_alt", "chr9_GL383539v1_alt",
    "chr9_GL383540v1_alt", "chr9_GL383541v1_alt", "chr9_GL383542v1_alt", "chr9_KI270823v1_alt", "chr10_GL383545v1_alt",
    "chr10_KI270824v1_alt", "chr10_GL383546v1_alt", "chr10_KI270825v1_alt", "chr11_KI270832v1_alt", "chr11_KI270830v1_alt",
    "chr11_KI270831v1_alt", "chr11_KI270829v1_alt", "chr11_GL383547v1_alt", "chr11_JH159136v1_alt", "chr11_JH159137v1_alt",
    "chr11_KI270827v1_alt", "chr11_KI270826v1_alt", "chr12_GL877875v1_alt", "chr12_GL877876v1_alt", "chr12_KI270837v1_alt",
    "chr12_GL383549v1_alt", "chr12_KI270835v1_alt", "chr12_GL383550v2_alt", "chr12_GL383552v1_alt", "chr12_GL383553v2_alt",
    "chr12_KI270834v1_alt", "chr12_GL383551v1_alt", "chr12_KI270833v1_alt", "chr12_KI270836v1_alt", "chr13_KI270840v1_alt",
    "chr13_KI270839v1_alt", "chr13_KI270843v1_alt", "chr13_KI270841v1_alt", "chr13_KI270838v1_alt", "chr13_KI270842v1_alt",
    "chr14_KI270844v1_alt", "chr14_KI270847v1_alt", "chr14_KI270845v1_alt", "chr14_KI270846v1_alt", "chr15_KI270852v1_alt",
    "chr15_KI270851v1_alt", "chr15_KI270848v1_alt", "chr15_GL383554v1_alt", "chr15_KI270849v1_alt", "chr15_GL383555v2_alt",
    "chr15_KI270850v1_alt", "chr16_KI270854v1_alt", "chr16_KI270856v1_alt", "chr16_KI270855v1_alt", "chr16_KI270853v1_alt",
    "chr16_GL383556v1_alt", "chr16_GL383557v1_alt", "chr17_GL383563v3_alt", "chr17_KI270862v1_alt", "chr17_KI270861v1_alt",
    "chr17_KI270857v1_alt", "chr17_JH159146v1_alt", "chr17_JH159147v1_alt", "chr17_GL383564v2_alt", "chr17_GL000258v2_alt",
    "chr17_GL383565v1_alt", "chr17_KI270858v1_alt", "chr17_KI270859v1_alt", "chr17_GL383566v1_alt", "chr17_KI270860v1_alt",
    "chr18_KI270864v1_alt", "chr18_GL383567v1_alt", "chr18_GL383570v1_alt", "chr18_GL383571v1_alt", "chr18_GL383568v1_alt",
    "chr18_GL383569v1_alt", "chr18_GL383572v1_alt", "chr18_KI270863v1_alt", "chr19_KI270868v1_alt", "chr19_KI270865v1_alt",
    "chr19_GL383573v1_alt", "chr19_GL383575v2_alt", "chr19_GL383576v1_alt", "chr19_GL383574v1_alt", "chr19_KI270866v1_alt",
    "chr19_KI270867v1_alt", "chr19_GL949746v1_alt", "chr20_GL383577v2_alt", "chr20_KI270869v1_alt", "chr20_KI270871v1_alt",
    "chr20_KI270870v1_alt", "chr21_GL383578v2_alt", "chr21_KI270874v1_alt", "chr21_KI270873v1_alt", "chr21_GL383579v2_alt",
    "chr21_GL383580v2_alt", "chr21_GL383581v2_alt", "chr21_KI270872v1_alt", "chr22_KI270875v1_alt", "chr22_KI270878v1_alt",
    "chr22_KI270879v1_alt", "chr22_KI270876v1_alt", "chr22_KI270877v1_alt", "chr22_GL383583v2_alt", "chr22_GL383582v2_alt",
    "chrX_KI270880v1_alt", "chrX_KI270881v1_alt", "chr19_KI270882v1_alt", "chr19_KI270883v1_alt", "chr19_KI270884v1_alt",
    "chr19_KI270885v1_alt", "chr19_KI270886v1_alt", "chr19_KI270887v1_alt", "chr19_KI270888v1_alt", "chr19_KI270889v1_alt",
    "chr19_KI270890v1_alt", "chr19_KI270891v1_alt", "chr1_KI270892v1_alt", "chr2_KI270894v1_alt", "chr2_KI270893v1_alt",
    "chr3_KI270895v1_alt", "chr4_KI270896v1_alt", "chr5_KI270897v1_alt", "chr5_KI270898v1_alt", "chr6_GL000251v2_alt",
    "chr7_KI270899v1_alt", "chr8_KI270901v1_alt", "chr8_KI270900v1_alt", "chr11_KI270902v1_alt", "chr11_KI270903v1_alt",
    "chr12_KI270904v1_alt", "chr15_KI270906v1_alt", "chr15_KI270905v1_alt", "chr17_KI270907v1_alt", "chr17_KI270910v1_alt",
    "chr17_KI270909v1_alt", "chr17_JH159148v1_alt", "chr17_KI270908v1_alt", "chr18_KI270912v1_alt", "chr18_KI270911v1_alt",
    "chr19_GL949747v2_alt", "chr22_KB663609v1_alt", "chrX_KI270913v1_alt", "chr19_KI270914v1_alt", "chr19_KI270915v1_alt",
    "chr19_KI270916v1_alt", "chr19_KI270917v1_alt", "chr19_KI270918v1_alt", "chr19_KI270919v1_alt", "chr19_KI270920v1_alt",
    "chr19_KI270921v1_alt", "chr19_KI270922v1_alt", "chr19_KI270923v1_alt", "chr3_KI270924v1_alt", "chr4_KI270925v1_alt",
    "chr6_GL000252v2_alt", "chr8_KI270926v1_alt", "chr11_KI270927v1_alt", "chr19_GL949748v2_alt", "chr22_KI270928v1_alt",
    "chr19_KI270929v1_alt", "chr19_KI270930v1_alt", "chr19_KI270931v1_alt", "chr19_KI270932v1_alt", "chr19_KI270933v1_alt",
    "chr19_GL000209v2_alt", "chr3_KI270934v1_alt", "chr6_GL000253v2_alt", "chr19_GL949749v2_alt", "chr3_KI270935v1_alt",
    "chr6_GL000254v2_alt", "chr19_GL949750v2_alt", "chr3_KI270936v1_alt", "chr6_GL000255v2_alt", "chr19_GL949751v2_alt",
    "chr3_KI270937v1_alt", "chr6_GL000256v2_alt", "chr19_GL949752v1_alt", "chr6_KI270758v1_alt", "chr19_GL949753v2_alt",
    "chr19_KI270938v1_alt", "chrM", "chrUn_KI270302v1", "chrUn_KI270304v1", "chrUn_KI270303v1", "chrUn_KI270305v1",
    "chrUn_KI270322v1", "chrUn_KI270320v1", "chrUn_KI270310v1", "chrUn_KI270316v1", "chrUn_KI270315v1", "chrUn_KI270312v1",
    "chrUn_KI270311v1", "chrUn_KI270317v1", "chrUn_KI270412v1", "chrUn_KI270411v1", "chrUn_KI270414v1", "chrUn_KI270419v1",
    "chrUn_KI270418v1", "chrUn_KI270420v1", "chrUn_KI270424v1", "chrUn_KI270417v1", "chrUn_KI270422v1", "chrUn_KI270423v1",
    "chrUn_KI270425v1", "chrUn_KI270429v1", "chrUn_KI270442v1", "chrUn_KI270466v1", "chrUn_KI270465v1", "chrUn_KI270467v1",
    "chrUn_KI270435v1", "chrUn_KI270438v1", "chrUn_KI270468v1", "chrUn_KI270510v1", "chrUn_KI270509v1", "chrUn_KI270518v1",
    "chrUn_KI270508v1", "chrUn_KI270516v1", "chrUn_KI270512v1", "chrUn_KI270519v1", "chrUn_KI270522v1", "chrUn_KI270511v1",
    "chrUn_KI270515v1", "chrUn_KI270507v1", "chrUn_KI270517v1", "chrUn_KI270529v1", "chrUn_KI270528v1", "chrUn_KI270530v1",
    "chrUn_KI270539v1", "chrUn_KI270538v1", "chrUn_KI270544v1", "chrUn_KI270548v1", "chrUn_KI270583v1", "chrUn_KI270587v1",
    "chrUn_KI270580v1", "chrUn_KI270581v1", "chrUn_KI270579v1", "chrUn_KI270589v1", "chrUn_KI270590v1", "chrUn_KI270584v1",
    "chrUn_KI270582v1", "chrUn_KI270588v1", "chrUn_KI270593v1", "chrUn_KI270591v1", "chrUn_KI270330v1", "chrUn_KI270329v1",
    "chrUn_KI270334v1", "chrUn_KI270333v1", "chrUn_KI270335v1", "chrUn_KI270338v1", "chrUn_KI270340v1", "chrUn_KI270336v1",
    "chrUn_KI270337v1", "chrUn_KI270363v1", "chrUn_KI270364v1", "chrUn_KI270362v1", "chrUn_KI270366v1", "chrUn_KI270378v1",
    "chrUn_KI270379v1", "chrUn_KI270389v1", "chrUn_KI270390v1", "chrUn_KI270387v1", "chrUn_KI270395v1", "chrUn_KI270396v1",
    "chrUn_KI270388v1", "chrUn_KI270394v1", "chrUn_KI270386v1", "chrUn_KI270391v1", "chrUn_KI270383v1", "chrUn_KI270393v1",
    "chrUn_KI270384v1", "chrUn_KI270392v1", "chrUn_KI270381v1", "chrUn_KI270385v1", "chrUn_KI270382v1", "chrUn_KI270376v1",
    "chrUn_KI270374v1", "chrUn_KI270372v1", "chrUn_KI270373v1", "chrUn_KI270375v1", "chrUn_KI270371v1", "chrUn_KI270448v1",
    "chrUn_KI270521v1", "chrUn_GL000195v1", "chrUn_GL000219v1", "chrUn_GL000220v1", "chrUn_GL000224v1", "chrUn_KI270741v1",
    "chrUn_GL000226v1", "chrUn_GL000213v1", "chrUn_KI270743v1", "chrUn_KI270744v1", "chrUn_KI270745v1", "chrUn_KI270746v1",
    "chrUn_KI270747v1", "chrUn_KI270748v1", "chrUn_KI270749v1", "chrUn_KI270750v1", "chrUn_KI270751v1", "chrUn_KI270752v1",
    "chrUn_KI270753v1", "chrUn_KI270754v1", "chrUn_KI270755v1", "chrUn_KI270756v1", "chrUn_KI270757v1", "chrUn_GL000214v1",
    "chrUn_KI270742v1", "chrUn_GL000216v2", "chrUn_GL000218v1", "chrX", "chrY", "chrY_KI270740v1_random"
], [
    248956422,133797422,135086622,100316,133275309,114364328,107043718,201709,211173,
    194050,191469,38115,39555,172810,43739,101991189,448248,90338345,1872759,83257441,
    185591,280839,112551,80373285,58617616,175055,32032,127682,66860,40176,42210,176043,
    40745,41717,242193529,64444167,46709983,50818468,150754,41543,179772,165050,42811,
    181920,103838,99375,73985,161471,153799,198295559,155397,190214555,209709,181538259,
    92689,170805979,159345973,145138636,138394717,40062,38054,176845,39050,354444,256271,
    109528,185285,182439,110268,366580,50258,911658,425601,165834,136240,70887,223625,
    120616,143390,133041,138019,110395,110099,96131,123821,174166,161578,173151,109187,
    224108,180671,173649,248252,113034,205312,162429,184404,220246,376187,111943,586476,
    158965,164536,119912,205944,244096,126136,179043,195710,82728,226852,164558,1612928,
    101241,172708,173459,131892,4672374,175808,152148,124736,870480,75005,185823,197536,
    271782,157952,209586,158166,119183,1111570,271455,126434,209988,145606,282736,292436,
    985506,300230,624492,141812,374415,133535,36640,158983,305841,132244,162988,71551,
    171286,60032,439082,179254,181496,309802,188315,210133,177092,296895,204059,154407,
    200998,191409,67707,186169,167313,408271,40090,120804,238139,169178,138655,152874,
    119498,184319,76061,56134,191684,180306,103832,169134,306913,37287,322166,1511111,
    180703,1351393,478999,263054,327382,296527,244917,388773,430880,134193,63982,232857,
    2659700,192462,89672,375691,391357,196688,2877074,278131,70345,133151,1821992,223995,
    235827,108763,90219,178921,111737,289831,164789,198278,104552,167950,159547,167999,
    61734,52969,385657,170222,188024,155864,43156,233762,987716,128386,118774,58661,
    183433,63917,166743,143900,201197,74653,116689,82692,259914,186262,304135,263666,
    101331,96924,162811,284869,144206,248807,170399,157053,171027,204239,209512,155532,
    170698,184499,170680,162212,214158,161218,162896,378547,1144418,130957,4795265,
    190869,136959,318687,106711,214625,572349,196384,5161414,137721,157099,325800,88070,
    1423190,174061,157710,729520,74013,274009,205194,170665,184516,190932,123111,170701,
    198005,282224,187935,189352,166540,555799,4604811,229282,218612,1064304,176103,
    186203,200773,170148,215732,170537,177381,163458,4677643,1091841,197351,4827813,
    1066390,164170,4606388,1002683,165607,4929269,987100,76752,796479,1066800,16569,2274,
    2165,1942,1472,21476,4416,1201,1444,2276,998,12399,37690,1179,2646,2489,1029,2145,
    2321,2140,2043,1445,981,1884,1361,392061,1233,1774,3920,92983,112505,4055,2415,2318,
    2186,1951,1300,22689,138126,5674,8127,6361,5353,3253,1899,2983,2168,993,91309,1202,
    1599,1400,2969,1553,7046,31033,44474,4685,4513,6504,6158,3041,5796,1652,1040,1368,
    2699,1048,1428,1428,1026,1121,1803,2855,3530,8320,1048,1045,1298,2387,1537,1143,1880,
    1216,970,1788,1484,1750,1308,1658,971,1930,990,4215,1136,2656,1650,1451,2378,2805,
    7992,7642,182896,179198,161802,179693,157432,15008,164239,210658,168472,41891,66486,
    198735,93321,158759,148850,150742,27745,62944,40191,36723,79590,71251,137718,186739,
    176608,161147,156040895,57227415,37240
]);


text:
using RecipesBase

@recipe f(::Type{InvariantDistribution}, ivd::InvariantDistribution) = ivd.dist


text:
module Types

export GraphQLSchema
export ExecutionArgs
export AbstractGraphQLServer, GraphQLServer, ServerConfiguration
export AbstractGraphQLClient, GraphQLClient, QueryOptions
export DEFAULT_QUERY, DEFAULT_TYPEDEFS

import Base.LibuvServer

#######
# Language Types
#######
struct GraphQLSchema
  noop::String
end

struct ExecutionArgs
  noop::String
end

#######
# Server Types
#######
DEFAULT_TYPEDEFS = """
type HealthCheck {
  ping: boolean
  status: string
}
type query {
  healthcheck: HealthCheck
}
"""
DEFAULT_PORT = 3968

Base.@kwdef struct ServerConfiguration
  typedefs::String = DEFAULT_TYPEDEFS
  port::Int64 = DEFAULT_PORT
end

abstract type AbstractGraphQLServer end
struct GraphQLServer <: AbstractGraphQLServer
  server::LibuvServer
  address::String
  config::ServerConfiguration
end

#######
# Client Types
#######
DEFAULT_QUERY = """
query Ping {
  healthcheck {
    ping
  }
}
"""

abstract type AbstractGraphQLClient end
Base.@kwdef struct GraphQLClient <: AbstractGraphQLClient
  url::String = "http://localhost:$DEFAULT_PORT"
end

Base.@kwdef struct QueryOptions
  query::String = DEFAULT_QUERY
  variables::Dict{String,Any} = Dict()
end

end


text:
module NashGames

include("api.jl")
include("problems.jl")

end


text:
roomTiles(layer::Layer, room::Maple.Room) = layerName(layer) == "fgTiles" ? room.fgTiles : room.bgTiles

function validTiles(layer::String, removeTemplate::Bool=true)
    meta = layer == "fgTiles" ? fgTilerMeta : bgTilerMeta

    res = collect(keys(meta.paths))
    push!(res, '0')

    if removeTemplate
        index = findfirst(isequal('z'), res)
        if index !== nothing
            deleteat!(res, index)
        end
    end

    sort!(res)

    return res
end

validTiles(layer::Layer, removeTemplate::Bool=true) = validTiles(layerName(layer), removeTemplate)

function tileNames(name::String)
    meta = name == "fgTiles" ? fgTilerMeta : bgTilerMeta
    res = Dict{Any, Any}(
        '0' => "Air",
        "Air" => '0'
    )

    for (id, path) in meta.paths
        name = humanizeVariableName(splitdir(path)[2])

        if startswith(name, "Bg ")
            name = name[4:end]
        end

        res[name] = id
        res[id] = name
    end

    return res
end

tileNames(layer::Layer) = tileNames(layerName(layer))

const nodeType = Tuple{Number, Number}
const edgeType = Tuple{nodeType, nodeType}

function shrinkMatrix(m::Array{Bool, 2})
    h, w = size(m)

    lx, ly = h, w
    hx, hy = 1, 1

    for i in 1:h, j in 1:w
        if m[i, j]
            lx = min(lx, j)
            ly = min(ly, i)

            hx = max(hx, j)
            hy = max(hy, i)
        end
    end

    return m[ly:hy, lx:hx], lx, ly
end

function nodesToBrushPixels(nodes::Array{nodeType, 1})
    # Take one pass through to make the output as small as possible
    # Lowest and highest
    lx = ly = typemax(Int)
    hx = hy = typemin(Int)

    for node in nodes
        x, y = node

        lx = min(lx, x)
        ly = min(ly, y)

        hx = max(hx, x)
        hy = max(hy, y)
    end

    w = hx - lx + 1
    h = hy - ly + 1

    # Generate the output
    res = fill(false, (h, w))

    for node in nodes
        x, y = node

        res[y - ly + 1, x - lx + 1] = true
    end

    return res, lx, ly
end

function brushEdges(a::AbstractArray{Bool})
    res = edgeType[]
    h, w = size(a)

    for y in 1:h, x in 1:w
        # Add edges where there are no adjacent value
        if a[y, x]
            if !get(a, (y - 1, x), false)
                push!(res, ((x + 1, y), (x, y)))
            end

            if !get(a, (y + 1, x), false)
                push!(res, ((x, y + 1), (x + 1, y + 1)))
            end

            if !get(a, (y, x - 1), false)
                push!(res, ((x, y + 1), (x, y)))
            end

            if !get(a, (y, x + 1), false)
                push!(res, ((x + 1, y), (x + 1, y + 1)))
            end
        end
    end

    return res
end

function connectEdges!(edges::Array{edgeType, 1})
    res = []

    while length(edges) > 0
        changed = true

        tail, head = pop!(edges)
        nodes = [tail, head]

        while changed
            changed = false

            for (i, edge) in enumerate(edges)
                tail, head = edge

                tailc = tail == nodes[end]
                headc = head == nodes[end]

                if tailc || headc
                    push!(nodes, tailc ? head : tail)
                    deleteat!(edges, i)

                    changed = true

                    break
                end
            end
        end

        push!(res, nodes)
    end

    return res
end

mutable struct Brush
    name::String
    pixels::AbstractArray{Bool}
    offset::Tuple{Number, Number}
    nodes::Array{Array{nodeType, 1}, 1}
    rotation::Number

    Brush(name::String, pixels::AbstractArray{Bool}, offset::Tuple{Number, Number}=(1, 1)) = new(name, pixels, offset, connectEdges!(brushEdges(pixels)), 0)
    Brush(name::String, pixels::Array{T, 2}, offset::Tuple{Number, Number}=(1, 1)) where T <: Integer = Brush(name, pixels .!= 0, offset)
end

function rotationOffset(ctx::CairoContext, brush::Brush)
    bh, bw = size(brush.pixels)
    rotation = brush.rotation

    if rotation == 0
        translate(ctx, 0, 0)

    elseif rotation == 1
        translate(ctx, 0, -bh + 1)

    elseif rotation == 2
        translate(ctx, -bw + 1, -bh + 1)

    elseif rotation == 3
        translate(ctx, -bw + 1, 0)
    end
end

function drawBrush(brush::Brush, layer::Layer, x::Number, y::Number, thickness::Int=2)
    ctx = getSurfaceContext(layer.surface)

    if ctx.ptr != C_NULL
        ox, oy = brush.offset

        Cairo.save(ctx)

        scale(ctx, 8, 8)
        translate(ctx, x + ox - 1.5, y + oy - 1.5)
        rotate(ctx, brush.rotation * pi / 2)
        translate(ctx, -1.5, -1.5)
        rotationOffset(ctx, brush)

        for nodes in deepcopy.(brush.nodes)
            # Fixes weird off by one error in Cairo paths
            # 1 / 8 is one "pixel", taking into account the (8, 8) scale
            # This should be good enough, but might look bad at non default stroke thickness

            extend = ceil(thickness / 2)
            nodes[end] = nodes[end] .+ (extend / 8, extend / 8) .* sign.(nodes[end] .- nodes[end - 1])
            drawLines(ctx, nodes, colors.brush_bc, filled=true, fc=colors.brush_fc, thickness=thickness)
        end

        restore(ctx)
    end
end

function applyBrush!(brush::Brush, tiles::Maple.Tiles, material::Char, x::Number, y::Number)
    brushPixels = rotr90(brush.pixels, brush.rotation)
    tilesData = tiles.data

    bh, bw = size(brushPixels)
    th, tw = size(tilesData)
    ox, oy = brush.offset

    for j in 1:bh, i in 1:bw
        ty, tx = j + y + oy - 2, i + x + ox - 2

        # Manual slice assignment
        # Brushes can be applied "off screen"
        if 1 <= ty <= th && 1 <= tx <= tw
            if brushPixels[j, i]
                tilesData[ty, tx] = material
            end
        end
    end
end

function applyBrush!(brush::Brush, tiles::Maple.Tiles, materials::Array{Char, 2}, x::Number, y::Number)
    brushPixels = rotr90(brush.pixels, brush.rotation)
    #materials = rotr90(materials, brush.rotation)
    tilesData = tiles.data

    bh, bw = size(brushPixels)
    th, tw = size(tilesData)
    ox, oy = brush.offset

    for j in 1:bh, i in 1:bw
        ty, tx = j + y - oy, i + x - ox

        # Manual slice assignment
        # Brushes can be applied "off screen"
        if 1 <= ty <= th && 1 <= tx <= tw
            if brushPixels[j, i]
                tilesData[ty, tx] = materials[j, i]
            end
        end
    end
end

function findFill(tiles::Array{Char, 2}, x::Number, y::Number)
    target = tiles[y, x]
    stack = Tuple{Number, Number}[(x, y)]

    h, w = size(tiles)

    res = fill(false, (h, w))

    while length(stack) > 0
        tx, ty = pop!(stack)

        if 1 <=tx <= w && 1 <= ty <= h && !res[ty, tx]
            if target == tiles[ty, tx]
                res[ty, tx] = true

                push!(stack, (tx - 1, ty))
                push!(stack, (tx, ty - 1))
                push!(stack, (tx + 1, ty))
                push!(stack, (tx, ty + 1))
            end
        end
    end

    return res
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
orig_input = readlines("day7.txt")
Init_state = hcat(split.(orig_input,",")...) .|> x -> parse(Int64,x)

# positions = sort(unique(Init_state))
un = unique(Init_state)
positions = sort(collect(minimum(un):maximum(un)))
all_dist = zeros(Int64, size(positions,1), 2)

all_dist[:,1] = positions

# for (ind, pos) in enumerate(positions)
#     all_dist[ind,2] = sum(abs.(Init_state .- all_dist[ind,1]))
# end

# minimum(all_dist[:,2])

minimum([sum(abs.(Init_state .- all_dist[ind,1])) for ind in (positions .+1)])


text:
@testset "replace!" begin
    y = [1, 2, 2, 2, 3]
    z1 = 1:5
    z2 = 0:3
    r1 = -2:2
    r2 = 2:5
    @test_throws DimensionMismatch _replace!(y, z1, r2)
    @test_throws DimensionMismatch _replace!(y, z2, r1)
    @test _replace!(deepcopy(y), z1, r1) == Base.replace!(deepcopy(y), (z1 .=> r1)...)
    @test _replace!(deepcopy(y), z2, r2) == Base.replace!(deepcopy(y), (z2 .=> r2)...)
end

@testset "softmax" begin
    X = rand(100,5)
    max_ = maximum(X, dims=2)
    exp_ = exp.(X .- max_)
    @test isapprox(MLJMultivariateStatsInterface.softmax(X), exp_ ./ sum(exp_, dims=2))
    X_ = MLJMultivariateStatsInterface.softmax!(X)
    @test X_ === X
    @test isapprox(X_, exp_ ./ sum(exp_, dims=2))
end


text:
using Documenter, ChristensenSet

makedocs(
  sitename="ChristensenSet.jl: Gallery",
  authors="Tomáš Kalvoda",
)


text:
using Documenter, FastaIO

makedocs(
    modules  = [FastaIO],
    format = Documenter.HTML(prettyurls = "--local" ∉ ARGS),
    sitename = "FastaIO.jl",
    pages    = Any[
        "Home" => "index.md",
       ]
    )

deploydocs(
    repo   = "github.com/carlobaldassi/FastaIO.jl.git",
)


text:
struct Field{T}
    name
    description
    computefn
    transformfn
    createfn
    formatfn
    filterfn
end

Base.show(io::IO, field::Field{T}) where T = print(io, "Field{", T, "}(name = \"", field.name, "\")")

"""
    struct Field(T, name; kwargs...)

A field defines what data will be stored in one column of a [`Registry`](#).

## Keyword arguments

- `description::String = ""`: More information on the field contents and how they
    can be used. May contain Markdown formatting.
- `optional = false`: Whether a registry entry can be entered without alue
    in this field.
- `default = missing`: The default value used if `optional === true`. If a
    non-`missing` value is passed, the value of `optional` is ignored and
    the field is treated as optional.
- `defaultfn = (row, key) -> default`: Function to compute dynamically compute
    a default value from an entry. If passed, the value of `default` is ignored.

## Examples

A simple `Field` with entries of type `String`:

```
Field(String, "My Field")
```
"""
function Field(
        T;
        name = "",
        description = "",
        default = missing,
        optional = !ismissing(default),
        defaultfn = optional ? ((row, key) -> default) : ((row, key) -> throw(RequiredKeyMissingError(key))),
        computefn = (row, key) -> haskey(row, key) ? get(row, key, nothing) : defaultfn(row, key),
        U = optional ? Union{typeof(default), T} : T,
        transformfn = x -> convert(U, x),
        containerfn = () -> U[],
        formatfn = identity,
        filterfn = nothing,
    )
    Field{T}(name, description, computefn, transformfn, containerfn, formatfn, filterfn)
end

struct RequiredKeyMissingError <: Exception
    key
end

Base.showerror(e::RequiredKeyMissingError) = print(io,
    "RequiredKeyMissingError: Key `", e.key, "` is not optional")

@testset "Field" begin
    @testset "default" begin
        field = Field(String, name = "Field")
        @test field.computefn((; field = 1,), :field) == 1
        @test_throws RequiredKeyMissingError field.computefn((;), :field)
    end

    @testset "optional" begin
        field = Field(String, name = "Field", optional = true)
        @test field.computefn((; field = 1,), :field) == 1
        @test field.computefn((;), :field) === missing

        field2 = Field(String, name = "Field", optional = true, default = nothing)
        @test field2.computefn((; field = 1,), :field) == 1
        @test field2.computefn((;), :field) === nothing

        field2 = Field(String, name = "Field", optional = true, defaultfn = (row, key) -> key)
        @test field2.computefn((; field = 1,), :field) == 1
        @test field2.computefn((;), :field) === :field
    end

    @testset "container" begin
        @test Field(String, name = "Field").createfn() isa Vector{String}
        @test Field(String, name = "Field", optional=true).createfn() isa Vector{Union{Missing, String}}
    end
end

struct Registry{F<:NamedTuple, N, S<:StructArray, D<:Dict}
    fields::F
    id::NTuple{N, Symbol}
    data::S
    index::D
    name
    description
    loadfn
    # hidecols: columns to be hidden from printing
end

"""
    struct RegistryEntry(row, registry)

An entry in a feature registry `registry` with values `row`.
Returned when indexing into a `FeatureRegistry`, i.e.
`registry[id]`.
"""
struct RegistryEntry
    row
    registry::Registry
end

Base.getproperty(entry::RegistryEntry, sym::Symbol) =
    getproperty(getfield(entry, :row), sym)

Base.propertynames(entry::RegistryEntry) =
    propertynames(getfield(entry, :row))

Base.NamedTuple(entry::RegistryEntry) = NamedTuple(getfield(entry, :row))

Base.keys(reg::Registry) = keys(getfield(reg, :index))

"""
    load(entry)
    load(registry[id])

Load a feature represented by a registry entry. The behavior depends on
the `Registry`.
"""
load(entry::RegistryEntry; kwargs...) = getfield(entry, :registry).loadfn(getfield(entry, :row); kwargs...)

getdata(registry::Registry) = getfield(registry, :data)
getfields(registry::Registry) = getfield(registry, :fields)

"""
    Registry(fields; kwargs...)

Create a feature registry with columns of [`Field`](#)s.

## Keyword arguments

- `name::String`: A descriptive name.
- `fields::NamedTuple`:
- `description::String = ""`: Description text for the registry. Shown when `info(registry)`
    is called.
- `loadfn = identity`: Function to apply over a row when [`load`](#) is called. For example,
    calling `load(registry["id"])` will call `loadfn`.
"""
function Registry(fields::NamedTuple; name = "", loadfn = identity, description = "", id = (:id,))
    id = id isa Symbol ? (id,) : id
    data = StructArray(NamedTuple(key => field.createfn() for (key, field) in pairs(fields)))
    index = Dict{Any, Int}()

    for id_ in id
        hasproperty(fields, id_) || throw(ArgumentError("ID field $id_ does not exist!"))
    end

    return Registry(fields, id, data, index, name, description, loadfn)
end

function makerow(fields, row::NamedTuple)
    NamedTuple(
        key => field.transformfn(field.computefn(row, key))
        for (key, field) in pairs(fields)
    )
end

struct DuplicateIDError <: Exception
    key
end

Base.showerror(io::IO, e::DuplicateIDError) = print(
    io, "DuplicateIDError: ID `", e.key, "` already registered!")

function Base.push!(registry::Registry, entry::NamedTuple)
    row = makerow(getfield(registry, :fields), entry)

    # check unique ID
    id = Tuple(row[id] for id in getfield(registry, :id))
    haskey(registry, id) && throw(DuplicateIDError(id))

    # enter row data and update index
    push!(getfield(registry, :data), row)
    index = getfield(registry, :index)
    index[id] = length(getfield(registry, :data))

    return row
end

# ### `Base` interface for `Registry`

function Base.getindex(registry::Registry, idx)
    return RegistryEntry(
        parent(getfield(registry, :data))[getfield(registry, :index)[_index_key(registry, idx)]],
        registry)
end

function Base.getindex(registry::Registry, ::Colon, sym::Symbol)
    return getproperty(getdata(registry), sym)
end

Base.length(registry::Registry) = length(registry.data)

Base.haskey(registry::Registry, idx) =
    haskey(getfield(registry, :index), _index_key(registry, idx))

Base.sort(registry::Registry, col::Symbol; rev = false) =
    withdata(
        registry,
        view(getdata(registry), sortperm(getproperty(getdata(registry), col); rev)))

Base.view(registry::Registry, inds...) = withdata(registry, view(getdata(registry), inds...))

_index_key(registry::Registry, nt::NamedTuple) =
    Tuple(nt[i] for i in getfield(registry, :id))

function _index_key(registry::Registry, id)
    length(getfield(registry, :id)) == 1 || throw(KeyError("Registry has a multi-column ID `$(getfield(registry, :id))`, please pass in a named tuple."))
    return (id,)
end

_index_key(::Registry, t::Tuple) = t

function withdata(registry::Registry, data)
    reg = Setfield.@set registry.data = data

    # recreate index
    index = Dict{Any, Int}()
    for (i, row) in enumerate(data)
        id = Tuple(row[id] for id in getfield(registry, :id))
        index[id] = i
    end
    reg = Setfield.@set reg.index = index

    return reg
end

# ### Tests for `Registry`

@testset "Registry" begin
    testregistry() = Registry((
        id = Field(String, name = "ID"),
    ), name = "Test registry")
    @test_nowarn testregistry()

    @testset "push!" begin
        registry = testregistry()
        @test_nowarn push!(registry, (; id = "myid"))
    end

    @testset "push! duplicate" begin
        registry = testregistry()
        push!(registry, (; id = "myid"))
        @test_throws DuplicateIDError push!(registry, (; id = "myid"))
    end

    @testset "getindex" begin
        @testset "get row" begin
            registry = testregistry()
            push!(registry, (; id = "myid"))
            @test_nowarn registry["myid"]
            @test NamedTuple(registry["myid"]) == (; id = "myid")
            @test_nowarn registry[("myid",)]
            @test_nowarn registry[(id = "myid",)]
        end

        @testset "get column" begin
            registry = testregistry()
            push!(registry, (; id = "myid"))
            @test_nowarn registry[:, :id]
            @test registry[:, :id] == ["myid"]
        end
    end

    @testset "length" begin
        registry = testregistry()
        @test length(registry) == 0
        push!(registry, (; id = "myid"))
        @test length(registry) == 1
    end
end


text:
using Random

state = [1,2,3,4,5,6,7,8,0]

function new_game()
    shuffle!(state)
end

function display()
    for j = 1:3
        for i = 1:3
            c = state[(j-1)*3+i]
            if c == 0
                print("  ")
            else
                print(c, " ")
            end
        end
        println()
    end
end

function user_input()
    print("w|a|s|d >> ")
    user = lowercase(strip(readline(stdin)))
    if user == "w"
        idx = findfirst(isequal(0), state)
        if idx < 7
            state[idx] = state[idx+3]
            state[idx+3] = 0
        end
    elseif user == "a"
        idx = findfirst(isequal(0), state)
        if !isnothing(findfirst(isequal(idx), [1,2,4,5,7,8]))
            state[idx] = state[idx+1]
            state[idx+1] = 0
        end
    elseif user == "s"
        idx = findfirst(isequal(0), state)
        if idx > 3
            state[idx] = state[idx-3]
            state[idx-3] = 0
        end
    elseif user == "d"
        idx = findfirst(isequal(0), state)
        if !isnothing(findfirst(isequal(idx), [2,3,5,6,8,9]))
            state[idx] = state[idx-1]
            state[idx-1] = 0
        end
    elseif user == "new"
        new_game()
    elseif user == "end"
        return false
    end

    if state == [1,2,3,4,5,6,7,8,0]
        println("You won!")
        display()
        return false
    end

    return true
end

new_game()
display()
while user_input()
    display()
end

max_stars_repo_name: aole/Games-in-Julia

max_issues_repo_name: aole/Games-in-Julia

max_forks_repo_name: aole/Games-in-Julia


text:
# t3r3
"""
    Fixed(body1::AbstractBody, body2; p1, p2, qoffset)

A fixed connection between two bodies.
"""
Fixed(body1::AbstractBody{T}, body2; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T})) where T =
    Translational3{T}(body1, body2; p1, p2), Rotational3{T}(body1, body2; qoffset)

# t2r3
"""
    Prismatic(body1, body2, axis; p1, p2, qoffset, spring, damper)

A prismatic joint between two bodies.
"""
Prismatic(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T}), spring = zero(T), damper = zero(T)) where T =
    Translational2{T}(body1, body2; p1, p2, axis, spring, damper), Rotational3{T}(body1, body2; qoffset, spring, damper)

# t1r3
"""
    Planar(body1, body2, axis; p1, p2, qoffset, spring, damper)

A planar joint between two bodies.
"""
Planar(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T}), spring = zero(T), damper = zero(T)) where T =
    Translational1{T}(body1, body2; p1, p2, axis, spring, damper), Rotational3{T}(body1, body2; qoffset, spring, damper)

# t0r3
"""
    FixedOrientation(body1, body2; qoffset, spring, damper)

Fixed orientation between two bodies (chicken's head).
"""
FixedOrientation(body1::AbstractBody{T}, body2; qoffset = one(UnitQuaternion{T}), spring = zero(T), damper = zero(T)) where T =
    Translational0{T}(body1, body2; spring, damper), Rotational3{T}(body1, body2; qoffset, spring, damper)

# t3r2
"""
    Revolute(body1, body2, axis; p1, p2, qoffset, spring, damper)

A revolute joint between two bodies (pin, continuous, hinge joint).
"""
Revolute(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T}), spring = zero(T), damper = zero(T)) where T =
    Translational3{T}(body1, body2; p1, p2, spring, damper), Rotational2{T}(body1, body2; axis, qoffset, spring, damper)

# t2r2
"""
    Cylindrical(body1, body2, axis; p1, p2, qoffset, spring, damper)

A cylindrical joint between two bodies.
"""
Cylindrical(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T}), spring = zero(T), damper = zero(T)) where T =
    Translational2{T}(body1, body2; p1, p2, axis, spring, damper), Rotational2{T}(body1, body2; axis, qoffset, spring, damper)

# t1r2
"""
    PlanarAxis(body1, body2, axis; p1, p2, qoffset, spring, damper)

A planar joint between two bodies with a rotation axis perpendicular to the plane (turtle bot).
"""
PlanarAxis(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T}), spring = zero(T), damper = zero(T)) where T =
    Translational1{T}(body1, body2; p1, p2, axis, spring, damper), Rotational2{T}(body1, body2; axis, qoffset, spring, damper)

# t0r2: something like a gyro?
# ?

# What is r1?
# t3r1: universal?
# ?
# t2r1: ?
# ?
# t1r1: ?
# ?
# t0r1: ?
# ?

# t3r0
"""
    Spherical(body1, body2; p1, p2, spring, damper)

A spherical joint between two bodies (ball-and-socket joint).
"""
Spherical(body1::AbstractBody{T}, body2; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T}), spring = zero(T), damper = zero(T)) where T =
    Translational3{T}(body1, body2; p1, p2, spring, damper), Rotational0{T}(body1, body2; qoffset, spring, damper)

# t2r0
"""
    CylindricalFree(body1, body2, axis; p1, p2, spring, damper)

A cylindrical joint between two bodies with unconstrained orientation (point-on-line).
"""
CylindricalFree(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), spring = zero(T), damper = zero(T)) where T =
    Translational2{T}(body1, body2; p1, p2, axis, spring, damper), Rotational0{T}(body1, body2; spring, damper)

# t1r0
"""
    PlanarFree(body1, body2, axis; p1, p2, spring, damper)

A planar joint between two bodies with unconstrained orientation.
"""
PlanarFree(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), spring = zero(T), damper = zero(T)) where T =
    Translational1{T}(body1, body2; p1, p2, axis, spring, damper), Rotational0{T}(body1, body2; spring, damper)

# t0r0
"""
    Floating(body1, body2; spring, damper)

An unconstrained connection between two bodies (connection between floating base and origin).
"""
Floating(body1::AbstractBody{T}, body2; spring = zero(T), damper = zero(T)) where T =
    Translational0{T}(body1, body2; spring, damper), Rotational0{T}(body1, body2; spring, damper)


text:
##  Modal analysis of the L-frame of Argyris.
# Eigenvector 1 frequency 11.2732  [Hz]
#   Eigenvector 2 frequency 30.5269  [Hz]
module argyr_l_frame_modal_examples

using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexStructures.FEMMCorotBeamModule
using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
mass = FEMMCorotBeamModule.mass
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using VisualStructures: plot_points, plot_nodes, plot_midline, render, plot_space_box, plot_solid, space_aspectratio, default_layout_3d
using PlotlyJS
using JSON

function argyr_l_frame_modal()
    # Parameters:
    E=71240.0;#MPa
    nu=0.31;# Poisson ratio
    rho=5e-9;
    b=3.0; h=30.0; L=240.0; # cross-sectional dimensions and length of each leg in millimeters
    # Choose the mass formulation:
    mass_type=2;

    # Reference frequencies
    reffs = [11.2732, 30.5269]
    neigvs = 2;

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [0.0, 1.0, 0.0])

    # Select the number of elements per leg.
    n=8;
    members = []
    push!(members, frame_member([0 0 L; L 0 L], n, cs))
    push!(members, frame_member([L 0 L; L 0 0], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    l1 = selectnode(fens; box = [0 0 0 0 L L], tolerance = L / 10000)
    for i in [1, 2, 3, 4, 5, 6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the global discrete system
    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);
    M = mass(femm, geom0, u0, Rfield0, dchi);

    # Solve the eigenvalue problem
    d,v,nconv = eigs(K, M; nev=2*neigvs, which=:SM, explicittransform=:none)
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Natural frequencies: $fs [Hz]")
    println("Reference: $reffs [Hz]")

    # Visualize vibration modes
    scattersysvec!(dchi, v[:, 1])
    update_rotation_field!(Rfield0, dchi)
    tbox = plot_space_box([[-L/2 -L/2 0]; [L/2 L/2 L]])
    plots = cat(tbox, plot_nodes(fens),
        plot_solid(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
        dims = 1)
    layout = default_layout_3d()
    layout[:scene][:aspectratio] = space_aspectratio(fens.xyz)
    layout[:scene][:aspectmode] = "manual"
    render(plots; layout = layout)

    return true
end # argyr_l_frame_modal

function argyr_l_frame_modal_anim()
    # Parameters:
    E=71240.0;#MPa
    nu=0.31;# Poisson ratio
    rho=5e-9;
    b=3.0; h=30.0; L=240.0; # cross-sectional dimensions and length of each leg in millimeters
    # Choose the mass formulation:
    mass_type=2;

    # Reference frequencies
    reffs = [11.2732, 30.5269]
    neigvs = 2;

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [0.0, 1.0, 0.0])

    # Select the number of elements per leg.
    n = 4;
    members = []
    push!(members, frame_member([0 0 L; L 0 L], n, cs))
    push!(members, frame_member([L 0 L; L 0 0], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    l1 = selectnode(fens; box = [0 0 0 0 L L], tolerance = L/10000)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the global discrete system
    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);
    M = mass(femm, geom0, u0, Rfield0, dchi);

    # Solve the eigenvalue problem
    d,v,nconv = eigs(K, M; nev=2*neigvs, which=:SM, explicittransform=:none)
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Natural frequencies: $fs [Hz]")
    println("Reference: $reffs [Hz]")

    # Visualize vibration modes
    mode = 2
    tbox = plot_space_box([[-L/2 -L/2 0]; [L/2 L/2 1.1*L]])
    tenv0 = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
    plots = cat(tbox, tenv0; dims = 1)
    layout = default_layout_3d()
    layout[:scene][:camera][:eye] = Dict(:x=>1.02, :y=> 1.02, :z=> 0.836)
    layout[:scene][:camera][:center] = Dict(:x=>0.058, :y=>0.065, :z=>-0.122)
    pl = render(plots; layout = layout, title = "Mode $(mode)")
    Rfield1 = deepcopy(Rfield0)
    scale = L/3  /  maximum(v[:, mode])
    for xscale in scale.*sin.(collect(0:1:72).*(2*pi/21))
        scattersysvec!(dchi, xscale.*v[:, mode])
        Rfield1 = deepcopy(Rfield0)
        update_rotation_field!(Rfield1, dchi)
        tenv1 = plot_solid(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield1.values, facecolor = "rgb(25, 255, 25)");
        plots = cat(tbox, tenv0, tenv1; dims = 1)
        react!(pl, plots, pl.plot.layout)
        sleep(0.115)
    end
    savejson(pl, "plots.json")

    return true
end # argyr_l_frame_modal_anim

function allrun()
    println("#####################################################")
    println("# argyr_l_frame_modal ")
    argyr_l_frame_modal()
    println("#####################################################")
    println("# argyr_l_frame_modal_anim ")
    argyr_l_frame_modal_anim()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
using .SampleChainsDynamicHMC
using Random

export sample

using ..Soss

using SampleChainsDynamicHMC:DynamicHMCConfig

function sample(rng::AbstractRNG,
    m::ConditionalModel,
    config::DynamicHMCConfig,
    nsamples::Int=1000,
    nchains::Int=4)

    ℓ(x) = Soss.logdensity(m, x)
    tr = xform(m)

    chains = newchain(rng, nchains, config, ℓ, tr)
    sample!(chains, nsamples - 1)
    return chains
end

function sample(
    m::ConditionalModel,
    config::DynamicHMCConfig,
    nsamples::Int=1000,
    nchains::Int=4)

    sample(Random.GLOBAL_RNG, m, config, nsamples, nchains)
end


text:
# This module solves the Ramsey problem for the quasi linear case
# We'll try to replicate figures from quant econ https://python-advanced.quantecon.org/opt_tax_recur.html#top

#----------------- CHECKING WITH QUANTECON ---------------------------------------------------#
## Parameters
β = 0.9;
σ = 2;
γ = 2;

## Utility functions
u(x) = x^(1 - σ) / (1 - σ);
v(x) = x^(1 + γ) / (1 + γ);
U(c, n) = u(c) - v(n);

## Differentials
Uc(c) = ForwardDiff.derivative(u, c);
Ucc(c) = ForwardDiff.derivative(Uc, c);

Un(n) = -ForwardDiff.derivative(v, n);
Unn(n) = ForwardDiff.derivative(Un, n);

## Defining states
g = [0.1; 0.1; 0.1; 0.1; 0.2; 0.1];
θ = ones(length(g));
s = [(1.0, g[i]) for i in 1:6];
Ps = [0.0 1.0 0.0 0.0 0.0 0.0;
      0.0 0.0 1.0 0.0 0.0 0.0;
      0.0 0.0 0.0 0.5 0.5 0.0;
      0.0 0.0 0.0 0.0 0.0 1.0;
      0.0 0.0 0.0 0.0 0.0 1.0;
      0.0 0.0 0.0 0.0 0.0 1.0];

include("NewtonRoot.jl");

## This function calculates Ramsey Allocation, given Φ and initial state s0
function RamseyAllocation(Φ::Real, s0::Tuple)

    function time1_allocation(s::Tuple)
        θ, g = s[1], s[2];
        n(c::Real) = (g + c) / θ;
        res1(c::Real) = (1 + Φ) * (Uc(c) * θ + Un(n(c))) + Φ * (c * Ucc(c) * θ + n(c) * Unn(n(c)));
        c = NewtonRoot(c -> res1(c), 0.5);

        return c, n(c)
    end

    c_time1 = [time1_allocation(s[i])[1] for i in 1:length(s)];
    n_time1 = [time1_allocation(s[i])[2] for i in 1:length(s)];

    function time0_allocation(s::Tuple)
        θ, g = s[1], s[2];
        n(c::Real) = (g + c) / θ;
        b0(c::Real) = 1.0;
        res0(c::Real) = (1 + Φ) * (Uc(c) * θ + Un(n(c))) + Φ * (c * Ucc(c) * θ + n(c) * Unn(n(c))) - Φ * Ucc(c) * b0(c) * θ
        c = NewtonRoot(c -> res0(c), 0.5)

        return c, n(c)
    end

    c_time0, n_time0 = time0_allocation(s0)[1], time0_allocation(s0)[2];

    return c_time1, n_time1, c_time0, n_time0

end

c_time1, n_time1, c_time0, n_time0 = RamseyAllocation(0.06175628494006816, s[1])

sum_mat = inv(1.0I - β * Ps) * (Uc.(c_time1) .* c_time1 .+ Un.(n_time1) .* n_time1);

ResImC = Uc(c_time0) - Uc(c_time0) * c_time0 - Un(n_time0) * n_time0 - β * sum([Ps[1,i]*sum_mat[i] for i in 1:length(sum_mat)]);

b_time1 = sum_mat ./ Uc.(c_time1);
b_time0 = 1.0;

τ_time1 = [1.0 + Un(n_time1[i]) / Uc(c_time1[i]) for i in 1:length(s)];

function search_Φ(Φ::Real, s0::Tuple)
    c_time1, n_time1, c_time0, n_time0 = RamseyAllocation(Φ, s0);

    i_s0 = findfirst(isequal(s0), s);

    sum_mat = inv(1.0I - β * Ps) * (Uc.(c_time1) .* c_time1 .+ Un.(n_time1) .* n_time1); #this sum is independent of s0

    implementCons = Uc(c_time0) - Uc(c_time0) * c_time0 - Un(n_time0) * n_time0 - β * sum([Ps[i_s0,i] * sum_mat[i] for i in 1:length(sum_mat)]);

    return implementCons

end

Φ_est = NewtonRoot(x -> search_Φ(x, s[1]), 0.1); # 0.06175628494006816

τ_time0 = 0.0;

## Plotting
sHist_h = [2, 3, 4, 6, 6, 6];
sHist_l = [2, 3, 5, 6, 6, 6];

c_sim_h = [c_time1[i] for i in sHist_h];
c_sim_l = [c_time1[i] for i in sHist_l];

n_sim_h = [n_time1[i] for i in sHist_h];
n_sim_l = [n_time1[i] for i in sHist_l];

b_sim_h = [b_time1[i] for i in sHist_h];
b_sim_l = [b_time1[i] for i in sHist_l];

τ_sim_h = [τ_time1[i] for i in sHist_h];
τ_sim_l = [τ_time1[i] for i in sHist_l];

c_plot_h = vcat(c_time0, c_sim_h)
c_plot_l = vcat(c_time0, c_sim_l)

n_plot_h = vcat(n_time0, n_sim_h)
n_plot_l = vcat(n_time0, n_sim_l)

b_plot_h = vcat(b_time0, b_sim_h)
b_plot_l = vcat(b_time0, b_sim_l)

τ_plot_h = vcat(τ_time0, τ_sim_h)
τ_plot_l = vcat(τ_time0, τ_sim_l)

plot(c_plot_h, label="g high")
plot!(c_plot_l, label = "g low")

plot(n_plot_h, label="g high")
plot!(n_plot_l, label = "g low")

plot(b_plot_l, label="g high")
plot!(b_plot_h, label = "g low")

plot(τ_plot_l, label="g high")
plot!(τ_plot_h, label = "g low")


text:
using GLMakie # openGL backend, uses GPu, very interactive, runs outside browser
# using WGLMakie # WebGL backend, interactive, runs inside browser
# using Cairo # 2d plots
using Makie.MakieLayout

# Plots latitudes and longitude of sphere.
function plotLatLon3d!(
        scene :: Union{Scene,LScene};
        lat_num :: Int64 = 10, # number of latitudes
        lon_num :: Int64 = 20, # number of longitude
        lat_den :: Int64 = 100, # number of points per latitude, den for density
        lon_den :: Int64 = 50 # number of points per longitude
) :: Nothing
    lats = [ [ sphericalToCartesian( SphericalCoord(theta, phi) ) for phi in LinRange(0, 2π, lat_den) ] for theta in LinRange(0, π, lat_num)]
    lons = [ [ sphericalToCartesian( SphericalCoord(theta, phi) ) for theta in LinRange(0, π, lon_den)] for phi in LinRange(0, 2π, lon_num)]

    # plot each line separately because plotting library connects consecutive ponts in the array in second argument using a line
    for l=lons
       lines!(scene, [l[j].x for j in 1:lon_den], [l[j].y for j in 1:lon_den], [l[j].z for j in 1:lon_den])
    end

    for l=lats
       lines!(scene, [l[j].x for j in 1:lat_den], [l[j].y for j in 1:lat_den], [l[j].z for j in 1:lat_den])
    end
    nothing
end

# Plots latitudes and longitude of sphere.
function plotLatLon2d!(
        scene :: Union{Scene,LScene};
        lat_num :: Int64 = 10, # number of latitudes
        lon_num :: Int64 = 20, # number of longitude
        lat_den :: Int64 = 100, # number of points per latitude, den for density
        lon_den :: Int64 = 50 # number of points per longitude
) :: Nothing
    lats = [ [ SphericalCoord(theta, phi) for phi in LinRange(0, 2π, lat_den) ] for theta in LinRange(0, π, lat_num)]
    lons = [ [ SphericalCoord(theta, phi) for theta in LinRange(0, π, lon_den)] for phi in LinRange(0, 2π, lon_num)]

    # plot each line separately because plotting library connects consecutive ponts in the array in second argument using a line
    for l=lons
       lines!(scene, [l[j].phi for j in 1:lon_den], [l[j].theta for j in 1:lon_den])
    end

    for l=lats
       lines!(scene, [l[j].phi for j in 1:lat_den], [l[j].theta for j in 1:lat_den])
    end
    nothing
end

# convert field specified in spherical coords to normal cartesian coordinates
function cartesianField( sphericalField )
    field = (x,y,z) ->  sphericalToCartesianVelocity(
        cartesianToSpherical( CartesianCoord( x,y,z ) ),
        sphericalField( cartesianToSpherical( CartesianCoord( x,y,z ) ) )
    )
    return( field )
end

function plotVelocityField3d!( threeDPlot::Union{Scene,LScene}, fieldGenerator, objects::Vector{<:Object}, lat_den::Int64, lon_den::Int64 )::Nothing
    sphericalField = @lift( p -> fieldGenerator( p, $objects ) )
    ∇f = @lift( cartesianField( $sphericalField ) )

    # create arrows plot for velocity vector field on a lot long grid
    gx = reduce(vcat,[[cospi(φ)*sinpi(θ) for φ in LinRange(0, 2π, round( Int, lat_den*sin(θ)))] for θ in LinRange(0, π, lon_den) ])
    gy = reduce(vcat,[[sinpi(φ)*sinpi(θ) for φ in LinRange(0, 2π, round( Int, lat_den*sin(θ)))] for θ in LinRange(0, π, lon_den) ])
    gz = reduce(vcat,[[cospi(θ)          for φ in LinRange(0, 2π, round( Int, lat_den*sin(θ)))] for θ in LinRange(0, π, lon_den) ])

    # we normalize the vectors because vectors can be too large at some places messing the entire plot
    n∇fxyz = @lift( map((x,y,z) -> normalize($∇f(x,y,z)), gx,gy,gz) )
    vx = @lift( map(v->v.x, $n∇fxyz) )
    vy = @lift( map(v->v.y, $n∇fxyz) )
    vz = @lift( map(v->v.z, $n∇fxyz) )
    arrows!(threeDPlot, gx, gy, gz, vx, vy, vz, arrowsize = 0.01, lengthscale = 0.05 )

    φs = LinRange(0, 2π, lat_den)
    θs = LinRange(0, π, lon_den)

    # plot the magnitude of velocity vector field on spherical surface
    sphereGrid = [ sphericalToCartesian( SphericalCoord(theta, phi) ) for phi in φs, theta in θs ]
    # TODO: we are recomputing ∇f on the sphere grid again, we already did that for a 'similar' grid when making arrows, can we reuse those computations
    colorsGrid = @lift( map( p->log(1 + norm($∇f(p.x,p.y,p.z))), sphereGrid ) )
    surface!(threeDPlot, map(s->s.x,sphereGrid), map(s->s.y,sphereGrid), map(s->s.z,sphereGrid), color=colorsGrid, colormap=:dense )
    nothing
end

function plotVelocityField2d!( twoDPlot::Union{Scene,LScene}, fieldGenerator, objects::Vector{<:Object}, lat_den::Int64, lon_den::Int64 )::Nothing
    # now we plot streams on theta-phi chart
    # phi is x axis and theta is y axis
    sphereCoordToPoint2f(vector::SphericalCoord, position::SphericalCoord ) = Point2f0( csc( position.theta)*vector.phi, vector.theta )
    # velocity func in polar coordinates
    # this gives (dϕ/dt, dθ/dt)
    ∇f_polar = @lift( (phi,theta) -> sphereCoordToPoint2f(
        fieldGenerator(
            SphericalCoord(theta,phi),
            $objects
        ),
        SphericalCoord(theta, phi)
    ) )

    φs = LinRange(0, 2π, lat_den)
    θs = LinRange(0, π, lon_den)

    # plot streams on theta phi chart
    streamplot!( twoDPlot, ∇f_polar, φs, θs, arrow_size=0.1 ) # streamplot is computationally expensive

    # plot vector magnitude
    # we don't need to use spherical coordinate norm as the spherical vector already include sin(theta) in e_phi component
    # TODO: we are recomputing ∇f on the sphere grid again, we already did this when making density surface plot on 3d chart, we can reuse the computations
    # TODO: in here vector norm should be on (sin(θ) dϕ/dt, dθ/dt), rather than (dϕ/dt, dθ/dt)
    colorsGridTwoD = @lift( [ log(1+norm($∇f_polar(phi,theta))) for phi in φs, theta in θs ] )
    heatmap!(twoDPlot, φs, θs, colorsGridTwoD, colormap=:dense)
    nothing
end

# Create a slider animation.
function animate(
    simObjects :: Vector{<:Vector{<:Object}}; #array of all time simuatled array of objects
    plotLatLon :: Bool = true,
    fieldGenerator = nothing, # a (point::SphericalCoord, objects::Vector{Object}} -> vector::SphericalCoord) function, if want to plot arrows and density on sphere
    resolution :: Tuple{Int64,Int64} = (1000,600),
    makeMovie :: Bool = true,
    cameraPosition :: Vec{3,Float32} = Vec3f0(-2, 0, 0),
    cameraLookingAt :: Vec{3,Float32} = Vec3f0(-1, 0, 0),
    frameRate ::Int64 = 1,
    movieName ::String = "animation.mp4",
) :: Nothing
    numFrames = length(simObjects)
    scene, layout = layoutscene(resolution = resolution) # create a new scene

    # a slider to move forward or backward in animation
    frameSlider = layout[1,:] = Slider(scene, range = LinRange(1, numFrames, numFrames), startvalue = 1)

    threeDPlot = layout[2,1] = LScene(scene) # full interactive shphere plot
    twoDPlot = layout[2,2] = LScene(scene) # theta phi chart

    cam3d!(threeDPlot.scene)
    update_cam!( threeDPlot.scene, cameraPosition, cameraLookingAt )

    # create an observable state which the plots will 'watch'
    # later we will simply update the state and scene will be automatically updated
    # objectState is dependent on slider so we 'lift' it to an observable
    sliderVal = frameSlider.value
    simObjectstate = @lift( simObjects[round( Int, $sliderVal)] )

    lat_den = 100 # number of points per latitude (= number of longitudes for grid below)
    lon_den = 50 # number of points per longitude (=number of longitude for grid below)

    if plotLatLon
        plotLatLon3d!(threeDPlot, lat_den = lat_den, lon_den = lon_den) # plot latitude and longitude
        plotLatLon2d!(twoDPlot, lat_den = lat_den, lon_den = lon_den)
    end

    # plot arrows and density for the fields if the fieldGenerator is provided
    if fieldGenerator !== nothing
        plotVelocityField3d!( threeDPlot, fieldGenerator, simObjectstate, lat_den, lon_den )
        plotVelocityField2d!( twoDPlot, fieldGenerator, simObjectstate, lat_den, lon_den )
    else
        φs = LinRange(0, 2π, lat_den)
        θs = LinRange(0, π, lon_den)
        sphereGrid = [ sphericalToCartesian( SphericalCoord(theta, phi) ) for phi in φs, theta in θs ]
        colorsGrid =  map( p->2, sphereGrid )
        surface!(threeDPlot, map(s->s.x,sphereGrid), map(s->s.y,sphereGrid), map(s->s.z,sphereGrid), color=colorsGrid, colormap=:dense )
    end

    # plot all simObjects
    for i in 1:length(simObjects[1])
        o = @lift( $simObjectstate[i] )
        plot3d!( threeDPlot, o )
        plot2d!( twoDPlot, o )
    end

    if makeMovie
        # create animation: update the state obervable node by moving slider, record, loop
        record(scene, movieName, 1:numFrames; framerate = frameRate) do frameNum
            set_close_to!(frameSlider, frameNum)
        end
    end

    # create an interactive simulation controlled by slider
    RecordEvents(scene, "output")
    nothing
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations3D(1.4)

initial_condition = initial_condition_convergence_test
boundary_condition = BoundaryConditionDirichlet(initial_condition)

solver = DGSEM(polydeg=4, surface_flux=flux_hll)

# Note that half the elements will have left-handed coordinates.
# This allows for easier coupling. A version with only right-handed coordinates
# is used in meshes/p4est_mesh.jl.
function cubed_sphere_mapping(xi, eta, zeta, inner_radius, thickness, direction)
  alpha = xi * pi/4
  beta = eta * pi/4

  # Equiangular projection
  x = tan(alpha)
  y = tan(beta)

  # Coordinates on unit cube per direction
  cube_coordinates = [SVector(-1, x, y),
                      SVector( 1, x, y),
                      SVector(x, -1, y),
                      SVector(x,  1, y),
                      SVector(x, y, -1),
                      SVector(x, y,  1)]

  # Radius on cube surface
  r = sqrt(1 + x^2 + y^2)

  # Radius of the sphere
  R = inner_radius + thickness * (0.5 * (zeta + 1))

  # Projection onto the sphere
  R / r * cube_coordinates[direction]
end

function cubed_sphere_mapping(inner_radius, thickness, direction)
  mapping(xi, eta, zeta) = cubed_sphere_mapping(xi, eta, zeta, inner_radius, thickness, direction)
end

mapping_as_string_(inner_radius, thickness, direction) = """
function cubed_sphere_mapping(xi, eta, zeta, inner_radius, thickness, direction)
  alpha = xi * pi/4
  beta = eta * pi/4

  x = tan(alpha)
  y = tan(beta)

  r = sqrt(1 + x^2 + y^2)

  R = inner_radius + thickness * (0.5 * (zeta + 1))

  # Cube coordinates per direction
  cube_coordinates = [SVector(-1, x, y),
                      SVector( 1, x, y),
                      SVector(x, -1, y),
                      SVector(x,  1, y),
                      SVector(x, y, -1),
                      SVector(x, y,  1)]

  R / r * cube_coordinates[direction]
end; function cubed_sphere_mapping(inner_radius, thickness, direction)
  mapping(xi, eta, zeta) = cubed_sphere_mapping(xi, eta, zeta, inner_radius, thickness, direction)
end; mapping = cubed_sphere_mapping($inner_radius, $thickness, $direction)
"""

cells_per_dimension = (3, 3, 2)

inner_radius = 1
thickness = 0.1

  # Illustration of the local coordinates of each face. ξ and η are the first
  # local coordinates of each face. The third local coordinate ζ is always
  # pointing outwards, which allows for easy coupling.
  #               ┌────────────────────────────────────────────────────┐
  #              ╱│                                                   ╱│
  #             ╱ │                  ┌───> ξ                         ╱ │
  #            ╱  │                 ╱                               ╱  │
  #           ╱   │                V       4 (+y)                  ╱   │
  #          ╱    │               η                               ╱    │
  #         ╱     │                                              ╱     │
  #        ╱      │                                             ╱      │
  #       ╱       │                                            ╱       │
  #      ╱        │                                           ╱        │
  #     ╱         │                    η    5 (-z)           ╱         │
  #    ╱          │                    ↑                    ╱          │
  #   ╱           │                    │                   ╱           │
  #  ╱            │                    └───> ξ            ╱            │
  # ┌────────────────────────────────────────────────────┐    2 (+x)   │
  # │             │                                      │             │
  # │             │                                      │      ξ      │
  # │             │                                      │      ↑      │
  # │    1 (-x)   │                                      │      │      │
  # │             │                                      │      │      │
  # │      ξ      │                                      │     ╱       │
  # │      ↑      │                                      │    V        │
  # │      │      │                                      │   η         │
  # │      │      └──────────────────────────────────────│─────────────┘
  # │     ╱      ╱         η   6 (+z)                    │            ╱
  # │    V      ╱          ↑                             │           ╱
  # │   η      ╱           │                             │          ╱
  # │         ╱            └───> ξ                       │         ╱
  # │        ╱                                           │        ╱
  # │       ╱                                            │       ╱ Global coordinates:
  # │      ╱                                             │      ╱        y
  # │     ╱                      ┌───> ξ                 │     ╱         ↑
  # │    ╱                      ╱                        │    ╱          │
  # │   ╱                      V      3 (-y)             │   ╱           │
  # │  ╱                      η                          │  ╱            └─────> x
  # │ ╱                                                  │ ╱            ╱
  # │╱                                                   │╱            V
  # └────────────────────────────────────────────────────┘            z

mesh1 = StructuredMesh(cells_per_dimension, cubed_sphere_mapping(inner_radius, thickness, 1), periodicity=false,
                       mapping_as_string=mapping_as_string_(inner_radius, thickness, 1))

semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition, solver,
                                     source_terms=source_terms_convergence_test,
                                     boundary_conditions=(
                                       x_neg=BoundaryConditionCoupled(3, (:begin, :i_forward, :j_forward), Float64),
                                       x_pos=BoundaryConditionCoupled(4, (:begin, :i_forward, :j_forward), Float64),
                                       y_neg=BoundaryConditionCoupled(5, (:begin, :i_forward, :j_forward), Float64),
                                       y_pos=BoundaryConditionCoupled(6, (:begin, :i_forward, :j_forward), Float64),
                                       z_neg=boundary_condition,
                                       z_pos=boundary_condition,
                                     ))

mesh2 = StructuredMesh(cells_per_dimension, cubed_sphere_mapping(inner_radius, thickness, 2), periodicity=false,
                   mapping_as_string=mapping_as_string_(inner_radius, thickness, 2))

semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition, solver,
                                     source_terms=source_terms_convergence_test,
                                     boundary_conditions=(
                                       x_neg=BoundaryConditionCoupled(3, (:end, :i_forward, :j_forward), Float64),
                                       x_pos=BoundaryConditionCoupled(4, (:end, :i_forward, :j_forward), Float64),
                                       y_neg=BoundaryConditionCoupled(5, (:end, :i_forward, :j_forward), Float64),
                                       y_pos=BoundaryConditionCoupled(6, (:end, :i_forward, :j_forward), Float64),
                                       z_neg=boundary_condition,
                                       z_pos=boundary_condition,
                                     ))

mesh3 = StructuredMesh(cells_per_dimension, cubed_sphere_mapping(inner_radius, thickness, 3), periodicity=false,
                   mapping_as_string=mapping_as_string_(inner_radius, thickness, 3))

semi3 = SemidiscretizationHyperbolic(mesh3, equations, initial_condition, solver,
                                     source_terms=source_terms_convergence_test,
                                     boundary_conditions=(
                                       x_neg=BoundaryConditionCoupled(1, (:begin, :i_forward, :j_forward), Float64),
                                       x_pos=BoundaryConditionCoupled(2, (:begin, :i_forward, :j_forward), Float64),
                                       y_neg=BoundaryConditionCoupled(5, (:i_forward, :begin, :j_forward), Float64),
                                       y_pos=BoundaryConditionCoupled(6, (:i_forward, :begin, :j_forward), Float64),
                                       z_neg=boundary_condition,
                                       z_pos=boundary_condition,
                                     ))

mesh4 = StructuredMesh(cells_per_dimension, cubed_sphere_mapping(inner_radius, thickness, 4), periodicity=false,
                   mapping_as_string=mapping_as_string_(inner_radius, thickness, 4))

semi4 = SemidiscretizationHyperbolic(mesh4, equations, initial_condition, solver,
                                     source_terms=source_terms_convergence_test,
                                     boundary_conditions=(
                                       x_neg=BoundaryConditionCoupled(1, (:end, :i_forward, :j_forward), Float64),
                                       x_pos=BoundaryConditionCoupled(2, (:end, :i_forward, :j_forward), Float64),
                                       y_neg=BoundaryConditionCoupled(5, (:i_forward, :end, :j_forward), Float64),
                                       y_pos=BoundaryConditionCoupled(6, (:i_forward, :end, :j_forward), Float64),
                                       z_neg=boundary_condition,
                                       z_pos=boundary_condition,
                                     ))

mesh5 = StructuredMesh(cells_per_dimension, cubed_sphere_mapping(inner_radius, thickness, 5), periodicity=false,
                   mapping_as_string=mapping_as_string_(inner_radius, thickness, 5))

semi5 = SemidiscretizationHyperbolic(mesh5, equations, initial_condition, solver,
                                     source_terms=source_terms_convergence_test,
                                     boundary_conditions=(
                                       x_neg=BoundaryConditionCoupled(1, (:i_forward, :begin, :j_forward), Float64),
                                       x_pos=BoundaryConditionCoupled(2, (:i_forward, :begin, :j_forward), Float64),
                                       y_neg=BoundaryConditionCoupled(3, (:i_forward, :begin, :j_forward), Float64),
                                       y_pos=BoundaryConditionCoupled(4, (:i_forward, :begin, :j_forward), Float64),
                                       z_neg=boundary_condition,
                                       z_pos=boundary_condition,
                                     ))

mesh6 = StructuredMesh(cells_per_dimension, cubed_sphere_mapping(inner_radius, thickness, 6), periodicity=false,
                   mapping_as_string=mapping_as_string_(inner_radius, thickness, 6))

semi6 = SemidiscretizationHyperbolic(mesh6, equations, initial_condition, solver,
                                     source_terms=source_terms_convergence_test,
                                     boundary_conditions=(
                                       x_neg=BoundaryConditionCoupled(1, (:i_forward, :end, :j_forward), Float64),
                                       x_pos=BoundaryConditionCoupled(2, (:i_forward, :end, :j_forward), Float64),
                                       y_neg=BoundaryConditionCoupled(3, (:i_forward, :end, :j_forward), Float64),
                                       y_pos=BoundaryConditionCoupled(4, (:i_forward, :end, :j_forward), Float64),
                                       z_neg=boundary_condition,
                                       z_pos=boundary_condition,
                                     ))

semi = SemidiscretizationCoupled((semi1, semi2, semi3, semi4, semi5, semi6))

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.6)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
module btcj

using HTTP
using JSON

token = Dict(
    "apikey" => "YOUR_WML_APIKEY",
    "grant_type" => "urn:ibm:params:oauth:grant-type:apikey"
)

request_token = HTTP.request("POST", "https://iam.cloud.ibm.com/identity/token",
    ["Content-Type" => "application/x-www-form-urlencoded"],
    HTTP.URIs.escapeuri(token))

access_token = JSON.parse(String(request_token.body))

submission = Dict(
    "email_addr" => "YOUR_EMAIL_HERE",
    "wml_url" => "WML_URL",
    "iam_token" => access_token["access_token"],
    "submit_confirmation" => false
)

submit = HTTP.request("POST", "http://172.21.86.186:5000/submit",
    ["Content-Type" => "application/json"],
    JSON.json(submission))

println(String(submit.body))

end # module


text:
using ConvDiffMIPDECO
using jInv.Mesh
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.LinearSolvers
using Statistics

function runInv()
domain = [0 3. -0.5 1.5]
n      = [60 40]
M      = getRegularMesh(domain,n)
display(M)

ftrue = getFEMsource2D(M)
#pFor  = getConvDiffFEMParam(M,sig=.01,Ainv = getMUMPSsolver())
 pFor  = getConvDiffFEMParam(M,sig=.01,Ainv = getJuliaSolver())
dobs, = getData(vec(ftrue),pFor)
dobs += 0.4*randn(size(dobs))*mean(abs.(vec(dobs))) # add noise

sigback     = 0.0
Wt          = ones(prod(n.+1))
pMis        = getMisfitParam(pFor,Wt,dobs,SSDFun);

# configure regularization
alpha	   	= 1e-2;
reg         = (m,mr,M,I=1.0) -> wTVReg(m,mr,M,eps=1e-6)
mref        = zeros(M.nc);
# configure optimization
HesPrec     = getSSORRegularizationPreconditioner(1.0,1e-15,50);
cgit       	= 5;
pcgTol     	= 1e-1;
maxIter    	= 5;
minUpdate 	= 1e-3;

# use identity model
modFun      = identityMod
boundsLow  	= 0*ones(prod(n));
boundsHigh 	= 1*ones(prod(n));

# # use level set approach
modFun      = boundMod
boundsLow  	= -Inf*ones(prod(n));
boundsHigh 	= Inf *ones(prod(n));

maxStep		= 0.1*maximum(boundsHigh);

pInv = getInverseParam(M,modFun,reg,alpha,mref,
         boundsLow,boundsHigh,maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
         minUpdate=minUpdate, maxIter = maxIter,HesPrec=HesPrec);

mc,Dc,flag,His = projGNCG(0*ftrue[:],pInv,pMis);
fc = pInv.modelfun(mc)[1];

println("time grad misfit: $(sum(His.timeGradMisfit))")
println("time GN solve:    $(sum(His.timeLinSol))")
println("time misfit       $(sum(His.timeMisfit))")
println("time reg          $(sum(His.timeReg))")
end

runInv()

# using ProfileView
# Profile.clear()
#
# @profile runInv()
#
# ProfileView.view()

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
module MeshFactory

using LinearAlgebra

include("ellipseBE.jl")

export ellipseBE

end


text:
module BatchUtils

using ...Concepts

abstract type BatchingStrategy end
abstract type SequentialScan           <: BatchingStrategy end
abstract type SampleWithReplacement    <: BatchingStrategy end
abstract type SampleWithoutReplacement <: BatchingStrategy end

export BatchFactory,
  next,
  reset,
  has_next,
  initialize,
  BatchingStrategy,
  SequentialScan

mutable struct BatchFactory{T<:Any}
  n::Optional{Int64}
  batch_size::Optional{Int64}
  cur_batch::Optional{Int64}

  function BatchFactory{S}(; size::Int64) where S<:BatchingStrategy
    return new{S}(nothing, size, nothing)
  end

  function BatchFactory{S}(n, batch_size) where S<:BatchingStrategy
    return new{S}(n, batch_size, 1)
  end
end

function initialize(obj::BatchFactory{T}, x::Array{S, 1}) where {T<:BatchingStrategy, S<:Any}
  obj.n = length(x)
  obj.cur_batch = 1
end

function next(obj::BatchFactory{SequentialScan})

  start_pos = obj.batch_size * (obj.cur_batch - 1)  + 1
  end_pos   = min(start_pos + obj.batch_size - 1, obj.n)
  obj.cur_batch = obj.cur_batch + 1
  return start_pos:end_pos
end

function Base.reset(obj::BatchFactory{SequentialScan})
  obj.cur_batch = 1
end

function has_next(obj::BatchFactory{SequentialScan})
  return obj.batch_size * (obj.cur_batch - 1)  + 1 <= obj.n
end

end # end of BatchUtils


text:
using Luxor
using Colors
using QHull: chull
using ConvexBodyProximityQueries: closest_points
using StaticArrays

function randdisk(sz)
    r = rand(sz)
    θ = 2*π*rand(sz)
    x = sqrt.(r).*cos.(θ)
    y = sqrt.(r).*sin.(θ)
    return hcat(x, y)
end

function randomcvx(sz::Int)
    pts = randdisk(sz)
    ch = chull(pts)
    sz = length(ch.vertices)
    return randomcvx(ch.points[ch.vertices, :], Val(sz))
end

function randomcvx(hull::Array, ::Val{N}) where {N}
    hulln = vcat(hull', ones(1, N))
    return SMatrix{3, N}(hulln)
end

logo = Movie(400, 283, "test")

function backdrop(scene, framenumber)
    background("white")
end

transform2(θ,c) = SMatrix{2,3}([cos(θ) -sin(θ) c[1]; sin(θ) cos(θ) c[2]])
polyA = randomcvx(35)
polyB = randomcvx(35)
polyC = randomcvx(35)

function drawpoly(polygon, θ, center)
    p = transform2(θ, center)*polygon
    plux = [Point(p[1:2,i]...) for i in 1:size(p, 2)]
    poly(plux, :fill)
    return p
end

function frame(scene, framenumber)
    global polyA, polyB, polyC
    centerA = SMatrix{2, 1}( 0.,-√2.)
    centerB = SMatrix{2, 1}( 2., √2.)
    centerC = SMatrix{2, 1}(-2., √2.)
    scale(50)

    # draw polygons
    sethue(Colors.RGB(96/255,173/255,81/255))
    pA = drawpoly(polyA, deg2rad(-framenumber/2), centerA)
    sethue(Colors.RGB(170/255,121/255,193/255))
    pB = drawpoly(polyB, deg2rad(framenumber), centerB)
    sethue(Colors.RGB(213/255,99/255,92/255))
    pC = drawpoly(polyC, deg2rad(framenumber), centerC)

    # get closest points between objects
    cp = map(p -> closest_points(p[1], p[2], SVector{2}(1., 1.)), [[pA, pB], [pB, pC], [pC, pA]])
    cp = map(p -> map(c -> Point(c.data...), p), cp)

    # draw line segments
    sethue(Colors.RGB(102/255,130/225,223/255))
    map(p -> line(p[1], p[2], :stroke), cp)

    # draw points
    sethue("blue")
    map(p -> map(c -> circle(c, 0.04, :fill), p), cp)
end

animate(logo, [
    Scene(logo, backdrop, 0:719),
    Scene(logo, frame, 0:719)
    ],
    creategif=true,
    framerate=30,
    pathname="logo.gif")


text:
using DashBootstrapComponents, DashHtmlComponents

button_group = html_div(
    [
        dbc_radioitems(
            id="radios",
            className="btn-group",
            labelClassName="btn btn-secondary",
            labelCheckedClassName="active",
            options=[
                Dict("label" => "Option 1", "value" => 1),
                Dict("label" => "Option 2", "value" => 2),
                Dict("label" => "Option 3", "value" => 3),
            ],
            value=1,
        ),
        html_div(id="output"),
    ],
    className="radio-group",
);

callback!(app, Output("output", "children"), Input("radios", "value")) do value
    return "Selected value: $value"
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
"""    dynamics(S, k; Δτ = 1e-7, t = 1e7, n = 128, rtol = sqrt(eps()), brtol = eps())

Computes the dynamical properties of a liquid `S` with structure factor
approximated over a finite grid.  In particular it computes the intermediate
scattering function `F` and the self intermediate scattering function `Fs` at
wavenumber `k` as functions of the correlation time, up to a time `t` (`1e7` by
default).  It also computes the memory function `ζ` of the Generalized Langevin
Equation, which can later be used to estimate the diffusion coefficient and the
mean square displacement as functions of the correlation time.

The algorithm is based on the theoretical model known as Self-consistent
Generalized Langevin Equation.  The precision of the method can be controlled
by the keyword arguments `n`, `Δτ` and `rtol`, which correspond to the number
of points in the time grid per 'time decade', the initial time grid spacing,
and the convergence of the relative tolerance for the memory function `ζ`,
respectively.

This function also returns the long-time-limit mobility `b` of the system up to
a relative tolerance specified by `brtol`.
"""
function dynamics(S, k; t = 1e7, Δτ = 1e-7, n = 128, rtol = sqrt(eps()), brtol = eps())
    # Number of time points for which the short-times approximation is used.
    @assert n ≥ (n₀ = 8)
    # Dynamical variables, memory kernel variables and auxiliar variables.
    dvars, kvars, auxvars = initialize_dynamics(S, n)

    # The current function is just an initialization routine,
    # the real work starts here.
    return dynamics!(dvars, kvars, auxvars, S, k, t, Δτ, n₀, n, rtol, brtol)
end

function dynamics!(dvars, kvars, auxvars, S, k, t, Δτ, n₀, n, rtol, brtol)
    # Use a short-times approximation initially.
    approximate!(dvars, kvars.svars, auxvars.Z, Δτ, n₀)
    # Fill the first time grid employing the SCGGLE equations.
    solve!(dvars, kvars, auxvars, Δτ, n₀ + 1, n, rtol)
    # Allocate the output with the solution for the first-decade time grid.
    #
    # TODO: This could be improved moving the following line to `dynamics`,
    #       but the size should be carefully allocated.
    output = DynamicsOutput(dvars, S.grid, k, Δτ, n)

    # As long as the target time `t` is not reached, keep doubling the time
    # grid spacing `Δτ` and solve for each new time point.  Since we decimate
    # the previous grid keeping the first-half information, we only need to
    # compute the dynamics for the second half of each new time grid.
    n₀ = div(n, 2) + 1
    b′ = one(output.b)
    while 2n * Δτ < t
        Δτ *= 2
        b′ = output.b
        decimate!(dvars)
        solve!(dvars, kvars, auxvars, Δτ, n₀, n, rtol)
        update!(output, dvars, S.grid, k, Δτ, n₀, n)
    end

    # Reaching the time `t` won't generally guarantee that the asymptotic value
    # for the mobility `b` has converged.  If the asymptotic value of the
    # memory kernel is finite and greater than zero, we know that `b` should
    # the be zero, so we use this fact.  Otherwise, we keep iterating until
    # convergence of `b` up to relative precision `brtol`, but no longer store
    # the solutions to the SCGLE.
    avars, akvars, D₀ = initialize_asymptotics(S)
    g = FixedPoint(auxvars.Z, akvars, D₀)
    asymptotics!(g, avars, last(dvars.ζ); rtol = rtol)
    asymptotic_mobility!(output.b, b′, avars.ζ∞, dvars, kvars, auxvars,
                         Δτ, n₀, n, rtol, brtol)

    return output
end


text:
"""
Uses UK Ordnance Survey "A Guide to Coordinate Systems in Great Britain"
https://www.ordnancesurvey.co.uk/documents/resources/guide-coordinate-systems-great-britain.pdf

Worked example (Airy1830)
E = 651409.903 # meters
N = 313177.270 # meters

"""
function OSENtoLLA(E,N)
    # Airy 1830
    a = 6377563.396
    b = 6356256.909
    e2 = (a^2 - b^2)/a^2
    E0 = 400000
    N0 = -100000
    F0 = 0.9996012717
    ϕ0 = 49 * (π/180)
    λ0 = -2 * (π/180)

    ϕ1 = (N - N0)/a*F0 + ϕ0
    n=(a-b)/(a+b)
    M_function(ϕ) = b*F0*( (1 + n + 5/4*n^2 + 5/4*n^3) * (ϕ-ϕ0)
        - (3*n + 3*n^2 + 21/8*n^3) * sin(ϕ-ϕ0) * cos(ϕ+ϕ0)
        + (15/8*n^2 + 15/8*n^3) * sin(2*(ϕ-ϕ0)) * cos(2*(ϕ+ϕ0))
        - 35/24*n^3 * sin(3*(ϕ-ϕ0))*cos(3*(ϕ+ϕ0)))
    residual(M) = (N-N0-M)
    M = M_function(ϕ1)

    while residual(M) > 0.00000001
        ϕ1 = (residual(M)/a*F0) + ϕ1
        M = M_function(ϕ1)
    end

    ϕ1 / (π/180)

    ν = a * F0 * abs(1-e2*sin(ϕ1)^2)^(-0.5) # \nu
    ρ = a * F0 * (1-e2) * (1-e2*sin(ϕ1)^2)^(-1.5)
    η2 = ν/ρ-1

    VII = tan(ϕ1)/(2*ρ*ν)
    VIII = tan(ϕ1)/(24*ρ*ν^3) * (5 + 3*tan(ϕ1)^2 + η2 - 9*(tan(ϕ1)^2)*η2)
    IX = tan(ϕ1)/(720*ρ*ν^5) * (61 + 90*tan(ϕ1)^2 + 45*(tan(ϕ1)^4))
    X = sec(ϕ1)/ν
    XI = sec(ϕ1)/(6*ν^3) * (ν/ρ + 2*tan(ϕ1)^2)
    XII = sec(ϕ1)/(120*ν^5) * (5 + 28*tan(ϕ1)^2 + 24*tan(ϕ1)^4)
    XIIA = sec(ϕ1)/(5040*ν^7) * (61 + 662*tan(ϕ1)^2 + 1320*tan(ϕ1)^4 + 720*tan(ϕ1)^6)

    ϕ = ϕ1 - VII*(E-E0)^2 + VIII*(E-E0)^4 - IX*(E-E0)^6
    λ = λ0 + X*(E-E0) - XI*(E-E0)^3 + XII*(E-E0)^5 - XIIA*(E-E0)^7

    return ϕ/(π/180), λ/(π/180)
end


text:
using FANN
using Base.Test

# one dimensional problem with many weigths, so it will overfit
net = ANN([1, 10, 1], [:sigmoid_symmetric, :linear], b=0.1)

# N uniformly spaced point over [0, 1]
N = 10
X = reshape([0:N-1], 1, N)/(N-1)
y = X + 0.1*randn(1, N)

# define training and validation sets
srand(1)
idx = randperm(N)
tidx = idx[1:N/2]
vidx = idx[N/2+1:end]
tset = DataSet(X[:, tidx], y[:, tidx])
vset = DataSet(X[:, vidx], y[:, vidx])

# train
train!(net, tset, vset,
	   max_epochs=100,
	   epochs_between_reports=100,
	   desired_error=1e-8, minratio=0.95)

# how to test this behaviour? I should return a full
# history of the training.... but it is overkill now
# but i can test it and it stops