text:
using SpecialFunctions, RecursiveArrayTools, DifferentialEquations, Plots
using ConservationLawsParticles

# model
V1(t, x) = 1 + sin(x)/2
V2(t, x) = -1 - cos(x)/2
Wₐ′(t,x) = sign(x) / (abs(x) + 1) + x^3/20
Wᵣ(t, x) = 1 / (abs(x) + 1)
mob1(ρ, σ) = max(1 - ρ - σ/2, 0)
mob2(ρ, σ) = max(1 - ρ/2 - σ, 0)
attr = SampledInteraction(Wₐ′)
rep = IntegratedInteraction(Wᵣ)

model = ParabolicModel(
    (V1, V2),
    ((attr, rep), (rep, attr)),
    (mob1, mob2),
    (SimpleDiffusion(1/16), SimpleDiffusion(1/8)))

mmodel = ParabolicModel(
    (V1, V2),
    ((attr, rep), (rep, attr)),
    (mob1, mob2),
    (MinDiffusion(1/16), MinDiffusion(1/8)))

# initial condition
n = 80
x0 = ArrayPartition(
    vcat(range(-5, -3, length=n÷2), range(-2, -1, length=n÷2)),
    vcat(range(1, 2, length=n÷2), range(2.5, 3, length=n÷4), range(3.5, 4, length=n÷4)))

# time span
tspan = (0., 12.0)

# ODE system for the particles
prob = ODEProblem(abstract_velocities_gen!, x0, tspan, model)
mprob = ODEProblem(abstract_velocities_gen!, x0, tspan, mmodel)

# solve it
abstol = reltol = 1e-7
@time sol = solve(prob, BS5(); abstol=abstol, reltol=reltol);
@time msol = solve(mprob, BS5(); abstol=abstol, reltol=reltol);

# plot the particle trajectories
plot(legend=false, title="Trajectories", xlabel="time", ylabel="position")
plot!(sol, vars=1:1:n; color=1)
plot!(sol, vars=n+1:1:2n; color=2)
plot!(msol, vars=1:1:n; color=1, ls=:dot)
plot!(msol, vars=n+1:1:2n; color=2, ls=:dot)
savefig("multi-spec-2.png")

# plot an animation of the density
anim = @animate for t in range(tspan...; length=100)
    plot(title="Density", xlabel="position", ylabel="density",
        xlims=(-6,5), ylims=(0,1), legend=false)
    densityplot!(sol(t).x[1]; color=1)
    densityplot!(sol(t).x[2]; color=2)
    densityplot!(msol(t).x[1]; color=1, ls=:dot)
    densityplot!(msol(t).x[2]; color=2, ls=:dot)
end
gif(anim, "multi-spec-2.gif")


text:
"""
Basic representation of the state of the age structured model:
u[wa_index,age_group,disease state]

This is row major unpacked so
first 1,...,n_wa entries are 0-4 year old susceptibles in wider areas 1,...,n_wa
then n_wa+1,...,2n_wa are 5-9 year old susceptibles in wider areas 1,...,n_wa
.
.
.
then n_wa*n_s + 1, ...,  n_wa*n_s + n_wa entries are 0-4 year old exposed in wider areas 1,...,n_wa
.
.
.

States:
1 -> S(usceptibles)
2 -> E(xposed/latent infected)
3 -> P(re-symptomatic infected)
4 -> A(symptomatic)
5 -> M(ild) symptomatics
6 -> first mild then eventually (se)V(ere) symptomatics
7 -> H(ospitalised)
8 -> Recovered
9 -> Cumulative P->A
10-> Cumulative P->M
11-> Cumulative P->V
12 -> Cumulative V->H

Events for each wider area and age group:

1-> S to E
2-> E to P
3-> P to A
4-> P to M
5-> P to V
6-> V to H
7-> M to R
8-> A to R

"""

dc_age = zeros(Int64,n_wa*n_a*n_s,n_ta*n*n_a)

function import_rate_mom(t,into_mom,global_prev)
    if t+1>min(length(into_mom),length(global_prev))
        t_int=min(length(into_mom),length(global_prev))
    else
        t_int=Int(floor(t))+1
    end
    return into_mom[t_int]*global_prev[t_int]
end

function import_rate_nai(t,into_nai,global_prev)
    if t+1>min(length(into_nai),length(global_prev))
        t_int=min(length(into_nai),length(global_prev))
    else
        t_int=Int(floor(t))+1
    end
    return into_nai[t_int]*global_prev[t_int]
end

# asymp_indices = zeros(Bool,n_wa,n_a,n_s)
# asymp_indices[:,:,3] .= true;
# f_asymp_indices = findall(asymp_indices[:])
# diseased_indices = zeros(Bool,n_wa,n_a,n_s)
# diseased_indices[:,:,4] .= true;
# f_diseased_indices = findall(diseased_indices[:])
# asymp_indices = 0;#free memory
# diseased_indices = 0;
"""
    function calculate_infection_rates!(u,p::CoVParameters_AS,t)

Inplace method for calculating the force of infection on susceptibles in each spatial and age group.
"""
function calculate_infection_rates!(u,p::CoVParameters_AS,t)
    I_P = @view u[:,:,3]
    I_A = @view u[:,:,4]
    I_M = @view u[:,:,5]
    I_V = @view u[:,:,6]
    mul!(p.Î,p.T,p.ϵ*I_P .+ p.ϵ*I_A .+ p.ϵ_D*I_M .+ p.ϵ_V*I_V)  #Local infecteds **if** everyone moved around
    p.Î[:,immobile_age_indices] .= p.ϵ*I_P[:,immobile_age_indices] .+ p.ϵ*I_A[:,immobile_age_indices] .+ p.ϵ_D*I_M[:,immobile_age_indices] .+ p.ϵ_V*I_V[:,immobile_age_indices]#This corrects for immobility
    mul!(p.λ_loc,p.β*p.c_t(t).*(p.Î ./p.N̂),p.M)#Local force of infection due to age-mixing --- M is in to (row), from (col) format
    mul!(p.λ,p.T',p.λ_loc)#this accounts for mobile susceptibles contracting away from home
    p.λ[:,immobile_age_indices] .= p.λ_loc[:,immobile_age_indices]#This corrects for immobility of susceptibles
    p.λ[ind_mombasa_as,:] .+= p.ext_inf_rate*import_rate_mom(t,p.into_mom,p.global_prev)
    p.λ[ind_nairobi_as,:] .+= p.ext_inf_rate*import_rate_nai(t,p.into_nai,p.global_prev)
    return nothing
end

"""
    function rates(out,u,p::CoVParameters_AS,t)

Inplace method for calculating the rate of each of the 8 events per spatial and age group.
    1-> S to E
    2-> E to P
    3-> P to A
    4-> P to M
    5-> P to V
    6-> V to H
    7-> M to R
    8-> A to R
The out vector is in linear index form.
"""
function rates(out,u,p::CoVParameters_AS,t)
    @unpack λ,γ,σ₁,σ₂,δ,τ,μ₁,χ,rel_detection_rate,clear_quarantine,hₐ = p
    calculate_infection_rates!(u,p,t)
    for k = 1:length(out)
        i,a,eventtype = Tuple(index_as_events[k])
        if eventtype ==1
            out[k] = χ[a]*λ[i,a]*u[i,a,1] #Transmission #S->E
        end
        if eventtype ==2
            out[k] = σ₁*u[i,a,2] #E->P
        end
        if eventtype ==3
            out[k] = σ₂*(1-(δ*rel_detection_rate[a]))*u[i,a,3] #P->A
        end
        if eventtype ==4
            out[k] = σ₂*δ*(1-hₐ[a])*rel_detection_rate[a]*u[i,a,3] #P->M
        end
        if eventtype ==5
            out[k] = σ₂*δ*hₐ[a]*rel_detection_rate[a]*u[i,a,3] #P->V
        end
        if eventtype ==6
            out[k] = τ*u[i,a,6] # V->H
        end
        if eventtype ==7
            out[k] = γ*u[i,a,5] # M->R
        end
        if eventtype ==8
            out[k] = γ*u[i,a,4] # A->R
        end
    end
end

"""
    function change_matrix(dc)

Inplace method for creating a change matrix, that is a matrix encoding the linear map between frequency of events
    occuring and the change this causes in the state of the epidemic. Both the list of events and the list of states
    are encoded in a linear index.
"""
function change_matrix(dc)
    d1,d2 = size(dc)
    for k = 1:d2 #loop over
        i,a,eventtype = Tuple(index_as_events[k]) #
        if eventtype ==1 #S->E
            ind_S = linear_as[i,a,1]
            ind_E = linear_as[i,a,2]
            dc[ind_S,k] = -1
            dc[ind_E,k] = 1
        end
        if eventtype ==2 #E->P
            ind_E = linear_as[i,a,2]
            ind_P = linear_as[i,a,3]
            dc[ind_E,k] = -1
            dc[ind_P,k] = 1
        end
        if eventtype ==3 # P->A
            ind_P = linear_as[i,a,3]
            ind_A = linear_as[i,a,4]
            ind_cumA = linear_as[i,a,9]
            dc[ind_P,k] = -1
            dc[ind_A,k] = 1
            dc[ind_cumA,k] = 1
        end
        if eventtype ==4 # P->M
            ind_P = linear_as[i,a,3]
            ind_M = linear_as[i,a,5]
            ind_cumM = linear_as[i,a,10]
            dc[ind_P,k] = -1
            dc[ind_M,k] = 1
            dc[ind_cumM,k] = 1
        end
        if eventtype ==5 # P->V
            ind_P = linear_as[i,a,3]
            ind_V = linear_as[i,a,6]
            ind_cumV = linear_as[i,a,11]
            dc[ind_P,k] = -1
            dc[ind_V,k] = 1
            dc[ind_cumV,k] = 1
        end
        if eventtype ==6# V->H
            ind_V = linear_as[i,a,6]
            ind_H = linear_as[i,a,7]
            ind_cumVH = linear_as[i,a,12]
            dc[ind_V,k] = -1
            dc[ind_H,k] = 1
            dc[ind_cumVH,k] = 1
        end
        if eventtype ==7 # M->R
            ind_M = linear_as[i,a,5]
            ind_R = linear_as[i,a,8]
            dc[ind_M,k] = -1
            dc[ind_R,k] = 1
        end
        if eventtype ==8 # A->R
            ind_A = linear_as[i,a,4]
            ind_R = linear_as[i,a,8]
            dc[ind_A,k] = -1
            dc[ind_R,k] = 1
        end
    end
end

"""
    function PP_drivers(dN::Vector{Int64},rates,p)
This method in-place generates the crude number of each type of event proposed by a Poisson process with
rates calculated by rates(out,u,p::CoVParameters_AS,t).
"""
function PP_drivers(dN::Vector{Int64},rates,p)
    for i = 1:length(dN)
        if rates[i] >= 0.
            dN[i] = rand(Poisson(p.dt*rates[i]))
        else
            dN[i] = 0
        end
    end
end

"""
    function max_change(out,u,p::CoVParameters_AS)
This method in-place modifies the number of each type of event proposed by the Poisson process
    so that non-negativity is respected.
"""
function max_change(out,u,p::CoVParameters_AS)
    @unpack δ,rel_detection_rate,hₐ = p
    for i = 1:n_wa,a = 1:n_a
        ind_trans = linear_as_events[i,a,1]
        ind_EP = linear_as_events[i,a,2]
        ind_PA = linear_as_events[i,a,3]
        ind_PM = linear_as_events[i,a,4]
        ind_PV = linear_as_events[i,a,5]
        ind_VH = linear_as_events[i,a,6]
        ind_MR = linear_as_events[i,a,7]
        ind_AR = linear_as_events[i,a,8]

        out[ind_trans] = min(out[ind_trans],u[i,a,1])
        out[ind_EP] = min(out[ind_EP],u[i,a,2])
        out[ind_VH] = min(out[ind_VH],u[i,a,6])
        out[ind_MR] = min(out[ind_MR],u[i,a,5])
        out[ind_AR] = min(out[ind_AR],u[i,a,4])

        #splitting events using multinomial sampling
        if out[ind_PA] + out[ind_PM] + out[ind_PV] > u[i,a,3] #More transitions than actual P so all P individuals transition
            rel_rate_each_event = [1-(δ*rel_detection_rate[a]),δ*(1-hₐ[a])*rel_detection_rate[a],δ*hₐ[a]*rel_detection_rate[a]]
            D = rand(Multinomial(u[i, a, 3], LinearAlgebra.normalize!(rel_rate_each_event,1) )) #Draw the states of the P individuals after transition
            out[ind_PA] = D[1]
            out[ind_PM] = D[2]
            out[ind_PV] = D[3]
        end
    end
end

"""
    nonneg_tauleap(du,u,p::CoVParameters_AS,t)

This performs one in-place simulation of a time step
"""
function nonneg_tauleap(du,u,p::CoVParameters_AS,t)
    @unpack dc,dN,poi_rates,du_linear = p
    rates(poi_rates,u,p,t) #calculate rates of underlying Poisson processes
    PP_drivers(dN,poi_rates,p)#Generate Poisson rvs with rates scaled by time step dt
    max_change(dN,u,p)#Cap the size of the Poisson rvs to maintain non-negativity
    mul!(du_linear,dc,dN)#Calculates the effect on the state in the inplace du vector
    du .= reshape(du_linear,n_wa,n_a,n_s)
    du .+= u #Calculates how the state should change
end

"""
    function ode_model(du,u,p::CoVParameters_AS,t)

This is the vector field of the related KenyaCoV ODE model
"""
function ode_model(du,u,p::CoVParameters_AS,t)
    @unpack λ,γ,σ₁,σ₂,δ,τ,μ₁,χ,rel_detection_rate,hₐ = p
    calculate_infection_rates!(u,p,t)
    for i = 1:n_wa,a = 1:n_a
        du[i,a,1] = (-1)*χ[a]*λ[i,a]*u[i,a,1]
        du[i,a,2] = χ[a]*λ[i,a]*u[i,a,1] - σ₁*u[i,a,2]
        du[i,a,3] = σ₁*u[i,a,2] - σ₂*u[i,a,3]
        du[i,a,4] = σ₂*(1-(δ*rel_detection_rate[a]))*u[i,a,3] - γ*u[i,a,4]
        du[i,a,5] = σ₂*δ*(1-hₐ[a])*rel_detection_rate[a]*u[i,a,3] - γ*u[i,a,5]
        du[i,a,6] = σ₂*δ*hₐ[a]*rel_detection_rate[a]*u[i,a,3] - τ*u[i,a,6]
        du[i,a,7] = τ*u[i,a,6]
        du[i,a,8] = γ*u[i,a,5] + γ*u[i,a,4]
        du[i,a,9] = σ₂*(1-(δ*rel_detection_rate[a]))*u[i,a,3]
        du[i,a,10] = σ₂*δ*(1-hₐ[a])*rel_detection_rate[a]*u[i,a,3]
        du[i,a,11] = σ₂*δ*hₐ[a]*rel_detection_rate[a]*u[i,a,3]
        du[i,a,12] = τ*u[i,a,6] #Same as being hosp. for this version
     end
end

ext: jl


text:
# Copyright (c) 2013: Iain Dunning, Miles Lubin, and contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module Ipopt

import Ipopt_jll
import MathOptInterface

const MOI = MathOptInterface

function __init__()
    global libipopt = Ipopt_jll.libipopt
    return
end

include("C_wrapper.jl")
include("MOI_wrapper.jl")

export IpoptProblem,
    CreateIpoptProblem,
    FreeIpoptProblem,
    AddIpoptStrOption,
    AddIpoptNumOption,
    AddIpoptIntOption,
    OpenIpoptOutputFile,
    SetIpoptProblemScaling,
    SetIntermediateCallback,
    IpoptSolve

end


text:
module WeakRefStrings

export WeakRefString, WeakRefStringArray, StringArray, StringVector

using Missings

########################################################################
# WeakRefString
########################################################################

"""
A custom "weakref" string type that only points to external string data.
Allows for the creation of a "string" instance without copying data,
which allows for more efficient string parsing/movement in certain data processing tasks.

**Please note that no original reference is kept to the parent string/memory, so `WeakRefString` becomes unsafe
once the parent object goes out of scope (i.e. loses a reference to it)**

Internally, a `WeakRefString{T}` holds:

  * `ptr::Ptr{T}`: a pointer to the string data (code unit size is parameterized on `T`)
  * `len::Int`: the number of code units in the string data

See also [`WeakRefStringArray`](@ref)
"""
struct WeakRefString{T} <: AbstractString
    ptr::Ptr{T}
    len::Int # of code units
end

WeakRefString(ptr::Ptr{T}, len) where {T} = WeakRefString(ptr, Int(len))
WeakRefString(t::Tuple{Ptr{T}, Int}) where {T} = WeakRefString(t...)
WeakRefString(a::Vector{T}) where T <: Union{UInt8,UInt16,UInt32} = WeakRefString(pointer(a), length(a))

const NULLSTRING = WeakRefString(Ptr{UInt8}(0), 0)
const NULLSTRING16 = WeakRefString(Ptr{UInt16}(0), 0)
const NULLSTRING32 = WeakRefString(Ptr{UInt32}(0), 0)
Base.zero(::Type{WeakRefString{T}}) where {T} = WeakRefString(Ptr{T}(0), 0)

import Base: ==
function ==(x::WeakRefString{T}, y::WeakRefString{T}) where {T}
    x.len == y.len && (x.ptr == y.ptr || ccall(:memcmp, Cint, (Ptr{T}, Ptr{T}, Csize_t),
                                           x.ptr, y.ptr, x.len) == 0)
end
function ==(x::String, y::WeakRefString{T}) where {T}
    sizeof(x) == y.len && ccall(:memcmp, Cint, (Ptr{T}, Ptr{T}, Csize_t),
                                 pointer(x), y.ptr, y.len) == 0
end
==(y::WeakRefString, x::String) = x == y

function Base.hash(s::WeakRefString{T}, h::UInt) where {T}
    h += Base.memhash_seed
    ccall(Base.memhash, UInt, (Ptr{T}, Csize_t, UInt32), s.ptr, s.len, h % UInt32) + h
end

function Base.show(io::IO, x::WeakRefString{T}) where {T}
    print(io, '"')
    print(io, string(x))
    print(io, '"')
    return
end
Base.print(io::IO, s::WeakRefString) = print(io, string(s))
Base.textwidth(s::WeakRefString) = textwidth(string(s))

chompnull(x::WeakRefString{T}) where {T} = unsafe_load(x.ptr, x.len) == T(0) ? x.len - 1 : x.len

Base.string(x::WeakRefString) = x == NULLSTRING ? "" : unsafe_string(x.ptr, x.len)
Base.string(x::WeakRefString{UInt16}) = x == NULLSTRING16 ? "" : String(transcode(UInt8, unsafe_wrap(Array, x.ptr, chompnull(x))))
Base.string(x::WeakRefString{UInt32}) = x == NULLSTRING32 ? "" : String(transcode(UInt8, unsafe_wrap(Array, x.ptr, chompnull(x))))

Base.convert(::Type{WeakRefString{UInt8}}, x::String) = WeakRefString(pointer(x), sizeof(x))
Base.convert(::Type{String}, x::WeakRefString) = convert(String, string(x))
Base.String(x::WeakRefString) = string(x)
Base.Symbol(x::WeakRefString{UInt8}) = ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int), x.ptr, x.len)

Base.pointer(s::WeakRefString) = s.ptr
Base.pointer(s::WeakRefString, i::Integer) = s.ptr + i - 1

Base.ncodeunits(s::WeakRefString{T}) where {T} = s.len
Base.codeunit(s::WeakRefString{T}) where {T} = T

@inline function Base.codeunit(s::WeakRefString, i::Integer)
    @boundscheck checkbounds(s, i)
    GC.@preserve s unsafe_load(pointer(s, i))
end

Base.thisind(s::WeakRefString, i::Int) = Base._thisind_str(s, i)
Base.nextind(s::WeakRefString, i::Int) = Base._nextind_str(s, i)
Base.isvalid(s::WeakRefString, i::Int) = checkbounds(Bool, s, i) && thisind(s, i) == i

Base.@propagate_inbounds function Base.iterate(s::WeakRefString, i::Int=firstindex(s))
    i > ncodeunits(s) && return nothing
    b = codeunit(s, i)
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
    return Base.next_continued(s, i, u)
end

function Base.next_continued(s::WeakRefString, i::Int, u::UInt32)
    u < 0xc0000000 && (i += 1; @goto ret)
    n = ncodeunits(s)
    # first continuation byte
    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); i += 1
@label ret
    return reinterpret(Char, u), i
end

########################################################################
# WeakRefStringArray
########################################################################

init(::Type{T}, rows) where {T} = fill(zero(T), rows)
init(::Type{Union{Missing, T}}, rows) where {T} = Vector{Union{Missing, T}}(undef, rows)

"""
A [`WeakRefString`](@ref) container.
Holds the "strong" references to the data pointed by its strings, ensuring that
the referenced memory blocks stay valid during `WeakRefStringArray` lifetime.

Upon indexing an elemnt in a `WeakRefStringArray`, the underlying `WeakRefString` is converted to a proper
Julia `String` type by copying the memory; this ensures safe string processing in the general case. If additional
optimizations are desired, the direct `WeakRefString` elements can be accessed by indexing `A.elements`, where
`A` is a `WeakRefStringArray`.
"""
struct WeakRefStringArray{T<:WeakRefString, N, U} <: AbstractArray{Union{String, U}, N}
    data::Vector{Any}
    elements::Array{Union{T, U}, N}

    WeakRefStringArray(data::Vector{Any}, A::Array{Union{T, Missing}, N}) where {T <: WeakRefString, N} =
        new{T, N, Missing}(data, A)
    WeakRefStringArray(data::Vector{Any}, A::Array{T, N}) where {T <: WeakRefString, N} =
        new{T, N, Union{}}(data, A)
end

WeakRefStringArray(data, ::Type{T}, rows::Integer) where {T} = WeakRefStringArray(Any[data], init(T, rows))
WeakRefStringArray(data, A::Array{T}) where {T <: Union{WeakRefString, Missing}} = WeakRefStringArray(Any[data], A)

wk(A, B::AbstractArray) = WeakRefStringArray(A.data, B)
wk(A, w::WeakRefString) = string(w)
wk(A, ::Missing) = missing

Base.size(A::WeakRefStringArray) = size(A.elements)
Base.getindex(A::WeakRefStringArray, I...) = wk(A, getindex(A.elements, I...))
Base.setindex!(A::WeakRefStringArray{T, N}, v::Missing, i::Int) where {T, N} = setindex!(A.elements, v, i)
Base.setindex!(A::WeakRefStringArray{T, N}, v::Missing, I::Vararg{Int, N}) where {T, N} = setindex!(A.elements, v, I...)
Base.setindex!(A::WeakRefStringArray{T, N}, v::WeakRefString, i::Int) where {T, N} = setindex!(A.elements, v, i)
Base.setindex!(A::WeakRefStringArray{T, N}, v::WeakRefString, I::Vararg{Int, N}) where {T, N} = setindex!(A.elements, v, I...)
Base.setindex!(A::WeakRefStringArray{T, N}, v::String, i::Int) where {T, N} = (push!(A.data, codeunits(v)); setindex!(A.elements, v, i))
Base.setindex!(A::WeakRefStringArray{T, N}, v::String, I::Vararg{Int, N}) where {T, N} = (push!(A.data, codeunits(v)); setindex!(A.elements, v, I...))
Base.resize!(A::WeakRefStringArray, i) = resize!(A.elements, i)

Base.push!(a::WeakRefStringArray{T, 1}, v::Missing) where {T} = (push!(a.elements, v); a)
Base.push!(a::WeakRefStringArray{T, 1}, v::WeakRefString) where {T} = (push!(a.elements, v); a)
function Base.push!(A::WeakRefStringArray{T, 1}, v::String) where T
    push!(A.data, codeunits(v))
    push!(A.elements, v)
    return A
end

function Base.append!(a::WeakRefStringArray{T, 1}, b::WeakRefStringArray{T, 1}) where {T}
    append!(a.data, b.data)
    append!(a.elements, b.elements)
    return a
end

function Base.vcat(a::WeakRefStringArray{T, 1}, b::WeakRefStringArray{T, 1}) where T
    WeakRefStringArray(Any[a.data, b.data], vcat(a.elements, b.elements))
end

Base.deleteat!(a::WeakRefStringArray{T, 1}, inds) where {T} = (deleteat!(a.elements, inds); return a)

########################################################################
# StringArray
########################################################################

const STR = Union{Missing, <:AbstractString}

"""
`StringArray{T,N}`

Efficient storage for N dimensional array of strings.

`StringArray` stores underlying string data for all elements of the array
in a single contiguous buffer. It maintains offset and length for each
element.

`T` can be `String`, `WeakRefString`, `Union{Missing, String}` or
`Union{Missing, WeakRefString}`. `getindex` will return this type although
all variants have the same storage format.

You can use `convert(StringArray{U}, ::StringArray{T})` to change the
element type (e.g. to `WeakRefString` for efficiency) without copying
the data.

# Example construction

```
julia> sa = StringArray(["x", "y"]) # from Array{String}
2-element WeakRefStrings.StringArray{String,1}:
 "x"
 "y"

julia> sa = StringArray{WeakRefString}(["x", "y"])
2-element WeakRefStrings.StringArray{WeakRefStrings.WeakRefString,1}:
 "x"
 "y"

julia> sa = StringArray{Union{Missing, String}}(["x", "y"]) # with Missing
2-element WeakRefStrings.StringArray{Union{Missings.Missing, String},1}:
 "x"
 "y"

julia> sa = StringArray{Union{Missing, String}}(2,2) # undef
2×2 WeakRefStrings.StringArray{Union{Missings.Missing, String},2}:
 #undef  #undef
 #undef  #undef
```
"""
struct StringArray{T, N} <: AbstractArray{T, N}
    buffer::Vector{UInt8}
    offsets::Array{UInt64, N}
    lengths::Array{UInt32, N}
end

const StringVector{T} = StringArray{T, 1}

const UNDEF_OFFSET = typemax(UInt64)
const MISSING_OFFSET = typemax(UInt64)-1

Base.size(a::StringArray) = size(a.offsets)
Base.IndexStyle(::Type{<:StringVector}) = IndexLinear()

# no-copy convert between eltypes
"""
`convert(StringArray{U}, A::StringArray{T})`

convert `A` to StringArray of another element type (`U`) without
copying the underlying data.
"""
function Base.convert(::Type{<:StringArray{T}}, x::StringArray{<:STR,N}) where {T, N}
    StringArray{T, ndims(x)}(x.buffer, x.offsets, x.lengths)
end

function StringArray{T, N}(dims::Tuple{Vararg{Integer}}) where {T,N}
    StringArray{T,N}(similar(Vector{UInt8}, 0), fill(UNDEF_OFFSET, dims), fill(zero(UInt32), dims))
end

function StringArray{T}(dims::Tuple{Vararg{Integer}}) where {T}
    StringArray{T,length(dims)}(dims)
end

function StringArray(dims::Tuple{Vararg{Integer}})
    StringArray{String}(dims)
end

(::Type{S})(dims::Vararg{Integer,N}) where {S<:StringArray{T,N}} where {T,N} = StringArray{T,N}(dims)
(::Type{<:StringArray})(dims::Integer...) = StringArray{String,length(dims)}(dims)

function Base.convert(::Type{<:StringArray{T}}, arr::AbstractArray{<:STR, N}) where {T,N}
    s = StringArray{T, N}(size(arr))
    @inbounds for i in eachindex(arr)
        if _isassigned(arr, i)
            s[i] = arr[i]
        else
            s.offsets[i] = UNDEF_OFFSET
            s.lengths[i] = 0
        end
    end
    s
end
Base.convert(::Type{StringArray}, arr::AbstractArray{T}) where {T<:STR} = StringArray{T}(arr)
Base.convert(::Type{StringArray{T, N} where T}, arr::AbstractArray{S}) where {S<:STR, N} = StringVector{S}(arr)
StringVector{T}() where {T} = StringVector{T}(Vector{UInt8}(0), UInt64[], UInt32[])
StringVector() = StringVector{String}()

(T::Type{<:StringArray})(arr::AbstractArray{<:STR}) = convert(T, arr)

_isassigned(arr, i...) = isassigned(arr, i...)
_isassigned(arr, i::CartesianIndex) = isassigned(arr, i.I...)
@inline Base.@propagate_inbounds function Base.getindex(a::StringArray{T}, i::Integer...) where T
    offset = a.offsets[i...]
    if offset == UNDEF_OFFSET
        throw(UndefRefError())
    end

    if Missing <: T && offset === MISSING_OFFSET
        return missing
    end

    convert(T, WeakRefString(pointer(a.buffer) + offset, a.lengths[i...]))
end

function Base.similar(a::StringArray, T::Type{<:STR}, dims::Tuple{Vararg{Int64, N}}) where N
    StringArray{T, N}(dims)
end

function Base.empty!(a::StringVector)
    empty!(a.buffer)
    empty!(a.offsets)
    empty!(a.lengths)
end

Base.copy(a::StringArray{T, N}) where {T,N} = StringArray{T, N}(copy(a.buffer), copy(a.offsets), copy(a.lengths))

@inline function Base.setindex!(arr::StringArray, val::WeakRefString, idx::Integer...)
    p = pointer(arr.buffer)
    if val.ptr <= p + sizeof(arr.buffer)-1 && val.ptr >= p
        # this WeakRefString points to data entirely within arr's buffer
        # don't add anything to the buffer in this case.
        # this optimization helps `permute!`
        arr.offsets[idx...] = val.ptr - p
        arr.lengths[idx...] = val.len
    else
        _setindex!(arr, val, idx...)
    end
end

@inline function Base.setindex!(arr::StringArray, val::STR, idx::Integer...)
    _setindex!(arr, val, idx...)
end

@inline function Base.setindex!(arr::StringArray, val::STR, idx::Integer)
    _setindex!(arr, val, idx)
end

function _setindex!(arr::StringArray, val::AbstractString, idx...)
    buffer = arr.buffer
    l = length(arr.buffer)
    resize!(buffer, l + sizeof(val))
    unsafe_copyto!(pointer(buffer, l+1), pointer(val,1), sizeof(val))
    arr.lengths[idx...] = sizeof(val)
    arr.offsets[idx...] = l
    val
end

function _setindex!(arr::StringArray, val::AbstractString, idx)
    buffer = arr.buffer
    l = length(arr.buffer)
    resize!(buffer, l + sizeof(val))
    unsafe_copyto!(pointer(buffer, l+1), pointer(val,1), sizeof(val))
    arr.lengths[idx] = sizeof(val)
    arr.offsets[idx] = l
    val
end

function _setindex!(arr::StringArray, val::Missing, idx)
    arr.lengths[idx] = 0
    arr.offsets[idx] = MISSING_OFFSET
end

function _setindex!(arr::StringArray, val::Missing, idx...)
    arr.lengths[idx...] = 0
    arr.offsets[idx...] = MISSING_OFFSET
end

function Base.resize!(arr::StringVector, len)
    l = length(arr)
    resize!(arr.offsets, len)
    resize!(arr.lengths, len)
    if l < len
        arr.offsets[l+1:len] .= UNDEF_OFFSET # undef
        arr.lengths[l+1:len] .= 0
    end
    arr
end

function Base.push!(arr::StringVector, val::AbstractString)
    l = length(arr.buffer)
    resize!(arr.buffer, l + sizeof(val))
    unsafe_copyto!(pointer(arr.buffer, l + 1), pointer(val,1), sizeof(val))
    push!(arr.offsets, l)
    push!(arr.lengths, sizeof(val))
    arr
end

function Base.push!(arr::StringVector, val::Missing)
    push!(arr.offsets, MISSING_OFFSET)
    push!(arr.lengths, 0)
end

function Base.deleteat!(arr::StringVector, idx::Integer)
    deleteat!(arr.lengths, idx)
    deleteat!(arr.offsets, idx)
    arr
end

function Base.permute!(arr::StringArray{String}, p::AbstractVector)
    permute!(convert(StringArray{WeakRefString{UInt8}}, arr), p)
    arr
end

function Base.vcat(a::StringVector{T}, b::StringVector{T}) where T
    StringVector{T}(vcat(a.buffer, b.buffer), vcat(a.offsets, b.offsets .+ length(a.buffer)), vcat(a.lengths, b.lengths))
end

function Base.append!(a::StringVector{T}, b::StringVector) where T
    append!(a.offsets, b.offsets .+ length(a.buffer))
    append!(a.buffer, b.buffer)
    append!(a.lengths, b.lengths)
    a
end

function Base.append!(a::StringVector{T}, b::AbstractVector) where T
    for x in b
        push!(a, x)
    end
end

end # module


text:
@assert begin
    local_interp=InterpolationFunction(
        :(a0+a1*xi+a2*eta),
        [:xi, :eta],
        [:a0, :a1, :a2],
        [
            [0, 0]
        ]
    )

    abs_interp=InterpolationFunction(
        :(a0+a1*x+a2*y),
        [:x, :y],
        [:a0, :a1, :a2],
        [
            [0, 0]
        ]
    )

    dmn=Domain(
        local_interp,
        [
            (0.0, 1.0),
            (0.0, :(1.0-xi))
        ],
        [
            [0.0, 0.0],
            [1.0, 0.0],
            [0.0, 1.0]
        ]
    )

    fin=Finel(
        [
            [1.0, 1.0],
            [2.0, 1.0],
            [2.0, 2.0]
        ],
        dmn
    )

    resd=Residual(
        [
            (
                :u, abs_interp, [
                    (:du!dx, [1, 0]),
                    (:du!dy, [0, 1])
                ]
            ),
            (
                :v, abs_interp, []
            )
        ],
        :(k*(u*(du!dx-du!dy)-v)),
        3;
        extra_args=[:k]
    )

    surfresd=Residual(
        [
            (
                :u, abs_interp, [
                    (:du!dx, [1, 0]),
                    (:du!dy, [0, 1])
                ]
            ),
            (
                :v, abs_interp, []
            )
        ],
        :([k*du!dx, k*du!dy]),
        3;
        extra_args=[:k]
    )

    resdfunc=get_punctual_resfunc(resd)

    isright=isapprox(1.0, resdfunc(1.0, 2.0, 3.0, -2.0, 1.0))

    varvals=[
        [0.0, 0.0, 1.0],
        [0.0, 1.0, 0.0]
    ]

    pls=vals_to_polys(
        fin,
        resd,
        varvals
    )

    isright=isright && isapprox(
        pls[1], [-1.0, 0.0, 1.0]
    )
    isright=isright && isapprox(
        pls[2], [0.0, 1.0, -1.0]
    )

    vs=polys_to_vars(resd, pls, [1.5, 1.5])

    isright=isright && isapprox(vs, [0.5, 0.0, 1.0, 0.0])

    absc=Vector{Float64}(undef, 2)
    get_abs_coords!(absc, fin, [0.5, 0.5])

    isright=isright && isapprox(
        absc,
        [2.0, 1.5]
    )

    funresd=get_volume_residual_function(resd)

    R=funresd(varvals, fin, -2.0)
    @time R=funresd(varvals, fin, -2.0)
    @time R=funresd(varvals, fin, -2.0)
    @time R=funresd(varvals, fin, -2.0)
    @time R=funresd(varvals, fin, -2.0)
    @time R=funresd(varvals, fin, -2.0)
    @time R=funresd(varvals, fin, -2.0)
    @time R=funresd(varvals, fin, -2.0)

    isright=isright && isapprox(R, 2.0/3)

    surf_funresd=get_surface_residual_function(surfresd)

    R=surf_funresd(varvals, fin, -2.0)
    @time R=surf_funresd(varvals, fin, -2.0)
    @time R=surf_funresd(varvals, fin, -2.0)
    @time R=surf_funresd(varvals, fin, -2.0)
    @time R=surf_funresd(varvals, fin, -2.0)
    @time R=surf_funresd(varvals, fin, -2.0)
    @time R=surf_funresd(varvals, fin, -2.0)
    @time R=surf_funresd(varvals, fin, -2.0)

    isright=isright && isapprox(R, 0.0; atol=1e-10)

    isright
end


text:
@testset "scatter1d assertions" begin
N=10
x=ones(2^N)
Q=[12, 1]
J=[8, 8]
σ0=0.1
# test difference in depth and number of filter banks for FilterBank1d
FilterBanks = get_FilterBanks(N, Q, J, σ0)
Transform = ScatteringTransform1d(3, Q, J)
@test_throws AssertionError scat=Scat.scatter1d(x, Transform, FilterBanks)

# test difference in depth and number of filter banks for FilterBank1dBlock
FilterBankBlocks=get_FilterBanks(N, Q, J, σ0, typ=FilterBank1dBlock)
@test_throws AssertionError scat=Scat.scatter1d(x, Transform, FilterBankBlocks)

# test transform depth
FilterBankBlocks=get_FilterBanks(N, [12,12,1], [8,8,8], σ0, typ=FilterBank1dBlock)
Transform = ScatteringTransform1d(3, Q, J)
@test_throws AssertionError scat=Scat.scatter1d(x, Transform, FilterBankBlocks)
end


text:
"""
The parser separates a raw HTTP Message into its component parts.

If the input data is invalid the Parser throws a `HTTP.ParseError`.

The `parse_*` functions processes a single element of a HTTP Message at a time
and return a `SubString` containing the unused portion of the input.

The Parser does not interpret the Message Headers. It is beyond the scope of the
Parser to deal with repeated header fields, multi-line values, cookies or case
normalization.

The Parser has no knowledge of the high-level `Request` and `Response` structs
defined in `Messages.jl`. However, the `Request` and `Response` structs must
have field names compatible with those expected by the `parse_status_line!` and
`parse_request_line!` functions.
"""
module Parsers

export Header, Headers,
       find_end_of_header, find_end_of_chunk_size, find_end_of_trailer,
       parse_status_line!, parse_request_line!, parse_header_field,
       parse_chunk_size,
       ParseError

include("parseutils.jl")

const emptyss = SubString("",1,0)
const emptyheader = emptyss => emptyss
const Header = Pair{SubString{String},SubString{String}}
const Headers = Vector{Header}

"""
Parser input was invalid.

Fields:
 - `code`, error code
 - `bytes`, the offending input.
"""
struct ParseError <: Exception
    code::Symbol
    bytes::SubString{String}
end

ParseError(code::Symbol, bytes="") =
    ParseError(code, first(split(String(bytes), '\n')))

# Regular expressions for parsing HTTP start-line and header-fields

"""
https://tools.ietf.org/html/rfc7230#section-3.1.1
request-line = method SP request-target SP HTTP-version CRLF
"""
const request_line_regex = r"""^
    (?: \r? \n) ?                       #    ignore leading blank line
    ([!#$%&'*+\-.^_`|~[:alnum:]]+) [ ]+ # 1. method = token (RFC7230 3.2.6)
    ([^.][^ \r\n]*) [ ]+                # 2. target
    HTTP/(\d\.\d)                       # 3. version
    \r? \n                              #    CRLF
"""x

"""
https://tools.ietf.org/html/rfc7230#section-3.1.2
status-line = HTTP-version SP status-code SP reason-phrase CRLF

See:
[#190](https://github.com/JuliaWeb/HTTP.jl/issues/190#issuecomment-363314009)
"""
const status_line_regex = r"""^
    [ ]?                                # Issue #190
    HTTP/(\d\.\d) [ ]+                  # 1. version
    (\d\d\d) .*                         # 2. status
    \r? \n                              #    CRLF
"""x

"""
https://tools.ietf.org/html/rfc7230#section-3.2
header-field = field-name ":" OWS field-value OWS
"""
const header_field_regex = r"""^
    ([!#$%&'*+\-.^_`|~[:alnum:]]+) :    # 1. field-name = token (RFC7230 3.2.6)
    [ \t]*                              #    OWS
    ([^\r\n]*?)                         # 2. field-value
    [ \t]*                              #    OWS
    \r? \n                              #    CRLF
    (?= [^ \t])                         #    no WS on next line
"""x

"""
https://tools.ietf.org/html/rfc7230#section-3.2.4
obs-fold = CRLF 1*( SP / HTAB )
"""
const obs_fold_header_field_regex = r"""^
    ([!#$%&'*+\-.^_`|~[:alnum:]]+) :    # 1. field-name = token (RFC7230 3.2.6)
    [ \t]*                              #    OWS
    ([^\r\n]*                           # 2. field-value
        (?: \r? \n [ \t] [^\r\n]*)*)    #    obs-fold
    [ \t]*                              #    OWS
    \r? \n                              #    CRLF
"""x

const empty_header_field_regex = r"^ \r? \n"x

# HTTP start-line and header-field parsing

"""
Arbitrary limit to protect against denial of service attacks.
"""
const header_size_limit = Int(0x10000)

"""
    find_end_of_header(bytes) -> length or 0

Find length of header delimited by `\\r\\n\\r\\n` or `\\n\\n`.
"""
function find_end_of_header(bytes::AbstractVector{UInt8}; allow_obs_fold=true)

    buf = 0xFFFFFFFF
    l = min(length(bytes), header_size_limit)
    i = 1
    while i <= l
        @inbounds x = bytes[i]
        if x == 0x0D || x == 0x0A
            buf = (buf << 8) | UInt32(x)
            # "Although the line terminator for the start-line and header
            #  fields is the sequence CRLF, a recipient MAY recognize a single
            #  LF as a line terminator"
            # [RFC7230 3.5](https://tools.ietf.org/html/rfc7230#section-3.5)
            buf16 = buf & 0xFFFF
            if buf == 0x0D0A0D0A || buf16 == 0x0A0A
                return i
            end
            # "A server that receives an obs-fold ... MUST either reject the
            # message by sending a 400 (Bad Request) ... or replace each
            # received obs-fold with one or more SP octets..."
            # [RFC7230 3.2.4](https://tools.ietf.org/html/rfc7230#section-3.2.4)
            if !allow_obs_fold && (buf16 == 0x0A20 || buf16 == 0x0A09)
                throw(ParseError(:HEADER_CONTAINS_OBS_FOLD, bytes))
            end
        else
            buf = 0xFFFFFFFF
        end
        i += 1
    end
    if i > header_size_limit
        throw(ParseError(:HEADER_SIZE_EXCEEDS_LIMIT))
    end

    return 0
end

"""
Parse HTTP request-line `bytes` and set the
`method`, `target` and `version` fields of `request`.
Return a `SubString` containing the header-field lines.
"""
function parse_request_line!(bytes::AbstractString, request)::SubString{String}
    re = request_line_regex
    if !exec(re, bytes)
        throw(ParseError(:INVALID_REQUEST_LINE, bytes))
    end
    request.method = group(1, re, bytes)
    request.target = group(2, re, bytes)
    request.version = VersionNumber(group(3, re, bytes))
    return nextbytes(re, bytes)
end

"""
Parse HTTP response-line `bytes` and set the
`status` and `version` fields of `response`.
Return a `SubString` containing the header-field lines.
"""
function parse_status_line!(bytes::AbstractString, response)::SubString{String}
    re = status_line_regex
    if !exec(re, bytes)
        throw(ParseError(:INVALID_STATUS_LINE, bytes))
    end
    response.version = VersionNumber(group(1, re, bytes))
    response.status = parse(Int, group(2, re, bytes))
    return nextbytes(re, bytes)
end

"""
Parse HTTP header-field.
Return `Pair(field-name => field-value)` and
a `SubString` containing the remaining header-field lines.
"""
function parse_header_field(bytes::SubString{String})::Tuple{Header,SubString{String}}
    # First look for: field-name ":" field-value
    re = header_field_regex
    if exec(re, bytes)
        return (group(1, re, bytes) => group(2, re, bytes)),
                nextbytes(re, bytes)
    end

    # Then check for empty termination line:
    re = empty_header_field_regex
    if exec(re, bytes)
        return emptyheader, nextbytes(re, bytes)
    end

    # Finally look for obsolete line folding format:
    re = obs_fold_header_field_regex
    if exec(re, bytes)
        unfold = SubString(replace(group(2, re, bytes), r"\r?\n"=>""))
        return (group(1, re, bytes) => unfold), nextbytes(re, bytes)
    end

    throw(ParseError(:INVALID_HEADER_FIELD, bytes))
end

# HTTP Chunked Transfer Coding

"""
Arbitrary limit to protect against denial of service attacks.
"""
const chunk_size_line_max = 64

const chunk_size_line_min = ncodeunits("0\r\n")

@inline function skip_crlf(bytes, i=1)
    if @inbounds bytes[i] == UInt('\r')
        i += 1
    end
    if @inbounds bytes[i] == UInt('\n')
        i += 1
    end
    return i
end

"""
Find `\\n` after chunk size in `bytes`.
"""
function find_end_of_chunk_size(bytes::AbstractVector{UInt8})
    l = length(bytes)
    if l < chunk_size_line_min
        return 0
    end
    if l > chunk_size_line_max
        l = chunk_size_line_max
    end
    i = skip_crlf(bytes)
    while i <= l
        if @inbounds bytes[i] == UInt('\n')
            return i
        end
        i += 1
    end
    return 0
end

"""
    find_end_of_trailer(bytes) -> length or 0

Find length of trailer delimited by `\\r\\n\\r\\n` (or starting with `\\r\\n`).
[RFC7230 4.1](https://tools.ietf.org/html/rfc7230#section-4.1)
"""
find_end_of_trailer(bytes::AbstractVector{UInt8}) =
    length(bytes) < 2 ? 0 :
    bytes[2] == UInt8('\n') ? 2 :
    find_end_of_header(bytes)

"""
Arbitrary limit to protect against denial of service attacks.
"""
const chunk_size_limit = typemax(Int32)

"""
Parse HTTP chunk-size.
Return number of bytes of chunk-data.

    chunk-size = 1*HEXDIG
[RFC7230 4.1](https://tools.ietf.org/html/rfc7230#section-4.1)
"""
function parse_chunk_size(bytes::AbstractVector{UInt8})::Int

    chunk_size = Int64(0)
    i = skip_crlf(bytes)

    while true
        x = Int64(unhex[@inbounds bytes[i]])
        if x == -1
            break
        end
        chunk_size = chunk_size * Int64(16) + x
        if chunk_size > chunk_size_limit
            throw(ParseError(:CHUNK_SIZE_EXCEEDS_LIMIT, bytes))
        end
        i += 1
    end
    if i > 1
        return Int(chunk_size)
    end

    throw(ParseError(:INVALID_CHUNK_SIZE, bytes))
end

const unhex = Int8[
        -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
    , 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,-1,-1,-1,-1,-1,-1
    ,-1,10,11,12,13,14,15,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,10,11,12,13,14,15,-1,-1,-1,-1,-1,-1,-1,-1,-1
    ,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
]

function __init__()
    # FIXME Consider turing off `PCRE.UTF` in `Regex.compile_options`
    # https://github.com/JuliaLang/julia/pull/26731#issuecomment-380676770
    Base.compile(status_line_regex)
    Base.compile(request_line_regex)
    Base.compile(header_field_regex)
    Base.compile(empty_header_field_regex)
    Base.compile(obs_fold_header_field_regex)
end

end # module Parsers


text:
module TestBiCGStabl

using IterativeSolvers
using Test
using Random
using LinearAlgebra

@testset ("BiCGStab(l)") begin

Random.seed!(1234321)
n = 20

@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A = rand(T, n, n) + 15I
    x = ones(T, n)
    b = A * x

    reltol = √eps(real(T))

    @testset "BiCGStab($l)" for l = (2, 4)
        # Solve without preconditioner
        x1, his1 = bicgstabl(A, b, l, max_mv_products = 100, log = true, reltol = reltol)
        @test isa(his1, ConvergenceHistory)
        @test norm(A * x1 - b) / norm(b) ≤ reltol

        # With an initial guess
        x_guess = rand(T, n)
        x2, his2 = bicgstabl!(x_guess, A, b, l, max_mv_products = 100, log = true, reltol = reltol)
        @test isa(his2, ConvergenceHistory)
        @test x2 == x_guess
        @test norm(A * x2 - b) / norm(b) ≤ reltol

        # The following tests fails CI on Windows and Ubuntu due to a
        # `SingularException(4)`
        if T == Float32 && (Sys.iswindows() || Sys.islinux())
            continue
        end
        # Do an exact LU decomp of a nearby matrix
        F = lu(A + rand(T, n, n))
        x3, his3 = bicgstabl(A, b, Pl = F, l, max_mv_products = 100, log = true, reltol = reltol)
        @test norm(A * x3 - b) / norm(b) ≤ reltol
    end
end

@testset "Termination criterion" begin
    for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = T[ 2 -1  0
              -1  2 -1
               0 -1  2]
        n = size(A, 2)
        b = ones(T, n)
        x0 = A \ b
        perturbation = 10 * sqrt(eps(real(T))) * T[(-1)^i for i in 1:n]

        # If the initial residual is small and a small relative tolerance is used,
        # many iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = bicgstabl!(x, A, b, log=true)
        @test 1 ≤ niters(ch) ≤ n÷2 # BiCGStab(2) makes more than one RHS evaluations per iteration

        # If the initial residual is small and a large absolute tolerance is used,
        # no iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = bicgstabl!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
        @test niters(ch) == 0
    end
end
end

end # module


text:
module Conversions

using ..WhiskerStructs

function hydro2topo(hydro::HydroWhisker)::TopoWhisker
    c₁ = 0.5 * hydro.T * (hydro.γ + hydro.Ac)
    c2 = 0.5 * hydro.T * (hydro.γ - hydro.Ac);
    M = hydro.T * hydro.λ * 0.5;
    ΔL₁ = 0.5 * hydro.ϵ * M;
    ΔL2 = M * (hydro.Phi - 1) + ΔL₁;
    a = sqrt(sqr(c₁) + sqr(ΔL₁));
    k = sqrt(sqr(c2) + sqr(ΔL2));
    α = asin(ΔL₁ / a);
    β = asin(ΔL2 / k);
    b = 0.5 * hydro.T * (1 - hydro.At);
    l = 0.5 * hydro.T * (1 + hydro.At);

    return TopoWhisker(a, b, k, l, α, β, M)
end

function topo2hyrdro(topo::TopoWhisker)::HydroWhisker

end

export hydro2topo

end


text:
function randclasses(classconfig, n, σ=10, q=2, seed=1618)
    Random.seed!(seed)

    nc = length(classconfig)

    configs = Vector{Float64}[]

    for i = 1:nc
        delta = rand(Gamma(q*classconfig[i].p, 1/q), classconfig[i].m)
        config = classconfig[i].start .+ cumsum(delta)
        push!(configs, config)
    end

    ln = length.(configs) .- 1
    class = rand(1:nc, n)
    left = [rand(1:ln[c]) for c in class]
    right = left .+ 1
    left_obs = [rand(Normal(configs[c][l], σ)) for (c, l) in zip(class, left)]
    right_obs = [rand(Normal(configs[c][l], σ)) for (c, l) in zip(class, right)]

    x = [left_obs' ; right_obs']

    obs   =  DataFrame(class=class, left=left, right=right, left_obs=left_obs, right_obs=right_obs)

    (configs=configs, n=n, obs=obs, classconfig=classconfig, x=x)

end

function simpletwoclass(n=1000, σ=10)
    classconfig = [(start=-900, p=180, m=6), (start=50, p=180, m=6)]
    randclasses(classconfig, n, σ)
end

###

# OK, here we go. Physical model for nimble, simplest version.
# 4 parameters: Position left side barrier, position right side barrier, a, b (you'll see), we can fix a subset for now
# Drawing data:

# Be able to draw from exponential with mean 1/b and add a
# 	From right side of barrier, add these random numbers until you reach the end of the domain
# 	From the left side of barrier, subtract these random numbers until you reach the end of the domain
# 	You now got positions of nucleosomes in a virtual cell
# 	Take difference of nearest neighbours to generate one virtual read.

# function rand_cellconfiguration(posleft=-50, posright=50, a=147, b=1/33, xrange=[-1000, 1000])

#     rightpositions = Float64[]
#     leftpositions  = Float64[]

#     cp = posright
#     while cp <= last(xrange)
#         cp = cp + rand(Exponential(1/b)) + a
#         push!(rightpositions, cp)
#     end

#     cp = posleft
#     while cp >= first(xrange)
#         cp = cp - rand(Exponential(1/b)) + a
#         push!(leftpositions, cp)
#     end

#     [leftpositions ; rightpositions]

# end

using Distributions


text:
procs = addprocs(1)

push!(LOAD_PATH,joinpath("..","src"))
push!(LOAD_PATH,joinpath("..","..","..","POMCP.jl","src"))
using Multilane
using POMCP
using POMDPToolbox
using GenerativeModels
using MCTS
using JLD

function POMCP.extract_belief(::POMDPToolbox.FastPreviousObservationUpdater{MLObs}, node::RootNode)
  rand(MersenneTwister(1),node.B)
end

POMCP.initialize_belief(u::FastPreviousObservationUpdater{MLObs}, o::Union{MLState,MLObs}) = o

POMCP.create_belief(u::FastPreviousObservationUpdater{MLObs}) = nothing

POMCP.extract_belief(::POMDPToolbox.FastPreviousObservationUpdater{MLObs}, node::BeliefNode) = node.label[2]

function write!(var_name::AbstractString, val::NoCrashStats, fname::AbstractString="results.jld")
  # XXX kinda hacky but AFAIK JLD has no append functionality
  # NOTE im not sure how to use this functionality AND do the addrequire stuff
  if !isfile(fname)
    d = Dict{ByteString,Any}()
  else
    d = load(fname)
  end
  if var_name in keys(d)
    warn("trying to overwrite value, adding spaghetti")
    var_name = string(var_name,"99999")
  end
  d[var_name] = val
  # TODO addrequire
  save(fname, d)
end

function test_solver(problems::Array{NoCrashMDP,1}, solver, solver_name::AbstractString, initial_states::Array{MLState,1})
  if typeof(solver) <: DPWSolver
    _ec = solver.exploration_constant
    best_score = -Inf
    for (i,ec) in enumerate(Float64[0.1,0.3,1.,3.,10.,30.,100.])
      solver.exploration_constant = ec
      _r = Multilane.evaluate_performance(problems, initial_states, solver)
      score = mean(_r)[3]
      if score > best_score
        r = _r
        best_score = score
      end
    end
    solver.exploration_constant = _ec
  else
    r = Multilane.evaluate_performance(problems, initial_states, solver)
  end
  # TODO hyperparameter optimization here, save best
  write!(solver_name, r)
  return r
end

function test_solver(problems::Array{NoCrashPOMDP,1}, solver, solver_name::AbstractString, initial_states::Array{MLState,1}, bu)
  r = Multilane.evaluate_performance(problems, initial_states, solver, bu)
  write!(solver_name, r)
  return r
end

function test_pomdp(problem::NoCrashPOMDP, solver, solver_name::AbstractString, initial_states, bu)
  r = NoCrashStats([test_run(problem, bu, s, solver) for s in initial_states], problem.rmodel)
  write!(solver_name, r)
  return r
end

nb_lanes = 4
pp = PhysicalParam(nb_lanes,lane_length=100.) #2.=>col_length=8\n",
_discount = 1.
nb_cars=10
dmodel = NoCrashIDMMOBILModel(nb_cars, pp)

N = 10

models = JLD.load("rmodels.jld")
rmodels = models["rmodels"]

# NOTE: random solver and heuristic solver are rmodel agnostic--won't get pareto curve
rmodel = rmodels[1]
mdp = NoCrashMDP(dmodel, rmodel, _discount);
isrng = MersenneTwister(123)
# Initial state is indifferent to rmodel
initial_states = MLState[initial_state(mdp, isrng) for i in 1:N]
problems = NoCrashMDP[mdp for i in 1:N]

#r_random = test_solver(problems, RandomSolver(), "random", initial_states)

#println("Random Average Reward: $(mean(r_random))")

#r_heur = test_solver(problems, SimpleSolver(), "simple_heuristic", initial_states)

#println("Heuristic Average Reward: $(mean(r_heur))")

# exploration constants for POMCP
ecs = [3., 5., 20., 50., 150., 300., 600., 1500., 3000., 5000.]

for (k,(rmodel,ec)) in enumerate(zip(rmodels,ecs))
  mdp = NoCrashMDP(dmodel, rmodel, _discount);

  problems = NoCrashMDP[mdp for i in 1:N]

  dpws = DPWSolver(depth=10,n_iterations=100) #k_action, alpha_action, k_state, alpha_state

  #r_dpw = test_solver(problems, dpws, string("dpw", k), initial_states)

  #println("DPW Average Reward: $(mean(r_dpw))")

  single_dpws = SingleBehaviorSolver(dpws, IDMMOBILBehavior("normal", 30.0, 10.0, 1))

  #r_single_dpw = test_solver(problems, single_dpws, string("single_dpw", k), initial_states)

  #println("Single Behavior DPW Average Reward: $(mean(r_single_dpw))")

  # POMDP stuff

  pomdp = NoCrashPOMDP(dmodel, rmodel, _discount)

  problems = NoCrashPOMDP[pomdp for i in 1:N]

  bu = ParticleUpdater(100, pomdp, MersenneTwister(555))

  pomcp = POMCPDPWSolver( k_observation=0.3,
                          alpha_observation=0.75,
                          k_action=0.5,
                          alpha_action=5.,
                          tree_queries=100,
                          c=ec,
                          eps=0.01)
                          #rollout_solver=SimpleSolver())

                          # Can't use simple solver! Issue is that somehow POMCPDPW is somehow generating an invalid action...

  #pomcp = POMCPSolver()
  r_pomcp = test_pomdp(pomdp, pomcp, string("pomcp", k), initial_states, bu)

  println("POMCP Avg Reward: $(mean(r_pomcp))")

  #r_despot = test_solver(problems, DespotSolver(pomdp,belief), string("dpw", k), initial_states)

  #println("DESPOT Avg Reward: $(mean(r_despot))")

end

rmprocs(procs)

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
struct TanYam7Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  k::rateType
  tab::TabType
end

u_cache(c::TanYam7Cache) = (c.utilde,c.atmp)
du_cache(c::TanYam7Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k)

function alg_cache(alg::TanYam7,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = TanYam7ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype); k2 = zero(rate_prototype) ; k3 = zero(rate_prototype); k4 = zero(rate_prototype)
  k5 = zero(rate_prototype); k6 = zero(rate_prototype) ; k7 = zero(rate_prototype); k8 = zero(rate_prototype)
  k9 = zero(rate_prototype); k10= zero(rate_prototype) ;
  utilde = similar(u,axes(u)); tmp = similar(u); atmp = similar(u,uEltypeNoUnits,axes(u)); k = zero(rate_prototype)
  TanYam7Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,utilde,tmp,atmp,k,tab)
end

alg_cache(alg::TanYam7,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = TanYam7ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct DP8Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  kupdate::rateType
  udiff::rateType
  bspl::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  dense_tmp5::rateType
  dense_tmp6::rateType
  dense_tmp7::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::DP8Cache) = (c.utilde,c.atmp)
du_cache(c::DP8Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13,c.k14,c.k15,c.k16,c.kupdate,c.udiff,c.bspl,c.dense_tmp3,c.dense_tmp4,c.dense_tmp5,c.dense_tmp6,c.dense_tmp7)

function alg_cache(alg::DP8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k1 = zero(rate_prototype); k2  = zero(rate_prototype); k3  = zero(rate_prototype);  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype); k6  = zero(rate_prototype); k7  = zero(rate_prototype);  k8 = zero(rate_prototype)
  k9 = zero(rate_prototype); k10 = zero(rate_prototype); k11 = zero(rate_prototype); k12 = zero(rate_prototype)
  kupdate = zero(rate_prototype); utilde = similar(u,axes(u));
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,axes(u))
  k13 = zero(rate_prototype)
  k14 = zero(rate_prototype)
  k15 = zero(rate_prototype)
  k16 = zero(rate_prototype)
  udiff = zero(rate_prototype)
  bspl = zero(rate_prototype)
  # dense_tmp1 = udiff
  # dense_tmp2 = bspl
  dense_tmp3 = zero(rate_prototype)
  dense_tmp4 = zero(rate_prototype)
  dense_tmp5 = zero(rate_prototype)
  dense_tmp6 = zero(rate_prototype)
  dense_tmp7 = zero(rate_prototype)
  tab = DP8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  DP8Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,kupdate,
           udiff,bspl,dense_tmp3,dense_tmp4,dense_tmp5,dense_tmp6,dense_tmp7,
           utilde,tmp,atmp,tab)
end

alg_cache(alg::DP8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = DP8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct TsitPap8Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  k::rateType
  tab::TabType
end

u_cache(c::TsitPap8Cache) = (c.utilde,c.atmp)
du_cache(c::TsitPap8Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13,c.k)

function alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = TsitPap8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype); k4 = zero(rate_prototype)
  k5 = zero(rate_prototype); k6 = zero(rate_prototype); k7 = zero(rate_prototype); k8 = zero(rate_prototype)
  k9 = zero(rate_prototype); k10 = zero(rate_prototype); k11 = zero(rate_prototype); k12 = zero(rate_prototype)
  k13 = zero(rate_prototype); utilde = similar(u,axes(u)); k = zero(rate_prototype)
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,axes(u));
  TsitPap8Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp,k,tab)
end

alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = TsitPap8ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))


text:
"""
Authors: Yaroslav Salii, 2020+
         Kara Ignatenko, 2021

Submodule responsible for producing the commute and air travel matrices.

travel.jl
2021-08-13 v.0.1: First modular version
2021-09-11 v.0.1.1 restore `talkDnsy` auxiliary function
"""
module Travel
using DataFrames
using FromFile
using Statistics,LinearAlgebra #for talkDnsy()

@from "airports.jl" using Airports: FlightMx
@from "aggregation.jl" using Aggregation: revexplPart

export mkCmtMx,
       mkPsgMx,
       talkDnsy

#--------COMMUTER---FLOW---MATRIX--------------#
"""
Calculate the commuter flow in both directions for all node pairs in `ns`,
given the commute table `wfs`. Returns a matrix `M` such that
`M[i1][i2]` is the total commuter flow from `ns[i1]` to `ns[i2].`
"""
function mkCmtMx(ns::DataFrame,wfs::DataFrame)
    if "Name" ∉ names(ns) || !(["ORG","DST","CMT"] ⊆ names(wfs))
        error("Wrong DF! Terminating.")
    end
    dim = nrow(ns)
    outM = fill(0.0,(dim,dim)) #commute matrix, node-to-node, ordered as `ns`
    ixs = zip(ns.Name, 1:dim) |> Dict #get index by name
    xis = zip(1:dim, ns.Name) |> Dict #get name by index
    for trip ∈ eachrow(wfs) #doggedly fill each trip into the matrix
        outM[ixs[trip.ORG],ixs[trip.DST]]=trip.CMT
    end #reflexive commute is removed beforehand

    return outM
end

"""
NB! `ns` MUST have [:Name]; `pns` MUST have [:Name]; `prt` :Name => [ns_row_indices]
Partition is reversed by a helper function, hence the need for `ns` and `pns` to have :Name.

Calculate the commuter flow in both directions for all node pairs in `ns`,
given the commute table `wfs`. Returns a matrix `M` such that
`M[i1][i2]` is the total commuter flow from `pns[i1]` to `pns[i2]`,
obtained by summing the flows from and to each node in the corresponding
partitions of `prt`.
"""
function mkCmtMx(ns::DataFrame,pns::DataFrame,prt::Dict,wfs::DataFrame)
    if "Name" ∉ names(ns) || !(["ORG","DST","CMT"] ⊆ names(wfs))
        error("Wrong DF! Terminating.")
    end
    dim = nrow(pns)
    ixs = zip(pns.Name, 1:dim) |> Dict # Name => index
   # xis = zip(1:dim, pns.Name) |> Dict # index => Name
    revp = revexplPart(prt,ns,ixs) #generate the reverse-exploded partition
    outM = fill(0.0,(dim,dim)) # part-to-part commute matrix, init to 0.0
    for trip ∈ eachrow(wfs)
        if revp[trip.ORG] ≠ revp[trip.DST] # assuming the commute was not reflexive,
            outM[revp[trip.ORG],revp[trip.DST]]+= trip.CMT #count it in outM[pfrom,pto]
        end
    end
    return outM
end

#------AIR---PASSENGER---FLOW---MATRIX----------------------#

"""
NB! `ns`:[:ID,:Pop,:IATA_Code,:shr]
NB! if provided, `pns` MUST have [:Name]; `prt` :Name => [ns_row_indices]

   mkPsgMx(ns, fmx[, pns, prt][, force_recompute=true|false])

Calculate the air passenger flow in both directions for all node pairs in `ns`,
given the flight matrix `fmx`.
- If `pns` and `prt` are not provided, returns a matrix `M` such that
  `M[i1][i2]` is the total air travel flow from `ns[i1]` to `ns[i2]`
- If `pns` and `prt` are provided, returns a matrix `M` such that
  `M[i1][i2]` is the total air travel flow from `pns[i1]` to `pns[i2]`,
  obtained by summing the flows from and to each node in the corresponding
  partitions of `prt`.
"""
function mkPsgMx(ns::DataFrame,
                 fmx::FlightMx,
                 pns::Union{DataFrame,Nothing}=nothing,
                 prt::Union{Dict, Nothing}=nothing;
                 force_recompute::Bool=false)
    #first determine if we're dealing with aggregated data
    if pns !== nothing && prt !== nothing
        aggregated = true
        #assuming force_recompute is disabled, try to detect by-AP aggregation..
        if !force_recompute && "IATA_Code" ∈ names(pns) && pns.IATA_Code == fmx.apCodes
            #...and if it is, just return the by-AP travel matrix
            println("By-AP aggregation engaged. Using AP-to-AP travel matrix directly.")
            return fmx.M
        end
        #the number of partitions is the size of our output matrix
        dim = length(prt)
        #build a dict mapping the name of a partition to its index in pns
        partIndexByName = Dict{String, Int}(
            partName => findfirst(isequal(partName), pns.Name)[1] for partName in keys(prt)
        )
        #build an array mapping each node in ns to the index of its partition in pns
        partIndexByNode = zeros(Int, nrow(ns))
        #prt contains the ns indices of each node in each partition
        #so we loop through prt assigning the corresponding part. index to each node index
        for partName ∈ keys(prt), nodeIndex ∈ prt[partName]
            partIndexByNode[nodeIndex] = partIndexByName[partName]
        end
    else
        aggregated = false
        #if we're dealing with non-aggregated tract data,
        #this is treated as a trivial case where each node is its own partition
        dim = nrow(ns)
        partIndexByNode = 1:dim
    end
    #attach partition indices to ns so that we remember where to save results
    insertcols!(ns, :partindex => partIndexByNode)
    #group the nodes into subDFs by their designated airport
    groupedbyAP = groupby(ns, :IATA_Code)
    #get the list of airports in the order they appear in `ns`
    retAPs = ns.IATA_Code |> unique
    #initialize output matrix
    outM = fill(0.0, (dim,dim))
    computePsgFlows!(outM, retAPs, fmx, groupedbyAP, accumulate=aggregated)
    select!(ns, Not(:partindex)) #drop the partindex column to keep ns clean
    return outM
end

"""
For each combination (i.e. unordered pair with distinct entries)
of groups in `groupedByAP`, and for each node pair within these combinations,
calculate the passenger flows for that pair and accumulate them in `outM`
at a location determined by the `partindex` column.
So M[i1, i2] is the total passenger flow from all nodes that have the `partindex` i1
to all nodes that have the `partindex` i2.

If `accumulate` is `false`, values in `outM` will be overwritten rather than accumulated,
which is useful for saving on memory accesses in the trivial case where each partition
has only a single node.
"""
@inline function computePsgFlows!(outM::Matrix{Float64},
                                  retAPs::AbstractArray{String},
                                  fmx::FlightMx,
                                  groupedByAP::GroupedDataFrame;
                                  accumulate::Bool=true)
    numAPs = length(retAPs)
    #for each combination of airports (AP1, AP2)...
    #(this weird iteration is to avoid repeating pairs)
    for i1 ∈ 1:numAPs, i2 ∈ (i1 + 1):numAPs
        # get IATA codes of AP1 and AP2
        iata1 = retAPs[i1]
        iata2 = retAPs[i2]
        # get AP1->AP2 and AP2->AP1 passenger numbers
        flow1to2 = fmx.M[fmx.ix[iata1], fmx.ix[iata2]]
        flow2to1 = fmx.M[fmx.ix[iata2], fmx.ix[iata1]]
        # skip over airport pairs that have no travel between them
        if flow1to2 == 0.0 && flow2to1 == 0.0
            continue
        end
        grp1 = groupedByAP[i1]
        grp2 = groupedByAP[i2]
        #pre-collect relevant columns to optimize performance
        indices1 = collect(grp1.partindex) :: Vector{Int}
        indices2 = collect(grp2.partindex) :: Vector{Int}
        shr1     = collect(grp1.shr)       :: Vector{Float64}
        shr2     = collect(grp2.shr)       :: Vector{Float64}
        #for each pair of nodes (n1 ∈ AP1, n2 ∈ AP2)...
        for n1 ∈ 1:nrow(grp1), n2 ∈ 1:nrow(grp2)
            #calculate psg flow n1->n2 and n2->n1 and store the result
            outi1 = indices1[n1]
            outi2 = indices2[n2]
            #do not count travel within the same subdivision
            if outi1 == outi2
                continue
            end
            result1to2 = psg(shr1[n1], shr2[n2], flow1to2)
            result2to1 = psg(shr2[n2], shr1[n1], flow2to1)
            #if we're accumulating, add the results to corresponding outM location...
            if accumulate
                outM[outi1, outi2] += result1to2
                outM[outi2, outi1] += result2to1
            #..otherwise overwrite
            else
                outM[outi1, outi2] = result1to2
                outM[outi2, outi1] = result2to1
            end
        end
    end
end

"""
Approximate passenger flow between two nodes, given both passenger shares
(`fromshr` and `toshr`) and the passenger flow between the corresponding
designated airports (`totalflow`)
"""
function psg(fromshr, toshr, totalflow)
    fromshr * toshr * totalflow
end

"""
Display (`println`) sparsity-related statistics for commute matrices.
Computes the input's `density`, the fraction of nonzero elements assuming
the diagonal may only have zeros, and also the minimum, median, and maximum
**nonzero** element of the matrix.
"""
function talkDnsy(M)
    dim = size(M,1)
    if ! (filter(x -> x > 0.0, diag(M)) |> isempty)
        println(stderr,"Warning: nonzero diagonal elements found.")
    end
    #nonzero entries
    nzs = filter(x -> x > 0.0, M); a= nzs |> length;
    dnsy = a/ dim^2; b = dim*(dim-1) #b is max no. nonzero elements (all 'cept diag)
    println("We've got ",a," ≠0 connections")
    println("With at most ",b, " ≠0 entries, we get the density ",dnsy)
    print("min ", nzs |> minimum, "   median ", nzs |> median, "   max ", nzs |> maximum)
    return(nnonzero = a, maxnonzero= b, density = dnsy, dim = size(M,1))
end

end #end module Oboe.Travel


text:
using DataFrames
using LinearAlgebra
using MixedModels
using MixedModelsSim
using PooledArrays
using Tables
using Test

@testset "factorproduct" begin
	items = (item = nlevels(50, 'I'), category = repeat(["ingroup", "outgroup"], inner=25))
	df = columntable(factorproduct(items, (subj = nlevels(100),)))
	@test length(df) == 3
	@test length(first(df)) == 5000
end

@testset "levels" begin
    levs = nlevels(100, 'k')
    @test length(levs) == 100
    @test length(unique(levs)) == 100
    @test first(levs) == "k001"
end

@testset "pooled!" begin
    subject = (subj = ["S1","S2","S3","S4","S5"], age=["O","O","O","Y","Y"]);

    @test_throws ArgumentError pooled!(subject)

    df = pooled!(DataFrame(subject))

    @test df.subj isa PooledArray
    @test df.age isa PooledArray
end

@testset "flatlowertri" begin
    l = LowerTriangular([1 0 0; 2 3 0; 4 5 6])
    @test all(flatlowertri(l) .== [1, 2, 4, 3, 5, 6])
end

@testset "update!" begin
    @testset "LMM" begin
        fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)),
                  MixedModels.dataset(:sleepstudy))
        update!(fm1; θ=[1.0, 0.0, 1.0])

        @test all(values(first(fm1.σs)) .== fm1.σ)
        @test only(VarCorr(fm1).σρ.subj.ρ) == 0.0
    end

    @testset "GLMM" begin
        cbpp = MixedModels.dataset(:cbpp)
        gm1 = fit(MixedModel, @formula((incid/hsz) ~ 1 + period + (1|herd)), cbpp, Binomial();
                wts=cbpp.hsz, fast=true)

        β = repeat([-1.], 4)
        θ = [0.5]

        β₀ = copy(fixef(gm1))
        update!(gm1; θ=θ)
        @test all(values(first(gm1.σs)) .== θ)
        @test all(gm1.β .== β₀)

        refit!(gm1, fast=false)
        β₀ = copy(fixef(gm1))
        update!(gm1; θ=θ)
        @test all(values(first(gm1.σs)) .== θ)
        @test all(gm1.β .== β₀)
    end

end

@testset "create_re" begin
    fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)),
              MixedModels.dataset(:sleepstudy))

    σs = values(first(fm1.σs)) ./ fm1.σ
    ρ = only(first(VarCorr(fm1).σρ).ρ)
    corr = [1 ρ; ρ 1]

    subjre = create_re(σs...; corrmat=corr)

    @test all(only(fm1.λ) .≈ subjre)
    @test create_re(σs...) == diagm([σs...])

    fm1unfitted = LinearMixedModel(@formula(reaction ~ 1 + days + (1 + days|subj)),
                                   MixedModels.dataset(:sleepstudy))

    update!(fm1unfitted, subjre)

    vcu, vc = VarCorr(fm1unfitted), VarCorr(fm1)

    @test all(values(vcu.σρ.subj.σ) .≈ values(vc.σρ.subj.σ))
    @test all(values(vcu.σρ.subj.ρ) .≈ values(vc.σρ.subj.ρ))
end