text:
# fft version
using FFTW

# *******************************************
#  Configuration
# *******************************************
"""
    ConfigFFT{T <: Union{Float64,Complex{Float64}},N}

Assembly all basic configurations of the Fourier-spectral method.
T represents the type of values of functions in 'real' space,
and N represents the dimension. The followings are all members:

    ngrids::NTuple{N,Int}
... Each element means the number of grids in each axis

    xranges::NTuple{N,NTuple{2,Float64}}
... Tuple of ((xmin, xmax), (ymin, ymax), ...)

    Xcoords::NTuple{N,Array{Float64,N}}
... Coordinates of 'real' space in each axis

    Kcoords::NTuple{N,Array{Complex{Float64},N}}
... Wavenumbers in each axis. It is used to calculate derivatives.

    P_fft, P_ifft, P_fftpad, P_ifftpad
... Operators used to accelerate fft and ifft

    cut_zigzag_mode::Bool
... The highest mode is suppressed or not for axes whose grids number is even.
If you use the easy constructor (explained below), it is set to be 'true' by default.
This flag is used in 'X_K' and 'K_X' function, the fourier transformation. If you
try to accelerate as much as possible and can manage the highest wavenumber mode,
set 'cut_zigzag_mode=false'.

For normal use, the easy constructor below is recommended for initialization.

    ConfigFFT(
        ngrids::NTuple{N,Int},
        xranges::NTuple{N,NTuple{2,Float64}};
        use_complex::Bool=false,
        cut_zigzag_mode::Bool=true)

The function in 'real' space is set to be real-valued (Float64) by default.
If you want to use complex-valued function, set 'use_complex=true'.

# Examples
1) 1-dimensional case. The number of grid is 64, and the range is 0 ≦ x < 10.

    julia> ngrids = (64,);

    julia> xranges = ((0., 10.),);

    julia> config = ConfigFFT(ngrids, xranges);

2) 2-dimensional case. The numbers of grids are 128 in x-axis and 256 in y-axis.
   The ranges in each axis are 0 ≦ x < 2π and 0 ≦ y < 4π.

    julia> ngrids = (128, 256);

    julia> xranges = ((0., 2π), (0., 4π));

    julia> config = ConfigFFT(ngrids, xranges);

3) 1-dimensional case. The number of grid is 64, and the range is 0 ≦ x < 30.
   The complex-valued functions in 'real' space is used.

   julia> ngrids = (64,);

   julia> xranges = ((0., 30.),);

   julia> config = ConfigFFT(ngrids, xranges, use_complex=true);
"""
struct ConfigFFT{T<:Union{Float64,Complex{Float64}},N}

    ngrids::NTuple{N,Int} # Number Of Grids

    xranges::NTuple{N,NTuple{2,Float64}} # tuple of (min, max) s
    Xcoords::NTuple{N,Array{Float64,N}}
    Kcoords::NTuple{N,Array{Complex{Float64},N}}

    P_fft::FFTW.cFFTWPlan{T1} where T1
    P_ifft::AbstractFFTs.ScaledPlan{T2} where T2
    P_fftpad::FFTW.cFFTWPlan{T3} where T3
    P_ifftpad::AbstractFFTs.ScaledPlan{T4} where T4

    cut_zigzag_mode::Bool

    # CONSTRUCTOR
    function ConfigFFT{T,N}(
            ngrids, xranges, Xcoords, Kcoords,
            P_fft, P_ifft, P_fftpad, P_ifftpad;
            cut_zigzag_mode=true
        ) where {T,N}

        warn_if_ngrid_is_not_power_of_2(ngrids)

        new(ngrids, xranges, Xcoords, Kcoords,
            P_fft, P_ifft, P_fftpad, P_ifftpad,
            cut_zigzag_mode)
    end

    # EASY CONSTRUCTOR
    function ConfigFFT(
            ngrids::NTuple{N,Int},
            xranges::NTuple{N,NTuple{2,Float64}};
            use_complex::Bool=false,
            cut_zigzag_mode::Bool=true
        ) where N

        T = use_complex ? Complex{Float64} : Float64

        carts = CartesianIndices(ngrids)
        Xcoords = Xcoordsgen(ngrids, xranges)
        Kcoords = Kcoordsgen(ngrids, xranges)

        P_fft = plan_fft(Xcoords[1])
        P_ifft = plan_ifft(Xcoords[1])

        pad_ngrids = @. ngrids ÷ 2 * 3
        P_fftpad = plan_fft(zeros(pad_ngrids))
        P_ifftpad = plan_ifft(zeros(pad_ngrids))

        return ConfigFFT{T,N}(
            ngrids, xranges, Xcoords, Kcoords,
            P_fft, P_ifft, P_fftpad, P_ifftpad,
            cut_zigzag_mode=cut_zigzag_mode)
    end

end

# --- helper function for constuctor ---
function Xcoordsgen(
        ngrids::NTuple{N,Int},
        xranges::NTuple{N,NTuple{2,Float64}}
    ) where N

    _Xcoordgen(axis) = Xcoordgen(axis, ngrids, xranges)
    return ntuple(_Xcoordgen, N)

end

function Xcoordgen(
        axis::Int,
        ngrids::NTuple{N,Int},
        xranges::NTuple{N,NTuple{2,Float64}}
    ) where N

    ngrid = ngrids[axis]
    xrange = xranges[axis]
    _Xcoordgen(indices) = (
        (  (indices[axis] - 1)*xrange[2]
         + (ngrid - indices[axis] + 1)*xrange[1] ) / ngrid
    )
    return _Xcoordgen.(CartesianIndices(ngrids))

end

function Kcoordsgen(
        ngrids::NTuple{N,Int},
        xranges::NTuple{N,NTuple{2,Float64}}
    ) where N

    _Kcoordgen(axis) = Kcoordgen(axis, ngrids, xranges)
    return ntuple(_Kcoordgen, N)

end

function Kcoordgen(
        axis::Int,
        ngrids::NTuple{N,Int},
        xranges::NTuple{N,NTuple{2,Float64}}
    ) where N

    ngrid = ngrids[axis]
    xrange = xranges[axis]
    _Kcoordgen(indices) = (
        kval(indices[axis], ngrid, xrange)
    )

    return _Kcoordgen.(CartesianIndices(ngrids))

end

function kval(index, ngrid, xrange)::Complex{Float64}

    index0 = index - 1  # start by 0
    if 2*index0 < ngrid
        return index0
    elseif 2*index0 == ngrid
        return 0
    elseif 2*index0 > ngrid
        return index0 - ngrid
    end

end

# *******************************************
#  XFunc
# *******************************************
"""
    XFunc{T,N} <: AbstractArray{T,N}

A struct of a function on 'real' space.
It consists of the following:

    vals::Array{T,N}
... Values of the function on grids

    config::ConfigFFT{T,N}
... Configuration of Fourier Transformation

For normal use, the easy constructor below is recommended for initialization.

    XFunc(
        vals::Array{Tc,N},
        config::ConfigFFT{T,N}
        ) where {N, Tc<:Number, T}

The type of elements of input 'vals' is converted to T.
Instead, you may generate a XFunc object, whose values are undef as

    XFunc(undef, config)

XFunc type object is transformed by 'K_X' to KFunc type object,
which represents a function in 'wavenumber' space.
See also the explanation of KFunc.

# Examples
Suppose 'config' be ConfigFFT{Float64,2} type object,
and c.ngrids = (2, 2). Then, for example,

    v = [1. 2.; 3. 4.]
    X_func = XFunc(v, config)

generates a XFunc object whose 'vals' is v.
If 'K_X' is applied to this function, say,

    K_func = K_X(X_func)

a new object 'K_func' is returned.
It is KFunc type, and its 'vals' is Fourier transform of v.
"""
mutable struct XFunc{T,N} <: AbstractArray{T,N}

    vals::Array{T,N}
    config::ConfigFFT{T,N}

    # CONSTRUCTOR
    function XFunc{T,N}(
            vals::Array{T,N},
            config::ConfigFFT{T,N}
        ) where {T,N}

        check_size_consistency(vals, config)
        return new(vals, config)

    end

    # EASY CONSTRUCTOR
    function XFunc(
            vals::Array{Tv,N},
            config::ConfigFFT{Tc,N}
        ) where {N, Tv<:Number, Tc}

        if Tv <: Complex && Tc <: Real
            return XFunc{Tc,N}(Tc.(real(vals)), config)
        else
            return XFunc{Tc,N}(Tc.(vals), config)
        end

    end

    # UNDEF CONSTRUCTOR
    function XFunc(
            undef::UndefInitializer,
            config::ConfigFFT{T,N}
        ) where {T,N}

        f_undef = Array{T,N}(undef, config.ngrids)
        return XFunc{T,N}(f_undef, config)

    end

end

Base.:size(f::XFunc) = size(f.vals)

Base.:getindex(f::XFunc, i::Int) = getindex(f.vals, i)

function Base.:getindex(
        f::XFunc{T,N}, I::Vararg{Int,N}
    ) where {T,N}

    getindex(f.vals, I...)
end

Base.:setindex!(f::XFunc, v, i::Int) = setindex!(f.vals, v, i)

function Base.:setindex!(
        f::XFunc{T,N}, v, I::Vararg{Int,N}
    ) where {T,N}

    setindex!(f.vals, v, I...)
end

Base.:copy(f::XFunc) = XFunc(copy(f.vals), f.config)

# *******************************************
#  KFunc
# *******************************************
"""
    KFunc{T,N} <: AbstractArray{T,N}

A struct of a function on 'wavenumber' space.
It consists of the following:

    vals::Array{Complex{Float64},N}
... Values of the function on grids

    config::ConfigFFT{T,N}
... Configuration of Fourier Transformation

For normal use, the easy constructor below is recommended for initialization.

    KFunc(
        vals::Array{Tc,N},
        config::ConfigFFT{T,N}
    ) where {N, Tc<:Number, T}

The type of elements of input 'vals' is converted to T.
Instead, you may generate a KFunc object, whose values are undef as

    KFunc(undef, config)

KFunc type object is transformed by 'X_K' to XFunc type object,
which represents a function in 'real' space.
See also the explanation of XFunc.

# Examples
Suppose 'config' be ConfigFFT{Float64,2} type object,
and c.ngrids = (2, 2). Then, for example,

    v = [1. 2.; 3. 4.]
    K_func = KFunc(v, config)

generates a XFunc object whose 'vals' is v.
If 'X_K' is applied to this function, say,

    X_func = X_K(K_func)

a new object 'X_func' is returned.
It is XFunc type, and its 'vals' is inverse-Fourier transform of v.
"""
mutable struct KFunc{T,N} <: AbstractArray{Complex{Float64},N}

    vals::Array{Complex{Float64},N}
    config::ConfigFFT{T,N}

    # CONSTRUCTOR
    function KFunc{T,N}(
            vals::Array{Complex{Float64},N},
            config::ConfigFFT{T,N}
        ) where {T,N}

        check_size_consistency(vals, config)
        return new(vals, config)

    end

    # EASY CONSTRUCTOR
    function KFunc(
            vals::Array{Tv,N},
            config::ConfigFFT{Tc,N}
        ) where {N, Tv<:Number, Tc}

        KFunc{Tc,N}(Complex{Float64}.(vals), config)
    end

    # UNDEF CONSTRUCTOR
    function KFunc(
            undef::UndefInitializer,
            config::ConfigFFT{T,N}
        ) where {T,N}

        f_undef = Array{Complex{Float64},N}(undef, config.ngrids)
        return KFunc{T,N}(f_undef, config)

    end

end

Base.:size(f::KFunc) = size(f.vals)

Base.:getindex(f::KFunc, i::Int) = getindex(f.vals, i)

function Base.:getindex(
        f::KFunc{T,N}, I::Vararg{Int,N}
    ) where {T,N}

    getindex(f.vals, I...)
end

Base.:setindex!(f::KFunc, v, i::Int) = setindex!(f.vals, v, i)

function Base.:setindex!(
        f::KFunc{T,N}, v, I::Vararg{Int,N}
    ) where {T,N}

    setindex!(f.vals, v, I...)
end

Base.:copy(f::KFunc) = KFunc(copy(f.vals), f.config)

# *******************************************
#  Binomial Operators
# *******************************************
const BINOP = (
    (:+, :.+), (:-, :.-), (:*, :.*),
    (:/, :./), (:\, :.\), (:^, :.^)
)

# +++++ XFunc +++++
Base.:+(f::XFunc) = f
Base.:-(f::XFunc) = XFunc(-f.vals, f.config)
Base.:^(f::XFunc, n::Integer) = XFunc(f.vals .^ n, f.config)
Base.:inv(f::XFunc) = f \ 1.0

for (op, opd) = BINOP
    @eval begin
        function Base.$op(f::XFunc, g::XFunc)
            check_config_consistency(f.config, g.config)
            return XFunc($opd(f.vals, g.vals), f.config)
        end

        Base.$op(f::XFunc, a::Number) = (
            XFunc($opd(f.vals, a), f.config)
        )
        Base.$op(a::Number, f::XFunc) = (
               XFunc($opd(a, f.vals), f.config)
        )
    end
end

# +++++ KFunc +++++
Base.:+(f::KFunc) = f
Base.:-(f::KFunc) = KFunc(-f.vals, f.config)
Base.:^(f::KFunc, n::Integer) = KFunc(f.vals .^ n, f.config)
Base.:inv(f::KFunc) = f \ 1.0

for (op, opd) = BINOP
    @eval begin
        function Base.$op(f::KFunc, g::KFunc)
            check_config_consistency(f.config, g.config)
            return KFunc($opd(f.vals, g.vals), f.config)
        end

        Base.$op(f::KFunc, a::Number) = (
            KFunc($opd(f.vals, a), f.config)
        )
        Base.$op(a::Number, f::KFunc) = (
            KFunc($opd(a, f.vals), f.config)
        )
    end
end

# *******************************************
#  Operators for Complex Numbers
# *******************************************
const OPERATOR = (
    :real, :imag, :conj
)

for op = OPERATOR
    @eval begin
        Base.$op(f::XFunc) = XFunc($op(f.vals), f.config)
        Base.$op(f::KFunc) = KFunc($op(f.vals), f.config)
    end
end

Base.:reim(f::XFunc) = (real(f), imag(f))
Base.:reim(f::KFunc) = (real(f), imag(f))

# *******************************************
#  Low/High-pass Filter
# *******************************************
# +++++ general pass filter +++++
function pass_K!(
        f::KFunc{T,N},
        slices::NTuple{N,UnitRange{Int}}
    ) where {T,N}

    vals = copy(f.vals)
    vals[slices...] .= 0.0 + 0.0im
    f.vals -= vals

end

# +++++ high-pass filter +++++
function highpass_K!(
        f::KFunc{T,N},
        min_nwaves::NTuple{N,Int}
    ) where {T,N}

    ngrids = f.config.ngrids
    if any(@. min_nwaves > ngrids ÷ 2)
        println("WARNING: all waves are suppressed")
        f.vals .= 0.
    else
        floors = tuple(ones(Int, N)...)
        ceils = ngrids
        min_indices = @. max(floors, floors + min_nwaves)
        max_indices = @. min(ceils, ceils + 1 - min_nwaves)

        slices = (:).(min_indices, max_indices)
        pass_K!(f, slices)
    end
end

function K_highpass_K(
        f::KFunc{T,N},
        min_nwaves::NTuple{N,Int}
    ) where {T,N}

    g = copy(f)
    highpass_K!(g, min_nwaves)
    return g

end

# +++++ low-pass filter +++++
function lowpass_K!(
        f::KFunc{T,N},
        max_nwaves::NTuple{N,Int}
    ) where {T,N}

    ngrids = f.config.ngrids
    f.vals .= fftshift(f.vals)

    center_indices = @. ngrids ÷ 2 + 1
    if any(max_nwaves .< 0)
        println("WARNING: all waves are suppressed")
        f.vals .= 0.
    else
        floors = tuple(ones(Int, N)...)
        ceils = ngrids
        min_indices = (
            @. max(floors, center_indices - max_nwaves)
        )
        max_indices = (
            @. min(ceils, center_indices + max_nwaves)
        )

        slices = (:).(min_indices, max_indices)
        pass_K!(f, slices)
    end

    f.vals .= ifftshift(f.vals)

end

function K_lowpass_K(
        f::KFunc{T,N},
        max_nwaves::NTuple{N,Int}
    ) where {T,N}

    g = copy(f)
    lowpass_K!(g, max_nwaves)
    return g

end

# *******************************************
#  Fourier Transformation
# *******************************************
function K_X(f::XFunc{T,N}) where {T,N}

    P = f.config.P_fft
    g = KFunc(P * f.vals, f.config)

    if f.config.cut_zigzag_mode
        ngrids = f.config.ngrids
        max_nwaves = @. ngrids ÷ 2 - 1
        lowpass_K!(g, max_nwaves)
    end

    return g
end

function X_K(f::KFunc{T,N}) where {T,N}

    P = f.config.P_ifft
    if T <: Real
        XFunc(real(P * f.vals), f.config)
    else
        XFunc(P * f.vals, f.config)
    end
end

# *******************************************
#  Differentiation
# *******************************************
# +++++ destructive +++++
function ∂Xaxis_K!(f::KFunc, axis::Int)

    config = f.config
    xrange = config.xranges[axis]
    xlen = xrange[2] - xrange[1]
    Kcoord = config.Kcoords[axis]

    f.vals .*= (2π*im/xlen)*Kcoord

end

# +++++ non-destructive +++++
function K_∂Xaxis_K(f::KFunc, axis::Int)
    g = copy(f)
    ∂Xaxis_K!(g, axis)
    return g
end

function X_∂Xaxis_X(f::XFunc, axis::Int)
    X_K(K_∂Xaxis_K(K_X(f), axis))
end

function K_laplacian_K(f::KFunc{T,N}) where {T,N}

    return sum(
        K_∂Xaxis_K(K_∂Xaxis_K(f, axis), axis)
        for axis = 1:N
    )

end

K_Δ_K = K_laplacian_K
X_laplacian_X = X_K ∘ K_laplacian_K ∘ K_X
X_Δ_X = X_laplacian_X

"""
K_laplainv_K returns g satisfying Δg = f
"""
function K_laplainv_K(f::KFunc{T,N}) where {T,N}

    c = f.config
    xlens = xlens_xranges(c.xranges)

    lap = - 4π^2 * sum(
        (x -> x .^ 2).(c.Kcoords) ./ (xlens .^ 2)
    )
    gvals = f.vals ./ lap

    zero_indiceses = (
        map(c.ngrids) do x
            if x % 2 == 0
                return (1, x÷2 + 1)
            else
                return (1,)
            end
        end
    )
    for indices in Iterators.product(zero_indiceses...)
        gvals[indices...] = 0.
    end

    return KFunc(gvals, c)

end

K_Δ⁻¹_K = K_laplainv_K
X_laplainv_X = X_K ∘ K_laplainv_K ∘ K_X
X_Δ⁻¹_X = X_laplainv_X

# *******************************************
#  De-aliased Products
# *******************************************
# 3/2-rule ... zero-padding
# 2/3-rule ... truncating

# de-aliased product by 3/2-rule (zero padding)
function K_prod_32_K_K(
        f::KFunc{T,N}, g::KFunc{T,N}
    ) where {T,N}

    check_config_consistency(f.config, g.config)

    fvals_pad = padding(f)
    gvals_pad = padding(g)

    P_fft = f.config.P_fftpad
    P_ifft = f.config.P_ifftpad
    if T <: Real
        fgvals_pad = P_fft * (
            real(P_ifft * fvals_pad) .* real(P_ifft * gvals_pad)
        )
    else
        fgvals_pad = P_fft * (
            (P_ifft * fvals_pad) .* (P_ifft * gvals_pad)
        )
    end
    fgvals = truncate(fgvals_pad, f.config)

    return KFunc(fgvals, f.config)

end

# dot-product
function K_dot_32_Ks_Ks(
        f::Array{KFunc{T,N},1},
        g::Array{KFunc{T,N},1}
    ) where {T,N}

    check_length_consistency(f, g)
    ndim = length(f)
    for dim = 1:ndim
        check_config_consistency(
            f[dim].config, g[dim].config
        )
        if dim < ndim
            check_config_consistency(
                f[dim].config, f[dim+1].config
            )
        end
    end

    config = f[1].config

    fvals_pad = padding.(f)
    gvals_pad = padding.(g)

    P_fft = config.P_fftpad
    P_ifft = config.P_ifftpad

    fgvals_pad = P_fft * sum(
        (P_ifft * fvals_pad[dim]) .* (P_ifft * gvals_pad[dim])
        for dim = 1:ndim
    )

    fgvals = truncate(fgvals_pad, config)
    return KFunc(fgvals, config)

end

# cross-product
function K3_cross_32_K3_K3(
        f::Array{KFunc{T,N},1},
        g::Array{KFunc{T,N},1}
    ) where {T,N}

    if !(length(f) == length(g) == 3)
        errmsg = (
            "lengths of arrays must be 3"
        )
        throw(ErrorException(errmsg))
    end

    config = f[1].config

    fvals_pad = padding.(f)
    gvals_pad = padding.(g)

    P_fft = config.P_fftpad
    P_ifft = config.P_ifftpad

    xfvals_pad, yfvals_pad, zfvals_pad = [
        P_ifft * fvals_pad[dim] for dim = 1:3
    ]
    xgvals_pad, ygvals_pad, zgvals_pad = [
        P_ifft * gvals_pad[dim] for dim = 1:3
    ]

    fgvals = [
        truncate( P_fft * (
            yfvals_pad .* zgvals_pad .- zfvals_pad .* ygvals_pad
        ), config),
        truncate( P_fft * (
            zfvals_pad .* xgvals_pad .- xfvals_pad .* zgvals_pad
        ), config),
        truncate( P_fft * (
            xfvals_pad .* ygvals_pad .- yfvals_pad .* xgvals_pad
        ), config),
    ]

    return [
        KFunc(fgvals[dim], config) for dim = 1:3
    ]

end

# padding helper functions
function padding(f::KFunc)
    config = f.config

    ngrids = config.ngrids
    pad_ngrids = to_pad_ngrids(ngrids)

    slices = slices_padded_core(ngrids, pad_ngrids)

    vals_shift = fftshift(f.vals)
    padded = zeros(Complex{Float64}, pad_ngrids)
    padded[slices...] = vals_shift
    padded .= ifftshift(padded)

    return padded
end

function truncate(
        padded::Array{Complex{Float64},N},
        config::ConfigFFT{T,N}
    ) where {T,N}

    ngrids = config.ngrids
    pad_ngrids = to_pad_ngrids(ngrids)

    slices = slices_padded_core(ngrids, pad_ngrids)
    vals = ifftshift(fftshift(padded)[slices...])

    N_origin = prod(ngrids)
    N_padded = prod(pad_ngrids)
    vals *= N_padded / N_origin

    return vals

end

# helper function
function slices_padded_core(
        ngrids::NTuple{N1,Int},
        pad_ngrids::NTuple{N2,Int}
    ) where {N1,N2}

    nshifts = @. pad_ngrids÷2 - ngrids÷2
    min_nwaves = @. nshifts + 1
    max_nwaves = @. nshifts + ngrids
    return (:).(min_nwaves, max_nwaves)

end

function to_pad_ngrids(ngrids::NTuple{N,Int}) where N
    @. ngrids ÷ 2 * 3
end

# de-aliased product by 2/3-rule (truncation)
function K_prod_23_K_K(f::KFunc, g::KFunc)

    check_config_consistency(f.config, g.config)
    config = f.config
    ngrids = config.ngrids
    max_nwaves = @. ngrids ÷ 3
    f_trunc = K_lowpass_K(f, max_nwaves)
    g_trunc = K_lowpass_K(g, max_nwaves)
    return K_X(X_K(f_trunc) * X_K(g_trunc))

end

# dot-product
function K_dot_23_Ks_Ks(
        f::Array{KFunc{T,N},1},
        g::Array{KFunc{T,N},1}
    ) where {T,N}

    check_length_consistency(f, g)
    ndim = length(f)
    for dim = 1:ndim
        check_config_consistency(
            f[dim].config, g[dim].config
        )
        if dim < ndim
            check_config_consistency(
                f[dim].config, f[dim+1].config
            )
        end
    end

    config = f[1].config
    ngrids = config.ngrids
    max_nwaves = ngrids .÷ 3
    f_trunc = [
        K_lowpass_K(f[dim], max_nwaves)
        for dim = 1:ndim
    ]
    g_trunc = [
        K_lowpass_K(g[dim], max_nwaves)
        for dim = 1:ndim
    ]
    return K_X(
        sum( X_K.(f_trunc) .* X_K.(g_trunc) )
    )

end

# cross-product
function K3_cross_23_K3_K3(
        f::Array{KFunc{T,N},1},
        g::Array{KFunc{T,N},1}
    ) where {T,N}

    if !(length(f) == length(g) == 3)
        errmsg = (
            "lengths of arrays must be 3"
        )
        throw(ErrorException(errmsg))
    end

    config = f[1].config
    ngrids = config.ngrids
    max_nwaves = ngrids .÷ 3
    X_xf, X_yf, X_zf = X_K.([
        K_lowpass_K(f[dim], max_nwaves) for dim = 1:3
    ])
    X_xg, X_yg, X_zg = X_K.([
        K_lowpass_K(g[dim], max_nwaves) for dim = 1:3
    ])

    return K_X.([
        X_yf * X_zg - X_zf * X_yg,
        X_zf * X_xg - X_xf * X_zg,
        X_xf * X_yg - X_yf * X_xg
    ])

end

# +++++ aliases +++++
# --- 3/2-rule ---
# \ast
∗(f, g) = K_prod_32_K_K(f, g)
# \cdot
⋅(f, g) = K_dot_32_Ks_Ks(f, g)
# \times
×(f, g) = K3_cross_32_K3_K3(f, g)

# --- 2/3-rule ---
# \circledast
⊛(f, g) = K_prod_23_K_K(f, g)
# \odot
⊙(f, g) = K_dot_23_Ks_Ks(f, g)
# \otimes
⊗(f, g) = K3_cross_23_K3_K3(f, g)

max_stars_repo_name: mat-der-D/SpectralMethodFFTW

max_issues_repo_name: mat-der-D/SpectralMethodFFTW

max_forks_repo_name: mat-der-D/SpectralMethodFFTW


text:
@testset "assemble_stiffness_matrix.jl" begin

end

ext: jl


text:
require("ArgParse")
require("Iterators")

using ArgParse, Iterators
import JSON
push!(LOAD_PATH, dirname(@__FILE__()))

aa3letter = {
  "A" => "ALA",
  "R" => "ARG",
  "N" => "ASN",
  "D" => "ASP",
  "C" => "CYS",
  "Q" => "GLN",
  "E" => "GLU",
  "G" => "GLY",
  "H" => "HIS",
  "I" => "ILE",
  "L" => "LEU",
  "K" => "LYS",
  "M" => "MET",
  "F" => "PHE",
  "P" => "PRO",
  "S" => "SER",
  "T" => "THR",
  "W" => "TRP",
  "Y" => "TYR",
  "V" => "VAL"
}

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--E-ONE", "-o"
            help = "converter for e-one energy term"
            action = :store_true
            #required = true
        "--input-file-phi"
            help = "file with rotamer angles"
            required = true
        "--input-file-energies"
            help = "file with rotamer energies"
            required = false
        "--output-file"
            help = ".json file name for saving results"
            required = false
    end

    return parse_args(s)
end
aaNames = ["GLY", "ALA", "SER", "CYS",
  "VAL", "THR", "ILE", "PRO",
  "MET", "ASP", "ASN", "LEU",
  "LYS", "GLU", "GLN", "ARG",
  "HIS", "PHE", "TYR", "TRP",
  "CYX"]
aaIndices = Dict(collect(zip(aaNames, 1:21)))

immutable E12Data
  aa1 :: String
  aa2 :: String
  n14 :: Int
  e214 :: Array{Number, 1}
end

function readRotamerLibraryInfo(angles_input_file_name, energies_input_file_name)
  angles_file = open(angles_input_file_name, "r")
  #println(readlines(input_file)[0])
  angles = Dict{Int, Array{Float64, 1}}()
  while !eof(angles_file)
    #println(readline(input_file))
    s = readline(angles_file)
    if (int(split(s)[2]) != 0)
      angles[int(split(s)[1])] = map(float, split(readline(angles_file)))
    end
    #return
  end
  close(angles_file)

  energies_file = open(energies_input_file_name, "r")
  by_aa = Dict{String, Dict{Int, Array{Float64, 1}}}()
  while !eof(energies_file)
    info = split(readline(energies_file))
    ibin = int(info[2])
    aa = aa3letter[info[1]]
    bdeg = int(info[4])
    if !haskey(by_aa, aa)
      by_aa[aa] = Dict{Int, Array{Float64, 1}}()
    end
    if int(info[3]) == 1
      by_aa[aa][ibin] = [
        - log(float(readline(energies_file))*bdeg) for i = 1 : bdeg
          ]
    else
      for i = 1:int(info[4])
        readline(energies_file)
      end
    end
  end
  close(energies_file)

  {"phi" => angles, "energies" => by_aa}
end

function main()
    parsed_args = parse_commandline()
    data = readRotamerLibraryInfo(parsed_args["input-file-phi"], parsed_args["input-file-energies"])
    #println(data)
    output_file = open(parsed_args["output-file"], "w")
    println(output_file, JSON.json(data, 1))
    close(output_file)
end

main()


text:
# Dantzig-Wolfe Reformulation and Column Generation
# Functions for Master Problem
# Edward J. Xu
# 2019.5.25
########################################################################################################################

function setModelMas(numQ, vecP, numSub, vecSenseP, dualPen, gurobi_env, whiSolver)
    vecDualGuessPi = 100 * ones(Float64, numQ)
    vecDualGuessKappa = 0 * ones(Float64, numSub)
    ## X1: initial matrix of extreme points
    # if whiSolver == 1
    #     modMas =  Model(solver = GurobiSolver(OutputFlag = 0, gurobi_env))
    # elseif whiSolver == 2
    #     modMas =  Model(solver = CplexSolver(CPX_PARAM_SCRIND=0))
    # else
    modMas =  Model(with_optimizer(GLPK.Optimizer))
    # end
    K = 1
    # In this case we do not use a starting set of extreme points.
    # we just use a dummy starting column
    @variable(modMas, vecLambda[1:K] >= 0 )
    @variable(modMas, 0 <= vecMuMinus[1:numQ] <= dualPen)
    @variable(modMas, - dualPen <= vecMuPlus[1:numQ] <= 0)
    @variable(modMas, 0 <= vecMuMinusConv[1:numSub] <= dualPen)
    @variable(modMas, - dualPen <= vecMuPlusConv[1:numSub] <= 0)
    # Remember to consider if we need to maximize or minimize
    @objective(
        modMas, Max, sum(- 1000 * vecLambda[j] for j = 1:K) +
        sum(vecDualGuessPi[i] * (vecMuMinus[i] + vecMuPlus[i]) for i = 1:numQ) +
        sum(vecDualGuessKappa[k] * (vecMuMinusConv[k] + vecMuPlusConv[k]) for k=1:numSub)
        )
    # @constraint(modMas, vecConsRef[i = 1: numQ], sum(vecP[i] * vecLambda[j] for j = 1:K) == vecP[i])
    JuMPConstraintRef = JuMP.ConstraintRef{Model, JuMP.GenericRangeConstraint{JuMP.GenericAffExpr{Float64,Variable}}}
    vecConsRef = Vector{JuMPConstraintRef}(undef, numQ)
    for i = 1:numQ
        if vecSenseP[i] == leq
            vecConsRef[i] = @constraint(modMas, sum(vecP[i] * vecLambda[j] for j = 1:K) +
                vecMuMinus[i] + vecMuPlus[i] <= vecP[i])
        elseif vecSenseP[i] == geq
            vecConsRef[i] = @constraint(modMas, sum(vecP[i] * vecLambda[j] for j = 1:K) +
                vecMuMinus[i] + vecMuPlus[i] >= vecP[i])
        else
            vecConsRef[i] = @constraint(modMas, sum(vecP[i] * vecLambda[j] for j = 1:K) +
                vecMuMinus[i] + vecMuPlus[i] == vecP[i])
        end
    end
    @constraint(modMas, vecConsConvex[k = 1: numSub], sum(vecLambda[j] for j = 1: K) +
        vecMuMinusConv[k] + vecMuPlusConv[k] == 1)
    return (modMas, vecConsRef, vecConsConvex, vecLambda, vecMuMinus, vecMuPlus, vecMuMinusConv, vecMuPlusConv)
end

function addColToMas(modMas::JuMP.Model, modSub::ModelSub, vec_x_result, vecConsRef, consConvex)
    (m, n) = size(modSub.mat_e)
    A0x = modSub.mat_e * vec_x_result
    vecConsTouched = ConstraintRef[]
    vals = Float64[]
    for i = 1: m
        # only insert non-zero elements (this saves memory and may make the modMas problem easier to optimize!)
        if A0x[i] != 0
            push!(vecConsTouched, vecConsRef[i])
            push!(vals, A0x[i])
        end
    end
    # add variable to convexity constraint.
    push!(vecConsTouched, consConvex)
    push!(vals, 1)
    objCoefNew = sum(modSub.vec_l[j] * vec_x_result[j] for j = 1: n)
    println("objCoefNew = $(objCoefNew).")
    @variable(
        modMas,
        lambdaNew >= 0,                     # New variable to be added
        objective = objCoefNew,             # cost coefficient of new varaible in the objective
        inconstraints = vecConsTouched,     # constraints to be modified
        coefficients = vals                 # the coefficients of the variable in those constraints
        )
    return (lambdaNew, objCoefNew)
end

function solveMas(modMas, vecConsRef, vecConsConvex)
    status = optimize!(modMas)
    if status != :Optimal
        throw("Error: Non-optimal modMas-problem status")
    end
    vecPi = dual(vecConsRef)
    vecKappa = dual(vecConsConvex)
    obj_master = objective_value(modMas)
    ## ensure that vecPi and vecKappa are row vectors
    vecPi = reshape(vecPi, 1, length(vecPi))
    vecKappa = reshape(vecKappa, 1, length(vecKappa))
    return (vecPi, vecKappa, obj_master)
end


text:
using LightGraphs.Experimental, LightGraphs.Experimental.ShortestPaths
using SparseArrays

import Base.==
function ==(a::ShortestPaths.AStarResult, b::ShortestPaths.AStarResult)
   return a.path == b.path && a.dist == b.dist
end
@testset "Shortest Paths" begin

    g1 = path_graph(5)
    g2 = path_digraph(5)
    s1 = shortest_paths(g1, 1) # bfs
    s2 = shortest_paths(g2, 1) # bfs
    @test paths(s1)[2] == paths(s1, 2) == [1, 2]
    @test paths(s2)[2] == paths(s2, 2) == [1, 2]
    @test dists(s1)[2] == dists(s1, 2) == 1
    @test dists(s2)[2] == dists(s2, 2) == 1
    @test !has_negative_weight_cycle(g1)
    @test !has_negative_weight_cycle(g2)

    @testset "AStar" begin
        g3 = path_graph(5)
        g4 = path_digraph(5)

        d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
        d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
        for g in testgraphs(g3), dg in testdigraphs(g4)
            y = @inferred(shortest_paths(g, 1, 4, d1, AStar()))
            @test y == @inferred(shortest_paths(dg, 1, 4, d1, AStar())) ==
            @inferred(shortest_paths(g, 1, 4, d2, AStar()))
            @test paths(y) == [y.path]
            @test_throws ArgumentError paths(y, 2)
            @test dists(y) == [[y.dist]]
            @test_throws ArgumentError dists(y, 2)

            z = @inferred(shortest_paths(dg, 4, 1, AStar()))
            @test isempty(z.path)
            @test isempty(paths(z)[1])
            @test dists(z)[1] == [typemax(Int)]
        end
    end

    @testset "BellmanFord" begin
	    g4 = path_digraph(5)

	    d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
	    d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
	    for g in testdigraphs(g4)
            y = @inferred(shortest_paths(g, 2, d1, BellmanFord()))
            z = @inferred(shortest_paths(g, 2, d2, BellmanFord()))
            @test y.dists == z.dists == [Inf, 0, 6, 17, 33]
            @test @inferred(paths(z))[2] == []
            @test @inferred(paths(z))[4] == paths(z, 4) == [2, 3, 4]
            @test @inferred(!has_negative_weight_cycle(g, BellmanFord()))
            @test @inferred(!has_negative_weight_cycle(g, d1, BellmanFord()))

            y = @inferred(shortest_paths(g, 2, d1, BellmanFord()))
            z = @inferred(shortest_paths(g, 2, d2, BellmanFord()))
            @test dists(y) == dists(z) == [Inf, 0, 6, 17, 33]
            @test @inferred(paths(z))[2] == []
            @test @inferred(paths(z))[4] == paths(z, 4) == [2, 3, 4]
            @test @inferred(!has_negative_weight_cycle(g, BellmanFord()))
            z = @inferred(shortest_paths(g, 2, BellmanFord()))
            @test dists(z) == [typemax(Int), 0, 1, 2, 3]
        end

        # Negative Cycle
        gx = complete_graph(3)
        for g in testgraphs(gx)
            d = [1 -3 1; -3 1 1; 1 1 1]
            @test_throws NegativeCycleError shortest_paths(g, 1, d, BellmanFord())
            @test has_negative_weight_cycle(g, d, BellmanFord())

            d = [1 -1 1; -1 1 1; 1 1 1]
            @test_throws NegativeCycleError shortest_paths(g, 1, d, BellmanFord())
            @test has_negative_weight_cycle(g, d, BellmanFord())
            @test has_negative_weight_cycle(g, d)
        end

        # Negative cycle of length 3 in graph of diameter 4
        gx = complete_graph(4)
        d = [1 -1 1 1; 1 1 1 -1; 1 1 1 1; 1 1 1 1]
        for g in testgraphs(gx)
            @test_throws NegativeCycleError shortest_paths(g, 1, d, BellmanFord())
            @test has_negative_weight_cycle(g, d)
        end
    end

    @testset "BFS" begin
        g1= path_graph(5)
        g2 = path_digraph(5)
        add_edge!(g1, 2, 5)
        add_edge!(g2, 2, 5)
        for g in testgraphs(g1, g2)
            d = shortest_paths(g, 1, Dijkstra())
            b = shortest_paths(g, 1, BFS())
            q = shortest_paths(g, 1, BFS(Base.Sort.QuickSort))
            @test dists(d) == dists(b) && paths(d) == paths(b)
            @test dists(b) == dists(q) && paths(b) == paths(q)
            d2 = shortest_paths(g, [1, 3], Dijkstra())
            b2 = shortest_paths(g, [1, 3], BFS(ShortestPaths.NOOPSort))
            q2 = shortest_paths(g, [1, 3], BFS(Base.Sort.MergeSort))
            @test dists(d2) == dists(b2) && paths(d2) == paths(b2)
            @test dists(b2) == dists(q2) && paths(b2) == paths(q2)
        end

        @test shortest_paths(g1, 1) isa ShortestPaths.BFSResult
        @test shortest_paths(g2, 1) isa ShortestPaths.BFSResult
        m1 = shortest_paths(g1, [1, 2])
        m2 = shortest_paths(g2, [1, 2])
        @test m1 isa ShortestPaths.BFSResult
        @test m2 isa ShortestPaths.BFSResult
        @test m1.dists == [0, 0, 1, 2, 1]
        @test m2.dists == [0, 0, 1, 2, 1]
    end

    @testset "DEsopoPape" begin
        g4 = path_digraph(5)
        d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
        d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
        @testset "generic tests: $g" for g in testdigraphs(g4)
            y = @inferred(shortest_paths(g, 2, d1, DEsopoPape()))
            z = @inferred(shortest_paths(g, 2, d2, DEsopoPape()))
            @test y.parents == z.parents == [0, 0, 2, 3, 4]
            @test y.dists == z.dists == [Inf, 0, 6, 17, 33]
        end

        gx = path_graph(5)
        add_edge!(gx, 2, 4)
        d = ones(Int, 5, 5)
        d[2, 3] = 100
        @testset "cycles: $g" for g in testgraphs(gx)
            z = @inferred(shortest_paths(g, 1, d, DEsopoPape()))
            @test z.dists == [0, 1, 3, 2, 3]
            @test z.parents == [0, 1, 4, 2, 4]
        end

        m = [0 2 2 0 0; 2 0 0 0 3; 2 0 0 1 2;0 0 1 0 1;0 3 2 1 0]
        G = SimpleGraph(5)
        add_edge!(G, 1, 2)
        add_edge!(G, 1, 3)
        add_edge!(G, 2, 5)
        add_edge!(G, 3, 5)
        add_edge!(G, 3, 4)
        add_edge!(G, 4, 5)

        @testset "more cycles: $g" for g in testgraphs(G)
            y = @inferred(shortest_paths(g, 1, m, DEsopoPape()))
            @test y.parents == [0, 1, 1, 3, 3]
            @test y.dists == [0, 2, 2, 3, 4]
        end

        G = SimpleGraph(5)
        add_edge!(G, 2, 2)
        add_edge!(G, 1, 2)
        add_edge!(G, 1, 3)
        add_edge!(G, 3, 3)
        add_edge!(G, 1, 5)
        add_edge!(G, 2, 4)
        add_edge!(G, 4, 5)
        m = [0 10 2 0 15; 10 9 0 1 0; 2 0 1 0 0; 0 1 0 0 2; 15 0 0 2 0]
        @testset "self loops: $g" for g in testgraphs(G)
            z = @inferred(shortest_paths(g, 1, m, DEsopoPape()))
            y = @inferred(shortest_paths(g, 1, m, Dijkstra()))
            @test isapprox(z.dists, y.dists)
        end

        G = SimpleGraph(5)
        add_edge!(G, 1, 2)
        add_edge!(G, 1, 3)
        add_edge!(G, 4, 5)
        inf = typemax(eltype(G))
        @testset "disconnected: $g" for g in testgraphs(G)
            z = @inferred(shortest_paths(g, 1, DEsopoPape()))
            @test z.dists == [0, 1, 1, inf, inf]
            @test z.parents == [0, 1, 1, 0, 0]
        end

        G = SimpleGraph(3)
        inf = typemax(eltype(G))
        @testset "empty: $g" for g in testgraphs(G)
            z = @inferred(shortest_paths(g, 1, DEsopoPape()))
            @test z.dists == [0, inf, inf]
            @test z.parents == [0, 0, 0]
        end

        @testset "random simple graphs" begin
            for i = 1:5
                nvg = Int(ceil(250*rand()))
                neg = Int(floor((nvg*(nvg-1)/2)*rand()))
                seed = Int(floor(100*rand()))
                g = SimpleGraph(nvg, neg; seed = seed)
                z = shortest_paths(g, 1, DEsopoPape())
                y = shortest_paths(g, 1, Dijkstra())
                @test isapprox(z.dists, y.dists)
            end
        end

        @testset "random simple digraphs" begin
            for i = 1:5
                nvg = Int(ceil(250*rand()))
                neg = Int(floor((nvg*(nvg-1)/2)*rand()))
                seed = Int(floor(100*rand()))
                g = SimpleDiGraph(nvg, neg; seed = seed)
                z = shortest_paths(g, 1, DEsopoPape())
                y = shortest_paths(g, 1, Dijkstra())
                @test isapprox(z.dists, y.dists)
            end
        end

        @testset "misc graphs" begin
            G = complete_graph(9)
            z = shortest_paths(G, 1, DEsopoPape())
            y = shortest_paths(G, 1, Dijkstra())
            @test isapprox(z.dists, y.dists)

            G = complete_digraph(9)
            z = shortest_paths(G, 1, DEsopoPape())
            y = shortest_paths(G, 1, Dijkstra())
            @test isapprox(z.dists, y.dists)

            G = cycle_graph(9)
            z = shortest_paths(G, 1, DEsopoPape())
            y = shortest_paths(G, 1, Dijkstra())
            @test isapprox(z.dists, y.dists)

            G = cycle_digraph(9)
            z = shortest_paths(G, 1, DEsopoPape())
            y = shortest_paths(G, 1, Dijkstra())
            @test isapprox(z.dists, y.dists)

            G = star_graph(9)
            z = shortest_paths(G, 1, DEsopoPape())
            y = shortest_paths(G, 1, Dijkstra())
            @test isapprox(z.dists, y.dists)

            G = wheel_graph(9)
            z = shortest_paths(G, 1, DEsopoPape())
            y = shortest_paths(G, 1, Dijkstra())
            @test isapprox(z.dists, y.dists)

            G = roach_graph(9)
            z = shortest_paths(G, 1, DEsopoPape())
            y = shortest_paths(G, 1, Dijkstra())
            @test isapprox(z.dists, y.dists)

            G = clique_graph(5, 19)
            z = shortest_paths(G, 1, DEsopoPape())
            y = shortest_paths(G, 1, Dijkstra())
            @test isapprox(z.dists, y.dists)
        end

        @testset "smallgraphs: $s" for s in [
            :bull, :chvatal, :cubical, :desargues,
            :diamond, :dodecahedral, :frucht, :heawood,
            :house, :housex, :icosahedral, :krackhardtkite, :moebiuskantor,
            :octahedral, :pappus, :petersen, :sedgewickmaze, :tutte,
            :tetrahedral, :truncatedcube, :truncatedtetrahedron,
            :truncatedtetrahedron_dir
         ]
            G = smallgraph(s)
            z = shortest_paths(G, 1, DEsopoPape())
            y = shortest_paths(G, 1, Dijkstra())
            @test isapprox(z.dists, y.dists)
        end

        @testset "errors" begin
            g = Graph()
            @test_throws DomainError shortest_paths(g, 1, DEsopoPape())
            g = Graph(5)
            @test_throws DomainError shortest_paths(g, 6, DEsopoPape())
        end
    end

    @testset "Dijkstra" begin
        g4 = path_digraph(5)
        d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
        d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))

        for g in testdigraphs(g4)
            y = @inferred(shortest_paths(g, 2, d1, Dijkstra()))
            z = @inferred(shortest_paths(g, 2, d2, Dijkstra()))

            @test y.parents == z.parents == [0, 0, 2, 3, 4]
            @test y.dists == z.dists == [Inf, 0, 6, 17, 33]

            y = @inferred(shortest_paths(g, 2, d1, Dijkstra(all_paths=true)))
            z = @inferred(shortest_paths(g, 2, d2, Dijkstra(all_paths=true)))
            @test z.predecessors[3] == y.predecessors[3] == [2]

            @test @inferred(paths(z)) == paths(y)
            @test @inferred(paths(z))[4] == paths(z, 4) == paths(y, 4) == [2, 3, 4]

            @test shortest_paths(g, 1, d1) isa ShortestPaths.DijkstraResult
        end

        gx = path_graph(5)
        add_edge!(gx, 2, 4)
        d = ones(Int, 5, 5)
        d[2, 3] = 100
        for g in testgraphs(gx)
            z = @inferred(shortest_paths(g, 1, d, Dijkstra()))
            @test z.dists == [0, 1, 3, 2, 3]
            @test z.parents == [0, 1, 4, 2, 4]
        end

        # small function to reconstruct the shortest path; I copied it from somewhere,
        # can't find the original source to give the credits
        # @Beatzekatze on github
        spath(target, dijkstraStruct, source) = target == source ? target : [spath(dijkstraStruct.parents[target], dijkstraStruct, source) target]
        spaths(ds, targets, source) = [spath(i, ds, source) for i in targets]

        G = Graph(4)
        add_edge!(G, 2, 1)
        add_edge!(G, 2, 3)
        add_edge!(G, 1, 4)
        add_edge!(G, 3, 4)
        add_edge!(G, 2, 2)
        w = [0. 3. 0. 1.;
            3. 0. 2. 0.;
            0. 2. 0. 3.;
            1. 0. 3. 0.]

        for g in testgraphs(G)
            ds = @inferred(shortest_paths(g, 2, w, Dijkstra()) )
          # this loop reconstructs the shortest path for vertices 1, 3 and 4
            @test spaths(ds, [1, 3, 4], 2) == Array[[2 1],
                                                [2 3],
                                                [2 1 4]]

        # here a selflink at source is introduced; it should not change the shortest paths
            w[2, 2] = 10.0
            ds = @inferred(shortest_paths(g, 2, w, Dijkstra()))
          # this loop reconstructs the shortest path for vertices 1, 3 and 4
            @test spaths(ds, [1, 3, 4], 2) == Array[[2 1],
                                                [2 3],
                                                [2 1 4]]
        end

        #615
        m = [0 2 2 0 0; 2 0 0 0 3; 2 0 0 1 2;0 0 1 0 1;0 3 2 1 0]
        G = SimpleGraph(5)
        add_edge!(G, 1, 2)
        add_edge!(G, 1, 3)
        add_edge!(G, 2, 5)
        add_edge!(G, 3, 5)
        add_edge!(G, 3, 4)
        add_edge!(G, 4, 5)
        for g in testgraphs(G)
            ds = @inferred(shortest_paths(g, 1, m, Dijkstra(all_paths=true)))
            @test ds.pathcounts   == [1, 1, 1, 1, 2]
            @test ds.predecessors == [[], [1], [1], [3], [3, 4]]
            @test ds.predecessors == [[], [1], [1], [3], [3, 4]]

            dm = @inferred(shortest_paths(g, 1, Dijkstra(all_paths=true, track_vertices=true)))
            @test dm.pathcounts       == [1, 1, 1, 1, 2]
            @test dm.predecessors     == [[], [1], [1], [3], [2, 3]]
            @test dm.closest_vertices == [1, 2, 3, 5, 4]
        end

        G = SimpleGraph(5)
        add_edge!(G, 1, 2)
        add_edge!(G, 1, 3)
        add_edge!(G, 4, 5)
        for g in testgraphs(G)
            dm = @inferred(shortest_paths(g, 1, Dijkstra(all_paths=true, track_vertices=true)))
            @test dm.closest_vertices == [1, 2, 3, 4, 5]
        end
    end

    @testset "FloydWarshall" begin
        g3 = path_graph(5)
        d = [0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]
        for g in testgraphs(g3)
            z = @inferred(shortest_paths(g, d, FloydWarshall()))
            @test z.dists[3, :][:] == [7, 6, 0, 11, 27]
            @test z.parents[3, :][:] == [2, 3, 0, 3, 4]

            @test @inferred(paths(z))[2][2] == []
            @test @inferred(paths(z))[2][4] == paths(z, 2)[4] == paths(z, 2, 4) == [2, 3, 4]
        end
        g4 = path_digraph(4)
        d = ones(4, 4)
        for g in testdigraphs(g4)
            z = @inferred(shortest_paths(g, d, FloydWarshall()))
            @test length(paths(z, 4, 3)) == 0
            @test length(paths(z, 4, 1)) == 0
            @test length(paths(z, 2, 3)) == 2
        end

        g5 = DiGraph([1 1 1 0 1; 0 1 0 1 1; 0 1 1 0 0; 1 0 1 1 0; 0 0 0 1 1])
        d = [0 3 8 0 -4; 0 0 0 1 7; 0 4 0 0 0; 2 0 -5 0 0; 0 0 0 6 0]
        for g in testdigraphs(g5)
            z = @inferred(shortest_paths(g, d, FloydWarshall()))
            @test z.dists == [0 1 -3 2 -4; 3 0 -4 1 -1; 7 4 0 5 3; 2 -1 -5 0 -2; 8 5 1 6 0]
        end

        @testset "paths infinite loop bug" begin
            g = SimpleGraph(2)
            add_edge!(g, 1, 2)
            add_edge!(g, 2, 2)
            @test paths(shortest_paths(g, FloydWarshall())) ==
                Vector{Vector{Int}}[[[], [1, 2]], [[2, 1], []]]

            g = SimpleDiGraph(2)
            add_edge!(g, 1, 1)
            add_edge!(g, 1, 2)
            add_edge!(g, 2, 1)
            add_edge!(g, 2, 2)
            @test paths(shortest_paths(g, FloydWarshall())) ==
                Vector{Vector{Int}}[[[], [1, 2]], [[2, 1], []]]
        end
        @testset "default for no source, no alg" begin
            g = path_graph(4)
            w = zeros(4,4)
            for i in 1:3
                w[i, i+1] = 1.0
                w[i+1, i] = 1.0
            end
            @test shortest_paths(g) isa ShortestPaths.FloydWarshallResult
            @test shortest_paths(g, w) isa ShortestPaths.FloydWarshallResult
        end
    end

    @testset "Johnson" begin
        g3 = path_graph(5)
        d = Symmetric([0 1 2 3 4; 1 0 6 7 8; 2 6 0 11 12; 3 7 11 0 16; 4 8 12 16 0])
        for g in testgraphs(g3)
            z = @inferred(shortest_paths(g, d, Johnson()))
            @test z.dists[3, :][:] == [7, 6, 0, 11, 27]
            @test z.parents[3, :][:] == [2, 3, 0, 3, 4]

            @test @inferred(paths(z))[2][2] == []
            @test @inferred(paths(z))[2][4] == paths(z, 2)[4] == paths(z, 2, 4) == [2, 3, 4]
        end

        g4 = path_digraph(4)
        for g in testdigraphs(g4)
            z = @inferred(shortest_paths(g, Johnson()))
            @test length(paths(z, 4, 3)) == 0
            @test length(paths(z, 4, 1)) == 0
            @test length(paths(z, 2, 3)) == 2
        end

        g5 = DiGraph([1 1 1 0 1; 0 1 0 1 1; 0 1 1 0 0; 1 0 1 1 0; 0 0 0 1 1])
        d = [0 3 8 0 -4; 0 0 0 1 7; 0 4 0 0 0; 2 0 -5 0 0; 0 0 0 6 0]
        for g in testdigraphs(g5)
            z = @inferred(shortest_paths(g, d, Johnson()))
            @test z.dists == [0 1 -3 2 -4; 3 0 -4 1 -1; 7 4 0 5 3; 2 -1 -5 0 -2; 8 5 1 6 0]
        end
    end

    @testset "SPFA" begin
        @testset "Generic tests for graphs" begin
            g4 = path_digraph(5)
            d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])

            for g in testdigraphs(g4)
                y = @inferred(shortest_paths(g, 2, d1, SPFA()))
                @test y.dists == [Inf, 0, 6, 17, 33]
            end

            @testset "Graph with cycle" begin
                gx = path_graph(5)
                add_edge!(gx, 2, 4)
                d = ones(Int, 5, 5)
                d[2, 3] = 100
                for g in testgraphs(gx)
                    z = @inferred(shortest_paths(g, 1, d, SPFA()))
                    @test z.dists == [0, 1, 3, 2, 3]
                end
            end

            m = [0 2 2 0 0; 2 0 0 0 3; 2 0 0 1 2; 0 0 1 0 1; 0 3 2 1 0]
            G = SimpleGraph(5)
            add_edge!(G, 1, 2)
            add_edge!(G, 1, 3)
            add_edge!(G, 2, 5)
            add_edge!(G, 3, 5)
            add_edge!(G, 3, 4)
            add_edge!(G, 4, 5)

            for g in testgraphs(G)
                y = @inferred(shortest_paths(g, 1, m, SPFA()))
                @test y.dists == [0, 2, 2, 3, 4]
            end
        end

        @testset "Graph with self loop" begin
            G = SimpleGraph(5)
            add_edge!(G, 2, 2)
            add_edge!(G, 1, 2)
            add_edge!(G, 1, 3)
            add_edge!(G, 3, 3)
            add_edge!(G, 1, 5)
            add_edge!(G, 2, 4)
            add_edge!(G, 4, 5)
            m = [0 10 2 0 15; 10 9 0 1 0; 2 0 1 0 0; 0 1 0 0 2; 15 0 0 2 0]
            for g in testgraphs(G)
                z = @inferred(shortest_paths(g, 1 , m, SPFA()))
                y = @inferred(dijkstra_shortest_paths(g, 1, m))
                @test isapprox(z.dists, y.dists)
            end
        end

        @testset "Disconnected graph" begin
            G = SimpleGraph(5)
            add_edge!(G, 1, 2)
            add_edge!(G, 1, 3)
            add_edge!(G, 4, 5)
            inf = typemax(eltype(G))
            for g in testgraphs(G)
                z = @inferred(shortest_paths(g, 1, SPFA()))
                @test z.dists == [0, 1, 1, inf, inf]
            end
        end

        @testset "Empty graph" begin
            G = SimpleGraph(3)
            inf = typemax(eltype(G))
            for g in testgraphs(G)
                z = @inferred(shortest_paths(g, 1, SPFA()))
                @test z.dists == [0, inf, inf]
            end
        end

        @testset "Random Graphs" begin
            @testset "Simple graphs" begin
                for i = 1:5
                    nvg = Int(ceil(250*rand()))
                    neg = Int(floor((nvg*(nvg-1)/2)*rand()))
                    seed = Int(floor(100*rand()))
                    g = SimpleGraph(nvg, neg; seed = seed)
                    z = shortest_paths(g, 1, SPFA())
                    y = dijkstra_shortest_paths(g, 1)
                    @test isapprox(z.dists, y.dists)
                end
            end

            @testset "Simple DiGraphs" begin
                for i = 1:5
                    nvg = Int(ceil(250*rand()))
                    neg = Int(floor((nvg*(nvg-1)/2)*rand()))
                    seed = Int(floor(100*rand()))
                    g = SimpleDiGraph(nvg, neg; seed = seed)
                    z = shortest_paths(g, 1, SPFA())
                    y = dijkstra_shortest_paths(g, 1)
                    @test isapprox(z.dists, y.dists)
                end
            end
        end

        @testset "Different types of graph" begin
            G = complete_graph(9)
            z = shortest_paths(G, 1, SPFA())
            y = dijkstra_shortest_paths(G, 1)
            @test isapprox(z.dists, y.dists)

            G = complete_digraph(9)
            z = shortest_paths(G, 1, SPFA())
            y = dijkstra_shortest_paths(G, 1)
            @test isapprox(z.dists, y.dists)

            G = cycle_graph(9)
            z = shortest_paths(G, 1, SPFA())
            y = dijkstra_shortest_paths(G, 1)
            @test isapprox(z.dists, y.dists)

            G = cycle_digraph(9)
            z = shortest_paths(G, 1, SPFA())
            y = dijkstra_shortest_paths(G, 1)
            @test isapprox(z.dists, y.dists)

            G = star_graph(9)
            z = shortest_paths(G, 1, SPFA())
            y = dijkstra_shortest_paths(G, 1)
            @test isapprox(z.dists, y.dists)

            G = wheel_graph(9)
            z = shortest_paths(G, 1, SPFA())
            y = dijkstra_shortest_paths(G, 1)
            @test isapprox(z.dists, y.dists)

            G = roach_graph(9)
            z = shortest_paths(G, 1, SPFA())
            y = dijkstra_shortest_paths(G, 1)
            @test isapprox(z.dists, y.dists)

            G = clique_graph(5, 19)
            z = shortest_paths(G, 1, SPFA())
            y = dijkstra_shortest_paths(G, 1)
            @test isapprox(z.dists, y.dists)

            @testset "Small Graphs" begin
                for s in [:bull, :chvatal, :cubical, :desargues,
                          :diamond, :dodecahedral, :frucht, :heawood,
                          :house, :housex, :icosahedral, :krackhardtkite, :moebiuskantor,
                          :octahedral, :pappus, :petersen, :sedgewickmaze, :tutte,
                          :tetrahedral, :truncatedcube, :truncatedtetrahedron, :truncatedtetrahedron_dir]
                    G = smallgraph(s)
                    z = shortest_paths(G, 1, SPFA())
                    y = dijkstra_shortest_paths(G, 1)
                    @test isapprox(z.dists, y.dists)
                end
            end
        end

        @testset "Normal graph" begin
            g4 = path_digraph(5)

            d1 = float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
            d2 = sparse(float([0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
            for g in testdigraphs(g4)
                y = @inferred(shortest_paths(g, 2, d1, SPFA()))
                z = @inferred(shortest_paths(g, 2, d2, SPFA()))
                @test y.dists == z.dists == [Inf, 0, 6, 17, 33]
                @test @inferred(!has_negative_weight_cycle(g, SPFA()))
                @test @inferred(!has_negative_weight_cycle(g, d1, SPFA()))

                y = @inferred(shortest_paths(g, 2, d1, SPFA()))
                z = @inferred(shortest_paths(g, 2, d2, SPFA()))
                @test y.dists == z.dists == [Inf, 0, 6, 17, 33]
                @test @inferred(!has_negative_weight_cycle(g, SPFA()))
                z = @inferred(shortest_paths(g, 2, SPFA()))
                @test z.dists == [typemax(Int), 0, 1, 2, 3]
            end
        end

        @testset "Negative Cycle" begin
            # Negative Cycle 1
            gx = complete_graph(3)
            for g in testgraphs(gx)
                d = [1 -3 1; -3 1 1; 1 1 1]
                @test_throws ShortestPaths.NegativeCycleError shortest_paths(g, 1, d, SPFA())
                @test has_negative_weight_cycle(g, d, SPFA())

                d = [1 -1 1; -1 1 1; 1 1 1]
                @test_throws ShortestPaths.NegativeCycleError shortest_paths(g, 1, d, SPFA())
                @test has_negative_weight_cycle(g, d, SPFA())
            end

            # Negative cycle of length 3 in graph of diameter 4
            gx = complete_graph(4)
            d = [1 -1 1 1; 1 1 1 -1; 1 1 1 1; 1 1 1 1]
            for g in testgraphs(gx)
                @test_throws ShortestPaths.NegativeCycleError shortest_paths(g, 1, d, SPFA())
                @test has_negative_weight_cycle(g, d, SPFA())
            end
        end

    end
end


text:
import ParserCombinator

function convert_results(;sub_path::String="", specific_run::String="")

    if specific_run != ""
        raw_data = load(joinpath("results", specific_run))
        data = raw_data[first(keys(raw_data))]
        filename = specific_run[1:first(findfirst(".jld2", specific_run))-1]

        if !in(filename, readdir("dataexchange"))
                mkdir(joinpath("dataexchange", filename))
        end

        CSV.write(
            joinpath("dataexchange", filename, "agent_log" * ".csv"),
            data.agent_log
        )

        CSV.write(
            joinpath("dataexchange", filename, "post_log" * ".csv"),
            data.post_log
        )

        for i in 1:length(data.graph_list)
            graphnr = lpad(
                string(i),
                length(string(length(data.graph_list))),
                "0"
            )

            savegraph(
                joinpath("dataexchange", filename, "graph_$graphnr.gml"),
                data.graph_list[i],
                GraphIO.GML.GMLFormat()
            )
        end
    else

        if !in("dataexchange", readdir())
                mkdir("dataexchange")
        end

        if sub_path != ""
            path = joinpath("results", sub_path)
        else
            path = "results"
        end

        for file in readdir(path)

            if !occursin("jld2", file)
                continue
            end

            raw_data = load(joinpath(path, file))
            data = raw_data[first(keys(raw_data))]
            filename = file[1:first(findfirst(".jld2", file))-1]

            if !in(filename, readdir("dataexchange"))
                    mkdir(joinpath("dataexchange", filename))
            end

            CSV.write(
                joinpath("dataexchange", filename, "agent_log" * ".csv"),
                data.agent_log
            )

            CSV.write(
                joinpath("dataexchange", filename, "post_log" * ".csv"),
                data.post_log
            )

            savegraph(
                joinpath("dataexchange", filename, "graph_final.gml"),
                data.final_state[1],
                GraphIO.GML.GMLFormat()
            )

            # for i in 1:length(data.graph_list)
            #     graphnr = lpad(
            #         string(i),
            #         length(string(length(data.graph_list))),
            #         "0"
            #     )
            #
            #     savegraph(
            #         joinpath("dataexchange", filename, "graph_$graphnr.gml"),
            #         data.graph_list[i],
            #         GraphIO.GML.GMLFormat()
            #     )
            # end
        end
    end
end


text:
# Start here to test the policy-compliance model
using Pkg
Pkg.activate(".")

# load other packages
using Agents

# Load submodules of the simulations
include(joinpath("..", "src", "PolicyCompliance.jl"))

include("../src/simulation/simulation.jl")

# run the actual simulation

abm = initialize_model()
(data, mdata) = runsimulation(abm)

#data_export tbd


text:
include("../../src/ToyModelsMD.jl")
using .ToyModelsMD
using DelimitedFiles
using Plots

function Afunc(system::AbstractSystem,mm::MixedModel,ind::Int64,dt::Float64)
  lam::Float64 = 2.0e3
  pot::Float64 = potential(system,mm.potentials[1])
  ind==1500 && system.x[2]>0.7 && pot<-145.0 ? lam : 0.0
end

function cbfunc(xstore::Array,system::AbstractSystem,mm::MixedModel,ind::Int64,dt::Float64)
  @. @views xstore[ind,:] .= copy(system.x)
end

### define system ###
dt = 0.0001
steps = 1500
tf = steps*dt

# mueller-brown potential
mb = MullerBrown()

# gaussian force model
Mx = 21
My = 21
Mt = 21
shape = (Mx,My,Mt,2)
# Define gaussian basis centers and variance
mux = zeros(Float64,Mx)
muy = zeros(Float64,My)
mut = zeros(Float64,Mt)
for j=1:Mx
  mux[j] = -1.5+(j-1)*(1.5-(-1.5))/(Mx-1)
end
for j=1:My
  muy[j] = -0.5+(j-1)*(2-(-0.5))/(My-1)
end
for j in 1:Mt
  mut[j] = 1+(j*steps/(Mt-1)) # 1,1+2sigamt,1+4sigmat,...,1+2(Mt-1)sigmat which is steps+1
  mut[j] = (j*steps/(Mt-1)) # 1,1+2sigamt,1+4sigmat,...,1+2(Mt-1)sigmat which is steps+1
  mut[j] *= dt
end
sigx = zeros(Float64,Mx)
sigy = zeros(Float64,My)
sigt = zeros(Float64,Mt)
sigx[:] .= ((1.5-(-1.5))/(Mx-1)/2.0)^2 #variance
sigy[:] .= ((2-(-0.5))/(My-1)/2.0)^2 #variance
sigt[:] .= ((steps*1.0/(Mt-1))/2.0)^2 #variance
sigt[:] .*= dt^2
# read in optimal coefficients
x = readdlm("coeff0x.txt")
y = readdlm("coeff0y.txt")
cs = zeros(Float64,shape)
for i in 1:21
  cs[:,:,i,1] .= reshape(x[:,i],(21,21))'
  cs[:,:,i,2] .= reshape(y[:,i],(21,21))'
end
cs = reshape(cs, (21*21*21,2))
# define model - yay!
gm = GaussianModel2D(Mx,My,Mt,mux,muy,mut,sigx,sigy,sigt,cs)

# mixed model
model = MixedModel(2,[mb,gm])

# thermostat and system
langevin = Langevin(1.0,1.0)
system = ThermostattedSystem(2,0.0,[0.63,0.03],zeros(Float64,2),langevin)

# define value baseline function
shape = (Mx,My,Mt)
vmux = zeros(Float64,Mx)
vmuy = zeros(Float64,My)
vmut = zeros(Float64,Mt)
for j=1:Mx
  vmux[j] = -1.5+(j-1)*(1.5-(-1.5))/(Mx-1)
end
for j=1:My
  vmuy[j] = -0.5+(j-1)*(2-(-0.5))/(My-1)
end
for j in 1:Mt
  vmut[j] = 1+(j-1)*steps*1.0/(Mt-1) # 1,1+2sigamt,1+4sigmat,...,1+2(Mt-1)sigmat which is steps+1
  vmut[j] *= dt
end
vsigx = zeros(Float64,Mx)
vsigy = zeros(Float64,My)
vsigt = zeros(Float64,Mt)
vsigx[:] .= ((1.5-(-1.5))/(Mx-1)/2.0)^2 #variance
vsigy[:] .= ((2-(-0.5))/(My-1)/2.0)^2 #variance
vsigt[:] .= ((steps*1.0/(Mt-1))/2.0)^2 #variance
vsigt[:] .*= dt^2
# read in optimal coefficients
v = readdlm("coeff0V.txt")
vcs = zeros(Float64,shape)
for i in 1:21
  vcs[:,:,i] .= reshape(v[:,i],(21,21))'
end
vcs = reshape(vcs, (21*21*21))
vbl = GaussianValueBaseline2D(Mx,My,Mt,mux,muy,mut,sigx,sigy,sigt,vcs)

# define callback function
mcvb = MCVBCallback(Afunc,gm,vbl)
# set up integrator and run!
integrator = StochasticEuler(system,model)
runtraj!(steps,dt,integrator,mcvb)
println(mcvb.dkl)
writedlm("gradf.txt",mcvb.gradf)
writedlm("gradv.txt",mcvb.gradv*dt)


text:
module ExamplePackage

end


text:
module PackageCompilerX

using Base: active_project
using Libdl: Libdl
using Pkg: Pkg

include("juliaconfig.jl")

const CC = (Sys.iswindows() ? `x86_64-w64-mingw32-gcc` : `gcc`)

function get_julia_cmd()
    julia_path = joinpath(Sys.BINDIR, Base.julia_exename())
    image_file = unsafe_string(Base.JLOptions().image_file)
    cmd = `$julia_path -J$image_file --color=yes --startup-file=no -Cnative`
end

# Returns a vector of precompile statemenets
function run_precompilation_script(project::String, precompile_file::String)
    tracefile = tempname()
    julia_code = """Base.__init__(); include($(repr(precompile_file)))"""
    run(`$(get_julia_cmd()) --project=$project --trace-compile=$tracefile -e $julia_code`)
    return tracefile
end

function create_object_file(object_file::String, packages::Union{Symbol, Vector{Symbol}};
                            project::String=active_project(),
                            precompile_execution_file::Union{String, Nothing}=nothing,
                            precompile_statements_file::Union{String, Nothing}=nothing)
    # include all packages into the sysimage
    packages = vcat(packages)
    julia_code = """
        if !isdefined(Base, :uv_eventloop)
            Base.reinit_stdio()
        end
        Base.__init__();
        """
    for package in packages
        julia_code *= "using $package\n"
    end

    # handle precompilation
    if precompile_execution_file !== nothing || precompile_statements_file !== nothing
        precompile_statements = ""
        if precompile_execution_file !== nothing
            tracefile = run_precompilation_script(project, precompile_execution_file)
            precompile_statements *= "append!(precompile_statements, readlines($(repr(tracefile))))\n"
        end
        if precompile_statements_file != nothing
            precompile_statements *= "append!(precompile_statements, readlines($(repr(precompile_statements_file))))\n"
        end

        precompile_code = """
            # This @eval prevents symbols from being put into Main
            @eval Module() begin
                PrecompileStagingArea = Module()
                for (_pkgid, _mod) in Base.loaded_modules
                    if !(_pkgid.name in ("Main", "Core", "Base"))
                        eval(PrecompileStagingArea, :(const \$(Symbol(_mod)) = \$_mod))
                    end
                end
                precompile_statements = String[]
                $precompile_statements
                for statement in sort(precompile_statements)
                    # println(statement)
                    try
                        Base.include_string(PrecompileStagingArea, statement)
                    catch
                        # See julia issue #28808
                        @error "failed to execute \$statement"
                    end
                end
            end # module
            """
        julia_code *= precompile_code
    end

    # finally, make julia output the resulting object file
    @debug "creating object file at $object_file"
    @info "PackageCompilerX: creating object file, this might take a while..."
    run(`$(get_julia_cmd()) --project=$project --output-o=$(object_file) -e $julia_code`)
end

default_sysimage_path() = joinpath(julia_private_libdir(), "sys." * Libdl.dlext)
backup_sysimage_path() = default_sysimage_path() * ".backup"

function create_sysimage(packages::Union{Symbol, Vector{Symbol}}=Symbol[];
                         sysimage_path::Union{String,Nothing}=nothing,
                         project::String=active_project(),
                         precompile_execution_file::Union{String, Nothing}=nothing,
                         precompile_statements_file::Union{String, Nothing}=nothing,
                         replace_default_sysimage::Bool=false)
    if sysimage_path === nothing && replace_default_sysimage == false
        error("`sysimage_path` cannot be `nothing` if `replace_default_sysimage` is `false`")
    end
    if sysimage_path === nothing
        sysimage_path = string(tempname(), ".", Libdl.dlext)
    end

    object_file = tempname() * ".o"
    create_object_file(object_file, packages; project=project, precompile_execution_file=precompile_execution_file,
                       precompile_statements_file=precompile_statements_file)
    @show sysimage_path
    create_sysimage_from_object_file(object_file, sysimage_path)
    if replace_default_sysimage
        if !isfile(backup_sysimage_path())
            cp(default_sysimage_path(), backup_sysimage_path())
            @debug "making a backup of sysimage"
        end
        @info "PackageCompilerX: default sysimage replaced, restart Julia for the new sysimage to be in effect"
        cp(sysimage_path, default_sysimage_path(); force=true)
    end
    # TODO: Remove object file
end

function create_sysimage_from_object_file(input_object::String, sysimage_path::String)
    julia_libdir = dirname(Libdl.dlpath("libjulia"))

    # Prevent compiler from stripping all symbols from the shared lib.
    # TODO: On clang on windows this is called something else
    if Sys.isapple()
        o_file = `-Wl,-all_load $input_object`
    else
        o_file = `-Wl,--whole-archive $input_object -Wl,--no-whole-archive`
    end
    extra = Sys.iswindows() ? `-Wl,--export-all-symbols` : ``
    run(`$CC -v -shared -L$(julia_libdir) -o $sysimage_path $o_file -ljulia $extra`)
    return nothing
end

function restore_default_sysimage()
    if !isfile(backup_sysimage_path())
        error("did not find a backup sysimage")
    end
    cp(backup_sysimage_path(), default_sysimage_path(); force=true)
    rm(backup_sysimage_path())
    @info "PackageCompilerX: default sysimage restored, restart Julia for the new sysimage to be in effect"
    return nothing
end

# This requires that the sysimage have been built so that there is a ccallable `julia_main`
# in Main.
function create_executable_from_sysimage(;sysimage_path::String,
                                          executable_path::String)
    flags = join((cflags(), ldflags(), ldlibs()), " ")
    flags = Base.shell_split(flags)
    wrapper = joinpath(@__DIR__, "embedding_wrapper.c")
     if Sys.iswindows()
        # Cannot create an executable without copying dlls on Windows...
        rpath = `` # functionality doesn't readily exist on this platform
    elseif Sys.isapple()
        rpath = `-Wl,-rpath,@executable_path`
    else
        rpath = `-Wl,-rpath,\$ORIGIN`
    end
    extra = Sys.iswindows() ? `-Wl,--export-all-symbols` : ``
    run(`$CC -DJULIAC_PROGRAM_LIBNAME=$(repr(sysimage_path)) -o $(executable_path) $(wrapper) $sysimage_path -O2 $rpath $flags $extra`)
    return nothing
end

function create_app(app_dir::String,
                    precompile_execution_file::Union{String, Nothing}=nothing,
                    precompile_statements_file::Union{String, Nothing}=nothing)

    project_path = Pkg.Types.projectfile_path(app_dir; strict=true)
    project_toml = Pkg.TOML.parsefile(project_path)
    app_name = get(() -> error("expected package to have a name entry"), project_toml, "name")
    sysimage_file = app_name * "." * Libdl.dlext
    cd(app_dir) do
        create_sysimage(Symbol(app_name); sysimage_path=sysimage_file, project = app_dir)
        create_executable_from_sysimage(;sysimage_path=sysimage_file, executable_path=app_name)
    end
end

end # module


text:
import ThreeJS
import FileIO: load
import Reactive: every

ascent = load("assets/ascent.png")
a = reinterpret(UInt8, ascent.data)

images = cycle((a, a[257:512, 257:512], a[292:375, 355:463]))

main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"nested-props")

    i = start(images)
    map(map(t -> next(images, i), every(1))) do pair
        a, i = pair
        ThreeJS.outerdiv() <<
        (ThreeJS.initscene() <<
            [
                ThreeJS.mesh(0.0, 0.0, 0.0) <<
                [
                    ThreeJS.plane(convert(Float64, size(a, 1)), convert(Float64, size(a, 2))),
                    ThreeJS.shadermaterial(open(readall, "assets/datatexture.vert", "r"), open(readall, "assets/datatexture.frag", "r")) <<
                        ThreeJS.datatexture("data", a; :minfilter => "LinearFilter")
                ],
                ThreeJS.camera(0.0, 0.0, 768.0)
            ]
        )
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 0d696b90-14b9-11eb-19c6-7110e9eded9e
begin
	using PotentialFlow
	using Plots
	using PlutoUI
end

# ╔═╡ 97a34c00-153a-11eb-0b35-e92e0aafac5e
md"""
### Flow over a cylinder
"""

# ╔═╡ 59e4bc40-14b9-11eb-3670-cdc8fa8f0ded
begin
	X = range(-3.5, 3.5, length=100)
	Y = range(-3.5, 3.5, length=100)
	Z = [x + im*y for y in Y, x in X]
end

# ╔═╡ 682ca32e-1568-11eb-35b2-15925e38f564
md"""
### Wall on the x-axis
"""

# ╔═╡ 7acb2d90-1568-11eb-2e53-f5b7e59f6c3d
md"""
Y location of vortexes: $(@bind Δy Slider(0.0:0.1:5.0,default = 1.0))
"""

# ╔═╡ b0349b90-1565-11eb-03e8-cf9127144400
begin
	points = Vortex.Point.([Δy*im, -Δy*im], [1, -1])
	streamlines(X, Y, points,
		color = cgrad([:black, :black]),
		ratio = 1,
		legend = false,
		markerstrokealpha = 0,
		grid = false,
		xlim = (-3.5, 3.5),
		ylim = (-3.5, 3.5))
end

# ╔═╡ 492f49b0-1568-11eb-28e1-0117b5154488
md"""
### Translate and Rotate the wall
"""

# ╔═╡ ba5a0330-156a-11eb-39db-85c8b88ad474
md"""
X location of vortexes: $(@bind Δx Slider(-2.0:0.1:2.0,default = 0.0))
"""

# ╔═╡ 748be780-1568-11eb-378b-75977d19d825
begin

	fs = Freestreams.Freestream(1.0)
	db = Doublets.Doublet(0.0im, π)

	point_vortex = Vortex.Point(Δx + 1.2im, 2π)
	image_vortex = Vortex.Point(conj(1/point_vortex.z), -2π)

	ψ_levels = collect(range(-5, 5, length=31))

	images = streamlines(X, Y, (fs,db,point_vortex,image_vortex),
		color = cgrad([:black, :black]),
		ratio = 1,
		levels = ψ_levels,
		legend = false,
		grid = false,
		xlim = (-3.5, 3.5),
		ylim = (-3.5, 3.5))

	plot!(cos.(0:0.1:2π), sin.(0:0.1:2π), linestyle = :dash, linewidth = 3, color = :red)

end

# ╔═╡ 27299cc0-1569-11eb-1e3d-f7188084f732
point_vortex.z

# ╔═╡ 24496030-1569-11eb-27d6-3d2fa7417c1a
conj(point_vortex.z)

# ╔═╡ bf594c00-1566-11eb-1acd-31fb57d91e48
md"""
### System of vortexes
"""

# ╔═╡ d4d166c0-1567-11eb-07d6-23315632cff4
md"""
Γ1: $(@bind Γ₁ Slider(-2π:π/10:2π,default = 1.0))
Γ2: $(@bind Γ₂ Slider(-2π:π/10:2π,default = 1.0))
Γ3: $(@bind Γ₃ Slider(-2π:π/10:2π,default = -1.0))
Γ4: $(@bind Γ₄ Slider(-2π:π/10:2π,default = -1.0))
"""

# ╔═╡ 2e843420-1538-11eb-14b2-d9ab38fa091b
begin

	zs = [1 + im, -1 + im, -1 - im, 1 - im]

	vs = Vortex.Point.(zs, [Γ₁, Γ₂, Γ₃, Γ₄])

	streamlines(X, Y, vs,
		color = cgrad([:black, :black]),
		ratio = 1,
		legend = false,
		grid = false,
		xlim = (-3.5, 3.5),
		ylim = (-3.5, 3.5))

	plot!(vs, markersize = 10, color = :RdBu)

end

# ╔═╡ be75d460-1567-11eb-39d5-0f62c8c97fd8

# ╔═╡ Cell order:
# ╟─97a34c00-153a-11eb-0b35-e92e0aafac5e
# ╠═0d696b90-14b9-11eb-19c6-7110e9eded9e
# ╠═59e4bc40-14b9-11eb-3670-cdc8fa8f0ded
# ╟─682ca32e-1568-11eb-35b2-15925e38f564
# ╟─7acb2d90-1568-11eb-2e53-f5b7e59f6c3d
# ╠═b0349b90-1565-11eb-03e8-cf9127144400
# ╠═492f49b0-1568-11eb-28e1-0117b5154488
# ╟─ba5a0330-156a-11eb-39db-85c8b88ad474
# ╠═748be780-1568-11eb-378b-75977d19d825
# ╠═27299cc0-1569-11eb-1e3d-f7188084f732
# ╠═24496030-1569-11eb-27d6-3d2fa7417c1a
# ╟─bf594c00-1566-11eb-1acd-31fb57d91e48
# ╠═d4d166c0-1567-11eb-07d6-23315632cff4
# ╠═2e843420-1538-11eb-14b2-d9ab38fa091b
# ╠═be75d460-1567-11eb-39d5-0f62c8c97fd8


text:
module DeepRL

greet() = print("Hello Deep Reinforcement Learning!")

using Reexport, Lazy
@reexport using RLCore

export ExperienceReplay,
    WeightedExperienceReplay, OnlineReplay,
    size, getindex, add!, sample

include("replay.jl")
include("ImageBuffer.jl")

export MountainCar, Atari
include("env.jl")

export ϵGreedy, sample
include("policy.jl")

export QLearning, DoubleQLearning, AuxQLearning, TDLearning
include("update.jl")

export DQNAgent, OnlineDQNAgent, AT_DQNAgent
include("agent.jl")

include("utils.jl")

end # module


text:
# Not really a "data" macro, but very convenient nonetheless
macro col(ex...)
    if length(ex) < 2
        error("Use @col(table, colname1, ...)")
    end

    if length(ex) == 2
        table = ex[1]
        colname = Expr(:quote,ex[2])

        return esc(:($table[Val{$colname}]))
    else
        table = ex[1]
        colnames = ex[2:end]

        return esc(:($table[Val{$colnames}]))
    end
end

# Select, rename and compute columns, returning a new table
macro select(x...)
    if length(x) == 0 # Is this an error or an empty table?
        return :(error("@select expects a table"))
    end
    table = x[1]
    exprs = Vector{Any}(length(x)-1)
    for i = 2:length(x)
        expr = x[i]

        if isa(expr, Symbol)
            # Straightforward extraction (by name)
            #exprs[i-1] = :( TypedTables.extractfield(TypedTables.index($(esc(table))),Val{$(Expr(:quote,expr))}) = $(esc(table))[Val{$(Expr(:quote,expr))}] )
            exprs[i-1] = :( $(esc(expr)) = $(esc(table))[Val{$(Expr(:quote, expr))}] )
        elseif isa(expr, Expr)
            if expr.head == :(=) || expr.head == :(kw)
                expr_left = expr.args[1]
                expr_right = expr.args[2]

                # Determine operation from RHS
                operation = :unknown
                if isa(expr_right, Symbol)
                    operation = :rename
                elseif isa(expr_right, Expr)
                    if expr_right.head == :(::)
                        return(:( error("Expected syntax like @select(table, col1, newname = col2, newcol::newtype = col1 -> f(col1))") ))
                    elseif expr_right.head == :(->)
                        operation = :compute
                    end
                end

                if operation == :unknown
                    return(:( error("Expected syntax like @select(table, col1, newname = col2, newcol::newtype = col1 -> f(col1))") ))
                elseif operation == :rename
                    if isa(expr_left,Symbol)
                        #exprs[i-1] = :( TypedTables.rename(TypedTables.extractfield(TypedTables.index($(esc(table))),Val{$(Expr(:quote,expr_right))}), Val{$(Expr(:quote,expr_left))}) = $(esc(table))[Val{$(Expr(:quote,expr_right))}] )
                        exprs[i-1] = :($(esc(expr_left)) = $(esc(table))[Val{$(Expr(:quote,expr_right))}])
                    elseif isa(expr_left, Expr) && expr_left.head == :(::)
                        expr_left.args[2] = :($(esc(expr_left.args[2])))
                        #exprs[i-1] = :($(macroexpand(:(TypedTables.@field($expr_left)))) = $(esc(table))[Val{$(Expr(:quote,expr_right))}] )
                        exprs[i-1] = :($(esc(expr_left)) = $(esc(table))[Val{$(Expr(:quote,expr_right))}])
                    else
                        return(:( error("Expected syntax like @select(table, col1, newname = col2, newcol::newtype = col1 -> f(col1))") ))
                    end
                elseif operation == :compute
                    # LHS must be new field
                    if isa(expr_left,Symbol)
                        lhs = esc(expr_left)
                    elseif isa(expr_left, Expr) && expr_left.head == :(::)
                        #expr_left.args[2] = :($(esc(expr_left.args[2])))
                        #lhs = :($(macroexpand(:(TypedTables.@field($expr_left)))))
                        lhs = esc(expr_left)
                    else
                        return(:( error("Expected syntax like @select(table, col1, newname = col2, newcol::newtype = col1 -> f(col1))") ))
                    end
                    #eltype = :( eltype($lhs) )

                    # RHS is (field_names) -> function. First identify the field names
                    fields = Vector{Symbol}()
                    if isa(expr_right.args[1], Symbol)
                        push!(fields, expr_right.args[1])
                    elseif isa(expr_right.args[1], Expr) && expr_right.args[1].head == :tuple
                        for j = 1:length(expr_right.args[1].args)
                            push!(fields, expr_right.args[1].args[j])
                        end
                    else
                        return(:( error("Expected syntax like @select(table, col1, newname = col2, newcol::newtype = col1 -> f(col1))") ))
                    end

                    # Given field names build the replacements
                    replacements = [:($(esc(table))[idx, Val{$(Expr(:quote,fields[i]))}]) for i = 1:length(fields)]

                    # Now replace the symbols in the function as necessary.
                    func = expr_right.args[2]
                    replace_symbols!(func, fields, replacements)

                    # Now make the expression
                    exprs[i-1] = :( $(lhs) = [$(func) for idx = 1:length($(esc(table)))] )
                end
            #elseif expr.head == :(::)
                # Straightforward extraction (by name and type)
                #exprs[i-1] = :($(macroexpand(:(TypedTables.@field($expr)))) = $(esc(table))[Val{$(Expr(:quote,expr.args[1]))}] )
            else
                return(:( error("Expected syntax like @select(table, col1, newname = col2, newcol::newtype = col1 -> f(col1))") ))
            end
        else
            return(:( error("Expected syntax like @select(table, col1, newname = col2, newcol::newtype = col1 -> f(col1))") ))
        end
    end

    return Expr(:macrocall, Symbol("@Table"), exprs...)
end

macro filter(x...)
    if length(x) == 0 # Is this an error or an empty table?
        return :(error("@filter expects a table"))
    elseif length(x) == 1
        return :($(esc(x[1])))
    end
    table = x[1]
    exprs = Vector{Any}(length(x)-1)
    for i = 1:length(x)-1
        expr = x[i+1]
        if !isa(expr,Expr) || !(expr.head == :(->))
            return :( error("Expected syntax like @filter(table, col1 -> col1 == 1, (col1, col2) -> col1 < col2 )"))
        end

        # Make list of symbols and replacements
        fields = Vector{Symbol}()
        if isa(expr.args[1], Symbol)
            push!(fields, expr.args[1])
        elseif isa(expr.args[1], Expr) && expr.args[1].head == :tuple
            for j = 1:length(expr.args[1].args)
                if !isa(expr.args[1].args[j], Symbol)
                    return :( error("Expected syntax like @filter(table, col1 -> col1 == 1, (col1, col2) -> col1 < col2 )"))
                end
                push!(fields, expr.args[1].args[j])
            end
        else
            return :( error("Expected syntax like @filter(table, col1 -> col1 == 1, (col1, col2) -> col1 < col2 )"))
        end

        # Given field names build the replacements
        replacements = [:($(esc(table))[idx, Val{$(Expr(:quote,fields[i]))}]) for i = 1:length(fields)]

        # Now make replacements in the boolean function
        func = expr.args[2]
        replace_symbols!(func, fields, replacements)

        # Now make the expression
        exprs[i] = func
    end
    test_expr = exprs[1]
    for i = 2:length(exprs)
        test_expr = :($test_expr && $(exprs[i]))
    end

    # Now build a few lines of code
    quote
        test = Bool[$test_expr for idx = 1:length($(esc(table)))]
        $(esc(table))[test]
    end
end

macro filter!(x...)
    if length(x) == 0 # Is this an error or an empty table?
        return :(error("@filter expects a table"))
    elseif length(x) == 1
        return :($(esc(x[1])))
    end
    table = x[1]
    exprs = Vector{Any}(length(x)-1)
    for i = 1:length(x)-1
        expr = x[i+1]
        if !isa(expr,Expr) || !(expr.head == :(->))
            return :( error("Expected syntax like @filter(table, col1 -> col1 == 1, (col1, col2) -> col1 < col2 )"))
        end

        # Make list of symbols and replacements
        fields = Vector{Symbol}()
        if isa(expr.args[1], Symbol)
            push!(fields, expr.args[1])
        elseif isa(expr.args[1], Expr) && expr.args[1].head == :tuple
            for j = 1:length(expr.args[1].args)
                if !isa(expr.args[1].args[j], Symbol)
                    return :( error("Expected syntax like @filter(table, col1 -> col1 == 1, (col1, col2) -> col1 < col2 )"))
                end
                push!(fields, expr.args[1].args[j])
            end
        else
            return :( error("Expected syntax like @filter(table, col1 -> col1 == 1, (col1, col2) -> col1 < col2 )"))
        end

        # Given field names build the replacements
        replacements = [:($(esc(table))[idx, Val{$(Expr(:quote,fields[i]))}]) for i = 1:length(fields)]

        # Now make replacements in the boolean function
        func = expr.args[2]
        replace_symbols!(func, fields, replacements)

        # Now make the expression
        exprs[i] = func
    end
    test_expr = exprs[1]
    for i = 2:length(exprs)
        test_expr = :($test_expr && $(exprs[i]))
    end

    # Now build a few lines of code
    # TODO can switch this to an anonymous function using filter!() for 0.5
    # TODO in general we should also consider a lazy view like filter()
    quote
        test_not = Bool[!$test_expr for idx = 1:length($(esc(table)))]
        deleteat!($(esc(table)), (:)[test_not])
    end
end

macro filter_mask(x...)
    if length(x) == 0 # Is this an error or an empty table?
        return :(error("@filter expects a table"))
    elseif length(x) == 1
        return :($(esc(x[1])))
    end
    table = x[1]
    exprs = Vector{Any}(length(x)-1)
    for i = 1:length(x)-1
        expr = x[i+1]
        if !isa(expr,Expr) || !(expr.head == :(->))
            return :( error("Expected syntax like @filter(table, col1 -> col1 == 1, (col1, col2) -> col1 < col2 )"))
        end

        # Make list of symbols and replacements
        fields = Vector{Symbol}()
        if isa(expr.args[1], Symbol)
            push!(fields, expr.args[1])
        elseif isa(expr.args[1], Expr) && expr.args[1].head == :tuple
            for j = 1:length(expr.args[1].args)
                if !isa(expr.args[1].args[j], Symbol)
                    return :( error("Expected syntax like @filter(table, col1 -> col1 == 1, (col1, col2) -> col1 < col2 )"))
                end
                push!(fields, expr.args[1].args[j])
            end
        else
            return :( error("Expected syntax like @filter(table, col1 -> col1 == 1, (col1, col2) -> col1 < col2 )"))
        end

        # Given field names build the replacements
        replacements = [:($(esc(table))[idx, Val{$(Expr(:quote,fields[i]))}]) for i = 1:length(fields)]

        # Now make replacements in the boolean function
        func = expr.args[2]
        replace_symbols!(func, fields, replacements)

        # Now make the expression
        exprs[i] = func
    end
    test_expr = exprs[1]
    for i = 2:length(exprs)
        test_expr = :($test_expr && $(exprs[i]))
    end

    # Now build a few lines of code
    quote # TODO can switch this to an anonymous function for 0.5
        Bool[$test_expr for idx = 1:length($(esc(table)))]
    end
end

function replace_symbols!(a::Expr, symbols::Vector{Symbol}, exprs::Vector)
    for i = 1:length(a.args)
		if isa(a.args[i], Expr) && a.args[i].head != :line && a.args[i].head != :.
		    replace_symbols!(a.args[i], symbols, exprs)
        elseif isa(a.args[i], Symbol)
            notfound = true
            for j = 1:length(symbols)
                if a.args[i] == symbols[j]
                    a.args[i] = exprs[j]
                    notfound = false
                    break
                end
            end
            # If it's not our symbol then it belongs to the caller's scope
            if notfound
                a.args[i] = :($(esc(a.args[i])))
            end
        elseif isa(a.args[i], Expr) && a.args[i].head == :.
            a.args[i] = :($(esc(a.args[i])))
        end
    end
end


text:
module StatWrapper

using JavaCall

const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const Comparable = @jimport java.lang.Comparable
const Comparator = @jimport java.util.Comparator
const Enum = @jimport java.lang.Enum
const Frequency = @jimport org.hipparchus.stat.Frequency
const JIterator = @jimport java.util.Iterator
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Locale = @jimport java.util.Locale
const LocalizedStatFormats = @jimport org.hipparchus.stat.LocalizedStatFormats
const LongFrequency = @jimport org.hipparchus.stat.LongFrequency
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const StatUtils = @jimport org.hipparchus.stat.StatUtils

export Frequency
export LocalizedStatFormats
export LongFrequency
export StatUtils
export add_value
export clear
export compare_to
export describe_constable
export entry_set_iterator
export equals
export geometric_mean
export get_class
export get_count
export get_cum_freq
export get_cum_pct
export get_declaring_class
export get_localized_string
export get_mode
export get_pct
export get_source_string
export get_sum_freq
export get_unique_count
export hash_code
export increment_value
export max
export mean
export mean_difference
export merge
export min
export mode
export name
export normalize
export notify
export notify_all
export ordinal
export percentile
export population_variance
export product
export sum
export sum_difference
export sum_log
export sum_sq
export to_string
export value_of
export values
export values_iterator
export variance
export variance_difference
export wait

include("CorrelationWrapper/CorrelationWrapper.jl")
include("DescriptiveWrapper/DescriptiveWrapper.jl")
include("FittingWrapper/FittingWrapper.jl")
include("InferenceWrapper/InferenceWrapper.jl")
include("IntervalWrapper/IntervalWrapper.jl")
include("RankingWrapper/RankingWrapper.jl")
include("RegressionWrapper/RegressionWrapper.jl")
include("frequency.jl")
include("localized_stat_formats.jl")
include("long_frequency.jl")
include("stat_utils.jl")

end


text:
using ClangCompiler
using Test

const CC = ClangCompiler

@testset "Diagnostic" begin
    @testset "Engine" begin
        id = CC.DiagnosticIDs()
        @test id.ptr != C_NULL

        opts = CC.DiagnosticOptions()
        @test opts.ptr != C_NULL

        diag = CC.DiagnosticsEngine(id, opts)
        @test diag.ptr != C_NULL
        dispose(diag)
    end

    @testset "Consumer" begin
        opts = CC.DiagnosticOptions()
        @test opts.ptr != C_NULL

        client = CC.IgnoringDiagConsumer()
        @test client.ptr != C_NULL

        diag = CC.DiagnosticsEngine(opts, client)
        @test diag.ptr != C_NULL
        dispose(diag)
    end

    diag = CC.DiagnosticsEngine()
    @test diag.ptr != C_NULL
    dispose(diag)
end


text:
import SMC
using PyPlot
import Distributions
import JSON
using DataFrames

@everywhere begin
    using SMC
    using Distributions
    include("dpmm.jl")
    include("smc_samplers.jl")
    include("../aide.jl")
end

include("estimation.jl")

# where to load and save data
const DATA_DIR = "data"

# load dataset
data = JSON.parsefile("$DATA_DIR/num_clusters_comparison_data.json")
data = convert(Vector{Float64}, data)
num_data = length(data)

# set prior
m_prior = Normal(0, 1) # mean cluster center
r_prior = Gamma(1, 1) # r * data-precision is the precision of cluster centers
nu_prior = Gamma(1, 1) # expected data-precision of cluster is nu/s
s_prior = Gamma(1, 1)
nign_params_prior = NIGNParamsPrior(m_prior, r_prior, nu_prior, s_prior)
alpha_prior = Gamma(1, 1)
hypers = Hyperparameters(nign_params_prior, alpha_prior)

num_sweeps_list = [0, 4]
num_particles_list= [1, 2, 4, 8, 10, 20, 40, 80, 100]

# the number of runs of inference to use for diagnostic is twice the number of
# AIDE replicates, because each AIDE replicate involves running the target
# inference sampler twice
num_aide_replicates = 62
num_diagnostic_replicates = num_aide_replicates * 2
println("generating sampler data...")
gold_standard_sampler = make_optimal_proposal_smc_sampler(data, hypers, 4, 100)
data = generate_sampler_data(data, hypers, num_aide_replicates, num_diagnostic_replicates,
                             gold_standard_sampler)
writetable("$DATA_DIR/sampler_data.csv", data)
println("done!")


text:
# Passing Julia callback functions to Python routines.
#
# Note that this will typically involve two functions: the
# desired Julia function/closure, and a top-level C-callable
# wrapper function used with PyCFunction_NewEx -- the latter
# is called from Python and calls the former as needed.

################################################################

# Define a Python method/function object from f(PyPtr,PyPtr)::PyPtr.
# Requires f to be a top-level function.
function pymethod(f::Function, name::AbstractString, flags::Integer)
    # Python expects the PyMethodDef structure to be a *constant*,
    # so we define an anonymous global to hold it.
    def = gensym("PyMethodDef")
    @eval const $def = PyMethodDef[PyMethodDef($name, $f, $flags)]
    PyObject(@pycheckn ccall((@pysym :PyCFunction_NewEx), PyPtr,
                             (Ptr{PyMethodDef}, Ptr{Void}, Ptr{Void}),
                             eval(def), C_NULL, C_NULL))
end

################################################################

# To pass an arbitrary Julia function to Python, we wrap it
# in a jl_Function Python class, where jl_Function.__call__
# executes the jl_Function_callback function in Julia, which
# in turn fetches the actual Julia function from the "f" attribute
# and calls it.

function jl_Function_call(self_::PyPtr, args_::PyPtr, kw_::PyPtr)
    args = PyObject(args_) # don't need pyincref because of finally clause below
    try
        f = unsafe_pyjlwrap_to_objref(self_)::Function
        if kw_ == C_NULL
            ret = PyObject(f(convert(PyAny, args)...))
        else
            kw = PyDict{Symbol,PyAny}(pyincref(kw_))
            kwargs = [ (k,v) for (k,v) in kw ]
            ret = PyObject(f(convert(PyAny, args)...; kwargs...))
        end
        return pystealref!(ret)
    catch e
        pyraise(e)
    finally
        args.o = PyPtr_NULL # don't decref
    end
    return PyPtr_NULL
end

# just string(Docs.doc(f)) doesn't work on VERSION < v"0.5"
function docstring(f::Function)
    buf = IOBuffer()
    show(buf, "text/plain", Docs.doc(f))
    return String(take!(buf))
end

# this function emulates standard attributes of Python functions,
# where possible.
function jl_Function_getattr(self_::PyPtr, attr__::PyPtr)
    attr_ = PyObject(attr__) # don't need pyincref because of finally clause below
    try
        f = unsafe_pyjlwrap_to_objref(self_)::Function
        attr = convert(String, attr_)
        if attr in ("__name__","func_name")
            return pystealref!(PyObject(string(f)))
        elseif attr in ("__doc__", "func_doc")
            return pystealref!(PyObject(docstring(f)))
        elseif attr in ("__module__","__defaults__","func_defaults","__closure__","func_closure")
            return pystealref!(PyObject(nothing))
        else
            # TODO: handle __code__/func_code (issue #268)

            # fallback: (probably raises AttributeError)
            return ccall(@pysym(:PyObject_GenericGetAttr), PyPtr, (PyPtr,PyPtr), self_, attr__)
        end
    catch e
        pyraise(e)
    finally
        attr_.o = PyPtr_NULL # don't decref
    end
    return PyPtr_NULL
end

function pycallback(f::Function)
    pyjlwrap_new(jl_FunctionType, f)
end

PyObject(f::Function) = pycallback(f)


text:
module OpenGL

import Images.imread

export glimread

# Going with a basic OpenGL RGB 1D-interleaved image-texture
function glimread(filename::String)
    img = imread(filename)
    w = size(img,2)
    h = size(img,3)
    img1D = Array(Uint8,w*h*3)
    for j=1:w
        for k=1:h
            for i=1:3
                img1D[3*w*(j-1)+3*(k-1)+i] = img[i,w-j+1,k]
            end
        end
    end
    return img1D, w, h
end

end


text:
#=
Given a list of integers S and a target number k, write a function that returns a subset of S that adds up to k. If such a subset cannot be made, then return null.

Integers can appear more than once in the list. You may assume all numbers in the list are positive.

For example, given S = [12, 1, 61, 5, 9, 2] and k = 24, return [12, 9, 2, 1] since it sums up to 24.
=#

function subarraywithsum(S::Array{Int, 1}, k::Int)
    sort!(S)
    subarray = copy(S)
    index = 0
    while sum(subarray) != k && index != length(S)
        if sum(subarray) - S[end-index] >= k
            deleteat!(subarray, findlast(isequal(S[end-index]), subarray))
        end
        index = index + 1
    end
    if sum(subarray) == k
        return sort!(subarray, rev = true)
    else
        return nothing
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem