text:
module Filthy

using StaticArrays
using GrowableArrays
using Parameters
using ForwardDiff
using MathProgBase
using LinearAlgebra

include("types.jl")
include("filter.jl")
include("mpb.jl")
include("methods.jl")

function SSMOptim(m, Z, CH, T, R, CQ, a, P, a1!, P1!, P1inf, Pstar, Y)
    t = prototypical(m)
    CS = ForwardDiff.pickchunksize(length(t))
    p = SSMOptimPar(Z, CH, T, R, CQ, a, P, a1!, P1!, P1inf, Pstar)
    jcfg = ForwardDiff.GradientConfig(nothing, t, ForwardDiff.Chunk(CS))
    hcfg = ForwardDiff.HessianConfig(nothing, t, ForwardDiff.Chunk(CS))
    M1 = map(a->Array{ForwardDiff.Dual{Nothing,Float64,CS},2}(size(a)...), (Z, CH, T, R, CQ))
    a1 = Array{ForwardDiff.Dual{Nothing,Float64,CS},1}(length(a))
    M2 = map(a->Array{ForwardDiff.Dual{Void,Float64,CS},2}(size(a)...), (P, P1inf, Pstar))
    d = SSMOptimParDual(M1..., a1, M2..., jcfg, hcfg)
    d = SSMOptim(m, p, d, Y)
    d
end

function get_parameters(s::SSMOptim, ::Type{T}) where T
    (s.p.Z::Array{T, 2},
     s.p.H::Array{T, 2},
     s.p.T::Array{T, 2},
     s.p.R::Array{T, 2},
     s.p.Q::Array{T, 2},
     s.p.a1::Array{T, 1},
     s.p.P1::Array{T, 2},
     s.p.P1inf::Array{T, 2},
     s.p.Pstar::Array{T, 2})
end

function get_parameters(s::SSMOptim, ::Type{T}) where T<:ForwardDiff.Dual
    (s.d.Z::Array{T, 2},
     s.d.H::Array{T, 2},
     s.d.T::Array{T, 2},
     s.d.R::Array{T, 2},
     s.d.Q::Array{T, 2},
     s.d.a1::Array{T, 1},
     s.d.P1::Array{T, 2},
     s.p.P1inf::Array{Float64, 2},
     s.p.Pstar::Array{Float64, 2})
end

function KFFiltered(a::T, P, Pinf::Nothing) where T<:Real
    KFFiltered([a], [P], [P], [P], [1], [0], BitArray([1]))
end

function KFFiltered(a::T, P, Pinf::Nothing) where T<:AbstractArray
    KFFiltered(GrowableArray(a), GrowableArray(P), [], [], [1], [1],  BitArray([1]))
end

function KFFiltered(a1, P1, P1inf::AbstractMatrix)
    ## Check positive element on the diagonal
    ## To Do: Check that P1inf is diagonal with diagonal entry > 0
    hasended = all(diag(P1inf) .== 0) ? 1 : 0
    a1′, P1′, Pstar = setupdiffusematrices(a1, P1, P1inf)
    KFFiltered(GrowableArray(a1′), GrowableArray(P1′), GrowableArray(Pstar), GrowableArray(P1inf), [1], [1],  BitArray([hasended]))
end

function KFFiltered(a1, P1, P1inf::AbstractFloat)
    ## Check positive element on the diagonal
    @assert P1inf > 0 "Initial value to diffuse part should be > 0"
    KFFiltered(GrowableArray(a1), GrowableArray(P1), GrowableArray(0.0), GrowableArray(P1inf), [1], [1],  BitArray([0]))
end

function KFCache(T::AbstractMatrix, p::Int, m::Int)
    KFCache(GrowableArray(zeros(p,p)),
            GrowableArray(zeros(p,p)),
            GrowableArray(zeros(p)),
            Array{Float64}(undef, p),
            Array{Float64}(undef, p),
            Array{Float64}(undef, p, p),  # C
            Array{Float64}(undef, m, p),  # M
            Array{Float64}(undef, p, p),  # ZM
            Array{Float64}(undef, m, p),  # TM
            Array{Float64}(undef, m, p),  # K
            Array{Float64}(undef, m, m),  # KZ
            Array{Float64}(undef, m, m),  # KZp
            Array{Float64}(undef, m, m),  # TPL
            Array{Float64}(undef, m, m),   # L
            Array{Float64}(undef, p),     # a′
            Array{Float64}(undef, m, m)   # P′
            )
end

function KFCache(T::AbstractFloat, p::Int, m::Int)
    @assert p==1 "Something wrong"
    KFCache(GrowableArray(zero(T)),
            GrowableArray(zero(T)),
            GrowableArray(zero(eltype(T))),
            Array{Float64}(undef, 1),
            Array{Float64}(undef, 1),
            Array{Float64}(undef, 1, 1),  # C
            Array{Float64}(undef, 1, 1),  # M
            Array{Float64}(undef, 1, 1),  # ZM
            Array{Float64}(undef, 1, 1),  # TM
            Array{Float64}(undef, 1, 1),  # K
            Array{Float64}(undef, 1, 1),  # KZ
            Array{Float64}(undef, 1, 1),  # KZp
            Array{Float64}(undef, 1, 1),  # TPL
            Array{Float64}(undef, 1, 1),  # L
            Array{Float64}(undef, 1),  # a
            Array{Float64}(undef, 1, 1)  # P
            )
end

function LinearStateSpace(Z::ZT, H::HT, T::TT, R::RT, Q::QT, a1::AZ, P1::PZ; P1inf::Union{Nothing, PZ}=nothing) where {ZT, HT, TT, RT, QT, AZ, PZ}
    @assert isa(P1, TT) "P1 must be of type $(typeof(TT))"
    if isa(Z, AbstractMatrix)
        @assert isa(a1, AbstractVector) "a1 must be of type Vector"
    elseif isa(Z, AbstractFloat)
        @assert isa(a1, AbstractFloat) "a1 must be of type AbstractFloat"
    end
    #checksizes(Z, H, T, R, Q, α0, P1)
    m = numstates(T)::Int64
    p = nummeasur(Z)::Int64
    params = KFParms(Z, H, T, R, Q)
    filter = KFFiltered(a1, P1, P1inf)
    smooth = KFSmoothed(Array{Float64}(undef, 1,1), Array{Float64}(undef, 1,1,1),
                        Array{Float64}(undef, 1,1), Array{Float64}(undef, 1,1,1))
    inival = KFInitVal(a1, P1)
    caches = KFCache(T, p, m)
    LinearStateSpace(params, filter, smooth, inival, caches, [0.0], Ref(1))
end

function smooth!(ss::LinearStateSpace)
    n = ss.t[]::Int64
    p, m, r = size(ss.p)
    @assert n > 1 "There is not data to smooth over"
    @unpack Z, H, T, R, Q = ss.p
    â  = Array{Float64}(m, n)
    V  = Array{Float64}(m, m, n)
    r  = Array{Float64}(p, n)
    N  = Array{Float64}(p, p, n)

    F⁻¹ = ss.c.F⁻¹
    Y  = ss.c.Y
    a  = ss.f.a
    P  = ss.f.P
    d  = ss.f.d[1]
    r[:, n] = 0.0
    N[:, :, n] = 0.0

    @inbounds for t in n:-1:d+1 #(d+1)
        ## Hoisting
        Finv = F⁻¹[t]
        PP = P[t-1]
        v = Y[t] .- Z*a[t-1]
        L = T .- T*P[t-1]*Z'*Finv*Z
        r[:, t-1] = Z'*Finv*v .+ L'*r[:, t]
        â[:, t] = a[t-1] .+ PP*r[:,t-1]
        N[:, :, t-1] = reshape(convert(Matrix, Z'*Finv*Z + L'*N[:,:,t]*L), (p, p, 1))
        V[:, :, t] = reshape(convert(Matrix, P[t-1] .- P[t-1]*N[:, :, t-1]*P[t-1]), (p, p, 1))
    end

    if d>=1
        rdz = Array{Float64}(p, d)
        rdo = Array{Float64}(p, d)
        Ndz = Array{Float64}(p, p, d)
        Ndo = Array{Float64}(p, p, d)
        Ndt = Array{Float64}(p, p, d)
        rdz[:, d] = r[:, d]
        rdo[:, d] = 0.0
        Ndz[:, :, d] = N[:, :, d]
        Ndo[:, :, d] = 0.0
        Ndt[:, :, d] = 0.0
    end

    ## Diffuse step
    for t in d:-1:2
        Pinf = ss.f.Pinf[t-1]
        Pstar = ss.f.Pstar[t]
        v = Y[t] .-Z*a[t-1]
        if ss.f.flagF[t] == 1
            Finf = Z*Pinf*Z'
            Finfinv = inv(Finf)
            Kinf = T*Pinf*Z*Finfinv
            Linf = T - Kinf*Z
            Fstar = Z*Pstar*Z' + H
            Kstar = (T*Pstar*Z' - Kinf*Fstar)*Finfinv
            F♯ = Kstar'Ndz[:,:,t]*Kstar - Finfinv*Fstar*Finfinv
            rdz[:,t-1] = Linf'rdz[:, t]
            rdo[:,t-1] = Z'*(Finfinv*v - Kstar'rdz[:,t]) + Linf'*rdo[:,t]
            Ndz[:,:,t-1] = Linf'*Ndz[:,:,t]*Linf
            A = Linf'Ndz[:,:,t]*Kstar*Z
            Ndo[:,:,t-1] = Z'Finfinv*Z + Linf'*Ndo[:, :, t]*Linf - A - A'
            A = Linf'Ndo[:,:,t]*Kstar*Z
            Ndt[:,:,t-1] = Z'*F♯*Z + Linf'*Ndt[:,:,t]*Linf - A - A'
            â[:,t] = a[t-1] .+ Pstar*rdz[:,t-1] .+ Pinf*rdo[:,t-1]
            A = Pinf*Ndo[:,:,t-1]*Pstar
            B = Pinf*Ndt[:,:,t-1]*Pinf
            V[:,:,t] = Pstar - Pstar*Ndz[:,:,t-1]*Pstar - (A + A') - B
        else
            Fstar = Z'*Pstar*Z + H
            Fstarinv = inv(Fstar)
            Lstar = T - Kstar*Z
            Kstar = T*Pstar*Z'*Fstarinv
            rdz[:,t-1] = Z'Fstarinv*v - Lstar'rdz[:, t]
            rdo[:,t-1] = T'rdo[:,t]
            Ndz[:,:,t-1] = Z'*Fstarinv*Z + Lstar'Ndz[:,:,t]*Lstar
            Ndo[:,:,t-1] = T'Ndo[:,:,t]*Lstar
            Ndt[:,:,t-1] = T*Ndt[:,:,t]*T
            â[:,t-1] = a[t] .+ Pstar*rdz[:,t] .+ Pinf*rdz[:,t-1]
            V[:,:,t] = Pstar - Pstar*Ndz[:,:,t-1]
        end
    end
    ss.s.r = r[:, 1:n-1]'
    ss.s.a = â'[2:end, :]
    ss.s.V = V[:, :, 2:end]
    ss.s.N = N[:, :, 1:n-1]
    nothing
end

function fastloglik(Z::AbstractMatrix, H, T, R, Q, a′, P′, Pinf′, Pstar′, Y::AbstractMatrix)
    n, p = size(Y)
    ll = zero(eltype(Z))
    j = 1; p = 1
    c = -p*n*log(2π)/2
    while (any(Pinf′ .> 0)) & (j < n)
        a′, P′, v, F, F⁻¹, Pinf′, Pstar′, flag = filterstep(a′, P′, Pinf′, Pstar′, Z, H, T, R, Q, Y[j,:])
        ll += loglik(v, F, F⁻¹, flag)
        j += 1
    end
    for t in j:n
        a′, P′, v, F, F⁻¹ = filterstep(a′, P′, Z, H, T, R, Q, Y[t, :])
        ll += loglik(v, F, F⁻¹)
    end
    ll +c
end

function fastloglikscalar(Z::FF, H::FF, T::FF, R::FF, Q::FF, a1::G, P1::G, P1inf::G, Pstar::G, y) where {FF<:Real, G<:AbstractFloat}
    ll = zero(eltype(Z));  n = length(y)
    a = a1; P = P1; j = 1; p = 1; c = -n/2*log(2π)
    while (P1inf > 0) & (j < n)
        a, P, v, F, F⁻¹, P1inf, Pstar, flag = filterstep(a, P, P1inf, Pstar, Z, H, T, R, Q, y[j])
        ll += loglik(v, F, F⁻¹, flag)
        j += 1
    end
    for t in j:n
        a, P, v, F, F⁻¹ = filterstep(a, P, Z, H, T, R, Q, y[t])
        ll += loglik(v, F, F⁻¹)
    end
    ll + c
end

function ssm_obj_fun(s::SSMOptim, theta)
    Z, H, T, R, Q, a1, P1, P1inf, Pstar = get_parameters(s, eltype(theta))
    ## H and Q are the chol
    remask!(s.m, (Z, H, T, R, Q), theta)
    try
        s.p.a1!(s, a1)
        s.p.P1!(s, P1)
        fastloglik(Z, H*H', T, R, Q*Q', a1, P1, P1inf, Pstar, s.Y)
    catch
        -Inf
    end
end

function ssm_obj_fun_scalar(s::OptimSSMScalar, theta)
    @unpack a1, P1, P1inf, Pstar, idx, par, y = s
    x = similar(theta, length(par))
    copy!(x, par)
    x[idx] = theta
    -fastloglikscalar(x[1], exp(x[2]), x[3], x[4], exp(x[5]), a1, P1, P1inf, Pstar, y)
end
"""
fit(Z::A, CH::A, T::A, R::A, CQ::A, a1::AZ, P1::PZ, P1inf::PW, Y, start, lower, upper) where {A<:AbstractMatrix, AZ<:Union{Function, AbstractVector}, PZ<:Union{Function, AbstractMatrix}, PW}

"""
function fit(::Type{LinearStateSpace}, Z::A, CH::A, T::A, R::A, CQ::A, a1!::Function,
            P1!::Function, P1inf::PINF, Y, x0, lower, upper, solver = IpoptSolver()) where {A<:AbstractMatrix, PINF}
    @assert islowertriangular(CH) "CH must be lower triangular"
    @assert islowertriangular(CQ) "CQ must be lower triangular"
    Z₀, CH₀, T₀, R₀, CQ₀ = map(obj-> convert(Matrix, obj), (copy(Z), copy(CH), copy(T), copy(R), copy(CQ)))
    mask = Masked((Z₀, CH₀, T₀, R₀, CQ₀))
    allfixed(mask) && (print_with_color(:cyan, "All State Space parameters are defined."); return nothing)
    ## Setup and do maximum likelihood
    ## Copy matrices to Matrix{Float64}
    a1m = Array{eltype(Z), 1}(size(T,1))
    P1m = Array{eltype(Z), 2}(size(R, 1), size(R, 1))
    a, P, Pstar = setupdiffusematrices(a1m, P1m, P1inf)
    #fitobj = OptimSSM(mask, Z₀, CH₀, T₀, R₀, CQ₀, a, P,  P1inf, Pstar, Y)
    s = SSMOptim(mask, Z₀, CH₀, T₀, R₀, CQ₀, a, P, a1!, P1!, P1inf, Pstar, Y)
    m = MathProgBase.NonlinearModel(solver)
    MathProgBase.loadproblem!(m, length(x0), 0, lower, upper, Array{Float64}(0), Array{Float64}(0),  :Max, s)
    MathProgBase.setwarmstart!(m, x0)
    MathProgBase.optimize!(m)
    sol = MathProgBase.getsolution(m)
    remask!(mask, (Z₀, CH₀, T₀, R₀, CQ₀), sol)
    ZZ = SMatrix{size(Z)...}(Z₀)
    HH = SMatrix{size(CH)...}(CH₀*CH₀')
    TT = SMatrix{size(T)...}(T₀)
    RR = SMatrix{size(R₀)...}(R₀)
    QQ = SMatrix{size(CQ)...}(CQ₀*CQ₀')
    aa = SVector{length(a)}(s.p.a1)
    PP = SMatrix{size(P)...}(s.p.P1)
    PPinf = SMatrix{size(P)...}(s.p.P1inf)
    LinearStateSpace(ZZ, HH, TT, RR, QQ, aa, PP, P1inf=PPinf)
end

function fit!(ss::LinearStateSpace, Y, x0, lower, upper, solver)
end

function fit(::Type{LinearStateSpace}, Z::A, H::A, T::A, R::B, Q::A, a1::AZ, P1::PZ, P1inf::PZ, Y, start) where {A<:AbstractFloat, B, AZ<:AbstractFloat, PZ<:AbstractFloat}
    @assert !((P1inf != zero(P1inf)) && (P1inf != one(P1inf))) "P1inf must be either 1.0 or 0.0"
    Pstar = one(P1inf)
    if P1inf == zero(P1inf)
        @assert P1 > 0 "P1 and P1inf cannot be both zero"
    elseif P1inf == one(P1inf)
        P1 = zero(P1inf)
        Pstar = zero(P1inf)
    end
    b = [Z, H, T, R, Q]
    idx = isnan.(b)
    b[idx] = start
    start[:] = log.(b[idx .& [false, true, false, false, true]])
    s = OptimSSMScalar(b, idx, a1, P1, P1inf, Pstar, y)
    f(x) = futil_scalar(s, x)
    d = Optim.OnceDifferentiable(f, start, autodiff = :forward)
    out = Optim.optimize(d, start, BFGS())
    sol = Optim.minimizer(out)
    b[idx] = sol
    sol[:] = exp.(b[idx .& [false, true, false, false, true]])
    b[idx] = sol
    out.minimizer[:] = sol
    out.initial_x[:] = sol
    LinearStateSpace(b[1], b[2], b[3], b[4], b[5], a1, P1, P1inf=P1inf)
end

function reset!(ss::LinearStateSpace)
    ## Clean filtered and smoothed states
    resize!(ss.f.a.data, 1)
    resize!(ss.f.P.data, 1)
    resize!(ss.f.Pstar.data, 1)
    resize!(ss.f.Pinf.data, 1)
    resize!(ss.f.flagF, 1)
    resize!(ss.c.Y.data, 1)
    resize!(ss.c.F.data, 1)
    resize!(ss.c.F⁻¹.data, 1)
    ss.f.d[1] = 1
    ss.f.hasdiffended = [false]
    ss.t[] = 1
    nothing
end

function simulate(cf::LinearStateSpace, nout)
    @unpack Z, H, T, R, Q = cf.p
    p, m = size(Z)
    Y = Array{Float64, 2}(p, nout+1)
    α = Array{Float64, 2}(m, nout+2)
    fH = chol(H)'
    fQ = chol(Q)'
    fP = chol(cf.i.P1)'
    α[:,1] = T*cf.i.a1 + R*fP*randn(m)
    for j in 1:nout+1
        Y[:, j] = Z*α[:,j] + fH*randn(p)
        α[:, j+1] = T*α[:,j] + R*fQ*randn(m)
    end
    (Y[:,1:end-1], α[:,2:end-1])
end

loglik(cf) = cf.loglik[]::Float64

export LinearStateSpace, simulate, filter!, smooth!, fit, reset!, get_parameters

end # module


text:
using Test
using LinearAlgebra, DomainSets, GridArrays, StaticArrays
using BasisFunctions

@testset begin
    F = Fourier(3) → -1..1
    S = SynthesisOperator(F)
    G = GridSampling(PeriodicEquispacedGrid(3,-1,1))
    @test Matrix(G*S)≈Matrix(G*F)
    # @test_throws MethodError SynthesisOperator(F,nothing)''

    g = interpolation_grid(Fourier(3) → -1..1)
    g1 = gram(Fourier(3) → -1..1,discretemeasure(map_grid(forward_map(g),FourierGrid(6))))
    g2 = gram(Fourier(3),discretemeasure(FourierGrid(6)))
    g3 = BasisFunctions.default_gram(Fourier(3),discretemeasure(FourierGrid(6)); warnslow=false)

    @test g1≈g2≈g3

    g1 = gram(Fourier(3) → -1..1)
    g2 = gram(Fourier(3))
    g3 = BasisFunctions.default_gram(Fourier(3); warnslow=false)

    @test g1≈g2≈g3

    @test weights(discretemeasure(subgrid(PeriodicEquispacedGrid(10,0,1)^2,(0.0..0.5)^2))) isa BasisFunctions.OuterProductArray

end

ext: jl


text:
using Oceananigans.Grids: topology

minimum_grid_spacing(Δx, TX          ) = minimum(Δx)
minimum_grid_spacing(Δx, ::Type{Flat}) = Inf

"Returns the time-scale for advection on a regular grid across a single grid cell."
function cell_advection_timescale(u, v, w, grid::RegularRectilinearGrid)

    umax = maximum(abs, u)
    vmax = maximum(abs, v)
    wmax = maximum(abs, w)

    topo = topology(grid)

    Δx = minimum_grid_spacing(grid.Δx, topo[1])
    Δy = minimum_grid_spacing(grid.Δy, topo[2])
    Δz = minimum_grid_spacing(grid.Δz, topo[3])

    return min(Δx/umax, Δy/vmax, Δz/wmax)
end

function cell_advection_timescale(u, v, w, grid::VerticallyStretchedRectilinearGrid)

    umax = maximum(abs, u)
    vmax = maximum(abs, v)
    wmax = maximum(abs, w)

    topo = topology(grid)

    Δx = minimum_grid_spacing(grid.Δx,    topo[1])
    Δy = minimum_grid_spacing(grid.Δy,    topo[2])
    Δz = minimum_grid_spacing(grid.Δzᵃᵃᶠ, topo[3])

    return min(Δx/umax, Δy/vmax, Δz/wmax)
end

cell_advection_timescale(model) =
    cell_advection_timescale(model.velocities.u.data.parent,
                             model.velocities.v.data.parent,
                             model.velocities.w.data.parent,
                             model.grid)


text:
using VegaLite
using Statistics

function build_plot_dataframe(agg::T) where {T}

    if T == Vector{Aggregate}
        # all fine
    elseif T == Aggregate
        agg = Aggregate[agg]
    else
        @error "Unsupported type - needs to be Array{Aggregate, 1}"
        return nothing
    end

    # Figure out the total number of entries required here:
    N = sum((p -> p.N).(agg))

    # Return empty DF if aggregate is empty
    if N == 0
        return DataFrame()
    end

    # Not sure if this is the best option:
    # create a dict with some keys and then cast it into
    # a DataFrame object. Pre-allocating a df and filling
    # it later on doesn't seem to be in the spirit of the
    # DataFrames package..

    df = DataFrame(
        index = fill(-1, N),
        instrument = fill("none", N),
        mode = fill("unknown", N),
        SZA = fill(-1.0, N),
        VZA = fill(-1.0, N),
        lon = fill(-999.99, N),
        lat = fill(-999.99, N),
        PPFD = fill(-999.99, N),
        )

    last_index = 1
    for i in 1:length(agg)
        df[last_index:last_index + agg[i].N - 1, :index] .= i
        df[last_index:last_index + agg[i].N - 1, :instrument] = get_instrument.(agg[i].scenes)
        df[last_index:last_index + agg[i].N - 1, :mode] = get_mode.(agg[i].scenes)
        df[last_index:last_index + agg[i].N - 1, :SZA] = get_sza.(agg[i].scenes)
        df[last_index:last_index + agg[i].N - 1, :VZA] = get_vza.(agg[i].scenes)
        df[last_index:last_index + agg[i].N - 1, :lon] = get_lon.(agg[i].scenes)
        df[last_index:last_index + agg[i].N - 1, :lat] = get_lat.(agg[i].scenes)
        df[last_index:last_index + agg[i].N - 1, :PPFD] = get_ppfd.(agg[i].scenes)

        last_index = last_index + agg[i].N
    end

    return df
end

function aggregate_overview_vl(agg::Vector{Aggregate})

    # Construct DataFrame that holds most needed into
    df = build_plot_dataframe(agg)

    # From that, let's construct a grouped dataframe with some aggregate (of aggregates)
    # values ..

    dfg = groupby(df, :index)

    # Containing the number of scenes per instrument mode
    mode_symbols = [Symbol(y) for y in unique(df.mode)]
    df_modes = combine(dfg, ([:mode => (x -> sum(x .== y)) => Symbol(y) for y in unique(df.mode)])...)

    # Containing the number of scenes per instrument
    instrument_symbols = [Symbol(y) for y in unique(df.instrument)];
    df_instrument = combine(dfg, ([:instrument => (x -> sum(x .== y)) => Symbol(y) for y in unique(df.instrument)])...);

    # Containing the mean locations of the aggregates
    df_locs = combine(dfg,
                      :lon => mean => :lon,
                      :lat => mean => :lat,
                      :index => length => :count)

    df_plot = innerjoin(df_modes, df_locs, df_instrument, on=:index)

    # Turn number of scenes per instrument into percentage
    for inst in instrument_symbols
        new_percentage = convert.(Int32, round.(100 .* df_plot[:, inst] ./ df_plot[:, "count"]))
        df_plot[:, inst] .= new_percentage
    end

    #################
    # Produce plots #
    #################

    # A set of bar plots showing the occurrence of different instrument
    # modes for the array of aggregates
    p1 = df_plot |> @vlplot(repeat=mode_symbols, columns=3) +
        @vlplot(mark={:line, "tooltip" = "content" => "data"},
                y={field={repeat=:repeat}, bin=false, type=:quantitative},
                x={"index", title="Aggregate #"})
    display(p1)

    # A bar plot showing the number of scenes per aggregate
    p2 = df_plot |> @vlplot() + @vlplot(
        mark={:bar, "tooltip" = "content" => "data"},
        x={"index:o", title="Aggregate #"}, y={"count:q", title="Number of scenes"}
    )
    display(p2)

    # Make a lon/lat scatter plot highlighting which instrument
    # was measuring where
    xrange = [minimum(df_plot.lon), maximum(df_plot.lon)]
    yrange = [minimum(df_plot.lat), maximum(df_plot.lat)]

    for inst in instrument_symbols
        _plot = df_plot |> @vlplot(
            mark={:circle, stroke=:black, strokeWidth=1.0, opacity=1.0, "tooltip" = "content" => "data"},
            title="Instrument fraction in aggregate",
            x={:lon, scale={domain=xrange}},
            y={:lat, scale={domain=yrange}}, color=inst, type="quantitative")

        display(_plot)
    end

end


text:
module TestGroupBy
    using TablesDemo
    using Base.Test

    tbl = Table(
        name = ["Niamh", "Roger", "Genevieve", "Aiden"],
        age = [27, 63, 26, 17],
        eye_color = ["green", "brown", "brown", "blue"]
    )
    _tbl = copy(tbl)

    group_indices = Dict{Any, Vector{Int}}([
        ("brown",false) => [3],
        ("green",true)  => [1],
        ("brown",true)  => [2],
        ("blue",false)  => [4],
    ])

    qrya = @query groupby(tbl, eye_color, age > 26)
    qryb = @query tbl |> groupby(eye_color, age > 26)

    tbl2a = collect(qrya)
    tbl2b = collect(qryb)
    @test isequal(tbl, _tbl)
    @test isequal(tbl2a, tbl2b)
    for (group_level, indices) in zip(
        [
            (Nullable("brown"), Nullable(false)),
            (Nullable("green"), Nullable(true)),
            (Nullable("brown"), Nullable(true)),
            (Nullable("blue"), Nullable(false))
        ],
        [
            [3],
            [1],
            [2],
            [4]
        ]
    )
        @test isequal(tbl2a.group_indices[group_level], indices)
    end

    @test isequal(tbl2a.groupbys, qrya.graph.args)

end


text:
"""Optimize a quantum control problem.

```julia
opt_result = optimize(problem; method=<method>, kwargs...)
```

optimizes towards a solution of given [`problem`](@ref ControlProblem) with the
given optimization `method`. Any keyword argument temporarily overrides the
corresponding keyword argument in `problem`.
"""
function optimize(problem::ControlProblem; method, kwargs...)
    if length(kwargs) > 0
        tempproblem = copy(problem)
        merge!(tempproblem.kwargs, kwargs)
        return optimize(tempproblem, method)
    else
        return optimize(problem, method)
    end
end

optimize(problem::ControlProblem, method::Symbol) = optimize(problem, Val(method))
#
# Note: Methods *must* be defined in the various optimization packages as e.g.
#
#   optimize(problem, method::Val{:krotov}) = optimize_krotov(problem)
#


text:
export Q_update!

"""
update the Q matrix in SparseMatrixCSC or CDS format. Q=sum_{i=1}^p AtA[i]rho[i]
if any of the rho[i] change
"""
function Q_update!(
                  Q            ::Union{Array{TF,2},SparseMatrixCSC{TF,TI}},
                  AtA          ::Union{Vector{Array{TF,2}},Vector{SparseMatrixCSC{TF,TI}}},
                  set_Prop,
                  rho          ::Vector{TF},
                  ind_updated  ::Vector{TI},
                  log_PARSDMM,
                  i            ::Integer,
                  Q_offsets=[] ::Vector{TI}
                  ) where {TF<:Real,TI<:Integer}

  if isempty(ind_updated) == false
    if typeof(Q)==SparseMatrixCSC{TF,TI}
      for ii=1:length(ind_updated)
          Q = Q + ( AtA[ind_updated[ii]] )*( rho[ind_updated[ii]]-log_PARSDMM.rho[i,ind_updated[ii]] );
      end
    else
      for ii=1:length(ind_updated)
        CDS_scaled_add!(Q,AtA[ind_updated[ii]],Q_offsets,set_Prop.AtA_offsets[ind_updated[ii]],rho[ind_updated[ii]] .- log_PARSDMM.rho[i,ind_updated[ii]] )
      end
    end
  end

return Q
end


text:
#=
	filename: BrunSan2010.jl
	author: Jonathan Payne

	cd("/Users/Juddy/GitHub/NYU/Projects/MacroFin/Julia")

	This file replicates Brunnermeir Sannikov in Julia

	This file computes the equilibrium in the setting of Brunnermeier and
	Sannikov under the assumptions that (1) experts can issue only debt and
	not equity and (2) the production set of experts is (a - iota, Phi(iota)
	- delta) and households, (a_ - iota, Phi(iota) - delta_), where the
	function [Phi iota] = investment(q) solves for Phi and iota

=#

#include("odeev.jl")
global qmax

using Roots

function rk4step(fdash, t, x, h)
    # Calculates the the nth step in the classical fourth order Runga-Kutta
    # procedure for solving a differential equation of the form:
    #   x'(t) = fdash(x, t),
    # Inputs:
    #   * fdash: user defined function
    #   * t: previous value of independent variable
    #   * x: previous value of dependent variable
    #   * h: increment
    # Outputs:
    #   * xnew: next step
    #   * tnew: new time
    k1 = h*fdash(x, t)
    k2 = h*fdash(x + k1/2, t + h/2)
    k3 = h*fdash(x + k2/2, t + h/2)
    k4 = h*fdash(x + k3, t + h)
    xnew = x + (k1 + 2*k2 + 2k3 + k4) / 6
    tnew = t + h
    return xnew, tnew
end

function rk4uniform(fdash, x0, tspan, N=100)
    # Approximate the solution of the initial value problem using the
    # classical fourth order Runge-Kutta method with uniform step sizes:
    #   * x'(t) = fdash(x, t),
    #   * x(t0) = x0
    # Inputs:
    #   * fdash: user defined function
    #   * t0: initial value of the independent variable
    #   * x0: initial value of the dependent variable(s) (column vector)
    #   * tf: final value of the independent variable
    #   * N: number of universally sized time steps to be taken to advance the
    #        the solution from t = t0 to t = tf
    # Outputs:
    #   * wi: vector / matrix containing the values of the approximate solution
    #         to the ode
    #   * ti: number of uniformly sized time steps to be taken to advance t
    t0, tf = tspan
    neqn = length(x0)
    ti = linspace(t0, tf, N+1)
    wi = zeros(neqn, N+1)
    wi[1:neqn, 1] = x0
    h = (tf - t0) / N
    for i = 1:N
        #t = t0; x = x0
        x0, t0 = rk4step(fdash, t0, x0, h)
        wi[1:neqn, i+1] = x0
    end
    return wi, ti
end

function odelay(func, y0, xspan,
                 events, TOLERANCE = 1e-6)#, fsolve_args=[], args...)
    #if fsolve_args == None
    #    fsolve_args = []
    #end
    N = length(xspan)
    x0 = xspan[1] # initial point

    X = [x0]
    Nf = length(y0)
    sol = zeros(N, Nf)
    sol[1,:] = float64(y0)
    TE, YE, IE = [], [], [] # to store where events occur

    # initial value of events
    e = zeros(length(events), N)
    for (i, event) in enumerate(events)
        e[i, 1], isterminal, direction = event(y0, x0)
    end

    # now we step through the integration
    for (i, x1) in enumerate(xspan[1:end-1])
        x2 = xspan[i + 1]
        f1 = sol[i, :]

        #fdash = func; x0 = f1'; tspan = [x1, x2]; N = 1
        f2, wz = rk4uniform(func, f1', [x1, x2], 1)#, args...)

        #X[i+1] = x2
        #println(X[i+1])
        X = push!(X, x2)
        sol[i+1,:] = f2[:, end]

        # check event functions. At each step we compute the event
        # functions, and check if they have changed sign since the
        # last step. If they changed sign, it implies a zero was
        # crossed.
        for (j, event) in enumerate(events)
            e[j, i + 1], isterminal, direction = event(sol[i + 1], X[i + 1])

            if ((e[j, i + 1] * e[j, i] < 0) || abs(e[j, i + 1]) < TOLERANCE || abs(e[j, i]) < TOLERANCE)
                xLt = X[end]       # Last point
                fLt = sol[i+1,:]
                eLt = e[j, i+1]

                function objective(x)
                    # evaluate ode from xLT to x
                    txspan = [xLt, x]
                    tempsol, wyz = rk4uniform(func, fLt', txspan, length(txspan))#, args...)
                    tsol = tempsol[:, end]
                    val, isterminal, direction = event(tsol, x)
                    return val[1]
                end
                xZ = fzero(objective, xLt)#, **fsolve_args)  # this should be the
                                              # value of x that makes
                                              # the event zero
                # now evaluate solution at this point, so we can
                # record the function values here.
                txspan = [xLt, xZ]
                tempsol, wyz = rk4uniform(func, fLt', txspan, length(txspan))#, args...)
                fZ = tempsol[:,end]

                vZ, isterminal, direction = event(fZ, xZ)

                COLLECTEVENT = false
                if direction == 0
                    COLLECTEVENT = true
                elseif (e[j, i + 1] > e[j, i] ) && direction == 1
                    COLLECTEVENT = true
                elseif (e[j, i + 1] < e[j, i] ) && direction == -1
                    COLLECTEVENT = true
                end

                if COLLECTEVENT
                    push!(TE, xZ)
                    push!(YE, fZ)
                    push!(IE, j)

                    if isterminal
                        X[end] = xZ[1]
                        sol[i+1,:] = fZ[1]
                        return X, sol[1:i+1,:], TE, YE, IE
                    end
                end

            end
        end
    end
    # at the end, return what we have
    return X, sol, TE, YE, IE

end

# Function finds the optimal investment rate iota and rate of capital creation
# Phi given the price q
function investment(q)
	# Input:
	# 	q = price of capital
	# Output:
	# 	Phi, iota
	Phi = (q - 1)/kappa
	iota = Phi + kappa*Phi^2/2
	return Phi, iota
end

a      = 0.11  # productivity of expert
a_     = 0.05  # productivity of household
rho    = 0.06  # expert discount rate
r      = 0.05  # riskfree rate of return
sigma  = 0.025 # volatility of capital ito process
delta  = 0.03  # expert discount rate
delta_ = 0.08  # household discount rate
kappa  = 10    # Adjustment cost parameter, such that investment
                # ...iota = Phi + kappa Phi^2/2 is required
                # ...to build capital Phi, is defined by investment(q)

# Some preparation: determine first best Q, and the worst Q
# This algorithm solves the following fixed point problem iteratively:
# q = max_{iota(q)} (a - iota)/(r - (Phi(iota(q)) - delta))
QL = 0 # initial lower bound on price
QR = 10 # initial upper bound on prices
for iter = 1:50
    qmax = (QL + QR)/2
    Phi, iota = investment(qmax)
    value = (a - iota)/(r + delta - Phi)
    if iota > a # iota is outside the bounds 0 < iota < a
        QR = qmax # iota is too large because the price is too large
    elseif value > qmax # qmax is to the right of fixed point q*
        QL = qmax # price too large
    elseif r + delta < Phi # Adjustment cost greater than output
        println("first-best price is infinite") # warn that first-best price is
        # infinite, but allow computation to continue
        QR = qmax  # without allowing q to grow above the level where
        # Phi > r + delta
    else
        QR = qmax
    end
end

# determine q_, the value of capital at eta = 0
QL = 0 # initial lower bound on price
QR = 10 # initial upper bound on prices
q_ = -1
for iter = 1:50
    q_ = (QL + QR)/2
    Phi, iota = investment(q_)
    value = (a_ - iota)/(r + delta_ - Phi)
    if iota > a_
        QR = q_ # price is too large
    elseif value < q_  # price still too large
        QR = q_
    else
        QL = q_
    end
end

# evntfcn
# ------------------------------
# [value,isterminal,direction] = evntfcn(eta, F) returns three
# variables used by ode45 to determine when to terminate integration
function event1(F, eta)
    value = qmax - F[3] # difference between qmax and q,
    isterminal = true # terminate computation in all three cases
    direction = 0 # event occurs whether we get there from above or below
    return value, isterminal, direction
end

function event2(F, eta)
    value = F[2] # first derivative of theta
    isterminal = true
    direction = 0
    return value, isterminal, direction
end

function event3(F, eta)
    value = F[4] # first derivative of q
    isterminal = true
    direction = 0
    return value, isterminal, direction
end
events = [event1, event2, event3]

# fnct
# ------------------------------
# This funciton implements proposition II.4. That is, it runs an algorithm
# for finding psi(eta) and coefficients of the ODEs for q(eta) and theta(eta)
# given (eta, q(eta), q'(eta), theta(eta), theta'(eta))
#
# Inputs:
# eta = scalar
# f = [theta, theta', q, q'], a 4x1 vector
# r, rho, a, a_, delta, delta_, sigma = parameters
#
# Outputs:
# fp = the derivative of f with respect to eta
# dyn = [psi, sigma_eta*eta, sigma_q, mu_eta*eta, mu_q, iota, leverage, rk, r_k]
#
# assumes that the production set of experts is (a - iota, Phi(iota) - delta)
# and households, (a_ - iota, Phi(iota) - delta_), where the function
# [Phi iota] = investment(q) solves for Phi and iota

function fnct(eta, f, r, rho, a, a_, delta, delta_, sigma)
    # search for psi between eta (lower bound) and
    # min(f(3)/f(4) + eta, 1) (upper bound)
    Phi, iota = investment(f[3])

    psi = 0
    sigma_eta_eta,sigma_q,sigma_theta,risk_premium,household_premium = 0,0,0,0,0
    psi_L = eta; psi_R = min(f[3]/f[4] + eta, 1)

    for n = 1:50
        psi = (psi_L + psi_R)/2
        amplification = 1 - f[4]/f[3]*(psi - eta)

        # VOLATILITY COMPUTATION
        sigma_eta_eta = sigma*(psi - eta)/amplification  # sigma_eta *times* eta
        sigma_q = sigma_eta_eta*f[4]/f[3]
        sigma_theta = sigma_eta_eta*f[2]/f[1]
        risk_premium = - sigma_theta*(sigma + sigma_q)

        household_premium = (a_ - a)/f[3] + delta - delta_ + risk_premium

        if household_premium > 0 # households want to hold more
            psi_R = psi
        else
            psi_L = psi
        end
    end
    (psi, sigma_eta_eta, sigma_q,
     sigma_theta, risk_premium, household_premium) = (
     float64(psi), float64(sigma_eta_eta), float64(sigma_q),
     float64(sigma_theta), float64(risk_premium), float64(household_premium))

    mu_q = r - (a - iota)/f[3] - Phi + delta - sigma*sigma_q + risk_premium
    mu_eta_eta = (-(psi - eta)*(sigma + sigma_q)*(sigma + sigma_q + sigma_theta)
                  + eta*(a - iota)/f[3] + eta*(1 - psi)*(delta_ - delta))
    qpp = 2*(mu_q*f[3] - f[4]*mu_eta_eta)/sigma_eta_eta^2
    thetapp = 2*((rho - r)*f[1] - f[2]*mu_eta_eta)/sigma_eta_eta^2

    fp = [f[2], thetapp, f[4], qpp];

    leverage = psi/eta
    rk = r + risk_premium
    r_k = r + household_premium

    dyn = [psi, sigma_eta_eta, sigma_q, mu_eta_eta, mu_q, iota, leverage, rk, r_k]
    return fp, dyn

end

# MAIN PART OF THE CODE:
# solves the system of ODE's starting with boundary conditions F0, i.e.
# theta(0) = 1, theta'(0) = #large negative number# and q(0) = q_;
# searching for q'(0) such that q'(eta*) and theta'(eta*) reach 0
# at the same boundary eta*

etaspan = linspace(0,1,50)
#F0 = [1 -1e+10 q_ 0]'  # [theta(0), theta'(0), q(0), q'(0)]
F0 = [1 -1000 q_ 0]'

# note that if theta(eta) satisfies the system of ODE's then so does
# const*theta(eta), for any positive constant const
# hence we set theta(0) = 1, and then normalize

odefun(f, eta) = fnct(eta, f, r, rho, a, a_, delta, delta_, sigma)[1]

QL = 0; QR = 1000#QR = 1e+15;
for iter = 1:50
    F0[4] = (QL + QR)/2  # this is q'(0)
    func = odefun; y0 = F0; xspan = etaspan
    [etaout,fout,TE,YE,IE] = odelay(odefun, F0, etaspan, events)
    if IE[1] == 3 # if q'(eta) has reached zero, we
        QL[1] = F0[4]  # increase q'(0)
    else        # if q(eta) reached qmax or theta'(0) reached 0
        QR[1] = F0[4]  # we reduce q'(0)
    end

end

# here we are basically done... let me just compute all other variables
# from the same function fnct
# N = length(etaout)
# dynout = zeros(N, 9)
# for n = 1:N
#     [fp dynout(n,:)] = fnct(etaout[n], fout[n,:], r, rho, a, a_, delta, delta_, sigma);
# end

# # normalize theta, to make sure that theta(eta*) = 1
# normalization = fout(N,1);
# fout(:,1:2) = fout(:,1:2)/normalization;

max_stars_repo_name: jepayne/NYU

max_issues_repo_name: jepayne/NYU

max_forks_repo_name: jepayne/NYU


text:
module Id

include("Sphere.jl")
import .Sphere: swal

import Polynomials: ChebyshevT
import HDF5: h5read

export set_gaussian!, set_qnm!

"""
Initial gaussian initial data for psi.

function set_gaussian(
      f,
      p,
      spin::Int64,
      mv::Int64,
      l_ang::Int64,
      ru::Float64,
      rl::Float64,
      width::Float64,
      amp::ComplexF64,
      cl::Float64,
      Rv::Vector{Float64},
      Yv::Vector{Float64}
   )
"""
function set_gaussian!(
      f,
      p,
      spin::Int64,
      mv::Int64,
      l_ang::Int64,
      ru::Float64,
      rl::Float64,
      width::Float64,
      amp::ComplexF64,
      cl::Float64,
      Rv::Vector{Float64},
      Yv::Vector{Float64}
   )
   @assert f.mv == mv
   @assert p.mv == mv

   nx, ny = f.nx, f.ny

   max_val = 0.0

   for j=1:ny
      for i=1:nx
         r = (cl^2)/Rv[i]

         bump = 0.0
         if ((r<ru) && (r>rl))
            bump = exp(-1.0*width/(r-rl))*exp(-2.0*width/(ru-r))
         end

         f.n[i,j]  = (((r-rl)/width)^2) * (((ru-r)/width)^2) * bump
         f.n[i,j] *= swal(spin,mv,l_ang,Yv[j])

         p.n[i,j] = 0.0

         max_val = max(abs(f.n[i,j]),max_val)
      end
   end

   ## rescale

   for j=1:ny
      for i=1:nx
         f.n[i,j] *= amp / max_val

         f.np1[i,j] = f.n[i,j]
         p.np1[i,j] = p.n[i,j]
      end
   end
   return nothing
end

"""
Initial qnm initial data for psi, read in
from HDF5 file.

function set_qnm(
      f,
      p,
      spin::Int64,
      mv::Int64,
      l_ang::Int64,
      ru::Float64,
      rl::Float64,
      width::Float64,
      amp::ComplexF64,
      cl::Float64,
      Rv::Vector{Float64},
      Yv::Vector{Float64}
   )
"""
function set_qnm(
      f,
      p,
      s::Int64,
      l::Int64,
      mv::Int64,
      n::Int64,
      a::Float64,
      amp::Float64,
      Rv::Vector{Float64},
      Yv::Vector{Float64}
   )
   @assert f.mv == mv
   @assert p.mv == mv
   nx, ny = f.nx, f.ny

   aval = round(digits=12,a)

   qnmpath = dirname(pwd())*"/qnm"
   h5f = h5read(
         qnmpath*"/prec1024_nr$(nr)_s$(s)_m$(mv)_n$(n).h5",
         "[a=$(aval),l=$(l)]"
        )
   rpoly = ChebyshevT(h5f["radial_coef"])
   lpoly = ChebyshevT(h5f["angular_coef"])
   lmin  = max(abs(s),abs(mv))
   lvals = [i+lmin for i in range(length(lpoly))]
   for j=1:ny
      for i=1:nx
         f.n[i,j]  = roply(Rv[i])
         f.n[i,j] *= sum(
            [lpoly[i]*swal(spin,mv,(i-1)+lmin,Yv[j])
             for i in 0:length(lpoly)
            ]
         )
      end
   end
   ## rescale
   for j=1:ny
      for i=1:nx
         f.n[i,j] *= amp / max_val
         f.np1[i,j] = f.n[i,j]
      end
   end
   ## p = f,t = -iωf
   omega = h5f["omega"]
   for j=1:ny
      for i=1:nx
         p.n[i,j]   = -im*omega*f.n[i,j]
         p.np1[i,j] = p.n[i,j]
      end
   end
   return nothing
end

end


text:
# module Grammars
module Parsimonious

include("Util.jl")
include("Nodes.jl")
include("Expressions.jl")

export govisit, debug_govisit, Grammar, lookup, rule_grammar, parse, unicode, match,  @grammar_mstr

# You have to import everything you want to overload and export
# import Expressions: parse, unicode, match

type Grammar
    rules::String
    default_rule::Expression
    exprs::Dict
end

# TODO: Grammar <: Associative
# for now, a lookup method
lookup(g::Grammar, key) = g.exprs[key]

# As an interface, calling with similar parameters in a different order could
# not be more fucking confusing. TODO: sort it out
function Grammar(rule_grammar::Grammar, rule_syntax::String;
    default_rule=nothing)
    exprs, first_rule = _expressions_from_rules(rule_grammar, rule_syntax)
    default_rule = is(default_rule,nothing) ? first_rule : exprs[default_rule]
    Grammar(rule_syntax, default_rule, exprs)
end

Grammar(rule_syntax::String) = Grammar(rule_grammar, rule_syntax)

function Grammar() # Bootstrapping
    exprs, default_rule = _expressions_from_rules(rule_syntax)
    Grammar(rule_syntax, default_rule, exprs)
end

parse(grammar::Grammar, text::String, pos=1) = parse(grammar.default_rule, text, pos)
match(grammar::Grammar, text::String, pos=1) = match(grammar.default_rule, text, pos)

rule_syntax = p"""
    # Ignored things (represented by _) are typically hung off the end of the
    # leafmost kinds of nodes. Literals like "/" count as leaves.

    rules = _ rule+
    rule = label equals expression
    equals = "=" _
    literal = spaceless_literal _

    # So you can't spell a regex like `~"..." ilm`:
    """ * """ # TODO: barf backslash handling
    spaceless_literal = ~"\\"[^\\"\\\\]*(?:\\\\.[^\\"\\\\]*)*\\""is /
                        ~"'[^'\\\\]*(?:\\\\.[^'\\\\]*)*'"is
    """ * p""" # /barf
    expression = ored / sequence / term
    or_term = "/" _ term
    ored = term or_term+
    sequence = term term+
    not_term = "!" term _
    lookahead_term = "&" term _
    term = not_term / lookahead_term / quantified / atom
    quantified = atom quantifier
    atom = reference / literal / regex / parenthesized
    regex = "~" spaceless_literal ~"[ilmsux]*"i _
    parenthesized = "(" _ expression ")" _
    quantifier = ~"[*+?]" _
    reference = label !equals

    # A subsequent equal sign is the only thing that distinguishes a label
    # (which begins a new rule) from a reference (which is just a pointer to a
    # rule defined somewhere else):
    label = ~"[a-zA-Z_][a-zA-Z_0-9]*" _

    # _ = ~"\s*(?:#[^\r\n]*)?\s*"
    _ = meaninglessness*
    meaninglessness = ~"\s+" / comment
    comment = ~"#[^\r\n]*"
"""

function boot_expressions()
    # Return a stripped-down rules which may be used to parse
    # rule_syntax to generate the complete rule_grammar
    comment = Regex(p"#[^\r\n]*", name="comment")
    meaninglessness = OneOf(Regex(p"\s+"), comment, name="meaninglessness")
    _ = ZeroOrMore(meaninglessness, name="_")
    equals = Sequence(Literal("="), _, name="equals")
    label = Sequence(Regex("[a-zA-Z_][a-zA-Z_0-9]*"), _, name="label")
    reference = Sequence(label, Not(equals), name="reference")
    quantifier = Sequence(Regex("[*+?]"), _, name="quantifier")
    spaceless_literal = Regex(p""" "[^"\\]*(?:\\.[^"\\]*)*" """[2:end-1],
        options="is",
        name="spaceless_literal")
    literal = Sequence(spaceless_literal, _, name="literal")
    regex = Sequence(Literal("~"),
                     literal,
                     Regex("[ilmsux]*", options="i"),
                     _,
                     name="regex")
    atom = OneOf(reference, literal, regex, name="atom")
    quantified = Sequence(atom, quantifier, name="quantified")

    term = OneOf(quantified, atom, name="term")
    not_term = Sequence(Literal("!"), term, _, name="not_term")
    term = OneOf(not_term, quantified, atom, name="term")

    sequence = Sequence(term, OneOrMore(term), name="sequence")
    or_term = Sequence(Literal("/"), _, term, name="or_term")
    ored = Sequence(term, OneOrMore(or_term), name="ored")
    expression = OneOf(ored, sequence, term, name="expression")
    rule = Sequence(label, equals, expression, name="rule")
    rules = Sequence(_, OneOrMore(rule), name="rules")

    # TODO: Is there a way to get a list of all of the local variables so
    # I don't have to make this dictionary manually?
    #exprnames = split("""
    #    comment meaninglessness _ equals label reference quantifier
    #    spaceless_literal literal regex atom quantified term not_term
    #    term sequence or_term ored expression rule rules""")
    #exprs = {n=>@deref(n) for n in exprnames}

    exprs ={
        "comment" => comment,
        "meaninglessness" => meaninglessness,
        "_" => _,
        "equals" => equals,
        "label" => label,
        "reference" => reference,
        "quantifier" => quantifier,
        "spaceless_literal" => spaceless_literal,
        "literal" => literal,
        "regex" => regex,
        "atom" => atom,
        "quantified" => quantified,
        "term" => term,
        "not_term" => not_term,
        "term" => term,
        "sequence" => sequence,
        "or_term" => or_term,
        "ored" => ored,
        "expression" => expression,
        "rule" => rule,
        "rules" => rules
        }
end

function boot_grammar()
    exprs = boot_expressions()
    Grammar("", exprs["rules"], exprs)
end

function _expressions_from_rules(grammar::Grammar, rules::String)
    tree = parse(grammar, rules)
    DEFUNCT_visit(RuleVisitor(), tree)
end

# make boot grammar
function _expressions_from_rules(rule_syntax::String)
    rules = boot_expressions()["rules"]
    rule_tree = parse(rules, rule_syntax)
    exprs, default_rule = DEFUNCT_visit(RuleVisitor(), rule_tree)
    return exprs, default_rule
end

type RuleVisitor <: NodeVisitor end

# generic_visit
DEFUNCT_visit(v::RuleVisitor, n::LeafNode) = n
DEFUNCT_visit(v::RuleVisitor, n::ParentNode, visited_children) = visited_children

# TODO: tuples not arrays for visited_children
DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:expression}, visited_children) = visited_children[1]
DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:term}, visited_children) = visited_children[1]
DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:atom}, visited_children) = visited_children[1]

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:parenthesized}, visited_children)
    left_paren, _1, expression, right_paren, _2 = visited_children
    expression
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:quantifier}, visited_children)
    symbol, _ = visited_children
    symbol
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:quantified}, visited_children)
    atom, quantifier = visited_children
    quantifier = nodetext(quantifier)[1]
    if quantifier == '?'
        return Optional(atom)
    elseif quantifier == '*'
        return ZeroOrMore(atom)
    elseif quantifier == '+'
        return OneOrMore(atom)
    end
    error("How is '" * string(quantifier) * "' a quantifier to you?")
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:lookahead_term}, visited_children)
    ampersand, term, _ = visited_children
    Lookahead(term)
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:not_term}, visited_children)
    exclamation, term, _ = visited_children
    Not(term)
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:rule}, visited_children)
    label, equals, expression = visited_children
    expression.name = label
    expression
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:sequence}, visited_children)
    term, terms = visited_children
    Sequence(term, terms...)
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:ored}, visited_children)
    term, terms, = visited_children
    OneOf(term, terms...)
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:or_term}, visited_children)
    slash, _, term = visited_children
    term
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:label}, visited_children)
    name, _ = visited_children
    nodetext(name)
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:reference}, visited_children)
    label, not_equals = visited_children
    LazyReference("", label)
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:regex}, visited_children)
    tilde, literal, flags, _ = visited_children
    pattern = literal.literal
    flags = lowercase(nodetext(flags))
    Parsimonious.Regex(pattern, options=flags)
end

function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:literal}, visited_children)
    spaceless_literal, _ = visited_children
    spaceless_literal
end

type UndefinedLabel <: Exception label end

showerror(io::IO, e::UndefinedLabel) = print(io,"The label $(e.label) was never defined.")

function _resolve_refs(rule_map, expr, unwalked_names, walking_names)
    #Return an expression with all its lazy references recursively
    #resolved.

    #Resolve any lazy references in the expression ``expr``, recursing into
    #all subexpressions. Populate ``rule_map`` with any other rules (named
    #expressions) resolved along the way. Remove from ``unwalked_names`` any
    #which were resolved.

    #:arg walking_names: The stack of labels we are currently recursing
    #    through. This prevents infinite recursion for circular refs.

    # If it's a top-level (named) expression and we've already walked it,
    # don't walk it again:
    exprname = expr.name
    if !isempty(exprname) && !in(exprname, unwalked_names)
        # unwalked_names started out with all the rule names in it, so, if
        # this is a named expr and it isn't in there, it must have been
        # resolved.
        return rule_map[exprname]
    # If not, resolve it:
    elseif isa(expr, LazyReference)
        label = expr.label
        if !in(label, walking_names)
            # We aren't already working on traversing this label:
            local reffed_expr
            try
                reffed_expr = rule_map[label]
            catch e
                if isa(e, KeyError)
                    # println(rule_map)
                    throw(UndefinedLabel(expr))
                end
                rethrow(e)
            end
            rule_map[label] = _resolve_refs(rule_map, reffed_expr, unwalked_names, Set(walking_names..., label))

            # If we recurse into a compound expression, the remove()
            # happens in there. But if this label points to a non-compound
            # expression like a literal or a regex or another lazy
            # reference, we need to do this here:
            delete!(unwalked_names, label)
        end
        return rule_map[label]
    else
        members = []
        try
            members = getfield(expr, :members)
        end
        if !isempty(members)
            expr.members = [_resolve_refs(rule_map, m, unwalked_names, walking_names) for m in members]
        end
        if exprname != ""
            delete!(unwalked_names, exprname)
        end
        return expr
    end
end

# return dictionary of expressions and a rule name
function DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:rules}, visited_children)
    _, rules = visited_children
    rule_map = {expr.name => expr for expr in rules}
    unwalked_names = Set(collect(keys(rule_map))...)
    while !isempty(unwalked_names)
        rule_name = first(unwalked_names)
        rule_map[rule_name] = _resolve_refs(rule_map, rule_map[rule_name], unwalked_names, Set(rule_name))
        delete!(unwalked_names, rule_name)
    end

    return rule_map, rules[1]
end

DEFUNCT_visit(v::RuleVisitor, n::ParentNode{:spaceless_literal}, _) = Literal(nodetext(n)[2:end-1])
DEFUNCT_visit(v::RuleVisitor, n::LeafNode{:spaceless_literal}) = Literal(nodetext(n)[2:end-1])

rule_grammar = Grammar()  # Level 1 -- Bootstrap grammar parses rule_symtax
rule_grammar = Grammar(rule_grammar, rule_syntax)  # Level 2 -- generated rule_grammar parses rule_syntax
rule_grammar = Grammar(rule_grammar, rule_syntax)  # Level 3 -- re-generated rule_grammar parses rule_syntax

# TODO: Stupid name
function unicode(g::Grammar)
    """Return a rule string that, when passed to the constructor, would
    reconstitute the grammar."""
# TODO: This got hacked up in debugging. It was prettier to start.
    sexprs = [as_rule(g.default_rule)]
    for expr in g.exprs
        if !is(expr[2], g.default_rule)
            push!(sexprs, as_rule(expr[2]))
        end
    end
    join([sexprs], "\n")
end

# cheap
cleanspace(s) = replace(replace(s, "\\n", "\n"), "\\t", "\t")
macro grammar_mstr(s) Grammar(rule_grammar, cleanspace(s)) end

end # module


text:
# ---
# title: 37. Sudoku Solver
# id: problem37
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, Backtracking
# link: <https://leetcode.com/problems/sudoku-solver/description/>
# hidden: true
# ---
#
# Write a program to solve a Sudoku puzzle by filling the empty cells.
#
# A sudoku solution must satisfy **all of  the following rules**:
#
#   1. Each of the digits `1-9` must occur exactly once in each row.
#   2. Each of the digits `1-9` must occur exactly once in each column.
#   3. Each of the digits `1-9` must occur exactly once in each of the 9 `3x3` sub-boxes of the grid.
#
# The `'.'` character indicates empty cells.
#
#
#
# **Example 1:**
#
# ![](https://upload.wikimedia.org/wikipedia/commons/thumb/f/ff/Sudoku-
# by-L2G-20050714.svg/250px-Sudoku-by-L2G-20050714.svg.png)
#
#
#
#     Input: board = [["5","3",".",".","7",".",".",".","."],["6",".",".","1","9","5",".",".","."],[".","9","8",".",".",".",".","6","."],["8",".",".",".","6",".",".",".","3"],["4",".",".","8",".","3",".",".","1"],["7",".",".",".","2",".",".",".","6"],[".","6",".",".",".",".","2","8","."],[".",".",".","4","1","9",".",".","5"],[".",".",".",".","8",".",".","7","9"]]
#     Output: [["5","3","4","6","7","8","9","1","2"],["6","7","2","1","9","5","3","4","8"],["1","9","8","3","4","2","5","6","7"],["8","5","9","7","6","1","4","2","3"],["4","2","6","8","5","3","7","9","1"],["7","1","3","9","2","4","8","5","6"],["9","6","1","5","3","7","2","8","4"],["2","8","7","4","1","9","6","3","5"],["3","4","5","2","8","6","1","7","9"]]
#     Explanation:  The input board is shown above and the only valid solution is shown below:
#
#     ![](https://upload.wikimedia.org/wikipedia/commons/thumb/3/31/Sudoku-by-L2G-20050714_solution.svg/250px-Sudoku-by-L2G-20050714_solution.svg.png)
#
#
#
#
# **Constraints:**
#
#   * `board.length == 9`
#   * `board[i].length == 9`
#   * `board[i][j]` is a digit or `'.'`.
#   * It is **guaranteed** that the input board has only one solution.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
T = Complex{Float64}
N = 51
sampling_rate = 100.
input_file_prefix = "data/"
input_sacpz = input_file_prefix * "SAC_PZs_NZ_CRLZ_HHZ"
input_sacpz_without_zeros = input_file_prefix * "SAC_PZs_NZ_HHZ_10"
expect_output_type = Vector{Complex{Float64}}

@testset "read response" begin
    @testset "read from SACPZ" begin
        freq_response = read_resp_from_sacpz(input_sacpz, sampling_rate, N)
        @test size(freq_response, 1) == N
        @test typeof(freq_response) == expect_output_type
    end

    @testset "read from SACPZ without zeros" begin
    end
end


text:
using IndexedTables
using SurrealNumbers
convert(SurrealFinite, 2) * convert(SurrealFinite, 1)

ns = [0; 1; 2; 3; 4; 5]
n = length(ns)
m = [100; 100; 100; 30; 1; 1; 1]
x = Array{SurrealFinite,1}(n)
t = Array{Float64,1}(n)
bytes = Array{Float64,1}(n)
gctime = Array{Float64,1}(n)
for i=1:n
    println(" 2 x ", ns[i])
    t[i] = 0
    b = 0
    g = 0
    mem = 0
    val = 0
    for j=1:m[i]
        val, time, b, g, mem = @timed convert(SurrealFinite, 2) * convert(SurrealFinite, ns[i])
        t[i] += time
    end
    x[i] = val
    bytes[i] = b
    gctime[i] = g
end
s = size.(x)
d = depth_av.(x)
f = n_zeros.(x)
t = t./m[1:n]
bytes = bytes./m[1:n]

t_p_ms = 1.0e6 * t ./ s

t1 = table(@NT( i=0:n-1,
                size=s,
                depth=d,
                leaves=f,
                time=t,
                time_per_n=t_p_ms,
                bytes=bytes);
           pkey = [:i])


text:
export costL2TV, costL2TVTV2, costL2TV2, costTV, costTV2,costIntrICTV12

@doc doc"""
    costIntrICTV12(M,f,u,v,α,β)

computes the intrinsic infimal convolution model, where the addition is replaced
by a mid point approach and the two functions involved are [`costTV2`](@ref)
and [`costTV`](@ref). The model reads

```math
E(u,v) =
  \frac{1}{2}\sum_{i\in\mathcal G}
    d_{\mathcal M}\bigl(g(\frac{1}{2},v_i,w_i),f_i\bigr)
  +\alpha\bigl( \beta\mathrm{TV}(v) + (1-\beta)\mathrm{TV}_2(w) \bigr).
```
"""
function costIntrICTV12(M::mT,f::P,u::P,v::P,α,β) where {mT <: Manifold, P <: MPoint}
  return 1/2*distance(M,midPoint(M,u,v,f),f)^2 + α*β*costTV(M,u) + α * (1-β) * costTV2(M,v)
end

@doc doc"""
    costL2TV(M,f,α,x)

compute the $\ell^2$-TV functional on the [`Power`](@ref) manifold `M` for given
(fixed) data `f` (on `M`), a nonnegative weight `α`, and evaluated at `x` (on `M`),
i.e.

```math
E(x) = d_{\mathcal M}^2(f,x) + \alpha \operatorname{TV}(x)
```

# See also
[`costTV`](@ref)
"""
costL2TV(M::Power,f::PowPoint,α::Number,x::PowPoint) =
  1/2*distance(M,f,x)^2 + α*costTV(M,x)

@doc doc"""
    costL2TVTV2(M,f,α,β,x)

compute the $\ell^2$-TV-TV2 functional on the [`Power`](@ref) manifold `M` for
given (fixed) data `f` (on `M`), nonnegative weight `α`, `β`, and evaluated
at `x` (on `M`), i.e.

```math
E(x) = d_{\mathcal M}^2(f,x) + \alpha\operatorname{TV}(x)
  + \beta\operatorname{TV}_2(x)
```

# See also
[`costTV`](@ref), [`costTV2`](@ref)
"""
costL2TVTV2(M::Power,f::PowPoint,α::Number,β::Number,x::PowPoint) =
  1/2*distance(M,f,x)^2 + α*costTV(M,x) + β*costTV2(M,x)

@doc doc"""
    costL2TV2(M,f,β,x)

compute the $\ell^2$-TV2 functional on the [`Power`](@ref) manifold `M`
for given data `f`, nonnegative parameter `β`, and evaluated at `x`, i.e.

```math
E(x) = d_{\mathcal M}^2(f,x) + \beta\operatorname{TV}_2(x)
```

# See also
[`costTV2`](@ref)
"""
costL2TV2(M::Power,f::PowPoint,β::Number,x::PowPoint) =
    1/2*distance(M,f,x)^2 + β*costTV2(M,x)

@doc doc"""
    costTV(M,x,p)

compute the $\operatorname{TV}^p$ functional for a tuple `pT` of [`MPoint`](@ref)s
on a [`Manifold`](@ref) `M`, i.e.

```math
E(x_1,x_2) = d_{\mathcal M}^p(x_1,x_2), \quad x_1,x_2\in\mathcal M
```

# See also
[`gradTV`](@ref), [`proxTV`](@ref)
"""
function costTV(M::mT,x::Tuple{P,P},p::Int=1) where {mT <: Manifold, P <: MPoint}
  return distance(M,x[1],x[2])^p
end
@doc doc"""
    costTV(M,x [,p=2,q=1])

compute the $\operatorname{TV}^p$ functional for data `x`on the [`Power`](@ref)
manifold `M`, i.e. $\mathcal M = \mathcal N^n$, where $n\in\mathbb N^k$ denotes
the dimensions of the data `x`.
Let $\mathcal I_i$ denote the forward neighbors, i.e. with $\mathcal G$ as all
indices from $\mathbf{1}\in\mathbb N^k$ to $n$ we have
$\mathcal I_i = \{i+e_j, j=1,\ldots,k\}\cap \mathcal G$.
The formula reads

```math
E^q(x) = \sum_{i\in\mathcal G}
  \bigl( \sum_{j\in \mathcal I_i} d^p_{\mathcal M}(x_i,x_j) \bigr)^{q/p}.
```

# See also
[`gradTV`](@ref), [`proxTV`](@ref)
"""
function costTV(M::Power, x::PowPoint, p::Int=1, q::Int=1)
  R = CartesianIndices(M.powerSize)
  d = length(M.powerSize)
  maxInd = last(R)
  cost = fill(0.,M.powerSize)
  for k in 1:d # for all directions
    ek = CartesianIndex(ntuple(i  ->  (i==k) ? 1 : 0, d) ) #k th unit vector
    for i in R # iterate over all pixel
      j = i+ek # compute neighbor
      if all( map(<=, j.I, maxInd.I)) # is this neighbor in range?
        cost[i] += costTV(M.manifold,(x[i],x[j]),p) # Compute TV on these
      end
    end
  end
  cost = (cost).^(1/p)
  if q > 0
    return sum(cost.^q)^(1/q)
  else
    return cost
  end
end
@doc doc"""
    costTV2(M,(x1,x2,x3) [,p=1])

compute the $\operatorname{TV}_2^p$ functional for the 3-tuple of points
`(x1,x2,x3)`on the [`Manifold`](@ref) `M`. Denote by

```math
  \mathcal C = \bigl\{ c\in \mathcal M \ |\ g(\tfrac{1}{2};x_1,x_3) \text{ for some geodesic }g\bigr\}
```

the set of mid points between $x_1$ and $x_3$. Then the functionr reads

$d_2^p(x_1,x_2,x_3) = \min_{c\in\mathcal C} d_{\mathcal M}(c,x_2).$

# See also
[`gradTV2`](@ref), [`proxTV2`](@ref)
"""
function costTV2(M::mT,pointTuple::Tuple{P,P,P},p::Int=1) where {mT <: Manifold, P <: MPoint}
  # note that here midPoint returns the closest to x2 from the e midpoints between x1 x3
  return 1/p*distance(M,midPoint(M,pointTuple[1],pointTuple[3]),pointTuple[2])^p
end
@doc doc"""
    costTV2(M,x [,p=1])

compute the $\operatorname{TV}_2^p$ functional for data `x` on the
[`Power`](@ref) manifoldmanifold `M`, i.e. $\mathcal M = \mathcal N^n$,
where $n\in\mathbb N^k$ denotes the dimensions of the data `x`.
Let $\mathcal I_i^{\pm}$ denote the forward and backward neighbors, respectively,
i.e. with $\mathcal G$ as all indices from $\mathbf{1}\in\mathbb N^k$ to $n$ we
have $\mathcal I^\pm_i = \{i\pm e_j, j=1,\ldots,k\}\cap \mathcal I$.
The formula then reads

```math
E(x) = \sum_{i\in\mathcal I,\ j_1\in \mathcal I^+_i,\ j_2\in \mathcal I^-_i}
d^p_{\mathcal M}(c_i(x_{j_1},x_{j_2}), x_i),
```

where $c_i(\cdot,\cdot)$ denotes the mid point between its two arguments that is
nearest to $x_i$.

# See also
[`gradTV2`](@ref), [`proxTV2`](@ref)
"""
function costTV2(M::Power, x::PowPoint, p::Int=1, Sum::Bool=true)
  R = CartesianIndices(M.powerSize)
  d = length(M.powerSize)
  minInd, maxInd = first(R), last(R)
  cost = fill(0.,M.powerSize)
  for k in 1:d # for all directions
    ek = CartesianIndex(ntuple(i  ->  (i==k) ? 1 : 0, d) ) #k th unit vector
    for i in R # iterate over all pixel
      jF = i+ek # compute forward neighbor
      jB = i-ek # compute backward neighbor
      if all( map(<=, jF.I, maxInd.I) ) && all( map(>=, jB.I, minInd.I)) # are neighbors in range?
        cost[i] += costTV2( M.manifold, (x[jB], x[i], x[jF]),p ) # Compute TV on these
      end
    end # i in R
  end # directions
  if p != 1
    cost = (cost).^(1/p)
  end
  if Sum
    return sum(cost)
  else
    return cost
  end
end


text:
mutable struct RBMSplitLUT{VT} <: NNLookUp
    θ::VT
end

lookup(net::RBMSplit) = RBMSplitLUT(net)
RBMSplitLUT(net::RBMSplit) = RBMSplitLUT(similar(net.b))

function set_lookup!(lt::RBMSplitLUT, net::RBMSplit, c::RBMSplitCache, σr_r, σc_r)
    σr = c.σr; copy!(σr, σr_r)
    σc = c.σc; copy!(σc, σc_r)

    mul!(lt.θ,    net.Wr, σr)
    mul!(c.θ_tmp, net.Wc, σc)
    lt.θ .+= net.b + c.θ_tmp

    return lt
end

function update_lookup!(lt::RBMSplitLUT, net::RBMSplit, c::RBMSplitCache, σr_r, σc_r,
                         changes_r, changes_c)
    σr = c.σr; copy!(σr, σr_r)
    σc = c.σc; copy!(σc, σc_r)

    θ = lt.θ

    for (i, nv) = changes_r
        θ .+= net.Wr[:,i] .* (nv - σr_r[i])
    end
    for (i, nv) = changes_c
        θ .+= net.Wc[:,i] .* (nv - σc_r[i])
    end
    return lt
end

function logψ(net::RBMSplit, c::RBMSplitCache, lt::RBMSplitLUT, σr_r, σc_r,
              changes_r, changes_c)
    σr = c.σr; copy!(σr, σr_r)
    σc = c.σc; copy!(σc, σc_r)
    θ  = c.θ;
    logℒθ = c.logℒθ
    ∂logℒθ = c.∂logℒθ

    θ .= lt.θ
    for (i, nv) = changes_r
        θ .+= net.Wr[:,i] .* (nv - σr_r[i])
        σr[i] = nv
    end
    for (i, nv) = changes_c
        θ .+= net.Wc[:,i] .* (nv - σc_r[i])
        σc[i] = nv
    end

    logℒθ .= logℒ.(θ)
    logψ = dot(σr, net.ar) + dot(σc, net.ac) + sum(logℒθ)
    return logψ
end

function logψ_and_∇logψ!(∇logψ, net::RBMSplit, c::RBMSplitCache, lt::RBMSplitLUT,
                         σr_r, σc_r, changes_r, changes_c)
    # copy the states to complex valued states for the computations.
    σr = c.σr; copy!(σr, σr_r)
    σc = c.σc; copy!(σc, σc_r)

    θ      = c.θ
    logℒθ  = c.logℒθ
    ∂logℒθ = c.∂logℒθ

    θ .= lt.θ
    for (i, nv) = changes_r
        θ .+= net.Wr[:,i] .* (nv - σr_r[i])
        σr[i] = nv
    end
    for (i, nv) = changes_c
        θ .+= net.Wc[:,i] .* (nv - σc_r[i])
        σc[i] = nv
    end

    logℒθ  .= logℒ.(θ)
    ∂logℒθ .= ∂logℒ.(θ)

    ∇logψ.ar .= σr
    ∇logψ.ac .= σc
    ∇logψ.b  .= ∂logℒθ
    ∇logψ.Wr .= ∂logℒθ .* transpose(σr)
    ∇logψ.Wc .= ∂logℒθ .* transpose(σc)

    logψ = dot(σr,net.ar) + dot(σc,net.ac) + sum(logℒθ)
    return logψ
end

function Δ_logψ(net::RBMSplit, c::RBMSplitCache, lt::RBMSplitLUT, σr_r, σc_r,
                     changes_r, changes_c)
    σr = c.σr; copy!(σr, σr_r)
    σc = c.σc; copy!(σc, σc_r)
    θ = c.θ
    logℒθ = c.logℒθ
    ∂logℒθ = c.∂logℒθ

    logvaldiff = zero(eltype(θ))

    if !isempty(changes_r) || !isempty(changes_c)
        θ .= lt.θ
        logℒθ .= logℒ.(θ)
        sumlogℒθ_old = sum(logℒθ)

        for (i, nv) = changes_r
            θ .+= net.Wr[:,i] .* (nv - σr_r[i])
            logvaldiff += net.ar[i] * (nv - σr_r[i])
        end
        for (i, nv) = changes_c
            θ .+= net.Wc[:,i] .* (nv - σc_r[i])
            logvaldiff += net.ac[i] * (nv - σc_r[i])
        end

        logℒθ .= logℒ.(θ)
        logvaldiff += sum(logℒθ) - sumlogℒθ_old
    end

    return logvaldiff
end

function Δ_logψ_and_∇logψ!(∇logψ, net::RBMSplit, c::RBMSplitCache, lt::RBMSplitLUT,
                         σr_r, σc_r, changes_r, changes_c)
    σr = c.σr; copy!(σr, σr_r)
    σc = c.σc; copy!(σc, σc_r)
    θ = c.θ;
    logℒθ  = c.logℒθ
    ∂logℒθ = c.∂logℒθ

    θ .= lt.θ

    logvaldiff = zero(eltype(θ))
    if !isempty(changes_r) || !isempty(changes_c)
        logℒθ .= logℒ.(θ)
        sumlogℒθ_old = sum(logℒθ)

        for (i, nv) = changes_r
            θ .+= net.Wr[:,i] .* (nv - σr_r[i])
            logvaldiff += net.ar[i] * (nv - σr_r[i])
            σr[i] = nv
        end
        for (i, nv) = changes_c
            θ .+= net.Wc[:,i] .* (nv - σc_r[i])
            logvaldiff += net.ac[i] * (nv - σc_r[i])
            σc[i] = nv
        end

        logℒθ .= logℒ.(θ)
        logvaldiff += sum(logℒθ) - sumlogℒθ_old
    end

    ∂logℒθ .= ∂logℒ.(θ)

    ∇logψ.ar .= σr
    ∇logψ.ac .= σc
    ∇logψ.b  .= ∂logℒθ
    ∇logψ.Wr .= ∂logℒθ .* transpose(σr)
    ∇logψ.Wc .= ∂logℒθ .* transpose(σc)

    return logvaldiff
end


text:
# parse readable numeric strings

parse_readable(::Type{T}, s::String, ch::Char) where {T <: Union{Integer, AbstractFloat}} =
    Base.parse(T, join(split(s,ch),""))

parse_readable(::Type{T}, s::String, ch1::Char, ch2::Char) where {T <: AbstractFloat} =
    Base.parse(T, join(split(s,(ch1,ch2)),""))

"""
how many times does char c occur in string s
"""
function count_char(s::String, c::Char)
    r = (c=='.') ? Regex("\\.") : Regex(string(c))
    return length( matchall(r,s) )
end


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

"""
    NearestMapping

A mapping method in which data points are assigned to their nearest
point in the domain.
"""
struct NearestMapping <: MappingMethod end

function Base.map(sdata, sdomain, targetvars, ::NearestMapping)
  @assert targetvars ⊆ name.(variables(sdata)) "target variables must be present in spatial data"

  # dictionary with mappings
  mappings = Dict(var => Dict{Int,Int}() for var in targetvars)

  # nearest neighbor search method
  neighbor = Vector{Int}(undef, 1)
  searcher = KNearestSearch(sdomain, 1)

  for ind in 1:nelements(sdata)
    # update datum coordinates
    pₒ = centroid(sdata, ind)

    # find nearest location in the domain
    search!(neighbor, pₒ, searcher)

    # save pair if there is data for variable
    for var in targetvars
      if !ismissing(sdata[var][ind])
        push!(mappings[var], neighbor[1] => ind)
      end
    end
  end

  mappings
end


text:
const ENUM_OUTER_RND = 1E-9
function f_init!(::RelaxMulEnum, g::DAT, b::RelaxCache{V,N,T}) where {V,N,T}
    b.use_apriori_mul = false
    fprop!(Relax(), g, b)
    xp = copy(val(b))
    vlbd = lbd(b)
    vubd = ubd(b)
    xl = copy(vlbd)
    xu = copy(vubd)
    vlbd .-= ENUM_OUTER_RND
    vubd .+= ENUM_OUTER_RND
    for k = 1:node_count(g)
        if !is_num(b, k)
            b._info[k].v = set(b, k)
        end
    end
    s = sparsity(g, 1)
    for i = 0:2^N-1
        s = last(bitstring(i), N)
        for (k,j) in enumerate(sparsity(g, 1))
            b.ic.v.x[j] = s[k] == '1' ? xl[j] : xu[j]
        end
        fprop!(Relax(), g, b)
        for k = node_count(g):-1:1
            if !is_num(b, k)
                b._info[k][i+1] = set(b, k)
            end
        end
    end
    b.ic.v.x .= xp
    vlbd .= xl
    vubd .= xu
    b.use_apriori_mul = true
    fprop!(Relax(), g, b)
    return
end

_cut_info(t::RelaxMulEnumInner, v, z, x) = z
_cut_info(t::RelaxMulEnum, v, z, x) = z

relax_info(s::RelaxMulEnumInner, n::Int, t::T) where T = MCBoxPnt{2^n,n,T}
relax_info(s::RelaxMulEnum, n::Int, t::T) where T = MCBoxPnt{2^n,n,T}

function estimator_extrema(x::MCBoxPnt{Q,N,T}, y::MCBoxPnt{Q,N,T}, s, dP) where {Q,N,T}
    xmax = maximum(cv, x.box)
    ymax = maximum(cv, y.box)
    xmin = minimum(cc, x.box)
    ymin = minimum(cc, y.box)
    return xmax, ymax, -xmin, -ymin
end

function estimator_under(xv, yv, x::MCBoxPnt{Q,N,T}, y::MCBoxPnt{Q,N,T}, s, dp, dP, p_rel, p_diam) where {Q,N,T}
    xv.cv, yv.cv, xv.cv_grad, yv.cv_grad
end

function estimator_over(xv, yv, x::MCBoxPnt{Q,N,T}, y::MCBoxPnt{Q,N,T}, s, dp, dP, p_rel, p_diam) where {Q,N,T}
    -xv.cc, -yv.cc, -xv.cc_grad, -yv.cc_grad
end

for d in ALL_ATOM_TYPES
    @eval function fprop!(t::RelaxMulEnum, v::Val{$d}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
        fprop!(Relax(), v, g, b, k)
    end
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
type Archive <: Associative
    filename::AbstractString
    bytes
    keys::Array{AbstractString,1}
    tempdir::AbstractString
    modified::Bool
end

function Archive(filename::AbstractString, bytes=nothing)

    # The command-line "zip" tool exepcts a file extension.
    @assert length(split(basename(filename),".")) > 1

    # Get list of files from archive...
    if isfile(filename)
        keys = [chomp(f) for f in readlines(`unzip -Z1 $filename`)]
        keys = filter(k->basename(k) != "", keys)
    else
        keys = AbstractString[]
    end

    Archive(filename, bytes, keys, mktempdir(), false)
end

# Open ZIP Archive from "filename".

open_zip(filename::AbstractString) = Archive(abspath(filename))
with_close(f::Function, io) = try f(io) finally close(io) end

# Open ZIP Archive from "bytes".

function open_zip(bytes::Array{UInt8,1})

    # Copy "bytes" into temporary .ZIP file...
    if !isempty(bytes)
        tmp, io = mktemp()
        write(io, bytes)
        close(io)
        filename = tmp * ".zip"
        mv(tmp, filename)
    else
        filename = tempname() * ".zip"
    end
    Archive(abspath(filename), bytes)
end

# "open_zip(archive) do..." syntax for exception safety.
open_zip(f::Function, archive) = with_close(f, open_zip(archive))

# Save changes to ZIP Archive.

function Base.close(z::Archive)

    # Add files from tempdir to the archive...
    if !isempty(readdir(z.tempdir))
        run(Cmd(`zip -q -r $(z.filename) .`, dir=z.tempdir))
        z.modified = true
    end

    # Write temporary file back to "bytes" array...
    if z.modified && z.bytes != nothing
        open(z.filename, "r") do io
            readbytes!(io, z.bytes, filesize(z.filename))
        end
        rm(z.filename)
    end

    rm(z.tempdir, recursive=true)
end

# Extract file from archive using Associative syntax: data = z[filename].

function Base.get(z::Archive, filename::AbstractString, default=nothing)

    @assert !ismatch(r"^/", filename)

    if !haskey(z, filename)
        return default
    end

    f = joinpath(z.tempdir, filename)
    if isfile(f)
        b = read(f)
    else
        b = read(`unzip -qc $(z.filename) $filename`)
    end
    return isvalid(ASCIIString, b) ? ASCIIString(b) :
           isvalid(UTF8String, b)  ? UTF8String(b)  : b
end

# Add file to archive using Associative syntax: z[filename] = data.

function Base.setindex!(z::Archive, data, filename::AbstractString)

    @assert !ismatch(r"^/", filename)

    # Write file to tempdir...
    mkpath(joinpath(z.tempdir, dirname(filename)))
    write(joinpath(z.tempdir, filename), data)

    # Add filename to "keys"...
    if !(filename in z.keys)
        push!(z.keys, filename)
    end

    return data
end

# Add contents of "dict" to archive.

function Base.merge!{T<:Associative}(z::Archive, dict::T)
    for (filename, data) in dict
        z[filename] = data
    end
end

# Read files from ZIP using iterator syntax.

Base.eltype(::Type{Archive}) =
    Tuple{AbstractString,Union{ByteString,Vector{UInt8},AbstractString}}

Base.keys(z::Archive) =   z.keys
Base.length(z::Archive) = length(z.keys)
Base.start(z::Archive) = start(z.keys)
Base.done(z::Archive, state) = done(z.keys, state)

function Base.next(z::Archive, state)
    (filename, state) = next(z.keys, state)
    return ((filename, get(z, filename)), state)
end

typealias FileOrBytes Union{AbstractString,Vector{UInt8}}

# Extract ZIP archive to "outputpath".
# Based on fhs/ZipFile.jl#16, thanks @timholy.

function unzip(archive::FileOrBytes, outputpath::AbstractString=pwd())
    open_zip(archive) do z
        run(Cmd(`unzip -qo $(z.filename)`, dir = outputpath))
    end
end

rm_archive(a::Array{UInt8,1}) = nothing
rm_archive(f::AbstractString) = !isfile(f) || rm(f)

# Create archive from filenames.

function create_zip(archive::FileOrBytes, files::Array)

    rm_archive(archive)
    open_zip(archive) do z
        run(`zip -q $(z.filename) $files`)
        z.modified = true
    end
end

# Create archive from "dict".

function create_zip{T<:Associative}(archive::FileOrBytes, dict::T)

    rm_archive(archive)
    open_zip(archive) do z
        for (filename, data) in dict
            z[filename] = data
        end
    end
end

# Create archive from (filename, data) tuples.

create_zip{T<:Tuple}(a::FileOrBytes, files::Array{T}) = create_zip(a, files...)

# Create archive from "files" and "data".

create_zip(a::FileOrBytes, files::Array, data::Array) = create_zip(a, zip(files, data)...)

# Create archive from filename => data pairs.

create_zip(archive::FileOrBytes, args...) = create_zip(archive, Dict(args))

# Use temporary memory buffer if "filename" or "io" are not provided.

function create_zip(arg, args...)

    buf = UInt8[]
    create_zip(buf, arg, args...)
    return buf
end


text:
"""
    Planet

Planetary parameters.
"""
module Planet

export molmass_dryair,
    R_d,
    kappa_d,
    cp_d,
    cv_d,
    ρ_cloud_liq,
    ρ_cloud_ice,
    molmass_water,
    molmass_ratio,
    R_v,
    cp_v,
    cp_l,
    cp_i,
    cv_v,
    cv_l,
    cv_i,
    T_freeze,
    T_min,
    T_max,
    T_icenuc,
    pow_icenuc,
    T_triple,
    T_0,
    LH_v0,
    LH_s0,
    LH_f0,
    e_int_v0,
    e_int_i0,
    press_triple,
    surface_tension_coeff,
    entropy_dry_air,
    entropy_water_vapor,
    entropy_reference_temperature,
    ρ_ocean,
    cp_ocean,
    planet_radius,
    day,
    Omega,
    grav,
    year_anom,
    orbit_semimaj,
    tot_solar_irrad,
    epoch,
    mean_anom_epoch,
    obliq_epoch,
    lon_perihelion_epoch,
    eccentricity_epoch,
    lon_perihelion,
    MSLP,
    T_surf_ref,
    T_min_ref

# Properties of dry air
""" Molecular weight dry air (kg/mol) """
function molmass_dryair end
""" Gas constant dry air (J/kg/K) """
function R_d end
""" Adiabatic exponent dry air """
function kappa_d end
""" Isobaric specific heat dry air """
function cp_d end
""" Isochoric specific heat dry air """
function cv_d end

# Properties of water
""" Density of liquid water (kg/m``^3``) """
function ρ_cloud_liq end
""" Density of ice water (kg/m``^3``) """
function ρ_cloud_ice end
""" Molecular weight (kg/mol) """
function molmass_water end
""" Molar mass ratio dry air/water """
function molmass_ratio end
""" Gas constant water vapor (J/kg/K) """
function R_v end
""" Isobaric specific heat vapor (J/kg/K) """
function cp_v end
""" Isobaric specific heat liquid (J/kg/K) """
function cp_l end
""" Isobaric specific heat ice (J/kg/K) """
function cp_i end
""" Isochoric specific heat vapor (J/kg/K) """
function cv_v end
""" Isochoric specific heat liquid (J/kg/K) """
function cv_l end
""" Isochoric specific heat ice (J/kg/K) """
function cv_i end
""" Freezing point temperature (K) """
function T_freeze end
""" Minimum temperature guess in saturation adjustment (K) """
function T_min end
""" Maximum temperature guess in saturation adjustment (K) """
function T_max end
""" Homogeneous nucleation temperature (K) """
function T_icenuc end
""" The assumed phase partition function power (-) """
function pow_icenuc end
""" Triple point temperature (K) """
function T_triple end
""" Reference temperature (K) """
function T_0 end
""" Latent heat vaporization at ``T_0`` (J/kg) """
function LH_v0 end
""" Latent heat sublimation at ``T_0`` (J/kg) """
function LH_s0 end
""" Latent heat of fusion at ``T_0`` (J/kg) """
function LH_f0 end
""" Specific internal energy of vapor at ``T_0`` (J/kg) """
function e_int_v0 end
""" Specific internal energy of ice at ``T_0`` (J/kg) """
function e_int_i0 end
""" Triple point vapor pressure (Pa) """
function press_triple end
""" Surface tension coefficient of water (J/m2) """
function surface_tension_coeff end

#=
The standard entropy value for dry air is computed based
on the reference data given in Lemmon et al. [2000], and
the standard entropy value for water vapor is based on the
reference data given in Chase [1998].
=#

""" Entropy of dry air J / (kg K)"""
function entropy_dry_air end
""" Entropy of water vapor J / (kg K)"""
function entropy_water_vapor end
""" Entropy reference temperature (K)"""
function entropy_reference_temperature end

# Properties of sea water
""" Reference density sea water (kg/m``^3``) """
function ρ_ocean end
""" Specific heat sea water (J/kg/K) """
function cp_ocean end

# Planetary parameters
""" Mean planetary radius (m) """
function planet_radius end
""" Length of day (s) """
function day end
""" Ang. velocity planetary rotation (1/s) """
function Omega end
""" Gravitational acceleration (m/s``^2``) """
function grav end
""" Length of anomalistic year (s) """
function year_anom end
""" ngth of semimajor orbital axis (m) """
function orbit_semimaj end
""" Total solar irradiance (W/m``^2``) """
function tot_solar_irrad end
""" Time of  epoch (J2000) (s) """
function epoch end
""" Mean anomaly at the epoch (radians) """
function mean_anom_epoch end
""" Orbital obliquity at the epoch (radians) """
function obliq_epoch end
""" Longitude of perihelion at the epoch (radians),
measured relative to vernal equinox (i.e., the longitude of perihelion is the angle
subtended at the Sun by the orbital arc from vernal equinox to perihelion). """
function lon_perihelion_epoch end
""" Orbital eccentricity at the epoch """
function eccentricity_epoch end
""" Longitude of perihelion (radians), measured relative to vernal equinox.
The calculation of the mean anomaly is formulated such that the vernal equinox
is fixed in the calendar. However, this requires tracking both the current longitude of perihelion
as well as the reference longitude of perihelion at the reference time (epoch). """
function lon_perihelion end
""" Mean sea level pressure (Pa) """
function MSLP end
""" Mean surface temperature (K) in reference state """
function T_surf_ref end
""" Minimum temperature (K) in reference state """
function T_min_ref end

end

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
module DocTestsTest
using Documenter, Test

println("="^50)
@info("Testing `doctest = :fix`")
mktempdir(@__DIR__) do dir
    srcdir = mktempdir(dir)
    builddir = mktempdir(dir)
    cp(joinpath(@__DIR__, "broken.md"), joinpath(srcdir, "index.md"))
    cp(joinpath(@__DIR__, "broken.jl"), joinpath(srcdir, "src.jl"))
    include(joinpath(srcdir, "src.jl"))
    @eval using .Foo
    # fix up
    makedocs(sitename="-", modules = [Foo], source = srcdir, build = builddir, doctest = :fix)
    # test that strict = true works
    makedocs(sitename="-", modules = [Foo], source = srcdir, build = builddir, strict = true)
    # also test that we obtain the expected output
    @test read(joinpath(srcdir, "index.md"), String) ==
          read(joinpath(@__DIR__, "fixed.md"), String)
    @test read(joinpath(srcdir, "src.jl"), String) ==
          read(joinpath(@__DIR__, "fixed.jl"), String)
end
@info("Done testing `doctest = :fix`")
println("="^50)

end # module


text:
@testset "SysEx events" begin
    validtestvalues = [
        # Structure: dT (variable length), begin byte (0xF0), length (variable length type), data, end byte (0xF7)
        # dT is a variable length value which is handled by readvariablelength and not readsysexevent. We don't need to test it here, which
        # is why it's set to 0 in each case.
        ([0x00, 0xF0, 0x05, 0x11, 0x21, 0x53, 0x1F, 0xF7], MIDI.SysexEvent(0, [0x11, 0x21, 0x53, 0x1F])),
        ([0x00, 0xF0, 0x07, 0x25, 0x61, 0x23, 0x5B, 0x00, 0x02, 0xF7], MIDI.SysexEvent(0, [0x25, 0x61, 0x23, 0x5B, 0x00, 0x02]))
    ]

    invalidtestvalues = [
        ([0x00, 0xF0, 0x05, 0x11, 0xF7, 0x53, 0x1F, 0xF7], ErrorException),
        ([0x00, 0xF0, 0x01, 0x11, 0xF7, 0x53, 0x1F, 0xF7], ErrorException)
    ]

    @testset "sysex event" begin
        @testset "it should correctly read a sysex event from raw data" begin
            for (input, output) in validtestvalues
                result = MIDI.readsysexevent(Int(input[1]), IOBuffer(input[2:length(input)]))
                @test result.dT == output.dT && result.data == output.data
            end
        end
    end

    @testset "it should correctly write a sysex event to raw data" begin
        for (output, input) in validtestvalues
            buf = IOBuffer()
            MIDI.writeevent(buf, input)
            @test take!(buf) == output
        end
    end

    @testset "it should fail to read a sysex event when given invalid data" begin
        for (input, errtype) in invalidtestvalues
            @test_throws errtype MIDI.readsysexevent(Int(input[1]), IOBuffer(input[2:length(input)]))
        end
    end
end


text:
# `target` will typically a `Function` or `Type` but could also be a function-like object
struct Patch{T}
    target::T
    alternate::Function
end

macro patch(expr::Expr)
    def = splitdef(expr)

    if haskey(def, :name) && haskey(def, :body)
        target = def[:name]
    elseif !haskey(def, :name)
        throw(ArgumentError("Function definition must be a named function"))
    else
        throw(ArgumentError("Function definition must not be an empty function"))
    end

    # Include the target function name in the patch to make stack traces easier to read.
    # If the provided target uses a fully-qualified reference we'll just extract the name
    # of the function (e.g `Base.foo` -> `foo`).
    target_name = if Meta.isexpr(target, :.)
        target.args[2].value
    else
        target
    end

    def[:name] = gensym(string(target_name, "_patch"))
    alternate = combinedef(def)

    # We need to evaluate the alternate function in the context of the `@patch` macro in
    # order to support closures.
    return esc(:(Mocking.Patch($target, $alternate)))
end

struct PatchEnv
    mapping::Dict{Any,Vector{Function}}
    debug::Bool
end

function PatchEnv(patches, debug::Bool=false)
    pe = PatchEnv(debug)
    apply!(pe, patches)
    return pe
end

PatchEnv(debug::Bool=false) = PatchEnv(Dict{Any,Vector{Function}}(), debug)

function apply!(pe::PatchEnv, p::Patch)
    alternate_funcs = get!(Vector{Function}, pe.mapping, p.target)
    push!(alternate_funcs, p.alternate)
    return pe
end

function apply!(pe::PatchEnv, patches)
    for p in patches
        apply!(pe, p)
    end
end

function apply(body::Function, pe::PatchEnv)
    prev_pe = get_active_env()
    set_active_env(pe)
    try
        return body()
    finally
        set_active_env(prev_pe)
    end
end

function apply(body::Function, patches; debug::Bool=false)
    return apply(body, PatchEnv(patches, debug))
end

const PATCH_ENV = Ref{PatchEnv}(PatchEnv())
set_active_env(pe::PatchEnv) = (PATCH_ENV[] = pe)
get_active_env() = PATCH_ENV[]


text:
#! /usr/bin/julia

# Rosetta Code, Textonyms

const tcode = (Regex=>Char)[r"A|B|C|Ä|Å|Á|Â|Ç" => '2',
                            r"D|E|F|È|Ê|É" => '3',
                            r"G|H|I|Í" => '4',
                            r"J|K|L" => '5',
                            r"M|N|O|Ó|Ö|Ô|Ñ" => '6',
                            r"P|Q|R|S" => '7',
                            r"T|U|V|Û|Ü" => '8',
                            r"W|X|Y|Z" => '9']

function tpad(str::IOStream)
    tnym = (String=>Array{String,1})[]
    for w in eachline(str)
        w = chomp(w)
        t = uppercase(w)
        for (k,v) in tcode
            t = replace(t, k, v)
        end
        t = replace(t, r"\D", '1')
        tnym[t] = [get(tnym, t, String[]), w]
    end
    return tnym
end

dname = "/usr/share/dict/american-english"
DF = open(dname, "r")
tnym = tpad(DF)
close(DF)

println("The character to digit mapping is done according to")
println("these regular expressions (following uppercase conversion):")
for k in sort(collect(keys(tcode)), by=x->tcode[x])
    println("    ", tcode[k], " -> ", k)
end
println("Unmatched non-digit characters are mapped to 1")

println()
print("There are ", sum(map(x->length(x), values(tnym))))
println(" words in ", dname)
println("  which can be represented by the digit key mapping.")
print("They require ", length(keys(tnym)))
println(" digit combinations to represent them.")
print(sum(map(x->length(x)>1, values(tnym))))
println(" digit combinations represent Textonyms.")

println()
println("The degeneracies of telephone key encodings are:")
println("  Words Encoded   Number of codes")
dgen = zeros(maximum(map(x->length(x), values(tnym))))
for v in values(tnym)
    dgen[length(v)] += 1
end
for (i, d) in enumerate(dgen)
    println(@sprintf "%10d  %15d" i d)
end

println()
dgen = length(dgen) - 2
println("Codes mapping to ", dgen, " or more words:")
for (k, v) in tnym
    dgen <= length(v) || continue
    println(@sprintf "%7s (%2d) %s" k length(v) join(v, ", "))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
###
### Bipartite graph utilities
###

"""
    maximal_matching(s::SystemStructure, eqfilter=eq->true, varfilter=v->true) -> Matching

Find equation-variable maximal bipartite matching. `s.graph` is a bipartite graph.
"""
BipartiteGraphs.maximal_matching(s::SystemStructure, eqfilter=eq->true, varfilter=v->true) =
    maximal_matching(s.graph, eqfilter, varfilter)

function error_reporting(state, bad_idxs, n_highest_vars, iseqs)
    io = IOBuffer()
    if iseqs
        error_title = "More equations than variables, here are the potential extra equation(s):\n"
        out_arr = equations(state)[bad_idxs]
    else
        error_title = "More variables than equations, here are the potential extra variable(s):\n"
        out_arr = state.fullvars[bad_idxs]
    end

    Base.print_array(io, out_arr)
    msg = String(take!(io))
    neqs = length(equations(state))
    if iseqs
        throw(ExtraEquationsSystemException(
            "The system is unbalanced. "
            * "There are $n_highest_vars highest order derivative variables "
            * "and $neqs equations.\n"
            * error_title
            * msg
        ))
    else
        throw(ExtraVariablesSystemException(
            "The system is unbalanced. "
            * "There are $n_highest_vars highest order derivative variables "
            * "and $neqs equations.\n"
            * error_title
            * msg
        ))
    end
end

###
### Structural check
###
function check_consistency(state::TearingState)
    fullvars = state.fullvars
    @unpack graph, var_to_diff = state.structure
    n_highest_vars = count(v->length(outneighbors(var_to_diff, v)) == 0, vertices(var_to_diff))
    neqs = nsrcs(graph)
    is_balanced = n_highest_vars == neqs

    if neqs > 0 && !is_balanced
        varwhitelist = var_to_diff .== nothing
        var_eq_matching = maximal_matching(graph, eq->true, v->varwhitelist[v]) # not assigned
        # Just use `error_reporting` to do conditional
        iseqs = n_highest_vars < neqs
        if iseqs
            eq_var_matching = invview(complete(var_eq_matching)) # extra equations
            bad_idxs = findall(isnothing, @view eq_var_matching[1:nsrcs(graph)])
        else
            bad_idxs = findall(isequal(unassigned), var_eq_matching)
        end
        error_reporting(state, bad_idxs, n_highest_vars, iseqs)
    end

    # This is defined to check if Pantelides algorithm terminates. For more
    # details, check the equation (15) of the original paper.
    extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist; map(collect, edges(var_to_diff))])
    extended_var_eq_matching = maximal_matching(extended_graph)

    unassigned_var = []
    for (vj, eq) in enumerate(extended_var_eq_matching)
        if eq === unassigned
            push!(unassigned_var, fullvars[vj])
        end
    end

    if !isempty(unassigned_var) || !is_balanced
        io = IOBuffer()
        Base.print_array(io, unassigned_var)
        unassigned_var_str = String(take!(io))
        errmsg = "The system is structurally singular! " *
                 "Here are the problematic variables: \n" *
                 unassigned_var_str
        throw(InvalidSystemException(errmsg))
    end

    return nothing
end

###
### BLT ordering
###

"""
    find_var_sccs(g::BipartiteGraph, assign=nothing)

Find strongly connected components of the variables defined by `g`. `assign`
gives the undirected bipartite graph a direction. When `assign === nothing`, we
assume that the ``i``-th variable is assigned to the ``i``-th equation.
"""
function find_var_sccs(g::BipartiteGraph, assign=nothing)
    cmog = DiCMOBiGraph{true}(g, Matching(assign === nothing ? Base.OneTo(nsrcs(g)) : assign))
    sccs = Graphs.strongly_connected_components(cmog)
    foreach(sort!, sccs)
    return sccs
end

function sorted_incidence_matrix(sys, val=true; only_algeqs=false, only_algvars=false)
    var_eq_matching, var_scc = algebraic_variables_scc(sys)
    s = structure(sys)
    @unpack fullvars, graph = s
    g = graph
    varsmap = zeros(Int, ndsts(graph))
    eqsmap = zeros(Int, nsrcs(graph))
    varidx = 0
    eqidx = 0
    for vs in scc, v in vs
        eq = var_eq_matching[v]
        if eq !== unassigned
            eqsmap[eq] = (eqidx += 1)
            varsmap[var] = (varidx += 1)
        end
    end
    for i in diffvars_range(s)
        varsmap[i] = (varidx += 1)
    end
    for i in dervars_range(s)
        varsmap[i] = (varidx += 1)
    end
    for i in 1:nsrcs(graph)
        if eqsmap[i] == 0
            eqsmap[i] = (eqidx += 1)
        end
    end

    I = Int[]
    J = Int[]
    for eq in 𝑠vertices(g)
        only_algeqs && (isalgeq(s, eq) || continue)
        for var in 𝑠neighbors(g, eq)
            only_algvars && (isalgvar(s, var) || continue)
            i = eqsmap[eq]
            j = varsmap[var]
            (iszero(i) || iszero(j)) && continue
            push!(I, i)
            push!(J, j)
        end
    end
    #sparse(I, J, val, nsrcs(g), ndsts(g))
    sparse(I, J, val)
end

###
### Structural and symbolic utilities
###

function find_eq_solvables!(state::TearingState, ieq; may_be_zero=false, allow_symbolic=false)
    fullvars = state.fullvars
    @unpack graph, solvable_graph = state.structure
    eq = equations(state)[ieq]
    term = value(eq.rhs - eq.lhs)
    to_rm = Int[]
    for j in 𝑠neighbors(graph, ieq)
        var = fullvars[j]
        isinput(var) && continue
        a, b, islinear = linear_expansion(term, var)
        a = unwrap(a)
        islinear || continue
        if a isa Symbolic
            allow_symbolic || continue
            add_edge!(solvable_graph, ieq, j)
            continue
        end
        (a isa Number) || continue
        if a != 0
            add_edge!(solvable_graph, ieq, j)
        else
            if may_be_zero
                push!(to_rm, j)
            else
                @warn "Internal error: Variable $var was marked as being in $eq, but was actually zero"
            end
        end
    end
    for j in to_rm
        rem_edge!(graph, ieq, j)
    end
end

function find_solvables!(state::TearingState; allow_symbolic=false)
    @assert state.structure.solvable_graph === nothing
    eqs = equations(state)
    graph = state.structure.graph
    state.structure.solvable_graph = BipartiteGraph(nsrcs(graph), ndsts(graph))
    for ieq in 1:length(eqs)
        find_eq_solvables!(state, ieq; allow_symbolic)
    end
    return nothing
end

# debugging use
function reordered_matrix(sys, torn_matching)
    s = TearingState(sys)
    complete!(s.structure)
    @unpack graph = s.structure
    eqs = equations(sys)
    nvars = ndsts(graph)
    max_matching = complete(maximal_matching(graph))
    torn_matching = complete(torn_matching)
    sccs = find_var_sccs(graph, max_matching)
    I, J = Int[], Int[]
    ii = 0
    M = Int[]
    solved = BitSet(findall(torn_matching .!== unassigned))
    for vars in sccs
        append!(M, filter(in(solved), vars))
        append!(M, filter(!in(solved), vars))
    end
    M = invperm(vcat(M, setdiff(1:nvars, M)))
    for vars in sccs
        e_solved = [torn_matching[v] for v in vars if torn_matching[v] !== unassigned]
        for es in e_solved
            isdiffeq(eqs[es]) && continue
            ii += 1
            js = [M[x] for x in 𝑠neighbors(graph, es) if isalgvar(s.structure, x)]
            append!(I, fill(ii, length(js)))
            append!(J, js)
        end

        e_residual = setdiff([max_matching[v] for v in vars if max_matching[v] !== unassigned], e_solved)
        for er in e_residual
            isdiffeq(eqs[er]) && continue
            ii += 1
            js = [M[x] for x in 𝑠neighbors(graph, er) if isalgvar(s.structure, x)]
            append!(I, fill(ii, length(js)))
            append!(J, js)
        end
    end
    # only plot algebraic variables and equations
    sparse(I, J, true)
end

"""
    uneven_invmap(n::Int, list)

returns an uneven inv map with length `n`.
"""
function uneven_invmap(n::Int, list)
    rename = zeros(Int, n)
    for (i, v) in enumerate(list)
        rename[v] = i
    end
    return rename
end

function torn_system_jacobian_sparsity(sys)
    state = get_tearing_state(sys)
    state isa TearingState || return nothing
    s = structure(sys)
    graph = state.structure.graph
    fullvars = state.fullvars

    states_idxs = findall(!isdifferential, fullvars)
    var2idx = Dict{Int,Int}(v => i for (i, v) in enumerate(states_idxs))
    I = Int[]; J = Int[]
    for ieq in 𝑠vertices(graph)
        for ivar in 𝑠neighbors(graph, ieq)
            nivar = get(var2idx, ivar, 0)
            nivar == 0 && continue
            push!(I, ieq)
            push!(J, nivar)
        end
    end
    return sparse(I, J, true)
end

###
### Nonlinear equation(s) solving
###

@noinline nlsolve_failure(rc) = error("The nonlinear solver failed with the return code $rc.")

function numerical_nlsolve(f, u0, p)
    prob = NonlinearProblem{false}(f, u0, p)
    sol = solve(prob, NewtonRaphson())
    rc = sol.retcode
    rc === :DEFAULT || nlsolve_failure(rc)
    # TODO: robust initial guess, better debugging info, and residual check
    sol.u
end


text:
# https://github.com/sglyon/PlotlyJS.jl

# ------------------------------------------------------------------------------

function plotlyjs_syncplot(plt::Plot{PlotlyJSBackend})
    plt[:overwrite_figure] && closeall()
    plt.o = PlotlyJS.plot()
    traces = PlotlyJS.GenericTrace[]
    for series_dict in plotly_series(plt)
        plotly_type = pop!(series_dict, :type)
        push!(traces, PlotlyJS.GenericTrace(plotly_type; series_dict...))
    end
    PlotlyJS.addtraces!(plt.o, traces...)
    layout = plotly_layout(plt)
    w, h = plt[:size]
    PlotlyJS.relayout!(plt.o, layout, width = w, height = h)
    return plt.o
end

# ------------------------------------------------------------------------------

for (mime, fmt) in (
    "application/pdf" => "pdf",
    "image/png"       => "png",
    "image/svg+xml"   => "svg",
    "image/eps"       => "eps",
)
    @eval _show(io::IO, ::MIME{Symbol($mime)}, plt::Plot{PlotlyJSBackend}) = PlotlyJS.savefig(io, plotlyjs_syncplot(plt), format = $fmt)
end

# Use the Plotly implementation for json and html:
_show(io::IO, mime::MIME"application/vnd.plotly.v1+json", plt::Plot{PlotlyJSBackend}) = plotly_show_js(io, plt)

html_head(plt::Plot{PlotlyJSBackend}) = plotly_html_head(plt)
html_body(plt::Plot{PlotlyJSBackend}) = plotly_html_body(plt)

_show(io::IO, ::MIME"text/html", plt::Plot{PlotlyJSBackend}) = write(io, embeddable_html(plt))

_display(plt::Plot{PlotlyJSBackend}) = display(plotlyjs_syncplot(plt))

function PlotlyJS.WebIO.render(plt::Plot{PlotlyJSBackend})
    return PlotlyJS.WebIO.render(plotlyjs_syncplot(plt))
end

function closeall(::PlotlyJSBackend)
    if !isplotnull() && isa(current().o, PlotlyJS.SyncPlot)
        close(current().o)
    end
end

Base.showable(::MIME"application/prs.juno.plotpane+html", plt::Plot{PlotlyJSBackend}) = true


text:
# Pair - collection of key and value

key, value = 1, 2;
p = Pair(key, value);
println(p);
println(typeof(p));

p1 = "Key" => 1;
println(p1);

# Dictionary - collection of key-value pairs
key1 = "Key1";

d1 = Dict(key1 => "Value1", 2 => 3, p1);
println(d1);

# Adding to Dictionary

d1[3] = "Key3" => 2;
println(d1);
d1[4] = "four";
println(d1);

# List Comprehension dictionary

d2 = Dict([i => i for i = 1:10]);
println(d2);

# Getting keys of dictionary

k1 = keys(d1);
println(k1);
println(typeof(k1));

# Check for a key
println(haskey(d1, "Key1"));

# Getting Values of dictionary
v1 = values(d1);
println(v1);
println(typeof(v1));

# Getting array from dictionary

a1 = collect(d1);
println(a1);

# Retrieve from dictionary
v2 = d1["Key1"];

# Merging dictionaries

d2 = Dict(1 => "one", 2 => "two");
d3 = Dict(2 => "one", 3 => "three"); # For conflicts, keeps most recent value
d4 = merge(d2, d3);
println(d4);

# Delete mapping and return dictionary

delete!(d1, "Key1");
println(d1);

# Delete mapping and return deleted KV pair, or default value
println(pop!(d1, "Key1", "Default"));

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming