text:
const PFB{T} = Matrix{T}          # polyphase filter bank

abstract type Filter end
abstract type FIRKernel{T} end

# Single rate FIR kernel
mutable struct FIRStandard{T} <: FIRKernel{T}
    h::Vector{T}
    hLen::Int
end

function FIRStandard(h::Vector)
    h    = reverse(h, dims=1)
    hLen = length(h)
    FIRStandard(h, hLen)
end

# Interpolator FIR kernel
mutable struct FIRInterpolator{T} <: FIRKernel{T}
    pfb::PFB{T}
    interpolation::Int
    Nϕ::Int
    tapsPerϕ::Int
    inputDeficit::Int
    ϕIdx::Int
    hLen::Int
end

function FIRInterpolator(h::Vector, interpolation::Integer)
    pfb           = taps2pfb(h, interpolation)
    tapsPerϕ, Nϕ  = size(pfb)
    interpolation = interpolation
    inputDeficit  = 1
    ϕIdx          = 1
    hLen          = length(h)

    FIRInterpolator(pfb, interpolation, Nϕ, tapsPerϕ, inputDeficit, ϕIdx,   hLen)
end

# Decimator FIR kernel
mutable struct FIRDecimator{T} <: FIRKernel{T}
    h::Vector{T}
    hLen::Int
    decimation::Int
    inputDeficit::Int
end

function FIRDecimator(h::Vector, decimation::Integer)
    h            = reverse(h, dims=1)
    hLen         = length(h)
    inputDeficit = 1
    FIRDecimator(h, hLen, decimation, inputDeficit)
end

# Rational resampler FIR kernel
mutable struct FIRRational{T}  <: FIRKernel{T}
    pfb::PFB{T}
    ratio::Rational{Int}
    Nϕ::Int
    ϕIdxStepSize::Int
    tapsPerϕ::Int
    ϕIdx::Int
    inputDeficit::Int
    hLen::Int
end

function FIRRational(h::Vector, ratio::Rational)
    pfb          = taps2pfb(h, numerator(ratio))
    tapsPerϕ, Nϕ = size(pfb)
    ϕIdxStepSize = mod(denominator(ratio), numerator(ratio))
    ϕIdx         = 1
    inputDeficit = 1
    hLen         = length(h)
    FIRRational(pfb, ratio, Nϕ, ϕIdxStepSize, tapsPerϕ, ϕIdx, inputDeficit, hLen)
end

FIRRational(h::Vector,ratio::Integer)=FIRRational(h,convert(Rational,ratio))

#
# Arbitrary resampler FIR kernel
#
# This kernel is different from the others in that it has two polyphase filtlter banks.
# The the second filter bank, dpfb, is the derivative of pfb. The purpose of this is to
# allow us to compute two y values, yLower & yUpper, whitout having to advance the input
# index by 1. It makes the kernel simple by not having to store extra state in the case
# when where's at the last polphase branch and the last available input sample. By using
# a derivitive filter, we can always compute the output in that scenario.
# See section 7.6.1 in [1] for a better explanation.

mutable struct FIRArbitrary{T} <: FIRKernel{T}
    rate::Float64
    pfb::PFB{T}
    dpfb::PFB{T}
    Nϕ::Int
    tapsPerϕ::Int
    ϕAccumulator::Float64
    ϕIdx::Int
    α::Float64
    Δ::Float64
    inputDeficit::Int
    xIdx::Int
    hLen::Int
end

function FIRArbitrary(h::Vector, rate_in::Real, Nϕ_in::Integer)
    rate         = convert(Float64, rate_in)
    Nϕ           = convert(Int, Nϕ_in)
    dh           = [diff(h); zero(eltype(h))]
    pfb          = taps2pfb(h,  Nϕ)
    dpfb         = taps2pfb(dh, Nϕ)
    tapsPerϕ     = size(pfb, 1)
    ϕAccumulator = 1.0
    ϕIdx         = 1
    α            = 0.0
    Δ            = Nϕ/rate
    inputDeficit = 1
    xIdx         = 1
    hLen         = length(h)
    FIRArbitrary(
        rate,
        pfb,
        dpfb,
        Nϕ,
        tapsPerϕ,
        ϕAccumulator,
        ϕIdx,
        α,
        Δ,
        inputDeficit,
        xIdx,
        hLen
    )
end

# FIRFilter - the kernel does the heavy lifting
mutable struct FIRFilter{Tk<:FIRKernel} <: Filter
    kernel::Tk
    history::Vector
    historyLen::Int
    h::Vector
end

# Constructor for single-rate, decimating, interpolating, and rational resampling filters
"""
    FIRFilter(h[, ratio])

Construct a stateful FIRFilter object from the vector of filter taps `h`.
`ratio` is an optional rational integer which specifies
the input to output sample rate relationship (e.g. `147//160` for
converting recorded audio from 48 KHz to 44.1 KHz).
"""
function FIRFilter(h::Vector, resampleRatio::Union{Integer,Rational} = 1)
    interpolation = numerator(resampleRatio)
    decimation    = denominator(resampleRatio)
    historyLen    = 0

    if resampleRatio == 1                                     # single-rate
        kernel     = FIRStandard(h)
        historyLen = kernel.hLen - 1
    elseif interpolation == 1                                 # decimate
        kernel     = FIRDecimator(h, decimation)
        historyLen = kernel.hLen - 1
    elseif decimation == 1                                    # interpolate
        kernel     = FIRInterpolator(h, interpolation)
        historyLen = kernel.tapsPerϕ - 1
    else                                                      # rational
        kernel     = FIRRational(h, resampleRatio)
        historyLen = kernel.tapsPerϕ - 1
    end

    history = zeros(historyLen)

    FIRFilter(kernel, history, historyLen, h)
end

# Constructor for arbitrary resampling filter (polyphase interpolator w/ intra-phase linear interpolation)
"""
    FIRFilter(h, rate[, Nϕ])

Returns a polyphase FIRFilter object from the vector of filter taps `h`.
`rate` is a floating point number that specifies the input to output
sample-rate relationship ``\\frac{fs_{out}}{fs_{in}}``. `Nϕ` is an
optional parameter which specifies the number of *phases* created from
`h`. `Nϕ` defaults to 32.
"""
function FIRFilter(h::Vector, rate::AbstractFloat, Nϕ::Integer=32)
    rate > 0.0 || error("rate must be greater than 0")
    kernel     = FIRArbitrary(h, rate, Nϕ)
    historyLen = kernel.tapsPerϕ - 1
    history    = zeros(historyLen)
    FIRFilter(kernel, history, historyLen, h)
end

# Constructor for a resampling FIR filter, where the user needs only to set the sampling rate
function FIRFilter(rate::AbstractFloat, Nϕ::Integer=32)
    h = resample_filter(rate, Nϕ)
    FIRFilter(h, rate)
end

function FIRFilter(rate::Union{Integer,Rational})
    h = resample_filter(rate)
    FIRFilter(h, rate)
end

#
# setphase! set's filter kernel phase index
#
function setphase!(kernel::FIRDecimator, ϕ::Real)
    ϕ >= zero(ϕ) || throw(ArgumentError("ϕ must be >= 0"))
    xThrowaway = round(Int, ϕ)
    kernel.inputDeficit += xThrowaway
    nothing
end

function setphase!(kernel::Union{FIRInterpolator, FIRRational}, ϕ::Real)
    ϕ >= zero(ϕ) || throw(ArgumentError("ϕ must be >= 0"))
    (ϕ, xThrowaway) = modf(ϕ)
    kernel.inputDeficit += round(Int, xThrowaway)
    kernel.ϕIdx = round(ϕ*(kernel.Nϕ) + 1.0)
    nothing
end

function setphase!(kernel::FIRArbitrary, ϕ::Real)
    ϕ >= zero(ϕ) || throw(ArgumentError("ϕ must be >= 0"))
    (ϕ, xThrowaway) = modf(ϕ)
    kernel.inputDeficit += round(Int, xThrowaway)
    kernel.ϕAccumulator = ϕ*(kernel.Nϕ) + 1.0
    kernel.ϕIdx         = round(kernel.ϕAccumulator)
    kernel.α            = modf(kernel.ϕAccumulator)[1]
    nothing
end

setphase!(self::FIRFilter, ϕ::Real) = setphase!(self.kernel, ϕ)

#
# reset! filter and its kernel to an initial state
#

# Generic case for FIRInterpolator and FIRStandard
function reset!(kernel::FIRKernel)
    kernel
end

function reset!(kernel::FIRRational)
    kernel.ϕIdx         = 1
    kernel.inputDeficit = 1
    kernel
end

function reset!(kernel::FIRDecimator)
    kernel.inputDeficit = 1
    kernel
end

function reset!(kernel::FIRArbitrary)
    kernel.ϕAccumulator = 1.0
    kernel.ϕIdx         = 1
    kernel.α            = 0.0
    kernel.inputDeficit = 1
    kernel.xIdx         = 1
    kernel
end

# For FIRFilter, set history vector to zeros of same type and required length
function reset!(self::FIRFilter)
    self.history = zeros(eltype(self.history), self.historyLen)
    reset!(self.kernel)
    self
end

#
# taps2pfb
#
# Converts a vector of coefficients to a matrix. Each column is a filter.
# NOTE: also flips the matrix up/down so computing the dot product of a
#       column with a signal vector is more efficient (since filter is convolution)
# Appends zeros if necessary.
# Example:
#   julia> taps2pfb( [1:9], 4 )
#   3x4 Array{Int64,2}:
#    9  0  0  0
#    5  6  7  8
#    1  2  3  4
#
#  In this example, the first phase, or ϕ, is [9, 5, 1].

function taps2pfb(h::Vector{T}, Nϕ::Integer) where T
    hLen     = length(h)
    tapsPerϕ = ceil(Int, hLen/Nϕ)
    pfbSize  = tapsPerϕ * Nϕ
    pfb      = Matrix{T}(undef, tapsPerϕ, Nϕ)
    hIdx     = 1

    for rowIdx in tapsPerϕ:-1:1, colIdx in 1:Nϕ
        tap = hIdx > hLen ? zero(T) : h[hIdx]
        @inbounds pfb[rowIdx,colIdx] = tap
        hIdx += 1
    end

    return pfb
end

#
# Calculates the resulting length of a multirate filtering operation, given a
#   FIRFilter{FIRRational} object and an input vector
#
# ( It's hard to explain how this works without a diagram )
#

function outputlength(inputlength::Integer, ratio::Union{Integer,Rational}, initialϕ::Integer)
    interpolation = numerator(ratio)
    decimation    = denominator(ratio)
    outLen        = ((inputlength * interpolation) - initialϕ + 1) / decimation
    ceil(Int, outLen)
end

function outputlength(kernel::FIRStandard, inputlength::Integer)
    inputlength
end

function outputlength(kernel::FIRInterpolator, inputlength::Integer)
    outputlength(inputlength-kernel.inputDeficit+1, kernel.interpolation, kernel.ϕIdx)
end

function outputlength(kernel::FIRDecimator, inputlength::Integer)
    outputlength(inputlength-kernel.inputDeficit+1, 1//kernel.decimation, 1)
end

function outputlength(kernel::FIRRational, inputlength::Integer)
    outputlength(inputlength-kernel.inputDeficit+1, kernel.ratio, 1)
end

function outputlength(kernel::FIRArbitrary, inputlength::Integer)
    ceil(Int, (inputlength-kernel.inputDeficit+1) * kernel.rate)
end

function outputlength(self::FIRFilter, inputlength::Integer)
    outputlength(self.kernel, inputlength)
end

#
# Calculates the input length of a multirate filtering operation,
# given the output length
#

function inputlength(outputlength::Int, ratio::Union{Integer,Rational}, initialϕ::Integer)
    interpolation = numerator(ratio)
    decimation    = denominator(ratio)
    inLen         = (outputlength * decimation + initialϕ - 1) / interpolation
    floor(Int, inLen)
end

function inputlength(kernel::FIRStandard, outputlength::Integer)
    outputlength
end

function inputlength(kernel::FIRInterpolator, outputlength::Integer)
    inLen = inputlength(outputlength, kernel.interpolation, kernel.ϕIdx)
    inLen += kernel.inputDeficit - 1
end

function inputlength(kernel::FIRDecimator, outputlength::Integer)
    inLen  = inputlength(outputlength, 1//kernel.decimation, 1)
    inLen += kernel.inputDeficit - 1
end

function inputlength(kernel::FIRRational, outputlength::Integer)
    inLen  = inputlength(outputlength, kernel.ratio, kernel.ϕIdx)
    inLen += kernel.inputDeficit - 1

end

# TODO: figure out why this fails. Might be fine, but the filter operation might not being stepping through the phases correcty.
function inputlength(kernel::FIRArbitrary, outputlength::Integer)
    inLen  = floor(Int, outputlength/kernel.rate)
    inLen += kernel.inputDeficit - 1
end

function inputlength(self::FIRFilter, outputlength::Integer)
    inputlength(self.kernel, outputlength)
end

#
# Calculates the delay caused by the FIR filter in # samples, at the input sample rate, caused by the filter process
#

function timedelay(kernel::Union{FIRRational, FIRInterpolator, FIRArbitrary})
    (kernel.hLen - 1)/(2.0*kernel.Nϕ)
end

function timedelay(kernel::Union{FIRStandard, FIRDecimator})
    (kernel.hLen - 1)/2
end

function timedelay(self::FIRFilter)
    timedelay(self.kernel)
end

#
# Single rate filtering
#

function filt!(buffer::AbstractVector{Tb}, self::FIRFilter{FIRStandard{Th}}, x::AbstractVector{Tx}) where {Tb,Th,Tx}
    kernel              = self.kernel
    history::Vector{Tx} = self.history
    bufLen              = length(buffer)
    xLen                = length(x)

    bufLen >= xLen || error("buffer length must be >= x length")

    h = kernel.h
    for i = 1:min(kernel.hLen-1, xLen)
        @inbounds buffer[i] = unsafe_dot(h, history, x, i)
    end
    for i = kernel.hLen:xLen
        @inbounds buffer[i] = unsafe_dot(h, x, i)
    end

    self.history = shiftin!(history, x)

    return xLen
end

function filt(self::FIRFilter{FIRStandard{Th}}, x::AbstractVector{Tx}) where {Th,Tx}
    bufLen         = outputlength(self, length(x))
    buffer         = Vector{promote_type(Th,Tx)}(undef, bufLen)
    samplesWritten = filt!(buffer, self, x)

    samplesWritten == bufLen || resize!(buffer, samplesWritten)

    return buffer
end

#
# Interpolation
#

function filt!(buffer::AbstractVector{Tb}, self::FIRFilter{FIRInterpolator{Th}}, x::AbstractVector{Tx}) where {Tb,Th,Tx}
    kernel              = self.kernel
    history::Vector{Tx} = self.history
    interpolation       = kernel.interpolation
    xLen                = length(x)
    bufLen              = length(buffer)
    bufIdx              = 0

    if xLen < kernel.inputDeficit
        self.history = shiftin!(history, x)
        kernel.inputDeficit -= xLen
        return bufIdx
    end

    inputIdx = kernel.inputDeficit
    bufLen >= outputlength(self, xLen) || error("length(buffer) must be >= interpolation * length(x)")

    while inputIdx <= xLen
        bufIdx += 1

        if inputIdx < kernel.tapsPerϕ
            accumulator = unsafe_dot(kernel.pfb, kernel.ϕIdx, history, x, inputIdx)
        else
            accumulator = unsafe_dot(kernel.pfb, kernel.ϕIdx, x, inputIdx)
        end

        buffer[bufIdx]          = accumulator
        (kernel.ϕIdx, inputIdx) = kernel.ϕIdx == kernel.Nϕ ? (1, inputIdx+1) : (kernel.ϕIdx+1, inputIdx)
    end

    kernel.inputDeficit = 1
    self.history        = shiftin!(history, x)

    return bufIdx
end

function filt(self::FIRFilter{FIRInterpolator{Th}}, x::AbstractVector{Tx}) where {Th,Tx}
    bufLen         = outputlength(self, length(x))
    buffer         = Vector{promote_type(Th,Tx)}(undef, bufLen)
    samplesWritten = filt!(buffer, self, x)

    samplesWritten == bufLen || resize!(buffer, samplesWritten)

    return buffer
end

#
# Rational resampling
#

function filt!(buffer::AbstractVector{Tb}, self::FIRFilter{FIRRational{Th}}, x::AbstractVector{Tx}) where {Tb,Th,Tx}
    kernel              = self.kernel
    history::Vector{Tx} = self.history
    xLen                = length(x)
    bufLen              = length(buffer)
    bufIdx              = 0

    if xLen < kernel.inputDeficit
        self.history = shiftin!(history, x)
        kernel.inputDeficit -= xLen
        return bufIdx
    end

    outLen = outputlength(xLen-kernel.inputDeficit+1, kernel.ratio, kernel.ϕIdx)
    bufLen >= outLen || error("buffer is too small")

    interpolation       = numerator(kernel.ratio)
    decimation          = denominator(kernel.ratio)
    inputIdx            = kernel.inputDeficit

    while inputIdx <= xLen
        bufIdx += 1

        if inputIdx < kernel.tapsPerϕ
            accumulator = unsafe_dot(kernel.pfb, kernel.ϕIdx, history, x, inputIdx)
        else
            accumulator = unsafe_dot(kernel.pfb, kernel.ϕIdx, x, inputIdx)
        end

        buffer[bufIdx]  = accumulator
        inputIdx       += div(kernel.ϕIdx + decimation - 1, interpolation)
        ϕIdx            = kernel.ϕIdx + kernel.ϕIdxStepSize
        kernel.ϕIdx     = ϕIdx > interpolation ? ϕIdx - interpolation : ϕIdx
    end

    kernel.inputDeficit = inputIdx - xLen
    self.history        = shiftin!(history, x)

    return bufIdx
end

function filt(self::FIRFilter{FIRRational{Th}}, x::AbstractVector{Tx}) where {Th,Tx}
    bufLen         = outputlength(self, length(x))
    buffer         = Vector{promote_type(Th,Tx)}(undef, bufLen)
    samplesWritten = filt!(buffer, self, x)

    samplesWritten == bufLen || resize!(buffer, samplesWritten)
    return buffer
end

#
# Decimation
#

function filt!(buffer::AbstractVector{Tb}, self::FIRFilter{FIRDecimator{Th}}, x::AbstractVector{Tx}) where {Tb,Th,Tx}
    kernel              = self.kernel
    bufLen              = length(buffer)
    xLen                = length(x)
    history::Vector{Tx} = self.history
    bufIdx              = 0

    if xLen < kernel.inputDeficit
        self.history = shiftin!(history, x)
        kernel.inputDeficit -= xLen
        return bufIdx
    end

    outLen              = outputlength(self, xLen)
    inputIdx            = kernel.inputDeficit

    nbufout = fld(xLen - inputIdx, kernel.decimation) + 1
    bufLen >= nbufout || throw(ArgumentError("buffer length insufficient"))

    while inputIdx <= xLen
        bufIdx += 1

        if inputIdx < kernel.hLen
            accumulator = unsafe_dot(kernel.h, history, x, inputIdx)
        else
            accumulator = unsafe_dot(kernel.h, x, inputIdx)
        end

        @inbounds buffer[bufIdx] = accumulator
        inputIdx                += kernel.decimation
    end

    kernel.inputDeficit = inputIdx - xLen
    self.history        = shiftin!(history, x)

    return bufIdx
end

function filt(self::FIRFilter{FIRDecimator{Th}}, x::AbstractVector{Tx}) where {Th,Tx}
    bufLen         = outputlength(self, length(x))
    buffer         = Vector{promote_type(Th,Tx)}(undef, bufLen)
    samplesWritten = filt!(buffer, self, x)

    samplesWritten == bufLen || resize!(buffer, samplesWritten)

    return buffer
end

#
# Arbitrary resampling
#
# Updates FIRArbitrary state. See Section 7.5.1 in [1].
#   [1] uses a phase accumilator that increments by Δ (Nϕ/rate)

function update(kernel::FIRArbitrary)
    kernel.ϕAccumulator += kernel.Δ

    if kernel.ϕAccumulator > kernel.Nϕ
        kernel.xIdx        += div(kernel.ϕAccumulator-1, kernel.Nϕ)
        kernel.ϕAccumulator = mod(kernel.ϕAccumulator-1, kernel.Nϕ) + 1
    end

    kernel.ϕIdx = floor(Int, kernel.ϕAccumulator)
    kernel.α    = kernel.ϕAccumulator - kernel.ϕIdx
end

function filt!(
    buffer::AbstractVector{Tb},
    self::FIRFilter{FIRArbitrary{Th}},
    x::AbstractVector{Tx}
) where {Tb,Th,Tx}
    kernel              = self.kernel
    pfb                 = kernel.pfb
    dpfb                = kernel.dpfb
    xLen                = length(x)
    bufIdx              = 0
    history::Vector{Tx} = self.history

    # Do we have enough input samples to produce one or more output samples?
    if xLen < kernel.inputDeficit
        self.history = shiftin!(history, x)
        kernel.inputDeficit -= xLen
        return bufIdx
    end

    # Skip over input samples that are not needed to produce output results.
    # We do this by seting inputIdx to inputDeficit which was calculated in the previous run.
    # InputDeficit is set to 1 when instantiation the FIRArbitrary kernel, that way the first
    #   input always produces an output.
    kernel.xIdx = kernel.inputDeficit

    while kernel.xIdx <= xLen
        bufIdx += 1

        if kernel.xIdx < kernel.tapsPerϕ
            yLower = unsafe_dot(pfb,  kernel.ϕIdx, history, x, kernel.xIdx)
            yUpper = unsafe_dot(dpfb, kernel.ϕIdx, history, x, kernel.xIdx)
        else
            yLower = unsafe_dot(pfb,  kernel.ϕIdx, x, kernel.xIdx)
            yUpper = unsafe_dot(dpfb, kernel.ϕIdx, x, kernel.xIdx)
        end

        # Used to have @inbounds. Restore @inbounds if buffer length
        # can be verified prior to access.
        buffer[bufIdx] = yLower + yUpper * kernel.α
        update(kernel)
    end

    kernel.inputDeficit = kernel.xIdx - xLen
    self.history        = shiftin!(history, x)

    return bufIdx
end

function filt(self::FIRFilter{FIRArbitrary{Th}}, x::AbstractVector{Tx}) where {Th,Tx}
    bufLen         = outputlength(self, length(x))
    buffer         = Vector{promote_type(Th,Tx)}(undef, bufLen)
    samplesWritten = filt!(buffer, self, x)

    samplesWritten == bufLen || resize!(buffer, samplesWritten)

    return buffer
end

#
# Stateless filt implementations
#

# Single-rate, decimation, interpolation, and rational resampling.
function filt(h::Vector, x::AbstractVector, ratio::Union{Integer,Rational})
    self = FIRFilter(h, ratio)
    filt(self, x)
end

# Arbitrary resampling with polyphase interpolation and two neighbor linear interpolation.
function filt(h::Vector, x::AbstractVector, rate::AbstractFloat, Nϕ::Integer=32)
    self = FIRFilter(h, rate, Nϕ)
    filt(self, x)
end

"""
    resample(x, rate[, coef])

Resample `x` at rational or arbitrary `rate`.
`coef` is an optional vector of FIR filter taps. If `coef`
is not provided, the taps will be computed using a Kaiser window.

Internally, `resample` uses a polyphase `FIRFilter` object,
but performs additional operations to make resampling a signal easier.
It compensates for the `FIRFilter`'s delay (ramp-up), and appends
zeros to `x`. The result is that when the input and output signals
are plotted on top of each other, they correlate very well, but one
signal will have more samples that the other.
"""
function resample(x::AbstractVector, rate::Real, h::Vector)
    self = FIRFilter(h, rate)

    # Get delay, in # of samples at the output rate, caused by filtering processes
    τ = timedelay(self)

    # Use setphase! to
    #   a) adjust the input samples to skip over before producing and output (integer part of τ)
    #   b) set the ϕ index of the PFB (fractional part of τ)
    setphase!(self, τ)

    # Calculate the number of 0's required
    outLen      = ceil(Int, length(x)*rate)
    reqInlen    = inputlength(self, outLen)
    reqZerosLen = reqInlen - length(x)
    xPadded     = [x; zeros(eltype(x), reqZerosLen)]

    filt(self, xPadded)
end

function resample(x::AbstractVector, rate::Real)
    h = resample_filter(rate)
    resample(x, rate, h)
end

"""
    resample(x::AbstractMatrix, rate::Real, h::Vector = resample_filter(rate); dims)

Resample a matrix `x` along dimension `dims`.
"""
function resample(x::AbstractMatrix, rate::Real, h::Vector = resample_filter(rate); dims)
    mapslices(x; dims=dims) do x
        resample(x, rate, h)
    end
end

#
# References
#

# [1] F.J. Harris, *Multirate Signal Processing for Communication Systems*. Prentice Hall, 2004
# [2] Dick, C.; Harris, F., "Options for arbitrary resamplers in FPGA-based modulators," Signals, Systems and Computers, 2004. Conference Record of the Thirty-Eighth Asilomar Conference on , vol.1, no., pp.777,781 Vol.1, 7-10 Nov. 2004
# [3] Kim, S.C.; Plishker, W.L.; Bhattacharyya, S.S., "An efficient GPU implementation of an arbitrary resampling polyphase channelizer," Design and Architectures for Signal and Image Processing (DASIP), 2013 Conference on, vol., no., pp.231,238, 8-10 Oct. 2013
# [4] Horridge, J.P.; Frazer, Gordon J., "Accurate arbitrary resampling with exact delay for radar applications," Radar, 2008 International Conference on , vol., no., pp.123,127, 2-5 Sept. 2008
# [5] Blok, M., "Fractional delay filter design for sample rate conversion," Computer Science and Information Systems (FedCSIS), 2012 Federated Conference on , vol., no., pp.701,706, 9-12 Sept. 2012


text:
using Random
using Test
using FFTW
using CUDA

Random.seed!(0)

FT = Float64
ε = eps(FT)

for N in (4, 5)
    A_cpu = rand(FT, (N, N, N))
    A_gpu = CuArray(A_cpu)

    B_cpu = FFTW.fft(A_cpu, 2)

    B_gpu = permutedims(A_gpu, (2, 1, 3))
    B_gpu = CUDA.CUFFT.fft(B_gpu, 1)
    B_gpu = permutedims(B_gpu, (2, 1, 3))

    @test isapprox.(B_cpu, Array(B_gpu), atol=4ε) |> all

    C_cpu = FFTW.ifft(B_cpu, 2)

    C_gpu = permutedims(B_gpu, (2, 1, 3))
    C_gpu = CUDA.CUFFT.ifft(C_gpu, 1)
    C_gpu = permutedims(C_gpu, (2, 1, 3))

    @test isapprox.(C_cpu, Array(C_gpu), atol=4ε) |> all

    @test isapprox.(C_cpu, A_cpu, atol=ε) |> all
    @test isapprox.(Array(C_gpu), Array(A_gpu), atol=ε) |> all
end

ext: jl


text:
# logistic regression on mnist data

using MultiStochGrad

using  LinearAlgebra, Test
using Random, Distributions

include("../src/applis/LinearRegression.jl")

# the . means we use the module LinearRegressionMod in the current module, other wise refer toa package!

function test_linear_regression_scsg()
    true_coefficients = [13.37, -4.2, 3.14]
    # sample noise according to a N(0,1)
    datas = Vector{Vector{Float64}}()
    values = Vector{Float64}()
    d = Normal()
    for i in 1:100
        v = Vector{Float64}()
        push!(v,1.)
        push!(v, rand())
        push!(v, rand())
        y = dot(true_coefficients, v) + rand(d)
        push!(datas, v)
        push!(values, y)
    end
    observations= Observations(datas,values)
    nbterms = length(observations.datas)
    # define our evaluations
    linear_reg_term_value = make_linear_reg_term_value(observations)
    term_function =  TermFunction{typeof(linear_reg_term_value)}(linear_reg_term_value, nbterms, Dims{1}(3))
    #
    linear_reg_term_gradient = make_linear_reg_term_gradient(observations)
    term_gradient = TermGradient{typeof(linear_reg_term_gradient)}(linear_reg_term_gradient ,nbterms, Dims{1}(3))
    evaluator = Evaluator{typeof(linear_reg_term_value),typeof(linear_reg_term_gradient)}(term_function,term_gradient)
    #
    scsg_pb = SCSG(0.1, 0.3, 1 , 0.7)
    # solve.
    nb_iter = 50
    position = fill(1., 3)
    initial_error = compute_value(evaluator, position)
    @info "\n initial error \n\n" initial_error
    position, value = minimize(scsg_pb, evaluator, nb_iter, position)
    @info "value , position " value  position
    @info "\n\n test_linear_regression_scsg done \n \n"
    value < 0.65 ? true : false
end

function test_linear_regression_svrg()
    true_coefficients = [13.37, -4.2, 3.14]
    # sample noise according to a N(0,1)
    datas = Vector{Vector{Float64}}()
    values = Vector{Float64}()
    d = Normal()
    for i in 1:100
        v = Vector{Float64}()
        push!(v,1.)
        push!(v, rand())
        push!(v, rand())
        y = dot(true_coefficients, v) + rand(d)
        push!(datas, v)
        push!(values, y)
    end
    observations= Observations(datas,values)
    nbterms = length(observations.datas)
    # define our evaluations
    linear_reg_term_value = make_linear_reg_term_value(observations)
    term_function =  TermFunction{typeof(linear_reg_term_value)}(linear_reg_term_value, nbterms, Dims{1}(3))
    #
    linear_reg_term_gradient = make_linear_reg_term_gradient(observations)
    term_gradient = TermGradient{typeof(linear_reg_term_gradient)}(linear_reg_term_gradient ,nbterms, Dims{1}(3))
    evaluator = Evaluator{typeof(linear_reg_term_value),typeof(linear_reg_term_gradient)}(term_function,term_gradient)
    #
    svrg_pb = SVRG(100, 0.2)
    # solve.
    nb_iter = 60
    position = fill(1., 3)
    initial_error = compute_value(evaluator, position)
    @info "\n initial error \n\n" initial_error
    position, value = minimize(svrg_pb, evaluator, nb_iter, position)
    @info "value , position ", value ,  position
    @info "\n\n test_linear_regression_svrg done \n \n"
    value < 0.65 ? true : false
end


text:
export overlap, overlaps, nucAttraction, nucAttractions, elecKinetic, elecKinetics,
       coreH, coreHij

@inline function oneBodyBFTensorCore(libcinFunc::Symbol,
                                     bf1::FloatingGTBasisFuncs, bf2::FloatingGTBasisFuncs,
                                     nuclei::Vector{String},
                                     nucleiCoords::Vector{<:AbstractArray};
                                     isGradient::Bool=false)
    env = Float64[]
    atm = Int32[]
    bas = Int32[]
    #=
    Add actual nuclei before bases can allow
        natm = length(nuclei)
    to save calculating overhead.
    Otherwise set
        natm = length(nuclei) + length(BasisSet)
    to include ghost atoms.
    =#
    addToDataChain!(env, atm, nuclei, nucleiCoords)
    addToDataChain!(env, atm, bas, bf1)
    addToDataChain!(env, atm, bas, bf2)

    cintFunc!(Val(libcinFunc),
              zeros(basisSize([bf1.subshell, bf2.subshell])..., 1+isGradient*2),
              [0,1], atm, length(nuclei), bas, 2, env)
end

"""

    oneBodyBFTensor(libcinFunc::Symbol, b1::AbstractGTBasisFuncs,
                    b2::AbstractGTBasisFuncs, nuclei::Array{String, 1}=String[],
                    nucleiCoords::Array{<:AbstractArray, 1}=Array[];
                    isGradient::Bool=false) ->
    Array{Float64, 3}

Core function for one-electron integrals.

`libcinFunc::Symbol` specifies the backend [libcint](https://github.com/sunqm/libcint)
function name, e.g. `"int1e_nuc_cart"` should be converted to `:int1e_nuc_cart` as the
input argument. If the integral does not need the information of nuclei and their
coordinates, those 2 arguments can be omitted. If the integral is a spacial gradient,
`isGradient` should be set to `true`.
"""
@inline function oneBodyBFTensor(libcinFunc::Symbol,
                                 b1::AbstractGTBasisFuncs,
                                 b2::AbstractGTBasisFuncs,
                                 nuclei::Vector{String}=String[],
                                 nucleiCoords::Vector{<:AbstractArray}=Array[];
                                 isGradient::Bool=false)
    f = @inline function (i,j)
        ints = oneBodyBFTensorCore(libcinFunc, i, j, nuclei, nucleiCoords; isGradient)
        ints[ijkIndex(i), ijkIndex(j), :]
    end
    sum([f(i,j) for i in unpackBasisFuncs(b1), j in unpackBasisFuncs(b2)])
end

@inline function oneBodyBSTensor(BasisSet::Vector{<:AbstractGTBasisFuncs},
                                 intFunc::F) where {F<:Function}
    subSize = basisSize(BasisSet) |> collect
    accuSize = vcat(0, accumulate(+, subSize))
    len = subSize |> sum
    nPage = (intFunc(BasisSet[1], BasisSet[1]) |> size)[3]
    tensor = cat(fill(-2*ones(len, len), nPage)..., dims=3)
    for i = 1:length(BasisSet), j = 1:i
        intTensor =  intFunc(BasisSet[i], BasisSet[j])
        for page = 1:nPage
            rowRange = accuSize[i]+1 : accuSize[i+1]
            colRange = accuSize[j]+1 : accuSize[j+1]
            int = intTensor[:, :, page]
            tensor[rowRange, colRange, page] = int
            tensor[colRange, rowRange, page] = int |> transpose
        end
    end
    tensor
end

@inline overlapCore(bf1::AbstractGTBasisFuncs, bf2::AbstractGTBasisFuncs) =
        oneBodyBFTensor(:cint1e_ovlp_cart, bf1, bf2)

"""

    overlap(fb1::AbstractGTBasisFuncs, fb2::AbstractGTBasisFuncs) ->
    Array{Float64, 2}

Return the orbital overlap matrix (an N×N `Matrix` where N is the number of spatial
orbitals) given 2 basis functions.
"""
overlap(bf1::AbstractGTBasisFuncs, bf2::AbstractGTBasisFuncs) =
dropdims(overlapCore(bf1, bf2), dims=3)

@inline overlapsCore(BSet::Vector{<:AbstractGTBasisFuncs}) =
        oneBodyBSTensor(BSet, overlapCore)

"""

    overlaps(BSet::Array{<:AbstractGTBasisFuncs, 1}) -> Array{Float64, 2}

Return the orbital overlap matrix (an N×N `Matrix` where N is the number of spatial
orbitals) given a basis set in the form of an `Array`.
"""
overlaps(BSet::Vector{<:AbstractGTBasisFuncs}) = dropdims(overlapsCore(BSet), dims=3)

@inline nucAttractionCore(bf1::AbstractGTBasisFuncs,
                          bf2::AbstractGTBasisFuncs,
                          nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) =
        oneBodyBFTensor(:cint1e_nuc_cart, bf1, bf2, nuc, nucCoords)

"""

    nucAttraction(fb1::AbstractGTBasisFuncs, fb2::AbstractGTBasisFuncs,
                  nuc::Array{String, 1}, nucCoords::Array{<:AbstractArray, 1}) ->
    Array{Float64, 2}

Return the nuclear attraction matrix (an N×N `Matrix` where N is the number of spatial
orbitals) given 2 basis functions, and the nuclei with their coordinates (in atomic unit).
"""
nucAttraction(bf1::AbstractGTBasisFuncs, bf2::AbstractGTBasisFuncs,
              nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) =
dropdims(nucAttractionCore(bf1, bf2, nuc, nucCoords), dims=3)

@inline nucAttractionsCore(BSet::Vector{<:AbstractGTBasisFuncs},
                           nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) =
        oneBodyBSTensor(BSet, @inline (bf1, bf2)->oneBodyBFTensor(:cint1e_nuc_cart,
                                                                  bf1, bf2, nuc, nucCoords))

"""

    nucAttractions(BSet::Array{<:AbstractGTBasisFuncs, 1}, nuc::Array{String, 1},
                   nucCoords::Array{<:AbstractArray, 1}) ->
    Array{Float64, 2}

Return the nuclear attraction matrix (an N×N `Matrix` where N is the number of spatial
orbitals) given a basis set in the form of an `Array`, and the nuclei with their
coordinates (in atomic unit).
"""
nucAttractions(BSet::Vector{<:AbstractGTBasisFuncs},
               nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) =
dropdims(nucAttractionsCore(BSet, nuc, nucCoords), dims=3)

@inline elecKineticCore(bf1::AbstractGTBasisFuncs,
                        bf2::AbstractGTBasisFuncs) =
        oneBodyBFTensor(:cint1e_kin_cart, bf1, bf2)

"""

    elecKinetic(fb1::AbstractGTBasisFuncs, fb2::AbstractGTBasisFuncs) ->
    Array{Float64, 2}

Return the electron kinetic energy matrix (an N×N `Matrix` where N is the number of spatial
orbitals) given 2 basis functions.
"""
elecKinetic(bf1::AbstractGTBasisFuncs, bf2::AbstractGTBasisFuncs) =
dropdims(elecKineticCore(bf1, bf2), dims=3)

@inline elecKineticsCore(BSet::Vector{<:AbstractGTBasisFuncs}) =
                     oneBodyBSTensor(BSet, elecKineticCore)

"""

    elecKinetics(BSet::Array{<:AbstractGTBasisFuncs, 1}) -> Array{Float64, 2}

Return the electron kinetic energy matrix (an N×N `Matrix` where N is the number of spatial
orbitals) given a basis set in the form of an `Array`.
"""
elecKinetics(BSet::Vector{<:AbstractGTBasisFuncs}) =
dropdims(elecKineticsCore(BSet), dims=3)

@inline coreHijCore(bf1::AbstractGTBasisFuncs, bf2::AbstractGTBasisFuncs,
                    nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) =
        elecKineticCore(bf1, bf2) + nucAttractionCore(bf1, bf2, nuc, nucCoords)

"""

    coreHij(fb1::AbstractGTBasisFuncs, fb2::AbstractGTBasisFuncs, nuc::Array{String, 1},
            nucCoords::Array{<:AbstractArray, 1}) ->
    Array{Float64, 2}

Return a matrix element or block of the core Hamiltonian (an N×N `Matrix` where N is the
number of spatial orbitals) given 2 basis functions.
"""
coreHij(bf1::AbstractGTBasisFuncs, bf2::AbstractGTBasisFuncs,
        nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) =
dropdims(coreHijCore(bf1, bf2, nuc, nucCoords), dims=3)

@inline coreHCore(BSet::Vector{<:AbstractGTBasisFuncs},
                  nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) =
        elecKineticsCore(BSet) + nucAttractionsCore(BSet, nuc, nucCoords)

"""

    coreH(BSet::Array{<:AbstractGTBasisFuncs, 1}, nuc::Array{String, 1},
          nucCoords::Array{<:AbstractArray, 1}) -> Array{Float64, 2}

Return the core Hamiltonian matrix (an N×N `Matrix` where N is the number of spatial
orbitals) given a basis set in the form of an `Array`.
"""
coreH(BSet::Vector{<:AbstractGTBasisFuncs},
      nuc::Vector{String}, nucCoords::Vector{<:AbstractArray}) =
dropdims(coreHCore(BSet, nuc, nucCoords), dims=3)


text:
using Luxor

struct Offset
    x::Int16
    y::Int16
end

Base.:+(o::Offset, p::Luxor.Point) = p + o
Base.:+(p::Luxor.Point, o::Offset) = p + (o.x, o.y)


text:
abstract type AbstractModelContainer end

mutable struct OptimizationContainer <: AbstractModelContainer
    JuMPmodel::JuMP.Model
    time_steps::UnitRange{Int}
    resolution::Dates.TimePeriod
    settings::Settings
    settings_copy::Settings
    variables::Dict{Symbol, AbstractArray}
    aux_variables::Dict{AuxVarKey, AbstractArray}
    constraints::Dict{Symbol, AbstractArray}
    cost_function::JuMP.AbstractJuMPScalar
    expressions::DenseAxisArrayContainer
    parameters::ParametersContainer
    initial_conditions::InitialConditions
    pm::Union{Nothing, PM.AbstractPowerModel}
    base_power::Float64
    solve_timed_log::Dict{Symbol, Any}

    function OptimizationContainer(
        sys::PSY.System,
        settings::Settings,
        jump_model::Union{Nothing, JuMP.Model},
    )
        resolution = PSY.get_time_series_resolution(sys)
        resolution = IS.time_period_conversion(resolution)
        use_parameters = get_use_parameters(settings)

        new(
            jump_model === nothing ? _make_jump_model(settings) :
            _prepare_external_jump_model!(jump_model, settings),
            1:1,
            resolution,
            settings,
            copy_for_serialization(settings),
            Dict{Symbol, AbstractArray}(),
            Dict{Symbol, AbstractArray}(),
            Dict{Symbol, AbstractArray}(),
            zero(JuMP.GenericAffExpr{Float64, JuMP.VariableRef}),
            DenseAxisArrayContainer(),
            ParametersContainer(),
            InitialConditions(use_parameters = use_parameters),
            nothing,
            PSY.get_base_power(sys),
            Dict{Symbol, Any}(),
        )
    end
end

function _validate_warm_start_support(JuMPmodel::JuMP.Model, warm_start_enabled::Bool)
    !warm_start_enabled && return warm_start_enabled
    solver_supports_warm_start =
        MOI.supports(JuMP.backend(JuMPmodel), MOI.VariablePrimalStart(), MOI.VariableIndex)
    if !solver_supports_warm_start
        solver_name = JuMP.solver_name(JuMPmodel)
        @warn("$(solver_name) does not support warm start")
    end
    return solver_supports_warm_start
end

function _finalize_jump_model!(JuMPmodel::JuMP.Model, settings::Settings)
    warm_start_enabled = get_warm_start(settings)
    solver_supports_warm_start = _validate_warm_start_support(JuMPmodel, warm_start_enabled)
    set_warm_start!(settings, solver_supports_warm_start)

    if get_optimizer_log_print(settings)
        JuMP.unset_silent(JuMPmodel)
        @debug "optimizer unset to silent"
    else
        JuMP.set_silent(JuMPmodel)
        @debug "optimizer set to silent"
    end
end

function _prepare_external_jump_model!(JuMPmodel::JuMP.Model, settings::Settings)
    parameters = get_use_parameters(settings)
    optimizer = get_optimizer(settings)
    if get_direct_mode_optimizer(settings)
        throw(
            IS.ConflictingInputsError(
                "Externally provided JuMP models are not compatible with the direct model keyword argument. Use JuMP.direct_model before passing the custom model",
            ),
        )
    end

    if parameters
        if !haskey(JuMPmodel.ext, :ParameterJuMP)
            @info("Model doesn't have Parameters enabled. Parameters will be enabled")
            PJ.enable_parameters(JuMPmodel)
            JuMP.set_optimizer(JuMPmodel, optimizer)
        end
    end
    _finalize_jump_model!(JuMPmodel, settings)
    return JuMPmodel
end

function _make_jump_model(settings::Settings)
    @debug "Instantiating the JuMP model"
    parameters = get_use_parameters(settings)
    optimizer = get_optimizer(settings)
    if get_direct_mode_optimizer(settings)
        JuMPmodel = JuMP.direct_model(MOI.instantiate(optimizer))
    elseif optimizer === nothing
        JuMPmodel = JuMP.Model()
        @debug "The optimization model has no optimizer attached"
    else
        JuMPmodel = JuMP.Model(optimizer)
    end
    parameters && PJ.enable_parameters(JuMPmodel)
    _finalize_jump_model!(JuMPmodel, settings)

    return JuMPmodel
end

function _make_container_array(parameters::Bool, ax...)
    if parameters
        return JuMP.Containers.DenseAxisArray{PGAE}(undef, ax...)
    else
        return JuMP.Containers.DenseAxisArray{GAE}(undef, ax...)
    end
    return
end

function _make_expressions_dict!(
    optimization_container::OptimizationContainer,
    bus_numbers::Vector{Int},
    ::Type{<:PM.AbstractPowerModel},
)
    settings = optimization_container.settings
    parameters = get_use_parameters(settings)
    time_steps = 1:get_horizon(settings)
    optimization_container.expressions = DenseAxisArrayContainer(
        :nodal_balance_active =>
            _make_container_array(parameters, bus_numbers, time_steps),
        :nodal_balance_reactive =>
            _make_container_array(parameters, bus_numbers, time_steps),
    )
    return
end

function _make_expressions_dict!(
    optimization_container::OptimizationContainer,
    bus_numbers::Vector{Int},
    ::Type{<:PM.AbstractActivePowerModel},
)
    settings = optimization_container.settings
    parameters = get_use_parameters(settings)
    time_steps = 1:get_horizon(settings)
    optimization_container.expressions = DenseAxisArrayContainer(
        :nodal_balance_active =>
            _make_container_array(parameters, bus_numbers, time_steps),
    )
    return
end

function optimization_container_init!(
    optimization_container::OptimizationContainer,
    ::Type{T},
    sys::PSY.System,
) where {T <: PM.AbstractPowerModel}
    @assert !(optimization_container.JuMPmodel === nothing)
    PSY.set_units_base_system!(sys, "SYSTEM_BASE")
    # The order of operations matter
    settings = get_settings(optimization_container)
    use_parameters = get_use_parameters(settings)
    use_forecasts = get_use_forecast_data(settings)

    if use_parameters
        if !use_forecasts
            throw(
                IS.ConflictingInputsError(
                    "enabling parameters without forecasts is not supported",
                ),
            )
        end
        set_use_parameters!(get_initial_conditions(optimization_container), use_parameters)
    end

    if get_initial_time(settings) == UNSET_INI_TIME
        set_initial_time!(settings, PSY.get_forecast_initial_timestamp(sys))
    end

    if use_forecasts
        if get_horizon(settings) == UNSET_HORIZON
            set_horizon!(settings, PSY.get_forecast_horizon(sys))
        end
        total_number_of_devices = length(get_available_components(PSY.Device, sys))
        optimization_container.time_steps = 1:get_horizon(settings)
        # The 10e6 limit is based on the sizes of the lp benchmark problems http://plato.asu.edu/ftp/lpcom.html
        # The maximum numbers of constraints and variables in the benchmark problems is 1,918,399 and 1,259,121,
        # respectively. See also https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2013/138847.pdf
        variable_count_estimate =
            length(optimization_container.time_steps) * total_number_of_devices
        if variable_count_estimate > 10e6
            @warn(
                "The estimated total number of variables that will be created in the model is $(variable_count_estimate). The total number of variables might be larger than 10e6 and could lead to large build or solve times."
            )
        end
    else
        set_horizon!(settings, 1)
    end

    bus_numbers = sort([PSY.get_number(b) for b in PSY.get_components(PSY.Bus, sys)])
    _make_expressions_dict!(optimization_container, bus_numbers, T)
    return
end

function has_initial_conditions(optimization_container::OptimizationContainer, key::ICKey)
    return has_initial_conditions(optimization_container.initial_conditions, key)
end

function iterate_initial_conditions(optimization_container::OptimizationContainer)
    return iterate_initial_conditions(optimization_container.initial_conditions)
end

function get_initial_conditions(
    optimization_container::OptimizationContainer,
    ::Type{T},
    ::Type{D},
) where {T <: InitialConditionType, D <: PSY.Device}
    return get_initial_conditions(optimization_container, ICKey(T, D))
end

function get_initial_conditions(optimization_container::OptimizationContainer, key::ICKey)
    return get_initial_conditions(optimization_container.initial_conditions, key)
end

function set_initial_conditions!(
    optimization_container::OptimizationContainer,
    key::ICKey,
    value,
)
    set_initial_conditions!(optimization_container.initial_conditions, key, value)
end

_variable_type(cm::OptimizationContainer) = JuMP.variable_type(cm.JuMPmodel)
model_time_steps(optimization_container::OptimizationContainer) =
    optimization_container.time_steps
model_resolution(optimization_container::OptimizationContainer) =
    optimization_container.resolution
model_has_parameters(optimization_container::OptimizationContainer) =
    get_use_parameters(optimization_container.settings)
model_uses_forecasts(optimization_container::OptimizationContainer) =
    get_use_forecast_data(optimization_container.settings)
model_initial_time(optimization_container::OptimizationContainer) =
    get_initial_time(optimization_container.settings)
# Internal Variables, Constraints and Parameters accessors
get_variables(optimization_container::OptimizationContainer) =
    optimization_container.variables
get_aux_variables(optimization_container::OptimizationContainer) =
    optimization_container.aux_variables
get_constraints(optimization_container::OptimizationContainer) =
    optimization_container.constraints
get_parameters(optimization_container::OptimizationContainer) =
    optimization_container.parameters
get_expression(optimization_container::OptimizationContainer, name::Symbol) =
    optimization_container.expressions[name]
get_initial_conditions(optimization_container::OptimizationContainer) =
    optimization_container.initial_conditions
get_PTDF(optimization_container::OptimizationContainer) =
    get_PTDF(optimization_container.settings)
get_settings(optimization_container::OptimizationContainer) =
    optimization_container.settings
get_jump_model(optimization_container::OptimizationContainer) =
    optimization_container.JuMPmodel
get_base_power(optimization_container::OptimizationContainer) =
    optimization_container.base_power

function get_variable(
    optimization_container::OptimizationContainer,
    var_type::AbstractString,
    ::Type{T},
) where {T <: PSY.Component}
    return get_variable(optimization_container, make_variable_name(var_type, T))
end

function get_variable(
    optimization_container::OptimizationContainer,
    ::Type{T},
    ::Type{U},
) where {T <: VariableType, U <: PSY.Component}
    return get_variable(optimization_container, make_variable_name(T, U))
end

function get_variable(
    optimization_container::OptimizationContainer,
    var_type::AbstractString,
)
    return get_variable(optimization_container, make_variable_name(var_type))
end

function get_variable(optimization_container::OptimizationContainer, update_ref::UpdateRef)
    return get_variable(optimization_container, update_ref.access_ref)
end

function get_variable(optimization_container::OptimizationContainer, name::Symbol)
    var = get(optimization_container.variables, name, nothing)
    if var === nothing
        @error "$name is not stored" sort!(get_variable_names(optimization_container))
        throw(IS.InvalidValue("variable $name is not stored"))
    end

    return var
end

function get_variable_names(optimization_container::OptimizationContainer)
    return collect(keys(optimization_container.variables))
end

function assign_variable!(
    optimization_container::OptimizationContainer,
    variable_type::AbstractString,
    ::Type{T},
    value,
) where {T <: PSY.Component}
    assign_variable!(optimization_container, make_variable_name(variable_type, T), value)
    return
end

function assign_variable!(
    optimization_container::OptimizationContainer,
    variable_type::AbstractString,
    value,
)
    assign_variable!(optimization_container, make_variable_name(variable_type), value)
    return
end

function _assign_container!(container::Dict, key, value)
    if haskey(container, key)
        @error "variable $key is already stored" sort!(collect(keys!(container)))
        throw(IS.InvalidValue("$key is already stored"))
    end
    container[key] = value
end

function assign_variable!(
    optimization_container::OptimizationContainer,
    name::Symbol,
    value,
)
    @debug "assign_variable" name
    _assign_container!(optimization_container.variables, name, value)
    return
end

function add_aux_var_container!(
    optimization_container::OptimizationContainer,
    var_key::AuxVarKey{<:AuxVariableType, <:PSY.Component},
    axs...;
    sparse = false,
)
    if sparse
        container = sparse_container_spec(Float64, axs...)
    else
        container = container_spec(Float64, axs...)
    end
    _assign_container!(optimization_container.aux_variables, var_key, container)
    return container
end

# TODO: Use this type of interface for regular vars later
function add_var_container!(
    optimization_container::OptimizationContainer,
    var_name::Symbol,
    axs...;
    sparse = false,
)
    if sparse
        container = sparse_container_spec(JuMP.VariableRef, axs...)
    else
        container = container_spec(JuMP.VariableRef, axs...)
    end
    assign_variable!(optimization_container, var_name, container)
    return container
end

function get_constraint(
    optimization_container::OptimizationContainer,
    constraint_type::AbstractString,
    ::Type{T},
) where {T <: PSY.Component}
    return get_constraint(optimization_container, make_constraint_name(constraint_type, T))
end

function get_constraint(
    optimization_container::OptimizationContainer,
    constraint_type::AbstractString,
)
    return get_constraint(optimization_container, make_constraint_name(constraint_type))
end

function get_constraint(optimization_container::OptimizationContainer, name::Symbol)
    var = get(optimization_container.constraints, name, nothing)
    if var === nothing
        @error "$name is not stored" sort!(get_constraint_names(optimization_container))
        throw(IS.InvalidValue("constraint $name is not stored"))
    end

    return var
end

function get_constraint_names(optimization_container::OptimizationContainer)
    return collect(keys(optimization_container.constraints))
end

function assign_constraint!(
    optimization_container::OptimizationContainer,
    constraint_type::AbstractString,
    ::Type{T},
    value,
) where {T <: PSY.Component}
    assign_constraint!(
        optimization_container,
        make_constraint_name(constraint_type, T),
        value,
    )
    return
end

function assign_constraint!(
    optimization_container::OptimizationContainer,
    constraint_type::AbstractString,
    value,
)
    assign_constraint!(optimization_container, make_constraint_name(constraint_type), value)
    return
end

function assign_constraint!(
    optimization_container::OptimizationContainer,
    name::Symbol,
    value,
)
    @debug "set_constraint" name
    optimization_container.constraints[name] = value
    return
end

function add_cons_container!(
    optimization_container::OptimizationContainer,
    cons_name::Symbol,
    axs...;
    sparse = false,
)
    if !haskey(optimization_container.constraints, cons_name)
        if sparse
            container = sparse_container_spec(JuMP.ConstraintRef, axs...)
        else
            container = container_spec(JuMP.ConstraintRef, axs...)
        end
        assign_constraint!(optimization_container, cons_name, container)
    else
        Throw(error("constraint container already exists: $cons_name"))
    end
    return container
end

function get_parameter_names(optimization_container::OptimizationContainer)
    return collect(keys(optimization_container.parameters))
end

function get_parameter_container(
    optimization_container::OptimizationContainer,
    name::AbstractString,
)
    return get_parameter_container(optimization_container, Symbol(name))
end

function get_parameter_container(
    optimization_container::OptimizationContainer,
    name::Symbol,
)
    container = get(optimization_container.parameters, name, nothing)
    if container === nothing
        @error "$name is not stored" sort!(get_parameter_names(optimization_container))
        throw(IS.InvalidValue("parameter $name is not stored"))
    end
    return container
end

function get_parameter_container(
    optimization_container::OptimizationContainer,
    name::Symbol,
    ::Type{T},
) where {T <: PSY.Component}
    return get_parameter_container(optimization_container, encode_symbol(T, name))
end

function get_parameter_container(
    optimization_container::OptimizationContainer,
    ref::UpdateRef,
)
    name = ref.access_ref
    dl = ref.data_label !== nothing ? ref.data_label : ""
    name = Symbol("$name" * "_" * dl)
    return get_parameter_container(optimization_container, name)
end

function get_parameter_array(optimization_container::OptimizationContainer, ref)
    return get_parameter_array(get_parameter_container(optimization_container, ref))
end

function assign_parameter!(
    optimization_container::OptimizationContainer,
    container::ParameterContainer,
)
    @debug "assign_parameter" container.update_ref
    name = container.update_ref.access_ref
    dl = container.update_ref.data_label !== nothing ? container.update_ref.data_label : ""
    name = Symbol(strip("$name" * "_" * dl))

    if haskey(optimization_container.parameters, name)
        @error "parameter $name is already stored" sort!(
            get_parameter_names(optimization_container),
        )
        throw(IS.InvalidValue("parameter $name is already stored"))
    end

    optimization_container.parameters[name] = container
    return
end

function add_param_container!(
    optimization_container::OptimizationContainer,
    param_reference::UpdateRef,
    axs...,
)
    container = ParameterContainer(
        param_reference,
        JuMP.Containers.DenseAxisArray{PJ.ParameterRef}(undef, axs...),
        fill!(JuMP.Containers.DenseAxisArray{Float64}(undef, axs...), NaN),
    )
    assign_parameter!(optimization_container, container)
    return container
end

function iterate_parameter_containers(optimization_container::OptimizationContainer)
    Channel() do channel
        for container in values(optimization_container.parameters)
            put!(channel, container)
        end
    end
end

function assign_expression!(
    optimization_container::OptimizationContainer,
    name::Symbol,
    value,
)
    @debug "set_expression" name
    optimization_container.expressions[name] = value
    return
end

function add_expression_container!(
    optimization_container::OptimizationContainer,
    exp_name::Symbol,
    axs...,
)
    container = JuMP.Containers.DenseAxisArray{JuMP.GenericAffExpr}(undef, axs...)
    assign_expression!(optimization_container, exp_name, container)
    return container
end

function is_milp(container::OptimizationContainer)
    type_of_optimizer = typeof(container.JuMPmodel.moi_backend.optimizer.model)
    supports_milp = hasfield(type_of_optimizer, :last_solved_by_mip)
    !supports_milp && return false
    return container.JuMPmodel.moi_backend.optimizer.model.last_solved_by_mip
end

function export_optimizer_stats(
    optimizer_stats::Dict{Symbol, Any},
    optimization_container::OptimizationContainer,
    path::String,
)
    optimizer_stats[:termination_status] =
        Int(JuMP.termination_status(optimization_container.JuMPmodel))
    optimizer_stats[:primal_status] =
        Int(JuMP.primal_status(optimization_container.JuMPmodel))
    optimizer_stats[:dual_status] = Int(JuMP.dual_status(optimization_container.JuMPmodel))

    if optimizer_stats[:primal_status] == MOI.FEASIBLE_POINT::MOI.ResultStatusCode
        optimizer_stats[:obj_value] = JuMP.objective_value(optimization_container.JuMPmodel)
    else
        optimizer_stats[:obj_value] = Inf
    end

    try
        optimizer_stats[:solve_time] =
            MOI.get(optimization_container.JuMPmodel, MOI.SolveTime())
    catch
        @warn("SolveTime() property not supported by the Solver")
        optimizer_stats[:solve_time] = NaN # "Not Supported by solver"
    end
    write_optimizer_stats(optimizer_stats, path)
    return
end

""" Exports the OpModel JuMP object in MathOptFormat"""
function serialize_optimization_model(
    optimization_container::OptimizationContainer,
    save_path::String,
)
    MOF_model = MOPFM(format = MOI.FileFormats.FORMAT_MOF)
    MOI.copy_to(MOF_model, JuMP.backend(optimization_container.JuMPmodel))
    MOI.write_to_file(MOF_model, save_path)
    return
end

function read_variables(optimization_container::OptimizationContainer)
    return Dict(
        k => axis_array_to_dataframe(v) for (k, v) in get_variables(optimization_container)
    )
end

function read_duals(optimization_container::OptimizationContainer)
    cons = get_constraint_duals(optimization_container.settings)
    return read_duals(optimization_container, cons)
end

function read_duals(op::OptimizationContainer, cons::Vector{Symbol})
    results_dict = Dict{Symbol, DataFrames.DataFrame}()
    isempty(cons) && return results_dict
    for c in cons
        v = get_constraint(op, c)
        results_dict[c] = axis_array_to_dataframe(v, [c])
    end
    return results_dict
end

function read_parameters(optimization_container::OptimizationContainer)
    # TODO: Still not obvious implementation since it needs to get the multipliers from
    # the system
    params_dict = Dict{Symbol, DataFrames.DataFrame}()
    parameters = get_parameters(optimization_container)
    (isnothing(parameters) || isempty(parameters)) && return params_dict
    for (k, v) in parameters
        !isa(v.update_ref, UpdateRef{<:PSY.Component}) && continue
        params_key_tuple = decode_symbol(k)
        params_dict_key = Symbol(params_key_tuple[1], "_", params_key_tuple[3])
        param_array = axis_array_to_dataframe(get_parameter_array(v))
        multiplier_array = axis_array_to_dataframe(get_multiplier_array(v))
        params_dict[params_dict_key] = param_array .* multiplier_array
    end
    return params_dict
end

function add_to_setting_ext!(
    optimization_container::OptimizationContainer,
    key::String,
    value,
)
    settings = get_settings(optimization_container)
    push!(get_ext(settings), key => value)
    @debug "Add to settings ext" key value
    return
end

function check_optimization_container(optimization_container::OptimizationContainer)
    valid = true
    # Check for parameter invalid values
    if model_has_parameters(optimization_container)
        for param_array in values(optimization_container.parameters)
            valid = !all(isnan.(param_array.multiplier_array.data))
        end
    end
    if !valid
        error("The model container has invalid values")
    end
    return
end

function get_problem_size(optimization_container::OptimizationContainer)
    model = optimization_container.JuMPmodel
    vars = JuMP.num_variables(model)
    cons = 0
    for (exp, c_type) in JuMP.list_of_constraint_types(model)
        cons += JuMP.num_constraints(model, exp, c_type)
    end
    return "The current total number of variables is $(vars) and total number of constraints is $(cons)"
end

function build_impl!(
    optimization_container::OptimizationContainer,
    template::OperationsProblemTemplate,
    sys::PSY.System,
)
    transmission = template.transmission
    # Order is required
    TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "Services" begin
        construct_services!(
            optimization_container,
            sys,
            template.services,
            template.devices,
        )
    end
    for device_model in values(template.devices)
        @debug "Building $(device_model.component_type) with $(device_model.formulation) formulation"
        TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "$(device_model.component_type)" begin
            construct_device!(optimization_container, sys, device_model, transmission)
            @debug get_problem_size(optimization_container)
        end
    end

    TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "$(transmission)" begin
        @debug "Building $(transmission) network formulation"
        construct_network!(optimization_container, sys, transmission, template)
        @debug get_problem_size(optimization_container)
    end

    for branch_model in values(template.branches)
        @debug "Building $(branch_model.component_type) with $(branch_model.formulation) formulation"
        TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "$(branch_model.component_type)" begin
            construct_device!(optimization_container, sys, branch_model, transmission)
            @debug get_problem_size(optimization_container)
        end
    end

    TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "Objective" begin
        @debug "Building Objective"
        JuMP.@objective(
            optimization_container.JuMPmodel,
            MOI.MIN_SENSE,
            optimization_container.cost_function
        )
    end
    @debug "Total operation count $(optimization_container.JuMPmodel.operator_counter)"

    check_optimization_container(optimization_container)
    return
end


text:
function acoeffs(ly::Float64,ny::Int,Ξ::Float64,τ::Float64=0.0;jw::Float64=0.05)
    ζjet = zeros(Float64,2*ny-1)
    Δθ::Float64 = jw
    κ::Float64 = τ == 0.0 ? 0.0 : 1.0/τ
    ζjet = [-κ*Ξ*tanh(-y/Δθ) for y in LinRange(-ly/2.0,ly/2.0,2*ny-1)]
    fftshift(fft(ζjet))
end

function bcoeffs(lx::Float64,ly::Float64,nx::Int,ny::Int,β::Float64,τ::Float64=0.0,νn::Float64=0.0)
    B = zeros(ComplexF64,2*ny-1,nx)
    α::Int = 2
    kxmax::Float64 = 2.0*Float64(pi)/lx*Float64(nx-1)
    kymax::Float64 = 2.0*Float64(pi)/ly*Float64(ny-1)
    γ::Float64 = τ == 0.0 ? 0.0 : 1.0/τ
    for m = 0:1:nx-1
        nmin = m == 0 ? 1 : -(ny-1)
        for n=nmin:1:ny-1
            kx::Float64 = 2.0*Float64(pi)*Float64(m)/lx
            ky::Float64 = 2.0*Float64(pi)*Float64(n)/ly
            B[n+ny,m+1] = -γ + im*β*kx/(kx^2 + ky^2) - νn*((kx^2 + ky^2)/(kxmax^2 + kymax^2))^(2*α)
        end
    end
    B
end

function ccoeffs(lx::Float64,ly::Float64,nx::Int,ny::Int)

    M::Int = nx - 1
    N::Int = ny - 1

    Cp = zeros(Float64,2*ny-1,nx,2*ny-1,nx)
    Cm = zeros(Float64,2*ny-1,nx,2*ny-1,nx)

    # ++ interactions note: +0 has only (0,+n)
    for m1=1:1:M
        for n1=-N:1:N
            for m2=0:1:min(m1,M-m1)

                n2min = m2 == 0 ? 1 : -N
                for n2=max(n2min,-N-n1):1:min(N,N-n1)

                    px::Float64 = 2.0*Float64(pi)/lx*Float64(m1)
                    py::Float64 = 2.0*Float64(pi)/ly*Float64(n1)
                    qx::Float64 = 2.0*Float64(pi)/lx*Float64(m2)
                    qy::Float64 = 2.0*Float64(pi)/ly*Float64(n2)

                    if m1 == m2
                        Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)/(px^2 + py^2)
                    else
                        Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))
                    end

                end
            end
        end
    end

    # +- interactions note: - includes (0,-n) because it is conj(0,n)
    for m1=1:1:M
        for n1=-N:1:N
            for m2=0:1:m1

                n2min = m2 == 0 ? 1 : -N
                n2max = m2 == m1 ? n1 - 1 : N
                for n2=max(n2min,n1-N):1:min(n2max,n1+N)

                    px::Float64 = 2.0*Float64(pi)/lx*Float64(m1)
                    py::Float64 = 2.0*Float64(pi)/ly*Float64(n1)
                    qx::Float64 = 2.0*Float64(pi)/lx*Float64(m2)
                    qy::Float64 = 2.0*Float64(pi)/ly*Float64(n2)

                    Cm[n2+ny,m2+1,n1+ny,m1+1] = (px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))

                end
            end
        end
    end
    Cp,Cm
end

function ccoeffs(lx::Float64,ly::Float64,nx::Int,ny::Int,Λ::Int)

    M::Int = nx - 1
    N::Int = ny - 1

    Cp = zeros(Float64,2*ny-1,nx,2*ny-1,nx)
    Cm = zeros(Float64,2*ny-1,nx,2*ny-1,nx)

    # L + L = L
    for m1=0:1:Λ
        n1min = m1 == 0 ? 1 : -N
        for n1=n1min:1:N
            for m2=0:1:min(m1,Λ-m1)

                n2min = m2 == 0 ? 1 : -N
                for n2=max(n2min,-N-n1):1:min(N,N-n1)

                    px::Float64 = 2.0*Float64(pi)/lx*Float64(m1)
                    py::Float64 = 2.0*Float64(pi)/ly*Float64(n1)
                    qx::Float64 = 2.0*Float64(pi)/lx*Float64(m2)
                    qy::Float64 = 2.0*Float64(pi)/ly*Float64(n2)

                    if m1 == m2
                        Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)/(px^2 + py^2)
                    else
                        Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))
                    end

                end
            end
        end
    end

    # L - L = L
    # note: -L should always include (0,-n)
    for m1=0:1:Λ
        n1min = m1 == 0 ? 1 : -N
        for n1=n1min:1:N

            for m2=0:1:m1

                n2min = m2 == 0 ? 1 : -N
                n2max = m2 == m1 ? n1 - 1 : N
                for n2=max(n2min,n1-N):1:min(n2max,n1+N)

                    px::Float64 = 2.0*Float64(pi)/lx*Float64(m1)
                    py::Float64 = 2.0*Float64(pi)/ly*Float64(n1)
                    qx::Float64 = 2.0*Float64(pi)/lx*Float64(m2)
                    qy::Float64 = 2.0*Float64(pi)/ly*Float64(n2)

                    Cm[n2+ny,m2+1,n1+ny,m1+1] = (px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))

                end
            end
        end
    end

    # H - H = L
    for m1=Λ+1:1:M
        for n1=-N:1:N
            for m2=max(Λ+1,m1-Λ):1:m1

                n2max = m2 == m1 ? n1 - 1 : N
                for n2=max(-N,n1-N):1:min(n2max,n1+N)

                    px::Float64 = 2.0*Float64(pi)/lx*Float64(m1)
                    py::Float64 = 2.0*Float64(pi)/ly*Float64(n1)
                    qx::Float64 = 2.0*Float64(pi)/lx*Float64(m2)
                    qy::Float64 = 2.0*Float64(pi)/ly*Float64(n2)

                    Cm[n2+ny,m2+1,n1+ny,m1+1] = (px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))

                end
            end
        end
    end

    # H + L = H
    for m1=Λ+1:1:M
        for n1=-N:1:N
            for m2=0:1:min(M-m1,Λ)

                n2min = m2 == 0 ? 1 : -N
                for n2=max(n2min,-N-n1):1:min(N,N-n1)

                    px::Float64 = 2.0*Float64(pi)/lx*Float64(m1)
                    py::Float64 = 2.0*Float64(pi)/ly*Float64(n1)
                    qx::Float64 = 2.0*Float64(pi)/lx*Float64(m2)
                    qy::Float64 = 2.0*Float64(pi)/ly*Float64(n2)

                    Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))

                end
            end
        end
    end

    # H - L = H
    # note: -L should always include (0,-n)
    for m1=Λ+1:1:M
        for n1=-N:1:N
            for m2=0:1:min(Λ,m1 - Λ - 1)

                n2min = m2 == 0 ? 1 : -N
                for n2=max(n2min,n1-N):1:min(N,n1+N)

                    px::Float64 = 2.0*Float64(pi)/lx*Float64(m1)
                    py::Float64 = 2.0*Float64(pi)/ly*Float64(n1)
                    qx::Float64 = 2.0*Float64(pi)/lx*Float64(m2)
                    qy::Float64 = 2.0*Float64(pi)/ly*Float64(n2)

                    Cm[n2+ny,m2+1,n1+ny,m1+1] = (px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))

                end
            end
        end
    end
    Cp,Cm
end


text:
@testset "lib.jl" begin
    @testset "accum" begin
        t1 = (a=1, b=2, c=3)
        t2 = (a=1, b=2)
        @test Zygote.accum(t1, t2) == (a = 2, b = 4, c = 3)
        @test_throws ArgumentError Zygote.accum(t2, t1)
    end
end


text:
using ImageContainers
using Test
using Base64, FileIO

function writebuf(c::ImageContainer, mime::AbstractString)
    buf = IOBuffer()
    show(buf, MIME(mime), c)
    take!(buf)
end

@testset "ImageContainers.jl" begin
    data = rand(UInt8, 4)
    b64data = base64encode(data)

    @test ImageContainer{:a, Int}(3).content == 3
    for (fmt, mime) in ImageContainers.plainmimes
        c = storeimage(fmt, data)
        @test c.content == data
        @test writebuf(c, mime) == data
    end
    for fmt in (:gif, :bmp)
        c = storeimage(fmt, data)
        @test c.content == data
        @test String(writebuf(c, "text/html")) ==
            ""
    end
    for fmt in (:mp4, :webm)
        for mime in ("text/html", "application/prs.juno.plotpane+html")
            c = storeimage(fmt, data)
            @test c.content == data
            @test String(writebuf(c, mime)) ==
                "<video controls autoplay loop " *
                "src=\"data:video/$fmt;base64,$b64data\" />"
        end
    end
    begin # Julia image
        c = storeimage(:jlc, data)
        buf = IOBuffer()
        save(format"BMP", buf, data)
        out1 = writebuf(c, "image/bmp") # Juno
        @test out1 == take!(buf)
        out2 = writebuf(c, "text/html") # Jupyter
        @test String(out2) ==
            ""
    end
end


text:
module HistoryTest

using MAFE
using Test
using DataFrames

@testset "Migration Histories" begin
    ctry = Ghana;
    df = read_migration_file(ctry);

    id = "G0003001";
    mh = make_migration_history(id, df);
    @test n_spells(mh) == 0

    id = "G0016001";
    mh = make_migration_history(id, df);
    @test n_spells(mh) == 2
    mig1 = get_spell(mh, 1);
    @test mig1.startYear == 1981
    @test mig1.endYear == 1984
    @test mig1.ctry == 99338

    @test find_country(mh, mig1.startYear, mig1.endYear) == mig1.ctry
    @test find_country(mh, 1940, 1950) == MAFE.homeCountry
    @test find_country(mh, 2000, 2005) == MAFE.homeCountry
    @test find_country(mh, 1980, 1984) == MAFE.unknownCountry
end

@testset "Job histories" begin
    ctry = Ghana;
    df = read_activity_file(ctry);

    id = "G0003001";
    jh = make_job_history(id, df);
    @test n_spells(jh) == 1
    job1 = get_spell(jh, 1);
    @test job1.startYear == 2008

    # Person with multiple spells
    id = "G0008001";
    jh = make_job_history(id, df);
    @test n_spells(jh) == 2
    job1 = get_spell(jh, 2);
    @test job1.startYear == 2003
end

end


text:
# helper

function createconnection(name, matrixtype)
	@eval begin
		type $name <: Connection
			w::$matrixtype
			pre::Neuron
			post::Neuron
			prename::Symbol
			postname::Symbol
		end
		export $name
	end
end

function updateconnection!(c::Connection)
end
export updateconnection!

# sparse connections

type SparseConnectionArray
	I::Array{Int64,1}	# post
	J::Array{Int64,1}	# pre
	V::Array{FloatXX,1} # value
	n::Int64
	m::Int64
	nnz::Int64
end
export SparseConnectionArray

import Base.A_mul_B!
function A_mul_B!(α::FloatXX,
				  c::SparseConnectionArray,
				  v::Array{FloatXX,1},
				  β::FloatXX,
				  res::Array{FloatXX, 1})
	@inbounds for i = 1:c.nnz
		res[c.I[i]] = α * c.V[i] * v[c.J[i]] + β * res[c.I[i]]
	end
	res
end
import Base.*
function *(c::SparseConnectionArray, v::Array{Float64,1})
	res = zeros(c.n)
	@inbounds for i = 1:c.nnz
		res[c.I[i]] += c.V[i] * v[c.J[i]]
	end
	res
end

function sampleconnections(n_ofpre, n_ofpost, nnz)
	c = Tuple[]
	while length(c) < nnz
		sample = (rand(1:n_ofpost), rand(1:n_ofpre))
		if !(sample in c); push!(c, sample); end
	end
	[c[i][1] for i in 1:nnz], [c[i][2] for i in 1:nnz]
end

using Distributions
function createsparseconnection(name)
	createconnection(name, SparseConnectionArray)
	@eval begin
		function $name(pre::Layer, post::Layer; p = .02)
			nnz = rand(Binomial(pre.n_of * post.n_of, p))
			I, J = sampleconnections(pre.n_of, post.n_of, nnz)
			$name(SparseConnectionArray(I, J, randn(FloatXX, nnz)/(10*p*sqrt(pre.n_of)),
										 post.n_of, pre.n_of, nnz),
				  pre.neurons,
				  post.neurons,
				  pre.name,
				  post.name)
		end
	end
end

createsparseconnection(:PlasticSparseConnection)
createsparseconnection(:StaticSparseConnection)

# dense connections
function createdenseconnection(name)
	createconnection(name, Array{FloatXX, 2})
	@eval begin
		function $name(pre::Layer, post::Layer)
			$name(randn(FloatXX, post.n_of, pre.n_of)/(10*sqrt(pre.n_of)),
				  pre.neurons,
				  post.neurons,
				  pre.name,
				  post.name)
		end
	end
end

createdenseconnection(:PlasticDenseConnection)
createdenseconnection(:StaticDenseConnection)

# One2OneConnection
type IdentityMatrix
end

createconnection(:One2OneConnection, IdentityMatrix)
function One2OneConnection(pre::Layer, post::Layer)
	One2OneConnection(IdentityMatrix(),
				      pre.neurons,
					  post.neurons,
					  pre.name,
					  post.name)
end

# TransposeDenseConnection
createconnection(:TransposeDenseConnection, Array{FloatXX, 2})
"""
	TransposeDenseConnection(c::PlasticDenseConnection, net::SimpleNetwork)

Create a connection from `c.postname` to `c.prename` and a pointer to `c.w`.
"""
function TransposeDenseConnection(c::PlasticDenseConnection, net::SimpleNetwork)
	TransposeDenseConnection(c.w,
						     net.layers[c.postname].neurons,
							 net.layers[c.prename].neurons,
							 c.postname,
							 c.prename)
end

# functions
"""
	weightedprerates!(v::Array{FloatXX, 1}, cons::Dict{Symbol, Connection}, scale::FloatXX)

Update `v = v + scale * ∑_{con ∈ cons} con.w * con.pre.outp`.

	weightedprerates!(v::Array{FloatXX, 1}, con::Connection, scale::FloatXX)

Update `v = v + scale * con.w * con.pre.outp`.
"""
function weightedprerates!(v::Array{FloatXX, 1},
						   con::Dict{Symbol, Connection},
						   scale::FloatXX)
	@inbounds for c in values(con)
		weightedprerates!(v, c, scale)
	end
end

function weightedprerates!(v, c::One2OneConnection, scale)
	BLAS.axpy!(scale, c.pre.outp, v);
end

function weightedprerates!(v, c::TransposeDenseConnection, scale)
	BLAS.gemv!('T', scale, c.w, c.pre.outp, FloatXX(1.), v);
end

function weightedprerates!(v, c::Union{PlasticSparseConnection,
									   StaticSparseConnection}, scale)
	A_mul_B!(scale, c.w, c.pre.outp, FloatXX(1.), v);
end

function weightedprerates!(v, c::Union{PlasticDenseConnection,
									   StaticDenseConnection}, scale)
	BLAS.gemv!('N', scale, c.w, c.pre.outp, FloatXX(1.), v);
end
export weightedprerates!

"""
	copyweights!(netfrom::SimpleNetwork, netto::SimpleNetwork)

Deepcopy weight parameters from `netfrom` to `netto`
"""
function copyweights!(netfrom::SimpleNetwork, netto::SimpleNetwork)
	for (ln, l) in netto.layers
        for (icn, ic) in l.inputconnections
            for (pren, pre) in ic
                pre.w = deepcopy(netfrom.layers[ln].inputconnections[icn][pren].w)
            end
        end
    end
end
export copyweights!

max_stars_repo_name: jbrea/NetSim

max_issues_repo_name: jbrea/NetSim

max_forks_repo_name: jbrea/NetSim


text:
using BinDeps
using CxxWrap

@BinDeps.setup

jlcxx_dir = joinpath(dirname(CxxWrap._l_jlcxx), "cmake", "JlCxx")

eigen_wrapper = library_dependency("eigen_wrapper", aliases=["libeigen_wrapper"])

cmake_prefix = jlcxx_dir

prefix=joinpath(BinDeps.depsdir(eigen_wrapper),"usr")
eigen_wrapper_srcdir = joinpath(BinDeps.depsdir(eigen_wrapper),"src","eigen_wrapper")
eigen_wrapper_builddir = joinpath(BinDeps.depsdir(eigen_wrapper),"builds","eigen_wrapper")
lib_prefix = @static is_windows() ? "" : "lib"
lib_suffix = @static is_windows() ? "dll" : (@static is_apple() ? "dylib" : "so")

makeopts = ["--", "-j", "$(Sys.CPU_CORES+2)"]

# Set generator if on windows
genopt = "Unix Makefiles"

build_type = get(ENV, "CXXWRAP_BUILD_TYPE", "Release")

eigen_steps = @build_steps begin
	`cmake -G "$genopt" -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_BUILD_TYPE="$build_type" -DCMAKE_PREFIX_PATH="$cmake_prefix" -DCMAKE_PROGRAM_PATH="$JULIA_HOME" $eigen_wrapper_srcdir`
	`cmake --build . --config $build_type --target install $makeopts`
end

# If built, always run cmake, in case the code changed
if isdir(eigen_wrapper_builddir)
  BinDeps.run(@build_steps begin
    ChangeDirectory(eigen_wrapper_builddir)
    eigen_steps
  end)
end

provides(BuildProcess,
  (@build_steps begin
    CreateDirectory(eigen_wrapper_builddir)
    @build_steps begin
      ChangeDirectory(eigen_wrapper_builddir)
      FileRule(joinpath(prefix,"lib", "$(lib_prefix)eigen_wrapper.$lib_suffix"),eigen_steps)
    end
  end),eigen_wrapper)

@BinDeps.install Dict([(:eigen_wrapper, :_l_eigen_wrap)])

@static if is_windows()
  if build_on_windows
    empty!(BinDeps.defaults)
    append!(BinDeps.defaults, saved_defaults)
  end
end