text:
module FiSUDiPPuS

using Optim: optimize, LBFGS, ParticleSwarm, Fminbox
import Optim
using Printf: @printf
using Images: imresize
using JLD2
using FileIO: load, save
using Dates: now
using PyPlot: figure, plot, legend, subplot

export runfit, viewsettings, plotresult, printresult

function runfit(options::String=joinpath(@__DIR__, "../data/default.jl");
        share_a::Bool=true, saveprefix="")
    include(options)
    runfit(settings; share_a=share_a, saveprefix=saveprefix)
end

function runfit(settings::Dict; share_a::Bool=true, saveprefix="")
    ### CONTENTS OF OPTIONS
    ## OPTIMIZATION
    # iterations
    # n_particles
    # lower
    # upper
    # start
    # time_limit
    ## DATA
    # filepath_ssp
    # filepath_ppp
    # roi_ssp
    # roi_ppp
    # size_ssp
    # size_ppp
    # reference_ssp
    # reference_ppp
    # bleach_weight
    # datatype
    # N
    ## SAVING
    # savedir

    ## Convert to correct types
    start = settings[:datatype].(settings[:start])
    lower = settings[:datatype].(settings[:lower])
    upper = settings[:datatype].(settings[:upper])
    bleach_weight = settings[:datatype](settings[:bleach_weight])

    # Process the data
    x, y = get_data(
        settings[:filepath_ssp],
        settings[:filepath_ppp],
        settings[:roi_ssp],
        settings[:roi_ppp],
        settings[:size_ssp],
        settings[:size_ppp],
        settings[:reference_ssp],
        settings[:reference_ppp,],
        bleach_weight,
        settings[:datatype],
    )

    costfun(p) = sum((model(x, p, settings[:N], bleach_weight; share_a=share_a) .- y).^2)

    if settings[:method] == :ParticleSwarm
        method = ParticleSwarm(
            lower=lower,
            upper=upper,
            n_particles=settings[:n_particles],
        )
    elseif settings[:method] == :LBFGS
        method = Fminbox(LBFGS())
    end

    result = optimize(
        costfun, lower, upper, start,
        method,
        Optim.Options(
            iterations=settings[:iterations],
            show_trace=false,
            time_limit=settings[:time_limit],
            )
        )

    save_result(settings[:savedir], x, y, result, settings;
                saveprefix=saveprefix)

    result
end

function viewsettings(file::String)
    include(file)   # contains the settings
    viewsettings(settings)
end

function viewsettings(s::Dict)
    function draw_roi(ax, roi)
        y0 = roi[1][1] - 1 # the plot starts at 0
        x0 = roi[2][1] - 1
        y1 = roi[1][end]
        x1 = roi[2][end]

        rectx = [
            x0 x1 x1 x0
            x1 x1 x0 x0
        ]
        recty = [
            y0 y0 y1 y1
            y0 y1 y1 y0
        ]

        ax.plot(rectx, recty, "r", linewidth=3)
    end

    x, y = get_data(s)

    # load the signal data to view the ROIs
    data_ssp = load(s[:filepath_ssp])
    data_ppp = load(s[:filepath_ppp])
    z_ssp = data_ssp["signal"]
    z_ppp = data_ppp["signal"]

    # Calculate start values
    ystart = startdata(settings)

    f1 = figure()
    plot(y, label="data")
    plot(ystart, label="initial")
    legend()

    f2 = figure()
    ax = subplot(111)
    ax.set_title("SSP Data ROI")
    ax.pcolormesh(z_ssp)
    draw_roi(ax, s[:roi_ssp])

    f3 = figure()
    ax = subplot(111)
    ax.set_title("PPP Data ROI")
    ax.pcolormesh(z_ppp)
    draw_roi(ax, s[:roi_ppp])

    f1, f2, f3
end

function startdata(settings::Dict)
    x, _ = get_data(settings)
    ystart = FiSUDiPPuS.model(
        x, settings[:datatype].(settings[:start]),
        settings[:N],
        settings[:datatype](settings[:bleach_weight])
        )
    ystart
end

function get_data(s::Dict)
    get_data(
        s[:filepath_ssp],
        s[:filepath_ppp],
        s[:roi_ssp],
        s[:roi_ppp],
        s[:size_ssp],
        s[:size_ppp],
        s[:reference_ssp],
        s[:reference_ppp],
        s[:bleach_weight],
        s[:datatype],
    )
end

function get_data(
        filepath_ssp,
        filepath_ppp,
        roi_ssp,
        roi_ppp,
        size_ssp,
        size_ppp,
        reference_ssp,
        reference_ppp,
        bleach_weight,
        datatype,
        )

    data_ssp = load(filepath_ssp)
    data_ppp = load(filepath_ppp)

    # time steps (this ends up being simply a unit range from
    # 1 to the size of the resized data)
    x1_ssp = 1:size_ssp[1]
    x1_ppp = 1:size_ppp[1]
    # wavenumbers
    x2_ssp = data_ssp["wavenumber"][roi_ssp[2]]
    x2_ppp = data_ppp["wavenumber"][roi_ppp[2]];
    # resize wavenumbers
    x2_ssp = imresize(x2_ssp, size_ssp[2])
    x2_ppp = imresize(x2_ppp, size_ppp[2])
    # signal
    y_ssp = data_ssp["signal"][roi_ssp...]
    y_ppp = data_ppp["signal"][roi_ppp...]

    # order the spectra from low to high delay in the roi
    p_ssp = sortperm(data_ssp["dltime"][roi_ssp[1]])
    p_ppp = sortperm(data_ppp["dltime"][roi_ppp[1]])
    # permute
    y_ssp = y_ssp[p_ssp,:]
    y_ppp = y_ppp[p_ppp,:]

    # resize signal
    y_ssp = imresize(y_ssp, size_ssp)
    y_ppp = imresize(y_ppp, size_ppp)

    x = build_independents(x1_ssp, x1_ppp, x2_ssp, x2_ppp)
    y = build_dependents(y_ssp, y_ppp, reference_ssp, reference_ppp, bleach_weight)

    # println("size ssp data: $(size(data_ssp["signal"])) | region of interest: $roi_ssp | resized to $size_ssp")
    # println("size ppp data: $(size(data_ppp["signal"])) | region of interest: $roi_ppp | resized to $size_ppp")

    datatype.(x), datatype.(y)
end

function get_results(filepath)
    data = load(filepath)
    x = data["x"]
    y = data["y"]
    r = data["result"]
    settings = data["settings"]

    x, y, r, settings
end

function plotresult(filepath::String)
    x,y,r,s = get_results(filepath)
    T = s[:datatype]

    ystart = model(
        x, T.(s[:start]),
        s[:N],
        T(s[:bleach_weight])
    )

    ymin = model(
        x, r.minimizer,
        s[:N],
        T(s[:bleach_weight])
    )

    f = figure()
    plot(y, label="data")
    plot(ystart, label="start")
    plot(ymin, label="minimized")
    legend()
    f
end

function printresult(r::String)
    data = load(r)
    result = data["result"]
    printresult(result)
end

function printresult(r)
    # Calculate the number of resonances by checking how many big-number
    # parameters there are.
    N = length(findall(x -> x > 1000, r.minimizer))

    # Get all the parameters from the result
    Assp = r.minimizer[   1: N]
    Appp = r.minimizer[ N+1:2N]
    ω    = r.minimizer[2N+1:3N]
    Γ    = r.minimizer[3N+1:4N]
    a_pow= r.minimizer[end-1]
    Δω   = r.minimizer[end]

    # Get the initial parameter
    Assp0 = r.initial_x[1:N]
    Appp0 = r.initial_x[N+1:2N]
    ω0    = r.initial_x[2N+1:3N]
    Γ0    = r.initial_x[3N+1:4N]
    a_pow0= r.initial_x[end-1]
    Δω0   = r.initial_x[end]

    n = ["", "Assp", "Appp", "ω", "Γ", "A₀ssp", "A₀ppp", "ω₀", "Γ₀"]
    @printf "Δω_ppp: %.3f (start: %.3f)\n" Δω Δω0
    @printf "a_pow:  %.3f (start: %.3f)\n" a_pow a_pow0

    # this prints the header
    for i = 1:length(n)
        @printf "%8.8s" n[i]
    end
    print("\n")

    # this prints the content of the table
    for i = 1:N
        @printf "%8.4u %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n" i Assp[i] Appp[i] ω[i] Γ[i] Assp0[i] Appp0[i] ω0[i] Γ0[i]
    end

    nothing
end

"""
sum frequency spectrum
"""
function sfspec(x, A, ω, Γ)
    y = 0.0 + 0.0im
    for i = 1:size(A,1)
        y += A[i] / (x - ω[i] - 1im * Γ[i]) #* exp(1im * φ[i])
    end
    abs2(y)
end

function model(x::Array{T,2}, p::Array{T,1}, N::Int, bleach_weight::T;
               share_a::Bool=true) where T

    t = @view x[:,1]
    wn = @view x[:,2]
    pol = @view x[:,3]
    dif = @view x[:,4]

    A_ssp = p[   1: N]
    A_ppp = p[ N+1:2N]
    ω     = p[2N+1:3N]
    Γ     = p[3N+1:4N]

    a_pow = p[end-1]    # power by which the a's for the ppp spectra differ from
                        # those of the ssp spectra
    ω_ppp = ω .+ p[end] # add offset parameter to resonance positions for
                        # ppp spectra

    a_idx = 4N+1

    y = Array{T,1}(undef, size(x,1))
    α = Array{T,1}(undef, N)
    a = p[a_idx:a_idx+N-1]
    # define our α
    # TODO: Make sure the ssp spectrum comes first!
    α = a .* A_ssp

    @inbounds for i ∈ eachindex(y)
        if pol[i] == 0
            #ssp spectrum
            # α .= a .* A_ssp

            if dif[i] == 1
                # difference ssp
                y[i] = sfspec(wn[i], α, ω, Γ)
                y[i] -= sfspec(wn[i], A_ssp, ω, Γ)
                y[i] *= bleach_weight
            else
                # plain spectrum ssp
                y[i] = sfspec(wn[i], α, ω, Γ)
            end
        else
            # ppp spectrum
            # if share_a
            #     @. α = (a ^ a_pow) * A_ppp
            # else
            #     @. α = a * A_ppp
            # end

            if dif[i] == 1
                # difference ppp
                y[i] = sfspec(wn[i], α, ω_ppp, Γ)
                y[i] -= sfspec(wn[i], A_ppp, ω_ppp, Γ)
                y[i] *= bleach_weight
            else
                # plain spectrum ppp
                y[i] = sfspec(wn[i], α, ω_ppp, Γ)
            end
        end
        if i != length(y)
            # check if we have to increase the index for
            # the `a`s
            if share_a && t[i+1] < t[i]
                # reset to start value if we change to
                # the next spectrum (no matter if its a plain or
                # a difference spectrum)
                # TODO: make sure the spectra are ordered
                # correctly timewise
                a_idx = 4N+1
                a .= p[a_idx:a_idx+N-1]
                # Since the α values are dependent on the a values we
                # have to change them too. Additionaly we have to check if
                # there's an ssp or a ppp spectrum coming up next.
                if pol[i+1] == 0
                    @. α = a * A_ssp
                else
                    @. α = a ^ a_pow * A_ppp
                end
            elseif t[i+1] != t[i] # just go to the next time step
                a_idx += N
                a .= p[a_idx:a_idx+N-1]
                if pol[i+1] == 0
                    @. α = a * A_ssp
                else
                    @. α = a ^ a_pow * A_ppp
                end
            end
            if !share_a && dif[i+1] > dif[i]
                # if the a values shall not be shared between the
                # spectra reset to the start index when we switch
                # to the difference spectra
                a_idx = 4N+1
                a .= p[a_idx:a_idx+N-1]
                if pol[i+1] == 0
                    @. α = a * A_ssp
                else
                    @. α = a ^ a_pow * A_ppp
                end
            end
        end
    end
    y
end

function convert_independents(x1, x2)
    a1 = eltype(x1)[]
    a2 = eltype(x2)[]
    for _x1 in x1, _x2 in x2
        push!(a1, _x1)
        push!(a2, _x2)
    end
    x = [a1 a2];
end

function build_independents(x1_ssp, x1_ppp, x2_ssp, x2_ppp)
    # we have the following order top to bottom:
    #    all ssp spectra
    #    all ppp spectra
    #    all ssp difference spectra
    #    all ppp difference spectra

    # the convert_independents function gives just a set of
    # time steps and wavenumbers for each polarization
    x_ssp = convert_independents(x1_ssp, x2_ssp)
    x_ppp = convert_independents(x1_ppp, x2_ppp)

    # we can add the polarization information and that we have spectra
    x_ssp = [x_ssp fill(0.0, size(x_ssp,1))]
    x_ppp = [x_ppp fill(1.0, size(x_ppp,1))]

    # to these we can now add if it's a plain spectrum or if we have
    # a difference spectrum
    x_ssp_spec = [x_ssp fill(0.0, size(x_ssp,1))]
    x_ppp_spec = [x_ppp fill(0.0, size(x_ppp,1))]
    x_ssp_diff = [x_ssp fill(1.0, size(x_ssp,1))]
    x_ppp_diff = [x_ppp fill(1.0, size(x_ppp,1))]

    # now we just have to stich everything together in the correct order
    # (vertical concatenation)
    [x_ssp_spec; x_ppp_spec; x_ssp_diff; x_ppp_diff]
end

function build_dependents(y_ssp, y_ppp, r_ssp::Int, r_ppp::Int, bleach_weight)
    # first the plain spectra (turn the matrix and reshape to 1d array)
    # then concatenate vertically
    y_spec = [y_ssp'[:]; y_ppp'[:]]
    # get the reference spectra
    y_ref_ssp = y_ssp[r_ssp,:]
    y_ref_ppp = y_ppp[r_ppp,:]
    # calculate the difference spectra and concatenate
    y_diff_ssp = vcat([y_ssp[i,:] .- y_ref_ssp for i=1:size(y_ssp,1)]...) .* bleach_weight
    y_diff_ppp = vcat([y_ppp[i,:] .- y_ref_ppp for i=1:size(y_ppp,1)]...) .* bleach_weight
    y_diff = [y_diff_ssp; y_diff_ppp]
    [y_spec; y_diff]
end

function save_result(savedir,x,y,r,settings; saveprefix="")
    method = typeof(r.method)
    !isdir(savedir) && mkdir(savedir)
    datestr = string(now())
    datestr = replace(datestr, ":" => ".")
    filename = saveprefix * "_" * datestr * ".jld2"
    savepath = joinpath(savedir, filename)

    d = Dict(
        "result" => r,
        "x" => x,
        "y" => y,
        "settings" => settings,
    )

    save(savepath, d)
    nothing
end

end # module


text:
abstract type AbstractMatrix end

struct SoA <: AbstractMatrix end
struct AoS <: AbstractMatrix end

"""
StateMatrix{F <: AbstractMatrix}(values::Matrix{T}) where {T <: AbstractFloat}
Holds a formated coordinate system. The format is given by the F type,
which can be SoA or AoS. The usage of this typed matrix allows the
correct retrieval of values in an uniform fashion, despite the format.
Example:
_coords = convert(Matrix{Float64}, [1 2 3; 4 5 6])
coords = LennardJones.StateMatrix{LennardJones.AoS, Float64}(_coords)
coords[1, 2] = 2.0
where 1 is the atom position and 2 is the desired coordinate;
"""
struct StateMatrix{F <: AbstractMatrix, T <: AbstractFloat}
    values::Matrix{T}
end

Base.getindex(M::StateMatrix{SoA, T}, i::Int, j::Int) where {T} = M.values[i, j]
Base.getindex(M::StateMatrix{AoS, T}, i::Int, j::Int) where {T} = M.values[j, i]
Base.getindex(M::StateMatrix{SoA, T}, i::Int) where {T} = M.values[i, :]
Base.getindex(M::StateMatrix{AoS, T}, i::Int) where {T} = M.values[:, i]
Base.setindex!(M::StateMatrix{SoA, T}, v, i::Int, j::Int) where {T <: AbstractFloat} = setindex!(M.values, v, i, j)
Base.setindex!(M::StateMatrix{AoS, T}, v, i::Int, j::Int) where {T <: AbstractFloat} = setindex!(M.values, v, j, i)
Base.length(M::StateMatrix) = max(size(M.values)...)

"""
    State(n::Int, coords::Matrix{T}, forces::Matrix{T}) where {T <: AbstractFloat}
    State(n::Int, coords::Matrix{T}) where {T <: AbstractFloat}
        Holds informating regarding the current coordinates and forces on a
        particle system with a given 'n' size. Input can be either in SoA or AoS
        format, in either case coordinates and forces can be retrieved and set
        using state.coords[i, j] (where i is the atom number and j is the
        coordinate position). If the 'forces' matrix is not given, a new Matrix
        is created with size equal to coords set to all zeros.
    Example:
        state = State(2, [1.0 2.0 3.0; 4.0 5.0 6.0], zeros(3, 2))
        state = State(2, [1.0 2.0 3.0; 4.0 5.0 6.0])
"""
mutable struct State{F <: AbstractMatrix, T <: AbstractFloat}
    n::Int
    coords::StateMatrix{F, T}
    forces::StateMatrix{F, T}

    function State(coords::Matrix{T}, forces::Matrix{T}) where {T <: AbstractFloat}
        if size(coords) != size(forces)
            error("Size of coords does not match size of forces")
        end
        if argmax(collect(size(coords))) == 2
            println("Creating state in AoS format")
            new{AoS, T}(max(size(coords)...), StateMatrix{AoS, T}(coords), StateMatrix{AoS, T}(forces))
        elseif argmax(collect(size(coords))) == 1
            println("Creating state in SoA format")
            new{SoA, T}(max(size(coords)...), StateMatrix{SoA, T}(coords), StateMatrix{SoA, T}(forces))
        else
            error("Could not infer the AoS or SoA format of input data.")
        end
    end
end

State(coords::Matrix{T}) where {T <: AbstractFloat} = State(coords, zeros(T, size(coords)))

ext: jl


text:
using CASInterface
using CCAS

const LIBCAS = @windows ? Pkg.dir("CCAS/libcas0.10.1/lib/libcas.dll") : @linux ? Pkg.dir("CCAS/libcas0.10.1/lib/libcas.so") : error("OS not supported")
const LIBCAS_CONFIG = Pkg.dir("CCAS/libcas0.10.1/parameters/0.10.1.standard.r15mtf.tcas.xa.config.txt")
const LIBCAS_R14 = @windows ? Pkg.dir("CCAS/libcas0.9.3/lib/libcas.dll") : @linux ? Pkg.dir("CCAS/libcas0.9.3/lib/libcas.so") : error("OS not supported")
const LIBCAS_CONFIG_R14 = Pkg.dir("CCAS/libcas0.9.3/parameters/0.9.3.standard.r14.xa.config.txt")

function set_input_coc!(input::Input)
    input.ownInput.dz = 0.0
    input.ownInput.z = 1665
    input.ownInput.psi = 0.0
    input.ownInput.h = 1665
    input.ownInput.modes = 0x1
    input.intruders[1].valid = true
    input.intruders[1].id = 100
    input.intruders[1].modes = 0x2
    input.intruders[1].chi = -1.2
    input.intruders[1].sr = 16500
    input.intruders[1].z = 2200
    input.intruders[1].cvc = 0x0
    input.intruders[1].vrc = 0x0
    input.intruders[1].vsb = 0x0
    input.intruders[1].equipage = EQUIPAGE_MODES
    input.intruders[1].quant = 25
    input.intruders[1].sensitivity_index = 0x0
    input.intruders[1].protection_mode = 0x0
    input
end

function set_input_ra!(input::Input, t)
    input.ownInput.dz = 0.0
    input.ownInput.z = 5000
    input.ownInput.psi = 0.0
    input.ownInput.h = 5000
    input.ownInput.modes = 0x1
    input.intruders[1].valid = true
    input.intruders[1].id = 100
    input.intruders[1].modes = 0x2
    input.intruders[1].chi = 0.0
    input.intruders[1].sr = 10000 - 500 * t
    input.intruders[1].z = 5000
    input.intruders[1].cvc = 0x0
    input.intruders[1].vrc = 0x0
    input.intruders[1].vsb = 0x0
    input.intruders[1].equipage = EQUIPAGE_TCAS
    input.intruders[1].quant = 25
    input.intruders[1].sensitivity_index = 0x0
    input.intruders[1].protection_mode = 0x0
    input
end

function timingtest(libcas::AbstractString=LIBCAS,
    libcas_config::AbstractString=LIBCAS_CONFIG; ntrials::Int64=10, nsteps::Int64=50,
    printout::Bool=true)

    println("===============")
    println("Timing test")
    println("===============")

    result = Dict{ASCIIString,Any}()
    result["reset"] = Array(Float64, ntrials)
    result["update"] = Array(Float64, (ntrials, nsteps))

    println(author())
    consts = Constants(25, libcas_config, 1)
    cas = CASShared(libcas, consts)
    println(version(cas))

    nintruders = max_intruders(cas)
    input = Input(nintruders)
    output = Output(nintruders)

    for i = 1:ntrials
        println("trial $i of $ntrials")
        tic();
        reset(cas)
        result["reset"][i] = toq();

        for t = 1:nsteps
            set_input_ra!(input, t)
            tic();
            update!(cas, input, output)
            result["update"][i,t] = toq();
        end
    end
    if printout
        println("Results")
        println("=======")
        for (item, times) in result
            println("$item: $(mean(times)) seconds")
        end
    end
    result
end

function estimate_cas_time(timingresult::Dict{ASCIIString,Any}, n_steps, n_iters, n_aircraft)
    reset_s = mean(timingresult["reset"])
    update_s = mean(timingresult["update"])
    cas_time_s = (reset_s + n_steps*update_s) * n_iters * n_steps * n_aircraft
    println("estimated cas time: $(cas_time_s) seconds, or $(cas_time_s/3600.0) hours")
    cas_time_s
end

function runtest(libcas::AbstractString=LIBCAS,
    libcas_config::AbstractString=LIBCAS_CONFIG)

    println(author())

    consts = Constants(25, libcas_config, 1)
    cas = CASShared(libcas, consts)
    println("cas handle = ", cas.handle)

    println(version(cas))
    println("max_intruders: ", max_intruders(cas))

    nintruders = max_intruders(cas)
    input = Input(nintruders)
    output = Output(nintruders)
    println("output handle = ", cas.output.handle)

    println("===============")
    println("COC test")
    println("===============")
    for i = 1:5
        println("\ni = ", i)
        reset(cas)

        for t = 1:1
            println("\nt = ", t)

            set_input_coc!(input)
            update!(cas, input, output)

            #xdump(input.ownInput)
            #xdump(input.intruders[1])
            #xdump(output)
            #xdump(output.intruders[1])
            println("[",output.cc,",",output.vc,",",output.ua,",",output.da,"]")
            println("target_rate=",output.target_rate)
            println("dh_min=",output.dh_min,", dh_max=",output.dh_max)

            errorMsg = error_msg(cas)
            println( errorMsg == nothing ? "No Errors" : errorMsg )
        end
    end

    println("")
    println("===============")
    println("active RA test")
    println("===============")
    for t = 1:20
        println("\nt = ", t)

        set_input_ra!(input, t)
        update!(cas, input, output)

        #xdump(input.ownInput)
        #xdump(input.intruders[1])
        #xdump(output)
        #xdump(output.intruders[1])
        println("own_out[",output.cc,",",output.vc,",",output.ua,",",output.da,"]")
        println("intr_out[",output.intruders[1].cvc,",",output.intruders[1].vrc,",",output.intruders[1].vsb,"]")
        println("target_rate=",output.target_rate)
        println("dh_min=",output.dh_min,", dh_max=",output.dh_max)

        errorMsg = error_msg(cas)
        println( errorMsg == nothing ? "No Errors" : errorMsg )
    end

    println()
    println("Done!")
end

max_stars_repo_name: sisl/CCAS

max_issues_repo_name: sisl/CCAS

max_forks_repo_name: sisl/CCAS


text:
include("./Algorithms/TEAVAR.jl")
include("./simulation.jl")
include("./util.jl")
include("./parsers.jl")

function findBeta(alpha, sigfigs, links, capacity, flows, demand, cutoff, T, Tf, k, scenarios, scenario_probs; cvar=false, allocations=nothing)

    step_size = 10.0^(-1 * sigfigs)
    low = 0
    high = 1 - step_size
    index = cvar ? 2 : 1

    while low <= high
        middle = round((low + high)/2, digits=sigfigs)
        val = allocations != nothing ?
            TEAVAR(links, capacity, flows, demand, middle, k, T, Tf, allocations=a)[index] :
            TEAVAR(links, capacity, flows, demand, middle, k, T, Tf, scenarios, scenario_probs)[index]
        if val <= alpha
            low = middle + step_size
            if (low == 1) break end
            val = allocations != nothing ?
                TEAVAR(links, capacity, flows, demand, low, k, T, Tf, allocations=a)[index] :
                TEAVAR(links, capacity, flows, demand, low, k, T, Tf, scenarios, scenario_probs)[index]
            if val > alpha
                return middle
            end
        else
            high = middle - step_size
            val = allocations != nothing ?
                TEAVAR(links, capacity, flows, demand, high, k, T, Tf, allocations=a)[index] :
                TEAVAR(links, capacity, flows, demand, high, k, T, Tf, scenarios, scenario_probs)[index]
            if val <= alpha
                return high
            end
        end
    end
    return high < .5 ? 0 : 1
end


text:
struct AHPQdata
    clusterdata
    qd::QuantizerData
    residuals::AbstractMatrix
    norms::AbstractArray
    config::NamedTuple
end

"""
    function builder(data::Matrix; T::Real, kwargs...)

Main constructor for building the (Anisotropic) Hierarchical Product Quantizer.\n
Leave out the keyword parameter (T::AbstractFloat) to train on Euclidean loss instead of Anisotropic loss

# Main Arguments
 * `Matrix` input data
 * `T::AbstractFloat` Threshold parameter of the Anisotropic loss function, set to 0 to use L2_loss

# Kwargs...
 * `T_preclustering::Int=-1` the `T` parameter for the preclustering loss function, set to -1 to use the same value as for `T`, 0 to use L2_loss
 * `n_codebooks::Int=n_dims/2` number of codebooks to train
 * `n_centers::Int=8` number of clusters/centers to train
 * `a::Int=√n_dp` number of clusters generated in preclustering step, set to 0 to turn off
 * `b::Int=(n_clusters_to_generate ÷ 5 + 2)` number of clusters pruned from preclustering step
 * `initialise_with_euclidean_loss::Bool=true` switches euclidean pre-training on/off
 * `inverted_index::Bool=true` switches to IVF methods for Euclidean quantization
 * `max_iter::Int=1000` the max number of iterations of the assignment_step-codebook_update loop
 * `stopcond<:AbstractFloat=1e-2` the stopcondition (l2-distance of codebook update) for the assignment_step-codebook_update loop
 * `verbose::Bool=false` switches training information updates on/off
 * `max_iter_assignments::Int=10` maximum iterations for the approximate assignment step of anisotorpic quantization
 * `optimisation::String="exact"` the optimisation method for the codebook udpate step, takes: ("exact", "nesterov")
 * `multithreading::Bool=false` switches multi-threading on/off
 * `GPU::Bool=false` switches optimization on GPU on/off (uses CUDA.jl)
 * `increment_steps::Int=floor(log10(training_points))-2)` number of incremental steps, use 0 to switch off incremental training
 * `reorder::Int=0` number of datapoints used for exact inner product reordering
 * `training_points::Int=ceil(n_dp/5))` number of training points to train on

 # Example
 ```julia
 traindata = rand(d, n)
 l2_searcher          = builder(traindata; n_codebooks=50,
                                           n_centers=16)
 anisotropic_searcher = builder(traindata, 0.2; n_codebooks=50,
                                                n_centers=16)
 ```
"""
function builder(data::Matrix; T::Real, kwargs...)
    # Step 0 - Retrieve configurations
    n_dims, n_dp = size(data)
    config, T = check_kwargs(kwargs, n_dp, n_dims, T)
    η = if T > 0 ComputeWeightsFromT(n_dims, T) else L2_loss() end
    η_cd = if (config.T_preclustering > 0) ComputeWeightsFromT(n_dims, config.T_preclustering) else L2_loss() end

    # Step 1 - Preclustering
    if config[:a] > 0
        if config[:verbose] @info("Generating precluster data...") end
        clusterdata =  GenerateClusterData(data, η_cd, config)

        if config[:verbose] @info("Computing residuals...") end
        for (i, center) in enumerate(eachcol(clusterdata.centers))
             data[:,clusterdata.assignments[i]] .-= center
        end
        norms = mapslices(norm, data, dims=1)
    else
        clusterdata = nothing
        norms = ones(n_dp)'
    end

    # Step 2 - Quantization
    codebook=0
    traindata = if config.training_points > 0 subsample(config.training_points, data./norms) else data./norms end
    if config[:increment_steps] > 0
        if config[:verbose] @info("Starting incremental initialisation...") end
        qd = incremental_quantization(traindata, η, config)
    else
        if (config[:initialise_with_l2_loss]) & (T > 0)
            if config[:verbose] @info("Starting L2 initialisation...") end
            qd = quantizer(traindata, L2_loss(), codebook, merge(config,(;verbose=false)))
            codebook = deepcopy(qd.C)
        end
        if config[:verbose] @info("Finished initialisation...") end
        qd = quantizer(traindata, η, codebook, config)
    end
    codebook = deepcopy(qd.C)

    ## Final Assignments on complete data set
    qd = GenerateQuantizerData(data, config.n_codebooks, config.n_centers)
    thread = if config.multithreading MultiThreaded() else SingleThreaded() end
    initialization!(data, qd, thread, codebook)
    assignment_step!(data, qd, BMatrix(), thread)
    if T>0 assignment_step!(data, qd, η, config[:max_iter_assignments], thread) end

    # Step 3 - Generate residuals for exact reordering
    for i in 1:n_dp
        data[:,i] .-= norms[i].*(qd.I.B[i]'qd.C)
    end

    return AHPQdata(clusterdata, qd, data, norms, config)
end


text:
# using Dolo
# import Dolo: Cartesian

module Temp

    using IterTools
    using StaticArrays
    const Point{d} = SVector{d, Float64}

    abstract type AbstractGrid{d} end

    # Cartesian Grid

    struct CGrid{d} <: AbstractGrid{d}
        min::Point{d}
        max::Point{d}
        n::SVector{d, Int64}
    end

    n_nodes(cg::CGrid{d}) where d = prod(cg.n)

    function iter(cg::CGrid{d}) where d
        # ranges = (range(e[1], e[2]; length=e[3]) for e in zip(cg.min, cg.max, cg.n))
        ranges = (range(cg.min[i], cg.max[i];length=cg.n[i]) for i=1:d)
        return ( SVector{d}(v...) for v in Iterators.product(ranges...) )
    end

    CGrid{d}(min::Vector, max::Vector, n::Vector) where d = CGrid(Point{d}(min...), Point{d}(max...), SVector(n...))

    struct PGrid{d, G1, G2} <: AbstractGrid{d}
        g1:: G1
        g2:: G2
    end

    ndims(pg::PGrid{d}) where d = d

    # that one looks efficient
    function iter(pg::PGrid{d, G1, G2}) where d where G1 where G2
        iters = (iter(pg.g1), iter(pg.g2))
        return ( ((k[1][1], k[2][1]), SVector(k[1][2]..., k[2][2]...)) for (k) in Iterators.product(iters...))
    end

    PGrid(g1::AbstractGrid{d1}, g2::AbstractGrid{d2}) where d1 where d2 = PGrid{d1+d2, typeof(g2), typeof(g2)}(g1, g2)
    n_nodes(pg::PGrid) = n_nodes(pg.g1)*n_nodes(pg.g2)

    struct PPGrid{T<:Tuple}
        grids::T
    end

    PPGrid(grids...) = PPGrid(grids)

    n_nodes(pg::PPGrid) = prod( (n_nodes(gr) for gr in pg.grids) )

    concat(x1::SVector{d1,Float64}, x2::SVector{d2,Float64}) where d1 where d2 = SVector{d1+d2, Float64}(x1...,x2...)
    concat(x1::SVector{d1,Float64}, x2::SVector{d2,Float64}, x3::SVector{d3,Float64}) where d1 where d2 where d3 = SVector{d1+d2+d3, Float64}(x1...,x2..., x3...)

    function iter(ppg::PPGrid)
        nodes = tuple( ([e for e in iter(gr)] for gr in ppg.grids)... )
        return ( concat(v...) for v in  Iterators.product(nodes...) )
    end

    import Base
    function Base.enumerate(ppg::PPGrid)
         return (
            ( tuple( (e[1] for e in v)...) , concat( (e[2] for e in v)...))
            for v in Iterators.product( (enumerate(iter(gr)) for gr in ppg.grids)... )
        )
    end

    # Base.getindex(ppg, inds...) = concat()

    # uses less memory but slower
    function iter2(ppg::PPGrid)
        return ( concat(v...) for v in Iterators.product( (iter(gr) for gr in ppg.grids) ... ) )
    end

    function test_iter_pg(pg::PGrid)
        k = 0.0
        for i in iter(pg)
            v = i[2]
            k += sum(v)
        end
        return k
    end

    function test_iter(cg::CGrid{d}) where d
        l = 0.0
        for v in iter(cg)
            l += sum(v)
        end
        return l
    end

    function test_iter_ppg(ppg::PPGrid)
        l = 0.0
        for v in iter(ppg)
            l += sum(v)
        end
        return l
    end

    function test_iter2_ppg(ppg::PPGrid)
        l = 0.0
        for v in iter2(ppg)
            l += sum(v)
        end
        return l
    end

    # PGrid(g1::AbstractGrid{d1}, g2::AbstractGrid{d2}) where d1 where d2 = PGrid{d1+d2, typeof(g2), typeof(g2)}(g1, g2)

end

import Main.Temp as M

g = M.CGrid{2}([0.0, 0.0], [1.0, 1.0], [5,5])

pg = M.PGrid(g,g)

pggg = M.PPGrid(g,g)

@time M.test_iter(g)
@time M.test_iter_pg(pg)
@time M.test_iter_ppg(pggg)
@time M.test_iter2_ppg(pggg)

@time M.test_iter(g)
@time M.test_iter_pg(pg)
@time M.test_iter_ppg(pggg)
@time M.test_iter2_ppg(pggg)

@code_warntype M.test_iter(g)

@code_warntype M.test_iter_ppg(pggg)
@code_warntype M.test_iter2_ppg(pggg)

function test()
    t = 0.0
    for k in M.enumerate(pggg)
        t += sum(k[2])

    end
    return t
end


text:
cd(dirname(@__FILE__))
include("./functions/Load_all_stuff.jl")

# */------------------------------------------------------------------------/* #
RMP = create_RMPo(ms,mp,ps,pp,hc,al)
solve_RMP!(RMP)
last_step_RMP!(RMP,hc,ps)
write_output(RMP,ms,ps,sd)
println("")
println(" finished - optimal solution found")
# */------------------------------------------------------------------------/* #

max_stars_repo_name: nimazzi/PlanModelDet

max_issues_repo_name: nimazzi/PlanModelDet

max_forks_repo_name: nimazzi/PlanModelDet


text:
using Plots, StatsBase, Statistics, LaTeXStrings, JLD, ProgressMeter, LsqFit

function InitialNetwork(dim, P)
    Network = sample([-1,0], Weights([1-P, P]),(dim,dim))
    return Network
end

function FindNeighbors(i, j, Network, dim)
    Numbereds = 0
    colorList = []
    for neighbor in ([0,-1], [-1,0])
        if j + neighbor[2] >= 1 && i + neighbor[1] >= 1 && Network[([i, j] + neighbor)...] > 0
            push!(colorList,Network[([i, j] + neighbor)...])
            Numbereds += 1
        end
    end
    return colorList, Numbereds #colorList contain the node and its neighbors
end

function RecursionMapping(x, L)
    while L[x] != L[L[x]]
        x = L[x]
    end
    return L[x]
end

@inline function HKNetworkDynamic(dim, P)
    Network = InitialNetwork(dim, P)
    S = []
    L = []
    k = 1

    for j in 1:dim
        for i in 1:dim
            if Network[i,j] != 0
                continue
            end
            NColor, NNum = FindNeighbors(i,j,Network,dim)
            if NNum == 0
                Network[i,j] = k
                push!(L, k)
                push!(S, 1)
                k += 1
            elseif NNum == 1
                Network[i,j] = RecursionMapping(NColor[1],L)
                S[RecursionMapping(NColor[1],L)] += 1
            elseif NNum == 2 && RecursionMapping(NColor[1],L) == RecursionMapping(NColor[2],L)
                Network[i,j] = RecursionMapping(NColor[1],L)
                S[RecursionMapping(NColor[1],L)] += 1
            else
                S[RecursionMapping(NColor[1],L)] += 1 + S[RecursionMapping(NColor[2],L)]
                S[RecursionMapping(NColor[2],L)] = 0
                Network[i,j] = RecursionMapping(NColor[1],L)
                L[NColor[2]] = RecursionMapping(NColor[1],L)
            end
        end
    end
    return Network, S, L
end

function JoinColors(Network, L)
    dim = size(Network)[1]
    for i in 1:dim
        for j in 1:dim
            if Network[i,j]!= -1
                Network[i,j] = RecursionMapping(Network[i,j],L)
            end
        end
    end
    return Network
end

function FindCorrelationLength(Network, S)
    if length(Set(S)) <= 2
        return 0.0
    end
    if findfirst(x->x==max(S...),S) ∈ intersect(Network[1,:],Network[end,:])
        BiggestFinite = findfirst(x->x==sort(S, rev=true)[2],S)
    else
        BiggestFinite = findfirst(x->x==max(S...),S)
    end
    BFiniteCluster = Tuple.(findall(x->x==BiggestFinite,Network))
    BFiniteClusterSize = length(BFiniteCluster)
    MassCenter = reduce(.+, BFiniteCluster) ./ BFiniteClusterSize
    return √(sum(pos -> reduce(.+, (pos .- MassCenter).^2),Tuple.(BFiniteCluster))/BFiniteClusterSize)
end

function ReturnXi(dim, p)
    Network, S, L = HKNetworkDynamic(dim, p)
    Network = JoinColors(Network, L)
    return  FindCorrelationLength(Network, S)
end

dimlist = [10,15,20,25,30,35,40,45,50,55,60,65,70,75,100,125,150]
# runlist = [10000,5000,2000,1000,500,200,100,100,100]
PList = hcat(0.52:0.001:0.62)
runnum = 100
progress = Progress(length(PList)*runnum*length(dimlist); showspeed=true)
# progress = Progress(length(PList)*sum(runlist); showspeed=true)
for n in 1:length(dimlist)
    dim = dimlist[n]
    AvgData = []
    STDData = []
    TopAvg = []
    for p in PList
        totalruns = []
        for i in 1:runnum
            push!(totalruns, ReturnXi(dim, p))
            next!(progress)
            update!(progress)
        end
        push!(TopAvg,mean(sort(totalruns,rev = true)[1:Int(runnum/5)]))
        push!(AvgData,mean(totalruns))
        push!(STDData,std(totalruns))
    end
    save("../../Data/Q2/Q2-$dim-fc.jld", "TopAvg", TopAvg, "AvgData", AvgData, "STDData", STDData)
end

Data = zeros(length(dimlist),length(PList),3)

for n in 1:length(dimlist)
    print("\r$n")
    dim = dimlist[n]
    data = load("../../Data/Q2/Q2-$dim-fc.jld")
    Data[n,:,1] = data["TopAvg"]
    Data[n,:,2] = data["AvgData"]
    Data[n,:,3] = data["STDData"]
end

PmaxList = []
for n in 1:length(dimlist)
    push!(PmaxList, PList[findfirst(x->x==max(Data[n,:,2]...),Data[n,:,2])])
end

@. modelP∞(x, p) = abs(x - p[1])^(-1.2)
@. modelν(x, p) = abs(x - 0.5927)^(-p[1])

xdata = PmaxList
ydata = dimlist
InitP∞ = [0.59]
Initν = [1.2]

ν = curve_fit(modelν, xdata, ydata, Initν).param[1]
P∞ = curve_fit(modelP∞, xdata, ydata, InitP∞).param[1]

P∞range = hcat(0.52:0.001:0.596)
νrange = abs.(P∞range .- P∞).^(-ν)

scatter(PmaxList,dimlist, label=L"Data\ Points", title=L"Curve\ fitting\ Plot\ (\nu = %$(round(ν,digits= 2)) ,\ Pc_{\infty} = %$(round(P∞,digits= 2)))",legend = 150)
plot!(P∞range,νrange,c= :purple,lw = 2,label = L"L=|Pc_{(L)} - Pc_{\infty}|^{-\nu}")
plot!([0.6,0.6],[0.0,200.0], linestyle = :dash,c = :black,lw = 2,
    label = L"X=Pc_{\infty}= %$(round(P∞,digits= 2))")
savefig("../../Figs/Q2/Q2-CF-res.pdf")


text:
# Main module, used to generate the lexer parser file (parser.yy.jl usually)

include("lexer.jl")
include("slr.jl")
include("lexerparser.inc.jl")
include("lexer_template.jl")

function generate(io, regs, tokens, prods, produce_rules)
    parserdir = dirname(@__FILE__)
    # Header
    src_header = "# This file is generated by the lexer-parser, do not modify"
    println(io, src_header)

    # Includes
    src_include = read("$parserdir/lexerparser.inc.jl", String)
    println(io, src_include)

    # Lexer
    generate_lexer(io, regs)

    # Parser
    generate_parser(io, tokens, prods, produce_rules)

    # Main function
    src_parse = read("$parserdir/lexerparser_template.jl", String)
    print(io, src_parse)
end


text:
# See imp.jl for the implementations of these functions
function listFunctions() end
function listMethods() end
function createDefaultArguments() end
function createDefaultObject() end
function invokeMethod() end
function invokeMethodOnAllProcs() end
function invokeMethodsOfModule() end
function getNonMissingTypeOfUnionType() end


text:
module TestUtils

using AbstractMCMC
using DynamicPPL
using LinearAlgebra
using Distributions
using Test

# A collection of models for which the mean-of-means for the posterior should
# be same.
@model function demo_dot_assume_dot_observe(
    x=[10.0, 10.0], ::Type{TV}=Vector{Float64}
) where {TV}
    # `dot_assume` and `observe`
    m = TV(undef, length(x))
    m .~ Normal()
    x ~ MvNormal(m, 0.25 * I)
    return (; m=m, x=x, logp=getlogp(__varinfo__))
end

@model function demo_assume_index_observe(
    x=[10.0, 10.0], ::Type{TV}=Vector{Float64}
) where {TV}
    # `assume` with indexing and `observe`
    m = TV(undef, length(x))
    for i in eachindex(m)
        m[i] ~ Normal()
    end
    x ~ MvNormal(m, 0.25 * I)

    return (; m=m, x=x, logp=getlogp(__varinfo__))
end

@model function demo_assume_multivariate_observe_index(x=[10.0, 10.0])
    # Multivariate `assume` and `observe`
    m ~ MvNormal(zero(x), I)
    x ~ MvNormal(m, 0.25 * I)

    return (; m=m, x=x, logp=getlogp(__varinfo__))
end

@model function demo_dot_assume_observe_index(
    x=[10.0, 10.0], ::Type{TV}=Vector{Float64}
) where {TV}
    # `dot_assume` and `observe` with indexing
    m = TV(undef, length(x))
    m .~ Normal()
    for i in eachindex(x)
        x[i] ~ Normal(m[i], 0.5)
    end

    return (; m=m, x=x, logp=getlogp(__varinfo__))
end

# Using vector of `length` 1 here so the posterior of `m` is the same
# as the others.
@model function demo_assume_dot_observe(x=[10.0])
    # `assume` and `dot_observe`
    m ~ Normal()
    x .~ Normal(m, 0.5)

    return (; m=m, x=x, logp=getlogp(__varinfo__))
end

@model function demo_assume_observe_literal()
    # `assume` and literal `observe`
    m ~ MvNormal(zeros(2), I)
    [10.0, 10.0] ~ MvNormal(m, 0.25 * I)

    return (; m=m, x=[10.0, 10.0], logp=getlogp(__varinfo__))
end

@model function demo_dot_assume_observe_index_literal(::Type{TV}=Vector{Float64}) where {TV}
    # `dot_assume` and literal `observe` with indexing
    m = TV(undef, 2)
    m .~ Normal()
    for i in eachindex(m)
        10.0 ~ Normal(m[i], 0.5)
    end

    return (; m=m, x=fill(10.0, length(m)), logp=getlogp(__varinfo__))
end

@model function demo_assume_literal_dot_observe()
    # `assume` and literal `dot_observe`
    m ~ Normal()
    [10.0] .~ Normal(m, 0.5)

    return (; m=m, x=[10.0], logp=getlogp(__varinfo__))
end

@model function _prior_dot_assume(::Type{TV}=Vector{Float64}) where {TV}
    m = TV(undef, 2)
    m .~ Normal()

    return m
end

@model function demo_assume_submodel_observe_index_literal()
    # Submodel prior
    m = @submodel _prior_dot_assume()
    for i in eachindex(m)
        10.0 ~ Normal(m[i], 0.5)
    end

    return (; m=m, x=[10.0], logp=getlogp(__varinfo__))
end

@model function _likelihood_dot_observe(m, x)
    return x ~ MvNormal(m, 0.25 * I)
end

@model function demo_dot_assume_observe_submodel(
    x=[10.0, 10.0], ::Type{TV}=Vector{Float64}
) where {TV}
    m = TV(undef, length(x))
    m .~ Normal()

    # Submodel likelihood
    @submodel _likelihood_dot_observe(m, x)

    return (; m=m, x=x, logp=getlogp(__varinfo__))
end

@model function demo_dot_assume_dot_observe_matrix(
    x=fill(10.0, 2, 1), ::Type{TV}=Vector{Float64}
) where {TV}
    m = TV(undef, length(x))
    m .~ Normal()

    # Dotted observe for `Matrix`.
    x .~ MvNormal(m, 0.25 * I)

    return (; m=m, x=x, logp=getlogp(__varinfo__))
end

const DEMO_MODELS = (
    demo_dot_assume_dot_observe(),
    demo_assume_index_observe(),
    demo_assume_multivariate_observe_index(),
    demo_dot_assume_observe_index(),
    demo_assume_dot_observe(),
    demo_assume_observe_literal(),
    demo_dot_assume_observe_index_literal(),
    demo_assume_literal_dot_observe(),
    demo_assume_submodel_observe_index_literal(),
    demo_dot_assume_observe_submodel(),
    demo_dot_assume_dot_observe_matrix(),
)

# TODO: Is this really the best/most convenient "default" test method?
"""
    test_sampler_demo_models(meanfunction, sampler, args...; kwargs...)

Test that `sampler` produces the correct marginal posterior means on all models in `demo_models`.

In short, this method iterators through `demo_models`, calls `AbstractMCMC.sample` on the
`model` and `sampler` to produce a `chain`, and then checks `meanfunction(chain)` against `target`
provided in `kwargs...`.

# Arguments
- `meanfunction`: A callable which computes the mean of the marginal means from the
  chain resulting from the `sample` call.
- `sampler`: The `AbstractMCMC.AbstractSampler` to test.
- `args...`: Arguments forwarded to `sample`.

# Keyword arguments
- `target`: Value to compare result of `meanfunction(chain)` to.
- `atol=1e-1`: Absolute tolerance used in `@test`.
- `rtol=1e-3`: Relative tolerance used in `@test`.
- `kwargs...`: Keyword arguments forwarded to `sample`.
"""
function test_sampler_demo_models(
    meanfunction,
    sampler::AbstractMCMC.AbstractSampler,
    args...;
    target=8.0,
    atol=1e-1,
    rtol=1e-3,
    kwargs...,
)
    @testset "$(nameof(typeof(sampler))) on $(m.name)" for model in DEMO_MODELS
        chain = AbstractMCMC.sample(model, sampler, args...; kwargs...)
        μ = meanfunction(chain)
        @test μ ≈ target atol = atol rtol = rtol
    end
end

"""
    test_sampler_continuous([meanfunction, ]sampler, args...; kwargs...)

Test that `sampler` produces the correct marginal posterior means on all models in `demo_models`.

As of right now, this is just an alias for [`test_sampler_demo_models`](@ref).
"""
function test_sampler_continuous(
    meanfunction, sampler::AbstractMCMC.AbstractSampler, args...; kwargs...
)
    return test_sampler_demo_models(meanfunction, sampler, args...; kwargs...)
end

function test_sampler_continuous(sampler::AbstractMCMC.AbstractSampler, args...; kwargs...)
    # Default for `MCMCChains.Chains`.
    return test_sampler_continuous(sampler, args...; kwargs...) do chain
        mean(Array(chain))
    end
end

end


text:
using Symbolics
using BenchmarkTools

@variables x[1:1]

m = max.(x, 0.0)
function height(x)
    [Symbolics.IfElse.ifelse(x[1] >= 0.0,
        Symbolics.IfElse.ifelse(x[1] >= 2.0, 2.0 * x[1], 1.0 * x[1]),
        0.0)]
end

m = height(x)

dm = Symbolics.jacobian(m, x)

m_func = eval(Symbolics.build_function(m, x)[1])
m_func! = eval(Symbolics.build_function(m, x)[2])
dm_func = eval(Symbolics.build_function(dm, x)[1])
dm_func! = eval(Symbolics.build_function(dm, x)[2])

a = [-1.0]

m_cache = zeros(length(m))
dm_cache = zeros(length(dm))

@benchmark m_func($a)
@benchmark m_func!($m_cache, $a)
@benchmark dm_func($a)
@benchmark dm_func!($dm_cache, $a)


text:
tron_zeros(S, n) = fill!(S(undef, Int64(n)), zero(eltype(S)))
tron_zeros(S, dims::Tuple) = fill!(S(undef, Int64(dims[1]), Int64(dims[1])), zero(eltype(S)))

export TronSparseMatrixCSC

abstract type AbstractTronMatrix end

struct TronSparseMatrixCSC{VI, VD} <: AbstractTronMatrix
    n::Int
    nnz::Int
    colptr::VI
    rowval::VI
    map::VI
    diag_vals::VD
    tril_vals::VD
end

mutable struct TronDenseMatrix{MD} <: AbstractTronMatrix
    n::Int      # the current dimension of submatrix
    max_n::Int  # the maximum dimension this matrix can hold
    vals::MD    # matrix entries
end

function TronDenseMatrix{MD}(n::Int) where {MD}
    @assert n >= 1
    vals = tron_zeros(MD, (n, n))
    return TronDenseMatrix{MD}(n, n, vals)
end

function TronDenseMatrix(A::TronDenseMatrix)
    return TronDenseMatrix{typeof(A.vals)}(A.n, A.max_n, copy(A.vals))
end

function TronDenseMatrix(I::VI, J::VI, V::Array, n) where {VI}
    @assert n >= 1
    @assert length(I) == length(J) == length(V)

    A = TronDenseMatrix{Array{Float64, 2}}(n, n, tron_zeros(Array{eltype(V)}, (n, n)))
    for i=1:length(I)
        @assert 1 <= I[i] <= n && 1 <= J[i] <= n && I[i] >= J[i]
        @inbounds A.vals[I[i], J[i]] += V[i]
    end

    return A
end

# Default constructor to allocate memory
function TronSparseMatrixCSC{VI, VD}(n::Int, nnz::Int) where {VI, VD}
    colptr = tron_zeros(VI, n+1)
    rowval = tron_zeros(VI, nnz)
    diag_vals = tron_zeros(VD, n)
    tril_vals = tron_zeros(VD, nnz)
    map = tron_zeros(VI, nnz)
    return TronSparseMatrixCSC{VI, VD}(
        n, nnz, colptr, rowval, map, diag_vals, tril_vals,
    )
end

# Constructor using sparsity pattern passed as input
function TronSparseMatrixCSC(I::VI, J::VI, V::VD, n) where {VI, VD}
    @assert length(I) == length(J) == length(V)
    # Create a CSC matrix without duplicates.
    nnz = length(I)
    # Temporary arrays for removal of duplicates
    csc_rows = zeros(Int, nnz)
    csc_cols = zeros(Int, n+1)
    imap = zeros(Int, nnz)
    w = zeros(Int, n)

    acol_ptr = tron_zeros(VI, n+1)
    arow_ind = tron_zeros(VI, nnz)
    diag_vals = tron_zeros(VD, n)
    tril_vals = tron_zeros(VD, nnz)
    map = tron_zeros(VI, nnz)

    # Count the number of entries of each column.
    for i in 1:nnz
        csc_cols[J[i]+1] += 1
    end

    # Compute the starting address (zero-based) of each column.
    for i in 2:n
        csc_cols[i+1] += csc_cols[i]
    end

    # Construct a CSC matrix.
    for i=1:nnz
        p = csc_cols[J[i]]
        csc_rows[p+1] = I[i]
        imap[p+1] = i
        csc_cols[J[i]] += 1
    end

    # Reset the starting address (one-based) of each column.
    for i=n+1:-1:2
        csc_cols[i] = csc_cols[i-1] + 1
    end
    csc_cols[1] = 1

    # Remove duplicates.
    fill!(w, -1)
    nz = 1
    for i=1:n
        p = nz
        for j=csc_cols[i]:csc_cols[i+1]-1
            r = csc_rows[j]
            if w[r] >= p  # already seen (r,i)
                if i == r # diagonal term
                    map[imap[j]] = -r
                else
                    map[imap[j]] = w[r]
                end
            else # first time seen (r,i)
                w[r] = nz
                if i == r # diagonal term
                    map[imap[j]] = -r
                else
                    map[imap[j]] = nz
                    arow_ind[nz] = r
                    nz = nz + 1
                end
            end
        end
        acol_ptr[i] = p
    end
    acol_ptr[n+1] = nz
    nnz_a = nz - 1

    return TronSparseMatrixCSC{VI, VD}(
        n, nnz_a, acol_ptr, arow_ind, map, diag_vals, tril_vals,
    )
end

Base.size(A::TronDenseMatrix) = (A.n, A.n)
Base.size(A::TronSparseMatrixCSC) = (A.n, A.n)

function Base.fill!(A::TronDenseMatrix, val)
    fill!(A.vals, val)
end

function Base.fill!(A::TronSparseMatrixCSC, val)
    fill!(A.diag_vals, val)
    fill!(A.tril_vals, val)
end

function transfer!(A::TronSparseMatrixCSC, rows, cols, values, nnz)
    for i in 1:nnz
        m = A.map[i]
        if m < 0
            A.diag_vals[-m] = values[i]
        else
            A.tril_vals[m] = values[i]
        end
    end
end
function transfer!(A::TronDenseMatrix, rows, cols, values, nnz)
    @inbounds for i in 1:nnz
        # It is assumed that rows[i] >= cols[i] for all i.
        A.vals[rows[i], cols[i]] += values[i]
    end
end

nrm2!(wa, A) = nrm2!(wa, A, A.n)

function nrm2!(wa, A::TronDenseMatrix, n)
    @inbounds for i=1:n
        wa[i] = A.vals[i,i]^2
    end
    @inbounds for j=1:n
        for i=j+1:n
            wa[j] += A.vals[i,j]^2
            wa[i] += A.vals[i,j]^2
        end
    end
    @inbounds for j=1:n
        wa[j] = sqrt(wa[j])
    end

    return
end

function nrm2!(wa, A::TronSparseMatrixCSC, n)
    @inbounds for i=1:n
        wa[i] = A.diag_vals[i]^2
    end
    @inbounds for j=1:n
        for i=A.colptr[j]:A.colptr[j+1]-1
            k = A.rowval[i]
            wa[j] = wa[j] + A.tril_vals[i]^2
            wa[k] = wa[k] + A.tril_vals[i]^2
        end
    end
    @inbounds for j=1:n
        wa[j] = sqrt(wa[j])
    end
end


text:
module HPF

# Import types from EchelleTools
import ..AbstractManifestFormatEntry
#import ..ManifestFormatEntry
import ..EchelleObservation, ..EchelleObservationSet, ..OrdersType
import ..AbstractInstrument, ..HPFType

include("constants.jl")

# Set paths # TODO: Move to config file
include("customizations.jl")
base_path = base_path_default

include("manifests.jl")
get_hpf_base_path = Manifests.get_base_path
set_hpf_base_path! = Manifests.set_base_path!
export get_hpf_base_path, set_hpf_base_path!
read_hpf_manifest_file = Manifests.read
read_hpf_manifest_files = Manifests.read
get_hpf_data_filenames = Manifests.get_data_filenames
export read_hpf_manifest_file, read_hpf_manifest_files, get_hpf_data_filenames

include("readdata.jl")
read_hpf_data = read_data
export read_hpf_data

end # module HPF


text:
function cmp(models::Vector{SampleModel}, type::Symbol)
    mnames = String[]
    lps = Matrix{Float64}[]
    for m in models
        nt = read_samples(m, :namedtuple)
        if :log_lik in keys(nt)
            append!(mnames, String.([m.name]))
            append!(lps, [Matrix(nt.log_lik')])
        else
            @warn "Model $(m.name) does not produce a log_lik matrix."
        end
    end
    (lps, type, mnames)
end

function cmp(m::Vector{Matrix{Float64}}, type::Symbol;
    mnames=String[])

    df = DataFrame()
    loo = Vector{Float64}(undef, length(m))
    loos = Vector{Vector{Float64}}(undef, length(m))
    pk = Vector{Vector{Float64}}(undef, length(m))
    for i in 1:length(m)
        loo[i], loos[i], pk[i] = psisloo(m[i])
    end
    ind = sortperm([-2loo[i][1] for i in 1:length(m)])

    mods = m[ind]
    loo = loo[ind]
    loos = loos[ind]
    pk = pk[ind]

    if length(mnames) > 0
        df.models = String.(mnames[ind])
    end

    df.PSIS = round.([-2loo[i] for i in 1:length(loo)], digits=1)
    df.SE = round.([sqrt(size(m[i], 2)*var2(-2loos[i])) for i in 1:length(m)],
        digits=2)

    dloo = zeros(length(m))
    for i in 2:length(m)
        dloo[i] = df[i, :PSIS] - df[1, :PSIS]
    end
    df.dPSIS = round.(dloo, digits=1)

    dse = zeros(length(m))
    for i in 2:length(m)
        diff = 2(loos[1] .- loos[i])
        dse[i] = √(length(loos[i]) * var2(diff))
    end
    df.dSE = round.(dse, digits=2)

    ps = zeros(length(m))
    lppds = zeros(length(m))
    for j in 1:length(m)
        n_sam, n_obs = size(mods[j])
        pd = zeros(length(m), n_obs)
        pd[j, :] = [var2(mods[j][:,i]) for i in 1:n_obs]
        ps[j] = sum(pd[j, :])
    end
    df.pPSIS = round.(ps, digits=2)

    weights = ones(length(m))
    sumval = sum([exp(-0.5df[i, :PSIS]) for i in 1:length(m)])
    for i in 1:length(m)
        weights[i] = exp(-0.5df[i, :PSIS])/sumval
    end
    df.weight = round.(weights, digits=2)
    df
end

compare([m8_1s, m8_2s], :waic) |> display

lps, type, mnames = cmp([m8_1s, m8_2s], :psis)

cmp(lps, type; mnames) |> display

md"
```

> PSIS(m8.1)
       PSIS    lppd  penalty  std_err
1 -188.6114 94.3057 2.755778 13.35391

> PSIS(m8.2)
       PSIS     lppd  penalty  std_err
1 -252.1136 126.0568 4.272533 15.29773

> compare( m8.1 , m8.2 )
       WAIC   SE dWAIC   dSE pWAIC weight
m8.2 -252.3 15.3   0.0    NA   4.3      1
m8.1 -188.7 13.3  63.5 15.15   2.7      0

> compare( m8.1 , m8.2, func=PSIS )
Some Pareto k values are high (>0.5). Set pointwise=TRUE to inspect individual points.
       PSIS    SE dPSIS   dSE pPSIS weight
m8.2 -252.1 15.28   0.0    NA   4.3      1
m8.1 -188.7 13.40  63.4 15.16   2.7      0
```
"

l, ls, pk = psisloo(m8_1s)
sum(ls)


text:
"""
    Quandl

The main access object for querying data from Quandl.com.

# Examples
```julia
quandl(TimeSeries("ML/BBY"))
quandl(TimeSeries("ML/BBY"); start_date = Date(2020,1,1), end_date = Date(2020,1,5))
quandl(TimeSeries("ML/BBY"); collapse = "weekly")
quandl(TimeSeries("ML/BBY"); transform = "diff")
quandl(TimeSeries("ML/BBY"); order = "asc")

quandl(Table("ETFG/FUND"), filters = [eq("ticker", "SPY")])
quandl(Table("ETFG/FUND"), filters = [eq("ticker", "SPY")], columns = ["ticker", "nav"])
quandl(Table("ETFG/FUND"), filters = [eq("ticker", "SPY"), gt("as_of_date", "2018-01-09")])
```
"""
struct Quandl
    api_key::String
end

Base.show(io::IO, q::Quandl) = print(io, "Quandl(api_key=******)")

"""
    TimeSeries

Represents source of a Quandl time series.
"""
struct TimeSeries
    code::String
end

"""
    Table

Represents source of a Quandl database source.
"""
struct Table
    code::String
end

"""
    FilterCondition{T}

Condition that can be used for filtering rows for a `Table` query.
"""
struct FilterCondition{T}
    op::String
    field::String
    value::T
end


text:
using Plots

using simulator.fresnelltools
include("../../src/utilities.jl")
using .utilities

n_fiber = 1.4
n_grating = 1

grating_width = 10e-9
gap_width = 400e-9

output(ν) = Grating(n_grating, n_fiber, grating_width, gap_width, 100, 2*π*ν/c_0) * [1, 0]

frequencies = 4e14:1e10:7e14
results = map(ν -> abs.(output(ν)).^2, frequencies)
transmission = map(pair -> pair[1], results)
reflection = map(pair -> pair[2], results)

plot(frequencies, transmission, label="transmission")
plot!(frequencies, reflection, label="reflection")


text:
using IteratorInterfaceExtensions
using TableTraitsUtils
import DataValues

IteratorInterfaceExtensions.isiterable(x::DataFrame) = true
TableTraits.isiterabletable(x::DataFrame) = true

function TableTraits.getiterator(df::DataFrame)
    col_names_raw = [c for c in df.columns]
    col_names = Symbol.(col_names_raw)
    column_data = [eltype(df[c]) == String ? [df[c].iloc[j] for j in 1:length(df)] : values(df[c]) for c in col_names_raw]
    return create_tableiterator(column_data, col_names)
end

TableTraits.supports_get_columns_copy_using_missing(df::DataFrame) = true

function TableTraits.get_columns_copy_using_missing(df::PandasLite.DataFrame)
    # return a named tuple of columns here
    col_names_raw = [c for c in df.columns]
    col_names = Symbol.(col_names_raw)
    cols = (Array(eltype(df[c]) == String ? [df[c].iloc[j] for j in 1:length(df)] : df[c]) for c in col_names_raw)
    return NamedTuple{tuple(col_names...)}(tuple(cols...))
end

function _construct_pandas_from_iterabletable(source)
    y = create_columns_from_iterabletable(source, errorhandling = :returnvalue)
    isnothing(y) && return nothing
    columns, column_names = y[1], y[2]
    cols = Dict{Symbol,Any}(i[1] => i[2] for i in zip(column_names, columns))
    for (k,v) in pairs(cols)
        if eltype(v) <: DataValues.DataValue
            T = eltype(eltype(v))
            if T <: AbstractFloat
                cols[k] = T[get(i, NaN) for i in v]
            elseif T <: Integer
                cols[k] = Float64[DataValues.isna(i) ? NaN : Float64(get(i)) for i in v]
            else
                throw(ArgumentError("Can't create a PandasLite.DataFrame from a source that has missing data."))
            end
        end
    end
    return DataFrame(cols)
end


text:
# Detect Elemental integer size
function ElIntType()
    using64 = Cint[0]
    ElError(ccall((:ElUsing64BitInt, libEl), Cuint, (Ptr{Cint},), using64))
    return using64[1] == 1 ? Int64 : Int32
end
const ElInt = ElIntType()

function ElCommType()
    sameSizeAsInt = Cint[0]
    ElError(ccall((:ElMPICommSameSizeAsInteger, libEl), Cuint, (Ptr{Cint},),
      sameSizeAsInt))
    return sameSizeAsInt[1] == 1 ? Cint : Ptr{Cvoid}
end
const ElComm = ElCommType()

function ElGroupType()
    sameSizeAsInt = Cint[0]
    ElError(ccall((:ElMPIGroupSameSizeAsInteger, libEl), Cuint, (Ptr{Cint},),
      sameSizeAsInt))
    return sameSizeAsInt[1] == 1 ? Cint : Ptr{Cvoid}
end
const ElGroup = ElGroupType()

# Detect Elemental Bool type
function ElBoolType()
    # NOTE: Returning Uint8 when C claims that sizeof(bool) is 1 byte leads
    #       to improperly passed structs to Elemental's C interface. This is
    #       worth investigating and might be an alignment issue.

    # warn("Hardcoding ElBool to Cint")
    boolsize = Ref(zero(Cuint))
    ElError(ccall((:ElSizeOfCBool, libEl), Cuint, (Ref{Cuint},), boolsize))
    return boolsize[] == 1 ? UInt8 : Cint
end
const ElBool = ElBoolType()

function ElBool(value::Bool)
    if value
      return ElBool(1)
    else
      return ElBool(0)
    end
end

using LinearAlgebra: BlasFloat, BlasReal, BlasComplex

const ElElementType = Union{ElInt,Float32,Float64,ComplexF32,ComplexF64}

const ElFloatType   = Union{Float32,Float64} # TODO: Maybe just use BlasReal here

abstract type ElementalMatrix{T} <: AbstractMatrix{T} end
eltype(A::ElementalMatrix{T}) where {T} = T

# Error is handled in error.jl as an Exception

@enum SortType UNSORTED DESCENDING ASCENDING
@enum Dist MC MD MR VC VR STAR CIRC
@enum Orientation NORMAL TRANSPOSE ADJOINT
@enum UpperOrLower LOWER UPPER
@enum LeftOrRight LEFT RIGHT
@enum UnitOrNonUnit NON_UNIT UNIT
@enum Pencil AXBX=1 ABX=2 BAX=3


text:
function copy_bibref_cli(argv=ARGS)

    # os
    _run_cli_os()

    # read ARGS
    argset = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table! argset begin
        "bibkey"
            help = "The entry identifier"
            arg_type = String
            required = true
        "--files", "-f"
            help = "Signal to include source files"
            action = :store_true
    end

    parsed_args = ArgParse.parse_args(argv, argset)
    bibkey = parsed_args["bibkey"]
    add_pdfs = parsed_args["files"]

    rem = _bibtorem(bibkey; add_pdfs)

    println("Rem at clipboard")
    println(rem)

    clipboard(rem)

    return nothing
end


text:
using Printf
using RecipesBase

function indices(compression :: Bool, x)

    step = compression ? 8 : 1
    1:step:length(x)

end

function indices(compression :: Int, x)

    step = compression
    1:step:length(x)

end

function solution_time(problem :: Problem, T)

    η, v, x, time = solution(problem; T = T)
    time

end
function solution_surface( problem :: Problem, T, x̃, interpolation, compression )

    η, v, x = solution(problem; x = x̃, T = T, interpolation = interpolation)
    i = indices(compression, x)
    x[i], η[i]

end

function solution_velocity( problem :: Problem, T, x̃, interpolation, compression )

    η, v, x = solution(problem; x = x̃, T = T, interpolation = interpolation)
    i = indices(compression, x)
    x[i], v[i]

end

function solution_fourier( problem :: Problem, T, x̃, interpolation, compression )

    η, v, x = solution(problem; x = x̃, T = T, interpolation = interpolation)
    fftη = fft(η)
    k = fftshift(Mesh(x).k)
    y = abs.(fftshift(fftη))
    i = indices(compression, x)
    k[i][y[i].!=0], y[i][y[i].!=0]

end

function difference( problems :: Vector{Problem}, T, x̃, interpolation )
    η1, v1, x1 = solution(problems[1]; x = x̃, T = T, interpolation = interpolation)
    η2, v2, x2 = solution(problems[2]; x = x̃, T = T, interpolation = interpolation)

    if !(x1 ≈ x2)

        if Symbol(typeof(problems[2].model)) !== :WaterWaves
            η2, v2, x2, t2 = solution(problems[2]; x = x1, T = T, interpolation = interpolation)

        elseif Symbol(typeof(problems[1].model)) !== :WaterWaves
            η1, v1, x1, t1 = solution(problems[1]; x = x2, T = T, interpolation = interpolation)

        else
            @warn("The difference is computed on different collocation points.")
            η2, v2, x2, t2 = solution(problems[2]; x = x1, T = T, interpolation = interpolation)
        end

    end
    η1, v1, η2, v2, x1

end

function difference_surface( problems :: Vector{Problem}, T, x̃, interpolation, compression )

    η1, v1, η2, v2, x = difference( problems :: Vector{Problem}, T, x̃, interpolation )
    i = indices(compression, x)
    x[i], η1[i].-η2[i]

end

function difference_velocity( problems :: Vector{Problem}, T, x̃, interpolation, compression )

    η1, v1, η2, v2, x = difference( problems :: Vector{Problem}, T, x̃, interpolation )
    i = indices(compression, x)
    x[i], v1[i].-v2[i]

end

function difference_fourier( problems :: Vector{Problem}, T, x̃, interpolation, compression )

    η1, v1, η2, v2, x = difference( problems :: Vector{Problem}, T, x̃, interpolation )
    if (x[2:end].-x[2]≈x[1:end-1].-x[1])
        k = fftshift(Mesh(x).k)

    else
        @warn("The FFT is (wrongly) computed on a non-regularly-spaced mesh.")
        k = fftshift(Mesh((xmin=x[1], xmax=x[end]+(x[end]-x[end-1]), N=length(x)) ).k)

    end
    y1 = fftshift(fft(η1))
    y2 = fftshift(fft(η2))
    y= abs.(y1.-y2)

    i = indices(compression, x)
    k[i][y[i].!=0], y[i][y[i].!=0]

end

@recipe function f(problem::Problem;
                        var = :surface,
                        T = nothing,
                        x = nothing,
                        interpolation = false,
                        compression = false)

	if var isa Array{Symbol}
		variables = var
	else
		variables = [var]
	end

    layout := (length(variables), 1)
	titlefontsize --> 10

    time = solution_time( problem, T )

	for (n, variable) in enumerate(variables)

        if n == 1
            string_title = @sprintf(" at t= %7.3f", time)
        else
            string_title = ""
        end

	    if variable == :surface

	    	@series begin

                title --> string("Surface deformation", string_title)
                label --> problem.label
                xguide --> "x"
                yguide --> "η"
	    		subplot := n

                solution_surface( problem, T, x, interpolation, compression )
	    	end

	    end

	    if variable == :velocity

	    	@series begin

                title --> string("Velocity",string_title)
                label --> problem.label
                xguide --> "x"
                yguide --> "v"
	    		subplot := n

                solution_velocity( problem, T, x, interpolation, compression )

	    	end

	    end

        if variable == :fourier || variable == :Fourier

	    	@series begin

                title --> string("Fourier coefficients (log scale)",string_title)
                label --> problem.label
                xguide --> "wavenumber"
                yguide --> "amplitude"
                yscale --> :log10
	    		subplot := n

                solution_fourier( problem, T, x, interpolation, compression )

	    	end

	    end

	end

end

@recipe function f(problems::Array{Problem};
                    var = :surface,
                    T = nothing,
                    x = nothing,
                    interpolation = false,
                    compression = false)

	if var isa Array{Symbol}
		variables = var
	else
		variables = [var]
	end

	layout := (length(variables), 1)
	titlefontsize --> 10

    time = solution_time( problems[1], T )

    for (n, variable) in enumerate(variables)

        if n == 1
            string_title = @sprintf(" at t= %7.3f", time)
        else
            string_title = ""
        end

        if variable == :surface

            for problem in problems

                @series begin

                    title --> string("Surface deformation", string_title)
                    label --> problem.label
                    xguide --> "x"
                    yguide --> "η"
                    subplot := n

                    solution_surface( problem, T, x, interpolation, compression )

                end

            end

        end

        if variable == :velocity

            for problem in problems

                @series begin

                    title --> string("Velocity", string_title)
                    label --> problem.label
                    xguide --> "x"
                    yguide --> "v"
                    subplot := n

                    solution_velocity( problem, T, x, interpolation, compression )

                end

            end

        end

        if variable in (:fourier , :Fourier)

            for problem in problems

                @series begin

                    title --> string("Fourier coefficients (log scale)", string_title)
                    label --> problem.label
                    xguide --> "wavenumber"
                    yguide --> "amplitude"
                    yscale --> :log10
                    subplot := n

                    solution_fourier( problem, T, x, interpolation, compression )

                end

            end
        end

	    if variable in ( :difference , :difference_surface)

            @series begin
                xguide --> "x"
                yguide --> "Δη"
	            label --> "$(problems[1].label) - $(problems[2].label)"
	            title --> string("Difference (surface deformation)", string_title)
                subplot := n

                difference_surface( [problems[1],problems[2]], T, x, interpolation, compression )

            end

	    end

        if variable == :difference_velocity

            @series begin
                xguide --> "x"
                yguide --> "Δv"
                label --> "$(problems[1].label) - $(problems[2].label)"
                title --> string("Difference (velocity)", string_title)
                subplot := n

                difference_velocity( [problems[1],problems[2]], T, x, interpolation, compression )

            end

        end

        if variable in ( :difference_fourier , :difference_Fourier )

            @series begin
                xguide --> "frequency"
                yguide --> "amplitude"
                label --> "$(problems[1].label) - $(problems[2].label)"
                title --> string("Difference (Fourier coefficients in log scale)", string_title)
                yscale --> :log10
                subplot := n

                difference_fourier( [problems[1],problems[2]], T, x, interpolation, compression )

            end

        end

        if variable in ( :differences, :differences_surface )

            pairs = [(i,j) for i in eachindex(problems) for j in 1:i-1]

            for (i, j) in pairs
                @series begin
                    xguide --> "x"
                    yguide --> "Δη"
                    label --> "$(problems[i].label) - $(problems[j].label)"
                    title --> string("Difference (surface deformation)", string_title)
                    subplot := n

                    difference_surface( [problems[i],problems[j]], T, x, interpolation, compression )

                end

            end

        end

        if variable == :differences_velocity

            pairs = [(i,j) for i in eachindex(problems) for j in 1:i-1]

            for (i, j) in pairs

                @series begin
                    xguide --> "x"
                    yguide --> "Δv"
                    label --> "$(problems[i].label) - $(problems[j].label)"
                    title --> string("Difference (velocity)", string_title)
                    subplot := n

                    difference_velocity( [problems[i],problems[j]], T, x, interpolation, compression )

                end

            end

        end

        if variable == :differences_fourier || variable == :differences_Fourier

            pairs = [(i,j) for i in eachindex(problems) for j in 1:i-1]

            for (i, j) in pairs

                @series begin
                    xguide --> "frequency"
                    yguide --> "amplitude"
                    label --> "$(problems[i].label) - $(problems[j].label)"
                    title --> string("Difference (Fourier coefficients in log scale)", string_title)
                    yscale --> :log10
                    subplot := n

                    difference_fourier( [problems[i],problems[j]], T, x, interpolation, compression )

                end

            end

        end

    end

end

@recipe function f(pairs::Union{Tuple{Problem, Problem},Array{Tuple{Problem, Problem}}};
                    var = :difference,
                    T = nothing,
                    x = nothing,
                    interpolation = false,
                    compression = false)

    if typeof(pairs) == Tuple{Problem, Problem} pairs = [pairs] end

    if var isa Array{Symbol}
		variables = var
	else
		variables = [var]
	end

	layout := (length(variables), 1)
	titlefontsize --> 10

    for problems in pairs

        time = solution_time( problems[1], T )

        for (n, variable) in enumerate(variables)

            if n == 1
                string_title = @sprintf(" at t= %7.3f", time)
            else
                string_title = ""
            end

            if variable in ( :difference  , :difference_surface, :surface, :surface_difference )

                @series begin
                    xguide --> "x"
                    yguide --> "Δη"
                    label --> "$(problems[1].label) - $(problems[2].label)"
                    title --> string("Difference (surface deformation)", string_title)
                    subplot := n

                    difference_surface( [problems[1],problems[2]], T, x, interpolation, compression )

                end

            end

            if variable in ( :difference_velocity, :velocity, :velocity_difference)

                @series begin
                    xguide --> "x"
                    yguide --> "Δv"
                    label --> "$(problems[1].label) - $(problems[2].label)"
                    title --> string("Difference (velocity)", string_title)
                    subplot := n

                    difference_velocity( [problems[1],problems[2]], T, x, interpolation, compression )

                end

            end

            if variable in ( :difference_fourier , :difference_Fourier , :fourier, :Fourier, :Fourier_difference, :fourier_difference )

                @series begin
                    xguide --> "frequency"
                    yguide --> "amplitude"
                    label --> "$(problems[1].label) - $(problems[2].label)"
                    title --> string("Difference (Fourier coefficients in log scale)", string_title)
                    yscale --> :log10
                    subplot := n

                    difference_fourier( [problems[1],problems[2]], T, x, interpolation, compression )

                end

            end

        end

    end

end


text:
module TwoLink
using RecipesBase

export torque, forward_kin, inverse_kin, inverse_kin_up, inverse_kin_down, traj, connect_points, acceleration, time_derivative, time_derivative!, inertia, ⊕
⊕(a::Tuple{Number,Number}, b::Tuple{Number,Number}) = (a[1]+b[1], a[2]+b[2])

const v1 = const v2 = 2
const m1 = const m2 = 0.2
const l1 = const l2 = 1
const k1 = const k2 = 0.1
const g = 1.82

function inertia(q2)
    x = cos(q2)*l1*l2*m2 + l2^2*m2
    [2*cos(q2)*l1*l2*m2+l2^2*m2+1    x;
    x                    l2^2*m2]
end
inertia(q1,q2) = inertia(q2)

signfunc(x) = tanh(100x) # A friendlier (smoother) version of sign(x)

"""
`τ1,τ2 = torque(q,qd,qdd)`\n
Inverse model
"""
function torque(q,qd,qdd)
    q1,q2     = q
    qd1,qd2   = qd
    qdd1,qdd2 = qdd

    c2  = cos(q2)
    s1  = sin(q1)
    s2  = sin(q2)
    s12 = sin(q1+q2)

    τ1 = m2*l2^2*(qdd1+qdd2) + m2*l1*l2*c2*(2qdd1+qdd2) +
    (m1+m2)*l1^2+qdd1 - m2*l1*l2*s2*qd2^2 - 2m2*l1*l2*s2*qd1*qd2 +
    g*l2*m2*s12 + (m1+m2)*l1*g*s1 + v1*qd1 +
    k1*signfunc(qd1)

    τ2 = m2*l1*l2*c2*qdd1 + m2*l1*l2*s2*qd1^2 + g*l2*m2*s12 +
    m2*l2^2*(qdd1 + qdd2) + v2*qd2 + k2*signfunc(qd2)

    (τ1,τ2)
end

function torque(q::AbstractMatrix,qd::AbstractMatrix,qdd::AbstractMatrix)
    tau = similar(q)
    for i = 1:size(q,2)
        tau1,tau2 = torque(q[:,i], qd[:,i],qdd[:,i])
        tau[1,i] = tau1
        tau[2,i] = tau2
    end
    tau
end

"""
`q̈ = acceleration(τ::Tuple, q::Tuple, qd::Tuple)`\n
`q̈1,q̈2 = acceleration(q1,q2,qd1,qd2,τ1,τ2)`\n
Model
"""
function acceleration(q1,q2,qd1,qd2,τ1,τ2)
    c2 = cos(q2)
    s1 = sin(q1)
    s2 = sin(q2)
    s12 = sin(q1+q2)

    qdd1 = (-l2*(g*l1*m1*s1 + g*l1*m2*s1 + g*l2*m2*s12 + k1*signfunc(qd1) + l1^2*m1 + l1^2*m2 - 2*l1*l2*m2*qd1*qd2*s2 - l1*l2*m2*qd2^2*s2 + qd1*v1 - τ1) + (c2*l1 + l2)*(g*l2*m2*s12 + k2*signfunc(qd2) + l1*l2*m2*qd1^2*s2 + qd2*v2 - τ2))/(l2*(2*c2*l1*l2*m2 + l2^2*m2 - m2*(c2*l1 + l2)^2 + 1))

    qdd2 = (l2*m2*(c2*l1 + l2)*(g*l1*m1*s1 + g*l1*m2*s1 + g*l2*m2*s12 + k1*signfunc(qd1) + l1^2*m1 + l1^2*m2 - 2*l1*l2*m2*qd1*qd2*s2 - l1*l2*m2*qd2^2*s2 + qd1*v1 - τ1) - (2*c2*l1*l2*m2 + l2^2*m2 + 1)*(g*l2*m2*s12 + k2*signfunc(qd2) + l1*l2*m2*qd1^2*s2 + qd2*v2 - τ2))/(l2^2*m2*(2*c2*l1*l2*m2 + l2^2*m2 - m2*(c2*l1 + l2)^2 + 1))

    qdd1,qdd2
end

function acceleration(τ::Tuple, q::Tuple, qd::Tuple)
    q1,q2   = q
    qd1,qd2 = qd
    τ1,τ2   = τ
    acceleration(q1,q2,qd1,qd2,τ1,τ2)
end

function acceleration(τ, q, q̇)
    q̈ = similar(q, dims=(2,size(q,2)))
    for i = 1:size(q,2)
        q̈1,q̈2 = acceleration(q[:,i]..., q̇[:,i]...,τ[:,i]...)
        q̈[1,i] = q̈1
        q̈[2,i] = q̈2
    end
    q̈
end

function time_derivative(state, τ)
    qdd1,qdd2 = acceleration(state[1],state[2],state[3],state[4],τ[1],τ[2])
    return [state[3],state[4], qdd1,qdd2]
end

function time_derivative!(state, τ, deriv)
    qdd1,qdd2 = acceleration(state[1],state[2],state[3],state[4],τ[1],τ[2])
    deriv[1] = state[3]
    deriv[2] = state[4]
    deriv[3] = qdd1
    deriv[4] = qdd2
end

# @syms v1 v2 k1 k2 m1 m2 l1 l2 g q1 q2 qd1 qd2 qdd1 qdd2 t1 t2 c1 c2 s1 s2 s12
# tau1 = m2*l2^2*(qdd1+qdd2) + m2*l1*l2*c2*(2qdd1+qdd2) +
# (m1+m2)*l1^2+qdd1 - m2*l1*l2*s2*qd2^2 - 2m2*l1*l2*s2*qd1*qd2 +
# g*l2*m2*s12 + (m1+m2)*l1*g*s1 + v1*qd1 +
# k1*sign(qd1)
# tau2 = m2*l1*l2*c2*qdd1 + m2*l1*l2*s2*qd1^2 + g*l2*m2*s12 +
# m2*l2^2*(qdd1 + qdd2) + v2*qd2 + k2*sign(qd2)
# sol = solve([t1-tau1,t2-tau2],[qdd1,qdd2])

function jacobian(q)
    q1,q2 = q[1],(q[1]+q[2])
    J = [l1*cos(q1)+l2*cos(q2) l2*cos(q2); l1*sin(q1)+l2*sin(q2) l2*sin(q2)]
end

function forward_kin(q)
    q1,q2 = q[1],(q[1]+q[2])
    p1 = (l1*sin(q1), -l1*cos(q1))
    p = p1 ⊕ (l2*sin(q2), -l2*cos(q2))
    p1, p
end

function forward_kin(jtraj::AbstractMatrix)
    N = size(jtraj,2)
    p = zeros(eltype(jtraj), 2, N)
    for i = 1:N
        q = jtraj[1:2,i]
        p[:,i] = forward_kin(q)[2] |> collect
    end
    p
end

"""
    (q11,q21), (q12,q22) = inverse_kin(p)
"""
function inverse_kin(p)
    x,y = p[1], p[2]
    a2 = x^2 + y^2
    a = sqrt(a2)
    α = acos((l1^2 + l2^2 - a2)/(2l1*l2))
    α2 = 2π-α
    q21 = π - α
    q22 = π - α2

    β = asin(sin(-α)/a*l2)
    β2 = -β
    q11 = asin(x/a)-β2
    q12 = asin(x/a)-β
    (q11,q21), (q12,q22)
end

function inverse_kin_up(p)
    q1,q2 = inverse_kin(p)
    return q1[2] > q2[2] ? q1 : q2
end

function inverse_kin_down(p)
    q1,q2 = inverse_kin(p)
    return q1[2] < q2[2] ? q1 : q2
end

"""
    p = inverse_kin(ctraj::Matrix, dir = :up)

`p` is of size N × 2
"""
function inverse_kin(ctraj::AbstractMatrix, dir = :up)
    if dir == :up
        kin_fun = inverse_kin_up
    elseif dir == :down
        kin_fun = inverse_kin_down
    else
        error("Symbol direction unknown")
    end
    N = size(ctraj,2)
    p = zeros(eltype(ctraj), 2, N)
    for i = 1:N
        q = ctraj[:,i]
        p[:,i] = kin_fun(q) |> collect
    end
    p
end

"""
    p, pd, pdd = traj(q0,q1,t)

Move from `q0` to `q1` during total time `t` using quadratic blends.
"""
function traj(q0,q1,t)
    tf = maximum(t)
    V = (q1-q0) ./tf * 1.5
    traj(q0,q1,t, V)
end

"""
    p, pd, pdd = traj(q0,q1,t, V)

Move from `q0` to `q1` with velocity `V` and total time `t` using quadratic blends.
"""
function traj(q0,q1,t, V)
    tf = maximum(t)
    V  = abs(V) * sign(q1-q0)
    if abs(V) < abs(q1-q0)/tf
        error("V too small")
    elseif abs(V) > 2*abs(q1-q0)/tf
        error("V too big")
    end

    if q0 == q1
        s = ones(eltype(q0),size(t)) * q0
        sd = zeros(eltype(q0),size(t))
        sdd = zeros(eltype(q0),size(t))
        return s, sd, sdd
    end

    tb = (q0 - q1 + V*tf)/V
    a = V/tb

    p = Array{eltype(q0)}(undef, size(t))
    pd = Array{eltype(q0)}(undef, size(t))
    pdd = Array{eltype(q0)}(undef, size(t))

    for (i,t) = enumerate(t)
        if t <= tb
            # initial blend
            p[i]   = q0 + a/2*t^2
            pd[i]  = a*t
            pdd[i] = a
        elseif t <= tf-tb
            # linear motion
            p[i]   = (q1+q0-V*tf)/2 + V*t
            pd[i]  = V
            pdd[i] = 0
        else
            # final blend
            p[i]   = q1 - a/2*tf^2 + a*tf*t - a/2*t^2
            pd[i]  = a*tf - a*t;
            pdd[i] = -a;
        end
    end
    return p, pd, pdd

end

"""
    p, pd, pdd = connect_points(points,ni)
Connect the points in matrix `points` ∈ R(ni+1 × 2) using third degree polinomials

# Example
```
traj = connect_points([0 1; 1 -0.5],ni)
q    = inverse_kin(traj[1],:up)'
qd   = centraldiff(q')'/sys.h
qdd  = centraldiff(qd')'/sys.h
u    = torque.(q,qd,qdd)
```
"""
function connect_points(points,ni)
    fx(i)    = traj(points[i,1],points[i+1,1],1:ni)
    fy(i)    = traj(points[i,2],points[i+1,2],1:ni)
    n        = size(points,2)
    x,xd,xdd = fx(1)
    y,yd,ydd = fy(1)
    p        = [x y]
    pd       = [xd yd]
    pdd      = [xdd ydd]
    for i = 2:n-1
        x,xd,xdd = fx(i)
        y,yd,ydd = fy(i)
        p        = vcat(p,[x y])
        pd       = vcat(pd,[xd yd])
        pdd      = vcat(pdd,[xdd ydd])
    end
    p', pd', pdd'
end

@userplot Plotworkspace
@recipe function plotworkspace(p::Plotworkspace, N=1000)
    q = [(2π*rand(),2π*rand()) for n in 1:N]
    p = getindex.(forward_kin.(q), 2)
    p1,p2 = getindex.(p,1), getindex.(p,2)
    seriesstyle := :scatter
    p1,p2
end

end


text:
x0 = 1.0
P0 = 1.0

Phi = -1.0
b = 0.0
Q = 1.0

yshadow = NaN
H = 1.0
R = 1.0
M = LinearHomogSystem(x0, P0, Phi, b, Q, yshadow, H, R)


text:
# Conversions

import SPTK: lpc2c!

function test_conversion_base(f::Function, f!::Function,
                              src_order, dst_order, args...;
                              kargs...)
    # expected dst length
    ka = Dict{Symbol,Any}(kargs)
    dst_length = haskey(ka, :dst_length) ? ka[:dst_length] : dst_order + 1

    srand(98765)
    src = rand(src_order + 1)
    dst_inplace = zeros(dst_length)
    @assert applicable(f, src, dst_order, args...)
    @assert applicable(f!, dst_inplace, src, args...)

    dst = f(src, dst_order, args...)
    f!(dst_inplace, src, args...)
    @test all(isfinite(dst))
    @test_approx_eq dst dst_inplace
end

function test_transform_base(f::Function, f!::Function, order, args...)
    srand(98765)
    src = rand(order + 1)
    dst_inplace = zeros(order + 1)
    @assert applicable(f, src, args...)
    @assert applicable(f!, dst_inplace, src, args...)

    dst = f(src, args...)
    f!(dst_inplace, src, args...)
    @test all(isfinite(dst))
    @test_approx_eq dst dst_inplace
end

function test_lpc2lsp(order)
    srand(98765)
    dummy_lpc = lpc(rand(512), order)
    lsp_inplace = zeros(order + 1)

    lsp = lpc2lsp(dummy_lpc)
    lpc2lsp!(lsp_inplace, dummy_lpc)
    @test all(isfinite(lsp))
    @test_approx_eq lsp lsp_inplace
end

function test_ic2ir_invertibility(order, len)
    srand(98765)
    dummy_ceps = rand(order+1)
    ir = c2ir(dummy_ceps, len)
    c = ic2ir(ir, order)
    @test_approx_eq c dummy_ceps
end

function test_gc2gc(src_order, dst_order, src_γ, dst_γ)
    srand(98765)
    src_ceps = rand(src_order + 1)
    dst_ceps_inplace = zeros(dst_order + 1)

    dst_ceps = gc2gc(src_ceps, src_γ, dst_order, dst_γ)
    gc2gc!(dst_ceps_inplace, dst_γ, src_ceps, src_γ)
    @test all(isfinite(dst_ceps))
    @test_approx_eq dst_ceps dst_ceps_inplace
end

function test_mgc2mgc(dst_order, dst_α, dst_γ, src_order, src_α, src_γ)
    srand(98765)
    src_ceps = rand(src_order + 1)
    dst_ceps_inplace = zeros(dst_order + 1)

    dst_ceps = mgc2mgc(src_ceps, src_α, src_γ, dst_order, dst_α, dst_γ)
    mgc2mgc!(dst_ceps_inplace, dst_α, dst_γ, src_ceps, src_α, src_γ)
    @test all(isfinite(dst_ceps))
    @test_approx_eq dst_ceps dst_ceps_inplace
end

function test_mgc2sp(order, α, γ, fftlen)
    srand(98765)
    dummy_ceps = rand(order + 1)
    sp = mgc2sp(dummy_ceps, α, γ, fftlen)
    @test length(sp) == fftlen>>1 + 1
    @test all(isfinite(sp))
end

function test_mgclsp2sp(order, α, γ, fftlen)
    srand(98765)
    dummy_lsp = rand(order + 1)
    sp = mgclsp2sp(dummy_lsp, α, γ, fftlen)
    @test length(sp) == fftlen>>1 + 1
    @test all(isfinite(sp))
end

println("-- test_lpc2c")
for src_order in [15, 20, 25, 30]
    for dst_order in [15, 20, 25, 30]
        test_conversion_base(lpc2c, lpc2c!, src_order, dst_order)
    end
end

println("-- test_lpc2lsp")
for order in [15, 20, 25, 30, 40, 50]
    println(" where order = $order")
    test_transform_base(lpc2lsp, lpc2lsp!, order)
    test_lpc2lsp(order)
end

let
    srand(98765)
    dummy_lpc = lpc(rand(512), 21)
    @test_throws ErrorException lpc2lsp(dummy_lpc, otype=2)
    @test_throws ErrorException lpc2lsp(dummy_lpc, otype=3)
    lsp1 = lpc2lsp(dummy_lpc, otype=2, fs=16000)
    lsp2 = lpc2lsp(dummy_lpc, otype=3, fs=16)
    @test_approx_eq lsp1 lsp2
    lsp3 = lpc2lsp(dummy_lpc, otype=3, fs=16, loggain=true)
    @test_approx_eq first(lsp3) log(first(lsp2))
end

let
    srand(98765)
    dummy_lpc = lpc(rand(512), 21)
    dummy_lsp_inplace = rand(50)
    @test_throws DimensionMismatch lpc2lsp!(dummy_lsp_inplace, dummy_lpc)
end

println("-- test_lpc2par")
for order in [15, 20, 25, 30]
    println(" where order = $order")
    test_transform_base(lpc2par, lpc2par!, order)
end

let
    @test_throws DimensionMismatch lpc2par!(ones(10), ones(9))
end

println("-- test_par2lpc")
for order in [15, 20, 25, 30]
    println(" where order = $order")
    test_transform_base(par2lpc, par2lpc!, order)
end

let
    @test_throws DimensionMismatch lpc2par!(ones(10), ones(9))
end

println("-- test_lsp2sp")
for order in [15, 20, 25, 30]
    for fftlen in [256, 512, 1024]
        println(" where order = $order, fftlen = $fftlen")
        test_conversion_base(lsp2sp, lsp2sp!, order, fftlen;
                             dst_length=fftlen>>1 + 1)
    end
end

let
    try lsp2sp(ones(20), 512) catch @test false end
    try lsp2sp!(ones(257), ones(20)) catch @test false end
    # invalid fftlen
    @test_throws ArgumentError lsp2sp(ones(20), 513)
    @test_throws ArgumentError lsp2sp!(ones(256), ones(20))
end

println("-- test_mc2b")
for order in [15, 20, 25, 30]
    for α in [0.35, 0.41, 0.5]
        println(" where order = $order, α = $α")
        test_transform_base(mc2b, mc2b!, order, α)
    end
end

let
    @test_throws DimensionMismatch mc2b!(ones(10), ones(9), 0.41)
end

println("-- test_b2mc")
for order in [15, 20, 25, 30]
    for α in [0.35, 0.41, 0.5]
        println(" where order = $order, α = $α")
        test_transform_base(b2mc, b2mc!, order, α)
    end
end

let
    @test_throws DimensionMismatch b2mc!(ones(10), ones(9), 0.41)
end

println("-- test_b2c")
for src_order in [15, 20, 25, 30]
    for dst_order in [15, 20, 25, 30]
        for α in [0.35, 0.41, 0.5]
            println(" where src_order = $src_order, dst_order = $dst_order, α = $α")
            test_conversion_base(b2c, b2c!, src_order, dst_order, α)
        end
    end
end

println("-- test_c2acr")
for src_order in [15, 20, 25, 30]
    for dst_order in [15, 20, 25, 30]
        for fftlen in [256, 512, 1024]
            println(" where src_order = $src_order, dst_order = $dst_order, fftlen = $fftlen")
            test_conversion_base(c2acr, c2acr!, src_order, dst_order, fftlen)
        end
    end
end

let
    # invalid fftlen
    @test_throws ArgumentError c2acr(ones(20), 19, 513)
end

println("-- test_c2ir")
for order in [15, 20, 25, 30]
    for len in [256, 512, 1024]
        println(" where order = $order, len = $len")
        test_conversion_base(c2ir, c2ir!, order, len; dst_length=len)
    end
end

println("-- test_ic2ir")
for len in [256, 512, 1024]
    for order in [15, 20, 25, 30]
        println(" where len = $len, order = $order")
        test_conversion_base(ic2ir, ic2ir!, len, order)
    end
end

println("-- test_ic2ir invertibility")
for order in [15, 20, 25, 30]
    for len in [256, 512, 1024]
        println(" where order = $order, len = $len")
        test_ic2ir_invertibility(order, len)
    end
end

println("-- test_gc2gc")
for src_order in [15, 20, 25, 30]
    for dst_order in [15, 20, 25, 30]
        for src_γ in [-1.0, -0.5, 0.0]
            for dst_γ in [-1.0, -0.5, 0.0]
                println(" where dst_order = $dst_order, dst_γ = $dst_γ, src_order = $src_order, src_γ = $src_γ")
                test_gc2gc(src_order, dst_order, src_γ, dst_γ)
            end
        end
    end
end

let
    @test_throws ArgumentError gc2gc!(ones(20), -0.1, ones(20), 0.1)
    @test_throws ArgumentError gc2gc!(ones(20), 0.1, ones(20), -0.1)
end

let
    srand(98765)
    order = 20
    src = rand(order + 1)
    @test_approx_eq gc2gc(src, 0.0, order, -0.1) gc2gc(src, 0.0, -0.1)
end

println("-- test_gnorm")
for order in [15, 20, 25, 30]
    for γ in [-1.0, -0.5, 0.0]
        println(" where order = $order, γ = $γ")
        test_transform_base(gnorm, gnorm!, order, γ)
    end
end

println("-- test_ignorm")
for order in [15, 20, 25, 30]
    for γ in [-1.0, -0.5, 0.0]
        println(" where order = $order, γ = $γ")
        test_transform_base(ignorm, ignorm!, order, γ)
    end
end

let
    c1 = ones(10)
    c2 = ones(11)
    @test_throws DimensionMismatch gnorm!(c1, c2, 0.0)
    @test_throws DimensionMismatch ignorm!(c1, c2, 0.0)
    @test_throws ArgumentError gnorm!(c1, c2, 0.1)
    @test_throws ArgumentError ignorm!(c1, c2, 0.1)
end

println("-- test_freqt and fqtr")
for src_order in [15, 20, 25, 30]
    for dst_order in [15, 20, 25, 30]
        for α in [0.35, 0.41, 0.5]
            println(" where dst_order = $dst_order, src_order = $src_order, α = $α")
            test_conversion_base(freqt, freqt!, src_order, dst_order, α)
            test_conversion_base(frqtr, frqtr!, src_order, dst_order, α)
        end
    end
end

let
    srand(98765)
    order = 20
    src = rand(order + 1)
    @test_approx_eq freqt(src, order, 0.41) freqt(src, 0.41)
    @test_approx_eq frqtr(src, order, 0.41) frqtr(src, 0.41)
end

println("-- test_mgc2mgc")
for dst_order in [15, 20, 25, 30]
    for dst_α in [-0.35, 0.41, 0.5]
        for src_γ in [-1.0, -0.5, 0.0]
            for src_order in [15, 20, 25, 30]
                for src_α in [0.35, 0.41, 0.5]
                    for dst_γ in [-1.0, -0.5, 0.0]
                        println(" where dst_order = $dst_order, dst_α = $dst_α, dst_γ = $dst_γ, src_order = $src_order, src_α = $src_α, dst_γ = $dst_γ")
                        test_mgc2mgc(dst_order, dst_α, dst_γ, src_order, src_α, src_γ)
                    end
                end
            end
        end
    end
end

let
    c1 = ones(20)
    c2 = copy(c1)
    @test_throws ArgumentError mgc2mgc!(c2, 0.0, -0.1, c1, 0.0, 1.0)
    @test_throws ArgumentError mgc2mgc!(c2, 0.0, 1.0, c1, 0.0, -0.1)
end

let
    srand(98765)
    order = 20
    src = rand(order + 1)
    @test_approx_eq mgc2mgc(src, 0.0, 0.0, order, 0.41, -0.1) mgc2mgc(src, 0.0, 0.0, 0.41, -0.1)
end

println("-- test_mgc2sp")
for order in [15, 20, 25, 30]
    for α in [0.35, 0.41, 0.5]
        for γ in [-1.0, -0.5, 0.0]
            for fftlen in [256, 512, 1024]
                println(" where order = $order, α = $α, γ = $γ, fftlen = $fftlen")
                test_mgc2sp(order, α, γ, fftlen)
            end
        end
    end
end

function mgc2sp_exceptions()
    try mgc2sp(zeros(20), 0.0, 0.0, 512) catch @test false end
    let sp = zeros(Complex{Cdouble}, 513)
        try mgc2sp!(sp, zeros(20), 0.0, 0.0) catch @test false end
    end
    # invalid fftlen
    @test_throws ArgumentError mgc2sp(zeros(20), 0.0, 0.0, 513)
    let sp = zeros(Complex{Cdouble}, 512)
        @test_throws ArgumentError mgc2sp!(sp, zeros(20), 0.0, 0.0)
    end

    # invalid γ
    @test_throws ArgumentError mgc2sp(zeros(20), 0.0, 1.0, 513)
end

mgc2sp_exceptions()

println("-- test_mgclsp2sp")

# TODO: Inf/-Inf will happens when γ = 0
for order in [15, 20, 25, 30]
    for α in [0.35, 0.41, 0.5]
        for γ in [-1.0, -0.5] # , 0.0]
            for fftlen in [256, 512, 1024]
                println(" where order = $order, α = $α, γ = $γ, fftlen = $fftlen")
                test_mgclsp2sp(order, α, γ, fftlen)
            end
        end
    end
end

let
    try mgclsp2sp(zeros(20), 0.0, 0.0, 512) catch @test false end
    try mgclsp2sp!(zeros(513), zeros(20), 0.0, 0.0) catch @test false end
    # invalid γ
    @test_throws ArgumentError mgclsp2sp!(zeros(513), zeros(20), 0.0, 1.0)
    # invalid fftlen
    @test_throws ArgumentError mgclsp2sp(zeros(20), 0.0, 0.0, 513)
    @test_throws ArgumentError mgclsp2sp!(zeros(512), zeros(20), 0.0, 0.0)
end

println("-- test_c2ndps")
for order in [15, 20, 25, 30]
    for fftlen in [256, 512, 1024]
        println(" where order = $order, fftlen = $fftlen")
        test_conversion_base(c2ndps, c2ndps!, order, fftlen,
                             dst_length=fftlen>>1 + 1)
    end
end

let
    try c2ndps(zeros(20), 512) catch @test false end
    try c2ndps!(zeros(257), zeros(20)) catch @test false end
    # invalid fftlen
    @test_throws ArgumentError c2ndps(zeros(20), 513)
    @test_throws ArgumentError c2ndps!(zeros(256), zeros(20))
end

println("-- test_ndps2c")
for order in [15, 20, 25, 30]
    for fftlen in [256, 512, 1024]
        println(" where order = $order, fftlen = $fftlen")
        test_conversion_base(ndps2c, ndps2c!, fftlen>>1, order)
    end
end

let
    try ndps2c(zeros(257), 20) catch @test false end
    # invalid fftlen
    @test_throws ArgumentError ndps2c(zeros(256), 20)
end


text:
# Unit tests for api/md_result.jl

struct MockResult1 <: MolecularDynamicsResult end
struct MockResult2 <: MolecularDynamicsResult end
Atomistic.kinetic_energy(::MockResult2, t::Integer = 0) = t * u"bohr"
Atomistic.potential_energy(::MockResult2, t::Integer = 0) = -2t * u"bohr"

@testset "api/md_result.jl" begin
    target1 = MockResult1()
    target2 = MockResult2()

    @test_throws Atomistic.UnimplementedError{MockResult1} get_system(target1)
    @test_throws Atomistic.UnimplementedError{MockResult1} get_system(target1, 1)
    @test_throws Atomistic.UnimplementedError{MockResult1} get_time_range(target1)
    @test_throws Atomistic.UnimplementedError{MockResult1} temperature(target1)
    @test_throws Atomistic.UnimplementedError{MockResult1} temperature(target1, 1)
    @test ismissing(reference_temperature(target1))
    @test_throws Atomistic.UnimplementedError{MockResult1} kinetic_energy(target1)
    @test_throws Atomistic.UnimplementedError{MockResult1} kinetic_energy(target1, 1)
    @test_throws Atomistic.UnimplementedError{MockResult1} Atomistic.potential_energy(target1)
    @test_throws Atomistic.UnimplementedError{MockResult1} Atomistic.potential_energy(target1, 1)
    @test total_energy(target2) == 0u"bohr"
    @test total_energy(target2, 1) == -1u"bohr"
    @test_throws Atomistic.UnimplementedError{MockResult1} rdf(target1)
    @test_throws Atomistic.UnimplementedError{MockResult1} rdf(target1, 0.5)
end


text:
# patching issues in base

export indextype

# https://github.com/JuliaLang/julia/issues/10178#issuecomment-74136186
println("  This warning is ok:")
function Core.Inference.func_for_method(m::Method, tt, env)
    if !m.isstaged
        return m.func.code
    else
        ## if !isleaftype(tt)
        ##     error("Need leaf-type in staged function")
        ## end
        f=ccall(:jl_instantiate_staged,Any,(Any,Any,Any),m,tt,env)
        return f.code
    end
end

# eltype for dicts
Base.eltype{K,V}(::Type{Associative{K,V}}) = V
Base.eltype{K}(::Type{Associative{K}}) = Any
Base.eltype(::Type{Associative}) = Any

indextype{K}(::Type{Associative{K}}) = K
indextype(::Type{Associative}) = Any
indextype(::Type{Associative}) = Any
indextype(t::DataType) = eltype(super(t))

indextype(::Any) = Int # fall back...

println("  endof ok-warnings.")

# iterate over Tuple. Remove after merge of https://github.com/JuliaLang/julia/pull/11547
Base.length{T<:Tuple}(t::Type{T}) = length(t.parameters)
Base.start{T<:Tuple}(::Type{T}) = 1
Base.next{T<:Tuple}(t::Type{T}, state) = (t.parameters[state], state+1)
Base.done{T<:Tuple}(t::Type{T}, state) = length(t)<state

# indexing into Tuple
Base.getindex{T<:Tuple}(t::Type{T}, i::Integer) = t.parameters[i]


text:
abstract type Schemes end
struct UpWind <: Schemes end
struct Second <: Schemes end
struct Third <: Schemes end
struct Fourth <: Schemes end

mutable struct FVMPDEOption

    Δt
    Δx
    Δs
    function HPDEOption(;kwargs...)
        Δt, Δx, Δs = nothing,nothing,nothing
        if haskey(kwargs,:dt)
            Δt = kwargs[:dt]
        end
        if haskey(kwargs,:dx)
            Δx = kwargs[:dx]
        end
        if haskey(kwargs,:ds)
            Δs = kwargs[:ds]
        end
        new(Δt, Δx, Δs)
    end
end

HPDEOption()
HPDEOption(dt=0.1)
HPDEOption(dx=0.1)
HPDEOption(ds=0.1)
HPDEOption(ds=0.1,dx=0.2)
HPDEOption(ds=0.1,dx=0.2,dt=0.3)


text:
# Author: Vincent Arel-Bundock
# Contact: varel@umich.edu
# License: BSD-3
# Original code from the python statsmodels project
# https://github.com/statsmodels/statsmodels

function qreg_coef(y, X::Matrix, q::Real, method::IRLS)
    n, p = size(X)
    xstar = copy(X)
    diff = Inf

    beta0 = zeros(p)
    beta  = zeros(p)
    xtx   = Array{Float64}(undef, p, p)
    xty   = Array{Float64}(undef, p)
    xbeta = Array{Float64}(undef, n)
    resid = Array{Float64}(undef, n)

    for itr in 1:method.maxiter
        if diff > method.tol
            copyto!(beta0, beta)

            mul!(xtx, xstar', X)
            mul!(xty, xstar', y)
            beta .= xtx \ xty
            mul!(xbeta, X, beta)

            @. resid = y - xbeta

            for i in 1:n
                if abs(resid[i]) < method.threshold
                    @inbounds resid[i] = sign(resid[i]) * method.threshold
                end
                if resid[i] < 0
                    @inbounds resid[i] = abs(q * resid[i])
                else
                    @inbounds resid[i] = abs((1 - q) * resid[i])
                end
            end

            xstar .= X ./ resid

            diff = norm(beta0 - beta, Inf)
        end
    end
    return beta
end


text:
# ============================================================
# General Structs
# ============================================================

"""
$(TYPEDEF)

Computational setup

# Fields

$(FIELDS)
"""
@kwdef struct Setup{S,I<:Integer,E<:AV,F<:Real,G<:Real} <: AbstractSetup
    matter::S = "gas"
    case::S = "dev"
    space::S = "1d0f0v"
    flux::S = "kfvs"
    collision::S = "bgk"
    nSpecies::I = 1
    interpOrder::I = 2
    limiter::S = "vanleer"
    boundary::E = ["fix", "fix"]
    cfl::F = 0.5
    maxTime::G = 0.1
end

function Setup(
    matter,
    case,
    space,
    flux,
    collision,
    ns,
    order,
    limiter,
    bc::T,
    cfl,
    time,
) where {T<:Union{AbstractString,Symbol}}
    boundary = begin
        if parse(Int, space[1]) == 1
            [bc, bc]
        elseif parse(Int, space[1]) == 2
            [bc, bc, bc, bc]
        end
    end

    return Setup{typeof(matter),typeof(ns),typeof(boundary),typeof(cfl),typeof(time)}(
        matter,
        case,
        space,
        flux,
        collision,
        ns,
        order,
        limiter,
        boundary,
        cfl,
        time,
    )
end