text:
mutable struct Green
    Tpair::Tuple{Int,Int}
    weight::Float64
    function Green(inT, outT)
        return new((inT, outT), 0.0)
    end
    function Green(tpair::Tuple{Int,Int})
        return new(tpair, 0.0)
    end
end

# add Tpairs to Green's function (in, out) or vertex4 (inL, outL, inR, outR)
function addTidx(obj, _Tidx)
    for (i, Tidx) in enumerate(obj.Tpair)
        if Tidx == _Tidx
            return i
        end
    end
    push!(obj.Tpair, _Tidx)
    push!(obj.weight, zero(eltype(obj.weight))) # add zero to the weight table of the object
    push!(obj.child, [])
    return length(obj.Tpair)
end

"""
    struct IdxMap{_Ver4}

    Map left vertex Tpair[lidx], right vertex Tpair[ridx], the shared Green's function G0 and the channel specific Green's function Gx to the top level 4-vertex Tpair[vidx]

# Arguments
- l::_Ver4 : left vertex
- r::_Ver4 : right vertex
- v::_Ver4 : composte vertex
- lidx::Int : left sub-vertex index
- ridx::Int : right sub-vertex index
- vidx::Int : composite vertex index
- G0::Green : shared Green's function index
- Gx::Green : channel specific Green's function index
"""
mutable struct IdxMap{_Ver4}
    l::_Ver4 #left vertex
    r::_Ver4 #right vertex
    v::_Ver4 #composte vertex
    lidx::Int # left sub-vertex index
    ridx::Int # right sub-vertex index
    vidx::Int # composite vertex index
    G0::Green # shared Green's function index
    Gx::Green # channel specific Green's function index
    node::Any # useful for constructing DiagTree
    function IdxMap(l::V, r::V, v::V, lidx, ridx, vidx, G0::Green, Gx::Green) where {V}
        return new{V}(l, r, v, lidx, ridx, vidx, G0, Gx, nothing)
    end
end

struct Bubble{_Ver4} # template Bubble to avoid mutually recursive struct
    id::Int
    chan::Int
    Lver::_Ver4
    Rver::_Ver4
    map::Vector{IdxMap{_Ver4}}

    function Bubble(ver4::_Ver4, chan::Int, oL::Int, level::Int, _id::Vector{Int}) where {_Ver4}
        # @assert chan in para.chan "$chan isn't a bubble channels!"
        @assert oL < ver4.loopNum "LVer loopNum must be smaller than the ver4 loopNum"

        idbub = _id[1] # id vector will be updated later, so store the current id as the bubble id
        _id[1] += 1

        oR = ver4.loopNum - 1 - oL # loopNum of the right vertex
        lLpidxOffset = ver4.loopidxOffset + 1
        rLpidxOffset = lLpidxOffset + oL
        LTidx = ver4.TidxOffset  # the first τ index of the left vertex
        TauNum = ver4.para.interactionTauNum # maximum tau number for each bare interaction
        RTidx = LTidx + (oL + 1) * TauNum   # the first τ index of the right sub-vertex

        if chan == T || chan == U
            LverChan = (level == 1) ? ver4.Fouter : ver4.F
            RverChan = (level == 1) ? ver4.Allouter : ver4.All
        elseif chan == S
            LverChan = (level == 1) ? ver4.Vouter : ver4.V
            RverChan = (level == 1) ? ver4.Allouter : ver4.All
        else
            error("chan $chan isn't implemented!")
        end

        # println("left ver chan: ", LsubVer, ", loop=", oL)
        # W = eltype(ver4.weight)
        Lver = _Ver4(ver4.para, LverChan, ver4.F, ver4.V, ver4.All;
            loopNum = oL, loopidxOffset = lLpidxOffset, tidxOffset = LTidx,
            level = level + 1, id = _id)

        Rver = _Ver4(ver4.para, RverChan, ver4.F, ver4.V, ver4.All;
            loopNum = oR, loopidxOffset = rLpidxOffset, tidxOffset = RTidx,
            level = level + 1, id = _id)

        @assert Lver.TidxOffset == ver4.TidxOffset "Lver Tidx must be equal to vertex4 Tidx! LoopNum: $(ver4.loopNum), LverLoopNum: $(Lver.loopNum), chan: $chan"

        ############## construct IdxMap ########################################
        map = []
        for (lt, LvT) in enumerate(Lver.Tpair)
            for (rt, RvT) in enumerate(Rver.Tpair)
                VerTidx = 0

                if chan == T
                    VerT = (LvT[INL], LvT[OUTL], RvT[INR], RvT[OUTR])
                    GTx = (RvT[OUTL], LvT[INR])
                elseif chan == U
                    VerT = (LvT[INL], RvT[OUTR], RvT[INR], LvT[OUTL])
                    GTx = (RvT[OUTL], LvT[INR])
                elseif chan == S
                    VerT = (LvT[INL], RvT[OUTL], LvT[INR], RvT[OUTR])
                    GTx = (LvT[OUTL], RvT[INR])
                else
                    throw("This channel is invalid!")
                end

                Gx = Green(GTx)
                push!(ver4.G[Int(chan)], Gx)

                GT0 = (LvT[OUTR], RvT[INL])
                G0 = Green(GT0)
                push!(ver4.G[1], G0)

                VerTidx = addTidx(ver4, VerT)
                for tpair in ver4.Tpair
                    @assert tpair[1] == ver4.TidxOffset "InL Tidx must be the same for all Tpairs in the vertex4"
                end

                ###### test if the internal + exteranl variables is equal to the total 8 variables of the left and right sub-vertices ############
                Total1 = vcat(collect(LvT), collect(RvT))
                Total2 = vcat(collect(GT0), collect(GTx), collect(VerT))
                # println(Total1)
                # println(Total2)
                @assert compare(Total1, Total2) "chan $(chan): G0=$GT0, Gx=$GTx, external=$VerT don't match with Lver4 $LvT and Rver4 $RvT"

                idxmap = IdxMap(Lver, Rver, ver4, lt, rt, VerTidx, G0, Gx)
                push!(map, idxmap)
                push!(ver4.child[VerTidx], idxmap)
            end
        end
        return new{_Ver4}(idbub, chan, Lver, Rver, map)
    end
end

"""
    function Ver4{W}(para::Para, loopNum = para.internalLoopNum, tidx = 1; chan = para.chan, F = para.F, V = para.V, level = 1, id = [1,]) where {W}

    Generate 4-vertex diagrams using Parquet Algorithm

#Arguments
- `para`: parameters. It should provide internalLoopNum, interactionTauNum, firstTauIdx
- `chan`: list of channels of the current 4-vertex.
- `F`   : channels of left sub-vertex for the particle-hole and particle-hole-exchange bubbles
- `V`   : channels of left sub-vertex for the particle-particle bubble
- `All`   : channels of right sub-vertex of all channels
- `Fouter`   : channels of left sub-vertex for the particle-hole and particle-hole-exchange bubbles, only take effect for the outermost bubble
- `Vouter`   : channels of left sub-vertex for the particle-particle bubble, only take effect for the outermost bubble
- `Allouter`   : channels of right sub-vertex of all channels
- `loopNum`: momentum loop degrees of freedom of the 4-vertex diagrams
- `tidx`: the first τ variable index. It will be the τ variable of the left incoming electron for all 4-vertex diagrams
- `level`: level in the diagram tree
- `id`: the first element will be used as the id of the Ver4. All nodes in the tree will be labeled in preorder depth-first search

#Remark:
- AbstractTrees interface is implemented for Ver4. So one can use the API in https://juliacollections.github.io/AbstractTrees.jl/stable/ to manipulate/print the tree structre of Ver4.
- There are three different methods to print/visualize the tree structre:
1) `print_tree(ver4::Ver4)` or `print_tree(bub::Bubble)` to print the tree to terminal. This function is provided by AbstractTrees API.
2) `newick(ver4::Ver4)` or `newick(bub::Bubble)` to serilize the tree to a newick format string. You may save the string to a text file, then visualize it with a newick format visualizer application.
3) `showTree(ver4::Ver4)` to visualize the tree using the python package ete3. You have to install ete3 properly to use this function.
"""
struct Ver4{W}
    para::Any
    chan::Vector{Int} # list of channels
    F::Vector{Int}
    V::Vector{Int}
    All::Vector{Int}
    Fouter::Vector{Int}
    Vouter::Vector{Int}
    Allouter::Vector{Int}

    ###### vertex topology information #####################
    id::Int
    level::Int

    #######  vertex properties   ###########################
    loopNum::Int
    loopidxOffset::Int # offset of the loop index from the first internal loop
    TidxOffset::Int # offset of the Tidx from the tau index of the left most incoming leg

    ######  components of vertex  ##########################
    G::SVector{16,Vector{Green}}  # large enough to host all Green's function
    bubble::Vector{Bubble{Ver4{W}}}

    ####### weight and tau table of the vertex  ###############
    Tpair::Vector{Tuple{Int,Int,Int,Int}}
    child::Vector{Vector{IdxMap{Ver4{W}}}}
    weight::Vector{W}

    function Ver4{W}(para, chan, F = [I, U, S], V = [I, T, U], All = union(F, V);
        loopNum = para.innerLoopNum, loopidxOffset = 0, tidxOffset = 0,
        Fouter = F, Vouter = V, Allouter = All,
        level = 1, id = [1,]
    ) where {W}

        @assert para.totalTauNum >= (loopNum + 1) * para.interactionTauNum "$para"

        if level > 1
            @assert Set(F) == Set(Fouter)
            @assert Set(V) == Set(Vouter)
            @assert Set(All) == Set(Allouter)
        end

        @assert (T in F) == false "F vertex is particle-hole irreducible, so that T channel is not allowed in F"
        @assert (S in V) == false "V vertex is particle-particle irreducible, so that S channel is not allowed in V"
        @assert (T in Fouter) == false "F vertex is particle-hole irreducible, so that T channel is not allowed in F"
        @assert (S in Vouter) == false "V vertex is particle-particle irreducible, so that S channel is not allowed in V"

        g = @SVector [Vector{Green}([]) for i = 1:16]
        ver4 = new{W}(para, chan, F, V, All, Fouter, Vouter, Allouter, id[1], level, loopNum, loopidxOffset, tidxOffset, g, [], [], [[],], [])
        id[1] += 1
        @assert loopNum >= 0

        if loopNum == 0
            tidx = tidxOffset
            # bare interaction may have one, two or four independent tau variables
            if para.interactionTauNum == 1  # instantaneous interaction
                addTidx(ver4, (tidx, tidx, tidx, tidx)) #direct instant intearction
                addTidx(ver4, (tidx, tidx, tidx, tidx)) #exchange instant interaction
            end
            if para.interactionTauNum == 2  # interaction with incoming and outing τ varibales
                addTidx(ver4, (tidx, tidx, tidx, tidx))  # direct and exchange instant interaction
                addTidx(ver4, (tidx, tidx, tidx + 1, tidx + 1))  # direct dynamic interaction
                addTidx(ver4, (tidx, tidx + 1, tidx + 1, tidx))  # exchange dynamic interaction
            end
            if para.interactionTauNum == 4  # interaction with incoming and outing τ varibales
                error("Not implemented!")
                # addTidx(ver4, (tidx, tidx + 1, tidx + 2, tidx + 3))  # direct dynamic interaction
                # addTidx(ver4, (tidx, tidx + 3, tidx + 2, tidx + 1))  # exchange dynamic interaction
            end
        else # loopNum>0
            for c in chan
                for ol = 0:loopNum-1
                    bubble = Bubble(ver4, c, ol, level, id)
                    if length(bubble.map) > 0  # if zero, bubble diagram doesn't exist
                        push!(ver4.bubble, bubble)
                    end
                end
            end
            # TODO: add envolpe diagrams
            # for c in II
            # end
            test(ver4) # more test
        end
        return ver4
    end
end

function compare(A, B)
    # check if the elements of XY are the same as Z
    XY, Z = copy(A), copy(B)
    for e in XY
        if (e in Z) == false
            return false
        end
        Z = (idx = findfirst(x -> x == e, Z)) > 0 ? deleteat!(Z, idx) : Z
    end
    return length(Z) == 0
end

function test(ver4)
    if length(ver4.bubble) == 0
        return
    end

    G = ver4.G
    for bub in ver4.bubble
        Lver, Rver = bub.Lver, bub.Rver
        @assert Rver.loopidxOffset + Rver.loopNum == ver4.loopidxOffset + ver4.loopNum "Rver loopidx offset = $(Rver.loopidxOffset) + loopNum = $(Rver.loopNum) is not equal to ver4 loopidx offset = $(ver4.loopidxOffset) + loopNum = $(ver4.loopNum)"
        @assert Lver.loopNum + Rver.loopNum + 1 == ver4.loopNum
        for map in bub.map
            LverT, RverT = collect(Lver.Tpair[map.lidx]), collect(Rver.Tpair[map.ridx]) # 8 τ variables relevant for this bubble
            G1T, GxT = collect(map.G0.Tpair), collect(map.Gx.Tpair) # 4 internal variables
            ExtT = collect(ver4.Tpair[map.vidx]) # 4 external variables
            @assert compare(vcat(G1T, GxT, ExtT), vcat(LverT, RverT)) "chan $(bub.chan): G1=$G1T, Gx=$GxT, external=$ExtT don't match with Lver4 $LverT and Rver4 $RverT"
        end
    end
end


text:
# TODO: Test
function prob_miss(y::AbstractFloat, beta::Vector{Float64})
  n = length(beta)
  x = sum((y .^ (0:n-1)) .* beta)
  return MCMC.sigmoid(x)
end

function assertValidBounds(y::Vector{Float64}, p::Vector{Float64})
  n = length(y)
  @assert length(y) == length(p) == n

  @assert 2 <= n <= 3

  if n == 3
    @assert p[1] < p[2] > p[3] # Quadratic missing mechanism
    @assert y[1] < y[2] < y[3]
  elseif n ==2
    @assert p[1] > p[2] # Linear missing mechanism
    @assert y[1] < y[2]
  else
    println("In assertValidBounds: This should never be printed!")
    @assert false
  end
end

function solveBeta(y::Vector{Float64}, p::Vector{Float64})
  assertValidBounds(y, p)

  n = length(y)
  Y = [fill(1.0, n) (y) (y .^ 2)]
  beta = Y \ MCMC.logit.(p)

  @assert all(abs.(MCMC.sigmoid.(Y * beta) - p) .< 1E-10)

  return beta
end

ext: jl

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch


text:
export RadialFlow

import Base: eltype, size, length

@doc raw"""
The RadialFlow function is defined as

```math
f({\bf{x}}) = {\bf{x}} + \frac{\beta({bf{z}} - {\bf{z}}_0)}{\alpha + |{\bf{z}} - {\bf{z}}_0|}
```

with input and output dimension ``D``. Here ``{\bf{x}}\in \mathbb{R}^D`` represents the input of the function. Furthermore ``{\bf{z}}_0\in \mathbb{R}^D``, ``\alpha\in \mathbb{R}`` and ``\beta\in\mathbb{R}`` represent the parameters of the function. The function contracts and expands the input space.

This function has been introduced in:

Rezende, Danilo, and Shakir Mohamed. "Variational inference with normalizing flows." _International conference on machine learning._ PMLR, 2015.
"""
mutable struct RadialFlow{T1, T2 <: Real} <: AbstractCouplingFlow
    z0      :: T1
    α       :: T2
    β       :: T2
    function RadialFlow(z0::T1, α::T2, β::T2) where { T1, T2 <: Real}
        @assert α > 0 "The parameter α should be larger than 0."
        return new{T1,T2}(z0, float(α), float(β))
    end
end
struct RadialFlowEmpty{N} <: AbstractCouplingFlowEmpty end
RadialFlowEmpty(dim::Int) = RadialFlowEmpty{dim}()
struct RadialFlowPlaceholder <: AbstractCouplingFlowPlaceholder end
@doc raw"""
The `RadialFlow()` function creates a planar flow object. Its dimensionality is automatically set when wrapping this object inside a model. Its parameters are initialized during compilation.
"""
RadialFlow() = RadialFlowPlaceholder()
function prepare(dim::Int, flow::RadialFlowPlaceholder)
    return RadialFlowEmpty(dim)
end

# compile placeholder
compile(f::RadialFlowEmpty{1})          = RadialFlow(randn(), rand(), randn())
compile(f::RadialFlowEmpty)             = RadialFlow(randn(getdim(f)), rand(), randn())
compile(f::RadialFlowEmpty{1}, params)  = RadialFlow(params[1], params[2], params[3])
compile(f::RadialFlowEmpty, params)     = RadialFlow(params[1:getdim(f)], params[getdim(f)+1], params[getdim(f)+2])

@doc raw"""
The `RadialFlow(dim::Int64)` function creates a mutable `RadialFlow` structure with parameters corresponding to input of dimensions `dim`. The parameters are each random sampled from a standard (multivariate) normal distribution.
"""
function RadialFlow(dim::Int64)
    return RadialFlow(randn(dim), rand(), randn())
end

# number of parameters
nr_params(flow::RadialFlow)      = 2 + length(flow.z0)
nr_params(flow::RadialFlowEmpty) = 2 + getdim(flow)

# get-functions for the RadialFlow structure.
getz0(f::RadialFlow)                      = return f.z0
getα(f::RadialFlow)                       = return f.α
getβ(f::RadialFlow)                       = return f.β
getall(f::RadialFlow)                     = return f.z0, f.α, f.β
getdim(f::RadialFlowEmpty{N}) where { N } = return N

# set-functions for the RadialFlow structure
function setz0!(f::RadialFlow{T1,T2}, z0::T1) where { T1, T2 <: Real}
    @assert length(f.z0) == length(z0) "The dimensionality of the current value of z0 and its new value do not match."
    f.z0 = z0
end

function setα!(f::RadialFlow{T1,T2}, α::T2) where { T1, T2 <: Real}
    f.α = α
end

function setβ!(f::RadialFlow{T1,T2}, β::T2) where { T1, T2 <: Real }
    f.β = β
end

# custom Base function for the RadialFlow structure
eltype(f::RadialFlow{T1,T2}) where { T1 <: Real, T2 <: Real}                 = promote_type(T1, T2)
eltype(f::RadialFlow{T1,T2}) where { T1 <: AbstractVector, T2 <: Real}       = promote_type(eltype(T1), T2)
eltype(::Type{RadialFlow{T1,T2}}) where { T1 <: Real, T2 <: Real}            = promote_type(T1, T2)
eltype(::Type{RadialFlow{T1,T2}}) where { T1 <: AbstractVector, T2 <: Real}  = promote_type(eltype(T1), T2)

size(f::RadialFlow{T1,T2}) where { T1 <: Real, T2 <: Real}                   = 1
size(f::RadialFlow{T1,T2}) where { T1 <: AbstractVector, T2 <: Real}         = length(f.z0)
size(f::RadialFlowEmpty{N}) where { N }                                      = return N

length(f::RadialFlow{T1,T2}) where { T1 <: Real, T2 <: Real}                 = 1
length(f::RadialFlow{T1,T2}) where { T1 <: AbstractVector, T2 <: Real}       = length(f.z0)
length(f::RadialFlowEmpty{N}) where { N }                                    = return N

# forward pass through the RadialFlow function (multivariate input)
function _forward(f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real }

    # fetch values
    z0, α, β = getall(f)

    # calculate result
    denominator = α + norm(input - z0) # Not sure whether this is the correct norm
    denominator /= β

    result = copy(input)
    result .-= z0
    result ./= denominator
    result .+= input

    # return result
    return result

end
forward(f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real } = _forward(f, input)
Broadcast.broadcasted(::typeof(forward), f::RadialFlow{T1,T2}, input::AbstractVector{T1}) where { T1, T2 <: Real } = broadcast(_forward, Ref(f), input)

# forward pass through the RadialFlow function (univariate input)
function _forward(f::RadialFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 <: Real }

    # fetch values
    z0, α, β = getall(f)

    # calculate result
    denominator = α + norm(input - z0) # Not sure whether this is the correct norm
    denominator /= β

    result = copy(input)
    result -= z0
    result /= denominator
    result += input

    # return result
    return result

end
forward(f::RadialFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 <: Real } = _forward(f, input)
Broadcast.broadcasted(::typeof(forward), f::RadialFlow{T1,T2}, input::AbstractVector{ <: Real }) where { T1 <: Real, T2 <: Real } = broadcast(_forward, Ref(f), input)

function _forward(f::RadialFlow{T1,T2}, input) where { T1 <: Real, T2 <: Real}
    # function when the input is an array with 1 element
    @assert length(input) == 1 "Something is wrong with the dimensionality of the input to the RadialFlow flow."
    return forward(f, input[1])
end
forward(f::RadialFlow{T1,T2}, input) where { T1 <: Real, T2 <: Real } = _forward(f, input)
Broadcast.broadcasted(::typeof(forward), f::RadialFlow{T1,T2}, input::AbstractVector) where { T1 <: Real, T2 <: Real } = broadcast(_forward, Ref(f), input)

# inplace forward pass through the RadialFlow function (multivariate input)
function forward!(output::T1, f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real }

    # check dimensionality
    @assert length(output) == length(input) "The length of the preallocated vector does not seem to match the length of the input vector."

    # fetch values
    z0, α, β = getall(f)

    # calculate result
    denominator = α + norm(input - z0) # Not sure whether this is the correct norm
    denominator /= β

    output .= input
    output .-= z0
    output ./= denominator
    output .+= input

end

# jacobian of the RadialFlow function (multivariate input)
function _jacobian(f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real}

    # fetch values
    z0, α, β = getall(f)

    # # calculate result
    diff = input - z0
    result = diff*diff'
    hi = α + norm(diff)
    βh = β/hi
    ζ  = -β / hi / hi / norm(diff)
    result .*= ζ
    @inbounds for k = 1:length(input)
        result[k,k] += 1
        result[k,k] += βh
    end

    # return result
    return result

end
jacobian(f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real } = _jacobian(f, input)
Broadcast.broadcasted(::typeof(jacobian), f::RadialFlow{T1,T2}, input::AbstractVector{T1}) where { T1, T2 <: Real } = broadcast(_jacobian, Ref(f), input)

function _jacobian(f::RadialFlow{T1,T2}, input) where { T1 <: Real, T2 <: Real }
    # function when the input is an array with 1 element
    @assert length(input) == 1 "Something is wrong with the dimensionality of the input to the RadialFlow flow."
    return jacobian(f, input[1])
end
jacobian(f::RadialFlow{T1,T2}, input) where { T1 <: Real, T2 <: Real } = _jacobian(f, input)
Broadcast.broadcasted(::typeof(jacobian), f::RadialFlow{T1,T2}, input::AbstractVector) where { T1 <: Real, T2 <: Real } = broadcast(_jacobian, Ref(f), input)

# jacobian of the RadialFlow function (univariate input)
function _jacobian(f::RadialFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 <: Real }

    # fetch values
    z0, α, β = getall(f)

    # calculate result (optimized)
    diff = input - z0
    h = 1 / (α + norm(diff))
    result = 1 + β*h - β * h * h /norm(diff) * diff * diff

    # return result
    return result

end
jacobian(f::RadialFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 <: Real } = _jacobian(f, input)
Broadcast.broadcasted(::typeof(jacobian), f::RadialFlow{T1,T2}, input::AbstractVector{ <: Real }) where { T1 <: Real, T2 <: Real } = broadcast(_jacobian, Ref(f), input)

# inplace jacobian of the RadialFlow function (multivariate input)
function jacobian!(output::AbstractMatrix{T2}, f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real}

    # check whether the dimensionality is correct
    @assert size(output) == (length(input), length(f.z0)) "The dimensionality of the preallocated jacobian matrix seems incorrect."

    # fetch values
    z0, α, β = getall(f)

    # # calculate result
    diff = input - z0
    for ku = 1:length(diff)
        for kw = 1:length(diff)
            output[ku,kw] = diff[ku]*diff[kw]
        end
    end
    hi = α + norm(diff)
    βh = β/hi
    ζ  = -β / hi / hi / norm(diff)
    output .*= ζ
    @inbounds for k = 1:length(input)
        output[k,k] += 1
        output[k,k] += βh
    end

end

# determinant of the jacobian of the RadialFlow function (multivariate input)
det_jacobian(f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real} = det(jacobian(f, input))

# determinant of the jacobian of the RadialFlow function (univariate input)
function det_jacobian(f::RadialFlow{T, T}, input::T) where { T <: Real }

    # fetch values
    z0, α, β = getall(f)

    # calculate result
    r = norm(input - z0)
    h = 1 / (α + r)
    result = ( 1 + β*h - β*h^2*r ) * ( 1 + β*h )^(length(input)-1)

    # return result
    return result

end

# extra utility function (multivariate)
inv_jacobian(f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real }       = inv(jacobian(f, input))
absdet_jacobian(f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real }    = abs(det_jacobian(f, input))
logabsdet_jacobian(f::RadialFlow{T1,T2}, input::T1) where { T1, T2 <: Real}  = logabsdet(jacobian(f, input))

# extra utility functions (univariate)
inv_jacobian(f::RadialFlow{T,T}, input::T) where { T <: Real }               = 1.0 / jacobian(f, input)
absdet_jacobian(f::RadialFlow{T,T}, input::T) where { T <: Real }            = abs(det_jacobian(f, input))
logabsdet_jacobian(f::RadialFlow{T,T}, input::T) where { T <: Real}          = log(absdet_jacobian(f, input))

function inv_jacobian(f::RadialFlow{T1,T2}, input::T3) where { T1 <: Real, T2 <: Real, T3 }
    # function when the input is an array with 1 element
    @assert length(input) == 1 "Something is wrong with the dimensionality of the input to the RadialFlow flow."
    return inv_jacobian(f, input[1])
end


text:
using DSP

"""
    sos = butter2sos(order, Fc, Fs, type=:lowpass)

Z-domain IIR Butterworth design of Lowpass/Highpass/Allpass filters.
`Fc` indicates the corner "-3 dB" natural frequency point, whereas `Fs` indicates
the discrete sample rate and `order` specifies the cascade filter order. The returned
`sos` SecondOrderSection type is compatible with DSP.jl.
"""
function butter2sos(order::Integer, Fc::Real, Fs::Real, type::Symbol=:lowpass)
    (Fc > 0 && Fs > 0) || throw(ArgumentError("Fc and/or Fs must be greater than 0 Hz (DC)."))
    (Fc <= Fs/2) || throw(ArgumentError("Upperband `Fhi` cannot exceed Nyquist."))
    (type ∈ [:lowpass, :highpass, :allpass]) || throw(ArgumentError("type must indicate ':lowpass', ':highpass', or ':allpass.'"))
    Nstages = ceil(Int, order/2)
    matptr = Vector{Float64}(undef, Nstages*6)
    if (type == :lowpass)
        ftype = Int(0)
    elseif (type == :highpass)
        ftype = Int(1)
    else
        ftype = Int(2)
    end
    ccall((:jl_butter, "butterlib"), Cvoid, (Cint, Cdouble, Cdouble, Cint, Ptr{Float64}), order, Fc, Fs, ftype, matptr)
    sosmat = Matrix{Float64}(transpose(reshape(matptr, (6, Nstages))))
    Vb = Vector{Biquad{:z, Float64}}(undef, Nstages)
    K = prod(sosmat[:,1]) # normalize by lead coefficient and extract gain, (whichever ordering.)
    kstg = 1
    for iStg = 1:Nstages
        kstg = sosmat[iStg,1]
        sosmat[iStg,1:3] ./= kstg
        Vb[iStg] = Biquad{:z, Float64}(sosmat[iStg,1], sosmat[iStg,2], sosmat[iStg,3], sosmat[iStg,5], sosmat[iStg,6])
    end
    !any(isnan, sosmat) || error("Filter specs ($(order), $(Fc), $(Fs), $(type)) produced a gain of zero.")
    SecondOrderSections(Vb, K)
end

"""
    sos = butterband(order, Fl, Fh, Fs, type=:bandpass)

Z-domain IIR Butterworth design of Bandpass/Bandstop filters.
`Fl` and `Fh` indicate the lower and upper natural frequency edges,
with `Fs` indicating the discrete sample rate, and `order` specifies
the cascaded Biquads filter order.
"""
function butterband(order::Integer, Fl::Real, Fh::Real, Fs::Real, type::Symbol=:bandpass)
    (Fs > 0) || throw(ArgumentError("Fs sampling rate must be greater than 0 Hz (DC)."))
    (order > 0) || throw(ArgumentError("Filter order must be greater than zero."))
    (type ∈ [:bandpass, :bandstop]) || throw(ArgumentError("type must indicate ':bandpass' or ':bandstop'."))
    if (Fl > Fh)
        Fl, Fh = Fh, Fl
    end
    ((Fl > 0 && Fl < Fs/2) && (Fh > 0 && Fh < Fs/2)) || throw(ArgumentError("Upper/lower bounds must be in range (0, Fs/2]."))
    matptr = Vector{Float64}(undef, order*6)
    ftype = (type == :bandpass) ? Int(0) : Int(1)
    ccall((:jl_butterband, "butterlib"), Cvoid, (Cint, Cdouble, Cdouble, Cdouble, Cint, Ptr{Float64}), order, Fl, Fh, Fs, ftype, matptr)
    sosmat = Matrix{Float64}(transpose(reshape(matptr, (6, order))))
    Vb = Vector{Biquad{:z, Float64}}(undef, order)
    kstg = 1
    K = prod(sosmat[:,1])
    for iStg = 1:order
        kstg = sosmat[iStg,1]
        sosmat[iStg,1:3] ./= kstg
        Vb[iStg] = Biquad{:z, Float64}(sosmat[iStg,1], sosmat[iStg,2], sosmat[iStg,3], sosmat[iStg,5], sosmat[iStg,6])
    end
    !any(isnan, sosmat) || error("Filter specs ($(order), $(Fl), $(Fh), $(Fs), $(type)) produced a gain of zero.")
    SecondOrderSections(Vb,K)
end


text:
using Millboard
using Test

Millboard.set_table_mode(:grid)
Base.eval(:(have_color = false))

board = [(+) (*); (-) (/)]
@test """
+---+---+---+
|   | 1 | 2 |
+===+===+===+
| 1 | + | * |
+---+---+---+
| 2 | - | / |
+---+---+---+""" == table(board) |> string

Base.eval(:(have_color = true))


text:
"""
test_args(args...; keyword = "keyword")

# Description
- Tests splat and keyword arguments

# Argument
- args... : Whatever your heart desires

# Keyword Arguments
- keyword: default value = "keyword"

# Return
- Nothing
"""
function test_args(args...; keyword = "keyword")
    println("The function argument is $args")
    println("The members are")
    for i in args
        println(i)
    end
    println("The keyword argument is ")
    println(keyword)
    return nothing
end

max_stars_repo_name: sandreza/Threadicare


text:
include("./Utils.jl")
using .Utils
using Plots

mycoolfunction()


text:
"""
    _can_propagate_constant(::NodeOp) -> Bool

Return true for ops which can propagate constant values if all their parents are constant.
"""
_can_propagate_constant(::NodeOp) = false

"""
    _propagate_constant_value(op::NodeOp{T}, parents::NTuple{N, Node}) -> T

Given that all parents are constants, get the value of the constant node we should output.
This assumes that `_can_propagate_constant(op)` is true.
"""
function _propagate_constant_value end

"""
    obtain_node(parents::NTuple{N,Node}, op::NodeOp) -> Node

Get a node for the given `op` and `parents`. If an equivalent node already exists in the
global identity map, use that one, otherwise create a new node, add to the identity map, and
return it.

# Constant propagation
If all `parents` are constant nodes, and `op` has a well-defined operation on constant
inputs, we will immediately perform the computation and return a constant node wrapping the
computed value.
"""
function obtain_node(parents::NTuple{N,Node}, op::NodeOp) where {N}
    if !isempty(parents) && _can_propagate_constant(op) && all(_is_constant, parents)
        # Constant propagation, since all inputs are constants & the op supports it.
        # Note that `constant` does, itself, call through to `obtain_node`, so the node will
        # be correctly registered with the id_map.
        return constant(_propagate_constant_value(op, parents))
    end

    return obtain_node!(global_identity_map(), parents, op)
end

"""
    ancestors(graph::AbstractGraph{T}, sources::AbstractVector{T}) -> Vector{T}

Given a (directed) graph, with edges from child to parent, find all ancestor vertices of the
given `sources`.

The result will be ordered such that the parents of any given vertex will always come before
the vertex itself.
"""
function ancestors(graph::AbstractGraph{T}, sources::AbstractVector{T}) where {T}
    isempty(sources) && throw(ArgumentError("Need at least 1 source"))

    seen = Set{T}()

    stack = Vector(sources)

    # Get the search data for the current source.
    while !isempty(stack)
        vertex = popfirst!(stack)
        # Either continue if we've seen this vertex before, or add it to the seen nodes.
        if vertex in seen
            continue
        else
            push!(seen, vertex)
        end
        append!(stack, outneighbors(graph, vertex))
    end

    # This won't be the most efficient way of ordering the results, especially when we only
    # care about a small subgraph.
    order = topological_sort_by_dfs(graph)
    reverse!(order)
    return [v for v in order if v in seen]
end

struct NodeIterator
    nodes::Vector{Node}
end

struct NodeIteratorState
    stack::Vector{Node}
    seen::Set{Node}
    NodeIteratorState(ni::NodeIterator) = new(ni.nodes, Set())
end

Base.IteratorSize(::Type{NodeIterator}) = Base.SizeUnknown()
Base.eltype(::Type{NodeIterator}) = Node
Base.iterate(ni::NodeIterator) = iterate(ni, NodeIteratorState(ni))
function Base.iterate(::NodeIterator, state::NodeIteratorState)
    # Find a node on the stack that we haven't seen already.
    local node
    while true
        isempty(state.stack) && return nothing
        node = popfirst!(state.stack)
        in(node, state.seen) || break
    end

    # We now know that this is a new node, so process it.
    append!(state.stack, parents(node))
    push!(state.seen, node)

    return (node, state)
end

iternodes(nodes::AbstractVector{<:Node}) = NodeIterator(Vector{Node}(nodes))

"""
    ancestors(nodes)

Get a list of all nodes in the graph defined by `nodes`, including all parents.
    * Every node in the graph will be visited exactly once.
    * The parents of any vertex will always come before the vertex itself.
"""
function ancestors(nodes::AbstractVector{<:Node})
    # Construct a LightGraphs representation of the whole node graph.
    node_to_vertex = Bijection(Dict(n => i for (i, n) in enumerate(iternodes(nodes))))

    # Initialising a SimpleDiGraph via an edge list is more efficient than calling add_edge!
    # repeatedly.
    #! format: off
    edges = Edge{Int64}[
        Edge(i, node_to_vertex[parent])
        for (node, i) in node_to_vertex
        for parent in parents(node)
    ]
    #! format: on
    light_graph = SimpleDiGraph(edges)

    # Suppose we have nodes with no parents on children in the graph - these will not show
    # up in the edge list, and may hence require manual addition.
    add_vertices!(light_graph, length(node_to_vertex) - nv(light_graph))

    vertices = [node_to_vertex[node] for node in nodes]
    ancestor_vertices = ancestors(light_graph, vertices)
    return [inverse(node_to_vertex, vertex) for vertex in ancestor_vertices]
end


text:
#=
For details, please visit:
https://computationalmindset.com/en/neural-networks/experiments-with-neural-odes-in-julia.html#exp1

Parametric system to train:
x' = a1 x + b1 y + c1 e^(-d1 t)
y' = a2 x + b2 y + c2 e^(-d2 t)
x(0) = 0
y(0) = 0

to generate dataset use the following param values
a1 = 1.11
b1 = 2.43
c1 = -3.66
d1 = 1.37
a2 = 2.89
b2 = -1.97
c2 = 4.58
d2 = 2.86

Analytical solution is:
https://www.wolframalpha.com/input/?i=x%27+%3D+1.11+x+%2B+2.43+y+%2B+-3.66+e%5E%28-1.37+t%29%3B+y%27+%3D+2.89+x+-1.97+y+%2B+4.58+e%5E+%28-2.86+t%29%3B+x%280%29%3D0%3B+y%280%29%3D0%3B
=#

using Flux, DiffEqFlux, DifferentialEquations, Plots

function parametric_ode_system!(du,u,p,t)
  x, y = u
  a1, b1, c1, d1, a2, b2, c2, d2 = p
  du[1] = dx = a1*x + b1*y + c1*exp(-d1*t)
  du[2] = dy = a2*x + b2*y + c2*exp(-d2*t)
end

true_params = [1.11, 2.43, -3.66, 1.37, 2.89, -1.97, 4.58, 2.86]

an_sol_x(t) =
  -1.38778e-17 * exp(-8.99002 * t) -
  2.77556e-17 * exp(-7.50002 * t) +
  3.28757 * exp(-3.49501 * t) -
  3.18949 * exp(-2.86 * t) +
  0.258028 * exp(-1.37 * t) -
  0.356108 * exp(2.63501 * t) +
  4.44089e-16 * exp(3.27002 * t) +
  1.11022e-16 * exp(4.76002 * t)
an_sol_y(t) =
  -6.23016 * exp(-3.49501 * t) +
  5.21081 * exp(-2.86 * t) +
  1.24284 * exp(-1.37 * t) -
  0.223485 * exp(2.63501 * t) +
  2.77556e-17 * exp(4.76002 * t)

tbegin=0.0
tend=1.5
tstep=0.01
trange = tbegin:tstep:tend
u0 = [0.0,0.0]
tspan = (tbegin,tend)
p = ones(8)

prob = ODEProblem(parametric_ode_system!, u0, tspan, p)

function net()
    solve(prob, Tsit5(), p=p, saveat=trange)
end

dataset_outs = [an_sol_x.(trange), an_sol_y.(trange)]
function loss_func()
  pred = net()
  sum(abs2, dataset_outs[1] .- pred[1,:]) +
  sum(abs2, dataset_outs[2] .- pred[2,:])
end

epochs = 1000
learning_rate = 0.05
data = Iterators.repeated((), epochs)
opt = ADAM(learning_rate)

callback_func = function ()
  loss_value = loss_func()
  println("Loss: ", loss_value)
end
fparams = Flux.params(p)
Flux.train!(loss_func, fparams, data, opt, cb=callback_func)

predict_prob = ODEProblem(parametric_ode_system!, u0, tspan, p)
predict_sol = solve(prob, Tsit5(), saveat=trange)
x_predict_sol = [u[1] for u in predict_sol.u]
y_predict_sol = [u[2] for u in predict_sol.u]

println("Learned parameters:", p)

plot(trange, dataset_outs[1],
    linewidth=2, ls=:dash,
    title="Neural ODEs to fit params",
    xaxis="t",
    label="dataset x(t)",
    legend=true)
plot!(trange, dataset_outs[2],
    linewidth=2, ls=:dash,
    label="dataset y(t)")
plot!(predict_sol.t, x_predict_sol,
    linewidth=1,
    label="predicted x(t)")
plot!(predict_sol.t, y_predict_sol,
    linewidth=1,
    label="predicted y(t)")

max_stars_repo_name: ettoremessina/differential-equations

max_issues_repo_name: ettoremessina/differential-equations

max_forks_repo_name: ettoremessina/differential-equations


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    DiscretizationMethod

A method for discretizing geometries into meshes.
"""
abstract type DiscretizationMethod end

"""
    discretize(geometry, method)

Discretize `geometry` with discretization `method`.
"""
function discretize end

discretize(box::Box, method::DiscretizationMethod) =
  discretize(boundary(box), method)

function discretize(polygon::Polygon{Dim,T}, method::DiscretizationMethod) where {Dim,T}
  # build bridges in case the polygon has holes,
  # i.e. reduce to a single outer boundary
  chain, dups = bridge(unique(polygon), width=2atol(T))

  # discretize using outer boundary
  mesh = discretize(chain, method)

  if isempty(dups)
    # nothing to be done, return mesh
    mesh
  else
    # remove duplicate vertices
    points = collect(vertices(mesh))
    for (i, j) in dups
      points[i] = centroid(Segment(points[i], points[j]))
    end
    repeated = sort(last.(dups))
    deleteat!(points, repeated)

    # adjust connectivities
    elems = elements(topology(mesh))
    twin  = Dict(reverse.(dups))
    rrep  = reverse(repeated)
    einds = map(elems) do elem
      inds = indices(elem)
      [get(twin, ind, ind) for ind in inds]
    end
    for inds in einds
      for r in rrep
        for i in 1:length(inds)
          inds[i] > r && (inds[i] -= 1)
        end
      end
    end
    connec = connect.(Tuple.(einds))

    # return mesh without duplicates
    SimpleMesh(points, connec)
  end
end

discretize(multi::Multi, method::DiscretizationMethod) =
  mapreduce(geom -> discretize(geom, method), merge, multi)

"""
    triangulate(geometry)

Discretize `geometry` of parametric dimension 2 into
triangles using an appropriate discretization method.
"""
triangulate(box::Box{2}) = discretize(box, Dehn1899())
triangulate(ngon::Ngon) = discretize(ngon, Dehn1899())
triangulate(poly::PolyArea) = discretize(poly, FIST())
triangulate(multi::Multi) = mapreduce(triangulate, merge, multi)

# ----------------
# IMPLEMENTATIONS
# ----------------

include("discretization/fist.jl")
include("discretization/dehn.jl")


text:
using Test
using NQCDynamics
using NQCDynamics.InitialConditions
using StatsBase: mean
using ComponentArrays
using Random
using FiniteDiff
Random.seed!(1)

ats = [3, 1, 40]
Ds = [1, 3, 3]
Ts = [0.5, 1, 10]
beads = [1, 16, 2]

@testset "Classical" begin

    @testset "get_proposal" begin
        for (natoms, DoFs, T) in zip(ats, Ds, Ts)
            sim = Simulation(Atoms(rand(natoms)), Harmonic(dofs=DoFs); temperature=T)
            proposals = ThermalMonteCarlo.get_proposal(sim, Dict(:X=>0.5), 1, 0)
            @test length(proposals.proposal) == natoms * DoFs
        end
    end

    @testset "Energy expectation" begin
        for (natoms, DoFs, T) in zip(ats, Ds, Ts)
            sim = Simulation(Atoms(rand(natoms)), Harmonic(dofs=DoFs); temperature=T)
            R0 = randn(size(sim))
            chain = ThermalMonteCarlo.run_advancedmh_sampling(sim, R0, 1e4, Dict(:X=>1); move_ratio=0.5)
            energy = Estimators.@estimate potential_energy(sim, chain)
            @test energy / (DoFs*natoms) ≈ T/2 rtol=1e-1
        end
    end

    @testset "HMC gradient" begin
        for (natoms, DoFs, T) in zip(ats, Ds, Ts)
            sim = Simulation(Atoms(rand(natoms)), Harmonic(dofs=DoFs); temperature=T)
            density = ThermalMonteCarlo.get_density_function(sim)
            r = randn(prod(size(sim)))
            grad = FiniteDiff.finite_difference_gradient(density, r)
            grad_func = ThermalMonteCarlo.get_deriv_function(sim)
            @test grad ≈ grad_func(r)[2]
        end
    end

    @testset "Energy expectation HMC" begin
        for (natoms, DoFs, T) in zip(ats, Ds, Ts)
            sim = Simulation(Atoms(rand(natoms)), Harmonic(dofs=DoFs); temperature=T)
            R0 = randn(size(sim))
            chain, stats = ThermalMonteCarlo.run_advancedhmc_sampling(sim, R0, 1e4; verbose=false)
            energy = Estimators.@estimate potential_energy(sim, chain)
            @test energy / (DoFs*natoms) ≈ T/2 rtol=1e-1
        end
    end

end

@testset "Ring polymer" begin

    @testset "get_proposal" begin
        for (natoms, DoFs, T, nbeads) in zip(ats, Ds, Ts, beads)
            sim = RingPolymerSimulation(Atoms(vcat([:C], fill(:H, natoms-1))), Harmonic(dofs=DoFs), nbeads; temperature=T, quantum_nuclei=[:H])
            proposals = ThermalMonteCarlo.get_proposal(sim, Dict(:H=>0.5, :C=>0.1), 1, 0)
            @test length(proposals.proposal) == natoms * DoFs * nbeads
        end
    end

    @testset "Energy expectation" begin
        for (natoms, DoFs, T, nbeads) in zip(ats, Ds, Ts, beads)
            sim = RingPolymerSimulation(Atoms(rand(natoms)), Harmonic(dofs=DoFs), nbeads; temperature=T)
            R0 = zeros(DoFs, natoms, nbeads)
            chain = ThermalMonteCarlo.run_advancedmh_sampling(sim, R0, 1e4, Dict(:X=>1); move_ratio=0.3, internal_ratio=0.95)
            energy = DynamicsUtils.classical_potential_energy.(sim, chain)
            @test mean(energy) / (DoFs*natoms*nbeads) ≈ nbeads*T/2 rtol=1e-1
        end
    end

    @testset "Quantum harmonic oscillator" begin
        ω = 100
        T = 1
        sim = RingPolymerSimulation(Atoms(1), Harmonic(ω=ω), 200; temperature=T)
        r = zeros(size(sim))
        chain = ThermalMonteCarlo.run_advancedmh_sampling(sim, r, 1e5, Dict(:X=>1/100); move_ratio=0.0, internal_ratio=0.95)

        exact = ω/2 * coth(ω / 2T)
        estimate = Estimators.@estimate total_energy(sim, chain)
        @test exact ≈ estimate rtol=1e-1
    end

    @testset "Free ring polymer radius_of_gyration" begin
        T = 3.0
        analytic_Rg = sqrt(1/12T)
        sim = RingPolymerSimulation(Atoms(1), Free(), 10; temperature=T)
        r = zeros(size(sim))
        chain = ThermalMonteCarlo.run_advancedmh_sampling(sim, r, 4e4, Dict(:X=>1); move_ratio=0.0, internal_ratio=0.5)
        estimate = Estimators.@estimate radius_of_gyration(sim, chain)
        @test analytic_Rg ≈ estimate[1] atol=1e-1
    end

end


text:
"""
    SplitQuaternion{T <: Real} = Hyperbolic{Binion{T}}

A split-quaternion is a hyperbolic Cayley-Dickson construct with complex numbers.
"""
const SplitQuaternion{T <: Real} = Hyperbolic{Binion{T}}

function SplitQuaternion(z::Binion{T}) where T <: Real
    SplitQuaternion{T}(z, zero(Binion{T}))
end

function SplitQuaternion(a::T, b::T, c::T, d::T) where T <: Real
    SplitQuaternion{T}(Binion{T}(a, b), Binion{T}(c, d))
end

function SplitQuaternion(a::T, b::T, c::T) where T <: Real
    SplitQuaternion{T}(Binion{T}(a, b), Binion{T}(c, zero(T)))
end

function SplitQuaternion(a::T, b::T) where T <: Real
    SplitQuaternion{T}(Binion{T}(a, b), zero(Binion{T}))
end

function SplitQuaternion(a::T) where T <: Real
    SplitQuaternion{T}(Binion{T}(a, zero(T)), zero(Binion{T}))
end

function SplitQuaternion(a::Real, b::Real)
    SplitQuaternion(promote(a, b)...)
end

function SplitQuaternion(a::Real, b::Real, c::Real)
    SplitQuaternion(promote(a, b, c)...)
end

function SplitQuaternion(a::Real, b::Real, c::Real, d::Real)
    SplitQuaternion(promote(a, b, c, d)...)
end

"""
    unreal(z::SplitQuaternion)

The unreal part of a split-quaternion is a 3-vector.
"""
function unreal(z::SplitQuaternion)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::SplitQuaternion)
    print(io, "[1: ")
    print(io, z.l.l)
    print(io, ", i: ")
    print(io, z.l.r)
    print(io, ", s: ")
    print(io, z.r.l)
    print(io, ", is: ")
    print(io, z.r.r)
    print(io, "]")
end

function random(::Type{SplitQuaternion{T}}) where T <: Real
    SplitQuaternion{T}(random(Binion{T}), random(Binion{T}))
end


text:
function cohort(data, dose, time, covariates...; dose_args = covariates, id = :ID)
    population = Subject[]
    for row in eachrow(data)
        values = map(covariate -> row[covariate], covariates)
        subject = Subject(;
            id = row[id],
            events = dose(map(covariate -> row[covariate], dose_args)...),
            time = time,
            covariates = NamedTuple{covariates}(values)
        )
        push!(population, subject)
    end
    return population
end

function nonmemify(population; output = missing, missingvalue = "0")
    df = DataFrame(population)
    rename!(df, Symbol.(uppercase.(String.(names(df)))))
    df = something.(df, missing)
    if missingvalue !== missing
        df = coalesce.(df, missingvalue)
    end
    if output !== missing
        CSV.write(output, df; missingstring = missingvalue)
    end
    return df
end


text:
using KernelAbstractions

function averageK(K, i, j, k)
    Kw = harmonicMean(K[i - 1, j,     k    ],     K[i,     j,     k    ])
    Ke = harmonicMean(K[i,     j,     k    ],     K[i + 1, j,     k    ])
    Ks = harmonicMean(K[i,     j - 1, k    ],     K[i,     j,     k    ])
    Kn = harmonicMean(K[i,     j,     k    ],     K[i,     j + 1, k    ])
    Kb = harmonicMean(K[i,     j,     k - 1],     K[i,     j,     k    ])
    Kt = harmonicMean(K[i,     j,     k    ],     K[i,     j,     k + 1])
    return Kw, Ke, Ks, Kn, Kb, Kt
end

function harmonicMean(v1, v2)
    return 2.0 * v1 * v2 / (v1 + v2)
end

@kernel function pressureKernel!(source, h, hⁿ⁺¹, Δx, Δy, Δz, Δt, K, SS)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        F = K[i,j,k] * (
            (h[i + 1,j,k] + h[i - 1,j,k] - 2.0 * h[i,j,k]) / (Δx * Δx) +
            (h[i,j + 1,k] + h[i,j - 1,k] - 2.0 * h[i,j,k]) / (Δy * Δy) +
            (h[i,j,k + 1] + h[i,j,k - 1] - 2.0 * h[i,j,k]) / (Δz * Δz)
            ) + source[i,j,k] / (Δx * Δy)

        hⁿ⁺¹[i,j,k] = h[i,j,k] + Δt * F * (1.0 / SS)
    end

end

@inline function pressureKernelAveraged!(i, j, k, source, h, u, v, w, hⁿ⁺¹, uⁿ⁺¹, vⁿ⁺¹, wⁿ⁺¹, Δx, Δy, Δz, Δt, K, SS)

    Kw, Ke, Ks, Kn, Kb, Kt = averageK(K, i, j, k)

    F = (1.0 / Δx * Δx) * (Ke * (h[i + 1,j  ,k  ] - h[i,j,k]) - Kw * (h[i,j,k] - h[i - 1,j  ,k  ])) +
        (1.0 / Δy * Δy) * (Kn * (h[i  ,j + 1,k  ] - h[i,j,k]) - Ks * (h[i,j,k] - h[i  ,j - 1,k  ])) +
        (1.0 / Δz * Δz) * (Kt * (h[i  ,j  ,k + 1] - h[i,j,k]) - Kb * (h[i,j,k] - h[i  ,j  ,k - 1])) +
        source[i,j,k] / (Δx * Δy)

    hⁿ⁺¹[i,j,k] = h[i,j,k] + Δt * F * (1.0 / SS)
end


text:
@testset "Algebraic Decision Diagrams" begin
    @info "ADD"
    ADD = AlgebraicDecisionDiagrams
    @show x1 = indicator(Int,1) # equivalent to x1 = Node(1,Base.zero(Int),Base.one(Int))
    @show one(x1) - x1
    @show ¬x1
    @show apply(x -> 1-x, x1)
    @show x1 + (¬x1)
    @show x1 * (¬x1)
    @show x2 = indicator(Int,2)
    @show x1 + x2
    @show x1 * x2
    @show apply(x -> x^2, x1 + x2)
    @show ⊕(x,y) = (x + y) % 2 # sum mod 2
    @show apply(⊕,x1,x2)
    @show x3 = indicator(Int,3)

    @info "Odd Parity Function"
    @show apply(⊕,x1,apply(⊕,x2,x3))
    @show apply(max, x1 * x2, x3)

    @info "Example from Fig. 1 in Bryant's paper"
    @show apply(max, (¬x1)*x2, x1*x3)

    @info "Example in Fig. 5 of Bryant's paper"
    A3 = Node(3,zero(Int),one(Int))
    @show A4 = ADD.reduce(Node(1,Node(2,zero(Int),A3),Node(2,A3,A3)))
    # for (i,α) in enumerate(x1*x2)
    #     println(i, "\n", α)
    # end
    println()

    @info "Directed Graph ADD"

    @info "Example in Fig 1. in Bahar et al. 1997's paper"
    x0 = indicator(Int,0)
    nx0 = indicator(Int,0,true)
    x1 = indicator(Int,2)
    nx1 = indicator(Int,2,true)
    y0 = indicator(Int,1)
    ny0 = indicator(Int,1,true)
    y1 = indicator(Int,3)
    ny1 = indicator(Int,3,true)
    @show graph = Terminal(2)*nx0*nx1*ny0*y1 + Terminal(2)*nx0*nx1*y0*ny1 + Terminal(2)*nx0*x1*ny0*y1 + Terminal(2)*nx0*x1*y0*ny1 + Terminal(4)*x0*nx1*y0*y1 + Terminal(4)*x0*x1*ny0*y1
    @info "Get scope"
    # collect nonterminal nodes
    nt = filter(n -> isa(n,Node), collect(graph))
    @show scope = mapreduce(index,union,nt) # 0,1,2,3
    @test length(intersect(scope,[0,1,2,3])) == 4
    println()

    @info "Matrix examples in page 9 of Bryant's paper"
    @show f = (¬x0)*(¬y0) + x0*y0
    @show g = Terminal(4)*nx0 + Terminal(2)*x0
    @show h = f + g
    println()
    @info "Restriction"
    @show restrict(h,0,false)
    @show h | (1 => true)
    @show (h | (index(x0) => true)) + (h | (index(x0) => false))
    @test value((h | (index(x0) => true)) + (h | (index(x0) => false))) == 7
    @show (h | (index(y0) => true)) + (h | (index(y0) => false))
    println()
    @info "marginalization"
    @show marginalize(h, index(x0))
    @show marginalize(h, index(y0))
    @show A4
    @show (A4 | (2 => true)) + (A4 | (2 => false))
    @show A5 = marginalize(A4, 2)
end # testset


text:
# zwart_powell - compare Matlab and Julia versions

using MATLAB
using Sinograms: zwart_powell
using LazyGrids: ndgrid_array
#using Plots
#using MIRTjim: jim

# setup matlab to use MIRT
if !@isdefined(irtdir)
	ENV["MATLAB_ROOT"] = "/Applications/freeware/matlab"

	irtdir = "/Users/fessler/src/matlab/alg"
	tmp = "addpath('$irtdir')"
	eval_string(tmp)
	mat"setup"
end

r = LinRange(-1, 1, 2001) * 2
phi = LinRange(0, 2π, 361)
jsino = zwart_powell.(r, phi')

rr, pp = ndgrid_array(r, phi) # for matlab...
@mput rr
@mput pp
mat"msino = ir_radon_zwart_powell(pp, rr)"
@mget msino

@assert jsino ≈ msino

#jim(jim(jsino, "julia"), jim(msino, "matlab"), jim(jsino-msino, "diff"))


text:
module SpringCollab2020DashSpring

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/dashSpring" DashSpring(x::Integer, y::Integer, playerCanUse::Bool=true)
@mapdef Entity "SpringCollab2020/wallDashSpringRight" DashSpringRight(x::Integer, y::Integer)
@mapdef Entity "SpringCollab2020/wallDashSpringLeft" DashSpringLeft(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Dash Spring (Up, Spring Collab 2020)" => Ahorn.EntityPlacement(
        DashSpring
    ),
    "Dash Spring (Left, Spring Collab 2020)" => Ahorn.EntityPlacement(
        DashSpringRight
    ),
    "Dash Spring (Right, Spring Collab 2020)" => Ahorn.EntityPlacement(
        DashSpringLeft
    ),
)

function Ahorn.selection(entity::DashSpring)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 6, y - 3, 12, 5)
end

function Ahorn.selection(entity::DashSpringLeft)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 1, y - 6, 5, 12)
end

function Ahorn.selection(entity::DashSpringRight)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 4, y - 6, 5, 12)
end

sprite = "objects/SpringCollab2020/dashSpring/00.png"

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DashSpring, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, -8)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DashSpringLeft, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 24, 0, rot=pi / 2)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DashSpringRight, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, -8, 16, rot=-pi / 2)

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# https://github.com/JuliaLang/julia/issues/4707

#----------------- Get the A matrix
function getDivGrad(n1,n2,n3)
        D1 = kron(speye(n3),kron(speye(n2),ddx(n1)))
        D2 = kron(speye(n3),kron(ddx(n2),speye(n1)))
        D3 = kron(ddx(n3),kron(speye(n2),speye(n1)))

        # DIV from faces to cell-centers
        Div = [D1 D2 D3]
        return Div*Div';
end

#----------------- 1D finite difference on staggered grid
function ddx(n)
    return spdiags(ones(n)*[-1 1],[0,1],n,n+1)
end

#------------- Build a diagonal matrix
function spdiags(B,d,m,n)
    d = d[:]
    p = length(d)

    len = zeros(Int, p+1, 1)
    for k = 1:p
        len[k+1] = len[k] + length(max(1,1-d[k]):min(m,n-d[k]))
    end

    a = zeros(Int, len[p+1], 3)
    for k = 1:p
        # Append new d[k]-th diagonal to compact form
        i = max(1,1-d[k]):min(m,n-d[k])
        a[(len[k]+1):len[k+1],:] = [i i+d[k] B[i+(m>=n)*d[k],k]]
    end

    return sparse(a[:,1],a[:,2],a[:,3],m,n)
end


text:
#
# Correctness Tests
#

using Base.Test
using IntertemporalChoiceHeuristics

my_tests = (
    "show.jl",
    "fits.jl",
    "gradients.jl",
    "optimization_algorithms.jl",
    "random_initial_parameters.jl",
    "consistency.jl",
)

println("Running tests:")

@testset "All tests" begin
    for my_test in my_tests
        @printf(" * %s\n", my_test)
        include(my_test)
    end
end


text:
using DelimitedFiles
using Statistics

open("data_3.dat") do f
   # line_number
  line = 1
  data = zeros(1000, 12)
  # read till end of file
  for j in 1:1000
     # read a new / next line for every itera
     s = readline(f)
     for ii in 1:12
         data[line, ii] = parse(Int, s[ii])
     end
     line += 1
  end

  jj = 0
  for ii in 1:12
      if ii == 1
         global dat_select_oxy = data
      else
         global dat_select_oxy = dat_select_oxy[ind_oxy, :]
      end
      if mean(dat_select_oxy[:,ii]) < 0.5
          oxygen_rating = 0.0
      else
          oxygen_rating = 1.0
      end

      global ind_oxy = findall(x -> x == oxygen_rating, dat_select_oxy[:, ii] )
      if length(ind_oxy) == 1 && jj == 0
          global dat_oxy = dat_select_oxy[ind_oxy, :]
          jj += 1
      end
  end

  jj = 0
  for ii in 1:12
      if ii == 1
          global dat_select_co2 = data
      else
          global dat_select_co2 = dat_select_co2[ind_co2, :]
      end

      # if 0 is more common than 1, then the mean is <= 0.5
      if mean(dat_select_co2[:,ii]) < 0.5
          co2_rating = 1.0
      else
          co2_rating = 0.0
      end
      global ind_co2 = findall(x -> x == co2_rating, dat_select_co2[:, ii] )
      if length(ind_co2) == 1 && jj == 0
          global dat_co2 = dat_select_co2[ind_co2, :]

          jj += 1
      end
   end
   println(dat_oxy)
   println(dat_co2)
   binary_o = 0
   for d in dat_oxy
      binary_o = binary_o * 2 + d
   end

   binary_c = 0
   for d in dat_co2
       binary_c = binary_c * 2 + d
   end
   println(binary_o, binary_c)
   println(binary_o * binary_c)

end


text:
function sbatch_spawn(tests::AbstractVector, objects::Dict;
                      batch_size=100,
                      time_per_batch="20:00",
                      data_dir=joinpath(get(ENV, "SCRATCH", tempdir()),
                                        string("sim_data_", Dates.format(Dates.now(),"u_d_HH_MM"))),
                      submit_command="sbatch",
                      template_name="sherlock.sh",
                      job_name=string("Multilane_", randstring())
                      )

    stats = setup_stats(tests, objects)
    # shuffle so that batches take approximately the same amount of time
    sort!(stats, cols=:uuid)
    objects["tests"] = Dict(t.key=>t for t in tests)

    try
        mkdir(data_dir)
    catch ex
        warn("When creating data_dir, got error $ex. Does the dir already exist?\nContinuing anyways.")
    end
    objectname = joinpath(data_dir, string("objects_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld"))
    JLD.save(objectname, objects)
    println("objects saved to $objectname")

    nb_sims = nrow(stats)
    nb_batches = cld(nb_sims, batch_size)
    @assert rem(nb_sims, batch_size) == 0
    results_file_list = []

    listname = joinpath(data_dir, string("list_\${SLURM_ARRAY_TASK_ID}_of_$(nb_batches).jld"))

    tpl = Mustache.parse(readstring(joinpath(Pkg.dir("Multilane"), "templates", template_name)))
    sbatch = Mustache.render(tpl,
                    job_name=job_name,
                    time=time_per_batch,
                    object_file_path=objectname,
                    list_file_path=listname,
                    nb_batches=nb_batches,
                    data_dir=data_dir)

    sbatchname = joinpath(data_dir, "$job_name.sbatch")
    open(sbatchname, "w") do f
        write(f, sbatch)
    end

    for i in 1:nb_batches
        println("\rpreparing job $i of $nb_batches")

        these_stats = stats[(i-1)*batch_size+1:i*batch_size, :]

        this_listname = replace(listname, "\${SLURM_ARRAY_TASK_ID}", i)
        JLD.save(this_listname, Dict("stats"=>these_stats))

        push!(results_file_list, joinpath(data_dir, string("results_$(i)_of_$(nb_batches).jld")))
    end

    cmd = `$submit_command $sbatchname`
    println("running $cmd ...")
    @time run(cmd)
    println("done")

    results_list_file = joinpath(data_dir, "results_list.txt")
    println("saving results file list to")
    println(results_list_file)
    open(results_list_file, "w") do opened
        for f in results_file_list
            write(opened, "$f\n")
        end
    end

    return results_file_list
end

function gather_results(results_file_list; save_file::Nullable=Nullable())
    results = JLD.load(first(results_file_list))
    stats = results["stats"]
    for (i,f) in enumerate(results_file_list[2:end])
        # print("\rmerging file $(i+1) of $(length(results_file_list))")
        # results = merge_results!(results, JLD.load(f), careful=false)
        tic()
        new_stats = JLD.load(f,"stats")
        t_load = toq()
        tic()
        stats = append!(stats, new_stats)
        t_append = toq()
        print(@sprintf("\rappended states %d of %d. t_load=%4.2f, t_append=%4.2f",
                       i+1,
                       length(results_file_list),
                       t_load,
                       t_append))
    end
    println()
    # results["stats"] = stats
    if !isnull(save_file)
        JLD.save(get(save_file), results)
    end
    return results
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
"""
    function from_jld2_file(filename::AbstractString,dataset::AbstractString="data")

Loads an object from a JLD2 file given its filename and, optionally, the dataset name.
The dataset specifies the name and the root location of the data inside the generated JLD2 file.
"""
function from_jld2_file(filename::AbstractString,dataset::AbstractString="data")
    load(filename)[dataset]
end

"""
    function from_jld2_file(::Type{T},filename::AbstractString,dataset::AbstractString="data") where T

Loads an object from a JLD2 file given its filename and, optionally, the dataset name.
The dataset specifies the name and the root location of the data inside the generated JLD2 file.
Checks if the returned object is of the expected `Type{T}`, if not return error.
"""
function from_jld2_file(::Type{T},filename::AbstractString,dataset::AbstractString="data") where T
    obj = from_jld2_file(filename,dataset)
    if isa(obj,T)
        return obj
    else
        error("Object from JLD2 file is of type ",typeof(obj)," instead of the expected ",T)
    end
end

"""
    function to_jld2_file(object,filename::AbstractString,dataset::AbstractString="data")

Stores an object to a JLD2 file given its filename and, optionally, the dataset name.
The dataset specifies the name and the root location of the data inside the JLD2 file.
"""
function to_jld2_file(object,filename::AbstractString,dataset::AbstractString="data")
    save(DataFormat{:JLD2}, filename, dataset, object)
end


text:
module CovarianceFunctions

export AbstractKernel, MercerKernel, StationaryKernel, IsotropicKernel, MulitKernel
export gramian, Gramian
export DotProductInput, IsotropicInput, GenericInput, InputTrait, input_trait

using LinearAlgebra
using SparseArrays
using StaticArrays
using FillArrays
using LazyArrays
using Base.Threads

using ForwardDiff
using DiffResults
using Functors
using SymEngine

using SpecialFunctions: gamma
using BesselK

using KroneckerProducts
using WoodburyFactorizations
using BlockFactorizations
using NearestNeighbors
using IterativeSolvers
using ToeplitzMatrices
using FFTW

# IDEA: AbstractKernel{T, IT}, where IT isa InputType
# then const IsotropicKernel{T} = AbstractKernel{T, IsotropicInput} ...
abstract type AbstractKernel{T} end
abstract type MercerKernel{T} <: AbstractKernel{T} end
abstract type StationaryKernel{T} <: MercerKernel{T} end
abstract type IsotropicKernel{T} <: StationaryKernel{T} end # ίσος + τρόπος (equal + way)

const default_tol = 1e-6 # default tolerance for matrix solves and products

# class of matrix-valued kernels for multi-output GPs
abstract type MultiKernel{T} <: AbstractKernel{T} end # MultiKernel
# compute the (i, j)th entry of k(x, y)
Base.getindex(K::MultiKernel, i, j) = (x, y) -> getindex(K(x, y), i, j)

function gramian end

# first, utility functions
include("util.jl")
include("givens.jl") # TODO: should be in base, makes givensAlgorithm differentiable
include("derivatives.jl")
include("lazy_linear_algebra.jl") # TODO: separate out into package
include("lazy_grid.jl")
include("toeplitz.jl") # special functions for toeplitz matrices -> put in ToeplitzMatrices.jl

# include all types of kernels
include("algebra.jl") # kernel operations
include("stationary.jl") # stationary mercer kernels
include("transformation.jl")
include("mercer.jl") # general mercer kernels

include("properties.jl")
include("gramian.jl")
include("sparse.jl")
include("barneshut.jl") # fast MVM algorithm for isotropic-kernel Gramians
include("taylor.jl")

# including multi-output kernels
include("gradient.jl")
include("gradient_algebra.jl")
include("hessian.jl")
# include("taylor_gradient.jl") # fast MVM algorithm for isotropic GradientKernel Gramians
include("separable.jl")

end # CovarianceFunctions


text:
# This file is a part of EuclidianNormalizingFlows.jl, licensed under the MIT License (MIT).

__precompile__(true)

"""
    EuclidianNormalizingFlows

Euclidian normalizing flows.
"""
module EuclidianNormalizingFlows

using LinearAlgebra
using Random
using Statistics

using ArgCheck
using ArraysOfArrays
using ChainRulesCore
using ChangesOfVariables
using Distributions
using DocStringExtensions
using ElasticArrays
using ForwardDiffPullbacks
using Functors
using InverseFunctions
using Optim
using Optimisers
using Parameters
using SpecialFunctions
using StatsBase
using ValueShapes

import Zygote
import ZygoteRules

using Distributions: log2π

import InverseFunctions.inverse
import ChainRulesCore.rrule
import ChangesOfVariables.with_logabsdet_jacobian

include("abstract_trafo.jl")
include("optimize_whitening.jl")
include("householder_trafo.jl")
include("scale_shift_trafo.jl")
include("center_stretch.jl")
include("johnson_trafo.jl")

end # module


text:
using DSGE, Test, JLD2, ModelConstructors, Dates
path = dirname(@__FILE__)

# Set up arguments
m = AnSchorfheide(testing = true)
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:forecast_horizons, 1)

system = JLD2.jldopen("$path/../reference/forecast_args.jld2","r") do file
    read(file, "system")
end
z0 = zeros(n_states_augmented(m))

# Read expected output
exp_states, exp_obs, exp_pseudo, exp_shocks =
    JLD2.jldopen("$path/../reference/forecast_out.jld2", "r") do file
        read(file, "exp_states"),
        read(file, "exp_obs"),
        read(file, "exp_pseudo"),
        read(file, "exp_shocks")
    end

# Without shocks
global states, obs, pseudo, shocks = forecast(m, system, z0; draw_shocks = false)

@testset "Testing forecasting without drawing shocks" begin
    @test @test_matrix_approx_eq exp_states states
    @test @test_matrix_approx_eq exp_obs    obs
    @test @test_matrix_approx_eq exp_pseudo pseudo
    @test @test_matrix_approx_eq exp_shocks shocks
end

# Supplying shocks
global states, obs, pseudo, shocks = forecast(m, system, z0; shocks = shocks)

@testset "Testing forecasting with pre-supplied shocks" begin
    @test @test_matrix_approx_eq exp_states states
    @test @test_matrix_approx_eq exp_obs    obs
    @test @test_matrix_approx_eq exp_pseudo pseudo
    @test @test_matrix_approx_eq exp_shocks shocks
end

# Draw normally distributed shocks
global states, obs, pseudo, shocks = forecast(m, system, z0; draw_shocks = true)

# Draw t-distributed shocks
m <= Setting(:forecast_tdist_shocks, true)
global states, obs, pseudo, shocks = forecast(m, system, z0; draw_shocks = true)
m <= Setting(:forecast_tdist_shocks, false)

# Enforce ZLB
ind_r = m.observables[:obs_nominalrate]
ind_r_sh = m.exogenous_shocks[:rm_sh]
zlb_value = forecast_zlb_value(m)
shocks = zeros(n_shocks_exogenous(m), forecast_horizons(m))
shocks[ind_r_sh, :] .= -10.

@testset "Ensure valid forecasting at the ZLB" begin
    global states, obs, pseudo, shocks = forecast(m, system, z0; shocks = shocks)
    @test all(x -> x < zlb_value, obs[ind_r, :])

    global states, obs, pseudo, shocks = forecast(m, system, z0; shocks = shocks, enforce_zlb = true)
    @test all(x -> abs(x - zlb_value) < 0.01, obs[ind_r, :])
    @test all(x -> x != -10.,                 shocks[ind_r_sh, :])
end

@testset "Test that forecasting and populating shocks under alternative policy works" begin
    m <= Setting(:regime_switching, true)
    setup_permanent_altpol!(m, AltPolicy(:taylor93, eqcond, solve))
    @test typeof(forecast(m, system, z0; draw_shocks = true)) == NTuple{4, Array{Float64, 2}}
    delete!(DSGE.get_settings(m), :regime_switching)
end

@testset "Enforce ZLB as a temporary alternative policy" begin

    # Set up model for forecast and permanent NGDP
    m = Model1002("ss10"; custom_settings = Dict{Symbol, Setting}(:add_altpolicy_pgap => Setting(:add_altpolicy_pgap, true)))
    m <= Setting(:date_forecast_start,  Date(2020, 6, 30))
    m <= Setting(:date_conditional_end, Date(2020, 6, 30))
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                                2 => Date(2020, 3, 31),
                                                3 => Date(2020, 6, 30)))
    m <= Setting(:forecast_horizons, 12)
    shocks = zeros(n_shocks_exogenous(m), forecast_horizons(m))
    shocks[m.exogenous_shocks[:b_sh], 1] = -0.8 # suppose massive negative shock to spreads -> MP should drop!
    m <= Setting(:gensys2, false)
    m <= Setting(:regime_switching, true)
    setup_regime_switching_inds!(m)
    m <= Setting(:pgap_value, 12.)
    m <= Setting(:pgap_type, :ngdp)
    setup_permanent_altpol!(m, DSGE.ngdp(); cond_type = :none)
    system = compute_system(m)
    z0 = zeros(n_states_augmented(m))

    # First test with unconditional
    ngdp_states, ngdp_obs, ngdp_pseudo, _ = forecast(m, system, z0; shocks = shocks, cond_type = :none)
    @test !all(ngdp_obs[m.observables[:obs_nominalrate], :] .> -1e-14)
    ngdp_states, ngdp_obs, ngdp_pseudo = forecast(m, z0, ngdp_states, ngdp_obs, ngdp_pseudo, shocks; cond_type = :none)
    @test all(ngdp_obs[m.observables[:obs_nominalrate], :] .> -1e-14)

    # Now test a conditional forecast with regime switching in the forecast
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                                2 => Date(2020, 3, 31),
                                                3 => Date(2020, 6, 30),
                                                4 => Date(2020, 9, 30)))
    m <= Setting(:forecast_horizons, 12)
    shocks = zeros(n_shocks_exogenous(m), forecast_horizons(m; cond_type = :full))
    shocks[m.exogenous_shocks[:b_sh], 1] = -0.8 # suppose massive negative shock to spreads -> MP should drop!
    m <= Setting(:replace_eqcond, false)
    m <= Setting(:gensys2, false)
    m <= Setting(:regime_switching, true)
    setup_regime_switching_inds!(m)
    m <= Setting(:pgap_value, 12.)
    m <= Setting(:pgap_type, :ngdp)
    setup_permanent_altpol!(m, DSGE.ngdp(); cond_type = :full)
    system = compute_system(m)
    z0 = zeros(n_states_augmented(m))

    ngdp_states, ngdp_obs, ngdp_pseudo, _ = forecast(m, system, z0; shocks = shocks, cond_type = :full)
    @test !all(ngdp_obs[m.observables[:obs_nominalrate], :] .> -1e-14)
    ngdp_states, ngdp_obs, ngdp_pseudo = forecast(m, z0, ngdp_states, ngdp_obs, ngdp_pseudo, shocks; cond_type = :full)
    @test all(ngdp_obs[m.observables[:obs_nominalrate], :] .> -1e-14)
end

nothing


text:
module AstroBase

using Reexport

@reexport using AstroTime

struct AstroException <: Exception
    msg::String
end

Base.showerror(io::IO, err::AstroException) = print(io, err.msg)

include("constants.jl")
include("stubs.jl")
include("util.jl")
include(joinpath("bodies", "Bodies.jl"))
include(joinpath("earth_attitude", "EarthAttitude.jl"))
include(joinpath("two_body", "TwoBody.jl"))
include(joinpath("ephemerides", "Ephemerides.jl"))
include(joinpath("n_body", "NBody.jl"))
include(joinpath("frames", "Frames.jl"))
include(joinpath("coords", "Coords.jl"))

@reexport using .Bodies
@reexport using .Ephemerides
@reexport using .EarthAttitude
@reexport using .TwoBody
@reexport using .Frames
@reexport using .Coords

end # module


text:
#=
HackerRank - Mean, median, mode
https://www.hackerrank.com/challenges/s10-basic-statistics/problem

=#
using Printf
n = parse(Int, readline(stdin))
array = sort(parse.(Int, split(readline(stdin))))
# Example from the problem
# array = parse.(Int, split("64630 11735 14216 99233 99233 14470 4978 73429 38120 51135 67060"))

n = size(array)[1]
# mean
mu = sum(array) / n
@printf("%.1f\n", mu)
# median
med = (n % 2 == 0) ? (array[n ÷ 2] + array[n ÷ 2 + 1]) / 2.0 : array[n ÷ 2 + 1]
@printf("%.1f\n", med)
# mode
count = Dict{Int, Int}()
for v in array
    if v in keys(count)
        count[v] += 1
    else
        count[v] = 1
    end
end
mode = sort(collect(count), lt=(l, r)-> l[2] == r[2] ? l[1] < r[1] : l[2] > r[2])[1][1]
@printf("%d\n", mode)


text:
@testset "simple" begin

    bc = rand(2,3)
    cd = rand(3,4)

    @matmul A[b,d] := sum(c) bc[b,c] * cd[c,d]
    @test A ≈ bc * cd

    B = similar(A)
    @matmul B[b,d] = sum(c) bc[b,c] * cd[c,d] # in-place
    @test B ≈ A

    ## matrix-vector
    c = randn(3)
    @matmul b[b] := sum(c) bc[b,c] * c[c]
    @test b ≈ bc * c

    @matmul c[c] := sum(b) b[b] * bc[b,c]
    @test c ≈ vec(b' * bc)

    b′ = similar(b)
    @matmul b′[b] = sum(c) bc[b,c] * c[c] # in-place
    @test b′ ≈ bc * c

    c′ = similar(c)
    @matmul c′[c] = sum(b) b[b] * bc[b,c]
    @test c′ ≈ c

    ## vector-vector
    @matmul z[] := sum(b) b[b] * b[b]
    @test z[] ≈ dot(b,b)
    @test z isa Array{Float64,0}

    z′ = similar(z)
    @matmul z′[] = sum(b) b[b] * b[b]
    @test z′ ≈ z

end
@testset "permutedims" begin

    bbb = rand(2,2,2);
    bbb2 = randn(2,2,2);

    @matmul A[i,j] := sum(k,k2) bbb[i,k,k2] * bbb2[j,k,k2]
    @reduce B[i,j] := sum(k,k2) bbb[i,k,k2] * bbb2[j,k,k2]
    @einsum C[i,j] := bbb[i,k,k2] * bbb2[j,k,k2]
    @test A ≈ B ≈ C

    D = similar(A);
    @matmul D[i,j] = sum(k,k2) bbb[i,k,k2] * bbb2[j,k,k2] # in-place
    @test A ≈ D

    @matmul A[i,j] := sum(k,k2) bbb[k2,i,k] * bbb2[j,k,k2]
    @reduce B[i,j] := sum(k,k2) bbb[k2,i,k] * bbb2[j,k,k2]
    @einsum C[i,j] := bbb[k2,i,k] * bbb2[j,k,k2]
    @test A ≈ B ≈ C

    D = similar(A);
    @matmul D[i,j] = sum(k,k2) bbb[k2,i,k] * bbb2[j,k,k2] # in-place
    @test A ≈ D

    cccc = rand(3,3,3,3);
    cccc2 = randn(3,3,3,3);

    @matmul A[i,j] := sum(x,y,z) cccc[x,y,j,z] * cccc2[z,i,y,x]
    @reduce B[i,j] := sum(x,y,z) cccc[x,y,j,z] * cccc2[z,i,y,x]
    @einsum C[i,j] := cccc[x,y,j,z] * cccc2[z,i,y,x]
    @test A ≈ B ≈ C

    D = similar(A);
    @matmul D[i,j] = sum(x,y,z) cccc[x,y,j,z] * cccc2[z,i,y,x]  # in-place
    @test A ≈ D

    @matmul A[i,j,k,l] := sum(x,y) cccc[i,x,y,k] * cccc2[l,j,y,x]
    @reduce B[i,j,k,l] := sum(x,y) cccc[i,x,y,k] * cccc2[l,j,y,x]
    @einsum C[i,j,k,l] := cccc[i,x,y,k] * cccc2[l,j,y,x]
    @test A ≈ B ≈ C

    D = similar(A);
    @matmul D[i,j,k,l] = sum(x,y) cccc[i,x,y,k] * cccc2[l,j,y,x]   # in-place
    @test A ≈ D

    @matmul A[i,j,k,l,m,n] := sum(x) cccc[x,k,i,m] * cccc2[l,x,j,n]
    @reduce B[i,j,k,l,m,n] := sum(x) cccc[x,k,i,m] * cccc2[l,x,j,n]
    @einsum C[i,j,k,l,m,n] := cccc[x,k,i,m] * cccc2[l,x,j,n]
    @test A ≈ B  ≈ C

    D = similar(A);
    @matmul D[i,j,k,l,m,n] = sum(x) cccc[x,k,i,m] * cccc2[l,x,j,n]  # in-place
    @test A ≈ D

    ## repeat with non-alphabeitcal LHS!
    @matmul A[j,i,l,k,n,m] := sum(x) cccc[x,k,i,m] * cccc2[l,x,j,n]
    @reduce B[j,i,l,k,n,m] := sum(x) cccc[x,k,i,m] * cccc2[l,x,j,n]
    @einsum C[j,i,l,k,n,m] := cccc[x,k,i,m] * cccc2[l,x,j,n]
    @test A ≈ B  ≈ C

    D = similar(A);
    @matmul D[j,i,l,k,n,m] = sum(x) cccc[x,k,i,m] * cccc2[l,x,j,n]  # in-place
    @test A ≈ D

    ## Mason Potter's example:
    # W = rand(2,2,2,2); M = rand(2,2,2);
    W = rand(6,7,4,5); M = rand(7,2,3);

    @reduce N[σ, b\a, b′\a′] := sum(σ′) W[σ,σ′,b,b′] * M[σ′,a,a′];
    @matmul N2[σ, b\a, b′\a′] := sum(σ′) W[σ,σ′,b,b′] * M[σ′,a,a′];
    @test N ≈ N2

    @reduce R[σ, b,a, b′\a′] := sum(σ′) W[σ,σ′,b,b′] * M[σ′,a,a′]  lazy;
    @matmul R2[σ, b,a, b′\a′] := sum(σ′) W[σ,σ′,b,b′] * M[σ′,a,a′];
    @test R ≈ R2
    # invperm((4,2,1,3)) == (3, 2, 4, 1)

    @reduce S[σ, b,a, b′,a′] := sum(σ′) W[σ,σ′,b,b′] * M[σ′,a,a′];
    @matmul S2[σ, b,a, b′,a′] := sum(σ′) W[σ,σ′,b,b′] * M[σ′,a,a′];
    @test S ≈ S2
    # invperm((4,2,1,3,5)) == (3, 2, 4, 1, 5)

    # ## in-place version
    N3 = similar(N); N4 = similar(N);
    @reduce N3[σ, b\a, b′\a′] = sum(σ′) W[σ,σ′,b,b′] * M[σ′,a,a′]  lazy;
    @matmul N4[σ, b\a, b′\a′] = sum(σ′) W[σ,σ′,b,b′] * M[σ′,a,a′];
    @test N ≈ N3 ≈ N4

    ## these were OK:
    M = rand(3,3,3);

    @reduce Q[σ, b\a, a′] := sum(σ′) M[σ,σ′,b] * M[σ′,a,a′];
    @matmul Q2[σ, b\a, a′] := sum(σ′) M[σ,σ′,b] * M[σ′,a,a′];
    @test Q ≈ Q2

    @reduce Q3[σ, b,a, a′] := sum(σ′) M[σ,σ′,b] * M[σ′,a,a′];
    @matmul    Q4[σ, b,a, a′] := sum(σ′) M[σ,σ′,b] * M[σ′,a,a′];
    @test Q3 ≈ Q4

end
@testset "in & out shapes" begin

    bc = rand(2,3)
    cd = rand(3,4)

    @matmul A[b,_,d] := sum(c) bc[b,c] * cd[c,d]
    @matmul A1[b,1,d] := sum(c) bc[b,c] * cd[c,d]
    @test size(A) == size(A1) == (2,1,4)

    B = similar(A)
    @matmul B[b,_,d] = sum(c) bc[b,c] * cd[c,d] # in-place
    @test B == A

    @matmul C[b] := sum(c) bc[b,c] * cd[c,3]
    @test C == bc * cd[:,3]

    D = similar(C)
    @matmul D[b] = sum(c) bc[b,c] * cd[c,3] # in-place
    @test D == C

end
@testset "recursion" begin

    A = rand(2,3); B = rand(3,4); C = rand(4,5);

    @einsum V[i,l] := A[i,j] * B[j,k] * C[k,l]
    @reduce V2[i,l] := sum(j,k) A[i,j] * B[j,k] * C[k,l]

    @reduce W[i,l] := sum(j) A[i,j] * @matmul [j,l] := sum(k) B[j,k] * C[k,l]

    # @matmul W2[i,l] := sum(j) A[i,j] * @matmul [j,l] := B[j,k] * C[k,l]
    # maybe I decided not to allow this for now.

    @test V ≈ V2 ≈ W #≈ W2

end


text:
struct MarkovModel
    background::Array{MixtureModel, 2}
    params::MixtureParams
    cut::MinCut
    function MarkovModel(G, kernels=5, smoothing=2.0)
        (n, m) = size(G)
        cut = MinCut(zeros(Float64, n, m),
                     VARIANCE_THRESHOLD,
                     smoothing * VARIANCE_THRESHOLD)
        M = [MixtureModel(kernels) for i in 1:n, j in 1:m]
        new(M, MixtureParams(), cut)
    end
end

function apply!(M::MarkovModel, image)
    update!(M.cut, bg_energy!.(M.background, M.params, image))
    increment_frames!(M.params)
    segment!(M.cut)
end

function apply_mrf!(cut, background, params, image, smoothing=2.0)
    energy = bg_energy!.(background, params, image)
    increment_frames!(params)
    reset!(cut, energy,
        VARIANCE_THRESHOLD,
        smoothing * VARIANCE_THRESHOLD
    )
    segment!(cut)
end


text:
import Base.show

export CustomParser, Quoted

abstract type AbstractToken{T} end
fieldtype(::AbstractToken{T}) where {T} = T
fieldtype(::Type{AbstractToken{T}}) where {T} = T
fieldtype(::Type{T}) where {T<:AbstractToken} = fieldtype(supertype(T))

"""
`tryparsenext{T}(tok::AbstractToken{T}, str, i, till, localopts)`

Parses the string `str` starting at position `i` and ending at or before position `till`. `localopts` is a [LocalOpts](@ref) object which contains contextual options for quoting and NA parsing. (see [LocalOpts](@ref) documentation)

`tryparsenext` returns a tuple `(result, nextpos)` where `result` is of type `Nullable{T}`, `Nullable{T}()` if parsing failed, non-null containing the parsed value if it succeeded. If parsing succeeded, `nextpos` is the position the next token, if any, starts at. If parsing failed, `nextpos` is the position at which the parsing failed.
"""
function tryparsenext end

## options passed down for tokens (specifically NAToken, StringToken)
## inside a Quoted token
"""
    LocalOpts

Options local to the token currently being parsed.
- `endchar`: Till where to parse. (e.g. delimiter or quote ending character)
- `spacedelim`: Treat spaces as delimiters
- `quotechar`: the quote character
- `escapechar`: char that escapes the quote
- `includequotes`: whether to include quotes while parsing
- `includenewlines`: whether to include newlines while parsing
"""
struct LocalOpts{T_ENDCHAR<:Union{Char,UInt8}, T_QUOTECHAR<:Union{Char,UInt8}, T_ESCAPECHAR<:Union{Char,UInt8}}
    endchar::T_ENDCHAR        # End parsing at this char
    spacedelim::Bool
    quotechar::T_QUOTECHAR       # Quote char
    escapechar::T_ESCAPECHAR      # Escape char
    includequotes::Bool   # Whether to include quotes in string parsing
    includenewlines::Bool # Whether to include newlines in string parsing
end

const default_opts = LocalOpts(UInt8(','), false, UInt8('"'), UInt8('"'), false, false)
# helper function for easy testing:
@inline function tryparsenext(tok::AbstractToken, str, opts::LocalOpts=default_opts)
    tryparsenext(tok, str, firstindex(str), lastindex(str), opts)
end

# fallback for tryparsenext methods which don't care about local opts
@inline function tryparsenext(tok::AbstractToken, str, i, len, locopts)
    tryparsenext(tok, str, i, len)
end

struct WrapLocalOpts{T, X<:AbstractToken} <: AbstractToken{T}
    opts::LocalOpts
    inner::X
end

WrapLocalOpts(opts, inner) = WrapLocalOpts{fieldtype(inner), typeof(inner)}(opts, inner)

@inline function tryparsenext(tok::WrapLocalOpts, str, i, len, opts::LocalOpts=default_opts)
    tryparsenext(tok.inner, str, i, len, tok.opts)
end

# needed for promoting guessses
struct Unknown <: AbstractToken{Missing} end
fromtype(::Type{Missing}) = Unknown()
const nullableNA = Nullable{Missing}(missing)
function tryparsenext(::Unknown, str, i, len, opts)
    nullableNA, i
end
show(io::IO, ::Unknown) = print(io, "<unknown>")
struct CustomParser{T, F} <: AbstractToken{T}
    f::Function
end

"""
    CustomParser(f, T)

Provide a custom parsing mechanism.

# Arguments:

- `f`: the parser function
- `T`: The type of the parsed value

The parser function must take the following arguments:
- `str`: the entire string being parsed
- `pos`: the position in the string at which to start parsing
- `len`: the length of the string the maximum position where to parse till
- `opts`: a [LocalOpts](@ref) object with options local to the current field.

The parser function must return a tuple of two values:

- `result`: A `Nullable{T}`. Set to `Nothing{T}()` if parsing must fail, containing the value otherwise.
- `nextpos`: If parsing succeeded this must be the next position after parsing finished, if it failed this must be the position at which parsing failed.
"""
CustomParser(f, T) = CustomParser{T,typeof(f)}(f)

show(io::IO, c::CustomParser{T}) where {T} = print(io, "{{custom:$T}}")

@inline function tryparsenext(c::CustomParser, str, i, len, opts)
    c.f(str, i, len, opts)
end

# Numberic parsing
"""
parse numbers of type T
"""
struct Numeric{T} <: AbstractToken{T}
    decimal::Char
    thousands::Char
end
show(io::IO, c::Numeric{T}) where {T} = print(io, "<$T>")

Numeric(::Type{T}, decimal='.', thousands=',') where {T} = Numeric{T}(decimal, thousands)
fromtype(::Type{N}) where {N<:Number} = Numeric(N)

### Unsigned integers

function tryparsenext(::Numeric{T}, str, i, len) where {T<:Signed}
    R = Nullable{T}
    @chk2 sign, i = tryparsenext_sign(str, i, len)
    @chk2 x, i = tryparsenext_base10(T, str, i, len)

    @label done
    return R(convert(T, sign*x)), i

    @label error
    return R(), i
end

@inline function tryparsenext(::Numeric{T}, str, i, len) where {T<:Unsigned}
    tryparsenext_base10(T,str, i, len)
end

@inline _is_e(str, i) = str[i]=='e' || str[i]=='E'

@inline _is_negative(str, i) = str[i]=='-'

@inline _is_positive(str, i) = str[i]=='+'

const pre_comp_exp_double = Double64[Double64(10.0)^i for i=0:308]

@inline function convert_to_double(f1::Int64, exp::Int)
    f = Float64(f1)
    r = f1 - Int64(f) # get the remainder
    x = Double64(f) + Double64(r)

    maxexp = 308
    minexp = -256

    if exp >= 0
        x *= pre_comp_exp_double[exp+1]
    else
        if exp < minexp # not sure why this is a good choice, but it seems to be!
            x /= pre_comp_exp_double[-minexp+1]
            x /= pre_comp_exp_double[-exp + minexp + 1]
        else
            x /= pre_comp_exp_double[-exp+1]
        end
    end
    return Float64(x)
end

@inline function tryparsenext(::Numeric{F}, str, i, len) where {F<:AbstractFloat}
    R = Nullable{F}

    y1 = iterate(str, i)
    y1===nothing && @goto error

    negate = false
    c = y1[1]
    if c=='-'
        negate = true
        i = y1[2]
    elseif c=='+'
        i = y1[2]
    end

    f1::Int64 = 0

    # read an integer up to the decimal point
    f1, rval1, idecpt = parse_uint_and_stop(str, i, len, f1)
    idecpt = read_digits(str, idecpt, len) # get any trailing digits
    i = idecpt

    ie = i
    frac_digits = 0

    # next thing must be dec pt.
    y2 = iterate(str, i)
    if y2!==nothing && y2[1]=='.'
        i =y2[2]
        f1, rval2, ie = parse_uint_and_stop(str, i, len, f1)
        # TODO This is incorrect for string types where a digit takes up
        # more than one codeunit, we need to return the number of digits
        # from parse_uint_and_stop instead. Ok for now because we are
        # not handling any such string types.
        frac_digits = ie - i

        ie = read_digits(str, ie, len) # get any trailing digits
    elseif !rval1 # no first number, and now no deciaml point => invalid
        @goto error
    end

    # Next thing must be exponent
    i = ie
    eval::Int32 = 0

    y3 = iterate(str, i)
    if y3!==nothing && _is_e(str, i)
        i = y3[2]

        y4 = iterate(str, i)
        if y4!==nothing
            enegate = false
            if _is_negative(str, i)
                enegate = true
                i = y4[2]
            elseif _is_positive(str, i)
                i = y4[2]
            end
        end
        eval, rval3, i = parse_uint_and_stop(str, i, len, eval)
        if enegate
            eval *= Int32(-1)
        end
    end

    exp = eval - frac_digits

    maxexp = 308
    minexp = -307

    if frac_digits <= 15 && -22 <= exp <= 22
        if exp >= 0
            f = F(f1)*10.0^exp
        else
            f = F(f1)/10.0^(-exp)
        end
    else
          f = convert_to_double(f1, exp)
    end

    if negate
        f = -f
    end

    @label done
    return R(convert(F, f)), i

    @label error
    return R(), i
end

struct Percentage <: AbstractToken{Float64}
end

const floatparser = Numeric(Float64)
function tryparsenext(::Percentage, str, i, len, opts)
    num, ii = tryparsenext(floatparser, str, i, len, opts)
    if isnull(num)
        return num, ii
    else
        # parse away the % char
        ii = eatwhitespaces(str, ii, len)
        y = iterate(str, ii)
        if y===nothing
            return Nullable{Float64}(), ii # failed to parse %
        else
            c = y[1]; k = y[2]
            if c != '%'
                return Nullable{Float64}(), ii # failed to parse %
            else
                return Nullable{Float64}(num.value / 100.0), k # the point after %
            end
        end
    end
end

"""
Parses string to the AbstractString type `T`. If `T` is `StrRange` returns a
`StrRange` with start position (`offset`) and `length` of the substring.
It is used internally by `csvparse` for avoiding allocating strings.
"""
struct StringToken{T} <: AbstractToken{T}
end

function StringToken(t::Type{T}) where T
    StringToken{T}()
end
show(io::IO, c::StringToken) = print(io, "<string>")

fromtype(::Type{S}) where {S<:AbstractString} = StringToken(S)

function tryparsenext(s::StringToken{T}, str, i, len, opts) where {T}
    inside_quoted_strong = Char(opts.endchar) == Char(opts.quotechar)
    escapecount = 0
    R = Nullable{T}
    p = ' '
    i0 = i
    if opts.includequotes
        y = iterate(str, i)
        if y!==nothing
            c = y[1]; ii = y[2]
            if c == Char(opts.quotechar)
                i = ii # advance counter so that
                       # the while loop doesn't react to opening quote
            end
        end
    end

    y2 = iterate(str, i)
    while y2!==nothing
        c = y2[1]; ii = y2[2]

        if inside_quoted_strong && p==Char(opts.escapechar)
            escapecount += 1
        end

        if opts.spacedelim && (c == ' ' || c == '\t')
            break
        elseif !opts.spacedelim && c == Char(opts.endchar)
            if inside_quoted_strong
                # this means we're inside a quoted string
                if Char(opts.quotechar) == Char(opts.escapechar)
                    # sometimes the quotechar is the escapechar
                    # in that case we need to see the next char
                    y3 = iterate(str, ii)
                    if y3===nothing
                        if opts.includequotes
                            i=ii
                        end
                        break
                    else
                        nxt = y3[1]; j = y3[2]
                        if nxt == Char(opts.quotechar)
                            # the current character is escaping the
                            # next one
                            i = j # skip next char as well
                            p = nxt
                            y2 = iterate(str, i)
                            continue
                        end
                    end
                elseif p == Char(opts.escapechar)
                    # previous char escaped this one
                    i = ii
                    p = c
                    y2 = iterate(str, i)
                    continue
                end
            end
            if opts.includequotes
                i = ii
            end
            break
        elseif (!opts.includenewlines && isnewline(c))
            break
        end
        i = ii
        p = c

        y2 = iterate(str, i)
    end

    return R(_substring(T, str, i0, i-1, escapecount, opts)), i
end

@inline function _substring(::Type{String}, str, i, j, escapecount, opts)
    if escapecount > 0
        buf = IOBuffer(sizehint=j-i+1-escapecount)
        cur_i = i
        c = str[cur_i]
        if opts.includequotes && c==Char(opts.quotechar)
            print(buf, c)
            cur_i = nextind(str, cur_i)
        end
        while cur_i <= j
            c = str[cur_i]
            if c == Char(opts.escapechar)
                next_i = nextind(str, cur_i)
                if next_i <= j && str[next_i] == Char(opts.quotechar)
                    print(buf, str[next_i])
                    cur_i = next_i
                else
                    print(buf, c)
                end
            else
                print(buf, c)
            end
            cur_i = nextind(str, cur_i)
        end
        return String(take!(buf))
    else
        return unsafe_string(pointer(str, i), j-i+1)
    end
end

@inline function _substring(::Type{T}, str, i, j, escapecount, opts) where {T<:SubString}
    escapecount > 0 && error("Not yet handled.")
    T(str, i, thisind(j))
end

fromtype(::Type{StrRange}) = StringToken(StrRange)

@inline function alloc_string(str, r::StrRange)
    unsafe_string(pointer(str, 1 + r.offset), r.length)
end

@inline function _substring(::Type{StrRange}, str, i, j, escapecount, opts)
    StrRange(i - 1, j - i + 1, escapecount)
end

@inline function _substring(::Type{<:WeakRefString}, str, i, j, escapecount, opts)
    escapecount > 0 && error("Not yet handled.")
    WeakRefString(convert(Ptr{UInt8}, pointer(str, i)), j - i + 1)
end

export Quoted

struct Quoted{T, S<:AbstractToken, T_QUOTECHAR<:Union{Char,UInt8}, T_ESCAPECHAR<:Union{Char,UInt8}} <: AbstractToken{T}
    inner::S
    required::Bool
    stripwhitespaces::Bool
    includequotes::Bool
    includenewlines::Bool
    quotechar::T_QUOTECHAR
    escapechar::T_ESCAPECHAR
end

function show(io::IO, q::Quoted)
    c = Char(q.quotechar)
    print(io, "$c")
    show(io, q.inner)
    print(io, "$c")
end

"""
`Quoted(inner::AbstractToken; <kwargs>...)`

# Arguments:
- `inner`: The token inside quotes to parse
- `required`: are quotes required for parsing to succeed? defaults to `false`
- `includequotes`: include the quotes in the output. Defaults to `false`
- `includenewlines`: include newlines that appear within quotes. Defaults to `true`
- `quotechar`: character to use to quote (default decided by `LocalOpts`)
- `escapechar`: character that escapes the quote char (default set by `LocalOpts`)
"""
function Quoted(inner::S,
    quotechar::T_QUOTECHAR, escapechar::T_ESCAPECHAR;
    required=false,
    stripwhitespaces=fieldtype(S)<:Number,
    includequotes=false,
    includenewlines=true) where {S<:AbstractToken,T_QUOTECHAR,T_ESCAPECHAR}

    T = fieldtype(S)
    Quoted{T,S,T_QUOTECHAR,T_ESCAPECHAR}(inner, required, stripwhitespaces, includequotes,
                includenewlines, quotechar, escapechar)
end

Quoted(t::Type, quotechar, escapechar; kwargs...) = Quoted(fromtype(t), quotechar, escapechar; kwargs...)

function tryparsenext(q::Quoted{T,S,T_QUOTECHAR,T_ESCAPECHAR}, str, i, len, opts) where {T,S,T_QUOTECHAR,T_ESCAPECHAR}
    y1 = iterate(str, i)
    if y1===nothing
        q.required && @goto error
        # check to see if inner thing is ok with an empty field
        @chk2 x, i = tryparsenext(q.inner, str, i, len, opts) error
        @goto done
    end
    c = y1[1]; ii = y1[2]
    quotestarted = false
    if Char(q.quotechar) == c
        quotestarted = true
        if !q.includequotes
            i = ii
        end

        if q.stripwhitespaces
            i = eatwhitespaces(str, i, len)
        end
    else
        q.required && @goto error
    end

    if quotestarted
        qopts = LocalOpts(q.quotechar, false, q.quotechar, q.escapechar,
                         q.includequotes, q.includenewlines)
        @chk2 x, i = tryparsenext(q.inner, str, i, len, qopts)
    else
        @chk2 x, i = tryparsenext(q.inner, str, i, len, opts)
    end

    if i > len
        if quotestarted && !q.includequotes
            @goto error
        end
        @goto done
    end

    if q.stripwhitespaces
        i = eatwhitespaces(str, i, len)
    end
    y2 = iterate(str, i)
    y2===nothing && error("Internal error.")
    c = y2[1]; ii = y2[2]

    if quotestarted && !q.includequotes
        c != Char(q.quotechar) && @goto error
        i = ii
    end

    @label done
    return Nullable{T}(x), i

    @label error
    return Nullable{T}(), i
end

## Date and Time
struct DateTimeToken{T,S<:DateFormat} <: AbstractToken{T}
    format::S
end

"""
    DateTimeToken(T, fmt::DateFormat)

Parse a date time string of format `fmt` into type `T` which is
either `Date`, `Time` or `DateTime`.
"""
DateTimeToken(T::Type, df::S) where {S<:DateFormat} = DateTimeToken{T, S}(df)
DateTimeToken(df::S) where {S<:DateFormat} = DateTimeToken{DateTime, S}(df)
fromtype(df::DateFormat) = DateTimeToken(DateTime, df)
fromtype(::Type{DateTime}) = DateTimeToken(DateTime, ISODateTimeFormat)
fromtype(::Type{Date}) = DateTimeToken(Date, ISODateFormat)

function tryparsenext(dt::DateTimeToken{T}, str, i, len, opts) where {T}
    R = Nullable{T}
    nt, i = tryparsenext_internal(T, str, i, len, dt.format, opts.endchar)
    if isnull(nt)
        return R(), i
    else
        return R(T(nt.value...)), i
    end
end

### Missing

const nastrings_upcase = ["NA", "NULL", "N/A","#N/A", "#N/A N/A", "#NA",
                          "-1.#IND", "-1.#QNAN", "-NaN", "-nan",
                          "1.#IND", "1.#QNAN", "N/A", "NA", "NaN", "nan"]

const NA_STRINGS = sort!(vcat(nastrings_upcase, map(lowercase, nastrings_upcase)))

struct NAToken{T, S<:AbstractToken} <: AbstractToken{T}
    inner::S
    emptyisna::Bool
    nastrings::Vector{String}
end

"""
`NAToken(inner::AbstractToken; options...)`

Parses a Nullable item.

# Arguments
- `inner`: the token to parse if non-null.
- `emptyisna`: should an empty item be considered NA? defaults to true
- `nastrings`: strings that are to be considered NA. Defaults to `$NA_STRINGS`
"""
function NAToken(
    inner::S
  ; emptyisna=true
  , nastrings=NA_STRINGS) where S

    T = fieldtype(inner)
    NAToken{UnionMissing{T}, S}(inner, emptyisna, nastrings)
end

function show(io::IO, na::NAToken)
    show(io, na.inner)
    print(io, "?")
end

function tryparsenext(na::NAToken{T}, str, i, len, opts) where {T}
    R = Nullable{T}
    i = eatwhitespaces(str, i, len)
    y1 = iterate(str,i)
    if y1===nothing
        if na.emptyisna
            @goto null
        else
            @goto error
        end
    end

    c = y1[1]; ii=y1[2]
    if (c == Char(opts.endchar) || isnewline(c)) && na.emptyisna
       @goto null
    end

    if isa(na.inner, Unknown)
        @goto maybe_null
    end
    @chk2 x,ii = tryparsenext(na.inner, str, i, len, opts) maybe_null

    @label done
    return R(convert(T, x)), ii

    @label maybe_null
    naopts = LocalOpts(opts.endchar, opts.spacedelim, opts.quotechar,
                       opts.escapechar, false, opts.includenewlines)
    @chk2 nastr, ii = tryparsenext(StringToken(WeakRefString{UInt8}), str, i, len, naopts)
    if !isempty(searchsorted(na.nastrings, nastr))
        i=ii
        i = eatwhitespaces(str, i, len)
        @goto null
    end
    return R(), i

    @label null
    return R(missing), i

    @label error
    return R(), i
end

fromtype(::Type{Union{Missing,T}}) where T = NAToken(fromtype(T))

struct SkipToken{S} <: AbstractToken{Nothing}
    inner::S
end

function tryparsenext(f::SkipToken, str, i, len, opts)
    x, ii = tryparsenext(f.inner, str, i, len, opts)

    if isnull(x)
        return Nullable{Nothing}(), ii
    else
        return Nullable{Nothing}(nothing), ii
    end
end

### Field parsing

abstract type AbstractField{T} <: AbstractToken{T} end # A rocord is a collection of abstract fields

struct Field{T,S<:AbstractToken} <: AbstractField{T}
    inner::S
    ignore_init_whitespace::Bool
    ignore_end_whitespace::Bool
    eoldelim::Bool
end

function Field(inner::S; ignore_init_whitespace=true, ignore_end_whitespace=true, eoldelim=false) where S
    T = fieldtype(inner)
    Field{T,S}(inner, ignore_init_whitespace, ignore_end_whitespace, eoldelim)
end

function Field(f::Field; inner=f.inner, ignore_init_whitespace=f.ignore_init_whitespace,
                  ignore_end_whitespace=f.ignore_end_whitespace,
                  eoldelim=f.eoldelim)
    T = fieldtype(inner)
    Field{T,typeof(inner)}(inner, ignore_init_whitespace,
                           ignore_end_whitespace, eoldelim)
end

function swapinner(f::Field, inner::AbstractToken;
        ignore_init_whitespace= f.ignore_end_whitespace
      , ignore_end_whitespace=f.ignore_end_whitespace
      , eoldelim=f.eoldelim
  )
    Field(inner;
        ignore_init_whitespace=ignore_end_whitespace
      , ignore_end_whitespace=ignore_end_whitespace
      , eoldelim=eoldelim
     )

end

function tryparsenext(f::Field{T}, str, i, len, opts) where {T}
    R = Nullable{T}
    i > len && @goto error
    if f.ignore_init_whitespace
        y1 = iterate(str, i)
        while y1!==nothing
            c = y1[1]; ii = y1[2]
            !isspace(c) && break
            i = ii
            y1 = iterate(str, i)
        end
    end
    @chk2 res, i = tryparsenext(f.inner, str, i, len, opts)

    if f.ignore_end_whitespace
        i0 = i
        y2 = iterate(str, i)
        while y2!==nothing
            c = y2[1]; ii = y2[2]
            !opts.spacedelim && Char(opts.endchar) == '\t' && c == '\t' && (i =ii; @goto done)
            !isspace(c) && c != '\t' && break
            i = ii
            y2 = iterate(str, i)
        end

        opts.spacedelim && i > i0 && @goto done
    end
    # todo don't ignore whitespace AND spacedelim

    if i > len
        if f.eoldelim
            @goto done
        else
            @goto error
        end
    end

    y3 = iterate(str, i)
    y3===nothing && error("Internal error.")
    c = y3[1]; ii = y3[2]
    opts.spacedelim && (isspace(c) || c == '\t') && (i=ii; @goto done)
    !opts.spacedelim && Char(opts.endchar) == c && (i=ii; @goto done)

    if f.eoldelim
        if c == '\r'
            i=ii
            y4 = iterate(str, i)
            if y4!==nothing
                c = y4[1]; ii = y4[2]
                if c == '\n'
                    i=ii
                end
            end
            @goto done
        elseif c == '\n'
            i=ii
            y5 = iterate(str, i)
            if y5!==nothing
                c = y5[1]; ii = y5[2]
                if c == '\r'
                    i=ii
                end
            end
            @goto done
        end
    end

    @label error
    return R(), i

    @label done
    return R(convert(T, res)), i
end


text:
using SimpleDiffEq, Test

@time begin
  @time @testset "Discrete Tests" begin include("discrete_tests.jl") end
  @time @testset "SimpleEM Tests" begin include("simpleem_tests.jl") end
  @time @testset "SimpleTsit5 Tests" begin include("simpletsit5_tests.jl") end
  @time @testset "SimpleATsit5 Tests" begin include("simpleatsit5_tests.jl") end
  @time @testset "GPUSimpleATsit5 Tests" begin include("gpusimpleatsit5_tests.jl") end
  @time @testset "SimpleRK4 Tests" begin include("simplerk4_tests.jl") end
end


text:
module PolyaGammaSamplers

using Distributions
using Random
using StatsFuns

export PolyaGammaPSWSampler

include("polyagammapswsampler.jl")

end