text:
# Define an abstract type to represent color difference metrics
abstract type DifferenceMetric end

abstract type EuclideanDifferenceMetric{T<:Color3} <: DifferenceMetric end

# TODO?: make the DifferenMetrics parametric, to preserve type-stability

struct DE_2000 <: DifferenceMetric
    kl::Float64
    kc::Float64
    kh::Float64
end
"""
    DE_2000()
    DE_2000(kl, kc, kh)

Construct a metric of the CIE Delta E 2000 recommendation, with weighting
parameters `kl`, `kc` and `kh` as provided for in the recommendation. When not
provided, these parameters default to `1`.
"""
DE_2000() = DE_2000(1,1,1)

struct DE_94 <: DifferenceMetric
    kl::Float64
    kc::Float64
    kh::Float64
end
"""
    DE_94()
    DE_94(kl, kc, kh)

Construct a metric of CIE Delta E 94 recommendation (1994), with weighting
parameters `kl`, `kc` and `kh` as provided for in the recommendation. When not
provided, these parameters default to `1`.
The `DE_94` is more perceptually uniform than the [`DE_AB`](@ref), but has some
non-uniformities resolved by the [`DE_2000`](@ref).
"""
DE_94() = DE_94(1,1,1)

struct DE_JPC79 <: DifferenceMetric

end
"""
    DE_JPC79()

Construct a metric using McDonald's "JP Coates Thread Company" color difference
formula.
"""
DE_JPC79()

struct DE_CMC <: DifferenceMetric
    kl::Float64
    kc::Float64
end
"""
    DE_CMC()
    DE_CMC(kl, kc)

Construct a metric using the CMC equation (CMC l:c), with weighting parameters
`kl` and `kc`. When not provided, these parameters default to `1`.
!!! note
    The `DE_CMC` is a quasimetric, i.e. violates symmetry. Therefore,
    `colordiff(a, b, metric=DE_CMC())` may not equal to
    `colordiff(b, a, metric=DE_CMC())`.
"""
DE_CMC() = DE_CMC(1,1)

struct DE_BFD <: DifferenceMetric
    wp::XYZ{Float64}
    kl::Float64
    kc::Float64
end
"""
    DE_BFD()
    DE_BFD([wp,] kl, kc)

Construct a metric using the BFD equation, with weighting parameters `kl` and
`kc`. Additionally, a whitepoint `wp` can be specified, because the BFD equation
must convert between `XYZ` and `Lab` during the computation. When not provided,
`kl` and `kc` default to `1`, and `wp` defaults to CIE D65 (`Colors.WP_D65`).
"""
DE_BFD() = DE_BFD(WP_DEFAULT,1,1)
DE_BFD(kl, kc) = DE_BFD(WP_DEFAULT,kl, kc)

struct DE_AB <: EuclideanDifferenceMetric{Lab}

end
"""
    DE_AB()

Construct a metric of the original CIE Delta E equation (ΔE*ab), or Euclidean
color difference equation in the `Lab` (CIELAB) colorspace.
"""
DE_AB()

struct DE_DIN99 <: EuclideanDifferenceMetric{DIN99}

end
"""
    DE_DIN99()

Construct a metric using Euclidean color difference equation applied in the
`DIN99` colorspace.
"""
DE_DIN99()

struct DE_DIN99d <: EuclideanDifferenceMetric{DIN99d}

end
"""
    DE_DIN99d()

Construct a metric using Euclidean color difference equation applied in the
`DIN99d` colorspace.
"""
DE_DIN99d()

struct DE_DIN99o <: EuclideanDifferenceMetric{DIN99o}

end
"""
    DE_DIN99o()

Construct a metric using Euclidean color difference equation applied in the
`DIN99o` colorspace.
"""
DE_DIN99o()

# Compute the mean of two hue angles
function mean_hue(h1, h2)
    if abs(h2 - h1) > 180
        if h1 + h2 < 360
            mh = (h1 + h2 + 360) / 2
        else
            mh = (h1 + h2 - 360) / 2
        end
    else
        mh = (h1 + h2) / 2
    end

    mh
end

# Color difference metrics
# ------------------------

# Evaluate the CIEDE2000 color difference formula, implemented according to:
#   Klaus Witt, CIE Color Difference Metrics, Colorimetry: Understanding the CIE
#   System. 2007
#
# Args:
#   a, b: Any two colors.
#
# Returns:
#   The CIEDE2000 color difference metric evaluated between a and b.
#

pow7(x) = (y = x*x*x; y*y*x)
pow7(x::Integer) = pow7(Float64(x))
const twentyfive7 = pow7(25)

# Delta E 2000
function _colordiff(ai::Color, bi::Color, m::DE_2000)
    # Ensure that the input values are in L*a*b* space
    a_Lab = convert(Lab, ai)
    b_Lab = convert(Lab, bi)

    # Calculate some necessary factors from the L*a*b* values
    mc = (chroma(a_Lab) + chroma(b_Lab))/2
    g = (1 - sqrt(pow7(mc) / (pow7(mc) + twentyfive7))) / 2

    # Convert to L*C*h, where the remainder of the calculations are performed
    a = convert(LCHab, Lab(a_Lab.l, a_Lab.a * (1 + g), a_Lab.b))
    b = convert(LCHab, Lab(b_Lab.l, b_Lab.a * (1 + g), b_Lab.b))

    # Calculate the delta values for each channel
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end
    # Calculate H*
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    # Calculate mean L* and C* values
    ml, mc = (a.l + b.l) / 2, (a.c + b.c) / 2

    # Calculate mean hue value
    if a.c * b.c == 0
        mh = a.h + b.h
    else
        mh = mean_hue(a.h, b.h)
    end

    # lightness weight
    mls = (ml - 50)^2
    sl = 1.0 + 0.015 * mls / sqrt(20 + mls)

    # chroma weight
    sc = 1 + 0.045mc

    # hue weight
    t = 1 - 0.17 * cosd(mh - 30) +
            0.24 * cosd(2mh) +
            0.32 * cosd(3mh + 6) -
            0.20 * cosd(4mh - 63)
    sh = 1 + 0.015 * mc * t

    # rotation term
    dtheta = 30 * exp(-((mh - 275)/25)^2)
    cr = 2 * sqrt(pow7(mc) / (pow7(mc) + twentyfive7))
    tr = -sind(2*dtheta) * cr

    # Final calculation
    sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/(m.kh*sh))^2 +
         tr * (dc/(m.kc*sc)) * (dh/(m.kh*sh)))
end

# Delta E94
function _colordiff(ai::Color, bi::Color, m::DE_94)

    a = convert(LCHab, ai)
    b = convert(LCHab, bi)

    # Calculate the delta values for each channel
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end
    # Calculate H*
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    # Calculate geometric mean of chroma
    mc = sqrt(a.c*b.c)

    # Lightness, hue, chroma correction terms
    sl = 1
    sc = 1+0.045*mc
    sh = 1+0.015*mc

    sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/(m.kh*sh))^2)

end

# Metric form of jpc79 color difference equation (mostly obsolete)
function _colordiff(ai::Color, bi::Color, m::DE_JPC79)

    # Convert directly into LCh
    a = convert(LCHab, ai)
    b = convert(LCHab, bi)

    # Calculate deltas in each direction
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end
    # Calculate H* from C*'s and h
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    #Calculate mean lightness
    ml = (a.l + b.l)/2
    mc = (a.c + b.c)/2
    # Calculate mean hue value
    if a.c * b.c == 0
        mh = a.h + b.h
    else
        mh = mean_hue(a.h, b.h)
    end

    # L* adjustment term
    sl = 0.08195*ml/(1+0.01765*ml)

    # C* adjustment term
    sc = 0.638+(0.0638*mc/(1+0.0131*mc))

    # H* adjustment term
    if (mc < 0.38)
        sh = sc
    elseif (mh >= 164 && mh <= 345)
        sh = sc*(0.56+abs(0.2*cosd(mh+168)))
    else
        sh = sc*(0.38+abs(0.4*cosd(mh+35)))
    end

    # Calculate the final difference
    sqrt((dl/sl)^2 + (dc/sc)^2 + (dh/sh)^2)

end

# Metric form of the cmc color difference
function _colordiff(ai::Color, bi::Color, m::DE_CMC)

    # Convert directly into LCh
    a = convert(LCHab, ai)
    b = convert(LCHab, bi)

    # Calculate deltas in each direction
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end
    # Calculate H* from C*'s and h
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    # Find the mean value of the inputs to use as the "standard"
    ml, mc = (a.l + b.l)/2, (a.c + b.c)/2 # TODO: use `a.l` and `a.c` instead

    # Calculate mean hue value
    if a.c * b.c == 0
        mh = a.h + b.h
    else
        mh = mean_hue(a.h, b.h)
    end

    # L* adjustment term
    if (a.l <= 16)
        sl = 0.511
    else
        sl = 0.040975*ml/(1+0.01765*ml)
    end

    # C* adjustment term
    sc = 0.0638*mc/(1+0.0131*mc)+0.638

    f = sqrt((mc^4)/(mc^4 + 1900))

    if (mh >= 164 && mh < 345)
        t = 0.56 + abs(0.2*cosd(mh+168))
    else
        t = 0.36 + abs(0.4*cosd(mh+35))
    end

    # H* adjustment term
    sh = sc*(t*f+1-f)

    # Calculate the final difference
    sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/sh)^2)

end

# The BFD color difference equation
function _colordiff(ai::Color, bi::Color, m::DE_BFD)
    # Currently, support for the `wp` argument of `convert` is limited.
    function to_xyz(c::Color, wp)
        c isa XYZ && return ai
        (c isa xyY || c isa LMS) && return convert(XYZ, c)
        (c isa Lab || c isa Luv) && return convert(XYZ, c, wp)
        c isa LCHuv && return convert(XYZ, convert(Luv, c), wp)
        convert(XYZ, convert(Lab, c), wp)
    end

    # We have to start back in XYZ because BFD uses a different L equation
    a_XYZ = to_xyz(ai, m.wp)
    b_XYZ = to_xyz(bi, m.wp)

    la = 54.6*log10(a_XYZ.y+1.5)-9.6
    lb = 54.6*log10(b_XYZ.y+1.5)-9.6

    # Convert into LCh with the proper white point
    a_Lab = convert(Lab, a_XYZ, m.wp)
    b_Lab = convert(Lab, b_XYZ, m.wp)

    # Substitute in the different L values into the L*C*h values
    a = LCHab(la, chroma(a_Lab), hue(a_Lab))
    b = LCHab(lb, chroma(b_Lab), hue(b_Lab))

    # Calculate deltas in each direction
    dl, dc, dh = (b.l - a.l), (b.c - a.c), (b.h - a.h)
    if a.c * b.c == 0
        dh = zero(dh)
    elseif dh > 180
        dh -= 360
    elseif dh < -180
        dh += 360
    end

    # Calculate H* from C*'s and h
    dh = 2 * sqrt(a.c * b.c) * sind(dh/2)

    # Find the mean value of the inputs to use as the "standard"
    ml, mc = (a.l + b.l)/2, (a.c + b.c)/2

    # Calculate mean hue value
    if a.c * b.c == 0
        mh = a.h + b.h
    else
        mh = mean_hue(a.h, b.h)
    end

    # Correction terms for a variety of nonlinearities in CIELAB.
    g = sqrt(mc^4/(mc^4 + 14000))

    t = 0.627 + 0.055*cosd(mh - 245) - 0.04*cosd(2*mh - 136) + 0.07*cosd(3*mh - 32) + 0.049*cosd(4*mh + 114) - 0.015*cosd(5*mh + 103)

    rc = sqrt(mc^6/(mc^6 + 7e7))

    rh = -0.26cosd(mh-308) - 0.379cosd(2*mh-160) - 0.636*cosd(3*mh - 254) + 0.226cosd(4*mh + 140) - 0.194*cosd(5*mh + 280)

    dcc = 0.035*mc/(1+0.00365*mc) + 0.521
    dhh = dcc*(g*t+1-g)
    rt = rc*rh

    # Final calculation
    sqrt((dl/m.kl)^2 + (dc/(m.kc*dcc))^2 + (dh/dhh)^2 + rt*((dc*dh)/(dcc*dhh)))
end

function _colordiff(ai::Color, bi::Color,
                    m::EuclideanDifferenceMetric{T}) where {T <: Color3}
    a, b = convert(T, ai), convert(T, bi)

    d1, d2, d3 = comp1(a) - comp1(b), comp2(a) - comp2(b), comp3(a) - comp3(b)

    sqrt(d1^2 + d2^2 + d3^2)
end

# Default to Delta E 2000
"""
    colordiff(a, b; metric=DE_2000())

Compute an approximate measure of the perceptual difference between colors `a`
and `b`. Optionally, a `metric` may be supplied, chosen among [`DE_2000`](@ref)
(the default), [`DE_94`](@ref), [`DE_JPC79`](@ref), [`DE_CMC`](@ref),
[`DE_BFD`](@ref), [`DE_AB`](@ref), [`DE_DIN99`](@ref), [`DE_DIN99d`](@ref) and
[`DE_DIN99o`](@ref).

The return value is a non-negative number in a type depending on the colors and
metric.

!!! note
    The supported metrics measure the difference within `Lab` or its variant
    colorspaces. When the input colors are not in the colorspace internally used
    by the metric, the colors (e.g. in `RGB`) are converted with the default
    whitepoint CIE D65 (`Colors.WP_D65`). If you want to use another whitepoint,
    convert the colors into the colorspace used by metric (e.g. `Lab` for
    [`DE_2000`](@ref)) in advance.
"""
colordiff(ai::Union{Number, Color},
          bi::Union{Number, Color};
          metric::DifferenceMetric=DE_2000()) = _colordiff(ai, bi, metric)
@deprecate colordiff(ai::Color, bi::Color, metric::DifferenceMetric) colordiff(ai, bi; metric=metric)

function colordiff(ai::Colorant, bi::Colorant; metric::DifferenceMetric=DE_2000())
    alpha(ai) == 1 && alpha(bi) == 1 && return _colordiff(color(ai), color(bi), metric)
    throw(ArgumentError("""
        cannot evaluate the difference in transparent colors.
          Their appearance depends on the backdrop."""))
end

_colordiff(ai::AbstractGray, bi::Number, metric::DifferenceMetric) = _colordiff(ai, Gray(bi), metric)
_colordiff(ai::Number, bi::AbstractGray, metric::DifferenceMetric) = _colordiff(Gray(ai), bi, metric)
_colordiff(ai::Number, bi::Number, metric::DifferenceMetric) = _colordiff(Gray(ai), Gray(bi), metric)


text:
using Test

@testset "Test LogisticRegression.jl" begin
    @testset "test_sigmoid.jl" begin
        include("test_sigmoid.jl")
    end
    @testset "test_propagate.jl" begin
        include("test_propagate.jl")
    end
    @testset "test_optimize.jl" begin
        include("test_optimize.jl")
    end
    @testset "test_predict.jl" begin
        include("test_predict.jl")
    end
end

ext: jl


text:
"""
    read_block_by_name(filename::String, blockname::String;
                                info::Info_Line=Info_Line(),
                                parttype::Int64=-1)

Reads a block in a snapshot with given name. Names are case sensitive.

# Examples
```jldoctest
julia> pos_info = Info_Line("POS", Float32, 1, [1, 1, 1, 1, 1, 1])
[...]
julia> gas_pos = read_block_by_name(filename, "POS", info=pos_info, parttype=0)
[...]
```
"""
function read_block_by_name(filename::String, blockname::String;
                            info::Info_Line=Info_Line(),
                            parttype::Int64=-1)

    # read header - super fast and needed for flexibility
    h = head_to_obj(filename)

    if parttype == -1
        parttypes = ["PartType0", "PartType1", "PartType2",
                     "PartType3", "PartType4", "PartType5"]
        d = Dict()
    end # if parttype

    blockname = strip(blockname)

    if info.block_name == ""
        info = read_info(filename)
        if info == 1
            if blockname == "MASS"
                info = Info_Line("MASS", Float32, 1, [0, 0, 0, 0, 0, 0])
            else
                error("No Info block in snapshot! Supply Info_Line type!")
            end
        else # info != 1
            for i ∈ 1:length(info)
                if info[i].block_name == blockname
                    info = info[i]
                    break
                end # if block found
            end # loop over info
            if isa(info, Array)
                if (blockname == "MASS")
                    info = Info_Line("MASS", Float32, 1, [0, 0, 0, 0, 0, 0])
                else
                    error("Block not present!")
                end
            end
        end # info == 1
    end

    f = open(filename)

    first_block = read(f, Int32)

    if first_block != Int32(8)
        error("Only snapshots of format 2 can be read by name!")
        return
    end

    while eof(f) != true
        name = Char.(read!(f, Array{Int8,1}(undef,4)))
        blname = String(name)
        blname = strip(blname)

        if blname == blockname
            break
        end
        p = position(f)
        seek(f,p+8)
        skipsize = read(f, Int32)
        p = position(f)
        seek(f,p+skipsize+8)
    end

    if eof(f) != true
        p = position(f)
        seek(f,p+8)

        blocksize = read(f,Int32)

        for i ∈ 1:length(info.is_present)
            p = position(f)

            if info.is_present[i] == Int32(1)
                if parttype == -1
                    d[parttypes[i]] = Dict()
                    d[parttypes[i]][blockname] = copy(transpose(
                                                        read!(f, Array{info.data_type,2}(undef,(info.n_dim,h.npart[i])))))
                else
                    if i == (parttype+1)

                        block = copy(transpose(
                                        read!(f, Array{info.data_type,2}(undef,(info.n_dim,h.npart[i])))))
                        close(f)
                        return block
                    else
                        seek(f, p + ( sizeof(info.data_type)*info.n_dim*h.npart[i] ))
                    end # if i == (parttype+1)
                end # parttype == -1
            end # info.is_present[i] == Int32(1)
        end # i ∈ 1:length(info.is_present)

        # fill in mass from header
        if blockname == "MASS"
            for i ∈ 1:length(info.is_present)
                if info.is_present[i] == Int32(0)
                    if parttype == -1
                        d[parttypes[i]] = Dict()
                        d[parttypes[i]][blockname] = Array{info.data_type,2}(undef,(info.n_dim,h.npart[i]))
                        d[parttypes[i]][blockname] .= h.massarr[i]
                    end # parttype == -1
                end # info block
            end # loop over i
        end
        close(f)

        return d

    else # blockname not found
        if blockname != "MASS"
            error("Block not present!")
        else
            if parttype == -1
                for i = 1:length(parttypes)
                    d[parttypes[i]] = Dict()
                    d[parttypes[i]][blockname] = Array{info.data_type,2}(undef,(info.n_dim,h.npart[i]))
                    d[parttypes[i]][blockname] .= h.massarr[i]
                end
            else
                block = Array{info.data_type,2}(undef,(info.n_dim,h.npart[parttype+1]))
                block .= h.massarr[parttype+1]
                return block
            end # parttype == -1
        end
    end # eof(f) != true

end

function snap_2_d_info(filename::String, d::Dict{Any,Any}, info::Array{Info_Line,1})

    f = open(filename)

    seek(f, 300)

    for i ∈ 1:length(d["Header"]["npart"])
        if d["Header"]["npart"][i] != Int32(0)
            d[d["Header"]["PartTypes"][i]] = Dict()
            d[d["Header"]["PartTypes"][i]]["MASS"] = Float32.(d["Header"]["massarr"][i] .* ones(d["Header"]["npart"][i]))
        end
    end

    for i ∈ 1:length(info)
        #println("Reading block: ", info[i].block_name)
        d = read_block_with_info(f, d, info[i])
        p = position(f)
        seek(f, p+24)
    end

    close(f)

    return d

end

function snap_2_d(filename::String, data::Dict{Any,Any})

    f = open(filename)

    seek(f, 296)

    N = sum(data["Header"]["npart"])
    skipsize = read(f, Int32)
    bit_size = Int64(skipsize/(3*N))

    if Int(bit_size) == 4
        dtype = Float32
        println("Reading single precision snapshot")
    elseif Int(bit_size) == 8
        dtype = Float64
        println("Reading double precision snapshot")
    else
        println("read error! neither 32 nor 64 bits data!")
        return -1
    end

    seek(f, 288)

    # set up dictionaries for particles
    for i in 1:length(data["Header"]["PartTypes"])
        if data["Header"]["npart"][i] != 0
            data[data["Header"]["PartTypes"][i]] = Dict()
        end
        if data["Header"]["massarr"][i] != 0
           data[data["Header"]["PartTypes"][i]]["MASS"] = dtype.(data["Header"]["massarr"][i] * ones(data["Header"]["npart"][i]))
        end
    end

    blockname = "POS"

    while (blockname != "INFO") #& (eof(f) != true)

        #println(blockname)
        #println(typeof(blockname))
        # skip identifiers
        p = position(f)

        seek(f, p+8)

        if eof(f) == true
            break
        end

        p, data = read_block(p, data, dtype, blockname, f, bit_size)

        seek(f,p+4)

        if eof(f) == true
            break
        end

        p = position(f)
        seek(f, p+4)
        # read blockname

        name = Char.(read!(f, Array{Int8,1}(undef,4)))
        blockname = String(name)

        blockname = String(strip(blockname))

        p = position(f)
        #seek(f, p+8)

    end

    close(f)

    return data

end

function read_block_with_info(f::IOStream, d::Dict{Any,Any}, info::Info_Line)

    parttypes = ["PartType0", "PartType1", "PartType2",
                 "PartType3", "PartType4", "PartType5"]

    for i ∈ 1:length(info.is_present)

        if info.is_present[i] == Int32(1)
            d[parttypes[i]][info.block_name] = copy(transpose(
                                         read!(f, Array{info.data_type,2}(undef,(info.n_dim,d["Header"]["npart"][i])))))
        end

    end

    return d

end

function read_block(p::Int64, data::Dict{Any,Any}, dtype::DataType, blockname::String,
                    f::IOStream, bit_size::Int64)

    file_curr = @__FILE__
    path_curr = file_curr[1:end-17]

    include(path_curr * "part_specific_fields.jl")

    N = sum(data["Header"]["npart"])

    skipsize = read(f, Int32)

    if blockname == "MASS"
        for i ∈ 1:length(data["Header"]["PartTypes"])

            n = data["Header"]["npart"][i]

            if data["Header"]["npart"][i] != Int32(0)

                if data["Header"]["massarr"][i] == Int32(0)
                    data[data["Header"]["PartTypes"][i]][blockname] = copy(transpose(read!(f, Array{dtype,2}(undef,(1,n)))))
                else
                    data[data["Header"]["PartTypes"][i]][blockname] =  Float32.(data["Header"]["massarr"][i] .* ones(n))
                end

            end
        end
    elseif blockname != "ID"

        gas_block = any(x->x==blockname,gas_arr)

        if gas_block == true

            #println("Reading gas block")

            n = Int64(data["Header"]["npart"][1])

            dim = Int(skipsize/(n*bit_size))

            #println("nr. of gas-particles: ", n)
            #println("reading dimensions: ", skipsize/(n*bit_size))

            data[data["Header"]["PartTypes"][1]][blockname] = copy(transpose(read!(f, Array{dtype,2}(undef,(dim,n)))))

        # check if bh
        elseif (any(x->x==blockname,bh_arr)) == true

            n = Int64(data["Header"]["npart"][6])

            #println("nr. of bh-particles: ", n)

            dim = Int(skipsize/(n*bit_size))

            #println("reading dimensions: ", skipsize/(n*bit_size))

            if blockname != "BHPC"
                data[data["Header"]["PartTypes"][6]][blockname] = copy(transpose(read!(f, Array{dtype,2}(undef,(dim,n)))))
            else
                data[data["Header"]["PartTypes"][6]][blockname] = Int.(copy(transpose(read!(f, Array{UInt32,2}(undef,(1,n))))))
            end

        elseif blockname == gas_bh

            n = Int64(data["Header"]["npart"][1])
            data[data["Header"]["PartTypes"][1]][blockname] = Int.(copy(transpose(read!(f, Array{UInt32,2}(undef,(1,n))))))

            n = Int64(data["Header"]["npart"][6])
            data[data["Header"]["PartTypes"][6]][blockname] = Int.(copy(transpose(read!(f, Array{UInt32,2}(undef,(1,n))))))

        else
            # read Blocks
            #println("Blockname: ", blockname)

            for i in 1:length(data["Header"]["PartTypes"])

                dim = (skipsize/(N*bit_size))

                dim = Int64(trunc(dim))

                if data["Header"]["npart"][i] != 0

                    n = Int64(data["Header"]["npart"][i])

                    data[data["Header"]["PartTypes"][i]][blockname] = copy(transpose(read!(f, Array{dtype,2}(undef,(dim,n)))))

                end

            end
        end

    else

        for i in 1:length(data["Header"]["PartTypes"])

            if data["Header"]["npart"][i] != 0

                n = Int64(data["Header"]["npart"][i])

                data[data["Header"]["PartTypes"][i]][blockname] = copy(transpose(read!(f, Array{UInt32,2}(undef,(1,n)))))

            end

        end

    end

    p = position(f)

    return p, data
end

function read_block_with_offset(filename::String, data_old, pos0::Integer, info::Info_Line,
                                offset::Integer, offset_key, n_to_read::Integer, part_per_key )

    # open the file
    f = open(filename)

    # number of bits in data_type
    len = sizeof(info.data_type) * info.n_dim

    # jump to position of particle type in relevant block
    seek(f, pos0+offset*len)

    # store position in file
    p = position(f)

    # allocate array to store data
    data =  Array{info.data_type,2}(undef, (n_to_read, info.n_dim))

    n_read = 1
    n_this_key = 0

    for i = 1:length(offset_key)

        # jump to start of key
        seek(f, p + len*offset_key[i])
        n_this_key += part_per_key[i]

        data[n_read:n_this_key, :] = copy(transpose(read!(f,
                    Array{info.data_type,2}(undef,(info.n_dim,part_per_key[i])))))

        n_read += part_per_key[i]

    end # for

    # close the file
    close(f)

    return [data_old; data]
end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
# the same with: Dict("A"=>1, "B"=>2)
D = Dict([("A", 1), ("B", 2)])

# has key
println(haskey(D, "A"))

# get(dict, key, default)
println(get(D, "A", 3))

# get!(dict, key, default)
# Return the value stored for the given key, or if no mapping for the key is present, store key => default, and return default.
println(get!(D, "C", 3))
println(get(D, "C", 5))

# delete!(dict, key)
# Delete the mapping for the given key in a collection, and return the collection.

# keys & values
println(keys(D))

# merge
D1 = Dict([("A", 1), ("B", 2)])
D2 = Dict([("C", 2), ("B", 3)])
println(merge(D1, D2))

# WeakKeyDict is a hash table implementation where the keys are weak references to objects, and thus may be garbage collected even when referenced in a hash table. Like Dict it uses hash for hashing and isequal for equality, unlike Dict it does not convert keys on insertion.

# ImmutableDict is a Dictionary implemented as an immutable linked list, which is optimal for small dictionaries that are constructed over many individual insertions Note that it is not possible to remove a value, although it can be partially overridden and hidden by inserting a new value with the same key

max_stars_repo_name: startprogress/Hello-Julia

max_issues_repo_name: startprogress/Hello-Julia

max_forks_repo_name: startprogress/Hello-Julia


text:
if Sys.iswindows()
# ccall Fortran function "mul(x, y)"
println(ccall(:mul, Float64, (Float64, Float64), 2.0, 3.0))

# ccall Fortran function "sub(x, y)"
println(ccall(:sub, Float64, (Float64, Float64), 10.0, 2.0))
end

max_stars_repo_name: yundantianchang/dynload-julia

max_issues_repo_name: yundantianchang/dynload-julia

max_forks_repo_name: yundantianchang/dynload-julia


text:
using DocOpt

doc="""Voice conversion based on differencial spectral compensation

Usage:
    diffvc.jl [options] <input_wav> <model_jld> <dst_wav>
    diffvc.jl --version
    diffvc.jl -h | --help

Options:
    -h --help        show this message
    --period=PERIOD  frame period in msec [default: 5.0]
    --order=ORDER    order of mel cepsrum [default: 40]
    --alpha=ALPHA    all-pass constant [default: 0.0]
    --trajectory     trajectory-based parameter conversion
    --T=t            maximum length for one-step trajectory conversion [default: 100]
"""

using VoiceConversion
using MelGeneralizedCepstrums
using WAV
using SynthesisFilters
using HDF5, JLD

import WORLD: synthesis
using WORLD

using Logging
@Logging.configure(level=DEBUG, output=STDOUT)

function main()
    args = docopt(doc, version=v"0.0.1")

    x, fs = wavread(args["<input_wav>"])
    @assert size(x, 2) == 1 "The input data must be monoral."
    x = map(Float64, vec(x))
    fs = Int(fs)
    @info("length of input signal is $(length(x)/fs) sec.")

    period = parse(Float64, args["--period"])
    order = parse(Int, args["--order"])
    α = parse(Float64, args["--alpha"])
    if α == 0.0
        α = mcepalpha(fs)
    end
    trajectory = args["--trajectory"]

    # Load mapping model
    gmm = load(args["<model_jld>"])
    if !gmm["diff"]
        @warn("The model doesn't seem to be trained on differencial features")
    end

    mapper = GMMMap(gmm["weights"], gmm["means"], gmm["covars"])
    if trajectory
        mapper = TrajectoryGMMMap(mapper, int(args["--T"]))
    end

    # shape (order+1, number of frames)
    elapsed_fe = @elapsed begin
        f0, timeaxis = dio(x, fs, DioOption(period=period))
        f0 = stonemask(x, fs, timeaxis, f0)
        spectrogram = cheaptrick(x, fs, timeaxis, f0)
        src = sp2mc(spectrogram, order, α)
    end
    @info("elapsed time in feature extraction is $(elapsed_fe) sec.")
    if trajectory
        # add delta feature
        src = [src[1,:]; push_delta(src[2:end,:])]
    end

    # remove power coef.
    src[1,:] = 0.0

    # Perform conversion
    elapsed_vc = @elapsed converted = vc(mapper, src)
    @info("elapsed time in conversion process is $(elapsed_vc) sec.")

    # Waveform synthesis using Mel-Log Spectrum Approximation filter
    mf = MLSADF(order, α)
    hopsize = round(Int, (fs / (1000 / period)))

    elapsed_syn = @elapsed y = synthesis!(mf, x, mc2b(converted, α), hopsize)
    @info("elapsed time in waveform moduration is $(elapsed_syn) sec.")

    wavwrite(map(Float64, y), args["<dst_wav>"], Fs=fs)
    @info("Dumped to ", args["<dst_wav>"])
end

main()


text:
module Example #hide
# # A first demo
# This example serves as a test since this project doesn't have a "testing" procedure yet. In this example,
# the equation ``u'(x) = 2`` with ``u(0) = 0`` is solved on the domain ``[0,1]`` using a backward finite
# difference scheme.
#
# In this example, PETSc classic method names are used. For more fancy names, check the fancy version.
#
# Note that the way we achieve things in the document can be highly improved and the purpose of this example
# is only demonstrate some method calls to give an overview.
#
# To run this example, execute : `mpirun -n your_favorite_positive_integer julia example1.jl`

# Import package
using PetscWrap

# Initialize PETSc. Command line arguments passed to Julia are parsed by PETSc. Alternatively, you can
# also provide "command line arguments by defining them in a string, for instance
# `PetscInitialize("-ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always")` or by providing each argument in
# separate strings : `PetscInitialize(["-ksp_monitor_short", "-ksp_gmres_cgs_refinement_type", "refine_always")`
PetscInitialize()

# Number of mesh points and mesh step
n = 11
Δx = 1. / (n - 1)

# Create a matrix and a vector (you can specify the MPI communicator if you want)
A = MatCreate()
b = VecCreate()

# Set the size of the different objects, leaving PETSC to decide how to distribute. Note that we should
# set the number of preallocated non-zeros to increase performance.
MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)
VecSetSizes(b, PETSC_DECIDE, n)

# We can then use command-line options to set our matrix/vectors.
MatSetFromOptions(A)
VecSetFromOptions(b)

# Finish the set up
MatSetUp(A)
VecSetUp(b)

# Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor.
# As in PETSc, the `rstart, rend = *GetOwnershipRange` methods indicate the first row handled by the local processor
# (starting at 0), and the last row (which is `rend-1`). This may be disturbing for non-regular PETSc users. Checkout
# the fancy version of this example for a more Julia-like convention.
b_start, b_end = VecGetOwnershipRange(b)

# Now let's build the right hand side vector. Their are various ways to do this, this is just one.
n_loc = VecGetLocalSize(b) # Note that n_loc = b_end - b_start...
VecSetValues(b, collect(b_start:b_end-1), 2 * ones(n_loc))

# And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices.
A_start, A_end = MatGetOwnershipRange(A)
for i in A_start:A_end-1
    MatSetValues(A, [i], [i-1, i], [-1. 1.] / Δx, INSERT_VALUES) # MatSetValues(A, I, J, V, INSERT_VALUES)
end

# Set boundary condition (only the proc handling index `0` is acting)
(b_start == 0) && VecSetValue(b, 0, 0.)

# Assemble matrices
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)
VecAssemblyBegin(b)
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)
VecAssemblyEnd(b)

# At this point, you can inspect `A` or `b` using a viewer (stdout by default), simply call
MatView(A)
VecView(b)

# Set up the linear solver
ksp = KSPCreate()
KSPSetOperators(ksp, A, A)
KSPSetFromOptions(ksp)
KSPSetUp(ksp)

# Solve the system. We first allocate the solution using the `VecDuplicate` method.
x = VecDuplicate(b)
KSPSolve(ksp, b, x)

# Print the solution
VecView(x)

# Access the solution (this part is under development), getting a Julia array; and then restore it
array, ref = VecGetArray(x) # do something with array
@show array
VecRestoreArray(x, ref)

# Free memory
MatDestroy(A)
VecDestroy(b)
VecDestroy(x)

# Finalize Petsc
PetscFinalize()

end #hide

max_stars_repo_name: bmxam/JuliaPetsc

max_issues_repo_name: bmxam/JuliaPetsc

max_forks_repo_name: bmxam/JuliaPetsc


text:
# this algorithms assumes that the initial-value problem is one-dimensional
# FIXME requires special scalar quadratic system type
function reach_homog_QINT(alg::QINT{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
              Δt0::TimeInterval=zeroI, kwargs...) where {N}

    n = statedim(ivp)
    n == 1 || throw(ArgumentError("this algorithm applies to one-dimensional " *
                    "systems, but this initial-value problem is $n-dimensional"))

    @unpack Δ, δ, θ, maxiter, approx_model = alg

    # normalize system to canonical form FIXME
    ivp_norm = _normalize(ivp)
    Ω0 = initial_state(ivp_norm)

    # discretize system
    #ivp_discr = discretize(ivp_norm, δ, approx_model)
    #Φ = state_matrix(ivp_discr)[1, 1]
    #Ω0 = initial_state(ivp_discr)
    #X = stateset(ivp_discr)

    # this algorithm requires Ω0 to be an interval
    #Ω0 = overapproximate(Ω0, Interval)

    # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
    got_homogeneous = !hasinput(ivp_discr)

    # preallocate output flowpipe
    #IT = typeof(Ω0)
    #N = eltype(Ω0)
    #F = Vector{ReachSet{N, IT}}(undef, NSTEPS)

    if got_homogeneous
        # FIXME consider Δt0
        F = reach_homog_QINT(a=a, b=b, c=c, X0=Ω0, T=T, Δ=Δ, δ=δ, θ=θ, maxiter=maxiter)
    else
        error("not implemented yet")
        #U = inputset(ivp_discr)
        #@assert isa(U, LazySet)
        #U = overapproximate(U, Interval) # TODO guarantee this on the discretization?
        #reach_inhomog_INT!(F, Ω0, Φ, NSTEPS, δ, X, U, Δt0)
    end

    return F
end


text:
#!/usr/bin/env julia

include("../../../../misc/StochasticArithmetic.jl/src/StochasticArithmetic.jl")
using Main.StochasticArithmetic
using Random

function twoSum(a, b)
  x = a + b
  z = x - a
  y = (a - (x-z)) + (b-z)
  return (x,y)
end

function split(a)
  # New methods should be added to compute the correct
  # constant for each FP representation
  #
  # factor = 1+2^s, where s is the mantissa length
  factor(x::Float64) = 134217729
  factor(x::Float32) = 4097

  c = factor(a) * a
  x = c - (c - a)
  y = a - x
  return (x, y)
end

function twoProd(a, b)
  x = a * b
  (a1, a2) = split(a)
  (b1, b2) = split(b)
  y = a2*b2 - (((x - a1*b1) - a2*b1) - a1*b2)
  return (x,y)
end

function twoDiv(a, b)
  x = a / b
  (p1, p2) = twoProd(x, b)
  y = (a-p1-p2)/b
  return (x, y)
end

function round_dir(op, a, b, mode)
  ops = Dict(:+ => twoSum,
             :* => twoProd,
             :/ => twoDiv)
  mode_pred   = Dict(:up   => (err -> (err>0)),
                     :down => (err -> (err<0)))
  mode_change = Dict(:up   => nextfloat,
                     :down => prevfloat)
  (res,err) = ops[op](a,b)
  if (err != 0 && mode_pred[mode](err))
      res = mode_change[mode](res)
  end

  return res
end

round_rnd = let rng = MersenneTwister()
  function(op, a, b)
      return round_dir(op, a, b, rand(rng, [:up, :down]))
  end
end

my_prod(x, y) = x * y
my_prod(x::Float64, y::Float64) = round_rnd(:*, x, y)
my_prod(x::Float32, y::Float32) = round_rnd(:*, x, y)

my_sum(x, y) = x + y
my_sum(x::Float64, y::Float64) = round_rnd(:+, x, y)
my_sum(x::Float32, y::Float32) = round_rnd(:+, x, y)

macro mca(defun)
    transform!(x) = Nothing()
    transform!(e::Expr) = begin
        for i in 1:length(e.args)
            if typeof(e.args[i]) == Expr && e.args[i].head == :(+=)
                ee = e.args[i]
                transform!.(ee.args)
                e.args[i] = :($(ee.args[1]) = $(ee.args[1]) + $(ee.args[2]))
            elseif typeof(e.args[i]) == Expr && e.args[i].head == :(*=)
                ee = e.args[i]
                transform!.(ee.args)
                e.args[i] = :($(ee.args[1]) = $(ee.args[1]) * $(ee.args[2]))
            elseif e.args[i] == :*
                e.args[i] = :my_prod
            elseif e.args[i] == :+
                e.args[i] = :my_sum
            end
            transform!(e.args[i])
        end
    end

    transform!(defun)
    println(defun)
    defun
end

struct Params{T}
  dt::T
  dx::T
  v::Vector{T}
end

function _unpack(p::Params{T}) where{T}
  return (p.dt,p.dx)
end
function export_data(filename::String,X::Vector{T},Y::Vector{U}) where{T,U}
  open(filename,"w") do f
    for i in 1:size(X,1)
      x = X[i]
      y = Y[i]
      write(f,"$x $y\n")
    end
  end
end
function export_data_tuple(filename::String,X::Vector{T},Y::Vector{U}) where{T,U}
  open(filename,"w") do f
    for i in 1:size(X,1)
      x = X[i]
      y1,y2 = Y[i]
      write(f,"$x $y1 $y2\n")
    end
  end
end

@mca function euler(un::Vector{T},L,tn::T,param::Params{T}) where{T}
  un1 = Vector{T}(undef,size(un))
  un1 = un - param.dt*L(un,tn,param)
  return un1
end

@mca function rk33(un::Vector{T},L,tn::T,param::Params{T}) where{T}
  dt = param.dt
  u1 = Vector{T}(undef,size(un))
  u2 = Vector{T}(undef,size(un))
  un1 = Vector{T}(undef,size(un))

  Lun = L(un,tn,param)
  for i in 1:size(un,1)
    u1[i] = un[i] - dt*Lun[i]
  end
  Lu1 = L(u1,tn+dt,param)
  for i in 1:size(un,1)
    u2[i] = 0.75*un[i] + 0.25*u1[i] - 0.25*(dt*Lu1[i])
  end
  Lu2 = L(u2,tn+0.5*dt,param)
  for i in 1:size(un,1)
    un1[i] = (1.0/3.0)*un[i] + (2.0/3.0)*u2[i] - (2.0/3.0)*(dt*Lu2[i])
  end
  return un1
end

@mca function cd2(un::Vector{T},t::T,param::Params{T}) where{T}
  dxu = Vector{T}(undef,size(un))

  dxu[1] = un[1+1]-un[size(un,1)]
  dxu[size(un,1)] = un[1]-un[size(un,1)-1]
  for i in 2:size(un,1)-1
    dxu[i] = un[i+1]-un[i-1]
  end

  return dxu/(2.0*param.dx)
end

@mca function flux(u::Vector{T},im2::Int,im1::Int,i::Int,ip1::Int,ip2::Int,ip3::Int) where{T}
  eps = 1e-6

  b0p = (13.0/12.0)*( u[im2] - 2.0*u[im1] + u[i]   )^2 + 0.25*( u[im2] - 4.0*u[im1] + 3.0*u[i] )^2
  b1p = (13.0/12.0)*( u[im1] - 2.0*u[i]   + u[ip1] )^2 + 0.25*( u[im1] - u[ip1] )^2
  b2p = (13.0/12.0)*( u[i]   - 2.0*u[ip1] + u[ip2] )^2 + 0.25*( 3.0*u[i] - 4.0*u[ip1] + u[ip2] )^2

  b0m = (13.0/12.0)*( u[ip1] - 2.0*u[ip2] + u[ip3] )^2 + 0.25*( 3.0*u[ip1] - 4.0*u[ip2] + u[ip3] )^2
  b1m = (13.0/12.0)*( u[i]   - 2.0*u[ip1] + u[ip2] )^2 + 0.25*( u[i] - u[ip2] )^2
  b2m = (13.0/12.0)*( u[im1] - 2.0*u[i]   + u[ip1] )^2 + 0.25*( u[im1] - 4.0*u[i] + 3.0*u[ip1] )^2

  a0p = 0.1/(eps+b0p)^2
  a1p = 0.6/(eps+b1p)^2
  a2p = 0.3/(eps+b2p)^2

  a0m = 0.1/(eps+b0m)^2
  a1m = 0.6/(eps+b1m)^2
  a2m = 0.3/(eps+b2m)^2

  sum_ap = a0p+(a1p+a2p)
  sum_am = a0m+(a1m+a2m)

  w0p = a0p/sum_ap
  w1p = a1p/sum_ap
  w2p = a2p/sum_ap

  w0m = a0m/sum_am
  w1m = a1m/sum_am
  w2m = a2m/sum_am

  uip12p = (w0p*( (2.0/6.0)*u[im2] + ((-7.0/6.0)*u[im1] + (11.0/6.0)*u[i]) ) + w1p*( (-1.0/6.0)*u[im1] + ((5.0/6.0)*u[i] + (2.0/6.0)*u[ip1]) )) + w2p*( (2.0/6.0)*u[i] + ((5.0/6.0)*u[ip1] + (-1.0/6.0)*u[ip2]) )
  uip12m = (w2m*( (-1.0/6.0)*u[im1] + ((5.0/6.0)*u[i] + (2.0/6.0)*u[ip1]) ) + w1m*( (2.0/6.0)*u[i] + ((5.0/6.0)*u[ip1] + (-1.0/6.0)*u[ip2]) )) + w0m*( (11.0/6.0)*u[ip1] + ((-7.0/6.0)*u[ip2] + (2.0/6.0)*u[ip3]) )

  return (uip12p,uip12m)
end

@mca function weno(un::Vector{T},t::T,param::Params{T}) where{T}
  dxu = Vector{T}(undef,size(un))
  vm = Vector{T}(undef,size(param.v))
  vp = Vector{T}(undef,size(param.v))
  for i in 1:size(param.v,1)
    vp[i] = max(param.v[i],0.0)
    vm[i] = min(param.v[i],0.0)
  end

  N = size(un,1)
  i=-1
  uim12p,uim12m = flux(un, ((i-2)+N)%N+1 , ((i-1)+N)%N+1 , (i+N)%N+1 , ((i+1)+N)%N+1 , (i+2)+1 , (i+3)+1 )
  for i in 0:2
    uip12p,uip12m = flux(un, ((i-2)+N)%N+1 , ((i-1)+N)%N+1 , (i+N)%N+1 , ((i+1)+N)%N+1 , (i+2)+1 , (i+3)+1 )
    dxu[i+1] = vp[i+1]*(uip12p-uim12p) + vm[i+1]*(uip12m-uim12m)
    uim12p = uip12p
    uim12m = uip12m
  end
  for i in 3:size(un,1)-4
    uip12p,uip12m = flux(un, (i-2)+1 , (i-1)+1 , (i)+1 , (i+1)+1 , (i+2)+1 , (i+3)+1 )
    dxu[i+1] = vp[i+1]*(uip12p-uim12p) + vm[i+1]*(uip12m-uim12m)
    uim12p = uip12p
    uim12m = uip12m
  end
  for i in size(un,1)-3:size(un,1)-1
    uip12p,uip12m = flux(un, (i-2)+1 , (i-1)+1 , (i)+1 , (i+1)%N+1 , (i+2)%N+1 , (i+3)%N+1 )
    dxu[i+1] = vp[i+1]*(uip12p-uim12p) + vm[i+1]*(uip12m-uim12m)
    uim12p = uip12p
    uim12m = uip12m
  end

  return dxu/param.dx
end

@mca function simu(time_method,space_method,Tf::T,u0::Vector{T},p::Params{T}) where{T}
  current_time = T(0.0)
  n_iter = 0
  un = Vector{T}(u0)
  while ( n_iter*p.dt < Tf )
    un = time_method(un,space_method,current_time,p)
    current_time += p.dt
    n_iter += 1
  end
  return un
end

function launcher(label::String,func,X::Vector{T},Tf::T,p::Params{T}) where{T}
  print("\033[94m>\033[0m $label\n")
  u0 = func(X)
  export_data("results/init_$label.dat",X,u0)

  print("  \033[95m+\033[0m cd2\n")
  un = @reliable_digits simu(rk33,cd2,Tf,Float64.(u0),p)
  export_data_tuple("results/end_cd2_$label.dat",X,un)

  print("  \033[95m+\033[0m weno\n")
  un = @reliable_digits simu(rk33,weno,Tf,Float64.(u0),p)
  export_data_tuple("results/end_weno_$label.dat",X,un)
end

function gate(X::Vector{T}) where{T}
  u = Vector{T}(undef,size(X))
  i = 1
  for x in X
    if x < 0.25 || x > 0.75
      u[i] = 0.0
    else
      u[i] = 1.0
    end
    i += 1
  end
  return u
end

function hat(X::Vector{T}) where{T}
  u = Vector{T}(undef,size(X))
  i = 1
  for x in X
    if x < 0.25 || x > 0.75
      u[i] = 0.0
    elseif x<0.5
      u[i] = 4*x-1
    else
      u[i] = -4*x+3
    end
    i += 1
  end
  return u
end

function sinus(X::Vector{T}) where{T}
  return sin.((2.0*π)*X)
end

function cosinus(X::Vector{T}) where{T}
  return cos.((2.0*π)*X)
end

function sinusHF(X::Vector{T}) where{T}
  return sin.(6*(2.0*π)*X)
end

function main()
  N = 100
  X = Vector{Float64}(Vector([0:N-1])[1]/N)
  dx = Float64(1.0/N)
  v = Vector{Float64}(fill(1.0,N))
  p = Params{Float64}(0.1*dx,dx,v)
  Tf = 1.0

  #--- cos -------------------------------
  launcher("cos",cosinus,X,Tf,p)
  #--- sin -------------------------------
  launcher("sin",sinus,X,Tf,p)
  #--- sinHF -------------------------------
  launcher("sinHF",sinusHF,X,Tf,p)
  #--- hat -------------------------------
  launcher("hat",hat,X,Tf,p)
  #--- gat -------------------------------
  launcher("gat",gate,X,Tf,p)
end

main()

max_stars_repo_name: kivvix/draft

max_issues_repo_name: kivvix/draft

max_forks_repo_name: kivvix/draft


text:
import Pkg
if VERSION < v"1.4"
   if get(Pkg.installed(), "Tables", v"0.0") < v"1.0"
      include("./dummy_tables.jl")
      @eval import .Tables
   else
      @eval import Tables
   end
else
   # no check for correct version of Tables in Julia 1.4 since
   # Pkg.dependencies() is "experimental" and Pkg.installed is deprecated
   try
      @eval import Tables
   catch ex
      include("./dummy_tables.jl")
      @eval import .Tables
   end
end

include("./RConnector.jl")
RConnector.serve(parse(Int, ARGS[1]), portfile = ARGS[2],
   multiclient = (ARGS[3] == "t"))


text:
export SolitaryWaveWhithamBoussinesq

"""
    `SolitaryWaveWhithamBoussinesq(param; kwargs...)`

Computes the Whitham-Boussinesq solitary wave with prescribed velocity.

# Argument
- `param :: NamedTuple`: parameters of the problem containing velocity `c` and dimensionless parameters `ϵ` and `μ`, and mesh size `L` and number of collocation points `N`;

## Keywords (optional)
- `guess :: Vector{Real}`: initial guess for the surface deformation (if not provided, the exact formula for SGN is used);
- `x₀ :: Real`: center of solitary wave (if guess is not provided);
- `model :: Real`: determines the model used (typically `1` or `1/2`, default is 1);
- `iterative :: Bool`: inverts Jacobian through GMRES if `true`, LU decomposition if `false`;
- `verbose :: Bool`: prints numerical errors at each step if `true`;
- `max_iter :: Int`: maximum number of iterations of the Newton algorithm;
- `tol :: Real`: general tolerance (default is `1e-10`);
- `ktol :: Real`: tolerance of the Krasny filter (default is `0`, i.e. no filtering);
- `gtol :: Real`: relative tolerance of the GMRES algorithm;
- `dealias :: Int`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `q :: Real`: Newton algorithm modified with
`u_{n+1}=q*u_{n+1}+(1-q)*u_n`
(default is `1`);
- `α :: Real`: adds `α` times spectral projection onto the Kernel to the Jacobian.
# Return values
`(η,u) :: Tuple{Vector{Float64},Vector{Float64}}`
with
- `η`: surface deformation;
- `u`: velocity.
"""
function SolitaryWaveWhithamBoussinesq(
        param :: NamedTuple;
        guess = zeros(0) :: Vector{Float64},
        x₀ = 0 :: Real,
        model = 1 :: Real,
        iterative = false :: Bool,
        verbose = true :: Bool,
        max_iter = 20 :: Int,
        tol = 1e-10 :: Real,
        ktol = 0 :: Real,
        gtol = 1e-10 :: Real,
        dealias = 0 :: Int,
        q=1 :: Real,
        α=0 :: Real)
        # A good guess for low velocities is
        #function sol(x,α)
        #	2*α*sech.(sqrt(3*2*α)/2*x).^2
        #end
        c = param.c
        ϵ = param.ϵ
        μ = param.μ

        mesh=Mesh(param)
        if guess == zeros(0)
                @info "Using the exact formula for the SGN solitary wave as initial guess"
                η = (c^2-1)/ϵ*sech.(sqrt(3*(c^2-1)/(c^2)/μ)/2*(mesh.x.-x₀)).^2
                h = 1 .+ param.ϵ*η
                u = c*η./h
                k = mesh.k
        	F₀ = sqrt(param.μ)*1im * k
        	DxF(v) = real.(ifft(F₀ .* fft(v)))
        	guess = u - 1/3 ./h .* (DxF(h.^3 .*DxF(u)))
        end

        k = mesh.k
        F₁ 	= tanh.(sqrt(μ)*abs.(k))./(sqrt(μ)*abs.(k))
        F₁[1] 	= 1
        F₂ = F₁.^(model)
        Dx       =  1im * sqrt(μ)*k

        if dealias == 0
                Π = ones(size(k))
        else
                Π = abs.(k) .< maximum(k) * (1-dealias/(dealias+2))
        end
        krasny(k) = (abs.(k).>= ktol ).*k
        krasny!(k) = k[abs.(k).< ktol ].=0

        function filt( v :: Vector{Float64} )
                real.(ifft(krasny(Π.*fft(v))))
        end
        function filt( v :: Vector{Complex{Float64}} )
                ifft(krasny(Π.*fft(v)))
        end
        function Four1( v  )
                ifft(F₁.*fft(v))
        end
        function Four2( v  )
                ifft(F₂.*fft(v))
        end
        function FF( v )
                w = Four2(v)
                return real.(-c^2*v .+ Four1(v) .+ c*ϵ/2*w.^2 .+ ϵ*Four2(c*w.*v .- ϵ/2*w.^3))
        end
        function Fabs( v )
                w = Four2(v)
                return real.(c^2*abs.(v) .+ abs.(Four1(v)) .+ c*ϵ/2*w.^2 .+ ϵ*Four2(c*abs.(w.*v) .+ ϵ/2*abs.(w).^3))
        end

        if iterative == false
                k = mesh.k
                x = mesh.x
                x₀ = mesh.x[1]
                FFT = exp.(-1im*k*(x.-x₀)');
                IFFT = exp.(1im*k*(x.-x₀)')/length(x);
                Id = Diagonal( ones(size(x)));
                M₁ = IFFT*(Diagonal( F₁)*FFT)
                M₂ = IFFT*(Diagonal( F₂)*FFT)
                function JacF( v ,dxv )
                        M₀ = Diagonal( Four2(v))
                        Symmetric(real.(-c^2 .* Id .+ M₁ .+ c*ϵ.*(M₀*M₂ .+ M₂*M₀) .+ ϵ*M₂*(Diagonal(c*v -3*ϵ/2*Four2(v).^2))*M₂ .+ α*dxv*dxv'))
                end
        else
                function JacFfast( v , dxv )
                        dF(φ) = real.(-c^2*φ+Four1(φ)+c*ϵ*Four2(v).*Four2(φ)+ϵ*Four2(c*Four2(v).*φ+(c*v -3*ϵ/2*(Four2(v).^2)).*Four2(φ))+ α*dot(dxv,φ)*dxv)
                        return LinearMap(dF, length(v); issymmetric=true, ismutating=false)
                end
        end

        flag=0
        iter = 0
        u = filt(guess)
        du = similar(u)
        fu = similar(u)
        dxu = similar(u)
        for i in range(0, stop=max_iter)
                dxu .= real.(ifft(Dx.*fft(u)))
                dxu ./= norm(dxu,2)
                fu .= FF(u)
                relerr = norm(fu,Inf)/norm(Fabs(u),Inf)
                abserr = norm(fu,Inf)
                if relerr < tol
    			@info string("Converged : relative error ",relerr," in ",i," steps\n")
    			break
    		elseif verbose == true
                        print(string("absolute error at step ",i,": ",abserr,"\n"))
                        print(string("relative error at step ",i,": ",relerr,"\n"))

    		end
                if i == max_iter
                        @warn string("The algorithm did not converge after ",i," steps: final relative error is ",relerr,"\n")
                        break
                end

                if iterative == false
                        du .= JacF(u,dxu) \ fu
                else
                        du .= gmres( JacFfast(u,dxu) , fu ; verbose = verbose, tol = gtol )
                end
    		u .-= q*filt(du)
        end
        return (c*u-ϵ/2*real.(Four2(u).^2),u,mesh)
end


text:
using DataStructures
using Mmap

ismissingtype(T) = Missing <: T
ismissingeltype(T) = missingtype(eltype(T))

const UnionMissing{T} = Union{Missing, T}

export csvread
const current_record = Ref{Any}()
const debug = Ref{Bool}(false)

const StringLike = Union{AbstractString, StrRange}

optionsiter(opts::AbstractDict) = opts
optionsiter(opts::AbstractVector) = enumerate(opts)

getbyheader(opts, header, i::Int) = opts[i]
getbyheader(opts, header, i::Symbol) = getbyheader(opts, header, string(i))
function getbyheader(opts, header, i::AbstractString)
    if !(i in header)
        throw(ArgumentError("Unknown column $i"))
    end
    getbyheader(opts, header, something(findfirst(isequal(i), header), length(header)+1))
end

function optionsiter(opts::AbstractDict, header)
    isempty(header) && return opts
    iter = Dict{Int,Any}()
    for (k, v) in opts
        i = try
            getbyheader(1:length(header), header, k)
        catch err
            if isa(err, ArgumentError)
                continue
            else
                rethrow(err)
            end
        end
        iter[i] = v
    end
    iter
end

optionsiter(opts::AbstractVector, header) = optionsiter(opts)

tofield(f::AbstractField, opts) = f
tofield(f::AbstractToken, opts) = Field(f)
tofield(f::StringToken, opts) = Field(Quoted(f))
tofield(f::Type, opts) = tofield(fromtype(f), opts)
tofield(f::Type{String}, opts) = tofield(fromtype(StrRange), opts)
tofield(f::DateFormat, opts) = tofield(DateTimeToken(DateTime, f), opts)

"""
    csvread(file::Union{String,IO}, delim=','; <arguments>...)

Read CSV from `file`. Returns a tuple of 2 elements:
1. A tuple of columns each either a `Vector`, or `StringArray`
2. column names if `header_exists=true`, empty array otherwise

# Arguments:

- `file`: either an IO object or file name string
- `delim`: the delimiter character
- `spacedelim`: (Bool) parse space-delimited files. `delim` has no effect if true.
- `quotechar`: character used to quote strings, defaults to `"`
- `escapechar`: character used to escape quotechar in strings. (could be the same as quotechar)
- `nrows`: number of rows in the file. Defaults to `0` in which case we try to estimate this.
- `skiplines_begin`: skips specified number of lines at the beginning of the file
- `header_exists`: boolean specifying whether CSV file contains a header
- `nastrings`: strings that are to be considered NA. Defaults to `TextParse.NA_STRINGS`
- `colnames`: manually specified column names. Could be a vector or a dictionary from Int index (the column) to String column name.
- `colparsers`: Parsers to use for specified columns. This can be a vector or a dictionary from column name / column index (Int) to a "parser". The simplest parser is a type such as Int, Float64. It can also be a `dateformat"..."`, see [CustomParser](@ref) if you want to plug in custom parsing behavior
- `type_detect_rows`: number of rows to use to infer the initial `colparsers` defaults to 20.
"""
csvread(file::String, delim=','; kwargs...) = _csvread_f(file, delim; kwargs...)[1:2]

function csvread(file::IOStream, delim=','; kwargs...)
    mmap_data = Mmap.mmap(file)
    _csvread(String(pointer(mmap_data), length(mmap_data)), delim; kwargs...)
end

function csvread(buffer::IO, delim=','; kwargs...)
    mmap_data = read(buffer)
    _csvread(String(mmap_data), delim; kwargs...)
end

function _csvread(str::AbstractString, delim=','; kwargs...)
    _csvread_internal(str, delim; kwargs...)[1:2]
end

function _csvread_f(file::AbstractString, delim=','; kwargs...)
    # Try to detect file extension for compressed files
    ext = last(split(file, '.'))

    if ext == "gz" # Gzipped
        return open(GzipDecompressorStream, file, "r") do io
            data = read(io)
            _csvread_internal(String(data), delim; filename=file, kwargs...)
        end
    else # Otherwise just try to read the file
        return open(file, "r") do io
            data = Mmap.mmap(io)
            _csvread_internal(String(data), delim; filename=file, kwargs...)
        end
    end
end

const ColsPool = OrderedDict{Union{Int, String}, AbstractVector}

function csvread(files::AbstractVector{T},
                 delim=','; kwargs...) where {T<:AbstractString}
    @assert !isempty(files)
    colspool = ColsPool()
    cols, headers, parsers, nrows = try
        _csvread_f(files[1], delim;
                   noresize=true,
                   colspool=colspool,
                   kwargs...)
    catch err
        println(stderr, "Error parsing $(files[1])")
        rethrow(err)
    end

    count = Int[nrows]
    prev = nrows
    for f in files[2:end]
        if !isempty(cols) && length(cols[1]) == nrows
            n = ceil(Int, nrows * sqrt(2))
            resizecols(colspool, n)
        end
        cols, headers, parsers, nrows = try
            _csvread_f(f, delim; rowno=nrows+1, colspool=colspool,
                       prevheaders=headers, noresize=true, prev_parsers=parsers, kwargs...)
        catch err
            println(stderr, "Error parsing $(f)")
            rethrow(err)
        end
        push!(count, nrows - prev)
        prev = nrows
    end

    resizecols(colspool, nrows)
    (values(colspool)...,), collect(keys(colspool)), count
end

# read CSV in a string
function _csvread_internal(str::AbstractString, delim=',';
                 spacedelim=false,
                 quotechar='"',
                 escapechar='"',
                 stringtype=String,
                 noresize=false,
                 rowno::Int=1,
                 prevheaders=nothing,
                 pooledstrings=nothing,
                 skiplines_begin=0,
                 samecols=nothing,
                 header_exists=true,
                 nastrings=NA_STRINGS,
                 colnames=String[],
                 #ignore_empty_rows=true,
                 colspool = ColsPool(),
                 nrows = !isempty(colspool) ?
                     length(first(colspool)[2]) : 0,
                 prev_parsers = nothing,
                 colparsers=[],
                 filename=nothing,
                 type_detect_rows=20)

    if pooledstrings === true
        warn("pooledstrings argument has been removed")
    end
    opts = LocalOpts(delim, spacedelim, quotechar, escapechar, false, false)
    len = lastindex(str)
    pos = firstindex(str)
    rowlength_sum = 0   # sum of lengths of rows, for estimating nrows
    lineno = 0

    y = iterate(str, pos)
    if y!==nothing
        c = y[1]; i = y[2]
        if c == '\ufeff'
            pos = i
        end
    end

    pos, lines = eatnewlines(str, pos)
    lineno += lines
    while lineno < skiplines_begin
        pos = getlineend(str, pos)
        y2 = iterate(str, pos)
        y2===nothing && error("Internal error.")
        pos = y2[2]
        pos, lines = eatnewlines(str, pos)
        lineno += lines
    end
    if header_exists
        merged_colnames, pos = readcolnames(str, opts, pos, colnames)
        lineno += 1
    else
        merged_colnames = colnames
    end
    merged_colnames = map(string, merged_colnames)

    if !issorted(nastrings)
        nastrings = sort(nastrings)
    end

    pos1 = pos

    if samecols === nothing
        canonnames = merged_colnames
    else
        canonnames = map(merged_colnames) do c
            canonical_name(samecols, c)
        end
    end

    if length(unique(canonnames)) != length(canonnames)
        error("""File has multiple column headers with the same name, specify `colnames` by hand
                 along with `header_exists=false`""")
    end

    # seed guesses using those from previous file
    guess, pos1 = guesscolparsers(str, canonnames, opts,
                                  pos, type_detect_rows, colparsers,
                                  nastrings, prev_parsers)
    if isempty(canonnames)
        canonnames = Any[1:length(guess);]
    end

    for (i, v) in enumerate(guess)
        c = get(canonnames, i, i)
        # Make column nullable if it's showing up for the
        # first time, but not in the first file
        if !(fieldtype(v) <: StringLike) && prev_parsers !== nothing && !haskey(colspool, c)
            v = isa(v, NAToken) ? v : NAToken(v)
        end
        p = tofield(v, opts)
        guess[i] = p
    end

    # the last field is delimited by line end
    if !isempty(guess)
        guess[end] = swapinner(guess[end], guess[end]; eoldelim = true)
        rec = Record((guess...,))
    else
        parsers = prev_parsers === nothing ? Dict() : copy(prev_parsers)
        rec = Record(())
        return (), String[], parsers, rowno-1
    end

    if isempty(canonnames)
        canonnames = Any[1:length(rec.fields);]
    end

    current_record[] = rec

    if nrows == 0
        # just an estimate, with some margin
        nrows = ceil(Int, pos1/max(1, lineno) * sqrt(2))
    end

    if isempty(colspool)
        # this is the first file, use nrows
        cols = makeoutputvecs(rec, nrows, stringtype)
        for (c2, h) in zip(cols, canonnames)
            colspool[h] = c2
        end
    else
        _cols = map(1:length(rec.fields)) do i
            c = get(canonnames, i, i)
            f = rec.fields[i]
            if haskey(colspool, c)
                if eltype(colspool[c]) == fieldtype(f) || (fieldtype(f) <: StrRange && eltype(colspool[c]) <: AbstractString)
                    return colspool[c]
                else
                    try
                        return colspool[c] = promote_column(colspool[c],
                                                            rowno-1,
                                                            fieldtype(f), stringtype)
                    catch err
                        error("Could not convert column $c of eltype $(eltype(colspool[c])) to eltype $(fieldtype(f))")
                    end
                end
            else
                return colspool[c] = makeoutputvec(f, nrows, stringtype)
            end
        end
        # promote missing columns to nullable
        missingcols = setdiff(collect(keys(colspool)), canonnames)
        for k in missingcols
            if !ismissingtype(eltype(colspool[k])) && !(eltype(colspool[k]) <: StringLike)
                colspool[k] = promote_column(colspool[k],
                                             rowno-1,
                                             UnionMissing{eltype(colspool[k])}, stringtype)
            end
        end
        cols = (_cols...,)
    end

    if any(c->length(c) != nrows, cols)
        resizecols(colspool, nrows)
    end

    finalrows = rowno
    @label retry
    try
        finalrows = parsefill!(str, opts, rec, nrows, cols, colspool,
                               pos, lineno, rowno, lastindex(str))
        if !noresize
            resizecols(colspool, finalrows)
        end
    catch err

        if !isa(err, CSVParseError)
            rethrow(err)
        end

        err.filename = filename

        if err.err_code == PARSE_ERROR

            rng = getlineat(str, err.fieldpos)
            f, l = first(rng), last(rng)
            field = rec.fields[err.colno]

            if l !== lastindex(str) && err.pos >= l && !field.eoldelim
                if fieldtype(field) <: AbstractString || fieldtype(field) <: StrRange
                    # retry assuming newlines can be part of the field
                    wopts = LocalOpts(opts.endchar, opts.spacedelim, opts.quotechar, opts.escapechar, opts.includequotes, true)
                    fieldsvec = Any[rec.fields...]
                    fieldsvec[err.colno] = swapinner(field, WrapLocalOpts(wopts, field.inner))
                    rec = Record((fieldsvec...,))
                    pos = first(rng)
                    rowno = err.rowno
                    lineno = err.lineno
                    @goto retry
                end
                println(stderr, "Expected another field on row $(err.rowno) (line $(err.lineno))")
                err.filename = filename
                rethrow(err)
            end

            # figure out a new token type for this column and the rest
            # it's very likely that a number of columns change type in a single row
            # so we promote all columns after the failed column
            failed_strs = quotedsplit(str[err.fieldpos:l], opts, true)

            if length(failed_strs) != length(cols[err.colno:end])
                fn = err.filename === nothing ? "" : "In $(err.filename) "
                warn("$(fn)line $(err.lineno) has $(length(err.colno) + length(failed_strs) - 1) fields but $(length(cols)) fields are expected. Skipping row.")
                pos = last(rng)+1
                rowno = err.rowno
                lineno = err.lineno+1
                @goto retry
            end
            promoted = map(failed_strs, err.colno:length(cols)) do s, colidx
                col = cols[colidx]
                f = rec.fields[colidx]
                name = get(canonnames, colidx, colidx)
                c = promote_field(s, f, col, err, nastrings, stringtype)
                colspool[name] = c[2]
                c
            end

            newfields = map(first, promoted)
            newcols = map(last, promoted)

            if field.inner == newfields[1].inner
                println(stderr, "Could not determine which type to promote column to.")
                rethrow(err)
            end

            fieldsvec = Any[rec.fields...]
            fieldsvec[err.colno:end] = newfields
            typeof(cols)
            colsvec = Any[cols...]
            colsvec[err.colno:end] = newcols

            rec = Record((fieldsvec...,))
            cols = (colsvec...,)
            rowno = err.rowno
            lineno = err.lineno
            pos = first(rng)
            @goto retry

        end

    end

    parsers = prev_parsers === nothing ? Dict() : copy(prev_parsers)
    for i in 1:length(rec.fields)
        name = get(canonnames, i, i)
        parsers[name] = rec.fields[i].inner
    end
    cols, canonnames, parsers, finalrows
end

function promote_field(failed_str, field, col, err, nastrings, stringtype)
    newtoken = guesstoken(failed_str, field.inner, nastrings)
    if newtoken == field.inner
        # no need to change
        return field, col
    end
    newcol = try
        promote_column(col,  err.rowno-1, fieldtype(newtoken), stringtype)
    catch err2
        Base.showerror(stderr, err2)
        rethrow(err)
    end
    swapinner(field, newtoken), newcol
end

function promote_column(col, rowno, T, stringtype, inner=false)
    if typeof(col) <: Array{Missing}
        if T <: StringLike
            arr = StringVector{stringtype}(length(col))
            for i = 1:rowno
                arr[i] = ""
            end
            return arr
        elseif ismissingtype(T)
            fill!(Array{UnionMissing{eltype(T)}}(undef, length(col)), missing) # defaults to fill missing
        else
            error("empty to non-nullable")
        end
    elseif ismissingtype(T)
        arr = convert(Array{UnionMissing{eltype(col)}}, col)
        for i=rowno+1:length(arr)
            # if we convert an Array{Int} to be missing-friendly, we will not have missing in here by default
            arr[i] = missing
        end
        return arr
    else
        newcol = Array{T, 1}(undef, length(col))
        copyto!(newcol, 1, col, 1, rowno)
        newcol
    end
end

function readcolnames(str, opts, pos, colnames)
    colnames_inferred = String[]

    len = lastindex(str)
    lineend = getlineend(str, pos, len)
    head = str[pos:lineend]

    colnames_inferred = quotedsplit(str, opts, false, pos, lineend)
    # TODO: unescape

    # set a subset of column names
    for (i, v) in optionsiter(colnames, colnames_inferred)
        colnames_inferred[i] = v
    end
    colnames_inferred, lineend+1
end

function guesscolparsers(str::AbstractString, header, opts::LocalOpts, pos::Int,
                       nrows::Int, colparsers, nastrings=NA_STRINGS, prevs=nothing)
    # Field type guesses
    guess = []
    prevfields = String[]

    givenkeys = !isempty(colparsers) ? first.(collect(optionsiter(colparsers, header))) : []
    for i2=1:nrows
        pos, _ = eatnewlines(str, pos)
        if pos > lastindex(str)
            break
        end

        lineend = getlineend(str, pos)

        fields = quotedsplit(str, opts, true, pos, lineend)

        if i2 == 1
            guess = Any[Unknown() for i3=1:length(fields)] # idk
            if prevs !== nothing && !isempty(header)
                # sometimes length(fields) can be != length(header).
                # this sucks!
                for i4 in 1:length(header)
                    i4 > length(fields) && break
                    guess[i4] = get(prevs, header[i4], Unknown())
                end
            end
        end

        # update guess
        for j in 1:length(guess)
            if j in givenkeys
                continue # user specified this
            end
            if length(fields) != length(guess)
                error("previous rows had $(length(guess)) fields but row $i2 has $(length(fields))")
            end
            try
                guess[j] = guesstoken(fields[j], guess[j], nastrings)
            catch err
                println(stderr, "Error while guessing a common type for column $j")
                println(stderr, "new value: $(fields[j]), prev guess was: $(guess[j])")
                if j > 1
                    println(stderr, "prev value: $(fields[j-1])")
                end

                rethrow(err)
            end
        end
        prevfields = fields
        pos = lineend+1
    end

    # override guesses with user request
    for (i, v) in optionsiter(colparsers, header)
        guess[i] = tofield(v, opts)
    end
    guess, pos
end

function parsefill!(str::AbstractString, opts, rec::RecN{N}, nrecs, cols, colspool,
                    pos, lineno, rowno, l=lastindex(str)) where {N}
    pos, lines = eatnewlines(str, pos)
    lineno += lines
    pos <= l && while true
        prev_j = pos
        lineno += lines
        res = tryparsesetindex(rec, str, pos, l, cols, rowno, opts)
        if !issuccess(res)
            pos, fieldpos, colno, err_code = geterror(res)
            throw(CSVParseError(err_code, str, rec, lineno+1, rowno,
                                colno, pos, fieldpos))
        else
            pos = value(res)
        end

        pos, lines = eatnewlines(str, pos)
        lineno += lines

        if pos > l
            return rowno
        end
        rowno += 1
        lineno += 1
        if rowno > nrecs
            # grow
            nrecs = ceil(Int, rowno * sqrt(2)) # updated estimate
            resizecols(colspool, nrecs)
        end
    end
    return rowno # finished before starting
end

function resizecols(colspool, nrecs)
    for (h, c) in colspool
        l = length(c)
        resize!(c, nrecs)
        if eltype(c) <: AbstractString
            # fill with blanks
            c[l+1:nrecs] .= ""
        elseif eltype(c) <: StrRange
            c[l+1:nrecs] .= StrRange(1,0)
        end
    end
end

function makeoutputvecs(rec, N, stringtype)
    map(f->makeoutputvec(f, N, stringtype), rec.fields)
end

function makeoutputvec(eltyp, N, stringtype)
    if fieldtype(eltyp) == Missing # we weren't able to detect the type,
                                   # all cells were blank
        Array{Missing}(undef, N)
    elseif fieldtype(eltyp) == StrRange
        StringVector{stringtype}(N)
    elseif ismissingtype(fieldtype(eltyp)) && fieldtype(eltyp) <: StrRange
        StringVector{Union{Missing, String}}(N)
    else
        Array{fieldtype(eltyp)}(undef, N)
    end
end

mutable struct CSVParseError <: Exception
    err_code
    location_display
    rec
    lineno
    rowno
    colno
    pos
    fieldpos
    charinline
    filename
end

function CSVParseError(e_code, str, rec, lineno, rowno, colno, pos, fieldpos)
    rng = getlineat(str, pos)
    charinline = pos - first(rng)
    CSVParseError(e_code, showerrorchar(str, pos, 100), rec, lineno, rowno, colno, pos, fieldpos, charinline, nothing)
end

function Base.showerror(io::IO, err::CSVParseError)
    if err.filename !== nothing
        print(io, "CSV parsing error in $(err.filename) ")
    else
        print(io, "CSV parsing error ")
    end

    println(io, "at line $(err.lineno) char $(err.charinline):")
    println(io, err.location_display)
    print(io, "column $(err.colno) is expected to be: ")
    print(io, string(err.rec.fields[err.colno]))
end

function showerrorchar(str, pos, maxchar)
    hmaxchar = round(Int, maxchar/2)
    rng = getlineat(str, pos)
    substr = strip(str[rng])
    pointer = String(['_' for i=1:(pos-first(rng)-1)]) * "^"
    if length(substr) > maxchar
        # center the error char
        lst = min(pos+ceil(Int, hmaxchar), last(rng))
        fst = max(first(rng), pos-hmaxchar)
        substr = "..." * strip(str[fst:lst]) * "..."
        pointer = String(['_' for i=1:(pos-fst+2)]) * "^"
    end
    substr * "\n" * pointer
end

function quotedsplit(str, opts, includequotes, i=firstindex(str), l=lastindex(str))
    strtok = Quoted(StringToken(String), required=false,
                    includequotes=includequotes)

    f = Field(strtok, eoldelim=true)
    strs = String[]
    if l == 0
        return strs
    end
    while i <= l # this means that there was an empty field at the end of the line
        @chk2 x, i = tryparsenext(f, str, i, l, opts)
        push!(strs, x)
    end
    y1 = iterate(str, prevind(str, i))
    y1===nothing && error("Internal error.")
    c = y1[1]; i = y1[2]
    if c == opts.endchar
        # edge case where there's a delim at the end of the string
        push!(strs, "")
    end

    return strs
    @label error
    error("Couldn't split line, error at char $i:\n$(showerrorchar(str, i, 100))")
end

function canonical_name(opts, name)
    for list in opts
        if name in list
            return first(list)
        end
    end
    return name
end


text:
# Reference implementations for testing
# =====================================

module Reference

using LightGraphs: nv, vertices, neighbors, has_edge, add_edge!
using SimpleWeightedGraphs: SimpleWeightedGraph
using StatsFuns: logsumexp
using Random: shuffle

# Types used throughout the reference implementations
#   vertex identifier: Int
#   edge weight: Float64
const Graph = SimpleWeightedGraph{Int,Float64}
const Partition = Vector{BitSet}

function louvain(adjmat::AbstractMatrix;
                 partition::Partition = create_singleton_partition(adjmat),
                 resolution::Float64 = 1.0)
    return louvain(Graph(adjmat), partition = partition, resolution = resolution)
end

function louvain(graph::Graph;
                 partition::Partition = create_singleton_partition(graph),
                 resolution::Float64 = 1.0)
    @assert resolution > 0
    check_partition(graph, partition)
    stack = Partition[]
    cardinality = fill(1, nv(graph))
    @label loop
    partition = move_nodes(graph, cardinality, partition, resolution)
    check_partition(graph, partition)
    push!(stack, partition)
    if length(partition) != nv(graph)
        H_base = H(graph, cardinality, partition, resolution)
        graph, cardinality = aggregate_graph(graph, cardinality, partition)
        partition = create_singleton_partition(graph)
        @assert H(graph, cardinality, partition, resolution) ≈ H_base
        check_partition(graph, partition)
        @goto loop
    end
    return (quality = H(graph, cardinality, partition, resolution),
            partition = flatten(stack),)
end

function move_nodes(graph::Graph, cardinality::Vector{Int},
                    partition::Partition, resolution::Float64)
    γ = resolution
    partition = map(copy, partition)
    @label loop
    H_init = H(graph, cardinality, partition, γ)
    for u in shuffle(vertices(graph))
        dst, gain = find_best_community(graph, cardinality, partition, u, γ)
        if gain > 0
            move_node!(partition, u => dst)
        end
    end
    if H(graph, cardinality, partition, γ) > H_init
        @goto loop
    end
    return partition
end

function leiden(adjmat::AbstractMatrix;
                partition::Partition = create_singleton_partition(adjmat),
                resolution::Float64 = 1.0)
    return leiden(Graph(adjmat), partition = partition, resolution = resolution)
end

function leiden(graph::Graph;
                partition::Partition = create_singleton_partition(graph),
                resolution::Float64 = 1.0)
    @assert resolution > 0
    check_partition(graph, partition)
    stack = Partition[]
    cardinality = fill(1, nv(graph))
    @label loop
    partition = move_nodes_fast(graph, cardinality, partition, resolution)
    check_partition(graph, partition)
    if length(partition) != nv(graph)
        partition′ = refine_partition(graph, cardinality, partition, resolution)
        push!(stack, partition′)
        check_partition(graph, partition′)
        H_base = H(graph, cardinality, partition′, resolution)
        graph, cardinality = aggregate_graph(graph, cardinality, partition′)
        @assert H(graph, cardinality, create_singleton_partition(graph), resolution) ≈ H_base
        partition = map(partition) do comm
            new = create_empty_community()
            for (u, comm′) in zip(vertices(graph), partition′)
                if issubset(comm′, comm)
                    push!(new, u)
                end
            end
            @assert !isempty(new)
            new
        end
        check_partition(graph, partition)
        @goto loop
    end
    push!(stack, partition)
    return (quality = H(graph, cardinality, partition, resolution),
            partition = flatten(stack),)
end

function move_nodes_fast(graph::Graph, cardinality::Vector{Int},
                         partition::Partition, resolution::Float64)
    γ = resolution
    partition = map(copy, partition)
    queue = shuffle(1:nv(graph))
    while !isempty(queue)
        u = popfirst!(queue)
        dst, gain = find_best_community(graph, cardinality, partition, u, γ)
        if gain > 0
            move_node!(partition, u => dst)
            dst = findfirst(comm -> u ∈ comm, partition)  # move_node! may move communities
            ns = [v for v in neighbors(graph, u) if v ∉ partition[dst]]
            append!(queue, setdiff(ns, queue))
        end
    end
    return partition
end

function refine_partition(graph::Graph, cardinality::Vector{Int},
                          partition::Partition, resolution::Float64)
    partition′ = create_singleton_partition(graph)
    for comm in partition
        partition′ = merge_nodes_subset(graph, cardinality, partition′, comm, resolution)
    end
    return partition′
end

function merge_nodes_subset(graph::Graph, cardinality::Vector{Int},
                            partition::Partition, subset::BitSet, γ::Float64)
    # well-connectedness
    function is_well_connected(u::Int)
        a = tally_weights(graph, create_singleton_community(u), setdiff(subset, u))
        b = cardinality[u]
        c = tally_cardinality(subset, cardinality)
        return a ≥ γ * b * (c - b)
    end
    function is_well_connected(comm::BitSet)
        a = tally_weights(graph, comm, setdiff(subset, comm))
        b = tally_cardinality(comm, cardinality)
        c = tally_cardinality(subset, cardinality)
        return a ≥ γ * b * (c - b)
    end
    θ = 0.01
    partition = map(copy, partition)
    for u in subset
        if !is_well_connected(u)
            continue
        end
        i = findfirst(comm -> u ∈ comm, partition)
        @assert !isnothing(i)
        if is_singleton(partition[i])
            logprobs = fill(-Inf, length(partition))
            for (i, comm) in enumerate(partition)
                if issubset(comm, subset) && is_well_connected(comm)
                    gain = ΔH(graph, cardinality, partition, u => i, γ)
                    logprobs[i] = gain ≥ 0 ? 1/θ * gain : -Inf
                end
            end
            if !isfinite(logsumexp(logprobs))
                break
            end
            move_node!(partition, u => sample(logprobs))
        end
    end
    return partition
end

function sample(logprobs::Vector{Float64})
    probs = exp.(logprobs .- logsumexp(logprobs))
    r = rand()
    i = 1
    cumprob = 0.0
    while i < lastindex(probs)
        cumprob += probs[i]
        if cumprob ≥ r
            break
        end
        i += 1
    end
    return i
end

# Find the best community to which `u` is assigned.
function find_best_community(graph::Graph, cardinality::Vector{Int},
                             partition::Partition, u::Int, γ::Float64)
    # new community?
    dst = lastindex(partition) + 1
    push!(partition, create_empty_community())
    maxgain = ΔH(graph, cardinality, partition, u => dst, γ)
    pop!(partition)
    # or existing communities?
    for i in 1:lastindex(partition)
        val = ΔH(graph, cardinality, partition, u => i, γ)
        if val > maxgain
            dst = i
            maxgain = val
        end
    end
    return dst, maxgain
end

function ΔH(graph::Graph, cardinality::Vector{Int},
            partition::Partition, move::Pair{Int,Int}, γ::Float64)
    u, dst = move
    partition′ = map(enumerate(partition)) do (i, comm)
        if i == dst
            # `u` can be in `comm`
            comm = copy(comm)
            push!(comm, u)
        elseif u ∈ comm
            comm = copy(comm)
            delete!(comm, u)
        end
        return comm
    end
    return H(graph, cardinality, partition′, γ) - H(graph, cardinality, partition, γ)
end

function H(graph::Graph, cardinality::Vector{Int}, partition::Partition, γ::Float64)
    result = 0.0
    for (i, comm) in enumerate(partition)
        weight = tally_weights(graph, partition, i, i)
        n = isempty(comm) ? 0 : sum(cardinality[u] for u in comm)
        result += weight - γ * div(n * (n - 1), 2)
    end
    return result
end

function move_node!(partition::Partition, move::Pair{Int,Int})
    u, dst = move
    src = findfirst(comm -> u ∈ comm, partition)
    @assert !isnothing(src)
    if src == dst
        return partition
    end
    if dst > lastindex(partition)
        # new community
        push!(partition, create_singleton_community(u))
    else
        # existing community
        push!(partition[dst], u)
    end
    delete!(partition[src], u)
    if isempty(partition[src])
        deleteat!(partition, src)
    end
    return partition
end

function aggregate_graph(graph::Graph, cardinality::Vector{Int}, partition::Partition)
    @assert all(!isempty, partition)
    nv = length(partition)
    aggregated = SimpleWeightedGraph(nv)
    for i in 1:nv, j in i:nv
        weight = tally_weights(graph, partition, i, j)
        @assert weight ≥ 0
        if weight > 0
            add_edge!(aggregated, i, j, weight)
        end
    end
    return aggregated, [sum(cardinality[u] for u in comm) for comm in partition]
end

function tally_weights(graph::Graph)
    n = nv(graph)
    weight = 0.0
    for i in 1:n, j in i:n
        weight += graph.weights[i,j]
    end
    return weight
end

function tally_weights(graph::Graph, partition::Partition, i::Int, j::Int)
    weight = 0.0
    if i == j
        comm = partition[i]
        for u in comm, v in comm
            if u ≤ v
                weight += graph.weights[u,v]
            end
        end
    else
        for u in partition[i], v in partition[j]
            weight += graph.weights[u,v]
        end
    end
    return weight
end

function tally_weights(graph::Graph, comm1::BitSet, comm2::BitSet)
    @assert isdisjoin(comm1, comm2)
    weight = 0.0
    for u in comm1, v in comm2
        if has_edge(graph, u, v)
            weight += graph.weights[u,v]
        end
    end
    return weight
end

function isdisjoin(set1::BitSet, set2::BitSet)
    for x in set1, y in set2
        if x == y
            return false
        end
    end
    return true
end

function tally_cardinality(comm::BitSet, cardinality::Vector{Int})
    return isempty(comm) ? 0 : sum(cardinality[u] for u in comm)
end

function flatten(partitions::Vector{Partition})
    result = map(copy, last(partitions))
    for i in lastindex(partitions):-1:2
        partition = partitions[i-1]
        for j in 1:lastindex(result)
            result[j] = mapfoldl(x -> partition[x], union, result[j],
                                 init = create_empty_community())
        end
    end
    sort!(result, by = length)
    return result
end

function create_singleton_partition(graph::Graph)
    return [create_singleton_community(u) for u in vertices(graph)]
end

function create_singleton_partition(adjmat::AbstractMatrix)
    return [create_singleton_community(u) for u in 1:size(adjmat, 1)]
end

function create_empty_community()
    return BitSet()
end

function create_singleton_community(u::Int)
    return BitSet(u)
end

function is_singleton(comm::BitSet)
    return length(comm) == 1
end

function check_partition(graph::Graph, partition::Partition)
    n = nv(graph)
    found = BitSet()
    for comm in partition
        if !isempty(comm ∩ found)
            throw(AssertionError("non-disjoint partition: $(collect(comm ∩ found))"))
        elseif maximum(comm) > n
            throw(AssertionError("out of bounds partition: $([u for u in comm if u > n])"))
        end
        union!(found, comm)
    end
    if length(found) != n
        throw(AssertionError("non-exhaustive partition: $(collect(setdiff(BitSet(1:n), found)))"))
    end
end

end


text:
using Mocking

batch_delete_image_1_patch = @patch function batch_delete_image(a...; b...)
    return Dict{String, Any}("imageIds" => Any[Dict{String, Any}("imageDigest" => "sha256:examplee6d1e504117a17000003d3753086354a38375961f2e665416ef4b1b2f", "imageTag" => "precise")], "failures" => Any[])
end

batch_get_image_1_patch = @patch function batch_get_image(a...; b...)
    return Dict{String, Any}("images" => Any[Dict{String, Any}("registryId" => "244698725403", "repositoryName" => "ubuntu", "imageManifest" => "{\n \"schemaVersion\": 1,\n \"name\": \"ubuntu\",\n \"tag\": \"precise\",\n...", "imageId" => Dict{String, Any}("imageDigest" => "sha256:example76bdff6d83a09ba2a818f0d00000063724a9ac3ba5019c56f74ebf42a", "imageTag" => "precise"))], "failures" => Any[])
end

create_repository_1_patch = @patch function create_repository(a...; b...)
    return Dict{String, Any}("repository" => Dict{String, Any}("repositoryArn" => "arn:aws:ecr:us-west-2:012345678901:repository/project-a/nginx-web-app", "registryId" => "012345678901", "repositoryName" => "project-a/nginx-web-app"))
end

delete_repository_1_patch = @patch function delete_repository(a...; b...)
    return Dict{String, Any}("repository" => Dict{String, Any}("repositoryArn" => "arn:aws:ecr:us-west-2:012345678901:repository/ubuntu", "registryId" => "012345678901", "repositoryName" => "ubuntu"))
end

delete_repository_policy_1_patch = @patch function delete_repository_policy(a...; b...)
    return Dict{String, Any}("registryId" => "012345678901", "repositoryName" => "ubuntu", "policyText" => "{ ... }")
end

describe_repositories_1_patch = @patch function describe_repositories(a...; b...)
    return Dict{String, Any}("repositories" => Any[Dict{String, Any}("repositoryArn" => "arn:aws:ecr:us-west-2:012345678910:repository/ubuntu", "registryId" => "012345678910", "repositoryName" => "ubuntu"), Dict{String, Any}("repositoryArn" => "arn:aws:ecr:us-west-2:012345678910:repository/test", "registryId" => "012345678910", "repositoryName" => "test")])
end

get_authorization_token_1_patch = @patch function get_authorization_token(a...; b...)
    return Dict{String, Any}("authorizationData" => Any[Dict{String, Any}("proxyEndpoint" => "https://012345678901.dkr.ecr.us-west-2.amazonaws.com", "authorizationToken" => "QVdTOkN...", "expiresAt" => "1470951892432")])
end

get_repository_policy_1_patch = @patch function get_repository_policy(a...; b...)
    return Dict{String, Any}("registryId" => "012345678901", "repositoryName" => "ubuntu", "policyText" => "{\n  \"Version\" : \"2008-10-17\",\n  \"Statement\" : [ {\n    \"Sid\" : \"new statement\",\n    \"Effect\" : \"Allow\",\n    \"Principal\" : {\n     \"AWS\" : \"arn:aws:iam::012345678901:role/CodeDeployDemo\"\n    },\n\"Action\" : [ \"ecr:GetDownloadUrlForLayer\", \"ecr:BatchGetImage\", \"ecr:BatchCheckLayerAvailability\" ]\n } ]\n}")
end

list_images_1_patch = @patch function list_images(a...; b...)
    return Dict{String, Any}("imageIds" => Any[Dict{String, Any}("imageDigest" => "sha256:764f63476bdff6d83a09ba2a818f0d35757063724a9ac3ba5019c56f74ebf42a", "imageTag" => "precise")])
end


text:
using Polyopt
using Base.Test

include("polyexpr.jl")
include("polyopt.jl")
include("pooling.jl")
include("bsos.jl")


text:
"""
    MSA

Describes parameters of an MSA.
"""
Base.@kwdef struct MSA
    id::String					# Pfam id (required)
    aln::String = "full"		# -Full
    format::String = "fasta"	# -Format:FASTA
    order::String = "t"			# -Order:Tree
    case::String = "l"			# -Sequence:Inserts lower case
    gaps::String = "default"	# -Gaps:Gaps as "." or "-"(mixed)

    function MSA(
        id::String, aln::String, format::String, order::String, case::String, gaps::String
    )
        if is_valid_pfam_id(id)
            return new(id, aln, format, order, case, gaps)
        else
            error("$id is not a valid Pfam code")
        end
    end
end

MSA(id::String; kwargs...) = MSA(; id=id, kwargs...)

function is_valid_pfam_id(id::String)
    return occursin(r"^PF\d{5}$"i, id)
end

function url(msa::MSA)
    return (
        "https://pfam.xfam.org/family/$(msa.id)/alignment/full/format?" *
        "format=$(msa.format)&" *
        "alnType=$(msa.aln)&" *
        "order=$(msa.order)&" *
        "case=$(msa.case)&" *
        "gaps=$(msa.gaps)&" *
        "download=1"
    )
end

function file(msa::MSA)
    return (
        "pfam_$(msa.id)_" *
        "format=$(msa.format)_" *
        "alnType=$(msa.aln)_" *
        "order=$(msa.order)_" *
        "case=$(msa.case)_" *
        "gaps=$(msa.gaps)"
    )
end

path(msa::MSA) = joinpath(pfam_msa_scratch, file(msa))

function download(msa::MSA)
    @info "Downloading $(msa.id)"
    #= PFAM has expired certificates. So we need this or `download` will fail.
    TODO: check if the certificates at PFAM got fixed =#
    withenv("JULIA_NO_VERIFY_HOSTS" => "pfam.xfam.org", ENV...) do
        Downloads.download(url(msa), path(msa); timeout = Inf)
    end
end

"""
    load(msa)

Loads an MSA.
"""
function load(msa::MSA)
    if !isfile(path(msa))
        download(msa)
    end

    df = DataFrame(uniprot_id = String[], start = Int[], stop = Int[], sequence = String[])
    FastaReader(path(msa)) do fr
        for (desc, seq) in fr
            id, start, stop = split(desc, ['>', '/', '-']; keepempty=false)
            push!(df, (id, parse(Int, start), parse(Int, stop), seq))
        end
    end

    return df
end

"""
    remove_inserts(sequence)

Remove inserts from a string sequence.
"""
function remove_inserts(sequence::String)
    return filter(c -> c == '-' || c ≠ lowercase(c), sequence)
end


text:
@testset "Node" begin
    # Test full construction
    node = Node{2}(12, Vec{2}((2.3, 4.5)), [1,7,9], 2)

    @test node.id == 12
    @test node.coords == [2.3, 4.5]
    @test node.dofs == [1,7,9]

    # Test DofObject interface
    @test dofs(node) == [1,7,9]
    @test processor_id(node) == 2
end


text:
module KelvinletsVideo
    include("tetrahedronInterpolator.jl")
    using Images, ProgressMeter, ImageView, LinearAlgebra
    export KelvinletsVideoObject, grab

    struct KelvinletsVideoObject
        a::Float64
        b::Float64
        c::Float64
        sizeX::Int64
        sizeY::Int64
        sizeZ::Int64
        video::AbstractArray{RGB{N0f8}, 3}
        function KelvinletsVideoObject(video::AbstractArray{RGB{N0f8}, 3},
                                        ν::Float64,
                                        μ::Float64
                )::KelvinletsObject

            a = 1 / (4pi * μ)
            b = a / (4(1 - ν))
            c = 2 / (3a- 2b)

            sizeY, sizeX, sizeZ = size(video)

            new(a, b, c, sizeX, sizeY, sizeZ, video)
        end
    end
    function __applyVariation__(object::KelvinletsVideoObject,
                                pressurePoint::Array{Int64, 1},
                                variationCalculator,
                                retardationFunction)::Array{RGB{N0f8}, 3}

        allΔ = zeros(object.sizeY, object.sizeX, object.sizeZ, 3)
        print("Applying Variation")
        @showprogress for i=1:object.sizeY
            for j=1:object.sizeX
                for k=1:object.sizeZ
                    Δ = variationCalculator([i, j, k])

                    dx1 = j - 1
                    dx2 = object.sizeX - j
                    dy1 = i - 1
                    dy2 = object.sizeY - i
                    dz1 = k - 1
                    dz2 = object.sizeZ - k

                    dx = min(dx1, dx2)
                    dy = min(dy1, dy2)
                    dz = min(dz1, dz2)

                    y = 2(object.sizeY/2 - (dy))/object.sizeY
                    x = 2(object.sizeX/2 - (dx))/object.sizeX
                    z = 2(object.sizeZ/2 - (dz))/object.sizeZ

                    Δ[1] *= retardationFunction(y)
                    Δ[2] *= retardationFunction(x)
                    Δ[3] *= retardationFunction(z)

                    Δ += [i, j, k]

                    allΔ[i, j, k, :] = Δ[:]

                end
            end
        end
        __interpolateVariation__(object, allΔ)
    end
    function __interpolateVariation__(object::KelvinletsVideoObject,
                                      allΔ::Array{Float64, 4})::Array{RGB{N0f8}, 3}

        interpVideo = fill(RGB{N0f8}(0, 0, 0), object.sizeY, object.sizeX, object.sizeZ)

        rasterize = function(A, B, C, D)
            colorA = object.video[A[1], A[2], A[3]]
            colorB = object.video[B[1], B[2], B[3]]
            colorC = object.video[C[1], C[2], C[3]]
            colorD = object.video[D[1], D[2], D[3]]
            TetrahedronInterpolator.rasterization(interpVideo,
                            allΔ[A[1], A[2], A[3], :],
                            allΔ[B[1], B[2], B[3], :],
                            allΔ[C[1], C[2], C[3], :],
                            allΔ[D[1], D[2], D[3], :],
                            colorA, colorB, colorC, colorD
            )
        end
        print("Applying Interpolation")
        @showprogress for i=1:object.sizeY-1
            for j=1:object.sizeX-1
                for k=1:object.sizeZ-1
                    #(0,0,0), (1,0,1), (1,0,0), (1,1,1)
                    rasterize([i, j, k],
                              [i+1, j, k+1],
                              [i+1, j, k],
                              [i+1, j+1, k+1]
                             )
                    #(0,0,0), (1,0,0), (1,1,0), (1,1,1)
                    rasterize([i, j, k],
                              [i+1, j, k],
                              [i+1, j+1, k],
                              [i+1, j+1, k+1]
                             )
                    #(0,0,0), (0,1,0), (1,1,0), (1,1,1)
                    rasterize([i, j, k],
                              [i, j+1, k],
                              [i+1, j+1, k],
                              [i+1, j+1, k+1]
                             )
                    #(0,0,0), (0,1,0), (0,1,1), (1,1,1)
                    rasterize([i, j, k],
                              [i, j+1, k],
                              [i, j+1, k+1],
                              [i+1, j+1, k+1]
                             )
                    #(0,0,0), (0,1,1), (0,0,1), (1,1,1)
                    rasterize([i, j, k],
                              [i, j+1, k+1],
                              [i, j, k+1],
                              [i+1, j+1, k+1]
                             )
                    #(0,0,0), (1,0,1), (0,0,1), (1,1,1)
                    rasterize([i, j, k],
                              [i+1, j, k+1],
                              [i, j, k+1],
                              [i+1, j+1, k+1]
                             )
                end
            end
        end
        interpVideo
    end
    function grab(object::KelvinletsVideoObject,
                  x0::Array{Int64, 1},
                  force::Array{Float64, 1},
                  ϵ::Float64)::Array{RGB{N0f8}, 3}
            grabFunc = function(x)
                r = x - x0
                rLength = norm(r)
                rϵ = sqrt(rLength^2 + ϵ^2)
                kelvinState = (((object.a - object.b)/rϵ) * I +
                                (object.b / rϵ^3) * r * r' +
                                (object.a / 2)*((ϵ^2) / rϵ^3) * I)

                object.c * ϵ * kelvinState * force
            end
            retardationFunc = x -> (cos(π*x) + 1)/2
            __applyVariation__(object, x0, grabFunc, retardationFunc)
        end
end


text:
using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking

Turing.setprogress!(false)

df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
df = df[df.age .>= 18, :]

@model function ppl4_1(heights)
    μ ~ Normal(178, 20)
    σ ~ Uniform(0, 50)
    heights .~ Normal(μ, σ)
end

m4_1t = ppl4_1(df.height)
precis(m4_1t)

res = quap(m4_1t)