text:
"""
    VehicleState
A default type to represent an agent physical state (position, velocity).
It contains the position in the global frame, Frenet frame and the longitudinal velocity
# constructors
    VehicleState(posG::VecSE2{Float64}, v::Float64)
    VehicleState(posG::VecSE2{Float64}, roadway::Roadway, v::Float64)
    VehicleState(posG::VecSE2{Float64}, lane::Lane, roadway::Roadway, v::Float64)
    VehicleState(posF::Frenet, roadway::Roadway, v::Float64)

# fields
- `posG::VecSE2{Float64}` global position
- `posF::Frenet` lane relative position
- `v::Float64` longitudinal velocity
"""
struct VehicleState
    posG::VecSE2{Float64} # global
    posF::Frenet # lane-relative frame
    v::Float64
end

VehicleState() = VehicleState(VecSE2(), NULL_FRENET, NaN)
VehicleState(posG::VecSE2{Float64}, v::Float64) = VehicleState(posG, NULL_FRENET, v)
VehicleState(posG::VecSE2{Float64}, roadway::Roadway, v::Float64) = VehicleState(posG, Frenet(posG, roadway), v)
VehicleState(posG::VecSE2{Float64}, lane::Lane, roadway::Roadway, v::Float64) = VehicleState(posG, Frenet(posG, lane, roadway), v)
VehicleState(posF::Frenet, roadway::Roadway, v::Float64) = VehicleState(get_posG(posF, roadway), posF, v)

Base.show(io::IO, s::VehicleState) = print(io, "VehicleState(", s.posG, ", ", s.posF, ", ", @sprintf("%.3f", s.v), ")")

function Base.write(io::IO, ::MIME"text/plain", s::VehicleState)
    @printf(io, "%.16e %.16e %.16e", s.posG.x, s.posG.y, s.posG.θ)
    @printf(io, " %d %.16e %d %d", s.posF.roadind.ind.i, s.posF.roadind.ind.t, s.posF.roadind.tag.segment, s.posF.roadind.tag.lane)
    @printf(io, " %.16e %.16e %.16e", s.posF.s, s.posF.t, s.posF.ϕ)
    @printf(io, " %.16e", s.v)
end
function Base.read(io::IO, ::MIME"text/plain", ::Type{VehicleState})
    tokens = split(strip(readline(io)), ' ')
    i = 0
    posG = VecSE2(parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]))
    roadind = RoadIndex(CurveIndex(parse(Int, tokens[i+=1]), parse(Float64, tokens[i+=1])),
                        LaneTag(parse(Int, tokens[i+=1]), parse(Int, tokens[i+=1])))
    posF = Frenet(roadind, parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]))
    v = parse(Float64, tokens[i+=1])
    return VehicleState(posG, posF, v)
end

"""
    Vec.lerp(a::VehicleState, b::VehicleState, t::Float64, roadway::Roadway)
Perform linear interpolation of the two vehicle states. Returns a VehicleState.
"""
function Vec.lerp(a::VehicleState, b::VehicleState, t::Float64, roadway::Roadway)
    posG = lerp(a.posG, b.posG, t)
    v = lerp(a.v, b.v, t)
    VehicleState(posG, roadway, v)
end

"""
    get_vel_s(s::VehicleState)
returns the longitudinal velocity (along the lane)
"""
get_vel_s(s::VehicleState) = s.v * cos(s.posF.ϕ) # velocity along the lane
"""
    get_vel_t(s::VehicleState)
returns the lateral velocity (⟂ to lane)
"""
get_vel_t(s::VehicleState) = s.v * sin(s.posF.ϕ) # velocity ⟂ to lane

"""
    move_along(vehstate::VehicleState, roadway::Roadway, Δs::Float64;
    ϕ₂::Float64=vehstate.posF.ϕ, t₂::Float64=vehstate.posF.t, v₂::Float64=vehstate.v)

returns a vehicle state after moving vehstate of a length Δs along its lane.
"""
function move_along(vehstate::VehicleState, roadway::Roadway, Δs::Float64;
    ϕ₂::Float64=vehstate.posF.ϕ, t₂::Float64=vehstate.posF.t, v₂::Float64=vehstate.v
    )

    roadind = move_along(vehstate.posF.roadind, roadway, Δs)
    try
        footpoint = roadway[roadind]
    catch
        println(roadind)
    end
    footpoint = roadway[roadind]
    posG = convert(VecE2, footpoint.pos) + polar(t₂, footpoint.pos.θ + π/2)
    posG = VecSE2(posG.x, posG.y, footpoint.pos.θ + ϕ₂)
    VehicleState(posG, roadway, v₂)
end

"""
    Vehicle
A specific instance of the Entity type defined in Records to represent Vehicles with
state `VehicleState` , definition `VehicleDef` and id `Int64`
"""
const Vehicle = Entity{VehicleState,VehicleDef,Int64}

# XXX Should this go in features
"""
    get_center(veh::Vehicle)
returns the position of the center of the vehicle
"""
get_center(veh::Vehicle) = veh.state.posG
"""
    get_footpoint(veh::Vehicle)
returns the position of the footpoint of the vehicle
"""
get_footpoint(veh::Vehicle) = veh.state.posG + polar(veh.state.posF.t, veh.state.posG.θ-veh.state.posF.ϕ-π/2)

"""
    get_front(veh::Vehicle)
returns the position of the front of the vehicle
"""
get_front(veh::Vehicle) = veh.state.posG + polar(veh.def.length/2, veh.state.posG.θ)

"""
    get_rear(veh::Vehicle)
returns the position of the rear of the vehicle
"""
get_rear(veh::Vehicle) = veh.state.posG - polar(veh.def.length/2, veh.state.posG.θ)


text:
# QN.jl
# Copyright 2017 Koki Kitai.
#
# Quasi-Newton methods
using Optim
using LineSearches
include("lnsrch.jl")

"""
    qnBalancing{T<:AbstractFloat}(P::Array{T,2}; ϵ=1e-6, max_iter=1000, show_trace=false, show_every=100)

Matrix balancing by E-projection with LBFGS.

# Parameters
* P - Input matrix
* ϵ - Threshold for the error of η-coordinate
* max_iter - Max number of iterations for optimization
* show_trace - Whether to print trace
* show_every - Trace information is printed every `show_every`'th iteration
"""
function qnBalancing{T<:AbstractFloat}(P::Array{T,2}; ϵ=1e-6, max_iter=1000, show_trace=false, show_every=100, limited=false)
    A = copy(P)
    ylen, xlen = size(A)
    cumθ_y, cumθ_x = zeros(ylen), zeros(xlen)
    η = zeros(xlen+ylen-1)
    β = genβ(A::Array{T,2})
    β = β[2:end]
    η_target = genTargetη(A, β)

    function f(θ)
        cumθ!(cumθ_x, cumθ_y, θ, β, size(A))
        local s = 0.0
        for j = 1:xlen, i = 1:ylen
            if P[i,j] > 1e-15 && A[i,j] > 1e-100
                A[i,j] = P[i,j] * exp(cumθ_y[i] + cumθ_x[j])
                if A[i,j] <= 1e-100
                    A[i,j] = 0.0
                end
                s += A[i,j]
            end
        end
        log(s) - η_target'θ
    end

    function g!(grad, θ)
        calcMarginals!(η, A)
        η ./= η[1]
        for i in 1:length(β)
            grad[i] = η[i+1] - η_target[i]
        end
    end

    function cb(x)
        local res = calcResidual(A, xlen/sum(A))
        return  res < ϵ
    end

    result = 0
    result = optimize((f, g!),
            zeros(length(β)), # initial θ
            if limited
                LBFGS(linesearch=Lnsrch(0.5, .0, 20, f), m=100)
                #LBFGS(linesearch=LineSearches.BackTracking(), m=100)
            else
                BFGS(linesearch=Lnsrch(0.5, .0, 10, f))
                #BFGS(linesearch=LineSearches.BackTracking())
            end,
            Optim.Options(
                    x_tol=NaN,
                    f_tol=NaN,
                    g_tol=NaN,
                    iterations=max_iter,
                    show_trace=show_trace,
                    show_every=show_every,
                    allow_f_increases=true,
                    callback=cb,
            ))
    return A./sum(A)
end

ext: jl


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

"""
    AbstractCoverer

A method for covering spatial objects (e.g. convex hull).
"""
abstract type AbstractCoverer end

"""
    cover(object, coverer)

Cover spatial `object` with `coverer` method.
"""
cover(object::AbstractSpatialObject, coverer::AbstractCoverer) =
  cover(domain(object), coverer)

"""
    cover(domain, coverer)

Cover spatial `domain` with `coverer` method.
"""
cover(domain::AbstractDomain, coverer::AbstractCoverer) = @error "not implemented"

#------------------
# IMPLEMENTATIONS
#------------------
include("covering/rectangle.jl")


text:
# ADMM

# made this to compare convergence with douglas rachford

# to be tested

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

#Base.@kwdef struct AdmmIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},F<:Function,G<:Function}
    #proxf!::F
#    f::F
#    g::G
#    x0::Tx #initial point
#    gamma::R #don't think gamma is needed
#end

Base.@kwdef struct AdmmIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},Tf,Tg}
    f::Tf = Zero()
    g::Tg = Zero()
    x0::Tx
    gamma::R
end

Base.IteratorSize(::Type{<:AdmmIteration}) = Base.IsInfinite()

Base.@kwdef struct AdmmState{Tx}
    x::Tx
    y::Tx = similar(x) #empty
    u::Tx = similar(x)
    dualres::Tx = similar(x)
end

function Base.iterate(iter::AdmmIteration, state::AdmmState = AdmmState(x=copy(iter.x0),y=copy(iter.x0),u=copy(iter.x0)),dualres=copy(iter.x0))
    prox!(state.x,iter.f,(state.y .- state.u), iter.gamma)
    prox!(state.y,iter.g, (state.x .+ state.u), iter.gamma)
    state.dualres .= state.x .- state.y
    state.u .+= state.dualres
    #iter.proxf!(state.xhat,state.uhat) # xhat^{k+1} = A(uhat^{k+1})
    #state.rhat .= 2 .*(state.xhat .- state.uhat)
    #iter.denoiser!(state.yhat,state.rhat) # yhat^{k+1} = B(2xhat^{k+1} - uhat^k)
    #state.res .= state.xhat .- state.yhat
    #state.uhat .+= (state.xhat - state.yhat) #uhat^{k+1} = uhat^k + (xhat^{k+1} - yhat^{k+1})
    return state, state
end


text:
export softmax_cross_entropy

# TODO: error on `scalar .+ vector`
# TODO: allow docstring
@i function softmax_cross_entropy(x, p, imax, xmax, Z, out::T) where T
    @anc logZ = zero(T)
    @anc yi = zero(T)
    # subtract maximum
    imax += argmax(x)  # trade off space of xmax to time
    xmax ⊕ x[imax]
    # accumulate exp(x) to Z, and finally get logZ
    for i=1:length(x)
        x[i] ⊖ xmax
        Z += exp(x[i])
    end
    logZ += log(Z)
    for i=1:length(x)
        yi ⊕ logZ
        yi ⊖ x[i]
        out += yi * p[i]
        yi ⊕ x[i]
        yi ⊖ logZ
    end
    logZ -= log(Z)
end

function _sce(x::AbstractArray{T,N}, p) where {T,N}
    x = x .- maximum(x; dims=N)  # avoid data overflow
    rho = exp.(x)
    Z = sum(rho; dims=N)
    return dropdims(sum((log.(Z) .- x) .* p; dims=N), dims=N)
end


text:
module EMLise
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("ise/ISE.jl")
end


text:
abstract type AbstractGate end
Base.adjoint(::AbstractGate) = error("not implemented")
matrix(::AbstractGate) = error("not implemented")

struct HGate <: AbstractGate end
const HMatrix = ComplexF64[1 1; 1 -1] ./ convert(ComplexF64, sqrt(2))
matrix(::HGate) = HMatrix

struct XGate <: AbstractGate end
const XMatrix = ComplexF64[0 1; 1 0]
matrix(::XGate) = XMatrix
X = XGate()

struct YGate <: AbstractGate end
const YMatrix = ComplexF64[0 -im; im 0]
matrix(::YGate) = YMatrix
Y = YGate()

struct ZGate <: AbstractGate end
const ZMatrix = ComplexF64[1 0; 0 1]
matrix(::ZGate) = ZMatrix
Z = ZGate()

struct SGate <: AbstractGate end
const SMatrix = ComplexF64[1 0; 0 im]
matrix(::SGate) = SMatrix
S = SGate()

struct SdagGate <: AbstractGate end
const SdagMatrix = ComplexF64[1 0; 0 -im]
matrix(::SdagGate) = SdagMatrix
Sdag = SdagGate()

struct TGate <: AbstractGate end
const TMatrix = ComplexF64[1 0; 0 Base.exp(im * π / 4)]
matrix(::TGate) = TMatrix
T = TGate()

struct TdagGate <: AbstractGate end
const TdagMatrix = ComplexF64[1 0; 0 Base.exp(-im * π / 4)]
matrix(::TdagGate) = TdagMatrix
Tdag = TGate()

struct RZGate
    θ::Vector{Float64}

    function RZGate(θ::Real)
        new([θ])
    end
end

function matrix(g::RZGate)
    m = zeros(ComplexF64, (2, 2))
    eθ = exp(im * g.θ[] / 2)
    m[1, 1] = conj(eθ)
    m[2, 2] = eθ
    m
end

max_stars_repo_name: isolatedinformation/JuliaGridSynth

max_issues_repo_name: isolatedinformation/JuliaGridSynth

max_forks_repo_name: isolatedinformation/JuliaGridSynth


text:
using FastTransforms, FillArrays, Test
import FastTransforms: clenshaw, clenshaw!, forwardrecurrence!, forwardrecurrence

@testset "clenshaw" begin
    @testset "Chebyshev T" begin
        for elty in (Float64, Float32)
            c = [1,2,3]
            cf = elty.(c)
            @test @inferred(clenshaw(c,1)) ≡ 1 + 2 + 3
            @test @inferred(clenshaw(c,0)) ≡ 1 + 0 - 3
            @test @inferred(clenshaw(c,0.1)) == 1 + 2*0.1 + 3*cos(2acos(0.1))
            @test @inferred(clenshaw(c,[-1,0,1])) == clenshaw!(c,[-1,0,1]) == [2,-2,6]
            @test clenshaw(c,[-1,0,1]) isa Vector{Int}
            @test @inferred(clenshaw(elty[],1)) ≡ zero(elty)

            x = elty[1,0,0.1]
            @test @inferred(clenshaw(c,x)) ≈ @inferred(clenshaw!(c,copy(x))) ≈
                @inferred(clenshaw!(c,x,similar(x))) ≈
                @inferred(clenshaw(cf,x)) ≈ @inferred(clenshaw!(cf,copy(x))) ≈
                @inferred(clenshaw!(cf,x,similar(x))) ≈ elty[6,-2,-1.74]

            @testset "Strided" begin
                cv = view(cf,:)
                xv = view(x,:)
                @test clenshaw!(cv, xv, similar(xv)) == clenshaw!(cf,x,similar(x))

                cv2 = view(cf,1:2:3)
                @test clenshaw!(cv2, xv, similar(xv)) == clenshaw([1,3], x)

                # modifies x and xv
                @test clenshaw!(cv2, xv) == xv == x == clenshaw([1,3], elty[1,0,0.1])
            end

            @testset "matrix coefficients" begin
                c = [1 2; 3 4; 5 6]
                @test clenshaw(c,0.1) ≈ [clenshaw(c[:,1],0.1), clenshaw(c[:,2],0.1)]
                @test clenshaw(c,[0.1,0.2]) ≈ [clenshaw(c[:,1], 0.1) clenshaw(c[:,2], 0.1); clenshaw(c[:,1], 0.2) clenshaw(c[:,2], 0.2)]
            end
        end
    end

    @testset "Chebyshev U" begin
        N = 5
        A, B, C = Fill(2,N-1), Zeros{Int}(N-1), Ones{Int}(N)
        @testset "forwardrecurrence!" begin
            @test @inferred(forwardrecurrence(N, A, B, C, 1)) == @inferred(forwardrecurrence!(Vector{Int}(undef,N), A, B, C, 1)) == 1:N
            @test forwardrecurrence!(Vector{Int}(undef,N), A, B, C, -1) == (-1) .^ (0:N-1) .* (1:N)
            @test forwardrecurrence(N, A, B, C, 0.1) ≈ forwardrecurrence!(Vector{Float64}(undef,N), A, B, C, 0.1) ≈
                    sin.((1:N) .* acos(0.1)) ./ sqrt(1-0.1^2)
        end

        c = [1,2,3]
        @test c'forwardrecurrence(3, A, B, C, 0.1) ≈ clenshaw([1,2,3], A, B, C, 0.1) ≈
            1 + (2sin(2acos(0.1)) + 3sin(3acos(0.1)))/sqrt(1-0.1^2)

        @testset "matrix coefficients" begin
            c = [1 2; 3 4; 5 6]
            @test clenshaw(c,A,B,C,0.1) ≈ [clenshaw(c[:,1],A,B,C,0.1), clenshaw(c[:,2],A,B,C,0.1)]
            @test clenshaw(c,A,B,C,[0.1,0.2]) ≈ [clenshaw(c[:,1], A,B,C,0.1) clenshaw(c[:,2], A,B,C,0.1); clenshaw(c[:,1], A,B,C,0.2) clenshaw(c[:,2], A,B,C,0.2)]
        end
    end

    @testset "Chebyshev-as-general" begin
        @testset "forwardrecurrence!" begin
            N = 5
            A, B, C = [1; fill(2,N-2)], fill(0,N-1), fill(1,N)
            Af, Bf, Cf = float(A), float(B), float(C)
            @test forwardrecurrence(N, A, B, C, 1) == forwardrecurrence!(Vector{Int}(undef,N), A, B, C, 1) == ones(Int,N)
            @test forwardrecurrence!(Vector{Int}(undef,N), A, B, C, -1) == (-1) .^ (0:N-1)
            @test forwardrecurrence(N, A, B, C, 0.1) ≈ forwardrecurrence!(Vector{Float64}(undef,N), A, B, C, 0.1) ≈ cos.((0:N-1) .* acos(0.1))
        end

        c, A, B, C = [1,2,3], [1,2,2], fill(0,3), fill(1,4)
        cf, Af, Bf, Cf = float(c), float(A), float(B), float(C)
        @test @inferred(clenshaw(c, A, B, C, 1)) ≡ 6
        @test @inferred(clenshaw(c, A, B, C, 0.1)) ≡ -1.74
        @test @inferred(clenshaw([1,2,3], A, B, C, [-1,0,1])) == clenshaw!([1,2,3],A, B, C, [-1,0,1]) == [2,-2,6]
        @test clenshaw(c, A, B, C, [-1,0,1]) isa Vector{Int}
        @test @inferred(clenshaw(Float64[], A, B, C, 1)) ≡ 0.0

        x = [1,0,0.1]
        @test @inferred(clenshaw(c, A, B, C, x)) ≈ @inferred(clenshaw!(c, A, B, C, copy(x))) ≈
            @inferred(clenshaw!(c, A, B, C, x, one.(x), similar(x))) ≈
            @inferred(clenshaw!(cf, Af, Bf, Cf, x, one.(x),similar(x))) ≈
            @inferred(clenshaw([1.,2,3], A, B, C, x)) ≈
            @inferred(clenshaw!([1.,2,3], A, B, C, copy(x))) ≈ [6,-2,-1.74]
    end

    @testset "Legendre" begin
        @testset "Float64" begin
            N = 5
            n = 0:N-1
            A = (2n .+ 1) ./ (n .+ 1)
            B = zeros(N)
            C = n ./ (n .+ 1)
            v_1 = forwardrecurrence(N, A, B, C, 1)
            v_f = forwardrecurrence(N, A, B, C, 0.1)
            @test v_1 ≈ ones(N)
            @test forwardrecurrence(N, A, B, C, -1) ≈ (-1) .^ (0:N-1)
            @test v_f ≈ [1,0.1,-0.485,-0.1475,0.3379375]

            n = 0:N # need extra entry for C in Clenshaw
            C = n ./ (n .+ 1)
            for j = 1:N
                c = [zeros(j-1); 1]
                @test clenshaw(c, A, B, C, 1) ≈ v_1[j] # Julia code
                @test clenshaw(c, A, B, C, 0.1) ≈  v_f[j] # Julia code
                @test clenshaw!(c, A, B, C, [1.0,0.1], [1.0,1.0], [0.0,0.0])  ≈ [v_1[j],v_f[j]] # libfasttransforms
            end
        end

        @testset "BigFloat" begin
            N = 5
            n = BigFloat(0):N-1
            A = (2n .+ 1) ./ (n .+ 1)
            B = zeros(N)
            C = n ./ (n .+ 1)
            @test forwardrecurrence(N, A, B, C, parse(BigFloat,"0.1")) ≈ [1,big"0.1",big"-0.485",big"-0.1475",big"0.3379375"]
        end
    end

    @testset "Int" begin
        N = 10; A = 1:10; B = 2:11; C = range(3; step=2, length=N+1)
        v_i = forwardrecurrence(N, A, B, C, 1)
        v_f = forwardrecurrence(N, A, B, C, 0.1)
        @test v_i isa Vector{Int}
        @test v_f isa Vector{Float64}

        j = 3
        @test clenshaw([zeros(Int,j-1); 1; zeros(Int,N-j)], A, B, C, 1) == v_i[j]
    end

    @testset "Zeros diagonal" begin
        N = 10; A = randn(N); B = Zeros{Int}(N); C = randn(N+1)
        @test forwardrecurrence(N, A, B, C, 0.1) == forwardrecurrence(N, A, Vector(B), C, 0.1)
        c = randn(N)
        @test clenshaw(c, A, B, C, 0.1) == clenshaw(c, A, Vector(B), C, 0.1)
    end
end


text:
module AIS

import DataFrames: DataFrame, select!
import CSV
import Geodesy: LLA
import ..ARLToolkit: dts

function read(filename; aoi=nothing, from=nothing, to=nothing, duration=nothing)
  to !== nothing && duration !== nothing && throw(ArgumentError("Both 'to' and 'duration' cannot be specificed"))
  to === nothing && duration !== nothing && (to = from + duration)
  df = CSV.read(filename, DataFrame; header=["yyyy", "mm", "dd", "HH", "MM", "SS", "mmsi", "lat", "lon", "sog", "cog", "heading"])
  df.SS .= min.(59, df.SS)
  df.time = dts.(df.yyyy, df.mm, df.dd, df.HH, df.MM, df.SS)
  df.location = LLA.(df.lat, df.lon)
  df = select!(df, [:time, :mmsi, :location])
  aoi === nothing || (df = filter(row -> row.location ∈ aoi, df))
  from === nothing || (df = filter(row -> row.time ≥ from, df))
  to === nothing || (df = filter(row -> row.time ≤ to, df))
  df
end

end # module


text:
"""
    RegularSpacing{T<:Real} <: AbstractVector{T}

`RegularSpacing(t0, Δt, N)` represents the same thing as `range(t0; step=Δt, length=N)`, but
has a different implementation that makes it possible to differentiate through with the
current version of `Zygote`. This data structure will be entirely removed once it's possible
to work with `StepRangeLen`s in `Zygote`.

Relevant issue: https://github.com/FluxML/Zygote.jl/issues/550
"""
struct RegularSpacing{T<:Real} <: AbstractVector{T}
    t0::T
    Δt::T
    N::Int
end

# Implements the AbstractArray interface.

Base.IndexStyle(::RegularSpacing) = Base.IndexLinear()

Base.size(x::RegularSpacing) = (x.N,)

Base.getindex(x::RegularSpacing, n::Int) = x.t0 + (n - 1) * x.Δt

Base.step(x::RegularSpacing) = x.Δt

ZygoteRules.@adjoint function (::Type{TR})(t0::T, Δt::T, N::Int) where {TR<:RegularSpacing, T<:Real}
    function pullback_RegularSpacing(Δ::TΔ) where {TΔ<:NamedTuple}
        return (
            hasfield(TΔ, :t0) ? Δ.t0 : nothing,
            hasfield(TΔ, :Δt) ? Δ.Δt : nothing,
            nothing,
        )
    end
    return RegularSpacing(t0, Δt, N), pullback_RegularSpacing
end


text:
using Pkg
Pkg.activate("./")
# load pkgs
using Ipopt, HarmonicPowerModels, PowerModels
using JuMP #adding avoids problems with Revise

# pkg const
const _PMs = PowerModels
const _HPM = HarmonicPowerModels

# set the solver
solver = Ipopt.Optimizer

# path to the data
path = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_2bus.m")

# load data
data  = _PMs.parse_file(path)
# sanity check
pm_fundamental = _PMs.instantiate_model(data, _PMs.IVRPowerModel, _PMs.build_opf_iv);
result_fundamental = optimize_model!(pm_fundamental, optimizer=solver)
result_fundamental["termination_status"]

# set up current ratings consistently
#assume minimum voltage
vmmin = 0.8
for (b,branch) in data["branch"]
    branch["c_rating"] = branch["rate_a"]/vmmin
end
for (d,load) in data["load"]
    load["c_rating"] = abs(load["pd"] + im* load["qd"])/vmmin
end
for (g,gen) in data["gen"]
    gen["c_rating"] = abs(gen["pmax"] + im* gen["qmax"])/vmmin
end

hdata = _HPM.replicate(data)
# hdata = _HPM.replicate(data, xfmr_exc=xfmr)

##
# solve the hopf
# result = run_hopf_iv(hdata, _PMs.IVRPowerModel, solver)
pm = _PMs.instantiate_model(hdata, _PMs.IVRPowerModel, _HPM.build_hopf_iv; ref_extensions=[_HPM.ref_add_xfmr!]);
result = optimize_model!(pm, optimizer=solver, solution_processors=[ _HPM.sol_data_model!])
_HPM.append_indicators!(result, hdata)
##

println("Harmonic 3")
_PMs.print_summary(result["solution"]["nw"]["2"])
println("Harmonic 1")
_PMs.print_summary(result["solution"]["nw"]["1"])
result["objective"]
result["termination_status"]


text:
using Flux
using BSON
using LinearAlgebra
using Statistics
using Random
using Distributions: Uniform
using Zygote: @adjoint

Random.seed!(1)

include("ch_ctc.jl")

const TRAINDIR = "train"

const EPOCHS = 300
const BATCH_SIZE = 1

const N_FILES = 400

losses = []

myinit(x) = Float32.(rand(Uniform(-0.1, 0.1), x))
myinit(x, y) = Float32.(rand(Uniform(-0.1, 0.1), x, y))

forward = LSTM(39, 200, init=myinit)
output = Dense(200, 62, initW=myinit)

m(x) = output.(forward.(x))

function loss(x, y)
  Flux.reset!(forward)
  yhat = m(x)
  yhat = reduce(hcat, yhat)
  l = c_ctc(yhat, y)
  addToGlobalLoss(l)
  return l
end

function readData(dataDir)
  fnames = open(readlines, "shuffled_names.txt")
  fnames = [x * ".bson" for x in fnames]
  fnames = fnames[1:N_FILES]
  Xs = []
  Ys = []

  for fname in fnames
    BSON.@load joinpath(dataDir, fname) x y
    x = [Float32.(x[i,:]) for i in 1:size(x,1)]
    push!(Xs, x)
    push!(Ys, y)
  end

  return (Xs, Ys)
end

"""
    lev(s, t)

Calculates the Levenshtein distance between `s` and `t`
"""
function lev(s, t)
  m = length(s)
  n = length(t)
  d = Array{Int}(zeros(m+1, n+1))

  for i=2:(m+1)
    @inbounds d[i, 1] = i-1
  end

  for j=2:(n+1)
    @inbounds d[1, j] = j-1
  end

  for j=2:(n+1)
    for i=2:(m+1)
      @inbounds if s[i-1] == t[j-1]
        substitutionCost = 0
      else
        substitutionCost = 1
      end
      @inbounds d[i, j] = min(d[i-1, j] + 1, # Deletion
                              d[i, j-1] + 1, # Insertion
                              d[i-1, j-1] + substitutionCost) # Substitution
    end
  end
  @inbounds return d[m+1, n+1]
end

"""
    collapse(seq)

Collapses `seq` into a version that has symbols collapsed into one repetition and removes all instances
of the blank symbol.
"""
function collapse(seq)
  s = [x for x in seq if x != 62]
  if isempty(s) return s end
  s = [seq[1]]
  for ch in seq[2:end]
    if ch != s[end] && ch != 62
      push!(s, ch)
    end
  end
  return s
end

"""
    per(x, y)

Compute the phoneme error rate of the model for input `x` and target `y`. The phoneme error rate
is defined as the Levenshtein distance between the labeling produced by running `x` through
the model and the target labeling in `y`, all divided by the length of the target labeling
in `y`
"""
function per(x, y)
  Flux.reset!(forward)
  yhat = m(x)
  yhat = reduce(hcat, yhat)
  yhat = mapslices(argmax, yhat, dims=1) |> vec |> collapse
  y = mapslices(argmax, y, dims=1) |> vec |> collapse
  return lev(yhat, y) / length(y)
end

function addToGlobalLoss(x)
	global losses
	push!(losses, x)
end

@adjoint function addToGlobalLoss(x)
	addToGlobalLoss(x)
	return nothing, () -> nothing
end

function main()
  println("Loading files")
  Xs, Ys = readData(TRAINDIR)
  data = collect(zip(Xs, Ys))

  println("Beginning training")

  opt = Momentum(1e-4)

  for i in 1:EPOCHS
    global losses
    losses = []
    println("Beginning epoch $i/$EPOCHS")
    Flux.train!(loss, params((forward, output)), data, opt)
    println("Calculating PER...")
    p = mean(map(x -> per(x...), data))
    println("PER: $(p*100)")
    println("Mean loss: ", mean(losses))
	  if p < 0.35 exit() end
  end
end

main()

max_stars_repo_name: maetshju/flux-ctc-grad

max_issues_repo_name: maetshju/flux-ctc-grad

max_forks_repo_name: maetshju/flux-ctc-grad


text:
#-------------canny.jl---------------------------------------------------------#
#
# Purpose: To implement Canny Edge Detection in Julia
#
#------------------------------------------------------------------------------#

# First, we need to use the appropriate packages
using Images  # All filtering and image stuff
using ImageView # imshow for showing images, must be run in REPL
using TestImages # Standart example images in julia

# simple implementation of canny edge detection using in-built Julia functions
function simple_edge_detection()
    img = testimage("fabio_color_256")
    img_edge = canny(img)
    save(string("fabio_edges.png"), img_edge)
end

simple_edge_detection()

max_stars_repo_name: vchuravy/skillpill-julia

max_issues_repo_name: oist/skillpill-julia

max_forks_repo_name: oist/skillpill-julia


text:
"""
All sampling functions have three versions:
1. new* : function to create a new variable, return params: 1) the new variable (through return value or reference) 2) the proposal probability
2. del* : function to remove the old variable, return the proposal probability
3. shift*  : shift the old variable to a new one, return params: 1) the new variable (through return value or reference) 2) the proposal probability
"""

# Before using the sampling functions, call setRand function to reset the random number generator
# function setRand(_rand)
#     global const rand=_rand
# end
import Random:GLOBAL_RNG

@inline newIdx(size, rng=GLOBAL_RNG) = rand(rng, 1:size), Float(size)
@inline delIdx(oldIdx, size) = Float(1.0) / size
@inline shiftIdx(oldIdx, size, rng=GLOBAL_RNG) = rand(rng, 1:size), Float(1.0)

@inline newTau(β, rng=GLOBAL_RNG) = rand(rng, Float) * β, β
@inline delTau(oldTau, β) = Float(1.0) / β
@inline function shiftTau(oldTau, β, TauScale=β/4.0, rng=GLOBAL_RNG)
    x = rand(rng, Float)
    newTau::Float = 0.0
    if x < 1.0 / 3
        newTau = oldTau + TauScale * (rand(rng, Float) - Float(0.5))
    elseif x < 2.0 / 3
        newTau = -oldTau
    else
        newTau = rand(rng, Float) * β
    end

    if newTau < Float(0.0)
        newTau += β
    elseif newTau > β
        newTau -= β
    end
    return newTau, Float(1.0)
end

@inline function newK!(newK, Kf, Kscale=Kf/2.0, rng=GLOBAL_RNG)::Float
    ############ Simple Way ########################
    # for i in 1:DIM
    #     newK[i] = Kf * (rand(rng) - 0.5) * 2.0
    # end
    # return (2.0 * Kf)^DIM
    ################################################

    Kamp = Kf + (2.0*rand(rng) - 1.0) * Kscale
    (Kamp <= 0.0) && (return 0.0)
    # Kf-dK<Kamp<Kf+dK
    ϕ = 2π * rand(rng, Float)
    if length(newK) == 3
        θ = π * rand(rng, Float)
        newK[1] = Kamp * cos(ϕ) * sin(θ)
        newK[2] = Kamp * sin(ϕ) * sin(θ)
        newK[3] = Kamp * cos(θ)
        return 2Kscale * 2π * π * (sin(θ) * Kamp^2)
        # prop density of Kamp in [Kf-dK, Kf+dK), prop density of ϕ
        # prop density of θ, Jacobian
    elseif length(newK)==2
        newK[1] = Kamp * cos(ϕ)
        newK[2] = Kamp * sin(ϕ)
        return 2Kscale * 2π * Kamp
        # prop density of Kamp in [Kf-dK, Kf+dK), prop density of ϕ, Jacobian
    else
        throw(DimensionMismatch("Dimension $(length(newK)) is not supported!"))
    end
end

@inline function delK(oldK, Kf, Kscale=Kf/2.0)::Float
    ############## Simple Way #########################
    # for i in 1:DIM
    #     if abs(oldK[i]) > Kf
    #         return 0.0
    #     end
    # end
    # return 1.0 / (2.0 * Kf)^DIM
    ####################################################

    Kamp = norm(oldK)
    if !(Kf - Kscale < Kamp < Kf + Kscale)
        return 0.0
    end
    # (Kamp < Kf - dK || Kamp > Kf + dK) && return 0.0
    if length(oldK) == 3
        sinθ = sqrt(oldK[1]^2 + oldK[2]^2) / Kamp
        (sinθ < eps(Float)) && (return 0.0)
        return 1.0 / (2Kscale * 2π * π * sinθ * Kamp^2)
    elseif length(oldK)==2
        return 1.0 / (2Kscale * 2π * Kamp)
    else
        throw(DimensionMismatch("Dimension $(length(newK)) is not supported!"))
    end
end

@inline function shiftK!(oldK, newK, Kf, Kscale=Kf/3.0, λ=1.5, rng=GLOBAL_RNG)::Float
    x = rand(rng)
    D=length(newK)
    if x < 1.0 / 3
        # dK = β > 1.0 ? Kf / β * 3.0 : Kf
        # newK .= oldK .+ (rand(rng, DIM) .- 0.5) .* dK
        for i in 1:D
            newK[i] = oldK[i] + (rand(rng) - 0.5) * Kscale
        end
        return 1.0
    elseif x < 2.0 / 3
        ratio = 1.0 / λ + rand(rng) * (λ - 1.0 / λ)
        newK .= oldK .* ratio
        if D==3
            return ratio
        elseif D==2
            return 1.0
        else
            throw(DimensionMismatch("Dimension $(length(newK)) is not supported!"))
        end
    else
        newK .= oldK .* (-1.0)
        return 1.0
    end
end

max_stars_repo_name: kunyuan/DiagramCalculator

max_issues_repo_name: kunyuan/DiagramCalculator

max_forks_repo_name: kunyuan/DiagramCalculator


text:
using NamedArrays, BenchmarkTools, Test

function f1()
    x = Vector{Int}(10^6)
    y = NamedArray(x)
    for i in 1:10^6, j in 1:10
        x[i] = i
    end
end

function f2()
    x = Vector{Int}(10^6)
    y = NamedArray(x)
    for i in 1:10^6, j in 1:10
        y[i] = i
    end
end

function f3()
    x = Vector{Int}(10^6)
    y = NamedArray(x, (NamedArrays.defaultnamesdict(length(x)),), (:A,))
    for i in 1:10^6, j in 1:10
        y[i] = i
    end
end

println("@benchmark f1()")


text:
#RSDeltaSigma: State space utilities
#-------------------------------------------------------------------------------
#SEARCH: even, odd

#==getorder
===============================================================================#
#Figure out nu, nq, and order of the system & ensure dimensions are ok
function get_nu_nq_order(u, ABCD, nlev)
	u = conv2seriesmatrix2D(u)
	nu = size(u,1)
	nq = length(nlev)
	if !(size(ABCD,2) > 2 && size(ABCD,2)==nu+size(ABCD,1))
		msg = "The ABCD argument does not have proper dimensions."
		throw(ArgumentError(msg))
	end
	order = size(ABCD,1)-nq
	return (nu, nq, order)
end

#==partitionABCD
===============================================================================#
"""`[A B C D] = partitionABCD(ABCD, m)`

Partition ABCD into A, B, C, D for an m-input state-space system.
"""
function partitionABCD(ABCD, m=nothing)
	n = 0
	if m!= nothing
		n = size(ABCD,2) - m
	else
		n = minimum(size(ABCD)) - 1
		m = size(ABCD,2) - n
	end
	r = size(ABCD,1) - n

	A = ABCD[1:n, 1:n]
	B = ABCD[1:n, n+1:n+m]
	C = ABCD[n+1:n+r, 1:n]
	D = ABCD[n+1:n+r, n+1:n+m]

	return (A, B, C, D)
end

#==stuffABCD
===============================================================================#
function _stuffABCD(::DS{:CRFB}, a,g,b,c)
	order = length(a)
	ABCD = zeros(order+1,order+2)
	(order2, isodd, Δodd) = orderinfo(order)
	Δeven = 1-Δodd

	#C=(0 0...c_n)
	#This is done as part of the construction of A, below
	#B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
	ABCD[:,order+1] = b'
	#B2 = -(a_1 a_2... a_n)
	ABCD[1:order,order+2] = -a'
	diagonal = 1:(order+2):order*(order+1)
	ABCD[diagonal] = ones(1,order)
	subdiag = diagonal[(1+Δeven):2:order] .+ 1
	ABCD[subdiag] = c[1+Δeven:2:order]
	supdiag = subdiag[(1+Δodd):length(subdiag)] .- 2
	ABCD[supdiag] = -g
	dly = (2+Δodd):2:order #row numbers of delaying integrators
	ABCD[dly,:] = ABCD[dly,:] .+ diagm(c[dly .- 1]) * ABCD[dly .- 1,:]

	return ABCD
end

#=
switch form
    case 'CRFF'
        %B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD(:,order+1) = b';
        %B2 = -(c_1 0... 0)
        ABCD(1,order+2) = -c(1);
        diagonal = 1:(order+2):order*(order+1);	% # of elements = order
        ABCD(diagonal) = ones(1,order);
        subdiag = diagonal(1:2:order-1)+1;
        ABCD(subdiag)= c(2:2:order);
        if even
            multg = 1:2:order;	% rows to have g*(following row) subtracted.
            ABCD(multg,:) = ABCD(multg,:) - diagm(g)*ABCD(multg+1,:);
        elseif order >= 3
            supdiag = diagonal(3:2:order)-1;
            ABCD(supdiag) = -g;
        end
        multc = 3:2:order;		% rows to have c*(preceding row) added.
        ABCD(multc,:) = ABCD(multc,:) + diagm(c(multc))*ABCD(multc-1,:);
        ABCD(order+1,1:2:order) = a(1:2:order);
        for i = 2:2:order
            ABCD(order+1,:) = ABCD(order+1,:) + a(i)*ABCD(i,:);
        end

    case 'CIFB'
        %C=(0 0...c_n)
        % This is done as part of the construction of A, below
        %B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD(:,order+1) = b';
        %B2 = -(a_1 a_2... a_n)
        ABCD(1:order,order+2) = -a';
        diagonal = 1:(order+2):order*(order+1);
        ABCD(diagonal) = ones(1,order);
        subdiag = diagonal(1:order)+1;
        ABCD(subdiag)= c;
        supdiag = diagonal((2+odd):2:order)-1;
        ABCD(supdiag) = -g;

    case 'CIFF'
        %B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD(:,order+1) = b';
        %B2 = -(c_1 0... 0)
        ABCD(1,order+2) = -c(1);
        diagonal = 1:(order+2):order*(order+1);
        ABCD(diagonal) = ones(1,order);
        subdiag = diagonal(1:order-1)+1;
        ABCD(subdiag)= c(2:end);
        %C = (a_1 a_2... a_n)
        ABCD(order+1,1:order) = a(1:order);
        supdiag = diagonal((2+odd):2:order)-1;
        ABCD(supdiag) = -g;

    case 'CRFBD'
        %C=(0 0...c_n)
        ABCD(order+1,order) = c(order);
        %B1 = (b_1 b_2... b_n), D=(b_n+1 0)
        ABCD(:,order+1) = b';
        %B2 = -(a_1 a_2... a_n)
        ABCD(1:order,order+2) = -a';
        diagonal = 1:(order+2):order*(order+1);
        ABCD(diagonal) = ones(1,order);
        dly = (1+odd):2:order;	% row numbers of delaying integrators
        subdiag = diagonal(dly)+1;
        ABCD(subdiag)= c(dly);
        supdiag = subdiag((1+odd):length(subdiag))-2;
        ABCD(dly,:) = ABCD(dly,:) - diagm(g)*ABCD(dly+1,:);
        if order>2
            coupl = 2+even:2:order;
            ABCD(coupl,:) = ABCD(coupl,:) + diagm(c(coupl-1))*ABCD(coupl-1,:);
        end

    case 'CRFFD'
        diagonal = 1:(order+2):order*(order+1);
        subdiag = diagonal(1:order-1)+1;
        supdiag = diagonal(2:order)-1;
        %B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD(:,order+1) = b';
        %B2 = -(c_1 0... 0)
        ABCD(1,order+2) = -c(1);
        ABCD(diagonal) = ones(1,order);
        multc = 2:2:order;		% rows to have c*(preceding row) added.
        if order>2
            ABCD(subdiag(2:2:end)) = c(3:2:end);
        end
        if even
            ABCD(supdiag(1:2:end)) = -g;
        else
            % subtract g*(following row) from the multc rows
            ABCD(multc,:) = ABCD(multc,:) - diagm(g)*ABCD(multc+1,:);
        end
        ABCD(multc,:) = ABCD(multc,:) + diagm(c(multc))*ABCD(multc-1,:);
        % C
        ABCD(order+1,2:2:order) = a(2:2:order);
        for i = 1:2:order
            ABCD(order+1,:) = ABCD(order+1,:) + a(i)*ABCD(i,:);
        end
        % The above gives y(n+1); need to add a delay to get y(n).
        % Do this by augmenting the states. Note: this means that
        % the apparent order of the NTF is one higher than it acually is.
        [A B C D] = partitionABCD(ABCD,2);
        A = [A zeros(order,1); C 0];
        B = [B; D];
        C = [zeros(1,order) 1];
        D = [0 0];
        ABCD = [A B; C D];

    case 'PFF'
        %B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD(:,order+1) = b';
        odd_1 = odd;		% !! Bold assumption !!
        odd_2 = 0;			% !! Bold assumption !!
        gc = g.* c(1+odd_1:2:end);
        theta = acos(1-gc/2);
        if odd_1
            theta0 = 0;
        else
            theta0 = theta(1);
        end
        order_2 = 2*length(find(abs(theta-theta0)>0.5));
        order_1 = order - order_2;
        %B2 = -(c_1 0...0 c_n 0...0)
        ABCD(1,order+2) = -c(1);
        ABCD(order_1+1,order+2) = -c(order_1+1);
        diagonal = 1:(order+2):order*(order+1);	% # of elements = order
        ABCD(diagonal) = ones(1,order);
        i = [1:2:order_1-1 order-order_2+1:2:order]
        subdiag = diagonal(i)+1;
        ABCD(subdiag)= c(i+1);
        if odd_1
            if order_1 >= 3
                supdiag = diagonal(3:2:order_1)-1;
                ABCD(supdiag) = -g(1:(order_1-1)/2);
            end
        else
            multg = 1:2:order_1;	% rows to have g*(following row) subtracted.
            ABCD(multg,:) = ABCD(multg,:) - diagm(g(1:order_1/2))*ABCD(multg+1,:);
        end
        if odd_2
            if order_2 >= 3
                supdiag = diagonal(order_1+2:2:order)-1;
                ABCD(supdiag) = -g(1:(order_1-1)/2);
            end
        else
            multg = order_1+1:2:order; % rows to have g*(following row) subtracted.
            gg = g((order_1-odd_1)/2+1:end);
            ABCD(multg,:) = ABCD(multg,:) - diagm(gg)*ABCD(multg+1,:);
        end
        % Rows to have c*(preceding row) added.
        multc = [3:2:order_1 order_1+3:2:order];
        ABCD(multc,:) = ABCD(multc,:) + diagm(c(multc))*ABCD(multc-1,:);
        % C portion of ABCD
        i = [1:2:order_1 order_1+1:2:order];
        ABCD(order+1,i) = a(i);
        for i = [2:2:order_1 order_1+2:2:order]
            ABCD(order+1,:) = ABCD(order+1,:) + a(i)*ABCD(i,:);
        end

    case 'Stratos'
        % code copied from case 'CIFF':
        %B1 = (b_1 b_2... b_n), D=(b_(n+1) 0)
        ABCD(:,order+1) = b';
        %B2 = -(c_1 0... 0)
        ABCD(1,order+2) = -c(1);
        diagonal = 1:(order+2):order*(order+1);
        ABCD(diagonal) = ones(1,order);
        subdiag = diagonal(1:order-1)+1;
        ABCD(subdiag)= c(2:end);
        % code based on case 'CRFF':
        multg = 1+odd:2:order-1;	% rows to have g*(following row) subtracted.
        ABCD(multg,:) = ABCD(multg,:) - diagm(g)*ABCD(multg+1,:);
        % code copied from case 'CIFF':
        %C = (a_1 a_2... a_n)
        ABCD(order+1,1:order) = a(1:order);

    case 'DSFB'
        I = eye(order);
        A0 = I;
        A1 = zeros(order,order);
        indices = 2:order+1:(order-1)*order;    % subdiagonal
        A1(indices) = c(1:order-1);
        indices = order+1 : 2*(order+1) : order*order-1; %supdiagonal
        A1(indices) = -g;
        C0 = I;
        C1 = I;
        D0 = zeros(order,2);
        if odd
            C1(order,order) = 0;
        else
            A1(order-1:order,order-1:order) = [ 0 -g(end)/2; c(end-1) 0];
            C0(order,order) = 0;
            D0(order-1,2) = g(end)/2;
        end
        M = inv(I - A1*C1);
        A = M * (A0 + A1*C0 );
        B = M * ( [b(1:order)' -a'] + A1*D0 );
        C = [zeros(1,order-1) c(order)];
        D = [b(order+1) 0];
        ABCD = [A B; C D];

    otherwise
        fprintf(1,'%s error: Form %s is not yet supported.\n',mfilename,form)
end
=#

"""`ABCD = stuffABCD(a,g,b,c,form='CRFB')`

Compute the ABCD matrix for the specified structure.
See realizeNTF.m for a list of supported structures.
`mapABCD` is the inverse function.
"""
function stuffABCD(a,g,b,c, form::Symbol=:CRFB)
	order = length(a)
	if length(b)==1
		b = [b zeros(1,order)]
	end

	return _stuffABCD(DS(form), a, g, b, c)
end

#==mapABCD
===============================================================================#
#Common implementation for CRFB, CIFB, & CRFBD:
function _mapABCD_CRFB_CIFB_CRFBD(ABCD, subdiag, supdiag, order::Int, form::Symbol)
	order2, isodd, Δodd = orderinfo(order)
	Δeven = 1 - Δodd
	c = ABCD[subdiag]
	g = -ABCD[supdiag]
	if :CRFB == form
		dly = (2+Δodd):2:order #row numbers of delaying integrators.
		ABCD[dly,:] = ABCD[dly,:] .- diagm(c[dly .- 1])*ABCD[dly .- 1,:]
	elseif :CRFBD == form
		dly = (1+Δodd):2:order #row numbers of delaying integrators.
		ABCD[dly,:] = ABCD[dly,:] .+ diagm(g)*ABCD[dly .+ 1,:]
		if order>2
			coupl = 2+Δeven:2:order
			ABCD[coupl,:] = ABCD[coupl,:] .- diagm(c[coupl .- 1])*ABCD[coupl .- 1,:]
		end
	end
	a = -ABCD[1:order,order+2]'
	b = ABCD[:,order+1]'
	return (a,g,b,c)
end

_mapABCD(::DS{:CRFB}, ABCD, diagonal, subdiag, supdiag, order::Int) =
	_mapABCD_CRFB_CIFB_CRFBD(ABCD, subdiag, supdiag, order, :CRFB)
_mapABCD(::DS{:CIFB}, ABCD, diagonal, subdiag, supdiag, order::Int) =
	_mapABCD_CRFB_CIFB_CRFBD(ABCD, subdiag, supdiag, order, :CIFB)
_mapABCD(::DS{:CRFBD}, ABCD, diagonal, subdiag, supdiag, order::Int) =
	_mapABCD_CRFB_CIFB_CRFBD(ABCD, subdiag, supdiag, order, :CRFBD)

#=
case 'CRFF'
    c = [-ABCD(1,order+2) ABCD(subdiag(1:end-1))];
    g = -ABCD(supdiag);
    if even
        multg = 1:2:order;	%Rows to have g*(following row) added.
        ABCD(multg,:) = ABCD(multg,:)+diag(g)*ABCD(multg+1,:);
    end
    multc=3:2:order;	%Rows to have c*(preceding row) subtracted.
    ABCD(multc,:)=ABCD(multc,:)-diag(c(multc))*ABCD(multc-1,:);
    a(2:2:order)=ABCD(order+1,2:2:order);
    for i=2:2:order                   %Recover a coeff.
        ABCD(order+1,:)=ABCD(order+1,:)-a(i)*ABCD(i,:);
    end
    a(1:2:order)=ABCD(order+1,1:2:order);
    b = ABCD(:,order+1)';

case 'CRFFD'
    % CRFFD has an extra order. Correct the common variables
    order = order-1;
	order2, isodd, Δodd = orderinfo(order)
    odd = rem(order,2); even = ~odd;
    diagonal = diagonal(1:order);
    subdiag = diagonal(1:order-1)+1;
    supdiag = diagonal(2+Δodd:2:order)-1;
    g = -ABCD(supdiag);
    c = [-ABCD(1,order+3) ABCD(subdiag)];
    a = zeros(1,order);
    for i=1:2:order                   %Recover the a coefficients
		a(i) = ABCD(order+1,i);
		ABCD(order+1,:) = ABCD(order+1,:) - a(i)*ABCD(i,:);
    end
    a(2:2:order) = ABCD(order+1,2:2:order);
    b = ABCD(1:order+1,order+2)';
    for i=2:2:order                   %Recover the b coefficients
		b(i) = b(i) - c(i) *b(i-1);
		if isodd
		    b(i) = b(i) + g(i/2) *b(i+1);
		end
    end
    yscale =  ABCD(order+2,order+1);
    a = a*yscale;
    b(end) = b(end)*yscale;

case {'CIFF','Stratos'}
    a = ABCD(order+1,1:order);
    c = [-ABCD(1,order+2) ABCD(subdiag(1:end-1))];
    g = -ABCD(supdiag);
    b = ABCD(:,order+1)';

otherwise
    fprintf(1,'%s error: Form %s is not yet supported.\n',mfilename,form)
end
=#

"""`(a,g,b,c) = mapABCD(ABCD, form=:CRFB)`

Compute the coefficients for the specified structure.
See `realizeNTF` for a list of supported structures.
`stuffABCD` is the inverse function.
"""
function mapABCD(ABCD, form::Symbol=:CRFB)
	order=size(ABCD,1)-1
	order2, isodd, Δodd = orderinfo(order)
	diagonal = 1:(order+2):order*(order+1)
	subdiag = diagonal[1:order] .+ 1
	supdiag = diagonal[2+Δodd:2:order] .- 1

	#Algorithm below modifies ABCDs copy so as not to modify original:
	ABCD = deepcopy(ABCD)

	#Do the mapping.
	return _mapABCD(DS(form), ABCD, diagonal, subdiag, supdiag, order)
end

#==scaleABCD
===============================================================================#
"""`(ABCDs,umax,S)=scaleABCD(ABCD; nlev=2,f=0,xlim=1,ymax=nlev+5,umax,N=10^5,N0=10)`

Scale the loop filter of a general delta-sigma modulator for dynamic range.

 - ABCD: The state-space description of the loop filter.
 - nlev: The number of levels in the quantizer.
 - xlim: A vector or scalar specifying the limit for each state variable.
 - ymax: The stability threshold. Inputs that yield quantizer inputs above ymax
   are considered to be beyond the stable range of the modulator.
 - umax: The maximum allowable input amplitude. umax is calculated if it
   is not supplied.
 - ABCDs: The state-space description of the scaled loop filter.
 - S: The diagonal scaling matrix, ABCDs = [S*A*Sinv S*B; C*Sinv D];
	xs = S*x;
"""
function scaleABCD(ABCD; nlev::Int=2, f=0, xlim=1, ymax=NaN, umax=NaN, N::Int=10^5, N0=10)
	N_sim=N
	order = size(ABCD,1)-1
	if isnan(ymax)
		ymax=nlev+5
	end
	if length(xlim)==1
		xlim = xlim * ones(1,order) #Convert scalar xlim to a vector
	end
	quadrature = !isreal(ABCD)

	#Make this function repeatable
	#rng_state = rng; rng('default') #Original code
	Random.seed!(0)

	#Envelope for smooth start-up
	raised_cosine = 0.5*(1 .- cos.(pi/N0 * (0:N0-1)))'
	if isnan(umax)
		#Simulate the modulator with DC or sine wave inputs to detect its stable
		#input range.
		#First get a rough estimate of umax.
		ulist =(0.1:0.1:1.0) * (nlev-1)
		umax = nlev-1
		N = 10^3
		u0 = [ exp.(2π*j*f*(-N0:-1))' .* raised_cosine exp.(2π*j*f*(0:N-1))' ] .+ 0.01*[1 j]*randn(2, N+N0)
		if !quadrature; u0 = real.(u0); end
		local simresult
		for u in ulist
			if !quadrature
				simresult = simulateDSM(u*u0,ABCD, nlev=nlev, trackmax=true)
			else
				simresult = simulateQDSM(u*u0,ABCD, nlev=nlev, trackmax=true)
			end
			if maximum(abs.(simresult.y)) > ymax
				umax = u #umax is the smallest input found which causes 'instability'
				break
			end
		end
		if umax == ulist[1]
			umaxstr = @sprintf("%.1f", umax)
			msg = "Modulator is unstable even with an input amplitude of $umaxstr."
			throw(msg)
		end
	end
	#More detailed simulation
	N = N_sim
	u0 = [ exp.(2π*j*f*(-N0:-1))' .* raised_cosine exp.(2π*j*f*(0:N-1))' ] + 0.01*[1 j]*randn(2, N+N0)
	if !quadrature; u0 = real.(u0); end
	maxima = zeros(1,order) .- 1
	ulist = range(0.7*umax, stop=umax, length=10)

	local simresult
	for u in ulist
		if !quadrature
			simresult = simulateDSM(u*u0,ABCD, nlev=nlev, trackmax=true)
		else
			simresult = simulateQDSM(u*u0,ABCD, nlev=nlev, trackmax=true)
		end
		if maximum(abs.(simresult.y)) > ymax
			break
		end
		#We need to do this at least once.
		umax = u #umax is the largest input which keeps |y| < ymax
		#maxima = maximum([maxima; simresult.xmax'], dims=1)
		maxima = max.(maxima, simresult.xmax')
	end

	#Scale the modulator so that all states are at most xlim.
	scale = maxima ./ xlim
	scale = scale[:] #Convert to simple vector
	S = diagm(1 ./ scale); Sinv = diagm(scale) #xs = S * x
	(A, B, C, D) = partitionABCD(ABCD)
	ABCDs = [S*A*Sinv S*B; C*Sinv D]

	Random.seed!() #original code restores state: rng(rng_state)
	return (ABCDs, umax, S)
end

#Last line


text:
using Makie
using Distributions: Normal, pdf

function main()
    μ = 0
    σ = 1
    d = Normal(μ, σ)
    xs = -3.:0.01:3.
    ys = pdf.(d, xs)
    lines(xs, ys)
end

function main()
    μ = 0
    σ = 1
    d = Normal(μ, σ)
    xs = -3.:0.01:3.
    pos = map(xs) do x
        Point2f0(x, pdf(d, x))
    end
    lines(pos)
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
    GeneralizedLinearMixedModel
Generalized linear mixed-effects model representation

Members:

- `LMM`: a [`LinearMixedModel`](@ref) - used for the random effects only.
- `dist`: a `UnivariateDistribution` - typically `Bernoulli()`, `Binomial()`, `Gamma()` or `Poisson()`.
- `link`: a suitable `GLM.Link` object
- `β`: the fixed-effects vector
- `θ`: covariance parameter vector
- `b`: similar to `u`, equivalent to `broadcast!(*, b, LMM.Λ, u)`
- `u`: a vector of matrices of random effects
- `u₀`: similar to `u`.  Used in the PIRLS algorithm if step-halving is necessary.
- `X`:
- `y`: the response vector
- `μ`: the mean vector
- `η`: the linear predictor
- `devresid`: vector of squared deviance residuals
- `offset`: offset₀ + `X * β`
- `offset₀`: prior offset; `T[]` is allowed
- `wrkresid`: vector of working residuals
- `wrkwt`: vector of working weights
- `wt`: vector of prior case weights, a value of `T[]` indicates equal weights.
"""

type GeneralizedLinearMixedModel{T <: AbstractFloat, D <: UnivariateDistribution, L <: Link} <: MixedModel
    LMM::LinearMixedModel{T}
    dist::D
    link::L
    β::Vector{T}
    θ::Vector{T}
    b::Vector{Matrix{T}}
    u::Vector{Matrix{T}}
    u₀::Vector{Matrix{T}}
    X::Matrix{T}
    y::Vector{T}
    μ::Vector{T}
    η::Vector{T}
    devresid::Vector{T}
    offset::Vector{T}
    offset₀::Vector{T}
    wrkresid::Vector{T}
    wrkwt::Vector{T}
    wt::Vector{T}
end

fixef(m::GeneralizedLinearMixedModel) = m.β

"""
    glmm(f::Formula, fr::ModelFrame, d::Distribution[, l::GLM.Link])

Return a `GeneralizedLinearMixedModel` object.

The value is ready to be `fit!` but has not yet been fit.
"""
function glmm(f::Formula, fr::AbstractDataFrame, d::Distribution, l::Link; wt=[], offset=[])
    if d == Binomial() && isempty(wt)
        d = Bernoulli()
    end
    wts = isempty(wt) ? ones(nrow(fr)) : Array(wt)
        # the weights argument is forced to be non-empty in the lmm as it will be used later
    LMM = lmm(f, fr; weights = wts)
    A, R, trms, u, y = LMM.A, LMM.R, LMM.trms, ranef(LMM), copy(model_response(LMM))
    wts = oftype(y, wts)
    kp1 = length(LMM.Λ) + 1
    X = trms[kp1]
            # zero the dimension of the fixed-effects in trms, A and R
    trms[kp1] = zeros(length(y), 0)
    LMM.wttrms[kp1] = trms[kp1]
    for i in 1:kp1
        qi = size(trms[i], 2)
        A[i, kp1] = zeros((qi, 0))
        R[i, kp1] = zeros((qi, 0))
    end
    qend = size(trms[end], 2)  # should always be 1 but no harm in extracting it
    A[kp1, end] = zeros((0, qend))
    R[kp1, end] = zeros((0, qend))
            # fit a glm to the fixed-effects only
    gl = glm(X, y, d, l; wts = wts)
    r = gl.rr
    res = GeneralizedLinearMixedModel(LMM, d, l, coef(gl), getθ(LMM), deepcopy(u), u, map(zeros, u),
        X, y, r.mu, r.eta, r.devresid, copy(r.eta), oftype(y, offset), r.wrkresid, r.wrkwts,
        oftype(y, wt))
    wrkresp!(trms[end], res)
    reweight!(LMM, res.wrkwt)
    res
end

glmm(f::Formula, fr::AbstractDataFrame, d::Distribution) = glmm(f, fr, d, GLM.canonicallink(d))

Base.logdet{T}(m::GeneralizedLinearMixedModel{T}) = logdet(m.LMM)

"""
    LaplaceDeviance{T,D}(m::GeneralizedLinearMixedModel{T,D})

Return the Laplace approximation to the deviance of `m`.

If the distribution `D` does not have a scale parameter the Laplace approximation
is defined as the squared length of the conditional modes, `u`, plus the determinant
of `Λ'Z'ZΛ + 1`, plus the sum of the squared deviance residuals.
"""
LaplaceDeviance{T}(m::GeneralizedLinearMixedModel{T}) =
    sum(m.devresid) + logdet(m) + mapreduce(sumabs2, +, m.u)

"""
    LaplaceDeviance!(m::GeneralizedLinearMixedModel)

Update `m.η`, `m.μ`, etc., install the working response and working weights in
`m.LMM`, update `m.LMM.A` and `m.LMM.R`, then evaluate `LaplaceDeviance`.
"""
function LaplaceDeviance!(m::GeneralizedLinearMixedModel)
    updateη!(m) |> updateμ!
    wrkresp!(m.LMM.trms[end], m)
    reweight!(m.LMM, m.wrkwt)
    LaplaceDeviance(m)
end

function StatsBase.loglikelihood{T,D}(m::GeneralizedLinearMixedModel{T,D})
    accum = zero(T)
    if D <: Binomial
        for (μ, y, n) in zip(m.μ, m.y, m.wt)
            accum += logpdf(D(round(Int, n), μ), round(Int, y * n))
        end
    else
        for (μ, y) in zip(m.μ, m.y)
            accum += logpdf(D(μ), y)
        end
    end
    accum - (mapreduce(sumabs2, + , m.u) + logdet(m)) / 2
end

lowerbd(m::GeneralizedLinearMixedModel) = vcat(fill(-Inf, size(m.β)), lowerbd(m.LMM))

function restoreX!(m::GeneralizedLinearMixedModel)
    if !isempty(m.LMM.R[end - 1, end - 1])
        return m
    end
    lm, X = m.LMM, m.X
    A, R, trms, k = lm.A, lm.R, lm.trms, length(lm.Λ)
    kp1 = k + 1
    trms[kp1] = X
    lm.wttrms[kp1] = copy(X)
    for i in 1 : kp1
        A[i, kp1] = trms[i]'X
        R[i, kp1] = copy(A[i, kp1])
    end
    A[kp1, end] = X'trms[end]
    R[kp1, end] = copy(A[kp1, end])
    reweight!(lm, m.wrkwt)
end

"""
    updateη!(m::GeneralizedLinearMixedModel)

Update the linear predictor, `m.η`, from the offset and the `B`-scale random effects.
"""
function updateη!(m::GeneralizedLinearMixedModel)
    η, lm, b, offset, u = m.η, m.LMM, m.b, m.offset, m.u
    Λ, trms = lm.Λ, lm.trms
    isempty(offset) ? fill!(η, 0) : copy!(η, offset)
    for i in eachindex(b)
        unscaledre!(η, trms[i], A_mul_B!(Λ[i], copy!(b[i], u[i])))
    end
    m
end

"""
    pirls!(m::GeneralizedLinearMixedModel)

Use Penalized Iteratively Reweighted Least Squares (PIRLS) to determine the conditional
modes of the random effects.
"""
function pirls!{T}(m::GeneralizedLinearMixedModel{T})
    iter, maxiter, obj = 0, 100, T(-Inf)
    u₀, u = m.u₀, m.u
    for j in eachindex(u)         # start from u all zeros
        copy!(u₀[j], fill!(u[j], 0))
    end
    obj₀ = LaplaceDeviance!(m)
    while iter < maxiter
        iter += 1
        ranef!(u, m.LMM, true)    # solve for new values of u
        obj = LaplaceDeviance!(m) # update GLM vecs and evaluate Laplace approx
        nhalf = 0
        while obj > obj₀
            nhalf += 1
            if nhalf > 10
                if iter < 2
                    throw(ErrorException("number of averaging steps > 10"))
                end
                break
            end
            for i in eachindex(u)
                ui = u[i]
                ui₀ = u₀[i]
                for j in eachindex(ui)
                    ui[j] += ui₀[j]
                    ui[j] *= 0.5
                end
            end
            obj = LaplaceDeviance!(m)
        end
        if isapprox(obj, obj₀; atol = 0.0001)
            break
        end
        for i in eachindex(u)
            copy!(u₀[i], u[i])
        end
        obj₀ = obj
    end
    obj
end

"""
    setβθ!{T}(m::GeneralizedLinearMixedModel{T}, v::Vector{T})

Set the parameter vector, `:βθ`, of `m` to `v`.

`βθ` is the concatenation of the fixed-effects, `β`, and the covariance parameter, `θ`.
"""
function setβθ!{T}(m::GeneralizedLinearMixedModel{T}, v::Vector{T})
    β, lm, offset, offset₀, X = m.β, m.LMM, m.offset, m.offset₀, m.X
    lb = length(β)
    copy!(β, Compat.view(v, 1 : lb))
    setθ!(m.LMM, copy!(m.θ, Compat.view(v, (lb + 1) : length(v))))
    BLAS.gemv!('N', one(T), X, β, one(T), isempty(offset₀) ? fill!(offset, 0) : copy!(offset, offset₀))
    m
end

sdest{T <: AbstractFloat}(m::GeneralizedLinearMixedModel{T}) = one(T)

"""
    fit!(m::GeneralizedLinearMixedModel[, verbose = false, optimizer=:LN_BOBYQA]])

Optimize the objective function for `m`
"""
function StatsBase.fit!(m::GeneralizedLinearMixedModel, verbose::Bool=false, optimizer::Symbol=:LN_BOBYQA)
    β, lm = m.β, m.LMM
    βθ = vcat(β, getθ(lm))
    opt = NLopt.Opt(optimizer, length(βθ))
    NLopt.ftol_rel!(opt, 1e-12)   # relative criterion on deviance
    NLopt.ftol_abs!(opt, 1e-8)    # absolute criterion on deviance
    NLopt.xtol_abs!(opt, 1e-10)   # criterion on parameter value changes
    NLopt.lower_bounds!(opt, vcat(fill!(similar(β), -Inf), lowerbd(lm)))
    feval = 0
    function obj(x::Vector{Float64}, g::Vector{Float64})
        if length(g) ≠ 0
            error("gradient not defined for this model")
        end
        feval += 1
        setβθ!(m, x) |> pirls!
    end
    if verbose
        function vobj(x::Vector{Float64}, g::Vector{Float64})
            if length(g) ≠ 0
                error("gradient not defined for this model")
            end
            feval += 1
            val = setβθ!(m, x) |> pirls!
            print("f_$feval: $(round(val,5)), [")
            showcompact(x[1])
            for i in 2:length(x) print(","); showcompact(x[i]) end
            println("]")
            val
        end
        NLopt.min_objective!(opt, vobj)
    else
        NLopt.min_objective!(opt, obj)
    end
    fmin, xmin, ret = NLopt.optimize(opt, βθ)
    ## very small parameter values often should be set to zero
#    xmin1 = copy(xmin)
#    modified = false
#    for i in eachindex(xmin1)
#        if 0. < abs(xmin1[i]) < 1.e-5
#            modified = true
#            xmin1[i] = 0.
#        end
#    end
#    if modified  # branch not tested
#        m[:θ] = xmin1
#        ff = objective(m)
#        if ff ≤ (fmin + 1.e-5)  # zero components if increase in objective is negligible
#            fmin = ff
#            copy!(xmin,xmin1)
#        else
#            m[:θ] = xmin
#        end
#    end
    m.LMM.opt = OptSummary(βθ,xmin,fmin,feval,optimizer)
    restoreX!(m)
    if verbose
        println(ret)
    end
    m
end

function VarCorr(m::GeneralizedLinearMixedModel)
    Λ, trms = m.LMM.Λ, m.LMM.trms
    VarCorr(Λ, [string(trms[i].fnm) for i in eachindex(Λ)],
        [trms[i].cnms for i in eachindex(Λ)], NaN)
end

function Base.show{T,D,L}(io::IO, m::GeneralizedLinearMixedModel{T,D,L}) # not tested
    println(io, "Generalized Linear Mixed Model fit by minimizing the Laplace approximation to the deviance")
    println(io, "  ", m.LMM.formula)
    println(io, "  Distribution: ", D)
    println(io, "  Link: ", L, "\n")
    println(io, string("  Deviance (Laplace approximation): ", @sprintf("%.4f", LaplaceDeviance(m))), "\n")

    show(io,VarCorr(m))
    gl = grplevels(m.LMM)
    print(io, "\n Number of obs: ", length(m.y), "; levels of grouping factors: ", gl[1])
    for l in gl[2:end]
        print(io, ", ", l)
    end
    println(io)
    println(io, "\nFixed-effects parameters:")
    show(io, coeftable(m))
end

varest{T <: AbstractFloat}(m::GeneralizedLinearMixedModel{T}) = one(T)


text:
@testset "LOAD" begin

  for i in 0:4
    # test loading pre-saved files
    bts = load("./test_pkl/builtin_type_p$(i).pkl")
    @test builtin_type_samples == bts

    # test loading from string
    sbts = loads(py"pystores(builtin_type_samples, $i)")
    @test builtin_type_samples == sbts
  end

  mktempdir() do path
    for i in 0:4
      # test loading from directly saved file
      file = joinpath(path, "pybuiltin_type_p$i.pkl")
      py"pystore($file, builtin_type_samples, $i)"
      @test builtin_type_samples == load(file)
    end
  end

end


text:
export padding

function padding(str::AbstractString, width::Int; background=' ')
    str_width = textwidth(str)
    (width < str_width) && (throw(BoundsError))

    str *= background^(width-str_width)

    return str
end

function padding_vertical(str::String, height::Int)
    str_height = length(findall("\n", str))+1
    (height < str_height) && (throw(BoundsError))

    str *= '\n'^(height-str_height)

    return str
end


text:
"""
    This function returns the hellinger distance calculated in n (10 by default) different intervals of a trajectory together with its error. The normalization factor associated with the distribution in `q` is passed as an argument.
    Example:
    ```
    julia> hellingerdistance(data, parameters, normalizationfactor, n)
    ```
    """
function hellingerdistance(data::Matrix{Float64}, parameters::Parameters, normalizationfactor::Float64, n = 10)

    beta = 1./parameters.T

    function jointdistribution(x::Vector{Float64})
        q,p,z = x
        rhoq = 1./normalizationfactor*exp(-beta*parameters.potential.f(q))
        rhop = sqrt(beta/(2.*pi))*exp(-beta*p.^2/2.)
        rhoz = parameters.thermo.distribution(z)

        rhoq*rhop*rhoz
    end

    time = zeros(n)
    hell = zeros(n)
    error = zeros(n)
    data = data[1:100:end,:] ##The 100 factor increases the sampling time by this factor
    len = size(data)[1] - 1

    for i in 1:n
        data1 = data[1:Int(len/n*i),:]
        time[i] = data1[:,1][end]
        dens = sm[:nonparametric][:KDEMultivariate](data = data1[:,2:4], var_type = "ccc")

        function kernel(x::Vector{Float64})
            dens[:pdf](x)[1]
        end

        k_2(x) = kernel(x)^0.5
        j_2(x) = jointdistribution(x)^0.5

        hellinger(x) = 2.*(k_2(x) - j_2(x))^2.

        q1, q2 = (minimum(data1[:,2]),maximum(data1[:,2]))
        p1, p2 = (minimum(data1[:,3]),maximum(data1[:,3]))
        z1, z2 = (minimum(data1[:,4]),maximum(data1[:,4]))
        hell[i], error[i] = hcubature(hellinger, [q1,p1,z1],[q2,p2,z2], reltol = 1.0e-2)
        println(i)
    end

    results = hcat(time, hell)
    hellinger_results = hcat(results, error)
end

max_stars_repo_name: dapias/STODs

max_issues_repo_name: dapias/STODs

max_forks_repo_name: dapias/STODs


text:
export paramscan

"""
    paramscan(parameters, initialize; kwargs...) → adf, mdf

Perform a parameter scan of a ABM simulation output by collecting data from all
parameter combinations into dataframes (one for agent data, one for model data).
The dataframes columns are both the collected data (as in [`run!`](@ref)) but also the
input parameter values used.

`parameters` is a dictionary with key type `Symbol` which contains various parameters that
will be scanned over (as well as other parameters that remain constant).
This function uses `DrWatson`'s [`dict_list`](https://juliadynamics.github.io/DrWatson.jl/dev/run&list/#DrWatson.dict_list)
convention. This means that every entry of `parameters` that is a `Vector`
contains many parameters and thus is scanned. All other entries of
`parameters` that are not `Vector`s are not expanded in the scan.

The second argument `initialize` is a function that creates an ABM and returns it.
It must accept keyword arguments which are the *keys* of the `parameters` dictionary.
Since the user decides how to use input arguments to make an ABM, `parameters` can be
used to affect model properties, space type and creation as well as agent properties,
see the example below.

## Keywords
The following keywords modify the `paramscan` function:

- `include_constants::Bool = false`: by default, only the varying parameters (Vector in
  `parameters`) will be included in the output `DataFrame`. If `true`, constant parameters
  (non-Vector in `parameteres`) will also be included.
- `parallel::Bool = false` whether `Distributed.pmap` is invoked to run simulations
  in parallel. This must be used in conjunction with `@everywhere` (see
  [Performance Tips](@ref)).

All other keywords are propagated into [`run!`](@ref).
Furthermore, `agent_step!, model_step!, n` are also keywords here, that are given
to [`run!`](@ref) as arguments. Naturally, `agent_step!, model_step!, n` and at least one
of `adata, mdata` are mandatory.
The `adata, mdata` lists shouldn't contain the parameters that are already in
the `parameters` dictionary to avoid duplication.

## Example
A runnable example that uses `paramscan` is shown in [Schelling's segregation model](@ref).
There, we define
```julia
function initialize(; numagents = 320, griddims = (20, 20), min_to_be_happy = 3)
    space = GridSpace(griddims, moore = true)
    properties = Dict(:min_to_be_happy => min_to_be_happy)
    model = ABM(SchellingAgent, space;
                properties = properties, scheduler = schedule_randomly)
    for n in 1:numagents
        agent = SchellingAgent(n, (1, 1), false, n < numagents / 2 ? 1 : 2)
        add_agent_single!(agent, model)
    end
    return model
end
```
and do a parameter scan by doing:
```julia
happyperc(moods) = count(moods) / length(moods)
adata = [(:mood, happyperc)]

parameters = Dict(
    :min_to_be_happy => collect(2:5), # expanded
    :numagents => [200, 300],         # expanded
    :griddims => (20, 20),            # not Vector = not expanded
)

adf, _ = paramscan(parameters, initialize; adata, agent_step!, n = 3)
```
"""
function paramscan(
    parameters::Dict,
    initialize;
    include_constants::Bool = false,
    parallel::Bool = false,
    agent_step! = dummystep,
    model_step! = dummystep,
    n = 1,
    kwargs...,
)

    if include_constants
        output_params = collect(keys(parameters))
    else
        output_params = [k for (k, v) in parameters if typeof(v) <: Vector]
    end

    combs = dict_list(parameters)

    if parallel
        all_data = ProgressMeter.@showprogress pmap(combs) do i
            run_single(i, output_params, initialize; agent_step!, model_step!, n, kwargs...)
        end
    else
        all_data = ProgressMeter.@showprogress map(combs) do i
            run_single(i, output_params, initialize; agent_step!, model_step!, n, kwargs...)
        end
    end

    df_agent = DataFrame()
    df_model = DataFrame()
    for (df1, df2) in all_data
        append!(df_agent, df1)
        append!(df_model, df2)
    end

    return df_agent, df_model
end

# This function is taken from DrWatson:
function dict_list(c::Dict)
    iterable_fields = filter(k -> typeof(c[k]) <: Vector, keys(c))
    non_iterables = setdiff(keys(c), iterable_fields)

    iterable_dict = Dict(iterable_fields .=> getindex.(Ref(c), iterable_fields))
    non_iterable_dict = Dict(non_iterables .=> getindex.(Ref(c), non_iterables))

    vec(map(Iterators.product(values(iterable_dict)...)) do vals
        dd = Dict(keys(iterable_dict) .=> vals)
        if isempty(non_iterable_dict)
            dd
        elseif isempty(iterable_dict)
            non_iterable_dict
        else
            merge(non_iterable_dict, dd)
        end
    end)
end

function run_single(
    param_dict::Dict,
    output_params::Vector{Symbol},
    initialize;
    agent_step! = dummystep,
    model_step! = dummystep,
    n = 1,
    kwargs...,
)
    model = initialize(; param_dict...)
    df_agent_single, df_model_single = run!(model, agent_step!, model_step!, n; kwargs...)
    output_params_dict = filter(j -> first(j) in output_params, param_dict)
    insertcols!(df_agent_single, output_params_dict...)
    insertcols!(df_model_single, output_params_dict...)
    return (df_agent_single, df_model_single)
end


text:
const LIBPATH = joinpath(@__DIR__, "../lib")
const AM_pat = r"([SPDFGHI]{1,2})\s+?(\d++)"
const prim_pat = r"([+-]?\d*?\.\d+[D+-]{0,2}\d\d)\s+?([+-]?\d*?\.\d+[D+-]{0,2}+\d\d)"
const prim_pat3 = r"([+-]?\d*?\.\d+[D+-]{0,2}\d\d)\s+?([+-]?\d*?\.\d+[D+-]{0,2}+\d\d)\s+?([+-]?\d*?\.\d+[D+-]{0,2}+\d\d)"
const AMDict = Dict(
        "S" => 0,
        "P" => 1,
        "D" => 2,
        "F" => 3,
        "G" => 4,
        "H" => 5,
        "I" => 6,
    )

"""
    GaussinBasis.read_basisset(bname::String, AtomSymbol::String)

Returns an array of BasisFunction objects given an atom (AtomicSymbol) and basis set name (bname).
"""
function read_basisset(bname::String, AtomSymbol::String)

    clean_bname = replace(bname, "*"=>"_st_")
    file_path = joinpath(LIBPATH, clean_bname*".gbs")

    if !(isfile(file_path))
        throw(ArgumentError("Basis set file for $bname was not found."))
    end

    info = extract_atom_from_bs(file_path, AtomSymbol)

    # Split info into basis
    BasisStrings = [info[1]*"\n"]
    for line in info[2:end]
        if occursin(AM_pat, line)
            push!(BasisStrings, line*"\n")
        else
            BasisStrings[end] *= line * "\n"
        end
    end

    out = BasisFunction[]

    for b in BasisStrings
        r = r"[SPDFGHI]{2}"
        if occursin(r, b)
            bf1, bf2 = two_basis_from_string(b)
            normalize_basisfunction!(bf1)
            normalize_basisfunction!(bf2)
            push!(out, bf1)
            push!(out, bf2)
        else
            bf = basis_from_string(b)
            normalize_basisfunction!(bf)
            push!(out, bf)
        end
    end

    return out
end

"""
    GaussinBasis.extract_atom_from_bs(file_path::String, AtomSymbol::String)

Extract the portion of the basis set file (file_path) relevant for the given atom in AtomSymbol.
"""
function extract_atom_from_bs(file_path::String, AtomSymbol::String)

    atom_pat = Regex("$(AtomSymbol)\\s+?0")
    # This flag indicates if the atom was found within the file
    flag_atom = false

    # Output is an array of strings (corresponding to the lines of the block matching the desired atom)
    out = String[]

    # Iterate through the lines of the file until the atom is found
    for line in eachline(file_path)

        if occursin(atom_pat, line)
            flag_atom = true
            continue
        end

        if flag_atom
            occursin(raw"****", line) ? break : nothing
            push!(out, String(line))
        end
    end

    if flag_atom
        return out
    else
        throw(FermiException("Atom $AtomSymbol not found in $file_path."))
    end
end

"""
    GaussinBasis.basis_from_string(bstring::String)

From a String block representing a Basis Function in gbs format, produces a BasisFunction object.
For the special case of two basis being described within the same block (e.g. SP blocks) see `two_basis_from_string`
"""
function basis_from_string(bstring::String)
    lines = split(strip(bstring), "\n")
    head = lines[1]
    m = match(AM_pat, head)

    if m === nothing
        throw(ArgumentError("cannot parse basis set file line: $head"))
    end

    AMsymbol = String(m.captures[1])
    l = AMDict[AMsymbol]
    nprim = parse(Int, m.captures[2])
    coef = zeros(nprim)
    exp = zeros(nprim)

    for i in eachindex(exp)
        m = match(prim_pat, lines[i+1])

        if m === nothing
            throw(ArgumentError("cannot parse basis set file line: $(lines[i+1])"))
        end

        e = replace(m.captures[1], "D"=>"E")
        c = replace(m.captures[2], "D"=>"E")
        coef[i] = parse(Float64, c)
        exp[i] = parse(Float64, e)
    end

    return BasisFunction(l, coef, exp)
end

"""
    GaussinBasis.two_basis_from_string(bstring::String)

From a String block representing two Basis Function (e.g. SP blocks) in gbs format, produces two BasisFunction objects.
For the case of a single basis being described within the block see `basis_from_string`
"""
function two_basis_from_string(bstring::String)
    lines = split(strip(bstring), "\n")
    head = lines[1]
    m = match(AM_pat, head)
    if m === nothing
        throw(ArgumentError("cannot parse basis set file line: $head"))
    end
    AMsymbol = String(m.captures[1])

    length(AMsymbol) == 2 || throw(FermiException("cannot extract two basis from $AMsymbol function"))

    l1 = AMDict[AMsymbol[1]*""]
    l2 = AMDict[AMsymbol[2]*""]

    nprim = parse(Int, m.captures[2])
    coef1 = zeros(nprim)
    coef2 = zeros(nprim)
    exp = zeros(nprim)

    for i in eachindex(exp)
        m = match(prim_pat3, lines[i+1])

        if m === nothing
            throw(ArgumentError("cannot parse basis set file line: $(lines[i+1])"))
        end

        e = replace(m.captures[1], "D"=>"E")
        c1 = replace(m.captures[2], "D"=>"E")
        c2 = replace(m.captures[3], "D"=>"E")
        coef1[i] = parse(Float64, c1)
        coef2[i] = parse(Float64, c2)
        exp[i] = parse(Float64, e)
    end

    return BasisFunction(l1, coef1, exp), BasisFunction(l2, coef2, exp)
end


text:
const PROPERTIES = [
    :acquisition_duration,
    :anatomical_landmark_coordinates,
    :anatomical_landmark_coordinate_system,
    :anatomical_landmark_coordinate_units,
    :anatomical_landmark_coordinate_system_description,
    :CogAtlasID,
    :CogPOID,
    :coil_combination_method,
    :contrast_bolus_ingredient,
    :delay_after_trigger,
    :delay_time,
    :device_serial_number,
    :dewar_position,
    :dwell_time,
    :echo_time,
    :effective_echo_spacing,
    :eeg_electrode_groups,
    :eeg_ground,
    :eeg_placement_scheme,
    :electrical_stimulation,
    :electrical_stimulation_parameters,
    :epoch_length,
    :fiducial_description,
    :flip_angle,
    :gradient_set_type,
    :head_coil_coordinates,
    :head_coil_coordinate_system,
    :head_coil_coordinate_system_description,
    :head_coil_coordinate_units,
    :institution_address,
    :institution_name,
    :institutional_department_name,
    :instructions,
    :intended_for,
    :inversion_time,
    :magnetic_field_strength,
    :manufacturer,
    :manufacturer_model_name,
    :matrix_coil_mode,
    :mr_transmit_coil_sequence,
    :multiband_acceleration_factor,
    :negative_contrast,
    :number_of_volumes_discarded_by_scanner,
    :number_of_volumes_discarded_by_user,
    :number_shots,
    :nonlinear_gradient_correction,
    :parallel_acquisition_technique,
    :parallel_reduction_factor,
    :partial_fourier,
    :partial_fourier_direction,
    :phase_encoding_direction,
    :power_line_frequency,
    :pulse_sequence,
    :pulse_sequence_details,
    :pulse_sequence_type,
    :receive_coil_active_elements,
    :receive_coil_name,
    :recording_duration,
    :recording_type,
    :repetition_time,
    :sampling_frequency,
    :scan_options,
    :scanning_sequence,
    :sequence_name,
    :sequence_varient,
    :subject_artefact_description,
    :slice_encoding_direction,
    :slice_timing,
    :software_versions,
    :start_time,
    :station_name,
    :task_description,
    :task_name,
    :total_readout_time,
    :volume_timing,
   ]

NO_SET_PROPERTIES = [:anatomical_landmark_coordinate_system,
                     :anatomical_landmark_coordinate_units,
                     :head_coil_coordinate_system,
                     :head_coil_coordinate_units,
                     :start_time,
                    ]

function generate_getproperty(props)
    e = :(return getindex(properties(x),  s))
    for p in props
        e = :(if $(Expr(:call, :(===), :s, QuoteNode(p)))
              return $p(x)
          else
              $e
          end
         )
    end
    @eval begin
        function neuroproperty(x, s::Symbol)
            Base.@_inline_meta
            $e
        end
    end
    return nothing
end
generate_getproperty(PROPERTIES)

function generate_setproperty(props, noprops)
    e = :(return setindex!(properties(x), val, s))
    for p in props
        if p in noprops
            continue
        else
            pp = Symbol(p, :!)
            e = :(if $(Expr(:call, :(===), :s, QuoteNode(p)))
                  return $pp(x, val)
              else
                  $e
              end
             )
        end
    end
    @eval begin
        function neuroproperty!(x, s::Symbol, val)
            Base.@_inline_meta
            $e
        end
    end
    return nothing
end
generate_setproperty(PROPERTIES, NO_SET_PROPERTIES)


text:
evaluate(H::AbstractHomotopy, x, t) = error("`evaluate` not implemented for $(typeof(H))<:AbstractHomotopy")
substitute(H::AbstractHomotopy, pair::Pair{Symbol,<:Number}) = error("`substitute` not implemented for $(typeof(H))<:AbstractHomotopy")
startsystem(H::AbstractHomotopy) = error("`startsystem` not implemented for $(typeof(H))<:AbstractHomotopy")
targetsystem(H::AbstractHomotopy) = error("`targetsystem` not implemented for $(typeof(H))<:AbstractHomotopy")
nvariables(H::AbstractHomotopy) = error("`nvars` not implemented for $(typeof(H))<:AbstractHomotopy")
removepoly(H::AbstractHomotopy, i::Int) = error("`removepoly` not implemented for $(typeof(H))<:AbstractHomotopy")
degrees(H::AbstractHomotopy) = error("`degrees` not implemented for $(typeof(H))<:AbstractHomotopy")
differentiate(H::AbstractHomotopy) = error("`differentiate` not implemented for $(typeof(H))<:AbstractHomotopy")
dt(H::AbstractHomotopy) = error("`dt` not implemented for $(typeof(H))<:AbstractHomotopy")

ishomogenous(H::AbstractHomotopy) = error("`ishomogenous` not implemented for $(typeof(H))<:AbstractHomotopy")
homogenize(H::AbstractHomotopy) = error("`homogenize` not implemented for $(typeof(H))<:AbstractHomotopy")
homogenized(H::AbstractHomotopy) = error("`homogenized` not implemented for $(typeof(H))<:AbstractHomotopy")

weylnorm(H::AbstractHomotopy, t::Float64) = error("`weyl_norm` not implemented for $(typeof(H)) for AbstractHomotopy")

@inline ∂t(H::AbstractHomotopy) = dt(H)
eltype(H::AbstractHomotopy{T}) where {T<:Number} = T
coefftype(H::AbstractHomotopy{T}) where {T<:Number} = T


text:
# TRAJECTORY CALCULATION #

const DEFAULT_DT = 0.001

"""
   compute_trajectory(u,v,X₀::Vector,Tmax[,Δt=0.001])

Calculate the trajectory of a tracer particle with initial location(s) `X₀`, which
can be specified as either a single vector `[x0,y0]` or a vector of vectors
for multiple tracer particles. The arguments
`u` and `v` are interpolated velocity field components, `Tmax` is the final
integration time, and `Δt` is the time step size. The output is the solution
structure for the `OrdinaryDiffEq` package (or, for multiple particles, a vector
of such solution structures).
"""
function compute_trajectory(ufield::AbstractInterpolation{T,2},
                            vfield::AbstractInterpolation{T,2},
                            X₀::Vector{S},Tmax::Real;Δt::Real=DEFAULT_DT) where {T,S<:Real}

  vfcn!(dR,R,p,t) = _vfcn_autonomous!(dR,R,p,t,ufield,vfield)

  sol = _solve_trajectory(vfcn!,X₀,Tmax,Δt)
  return sol

end

function compute_trajectory(ufield::AbstractInterpolation{T,2},vfield::AbstractInterpolation{T,2},
   pts::Vector{Vector{S}},Tmax;Δt=DEFAULT_DT) where {T,S<:Real}

  sol_array = ODESolution[]
  for X₀ in pts
    sol = compute_trajectory(ufield,vfield,X₀,Tmax,Δt=Δt)
    push!(sol_array,sol)
  end
  return sol_array

end

function _solve_trajectory(vfcn!,u0,Tmax,Δt)
  Path = ODEProblem(vfcn!,u0,(0.0,Tmax))
  sol = solve(Path,Tsit5(), dt = Δt, maxiters = 1e8, adaptive = false, dense = false)
end

function _vfcn_autonomous!(dR,R,p,t,u,v)
 dR[1] = u(R[1],R[2])
 dR[2] = v(R[1],R[2])

 return dR
end


text:
#' @exported
#' @description
#'
#' Returns a string summary of an AbstractDataFrame in a standardized
#' form. For example, a standard DataFrame with 10 rows and 5 columns
#' will be summarized as "10x5 DataFrame".
#'
#' @param adf::AbstractDataFrame The AbstractDataFrame to be summarized.
#'
#' @returns res::UTF8String The summary of `adf`.
#'
#' @examples
#'
#' summary(DataFrame(A = 1:10))
function Base.summary(adf::AbstractDataFrame) # -> UTF8String
    nrows, ncols = size(adf)
    return utf8(@sprintf "%dx%d %s" nrows ncols typeof(adf))
end

#' @description
#'
#' Determine the number of UTF8 characters that would be used to
#' render a value.
#'
#' @param x::Any A value whose string width will be computed.
#'
#' @returns w::Int The width of the string.
#'
#' @examples
#'
#' ourstrwidth("abc")
#' ourstrwidth(10000)
begin
    local io = IOBuffer(Array(Uint8, 80), true, true)
    global ourstrwidth
    function ourstrwidth(x::Any) # -> Int
        truncate(io, 0)
        ourshowcompact(io, x)
        return position(io)
    end
    ourstrwidth(x::String) = strwidth(x) + 2 # -> Int
    ourstrwidth(s::Symbol) = int(ccall(:u8_strwidth,
                                       Csize_t,
                                       (Ptr{Uint8}, ),
                                       convert(Ptr{Uint8}, s)))
end

#' @description
#'
#' Render a value to an IO object in a compact format. Unlike
#' Base.showcompact, we render strings without surrounding quote
#' marks.
#'
#' @param io::IO An IO object to be printed to.
#' @param x::Any A value to be printed.
#'
#' @returns x::Nothing A `nothing` value.
#'
#' @examples
#'
#' ourshowcompact(STDOUT, "abc")
#' ourshowcompact(STDOUT, 10000)
ourshowcompact(io::IO, x::Any) = showcompact(io, x) # -> Nothing
ourshowcompact(io::IO, x::String) = showcompact(io, x) # -> Nothing
ourshowcompact(io::IO, x::Symbol) = print(io, x) # -> Nothing

#' @description
#'
#' Calculates, for each column of an AbstractDataFrame, the maximum
#' string width used to render either the name of that column or the
#' longest entry in that column -- among the rows of the AbstractDataFrame
#' will be rendered to IO. The widths for all columns are returned as a
#' vector.
#'
#' NOTE: The last entry of the result vector is the string width of the
#'       implicit row ID column contained in every AbstractDataFrame.
#'
#' @param adf::AbstractDataFrame The AbstractDataFrame whose columns will be
#'        printed.
#' @param rowindices1::AbstractVector{Int} A set of indices of the first
#'        chunk of the AbstractDataFrame that would be rendered to IO.
#' @param rowindices2::AbstractVector{Int} A set of indices of the second
#'        chunk of the AbstractDataFrame that would be rendered to IO. Can
#'        be empty if the AbstractDataFrame would be printed without any
#'        ellipses.
#' @param rowlabel::String The label that will be used when rendered the
#'        numeric ID's of each row. Typically, this will be set to "Row #".
#'
#' @returns widths::Vector{Int} The maximum string widths required to render
#'          each column, including that column's name.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "yy", "z"])
#' maxwidths = getmaxwidths(df, 1:1, 3:3, "Row #")
function getmaxwidths(adf::AbstractDataFrame,
                      rowindices1::AbstractVector{Int},
                      rowindices2::AbstractVector{Int},
                      rowlabel::Symbol) # -> Vector{Int}
    ncols = size(adf, 2)
    cnames = names(adf)
    maxwidths = Array(Int, ncols + 1)

    # TODO: Move this definition somewhere else
    NAstrwidth = 2
    undefstrwidth = ourstrwidth(Base.undef_ref_str)

    for j in 1:ncols
        # (1) Consider length of column name
        maxwidths[j] = ourstrwidth(cnames[j])

        # (2) Consider length of longest entry in that column
        col = adf[j]
        for indices in (rowindices1, rowindices2)
            for i in indices
                if isna(col, i)
                    maxwidths[j] = max(maxwidths[j], NAstrwidth)
                else
                    try
                        maxwidths[j] = max(maxwidths[j], ourstrwidth(col[i]))
                    catch
                        maxwidths[j] = max(maxwidths[j], undefstrwidth)
                    end
                end
            end
        end
    end

    rowmaxwidth1 = isempty(rowindices1) ? 0 : ndigits(maximum(rowindices1))
    rowmaxwidth2 = isempty(rowindices2) ? 0 : ndigits(maximum(rowindices2))

    maxwidths[ncols + 1] = max(max(rowmaxwidth1,
                                   rowmaxwidth2),
                               ourstrwidth(rowlabel))

    return maxwidths
end

#' @description
#'
#' Given the maximum widths required to render each column of an
#' AbstractDataFrame, this returns the total number of UTF8 characters
#' that would be required to render an entire row to an IO system.
#'
#' NOTE: This width includes the whitespace and special characters used to
#'       pretty print the AbstractDataFrame.
#'
#' @param maxwidths::Vector{Int} The maximum width needed to render each
#'        column of an AbstractDataFrame.
#'
#' @returns totalwidth::Int The total width required to render a complete row
#'          of the AbstractDataFrame for which `maxwidths` was computed.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "yy", "z"])
#' maxwidths = getmaxwidths(df, 1:1, 3:3, "Row #")
#' totalwidth = getprintedwidth(maxwidths))
function getprintedwidth(maxwidths::Vector{Int}) # -> Int
    # Include length of line-initial |
    totalwidth = 1
    for i in 1:length(maxwidths)
        # Include length of field + 2 spaces + trailing |
        totalwidth += maxwidths[i] + 3
    end
    return totalwidth
end

#' @description
#'
#' When rendering an AbstractDataFrame to a REPL window in chunks, each of
#' which will fit within the width of the REPL window, this function will
#' return the indices of the columns that should be included in each chunk.
#'
#' NOTE: The resulting bounds should be interpreted as follows: the
#'       i-th chunk bound is the index MINUS 1 of the first column in the
#'       i-th chunk. The (i + 1)-th chunk bound is the EXACT index of the
#'       last column in the i-th chunk. For example, the bounds [0, 3, 5]
#'       imply that the first chunk contains columns 1-3 and the second chunk
#'       contains columns 4-5.
#'
#' @param maxwidths::Vector{Int} The maximum width needed to render each
#'        column of an AbstractDataFrame.
#' @param splitchunks::Bool Should the output be split into chunks at all or
#'        should only one chunk be constructed for the entire
#'        AbstractDataFrame?
#'
#' @returns chunkbounds::Vector{Int} The bounds of each chunk of columns.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "yy", "z"])
#' maxwidths = getmaxwidths(df, 1:1, 3:3, "Row #")
#' chunkbounds = getchunkbounds(maxwidths, true)
function getchunkbounds(maxwidths::Vector{Int},
                        splitchunks::Bool) # -> Vector{Int}
    ncols = length(maxwidths) - 1
    rowmaxwidth = maxwidths[ncols + 1]
    availablewidth = Base.tty_cols()
    if splitchunks
        chunkbounds = [0]
        # Include 2 spaces + 2 | characters for row/col label
        totalwidth = rowmaxwidth + 4
        for j in 1:ncols
            # Include 2 spaces + | character in per-column character count
            totalwidth += maxwidths[j] + 3
            if totalwidth > availablewidth
                push!(chunkbounds, j - 1)
                totalwidth = rowmaxwidth + 4 + maxwidths[j] + 3
            end
        end
        push!(chunkbounds, ncols)
    else
        chunkbounds = [0, ncols]
    end
    return chunkbounds
end

#' @description
#'
#' Render a subset of rows and columns of an AbstractDataFrame to an
#' IO system. For chunked printing, this function is used to print a
#' single chunk, starting from the first indicated column and ending with
#' the last indicated column. Assumes that the maximum string widths
#' required for printing have been precomputed.
#'
#' @param io::IO The IO system to which `adf` will be printed.
#' @param adf::AbstractDataFrame An AbstractDataFrame.
#' @param rowindices::AbstractVector{Int} The indices of the subset of rows
#'        that will be rendered to `io`.
#' @param maxwidths::Vector{Int} The pre-computed maximum string width
#'        required to render each column.
#' @param leftcol::Int The index of the first column in a chunk to be
#'        rendered.
#' @param rightcol::Int The index of the last column in a chunk to be
#'        rendered.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' showrowindices(STDOUT, df, 1:2, [1, 1, 5], 1, 2)
function showrowindices(io::IO,
                        adf::AbstractDataFrame,
                        rowindices::AbstractVector{Int},
                        maxwidths::Vector{Int},
                        leftcol::Int,
                        rightcol::Int) # -> Nothing
    rowmaxwidth = maxwidths[end]

    for i in rowindices
        # Print row ID
        @printf io "| %d" i
        padding = rowmaxwidth - ndigits(i)
        for itr in 1:padding
            write(io, ' ')
        end
        print(io, " | ")
        # Print DataFrame entry
        for j in leftcol:rightcol
            strlen = 0
            try
                strlen = ourstrwidth(adf[i, j])
                ourshowcompact(io, adf[i, j])
            catch
                strlen = ourstrwidth(Base.undef_ref_str)
                ourshowcompact(io, Base.undef_ref_str)
            end
            padding = maxwidths[j] - strlen
            for itr in 1:padding
                write(io, ' ')
            end
            if j == rightcol
                if i == rowindices[end]
                    print(io, " |")
                else
                    print(io, " |\n")
                end
            else
                print(io, " | ")
            end
        end
    end
    return
end

#' @description
#'
#' Render a subset of rows (possibly in chunks) of an AbstractDataFrame to an
#' IO system. Users can control
#'
#' NOTE: The value of `maxwidths[end]` must be the string width of
#' `rowlabel`.
#'
#' @param io::IO The IO system to which `adf` will be printed.
#' @param adf::AbstractDataFrame An AbstractDataFrame.
#' @param rowindices1::AbstractVector{Int} The indices of the first subset
#'        of rows to be rendered.
#' @param rowindices2::AbstractVector{Int} The indices of the second subset
#'        of rows to be rendered. An ellipsis will be printed before
#'        rendering this second subset of rows.
#' @param maxwidths::Vector{Int} The pre-computed maximum string width
#'        required to render each column.
#' @param splitchunks::Bool Should the printing of the AbstractDataFrame
#'        be done in chunks? Defaults to `false`.
#' @param rowlabel::Symbol What label should be printed when rendering the
#'        numeric ID's of each row? Defaults to `"Row #"`.
#' @param displaysummary::Bool Should a brief string summary of the
#'        AbstractDataFrame be rendered to the IO system before printing the
#'        contents of the renderable rows? Defaults to `true`.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' showrows(STDOUT, df, 1:2, 3:3, [1, 1, 5], false, "Row #", true)
function showrows(io::IO,
                  adf::AbstractDataFrame,
                  rowindices1::AbstractVector{Int},
                  rowindices2::AbstractVector{Int},
                  maxwidths::Vector{Int},
                  splitchunks::Bool = false,
                  rowlabel::Symbol = symbol("Row #"),
                  displaysummary::Bool = true) # -> Nothing
    ncols = size(adf, 2)
    cnames = names(adf)

    if displaysummary
        println(io, summary(adf))
    end

    if isempty(rowindices1)
        return
    end

    rowmaxwidth = maxwidths[ncols + 1]
    chunkbounds = getchunkbounds(maxwidths, splitchunks)
    nchunks = length(chunkbounds) - 1

    for chunkindex in 1:nchunks
        leftcol = chunkbounds[chunkindex] + 1
        rightcol = chunkbounds[chunkindex + 1]

        # Print table bounding line
        write(io, '|')
        for itr in 1:(rowmaxwidth + 2)
            write(io, '-')
        end
        write(io, '|')
        for j in leftcol:rightcol
            for itr in 1:(maxwidths[j] + 2)
                write(io, '-')
            end
            write(io, '|')
        end
        write(io, '\n')

        # Print column names
        @printf io "| %s" rowlabel
        padding = rowmaxwidth - ourstrwidth(rowlabel)
        for itr in 1:padding
            write(io, ' ')
        end
        @printf io " | "
        for j in leftcol:rightcol
            s = cnames[j]
            ourshowcompact(io, s)
            padding = maxwidths[j] - ourstrwidth(s)
            for itr in 1:padding
                write(io, ' ')
            end
            if j == rightcol
                print(io, " |\n")
            else
                print(io, " | ")
            end
        end

        # Print main table body, potentially in two abbreviated sections
        showrowindices(io,
                       adf,
                       rowindices1,
                       maxwidths,
                       leftcol,
                       rightcol)
        if !isempty(rowindices2)
            print(io, "\n⋮\n")
            showrowindices(io,
                           adf,
                           rowindices2,
                           maxwidths,
                           leftcol,
                           rightcol)
        end

        # Print newlines to separate chunks
        if chunkindex < nchunks
            print(io, "\n\n")
        end
    end

    return
end

#' @exported
#' @description
#'
#' Render an AbstractDataFrame to an IO system. The specific visual
#' representation chosen depends on the width of the REPL window
#' from which the call to `show` derives. If the DataFrame could not
#' be rendered without splitting the output into chunks, a summary of the
#' columns is rendered instead of rendering the raw data. This dynamic
#' response to screen width can be configured using the argument
#' `splitchunks`.
#'
#' @param io::IO The IO system to which `adf` will be printed.
#' @param adf::AbstractDataFrame An AbstractDataFrame.
#' @param splitchunks::Bool Should the printing of the AbstractDataFrame
#'        be done in chunks? Defaults to `false`.
#' @param rowlabel::Symbol What label should be printed when rendering the
#'        numeric ID's of each row? Defaults to `"Row #"`.
#' @param displaysummary::Bool Should a brief string summary of the
#'        AbstractDataFrame be rendered to the IO system before printing the
#'        contents of the renderable rows? Defaults to `true`.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' show(STDOUT, df, false, "Row #", true)
function Base.show(io::IO,
                   adf::AbstractDataFrame,
                   splitchunks::Bool = false,
                   rowlabel::Symbol = symbol("Row #"),
                   displaysummary::Bool = true) # -> Nothing
    nrows = size(adf, 1)
    availableheight = Base.tty_rows() - 5
    nrowssubset = fld(availableheight, 2)
    bound = min(nrowssubset - 1, nrows)
    if nrows <= availableheight
        rowindices1 = 1:nrows
        rowindices2 = 1:0
    else
        rowindices1 = 1:bound
        rowindices2 = max(bound + 1, nrows - nrowssubset + 1):nrows
    end
    maxwidths = getmaxwidths(adf, rowindices1, rowindices2, rowlabel)
    width = getprintedwidth(maxwidths)
    if width > Base.tty_cols() && !splitchunks
        showcols(io, adf)
    else
        showrows(io,
                 adf,
                 rowindices1,
                 rowindices2,
                 maxwidths,
                 splitchunks,
                 rowlabel,
                 displaysummary)
    end
    return
end

#' @exported
#' @description
#'
#' Render an AbstractDataFrame to STDOUT with or without chunking. See
#' other `show` documentation for details. This is mainly used to force
#' showing the AbstractDataFrame in chunks.
#'
#' @param adf::AbstractDataFrame An AbstractDataFrame.
#' @param splitchunks::Bool Should the printing of the AbstractDataFrame
#'        be done in chunks? Defaults to `false`.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' show(df, true)
function Base.show(adf::AbstractDataFrame,
                   splitchunks::Bool = false) # -> Nothing
    return show(STDOUT, adf, splitchunks)
end

#' @exported
#' @description
#'
#' Render a DataFrameRow to an IO system. Each column of the DataFrameRow
#' is printed on a separate line.
#'
#' @param io::IO The IO system where rendering will take place.
#' @param r::DataFrameRow The DataFrameRow to be rendered to `io`.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' for r in eachrow(df)
#'     show(STDOUT, r)
#' end
function Base.show(io::IO, r::DataFrameRow)
    labelwidth = mapreduce(length, max, names(r)) + 2
    @printf("DataFrameRow (row %d)\n", r.row)
    for (label, value) in r
        println(io, rpad(label, labelwidth, ' '), value)
    end
end

#' @exported
#' @description
#'
#' Render a DataFrameRow to STDOUT. See other `show` documentation for
#' details.
#'
#' @param r::DataFrameRow The DataFrameRow to be rendered to `io`.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' for r in eachrow(df)
#'     show(r)
#' end
Base.show(row::DataFrameRow) = show(STDOUT, row)

#' @exported
#' @description
#'
#' Render all of the rows of an AbstractDataFrame to an IO system. See
#' `show` documentation for details.
#'
#' @param io::IO The IO system to which `adf` will be printed.
#' @param adf::AbstractDataFrame An AbstractDataFrame.
#' @param splitchunks::Bool Should the printing of the AbstractDataFrame
#'        be done in chunks? Defaults to `false`.
#' @param rowlabel::Symbol What label should be printed when rendering the
#'        numeric ID's of each row? Defaults to `"Row #"`.
#' @param displaysummary::Bool Should a brief string summary of the
#'        AbstractDataFrame be rendered to the IO system before printing the
#'        contents of the renderable rows? Defaults to `true`.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' showall(STDOUT, df, false, "Row #", true)
function Base.showall(io::IO,
                      adf::AbstractDataFrame,
                      splitchunks::Bool = false,
                      rowlabel::Symbol = symbol("Row #"),
                      displaysummary::Bool = true) # -> Nothing
    rowindices1 = 1:size(adf, 1)
    rowindices2 = 1:0
    maxwidths = getmaxwidths(adf, rowindices1, rowindices2, rowlabel)
    width = getprintedwidth(maxwidths)
    showrows(io,
             adf,
             rowindices1,
             rowindices2,
             maxwidths,
             splitchunks,
             rowlabel,
             displaysummary)
    return
end

#' @exported
#' @description
#'
#' Render all of the rows of an AbstractDataFrame to STDOUT. See
#' `showall` documentation for details.
#'
#' @param adf::AbstractDataFrame An AbstractDataFrame.
#' @param splitchunks::Bool Should the printing of the AbstractDataFrame
#'        be done in chunks? Defaults to `false`.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' showall(df, true)
function Base.showall(adf::AbstractDataFrame,
                      splitchunks::Bool = false) # -> Nothing
    showall(STDOUT, adf, splitchunks)
    return
end

#' @description
#'
#' Render a summary of the column names, column types and column missingness
#' count.
#'
#' @param io::IO The `io` to be rendered to.
#' @param adf::AbstractDataFrame An AbstractDataFrame.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' showcols(df, true)
function showcols(io::IO, adf::AbstractDataFrame) # -> Nothing
    println(io, summary(adf))
    metadata = DataFrame(Name = names(adf),
                         Eltype = eltypes(adf),
                         Missing = colmissing(adf))
    showall(io, metadata, true, symbol("Col #"), false)
    return
end

showcols(adf::AbstractDataFrame) = showcols(STDOUT, adf) # -> Nothing

#' @exported
#' @description
#'
#' Print an AbstractDataFrame to an IO system with an added newline.
#'
#' @param io::IO The `io` system to be rendered to.
#' @param adf::AbstractDataFrame An AbstractDataFrame.
#'
#' @returns o::Nothing A `nothing` value.
#'
#' @examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' print(STDOUT, df)
# TODO: Determine if this method is strictly necessary.
function Base.print(io::IO, adf::AbstractDataFrame)
    show(io, adf)
    print(io, '\n')
end