text:
@inline function eatwhitespaces(str::Union{VectorBackedUTF8String, String}, i=1, len=lastindex(str))
    while i<=len
        @inbounds b = codeunit(str, i)

        if b==0x20 # This is ' '
            i += 1
        else
            break
        end
    end
    return i
end

@inline function eatnewlines(str::Union{VectorBackedUTF8String, String}, i=1, len=lastindex(str))
    count = 0
    while i<=len
        @inbounds b = codeunit(str, i)
        if b == 0xd # '\r'
            i += 1
            if i<=len
                @inbounds b = codeunit(str, i)
                if b == 0xa # '\n'
                    i += 1
                end
            end
            count += 1
        elseif b == 0xa
            i += 1
            if i<=len
                @inbounds b = codeunit(str, i)
                if b == 0xd
                    i += 1
                end
            end
            count += 1
        else
            break
        end
    end

    return i, count
end

@inline function tryparsenext_base10_digit(T,str::Union{VectorBackedUTF8String, String},i, len)
    i > len && @goto error
    @inbounds b = codeunit(str,i)
    diff = b-0x30
    diff >= UInt8(10) && @goto error
    return convert(T, diff), i+1

    @label error
    return nothing
end

@inline _isdigit(b::UInt8) = ( (0x30 ≤ b) & (b ≤ 0x39) )

@inline function parse_uint_and_stop(str::Union{VectorBackedUTF8String, String}, i, len, n::T) where {T <: Integer}
    ten = T(10)
    # specialize handling of the first digit so we can return an error
    max_without_overflow = div(typemax(T)-9,10) # the larg
    i <= len || return n, false, i
    @inbounds b = codeunit(str, i)
    diff = b-0x30
    if diff < UInt8(10) && n <= max_without_overflow
        n *= ten
        n += T(diff)
    else
        return n, false, i
    end
    i += 1

    while i <= len && n <= max_without_overflow
        @inbounds b = codeunit(str, i)
        diff = b-0x30
        if diff < UInt8(10)
            n *= ten
            n += T(diff)
        else
            return n, true, i
        end
        i += 1
    end
    return n, true, i
end

@inline function read_digits(str::Union{VectorBackedUTF8String, String}, i, len)
    # slurp up extra digits
    while i <= len
        @inbounds b = codeunit(str, i)
        if !_isdigit(b) # do nothing
            return i
        end
        i += 1
    end
    return i
end

@inline function _is_e(str::Union{VectorBackedUTF8String, String}, i)
    @inbounds b = codeunit(str, i)
    return  (b==0x65) | (b==0x45)
end

@inline function _is_negative(str::Union{VectorBackedUTF8String, String}, i)
    @inbounds b = codeunit(str, i)
    return b==0x2d
end

@inline function _is_positive(str::Union{VectorBackedUTF8String, String}, i)
    @inbounds b = codeunit(str, i)
    return b==0x2b
end

const pre_comp_exp = Float64[10.0^i for i=0:22]

@inline function tryparsenext(::Numeric{F}, str::Union{VectorBackedUTF8String, String}, i, len) where {F<:AbstractFloat}
    R = Nullable{F}

    i>len && @goto error

    negate = false
    @inbounds b = codeunit(str, i)
    if b==0x2d # '-'
        negate = true
        i += 1
    elseif b==0x2b # '+'
        i +=1
    end

    f1::Int64 = 0

    # read an integer up to the decimal point
    f1, rval1, idecpt = parse_uint_and_stop(str, i, len, f1)
    idecpt = read_digits(str, idecpt, len) # get any trailing digits
    i = idecpt

    ie = i
    frac_digits = 0

    # next thing must be dec pt.
    if i <= len && @inbounds(codeunit(str, i)) == 0x2e # Check for '.'
        i += 1
        f1, rval2, ie = parse_uint_and_stop(str, i, len, f1)
        frac_digits = ie - i

        ie = read_digits(str, ie, len) # get any trailing digits
    elseif !rval1 # no first number, and now no deciaml point => invalid
        @goto error
    end

    # Next thing must be exponent
    i = ie
    eval::Int32 = 0

    if i <= len && _is_e(str, i)
        i += 1

        enegate = false
        if i<=len
            if _is_negative(str, i)
                enegate = true
                i += 1
            elseif _is_positive(str, i)
                i += 1
            end
        end
        eval, rval3, i = parse_uint_and_stop(str, i, len, eval)
        if enegate
            eval *= Int32(-1)
        end
    end

    exp = eval - frac_digits

    maxexp = 308
    minexp = -307

    if frac_digits <= 15 && -22 <= exp <= 22
        if exp >= 0
            f = F(f1)*pre_comp_exp[exp+1]
        else
            f = F(f1)/pre_comp_exp[-exp+1]
        end
    else
          f = convert_to_double(f1, exp)
    end

    if negate
        f = -f
    end

    @label done
    return R(convert(F, f)), i

    @label error
    return R(), i
end

function tryparsenext(f::Field{T}, str::Union{VectorBackedUTF8String, String}, i, len, opts::LocalOpts{T_ENDCHAR}) where {T, T_ENDCHAR<:UInt8}
    R = Nullable{T}
    i > len && @goto error
    if f.ignore_init_whitespace
        i = eatwhitespaces(str, i, len)
    end
    @chk2 res, i = tryparsenext(f.inner, str, i, len, opts)

    if f.ignore_end_whitespace
        i0 = i

        while i<=len
            @inbounds b = codeunit(str, i)

            !opts.spacedelim && opts.endchar == 0x09 && b == 0x09 && (i = i+1; @goto done) # 0x09 is \t

            b!=0x20 && b!=0x09 && break
            i=i+1
        end

        opts.spacedelim && i > i0 && @goto done
    end
    # todo don't ignore whitespace AND spacedelim

    if i > len
        if f.eoldelim
            @goto done
        else
            @goto error
        end
    end

    i>len && error("Internal error.")
    @inbounds b = codeunit(str, i)
    opts.spacedelim && (b!=0x20 || b!=0x09) && (i+=1; @goto done)
    !opts.spacedelim && opts.endchar == b && (i+=1; @goto done)

    if f.eoldelim
        if b == 0x0d # '\r'
            i+=1
            if i<=len
                @inbounds b = codeunit(str, i)
                if b == 0x0a # '\n'
                    i+=1
                end
            end
            @goto done
        elseif b == 0x0a # '\n'
            i+=1
            if i<=len
                @inbounds b = codeunit(str, i)
                if b == 0x0d # '\r'
                    i+=1
                end
            end
            @goto done
        end
    end

    @label error
    return R(), i

    @label done
    return R(convert(T, res)), i
end

function tryparsenext(q::Quoted{T,S,<:UInt8,<:UInt8}, str::Union{VectorBackedUTF8String, String}, i, len, opts::LocalOpts{<:UInt8,<:UInt8,<:UInt8}) where {T,S}
    if i>len
        q.required && @goto error
        # check to see if inner thing is ok with an empty field
        @chk2 x, i = tryparsenext(q.inner, str, i, len, opts) error
        @goto done
    end
    @inbounds b = codeunit(str, i)
    ii = i+1
    quotestarted = false
    if q.quotechar == b
        quotestarted = true
        if !q.includequotes
            i = ii
        end

        if q.stripwhitespaces
            i = eatwhitespaces(str, i, len)
        end
    else
        q.required && @goto error
    end

    if quotestarted
        qopts = LocalOpts(q.quotechar, false, q.quotechar, q.escapechar,
                         q.includequotes, q.includenewlines)
        @chk2 x, i = tryparsenext(q.inner, str, i, len, qopts)
    else
        @chk2 x, i = tryparsenext(q.inner, str, i, len, opts)
    end

    if i > len
        if quotestarted && !q.includequotes
            @goto error
        end
        @goto done
    end

    if q.stripwhitespaces
        i = eatwhitespaces(str, i, len)
    end
    i>len && error("Internal error.")
    @inbounds b = codeunit(str, i)
    ii = i + 1

    if quotestarted && !q.includequotes
        b != q.quotechar && @goto error
        i = ii
    end

    @label done
    return Nullable{T}(x), i

    @label error
    return Nullable{T}(), i
end

@inline function isnewline(b::UInt8)
    b == UInt8(10) || b == UInt8(13)
end

function tryparsenext(s::StringToken{T}, str::Union{VectorBackedUTF8String, String}, i, len, opts::LocalOpts{<:UInt8,<:UInt8,<:UInt8}) where {T}
    len = ncodeunits(str)
    inside_quoted_strong = opts.endchar == opts.quotechar
    escapecount = 0
    R = Nullable{T}
    p = UInt8(0)
    i0 = i
    if opts.includequotes
        if i<=len
            @inbounds b = codeunit(str, i)
            if b==opts.quotechar
                # advance counter so that
                # the while loop doesn't react to opening quote
                i += 1
            end
        end
    end

    while i<=len
        @inbounds b = codeunit(str, i)
        ii = i + 1

        if inside_quoted_strong && p==opts.escapechar
            escapecount += 1
        end

        if opts.spacedelim && (b == UInt8(32) || b == UInt8(9)) # 32 = ' ' and 9 = '\t'
            break
        elseif !opts.spacedelim && b == opts.endchar
            if inside_quoted_strong
                # this means we're inside a quoted string
                if opts.quotechar == opts.escapechar
                    # sometimes the quotechar is the escapechar
                    # in that case we need to see the next char
                    if ii > len
                        if opts.includequotes
                            i=ii
                        end
                        break
                    else
                        @inbounds next_b = codeunit(str, ii)
                        if next_b == opts.quotechar
                            # the current character is escaping the
                            # next one
                            i = ii + 1 # skip next char as well
                            p = next_b
                            continue
                        end
                    end
                elseif p == opts.escapechar
                    # previous char escaped this one
                    i = ii
                    p = b
                    continue
                end
            end
            if opts.includequotes
                i = ii
            end
            break
        elseif (!opts.includenewlines && isnewline(b))
            break
        end
        i = ii
        p = b
    end

    return R(_substring(T, str, i0, i-1, escapecount, opts)), i
end

@inline function _substring(::Type{String}, str::Union{VectorBackedUTF8String, String}, i, j, escapecount, opts::LocalOpts{<:UInt8,<:UInt8,<:UInt8})
    if escapecount > 0
        buffer = Vector{UInt8}(undef, j-i+1-escapecount)
        cur_i = i
        cur_buffer_i = 1
        @inbounds c = codeunit(str, cur_i)
        if opts.includequotes && c==opts.quotechar
            @inbounds buffer[cur_buffer_i] = c
            cur_i += 1
            cur_buffer_i += 1
        end
        while cur_i <= j
            @inbounds c = codeunit(str, cur_i)
            if c == opts.escapechar
                next_i = cur_i + 1
                if next_i <= j
                    @inbounds next_c = codeunit(str, next_i)
                    if next_c == opts.quotechar
                        @inbounds buffer[cur_buffer_i] = next_c
                        cur_buffer_i += 1
                        cur_i = next_i
                    end
                else
                    @inbounds buffer[cur_buffer_i] = c
                    cur_buffer_i += 1
                end
            else
                @inbounds buffer[cur_buffer_i] = c
                cur_buffer_i += 1
            end
            cur_i += 1
        end
        return String(buffer)
    else
        return unsafe_string(pointer(str, i), j-i+1)
    end
end


text:
module TestOptics

using Accessors
using Accessors: mapproperties
using Test
import ConstructionBase

@testset "mapproperties" begin
    res = @inferred mapproperties(x->2x, (a=1, b=2))
    @test res === (a=2, b=4)
    @test NamedTuple() === @inferred mapproperties(cos, NamedTuple())
    struct S{A,B}
        a::A
        b::B
    end
    res = @inferred mapproperties(x->2x, S(1, 2.0))
    @test res === S(2, 4.0)

    # overloading
    struct AB
        a::Int
        b::Int
        _checksum::UInt
        AB(a,b) = new(a,b,hash(a,hash(b)))
    end

    ConstructionBase.getproperties(o::AB) = (a=o.a, b=o.b)
    ConstructionBase.setproperties(o::AB, patch::NamedTuple) = AB(patch.a, patch.b)
    ab = AB(1,2)
    ab2 = @inferred mapproperties(x -> 2x, ab)
    @test ab2 === AB(2,4)
end

@testset "Properties" begin
    pt = (x=1, y=2, z=3)
    @test (x=0, y=1, z=2) === @set pt |> Properties() -= 1
    @inferred modify(x->x-1, pt, Properties())

    # custom struct
    struct Point{X,Y,Z}
        x::X; y::Y; z::Z
    end
    pt = Point(1f0, 2e0, 3)
    pt2 = @inferred modify(x->2x, pt, Properties())
    @test pt2 === Point(2f0, 4e0, 6)
end

@testset "Elements" begin

    @test [0,0,0] == @set 1:3 |> Elements() = 0

    arr = 1:3
    @test 2:4 == (@set arr |> Elements() += 1)
    @test map(cos, arr) == modify(cos, arr, Elements())

    @test modify(cos, (), Elements()) === ()

    @inferred modify(cos, arr, Elements())
    @inferred modify(cos, (), Elements())
end

@testset "Recursive" begin
    obj = (a=1, b=(1,2), c=(A=1, B=(1,2,3), D=4))
    rp = Recursive(x -> !(x isa Tuple), Properties())
    @test modify(collect, obj, rp) == (a = 1, b = [1, 2], c = (A = 1, B = [1, 2, 3], D = 4))

    arr = [1,2,[3,4], [5, 6:7,8, 9,]]
    oc = Recursive(x -> x isa AbstractArray, Elements())
    expected = [0,1,[2,3], [4, 5:6,7, 8,]]
    @test modify(x-> x - 1, arr, oc) == expected
end

@testset "If" begin
    @test 10 === @set(1 |> If(>=(0)) = 10)
    @test -1 === @set(-1 |> If(>=(0)) = 10)
    @inferred set(1, If(iseven), 2)
    @inferred modify(x -> 0, 1, If(iseven))

    arr = 1:6
    @test [1, 0, 3, 0, 5, 0] == @set(arr |> Elements() |> If(iseven) = 0)
    @inferred modify(x -> 0, arr, @optic _ |> Elements() |> If(iseven))
end

end#module

ext: jl


text:
leaderuri(cli::Client) = "http://$(cli.host):$(cli.port)/$(cli.version)/stats/leader"

"""
    leader(cli) -> Dict

Fetchers the stats for the leader of the etcd cluster.
See the [etcd API](https://github.com/coreos/etcd/blob/master/Documentation/v2/api.md#leader-statistics)
for more details.
"""
function leader(cli::Client)
    leader_id = request(Requests.get, leaderuri(cli), Dict())["leader"]
    return members(cli)[leader_id]
end


text:
include("activate.jl")

import LinearAlgebra

n = Int(1e5)
k = 10000

xpi = rand(n);
total = sum(xpi);
const xp = xpi ./ total;

f(x) = LinearAlgebra.norm(x - xp)^2

function grad!(storage, x)
    storage .= 2 * (x - xp)
    return nothing
end

# better for memory consumption as we do coordinate-wise ops

function cf(x, xp)
    return @. LinearAlgebra.norm(x - xp)^2
end

# lmo = FrankWolfe.KSparseLMO(100, 1.0)
lmo = FrankWolfe.LpNormLMO{Float64,1}(1.0)
# lmo = FrankWolfe.ProbabilitySimplexOracle(1.0);
# lmo = FrankWolfe.UnitSimplexOracle(1.0);
x00 = FrankWolfe.compute_extreme_point(lmo, zeros(n))
# print(x0)

gradient = similar(x00)

FrankWolfe.benchmark_oracles(
    f,
    grad!,
    () -> randn(n),
    lmo;
    k=100,
)

# 1/t *can be* better than short step

println("\n==> Short Step rule - if you know L.\n")

x0 = copy(x00)
@time x, v, primal, dual_gap, trajectory_shortstep = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

println("\n==> Short Step rule with momentum - if you know L.\n")

x0 = copy(x00)

@time x, v, primal, dual_gap, trajectoryM = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.blas,
    verbose=true,
    trajectory=true,
    momentum=0.9,
);

println("\n==> Adaptive if you do not know L.\n")

x0 = copy(x00)

@time x, v, primal, dual_gap, trajectory_adaptive = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    L=100,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

println("\n==> Agnostic if function is too expensive for adaptive.\n")

x0 = copy(x00)

@time x, v, primal, dual_gap, trajectory_agnostic = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Agnostic(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

data = [trajectory_shortstep, trajectory_adaptive, trajectory_agnostic, trajectoryM]
label = ["short step" "adaptive" "agnostic" "momentum"]

FrankWolfe.plot_trajectories(data, label)


text:
module ODEToolsTests

using Test

@testset "DiffOperators" begin include("DiffOperatorsTests.jl") end

@testset "ODEOperators" begin include("ODEOperatorsTests.jl") end

# @testset "ODESolvers" begin include("ODESolversTests.jl") end

end # module


text:
using EarCut, GeometryTypes
using Test

# just check basic wrapper functionality. Check for correct size and that data is not garbage
for typ in (Float32, Float64, Int32, Int64)
    @testset "Testing polygon with element type $typ" begin
        usborder = [
          Point{2, typ}[
            (1130,1713),(1131,1710),(1137,1731),(1133,1752),(1125,1753),(1118,1742),(1110,1717),(1105,1718),(1108,1704),(1096,1691),(1077,1694),(1067,1683),(1019,1687),(1031,1689),(1031,1704),(1022,1696),(1022,1702),(1010,1700),(1003,1692),(998,1696),(980,1690),(970,1698),(966,1694),(966,1702),(938,1718),(943,1742),(920,1736),(916,1721),(894,1693),(884,1691),(872,1703),(837,1667),(785,1672),(743,1654),(715,1656),(699,1636),(676,1628),(654,1587),(656,1583),(660,1588),(657,1579),(649,1580),(633,1547),(637,1529),(631,1507),(638,1454),(647,1454),(637,1452),(639,1441),(635,1442),(629,1417),(651,1421),(647,1434),(655,1428),(650,1440),(656,1434),(654,1423),(651,1420),(653,1419),(651,1407),(965,1407),(966,1400),(972,1411),(1008,1423),(1043,1419),(1083,1442),(1086,1450),(1091,1448),(1109,1468),(1114,1496),(1102,1520),(1107,1525),(1149,1508),(1147,1498),(1152,1495),(1174,1495),(1195,1474),(1242,1470),(1260,1433),(1277,1440),(1277,1462),(1286,1476),(1274,1484),(1265,1480),(1243,1503),(1240,1516),(1252,1526),(1238,1529),(1236,1523),(1234,1530),(1218,1531),(1206,1540),(1205,1554),(1195,1567),(1188,1556),(1194,1574),(1185,1590),(1187,1581),(1179,1567),(1185,1557),(1176,1562),(1180,1579),(1179,1585),(1170,1577),(1180,1593),(1169,1590),(1183,1596),(1186,1607),(1175,1605),(1183,1613),(1182,1618),(1171,1615),(1179,1624),(1167,1626),(1145,1650),(1132,1659),(1128,1656),(1121,1675),(1131,1708),(1129,1710),(1130,1713)],
          Point{2, typ}[(654,1419),(653,1419),(654,1423),(656,1425),(654,1419)]
        ]
        triangles = triangulate(usborder)
        @test length(triangles) == 120
        @test first(triangles) == GLTriangle(1, 2, 3)
        @test last(triangles)  == GLTriangle(124, 58, 123)
    end
end


text:
using Oceananigans

# Idealized convection in a rectilinear in x-z
grid = RectilinearGrid(size=(64, 64), x=(0, 128), z=(-64, 0), halo=(3, 3),
                       topology=(Periodic, Flat, Bounded))

# Build a model
top_buoyancy_bc = FluxBoundaryCondition(1e-8)
buoyancy_bcs = FieldBoundaryConditions(top=top_buoyancy_bc)

closure = ScalarDiffusivity(ν=1e-4, κ=1e-4)

model = NonhydrostaticModel(grid = grid,
                            advection = WENO5(),
                            tracers = :b,
                            closure = closure,
                            buoyancy = BuoyancyTracer(),
                            boundary_conditions = (; b=buoyancy_bcs))

N² = 1e-6
bᵢ(x, y, z) = N² * z
ϵ(x, y, z) = 1e-6 * randn()
set!(model, b=bᵢ, u=ϵ, w=ϵ)

# Set up simulation
simulation = Simulation(model, Δt=2minutes, stop_time=2hours)

progress(sim) = @printf("Iter: %d, time: %s \n", iteration(sim), prettytime(sim))
simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

run!(simulation)


text: #meaningless comment for test purposes


text:
using Pkg
println("Activating environment in $(pwd())...")
Pkg.activate(".")
println("Installing packages...")
flush(stdout)
Pkg.instantiate()
Pkg.precompile()

println("Done!")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DiskAccessProperties <: SwaggerModel
    privateEndpointConnections::Any # spec type: Union{ Nothing, Vector{PrivateEndpointConnection} } # spec name: privateEndpointConnections
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    timeCreated::Any # spec type: Union{ Nothing, DateTime } # spec name: timeCreated

    function DiskAccessProperties(;privateEndpointConnections=nothing, provisioningState=nothing, timeCreated=nothing)
        o = new()
        validate_property(DiskAccessProperties, Symbol("privateEndpointConnections"), privateEndpointConnections)
        setfield!(o, Symbol("privateEndpointConnections"), privateEndpointConnections)
        validate_property(DiskAccessProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(DiskAccessProperties, Symbol("timeCreated"), timeCreated)
        setfield!(o, Symbol("timeCreated"), timeCreated)
        o
    end
end # type DiskAccessProperties

const _property_map_DiskAccessProperties = Dict{Symbol,Symbol}(Symbol("privateEndpointConnections")=>Symbol("privateEndpointConnections"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("timeCreated")=>Symbol("timeCreated"))
const _property_types_DiskAccessProperties = Dict{Symbol,String}(Symbol("privateEndpointConnections")=>"Vector{PrivateEndpointConnection}", Symbol("provisioningState")=>"String", Symbol("timeCreated")=>"DateTime")
Base.propertynames(::Type{ DiskAccessProperties }) = collect(keys(_property_map_DiskAccessProperties))
Swagger.property_type(::Type{ DiskAccessProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DiskAccessProperties[name]))}
Swagger.field_name(::Type{ DiskAccessProperties }, property_name::Symbol) =  _property_map_DiskAccessProperties[property_name]

function check_required(o::DiskAccessProperties)
    true
end

function validate_property(::Type{ DiskAccessProperties }, name::Symbol, val)
end


text:
#######################################
#
# dat file
#
#######################################

"""
Read dat files

#### Arguments
* `fname`: Name or path for the dat file

#### Returns
* `x`: Range of x values
* `y`: Range of y values
* `z`: Range of z values
* `complete_data`: Array (x × y × z x t)
* `sample_times`

#### References
File specs were taken from [fieldtrip](https://github.com/fieldtrip/fieldtrip/blob/1cabb512c46cc70e5b734776f20cdc3c181243bd/external/besa/readBESAimage.m)
"""
function read_dat(fname::AbstractString)
    @info("Reading dat file = $fname")

    read_dat(open(fname, "r"))
end

function read_dat(fid::IO)

    if isa(fid, IOBuffer)
        fid.ptr = 1
    end

    # Ensure we are reading version 2
    versionAbstractString = match(r"(\S+):(\d.\d)", readline(fid))
    version = float(versionAbstractString.captures[2])
    @debug("Version = $version")

    # Use @assert here?
    if version != 2
        Logging.warn("Unknown dat file version")
        return
    end

    # Header info
    readline(fid) # Empty line
    data_file = readline(fid)
    condition = readline(fid)
    typeline  = readline(fid)

    # Types of data that can be stored
    if search(typeline, "Method") != 0:-1  # TODO: change to imatch
        @debug("File type is Method")

        image_type = typeline[21:end]
        image_mode = "Time"
        regularization = readline(fid)[21:end-1]

        #=Latency=#  # TODO: Fix for latencies. See fieldtrip

        # Units
        units          = readline(fid)[3:end-1]

        @debug("Regularisation = $regularization")
        @debug("Units = $units")
    elseif search(typeline, "MSBF") != 0:-1

        image_mode = "Single Time"
        image_type = "Multiple Source Beamformer"
        units = condition[3:end-1]
        regularization = "None"

        Logging.warn("MSBF type under development")
    elseif search(typeline, "MSPS") != 0:-1
        Logging.warn("MSPS type not implemented yet")
    elseif search(typeline, "Sens") != 0:-1
        Logging.warn("Sens type not implemented yet")
    else
        Logging.warn("Unknown type")
    end

    readline(fid) # Empty line
    description = readline(fid)

    # Read in the dimensions
    regexp = r"[X-Z]:\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+)"
    xrange = match(regexp, readline(fid))
    x = range(float(xrange.captures[1]), stop=float(xrange.captures[2]), length=parse(Int, xrange.captures[3]))
    yrange = match(regexp, readline(fid))
    y = range(float(yrange.captures[1]), stop=float(yrange.captures[2]), length=parse(Int, yrange.captures[3]))
    zrange = match(regexp, readline(fid))
    z = range(float(zrange.captures[1]), stop=float(zrange.captures[2]), length=parse(Int, zrange.captures[3]))

    empty       = readline(fid)

    # Variables to fill
    t = 1
    complete_data = Array{Float64}((length(x), length(y), length(z), t))
    sample_times  = Float64[]

    description = readline(fid)
    if search(description, "Sample") != 0:-1

        #
        # 4D file
        #

        s = match(r"Sample \d+, (-?\d+.\d+) ms", description)
        push!(sample_times, float(s.captures[1]))

        file_still_going = true
        while file_still_going
            for zind = 1:length(z)
                readline(fid)           # Z: z

                for yind = 1:length(y)
                    d = readline(fid)       # values
                    m = matchall(r"(-?\d+.\d+)", d)
                    complete_data[:, yind, zind, t] = float(m)
                end

                readline(fid)           # blank or dashed
            end

            if eof(fid)
                file_still_going = false
            else
                t += 1
                s = readline(fid)               # Sample n, t.tt ms
                s = match(r"Sample \d+, (-?\d+.\d+) ms", s)
                push!(sample_times, float(s.captures[1]))

                # There is no nice way to grow a multidimensional array
                temp = complete_data
                complete_data = Array{Float64}((length(x), length(y), length(z), t))
                complete_data[:,:,:,1:t-1] = temp
            end
        end
    else

        #
        # 3D file
        #

        file_still_going = true
        idx = 1
        while file_still_going
            for zind = 1:length(z)

                for yind = 1:length(y)
                    d = readline(fid)       # values
                    m = matchall(r"(-?\d+.\d+)", d)
                    complete_data[:, yind, zind, t] = float(m)
                end

                d = readline(fid)           # blank or dashed
                d = readline(fid)           # blank or dashed
            end

            if eof(fid)
                file_still_going = false
            else
                println("File should have finished")
            end
        end

        sample_times = [0]
    end

    close(fid)

    return x, y, z, complete_data, sample_times
end

"""
Write dat file
"""
function write_dat(fname::AbstractString,
X::AbstractVector, Y::AbstractVector, Z::AbstractVector,
S::Array{DataT,4}, T::AbstractVector;
data_file::AbstractString="NA", condition::AbstractString="NA", method::AbstractString="NA", regularization::AbstractString="NA",
units::AbstractString="NA") where DataT <: AbstractFloat

    if size(S,1) != length(X); Logging.warn("Data and x sizes do not match"); end
    if size(S,2) != length(Y); Logging.warn("Data and y sizes do not match"); end
    if size(S,3) != length(Z); Logging.warn("Data and z sizes do not match"); end
    if size(S,4) != length(T); Logging.warn("Data and t sizes do not match"); end

    @info("Saving dat to $fname")

    open(fname, "w") do fid
        @printf(fid, "BESA_SA_IMAGE:2.0\n")
        @printf(fid, "\n")
        @printf(fid, "Data file:          %s\n", data_file)
        @printf(fid, "Condition:          %s\n", condition)
        @printf(fid, "Method:             %s\n", method)
        @printf(fid, "Regularization:     %s\n", regularization)
        @printf(fid, "  %s\n",                   units)

        @printf(fid, "\n")
        @printf(fid, "Grid dimensions ([min] [max] [nr of locations]):\n")
        @printf(fid, "X:  %2.6f %2.6f %d\n", minimum(X), maximum(X), length(X))
        @printf(fid, "Y:  %2.6f %2.6f %d\n", minimum(Y), maximum(Y), length(Y))
        @printf(fid, "Z:  %2.6f %2.6f %d\n", minimum(Z), maximum(Z), length(Z))
        @printf(fid, "==============================================================================================\n")

        for t = 1:size(S, 4)
            @printf(fid, "Sample %d, %1.2f ms\n", t-1, T[t])
            for z = 1:size(S, 3)
                @printf(fid, "Z: %d\n", z-1)
                for y = 1:size(S,2)
                    for x = 1:size(S,1)
                        @printf(fid, "%2.10f ", S[x,y,z,t])
                    end
                    @printf(fid, "\n")
                end
                @printf(fid, "\n")
            end
        end
        @debug("File successfully written")
    end
end


text:
# This file comtains some useful functions to import logfiles for use in either julia or matlab.
# In julia, the data object from orcalog2mat can be used together with the function getdata
# In matlab, the function csv2mat can be used to convert the textfile to a .mat-file which is very fast to read

using DelimitedFiles

"""
    csv2dict(filename, destination="log.mat"; startline=0, writeNames = true, subsample = 1, separator = ',')

convert a csv table file to a dict
"""
function csv2dict(filename, destination="log.mat"; startline=0, writeNames = true, subsample = 1, separator = ',')

    # Read the csv file into a matrix
    println("Reading file $filename of size $(filesize(filename)/1_000_000) MB")
    df,n = readdlm(filename, separator, header=true, skipstart=startline)
    println("Preprocessing")
    (rows, cols) = size(df)
    # replace NA values with 0
    for  j = 1:cols, i = 1:rows
        ismissing(df[i,j]) && (df[i,j] = 0)
    end

    # Remove strange characters from variable names
    replacechars = [' ', ':', '.', '[', ']']
    stripchars = ['_']
    for (i,symbol) in enumerate(n)
        n[i] = (replace(string(symbol),replacechars => "_"))
        n[i] = (strip(n[i],stripchars))
        writeNames && println(n[i])
    end

    # Write the results to a .mat-file
    written = text = 0
    println("Creating save dict")
    d = Dict{String, Vector{<:Real}}()
    for i in 1:cols
        if !isa(df[1,i],Number)
            text += 1
            continue
        end
        written += 1
        a = replace(convert(Array,df[1:subsample:end,i]), missing=>0)
        d[string(n[i])] = a
    end
    return d
end

"""`findData(pattern, data::Dict)`"""
function findData(pattern, data::AbstractDict)
    inds = falses(length(keys(data)))
    pattern isa Regex || (pattern = Regex(pattern))
    for (i,key) in enumerate(keys(data))
        if match(pattern, key) != nothing
            inds[i] = true
        end
    end
    return inds
end

"""
    `array = getdata(pattern, data::Dict, ds=1; removeNaN = false)`

pattern is a string representation of a regex used to match the data
e.g. \"posRawAbs\"
"""
function getdata(pattern, data::AbstractDict, ds=1; removeNaN = false)
    inds = findData(pattern, data)
    if sum(inds) < 1
        @warn("No data found using $pattern, available headers are")
        foreach(println, keys(data))
        return
    end
    names = collect(keys(data))[inds]
    function lt(x,y) # String comparison function, this soring places 2 in front of 10
        if length(x) < length(y)
            return true
        elseif length(x) > length(y)
            return false
        else
            return x < y
        end
    end
    sort!(names,lt=lt)
    retmat = Matrix{Float64}(undef, ceil(Int,length(data[names[1]])/ds),length(names))
    for (j,key) in enumerate(names)
        for (i,k)= enumerate(1:ds:length(data[key]))
            retmat[i,j] = data[key][k]
        end
    end

    if removeNaN
        retmat[any(isnan(retmat),2),:] = []
    end
    return retmat
end

max_forks_repo_name: baggepinnen/Robotlib


text:
module GridapFSI

using TimerOutputs

export main

include("WeakForms/WeakForms.jl")
include("Drivers/Drivers.jl")

using GridapFSI.Drivers

function main(;problemName::String="elasticFlag",kwargs...)

    reset_timer!()

    # Setup problem
    problem = get_problem(problemName,kwargs)

    # Execute driver
    outputs = execute(problem; kwargs...)

    print_timer()
    println()

    return outputs

end

end # module


text:
# Parts of this code were taken / derived from NetworkX. See LICENSE for
# licensing details.

using Base.Threads

"""
    pagerank(g, α=0.85, n=100, ϵ=1.0e-6)

Calculate the [PageRank](https://en.wikipedia.org/wiki/PageRank) of the
graph `g` parameterized by damping factor `α`, number of iterations
`n`, and convergence threshold `ϵ`. Return a vector representing the
centrality calculated for each node in `g`, or an error if convergence
is not reached within `n` iterations.
"""
function pagerank(
    g::AbstractGraph{U},
    α=0.85,
    n::Integer=100,
    ϵ=1.0e-6
    ) where U <: Integer
    α_div_outdegree = Vector{Float64}(undef,nv(g))
    dangling_nodes = Vector{U}()
    for v in vertices(g)
        if outdegree(g, v) == 0
            push!(dangling_nodes, v)
        end
        α_div_outdegree[v] = (α/outdegree(g, v))
    end
    N = Int(nv(g))
    # solution vector and temporary vector
    x = fill(1.0 / N, N)
    xlast = copy(x)
    for _ in 1:n
        dangling_sum = 0.0
        for v in dangling_nodes
            dangling_sum += x[v]
        end
        # flow from teleprotation
        for v in vertices(g)
            xlast[v] = (1 - α + α * dangling_sum) * (1.0 / N)
        end
        # flow from edges

        for v in vertices(g)
            for u in inneighbors(g, v)
                xlast[v] += (x[u] * α_div_outdegree[u])
            end
        end
        # l1 change in solution convergence criterion
        err = 0.0
        for v in vertices(g)
            err += abs(xlast[v] - x[v])
            x[v] = xlast[v]
        end
        if (err < N * ϵ)
            return x
        end
    end
    error("Pagerank did not converge after $n iterations.") # TODO 0.7: change to InexactError with appropriate msg.
end


text:
function _build_reified_constraint(
    _error::Function,
    variable::JuMP.AbstractVariableRef,
    constraint::JuMP.ScalarConstraint,
    ::Type{<:CS.ReifiedSet{A}},
) where {A}
    S = typeof(JuMP.moi_set(constraint))
    set = ReifiedSet{A,S}(JuMP.moi_set(constraint), 2)
    return JuMP.VectorConstraint([variable, JuMP.jump_function(constraint)], set)
end

function _build_reified_constraint(
    _error::Function,
    variable::JuMP.AbstractVariableRef,
    jump_constraint::JuMP.VectorConstraint,
    ::Type{<:CS.ReifiedSet{A}},
) where {A}
    S = typeof(jump_constraint.set)
    set = CS.ReifiedSet{A,S}(jump_constraint.set, 1 + length(jump_constraint.func))
    if jump_constraint.func isa Vector{VariableRef}
        vov = JuMP.VariableRef[variable]
    else
        vov = JuMP.AffExpr[variable]
    end
    append!(vov, jump_constraint.func)
    return JuMP.VectorConstraint(vov, set)
end

function _reified_variable_set(::Function, variable::Symbol)
    return variable, ReifiedSet{MOI.ACTIVATE_ON_ONE}
end

function _reified_variable_set(_error::Function, expr::Expr)
    if expr.args[1] == :¬ || expr.args[1] == :!
        if length(expr.args) != 2
            _error("Invalid binary variable expression `$(expr)` for reified constraint.")
        end
        return expr.args[2], ReifiedSet{MOI.ACTIVATE_ON_ZERO}
    else
        return expr, ReifiedSet{MOI.ACTIVATE_ON_ONE}
    end
end

function JuMP.parse_constraint_head(_error::Function, ::Val{:(:=)}, lhs, rhs)
    variable, S = _reified_variable_set(_error, lhs)
    if !JuMP.isexpr(rhs, :braces) || length(rhs.args) != 1
        _error("Invalid right-hand side `$(rhs)` of reified constraint. Expected constraint surrounded by `{` and `}`.")
    end
    rhs_con = rhs.args[1]
    rhs_vectorized, rhs_parsecode, rhs_buildcall =
        JuMP.parse_constraint_expr(_error, rhs_con)

    # TODO implement vectorized version
    vectorized = false
    if rhs_vectorized
        _error("`$(rhs)` should be non vectorized. There is currently no vectorized support for reified constraints. Please open an issue at ConstraintSolver.jl")
    end

    buildcall = :($(esc(:(CS._build_reified_constraint)))(
        $_error,
        $(esc(variable)),
        $rhs_buildcall,
        $S,
    ))
    return vectorized, rhs_parsecode, buildcall
end


text:
export CalibratedArraySwitches, CalibratedArray
export validate_ca

function CalibratedArraySwitches(objId :: ObjectId,
    defaultValueM :: Array{T1, N}, lbM :: Array{T1, N}, ubM :: Array{T1, N},
    isCalInM) where {T1, N}

    isCalM = Array{Bool, N}(isCalInM);
    @assert size(defaultValueM) == size(lbM) == size(ubM) == size(isCalM);

    doCal = true;
    valueV = defaultValueM[isCalM];
    if all(.!isCalM)
        doCal = false;
    end
    p = Param(:calValueV, "Calibrated values", "calValueV",
        valueV, valueV, lbM[isCalM], ubM[isCalM], doCal);
    pvec = ParamVector(objId, [p]);
    return CalibratedArraySwitches(pvec, defaultValueM, lbM, ubM, isCalM)
end

"""
	$(SIGNATURES)

Array with a mix of calibrated and fixed elements.
"""
function CalibratedArray(switches :: CalibratedArraySwitches{T1, N}) where {T1, N}
    calValueV = switches.defaultValueM[switches.isCalM];
    valueM = copy(switches.defaultValueM);
    objId = get_object_id(switches);
    return CalibratedArray(objId, switches, calValueV, valueM)
end

function make_test_calibrated_array_switches(N :: Integer; allFixed = false)
    rng = MersenneTwister(23);
    objId = ObjectId(:test);
    sizeV = 2 : (1 + N);
    sz = (sizeV..., );

    T1 = Float64;
    defaultValueM = randn(rng, T1, sz);
    lbM = defaultValueM .- one(T1);
    ubM = defaultValueM .+ one(T1);
    if allFixed
        isCalM = trues(size(lbM));
    else
        isCalM = defaultValueM .> 0.0;
    end

    switches = CalibratedArraySwitches(objId, defaultValueM, lbM, ubM, isCalM);
    return switches
end

## -------------  Basics

Lazy.@forward CalibratedArray.switches (
    Base.size
);

has_pvector(::CalibratedArraySwitches{T1, N}) where {T1, N} = true;

# ModelObjectsLH.get_object_id(switches :: CalibratedArraySwitches{T1, N}) where {T1, N} = get_object_id(switches.pvec);
# get_pvector(switches :: CalibratedArraySwitches{T1, N}) where {T1, N} =
#     switches.pvec;
Base.size(switches :: CalibratedArraySwitches{T1, N}) where {T1, N} =
    size(switches.defaultValueM);

## --------  Validate

function validate_ca(ca :: CalibratedArray{T1, N}) where {T1, N}
    switches = ca.switches;
    isValid = true;
    isValid = isValid  && (size(switches.defaultValueM) == size(switches.lbM) == size(switches.ubM) == size(switches.isCalM));
    isValid = isValid  &&  check_values(ca);
    return isValid
end

function check_values(ca :: CalibratedArray{T1, N}) where {T1, N}
    isValid = check_calibrated_values(ca)  &&  check_fixed_values(ca);
    return isValid
end

function check_calibrated_values(ca :: CalibratedArray{T1, N}) where {T1, N}
    idxV = map_indices(ca);
    isempty(idxV)  &&  return true;

    calValueV = ca.valueM[idxV];
    isValid = isapprox(calValueV, ca.calValueV);

    pvec = get_pvector(ca);
    isValid = isValid  &&  isapprox(calValueV, param_value(pvec, :calValueV));
    return isValid
end

function check_fixed_values(ca :: CalibratedArray{T1, N}) where {T1, N}
    valueM = values(ca);
    isValid = true;
    for (j, isCal) in enumerate(ca.switches.isCalM)
        if !isCal
            isValid = isValid && (valueM[j] ≈ ca.switches.defaultValueM[j]);
        end
    end
    return isValid
end

## ------------  Retrieve

# change +++: retrieve with array interface; without allocating
function values(ca :: CalibratedArray{T1, N}) where {T1, N}
    update_values!(ca);
    return ca.valueM;
end

# Update calibrated array values in object from calibrated vector values.
function update_values!(ca :: CalibratedArray{T1, N}) where {T1, N}
    ca.valueM .= ca.switches.defaultValueM;
    idxV = map_indices(ca);
    for (j, idx) in enumerate(idxV)
        ca.valueM[idx] = ca.calValueV[j];
    end
end

# Linear index into arrays for each calibrated value.
function map_indices(ca :: CalibratedArray{T1, N}) where {T1, N}
    nCal = sum(ca.switches.isCalM);
    if nCal > 0
        idxV = zeros(Int, nCal);
        j2 = 0;
        for (j, isCal) in enumerate(ca.switches.isCalM)
            if isCal
                j2 += 1;
                idxV[j2] = j;
            end
        end
        @assert j2 == length(ca.calValueV);
    else
        idxV = Vector{Int}();
    end
    return idxV
end

# ---------------


text:
"""
    rlap(A)

Compute renormalized graph Laplacian for adjacency matrix A.
"""
function rlap(A)
    W = diagm(A * ones(size(A, 1)))
    L̃ = W \ (W \ A)'
    W̃ = diagm(L̃ * ones(size(L̃, 1)))

    println("Laplacian sparsity = $(sparsity(W̃ \ (W̃ - L̃)))")
    return sparse(W̃ \ (W̃ - L̃))
end

"""
    gaussian_similarity(X, σ=1)

Compute similarity using Gaussian kernel. Enforces sparsity by setting the similarity of
distant points to zero.
"""
function gaussian_similarity(X, σ=1)
    𝒦 = pairwise(Euclidean(), X, dims=2)
    n = size(X, 2)
    for i in 1:n
        for j in i+1:n
            𝒦[j, i] = 𝒦[j, i] <= 3σ ? exp(-(𝒦[j,i]/σ)^2) : 0
            𝒦[i, j] = 𝒦[j, i]
        end
    end
    𝒦[diagind(𝒦)] .= 1

    s = sparsity(𝒦)
    if s < 0.8
        println("Warning. sparsity = $s")
    end

    return 𝒦
end


text:
include("common.jl")
@testset "Dependent delays" begin
    alg = MethodOfSteps(BS3())
    dde_int = init(prob_dde_1delay, alg)
    sol = solve!(dde_int)

    @testset "reference" begin
        @test sol.errors[:l∞] < 3.0e-5
        @test sol.errors[:final] < 2.1e-5
        @test sol.errors[:l2] < 1.2e-5
    end

    @testset "constant lag function" begin
        # constant delay specified as lag function
        prob2 = DDEProblem(DiffEqProblemLibrary.DDEProblemLibrary.ff_1delay,
                           [1.0],
                           (p, t) -> [0.0], (0., 10.),
                           dependent_lags = [(u,p,t) -> 1])
        dde_int2 = init(prob2, alg)
        sol2 = solve!(dde_int2)

        @test dde_int.tracked_discontinuities == dde_int2.tracked_discontinuities

        # worse than results above with constant delays specified as scalars
        @test sol2.errors[:l∞] < 4.2e-5
        @test sol2.errors[:final] < 2.2e-5
        @test sol2.errors[:l2] < 1.7e-5

        # simple convergence tests
        @testset "convergence" begin
            sol3 = solve(prob2, alg, abstol=1e-9, reltol=1e-6)

            @test sol3.errors[:l∞] < 7.5e-8
            @test sol3.errors[:final] < 4.6e-8
            @test sol3.errors[:l2] < 3.9e-8

            sol4 = solve(prob2, alg, abstol=1e-13, reltol=1e-13)

            @test sol4.errors[:l∞] < 6.9e-11
            @test sol4.errors[:final] < 1.1e-11
            @test sol4.errors[:l2] < 6.8e-12 # 6.7e-12, relaxed for Win32
        end
    end

    # without any delays specified is worse
    @testset "without delays" begin
        prob2 = DDEProblem(DiffEqProblemLibrary.DDEProblemLibrary.ff_1delay,
                           [1.0],
                           (p, t) -> [0.0], (0., 10.))
        sol2 = solve(prob2, alg)

        @test sol2.errors[:l∞] > 1e-3
        @test sol2.errors[:final] > 4e-5
        @test sol2.errors[:l2] > 7e-4
    end
end


text:
"""
Momentum-based utterance selector, after Stadler et al. (2016).

# Reference

Stadler, K., Blythe, R. A., Smith, K. & Kirby, S. (2016) Momentum in language
change: a model of self-actuating s-shaped curves. *Language Dynamics and Change*,
6, 171–198. <https://doi.org/10.1163/22105832-00602005>
"""
mutable struct MomentumSelector <: AbstractUtteranceSelector
  replicator::BinaryVariable{Float64}
  gamma_replicator::BinaryVariable{Float64}
  alpha::Float64
  gamma::Float64
  b::Float64
  T::Int
  mmax::Float64
  meta::Any
end

"""
    MomentumSelector(alpha::Float64, gamma::Float64, b::Float64, T::Int, x0::Float64)

Construct a `MomentumSelector` with learning rate `alpha`, short-term memory
smoothing factor `gamma`, momentum bias `b`, number of utterances per interaction `T`,
and initial frequency `x0`.
"""
function MomentumSelector(alpha::Float64, gamma::Float64, b::Float64, T::Int, x0::Float64)
  tmmax = log(alpha/gamma)/(alpha - gamma)
  mmax = exp(-alpha*tmmax) - exp(-gamma*tmmax)
  MomentumSelector(BinaryVariable{Float64}(x0), BinaryVariable{Float64}(x0), alpha, gamma, b, T, mmax, nothing)
end

"""
    act!(x::MomentumSelector, y::MomentumSelector)

Let one `MomentumSelector` act upon another, i.e. `x` speaks and `y` listens.
"""
function act!(x::MomentumSelector, y::MomentumSelector)
  normalized_momentum = (y.gamma_replicator.x - y.replicator.x)/y.mmax
  observed_frequency = rand(Distributions.Binomial(x.T, x.replicator.x))/x.T

  transformed_frequency = 0.0

  if observed_frequency == 1.0
    transformed_frequency = 1.0
  elseif observed_frequency != 0.0
    transformed_frequency = observed_frequency + y.b*normalized_momentum
  end

  if transformed_frequency > 1.0
    transformed_frequency = 1.0
  elseif transformed_frequency < 0.0
    transformed_frequency = 0.0
  end

  y.replicator.x = y.alpha*transformed_frequency + (1-y.alpha)*y.replicator.x
  y.gamma_replicator.x = y.gamma*transformed_frequency + (1-y.gamma)*y.gamma_replicator.x
end


text:
@testset "Combined Manifolds" begin
using Random: seed!

M1 = Euclidean(3)
x1 = RnPoint([1.,2.,0.])
y1 = RnPoint([2.,3.,0.])
ξ1 = log(M1,x1,y1)
η1 = RnTVector([1.,0.,0.])

M2 = Sphere(2)
x2 = SnPoint([1.,0.,0.])
y2 = SnPoint([0.,0.,1.])
ξ2 = log(M2,x2,y2)
η2 = SnTVector(1/sqrt(2)*[0., 1.,1.])

#
# Prod
M = Product([M1, M2])
x = ProdPoint([x1,x2])
y = ProdPoint([y1,y2])
ξ = ProdTVector([ξ1,ξ2])
η = ProdTVector([η1, η2])

@test getValue(x) == [x1,x2]
@test getValue(ξ) == [ξ1, ξ2]

@test distance(M,x,y) == sqrt(distance(M1,x1,y1)^2 + distance(M2,x2,y2).^2)
@test dot(M,x,ξ,η) == dot(M1,x1,ξ1,η1) + dot(M2,x2,ξ2,η2)
@test distance( M, exp(M,x,ξ), ProdPoint([exp(M1,x1,ξ1),exp(M2,x2,ξ2)]) ) == 0
@test norm( M, x, log(M,x,y) - ProdTVector([log(M1,x1,y1), log(M2,x2,y2)]) ) ≈ 0
@test manifoldDimension(x) == 5
@test manifoldDimension(M) == 5
@test norm(M,y, parallelTransport(M,x,y,ξ) - ProdTVector([
    parallelTransport(M1,x1,y1,ξ1), parallelTransport(M2,x2,y2,ξ2) ]) ) ≈ 0

@test typeofMPoint(ξ) == ProdPoint{Array{MPoint,1}}
@test typeofTVector(x) == ProdTVector{Array{TVector,1}}
@test typeofMPoint(ProdTVector{Array{TVector,1}}) == ProdPoint{Array{MPoint,1}}
@test typeofTVector(ProdPoint{Array{MPoint,1}}) == ProdTVector{Array{TVector,1}}
@test validateMPoint(M, randomMPoint(M) )
@test validateTVector(M,x, randomTVector(M,x) )

@test_throws DomainError validateMPoint(Product([M1]),x)
@test_throws DomainError validateTVector(M,x,ProdTVector([ξ1]))

@test typicalDistance(M) ≈ sqrt( 2* (typicalDistance(M1).^2 + typicalDistance(M2).^2) )

@test zeroTVector(M,x) == ProdTVector([zeroTVector(M1,x1), zeroTVector(M2,x2)])

@test "$M" == "The Product Manifold of [ $(M1), $(M2) ]"
@test "$x" == "Prod[ $(x1), $(x2) ]"
@test "$ξ" == "ProdT[ $(ξ1), $(ξ2) ]"
#
# Power
M = Power(M1,2)
@test M.powerSize == (2,)
@test Power(M1,(3,2)).powerSize == (3,2)

x = PowPoint([x1,x1])
y = PowPoint([y1,y1])
ξ = PowTVector([ξ1,ξ1])
η = PowTVector([η1,η1])
#
@test size(x) == size(x.value)
@test distance(M, copy(x), x) == 0
@test norm(M,x, ξ - copy(ξ)) == 0
@test ndims(x) == ndims(x.value)
@test getValue(x) == [x1,x1]
@test getValue(ξ) == [ξ1,ξ1]
#
x[2] = RnPoint([1., 2., 0.])
@test getValue(x[2]) == [1., 2., 0.]
ξ[2] = RnTVector([1., 2., 0.])
@test getValue(ξ[2]) == [1., 2., 0.]
@test PowPoint([x1 x1; x1 x1])[1,2] == x1
x[2] = x1
ξ[2] = ξ1
xG = PowPoint([x1 x1;y1 y1])
@test xG[1:2,2].value == [x1,y1]
@test ndims(xG) == 2
ξG = PowTVector([ξ1 ξ1;η1 η1])
@test ξG[1:2,2].value == [ξ1,η1]
@test ndims(ξG) == 2
@test size(repeat(x, outer=[2,1])) == (4,1)
@test size(repeat(x, 2,1)) == (4,1)
@test size(repeat(x, inner=[1,2])) == (2,2)
@test size(repeat(ξ, outer=[2,1])) == (4,1)
@test size(repeat(ξ, 2,1)) == (4,1)
@test size(repeat(ξ, inner=[1,2])) == (2,2)
#
@test getValue.(getValue( adjointJacobiField(M, x, y, 0.5, ξ))) == getValue.(
        adjointJacobiField.(Ref(M1), [x1,x1], [y1,y1], 0.5, [ξ1,ξ1]) )
@test getValue.(getValue( jacobiField(M, x, y, 0.5, ξ))) ==
        getValue.(jacobiField.(Ref(M1), [x1,x1], [y1,y1], 0.5, [ξ1,ξ1]) )
#
@test distance(M,x,y) ≈ sqrt(2)*distance(M1,x1,y1)
@test dot(M,x,ξ,η) ≈ 2*dot(M1,x1,ξ1,η1)
#
@test getValue(exp(M, x, ξ)) == [exp(M1, x1, ξ1), exp(M1, x1, ξ1)]
@test getValue(log(M, x, y)) == [log(M1, x1, y1), log(M1, x1, y1)]
#
@test manifoldDimension(x) == 2*manifoldDimension(x1)
@test manifoldDimension(M) == 2*manifoldDimension(M1)
@test norm(M,x,ξ) ≈ sqrt(2)*norm(M1,x1,ξ1)
@test ndims(x) == 1
@test ndims(ξ) == 1
@test getValue( parallelTransport(M,x,y,ξ) ) == [parallelTransport(M1, x1, y1, ξ1),parallelTransport(M1, x1, y1, ξ1)]
@test validateMPoint(M,randomMPoint(M))
@test validateTVector(M,x,randomTVector(M,x))
@test validateTVector(M,x,randomTVector(M,x,:Gaussian))

@test typeofMPoint(ξ) == PowPoint{RnPoint{Float64},1}
@test typeofTVector(x) == PowTVector{RnTVector{Float64},1}
@test typeofTVector(PowPoint{RnPoint{Float64},1}) == PowTVector{RnTVector{Float64},1}
@test typeofMPoint(PowTVector{RnTVector{Float64},1}) == PowPoint{RnPoint{Float64},1}

@test_throws DomainError validateMPoint(Power(M,3),x)
@test_throws DomainError validateTVector(M,x,PowTVector([ξ1]))

@test typicalDistance(M) ≈ sqrt( 2 ) * typicalDistance(M1)

@test zeroTVector(M,x) == PowTVector([zeroTVector(M1,x1), zeroTVector(M1,x1)])

a0 = RnTVector([0.,0.,0.]); a1 = RnTVector([1.,0.,0.]);
a2 = RnTVector([0.,1.,0.]); a3 = RnTVector([0.,0.,1.])
ΞT = [PowTVector([a1,a0]), PowTVector([a2,a0]), PowTVector([a3,a0]),
    PowTVector([a0,a1]), PowTVector([a0,a2]), PowTVector([a0,a3]) ]
κT = zeros(6)
Ξ,κ = tangentONB(M,x,y)

@test sum( abs.(κ-κT)) == 0
@test sum( norm.(Ref(M), Ref(x),  ΞT .- Ξ ) ) == 0
# show
@test "$x" == "Pow[Rn([1.0, 2.0, 0.0]), Rn([1.0, 2.0, 0.0])]"
@test "$ξ" == "PowT[RnT([1.0, 1.0, 0.0]), RnT([1.0, 1.0, 0.0])]"
@test "$M" == "The Power Manifold of The 3-dimensional Euclidean space of size (2,)."
#
# Tangent Bundle
M = TangentBundle(M1)
@test getBase(M) == M1
x = TBPoint(x1,ξ1)
y = TBPoint(y1,η1)
xT = TBPoint( (x1,ξ1) )
@test x == xT
@test getValue(x) == (x1,ξ1)
xE = MPointE(x)
@test getBase(x) == x1
@test getTangent(x) == ξ1
@test getBase(xE) == MPointE(x1)
@test getTangent(xE) == TVectorE(ξ1,x1)

ξ = TBTVector(ξ1,η1)
ξT = TBTVector( (ξ1,η1) )
@test ξ == ξT
@test getValue(ξ) == (ξ1,η1)
ξE = TVectorE(ξ,x)
@test getBase(ξ) == ξ1
@test getTangent(ξ) == η1
@test getBase(ξE) == TVectorE(ξ1,x1)
@test getTangent(ξE) == TVectorE(η1,x1)

@test 2*ξ == TBTVector(2*ξ1, 2*η1)
@test ξ+ξ == TBTVector(ξ1+ξ1,η1+η1)
@test ξ-ξ == TBTVector(ξ1-ξ1,η1-η1)
@test -ξ == TBTVector(-ξ1,-η1)
@test +ξ == TBTVector(ξ1,η1)

@test distance(M,x,y) ≈ sqrt(3)
@test norm(M,x,ξ) ≈ sqrt(dot(M,x,ξ,ξ))
@test norm(M,x,ξ) ≈ sqrt(norm(M1,x1,ξ1)^2 + norm(M1,x1,η1))
@test exp(M,x,ξ) == TBPoint(exp(M1,x1,ξ1), ξ1+η1)
@test log(M,x,y) ==  TBTVector(log(M1,x1,y1),ξ1-η1)
@test manifoldDimension(M) == 2*manifoldDimension(M1)
@test manifoldDimension(x) == 2*manifoldDimension(x1)
@test parallelTransport(M,x,y,ξ) == ξ
@test validateMPoint(M,randomMPoint(M))
@test validateTVector(M,x,randomTVector(M,x))

@test typeofMPoint(ξ) == TBPoint{RnPoint{Float64},RnTVector{Float64}}
@test typeofTVector(x) == TBTVector{RnTVector{Float64}}

Ξ, κ = tangentONB(M,x,y)
@test κ == zeros(6)
@test Ξ == [TBTVector(a1,a0), TBTVector(a2,a0), TBTVector(a3,a0), TBTVector(a0,a1), TBTVector(a0,a2), TBTVector(a0,a3)]

@test typicalDistance(M) == sqrt( typicalDistance(M1)^2 + manifoldDimension(M1))
@test zeroTVector(M,x) == TBTVector(a0,a0)

@test "$M" == "The Tangent bundle of <$(M.manifold)>"
@test "$x" == "TB(Rn([1.0, 2.0, 0.0]), RnT([1.0, 1.0, 0.0]))"
@test "$ξ" == "TBT(RnT([1.0, 1.0, 0.0]), RnT([1.0, 0.0, 0.0]))"
end


text:
#=
@doc """
    cmp(op1, op2)

Compare `op1` and `op2`. Return a positive value if `op1>op2`, zero if
`op1=op2`, and a negative value if `op1<op2`. Both `op1` and `op2` are
considered to their full own precision, which may differ. If one of
the operands is NaN, set the `erange` flag and return zero.

Note: These functions may be useful to distinguish the three possible
cases. If you need to distinguish two cases only, it is recommended to
use the predicate functions (e.g., `equal_p` for the equality); they
behave like the IEEE 754 comparisons, in particular when one or both
arguments are NaN. But only floating-point numbers can be compared
(you may need to do a conversion first).
"""
=#

#=
e has type mpft_exp_t which must be at least 32 bits and it is defined
in mpfr.h according to the value of _MPFR_EXP_FORMAT. I set to it
Cint.
=#
export cmp_2exp
for (specC, tJ, tC) in ((:ui, CulongMax,  Culong),
                        (:si, ClongMax,   Clong))
    @eval begin
        function cmp_2exp(op1::BigFloat, op2::$tJ, e)
            isnan(op1) && throw(DomainError("op1 is NaN"))
            ccall(($(string("mpfr_cmp_", specC, "_2exp")), :libmpfr), Int32,
                  (Ref{BigFloat}, $tC, Cint), op1, op2, Int32(e))
        end
    end
end

@doc """
    cmp_2exp(op1, op2, e)

Compare `op1` and `op2 × 2^e`: returns a positive value if `op1>op2 ×
2^e`, zero if `op1=op2 × 2^e` or a negative number if `op1<op2 ×
2^e`. Throws an exception if `op1` is NaN.
""" cmp_2exp

for fJ = (:number_p, :regular_p)
    @eval begin
        export $fJ
        ($fJ)(op::BigFloat) =
            ccall(($(string(:mpfr_, fJ)), :libmpfr), Int32, (Ref{BigFloat},), op) != 0
    end
end

@doc """
    number_p(op)

Return true if `op` is an ordinary number (i.e., neither NaN nor an
infinity). Return false otherwise.
""" number_p

@doc """
    regular_p(op)

Return true if `op` is a regular number (i.e., neither NaN, nor an
infinity nor zero). Return false otherwise.
""" regular_p

export cmpabs
"""
    cmpabs(op1, op2)

Compare `|op1|` and `|op2|`. Return a positive value if `|op1|>|op2|`,
zero if `|op1|=|op2|`, and a negative value if `|op1|<|op2|`. Raise a
DomainError exception if one of the operands is NaN.
"""
function cmpabs(op1::BigFloat, op2::BigFloat)
    isnan(op1) || isnan(op2) && throw(DomainError("Operands must not be NaNs"))
    ccall((:mpfr_cmpabs, "libmpfr"), Int32, (Ref{BigFloat}, Ref{BigFloat}),
          op1, op2)
end

export unordered_p
"""
    unordered_p(op1, op2)

Return true if `op1` or `op2` is a NaN (i.e., they cannot be
compared), false otherwise.
"""
unordered_p(op1::BigFloat, op2::BigFloat) =
    ccall((:mpfr_unordered_p, "libmpfr"), Int32, (Ref{BigFloat}, Ref{BigFloat}),
          op1, op2) != 0


text:
# TO DO: revise in light of how gigaword no longer a module
include(joinpath(@__DIR__, "../src/gigaword_64k.jl"))
using Test, Revise

const path_test_data = joinpath(@__DIR__, "test_data")
const path_test_afp = joinpath(path_test_data,"afp")
const path_test_cna = joinpath(path_test_data,"cna")

const path_mtdoc_with_text_tags = joinpath(path_test_data,"afp/empty_doc_with_empty_text")
const path_empty_file = joinpath(path_test_data, "afp/empty_file")
const path_onedoc = joinpath(path_test_data, "afp/one_doc")
const path_twodocs = joinpath(path_test_data, "afp/two_docs")

const path_big_file = joinpath(path_test_data, "afp/afp_eng_200304")

const path_test_output = joinpath(path_test_data, "test_output")

# ===============
# TESTS FOR UTILS
# ===============
#using gigaword_64k: is_wanted

@testset "is_wanted" begin
    @test is_wanted("abc") == true
    @test is_wanted("aBc") == true

    # Filters out punctuation
    @test is_wanted("!") == false

    # Filters out stuff that's just fancy numbers and not v long
    @test is_wanted("3.77") == false

    # note
    @test is_wanted("") == true # but won't usually happen
end

@testset "lzfoldl" begin
    @test lzfoldl((x,y)->x+y, list(1, 2, 3)) == 6

    @test lzfoldl((x,y)->x-y, list(1, 2, 3)) == -4
    # this is a good test case for foldl, since foldr would return -2
end

# ===================
# TESTS FOR MAIN FTNS
# ===================

@testset "acc_to_df" begin
    # empty case
    mt_acc = counter(String)
    df_from_mt_acc = acc_to_df(mt_acc)
    @test size(df_from_mt_acc) == (0, 2)

    # case with three entries
    acc_t = counter(String)
    acc_t["a"] = 1; acc_t["b"] = 2; acc_t["c"] = 7
    df_at = acc_to_df(acc_t)

    @test df_at[:word] == ["c", "b", "a"]
    @test df_at[:freq] == [7, 2, 1]
end

@testset "df_to_acc" begin
    mt_acc = counter(String)
    df_from_mt_acc = acc_to_df(mt_acc)
    @test df_to_acc(df_from_mt_acc) == mt_acc

    # case with three entries
    acc_t = counter(String)
    acc_t["a"] = 1; acc_t["b"] = 2; acc_t["c"] = 7
    df_at = acc_to_df(acc_t)

    @test df_to_acc(df_at) == acc_t
end

#using gigaword_64k: year_from_fnm
@test year_from_fnm("wpb_eng_201012") == "2010"

#using gigaword_64k: pick_elts
@test pick_elts([1,2], 3) == [1,2]
@test pick_elts([1,2, 3], 3) == [1,2, 3]

# for count_words_from_file

#using gigaword_64k: count_words_from_file, most_common

@testset "count_words_from_file: empty docs" begin
    empty_f = parsehtml(read(path_empty_file, String))
    @test count_words_from_file(empty_f) == counter(String)
    @test_throws BoundsError most_common(count_words_from_file(empty_f), 1)

    empty_f_with_txt_tags = parsehtml(read(path_mtdoc_with_text_tags, String))
    @test count_words_from_file(empty_f_with_txt_tags) == counter(String)
    @test_throws BoundsError most_common(count_words_from_file(empty_f_with_txt_tags), 1)
end

@testset "count_words_from_file: simple" begin
    doc = parsehtml(read(path_onedoc, String))
    top_8 = most_common(count_words_from_file(doc), 8)
    @test top_8 == [("trade" => 4),
                    ("and" => 2),
                    ("the" => 2),
                    ("japan" => 2),
                    ("on" => 2),
                    ("vaile" => 1),
                    ("eu" => 1),
                    ("global" => 1)]
end

@testset "process_part_of_tree" begin
    path_tmp_test_output = joinpath(path_test_output, "tmp")
    process_part_of_tree(path_test_cna, path_tmp_test_output, 4)

    jdf_cna_97_10 = @_ JDF.load(joinpath(path_tmp_test_output,
        "cna_eng_199710.jdf")) |> DataFrame
    acc_df_cna_97_10 = read_and_wc(joinpath(path_test_cna, "cna_eng_199710")) |> acc_to_df

    @test jdf_cna_97_10 == acc_df_cna_97_10

    # to eyeball: case where there's more than 4 files per year
    path_manyfiles = joinpath(path_test_data, "manyfiles")
    process_part_of_tree(path_manyfiles, path_tmp_test_output, 4)
    # note that this will delete everything that was in the tmp test output directory before outputting the stuff

    many_files_dir = FileTree(path_tmp_test_output)
    @test length(files(many_files_dir[glob"*.jdf"])) == 24

    rm(path_tmp_test_output, recursive=true, force=true)
end

#    jdf_to_df(joinpath(path_test_data_reds, "a", "df1") * ".jdf")

@testset "reduce_files" begin
    # 1. Make the test data
    path_test_data_reds = joinpath(path_test_data, "test_data_for_reductions")

    a_1 = counter(String)
    a_1["w1"] = 1; a_1["w2"] = 3
    df_a1 = acc_to_df(a_1)

    a_2 = counter(String)
    a_2["w1"] = 1; a_2["w2"] = 2; a_2["w3"] = 7
    df_a2 = acc_to_df(a_2)

    a_3 = counter(String)
    a_3["w1"] = 1
    df_a3 = acc_to_df(a_3)

    function dmapper(dirname)
        for (i, df) in enumerate([df_a1, df_a2, df_a3])
            fname = joinpath(path_test_data_reds, dirname, "df") * string(i) * ".jdf"
            savejdf(fname, df)
        end
    end
    map(dmapper, ["a", "b", "c"])

    # 2. Reduce the files

    ## Testing reduce_jdfs
    final_to_test = reduce_jdfs(path_test_data_reds)
    correct_final_acc = counter(String)
    correct_final_acc["w1"] = 9; correct_final_acc["w2"] = 15; correct_final_acc["w3"] = 21
    @test final_to_test == correct_final_acc

    ## Testing reduce_jdfs_and_save
    reduce_jdfs_and_save(path_test_data_reds, path_test_data; ns_tuple = (2, 11))

    top2_loaded_to_test = joinpath(path_test_data, "final_df_2.jdf") |> jdf_to_df |> df_to_acc
    top3_loaded_to_test = joinpath(path_test_data, "final_df_3.jdf") |> jdf_to_df |> df_to_acc

    @test acc_to_df(top2_loaded_to_test) == acc_to_df(nlargest(correct_final_acc, 2))
    @test top3_loaded_to_test == correct_final_acc

    rm(path_test_data_reds, recursive=true, force=true)

end

    # dirnames = name.(dirs(testft))
    # jdf_dirs = Set(filter(s->endswith(s, "jdf"), dirnames))
    # jdf_tree = filter(x->name(x) ∈ jdf_dirs, testft, dirs=true)

    # tree_of_dfs = FileTrees.load(testft[glob"*/*.jdf"], dirs = true) do dir
    #        @_ dir |> string(FileTrees.path(__)) |> JDF.load |> DataFrame
    #     end

    # 3. Test the reduction code
    # function red_inner_dir(subtree)
    #     reducevalues(subtree, init=nothing) do df1, df2
    #         if df1 isa DataFrame && df2 isa DataFrame
    #             @info "first branch"
    #             (acc1, acc2) = df_to_acc.((df1, df2))
    #             merge!(acc1, acc2)
    #             @info acc1
    #             return acc1

    #         elseif df1 isa DataFrame && !isa(df2, DataFrame)
    #             return df_to_acc(df1)

    #         elseif !isa(df1, DataFrame) && isa(df2, DataFrame)
    #             return df_to_acc(df2)

    #         else
    #             @info inputs to reducevalues were not dfs
    #             return nothing
    #         end
    #     end
    # # end
    # reduced_subdirs = mapsubtrees(red_inner_dir, testft, glob"*")

#=
Note on the many_files_dir test

julia> many_files_dir[glob"*.jdf"]
/gpfs/scratch/yh31/projects/gigaword_64k/test/test_data/test_output/tmp/
├─ cna_eng_200404.jdf/
│  ├─ freq
│  ├─ metadata.jls
│  └─ word
├─ cna_eng_200406.jdf/
│  ├─ freq
│  ├─ metadata.jls
│  └─ word
├─ cna_eng_200409.jdf/
│  ├─ freq
│  ├─ metadata.jls
│  └─ word
├─ cna_eng_200410.jdf/
│  ├─ freq
│  ├─ metadata.jls
│  └─ word
├─ cna_eng_200507.jdf/
│  ├─ freq
│  ├─ metadata.jls
│  └─ word
├─ cna_eng_200508.jdf/
│  ├─ freq
│  ├─ metadata.jls
│  └─ word
├─ cna_eng_200509.jdf/
│  ├─ freq
│  ├─ metadata.jls
│  └─ word
└─ cna_eng_200510.jdf/
   ├─ freq
   ├─ metadata.jls
   └─ word
=#

#=
jdf_cna_200403 = @_ JDF.load(joinpath(path_tmp_test_output,
        "cna_eng_199710.jdf")) |> DataFrame(__)
=#

# testing with mt-ish files
#=
tree_of_mts = FileTrees.load(data_tree[glob"empty*"]; lazy = true) do file
    @_ file |>
    string(FileTrees.path(__)) |>
    read(__, String) |>
    parsehtml(__) |>
    count_words_from_file(__)
end
=#


text:
using Base.Test
using Fuji
using Requests

import Fuji.post, Fuji.put, Fuji.patch, Fuji.delete

get(endpoint) = Requests.get("http://localhost:8000$endpoint")
post_req(endpoint) = Requests.post("http://localhost:8000$endpoint")
put_req(endpoint) = Requests.put("http://localhost:8000$endpoint")
patch_req(endpoint) = Requests.patch("http://localhost:8000$endpoint")
delete_req(endpoint) = Requests.delete("http://localhost:8000$endpoint")

@spawn Fuji.start()

@testset "simple route" begin
    @test size(Fuji.server.routes, 1) == 0

    route("/hi") do req, res
        "hi"
    end

    route("/histatus") do req, res
        204
    end

    @test size(Fuji.server.routes, 1) == 2
    @test Fuji.server.routes[1].endpoint == "/hi"
    @test Fuji.server.routes[2].endpoint == "/histatus"

    @test get("/hi").status == 200
    @test get("/histatus").status == 204
end

@testset "routes with other http methods" begin
    post("/post") do req, res
        "hi"
    end

    @test post_req("/post").status == 200

    put("/put") do req, res
        "hi"
    end

    @test put_req("/put").status == 200

    patch("/patch") do req, res
        "hi"
    end

    @test patch_req("/patch").status == 200

    delete("/delete") do req, res
        "hi"
    end

    @test delete_req("/delete").status == 200
end

@testset "adding two of the same route" begin
    @test size(Fuji.server.routes, 1) == 6

    route("/hi") do req, res
        "hello"
    end

    @test size(Fuji.server.routes, 1) == 6
    @test length(get("/hi").data) == 5
end

@testset "adding a route with the same endpoint but a different method" begin
    @test size(Fuji.server.routes, 1) == 6

    post("/hi") do req, res
        "hello"
    end

    @test size(Fuji.server.routes, 1) == 7
end

@testset "removing a route" begin
    @test size(Fuji.server.routes, 1) == 7

    r = route("/toremove") do req, res
        "hello"
    end

    @test size(Fuji.server.routes, 1) == 8

    unroute(r)

    @test size(Fuji.server.routes, 1) == 7
end

@testset "named parameter" begin
    route("/hello/:name") do req, res
        string("hello, ", req.params["name"], "!")
    end

    @test get("/hello/:name").status == 404
    @test get("/hello").status == 404
    @test get("/hello/jack/asdf").status == 404
    @test get("/hello/ja;ck").status == 404

    @test get("/hello/jack").status == 200
    @test length(get("/hello/jack").data) == 12

    @test get("/hello/j-a_ck0-1").status == 200
end

@testset "splat parameter" begin
    route("/hi/*") do req, res
        string("hi, ", req.splat[1], "!")
    end

    @test get("/hi/:name").status == 404
    @test get("/hi/").status == 404
    @test get("/hi/jack/asdf").status == 404
    @test get("/hi/ja;ck").status == 404

    @test get("/hi/jack").status == 200
    @test length(get("/hi/jack").data) == 9

    @test get("/hi/j-a_ck0-1").status == 200
end

@testset "before filter" begin
    x = 0

    before() do req, res
        x += 1
    end

    route("/before") do req, res
        string(x)
    end

    @test get("/before").data[1] == 0x31 # 1
    @test get("/something").status == 404

    @test x == 2
end

@testset "after filter" begin
    x = 0

    after() do req, res
        x += 1
    end

    route("/after") do req, res
        string(x)
    end

    @test get("/after").data[1] == 0x30 # 0
end

@testset "query parameters" begin
    route("/hello") do req, res
        # if "name" was in the request, return it, otherwise return an empty 204
        try
            req.query_params["name"]
        catch
            204
        end
    end

    @test length(get("/hello?name=jack").data) == 4
    @test get("/hello?name=jack").status == 200
    @test get("/hello?").status == 204
    @test get("/hello?test&hello").status == 204
    @test get("/hello?t=4&x=3&name=2").data[1] == 0x32 # 2
end

@testset "templating" begin
    route("/template") do req, res
        render_template("template.html", Dict("name" => "Jack"))
    end

    @test get("/template").status == 200

    f = open("template_expected.html")
    contents = readstring(f)
    close(f)

    @test contents == render_template("template.html", Dict("name" => "Jack"))
end


text:
function mark_innerrectangle_perimetr!(r::Robot)
    num_steps=fill(0,3) # - вектор-столбец из 3-х нулей
    for (i,side) in enumerate((Sud,West,Sud))
        num_steps[i]=moves!(r,side)
    end
    #УТВ: Робот - в Юго-западном углу внешней рамки

    side = find_border!(r,Ost,side)
    #УТВ: Робот - у западной границы внутренней перегородки

    mark_innerrectangle_perimetr!(r,side)
    #УТВ: Робот - снова у западной границы внутренней прямоугольной перегородки

    moves!(r,Sud)
    moves!(r,West)
    #УТВ: Робот - в Юго-западном улу внешней рамки

    for (i,side) in enumerate((Nord,Ost,Nord))
        moves!(r,side, num_steps[i])
    end
    #УТВ: Робот - в исходном положении
end

function mark_innerrectangle_perimetr!(r::Robot, side::HorizonSide)
    direction_of_movement, direction_to_border = get_directions(side)
    for i ∈ 1:4
        putmarkers!(r, direction_of_movement[i], direction_to_border[i])
    end
end

get_directions(side::HorizonSide) =
    if side == Nord
    # - обход будет по часовой стрелке
        return (Nord,Ost,Sud, West), (Ost,Sud,West,Nord)
    else # - обход будет против часовой стрелки
        return (Sud,Ost,Nord,West), (Ost,Nord,West,Sud)
    end


text:
using StippleCharts

vector = [1, 2, 3, 4];

PlotData(vector)


text:
# center is (0,0,0)
# not working for general case
function gen_dr( r, LatVecs )
    #
    v1 = LatVecs[:,1]
    v2 = LatVecs[:,2]
    v3 = LatVecs[:,3]
    #
    Npoints = size(r)[2]
    dr = Array{Float64}(Npoints)
    #
    Lx = maximum( abs(LatVecs[:,1]) )
    Ly = maximum( abs(LatVecs[:,2]) )
    Lz = maximum( abs(LatVecs[:,3]) )
    @printf("VecLen: (%18.10f,%18.10f,%18.10f)\n", Lx,Ly,Lz)
    # r are guaranted to contained within the unit cell
    for ip=1:Npoints

        x = r[1,ip] < 0.0 ? r[1,ip] + Lx : r[1,ip]
        y = r[2,ip] < 0.0 ? r[2,ip] + Ly : r[2,ip]
        z = r[3,ip] < 0.0 ? r[3,ip] + Lz : r[3,ip]

        dx = x > 0.5*Lx ? Lx-x : x
        dy = y > 0.5*Ly ? Ly-y : y
        dz = z > 0.5*Lz ? Lz-z : z

        dr[ip] = sqrt(dx^2 + dy^2 + dz^2)
        @printf("\nWrapped: (%10.5f,%10.5f,%10.5f)\n", x, y, z);
        @printf("%8d (%10.5f,%10.5f,%10.5f) (%10.5f,%10.5f,%10.5f): %18.10f\n",
                    ip, r[1,ip], r[2,ip], r[3,ip],
                    dx, dy, dz, dr[ip] )
    end
    return dr
end


text:
mutable struct Frame{E}
    entities::Vector{E} # NOTE: I tried StaticArrays; was not faster
    n::Int
end
function Frame(arr::AbstractVector{E}; capacity::Int=length(arr)) where {E}
    capacity ≥ length(arr) || error("capacity cannot be less than entitiy count! (N ≥ length(arr))")
    entities = Array{E}(undef, capacity)
    copyto!(entities, arr)
    return Frame{E}(entities, length(arr))
end
function Frame(::Type{E}, capacity::Int=100) where {E}
    entities = Array{E}(undef, capacity)
    return Frame{E}(entities, 0)
end

Base.show(io::IO, frame::Frame{E}) where {E}= @printf(io, "Frame{%s}(%d entities)", string(E), length(frame))

capacity(frame::Frame) = length(frame.entities)
Base.length(frame::Frame) = frame.n
Base.getindex(frame::Frame, i::Int) = frame.entities[i]
Base.eltype(frame::Frame{E}) where {E} = E

Base.lastindex(frame::Frame) = frame.n
function Base.setindex!(frame::Frame{E}, entity::E, i::Int) where {E}
    frame.entities[i] = entity
    return frame
end
function Base.empty!(frame::Frame)
    frame.n = 0
    return frame
end
function Base.deleteat!(frame::Frame, entity_index::Int)
    for i in entity_index : frame.n - 1
        frame.entities[i] = frame.entities[i+1]
    end
    frame.n -= 1
    frame
end

function Base.iterate(frame::Frame{E}, i::Int=1) where {E}
    if i > length(frame)
        return nothing
    end
    return (frame.entities[i], i+1)
end

function Base.copyto!(dest::Frame{E}, src::Frame{E}) where {E}
    for i in 1 : src.n
        dest.entities[i] = src.entities[i]
    end
    dest.n = src.n
    return dest
end
Base.copy(frame::Frame{E}) where {E} = copyto!(Frame(E, capacity(frame)), frame)

function Base.push!(frame::Frame{E}, entity::E) where {E}
    frame.n += 1
    frame.entities[frame.n] = entity
    return frame
end

####

const EntityFrame{S,D,I} = Frame{Entity{S,D,I}}
EntityFrame(::Type{S},::Type{D},::Type{I}) where {S,D,I} = Frame(Entity{S,D,I})
EntityFrame(::Type{S},::Type{D},::Type{I},N::Int) where {S,D,I} = Frame(Entity{S,D,I}, N)

Base.in(id::I, frame::EntityFrame{S,D,I}) where {S,D,I} = findfirst(id, frame) != nothing
function Base.findfirst(id::I, frame::EntityFrame{S,D,I}) where {S,D,I}
    for entity_index in 1 : frame.n
        entity = frame.entities[entity_index]
        if entity.id == id
            return entity_index
        end
    end
    return nothing
end
function id2index(frame::EntityFrame{S,D,I}, id::I) where {S,D,I}
    entity_index = findfirst(id, frame)
    if entity_index === nothing
        throw(BoundsError(frame, [id]))
    end
    return entity_index
end
get_by_id(frame::EntityFrame{S,D,I}, id::I) where {S,D,I} = frame[id2index(frame, id)]
function get_first_available_id(frame::EntityFrame{S,D,I}) where {S,D,I}
    ids = Set{I}(entity.id for entity in frame)
    id_one = one(I)
    id = id_one
    while id ∈ ids
        id += id_one
    end
    return id
end
function Base.push!(frame::EntityFrame{S,D,I}, s::S) where {S,D,I}
    id = get_first_available_id(frame)
    entity = Entity{S,D,I}(s, D(), id)
    push!(frame, entity)
end

Base.delete!(frame::EntityFrame{S,D,I}, entity::Entity{S,D,I}) where {S,D,I} = deleteat!(frame, findfirst(entity.id, frame))
function Base.delete!(frame::EntityFrame{S,D,I}, id::I) where {S,D,I}
    entity_index = findfirst(id, frame)
    if entity_index != nothing
        deleteat!(frame, entity_index)
    end
    return frame
end

###

function Base.write(io::IO, mime::MIME"text/plain", frames::Vector{EntityFrame{S,D,I}}) where {S,D,I}
    println(io, length(frames))
    for frame in frames
        println(io, length(frame))
        for entity in frame
            write(io, mime, entity.state)
            print(io, "\n")
            write(io, mime, entity.def)
            print(io, "\n")
            write(io, mime, entity.id)
            print(io, "\n")
        end
    end
end
function Base.read(io::IO, mime::MIME"text/plain", ::Type{Vector{EntityFrame{S,D,I}}}) where {S,D,I}

    n = parse(Int, readline(io))
    frames = Array{EntityFrame{S,D,I}}(undef, n)

    for i in 1 : n
        m = parse(Int, readline(io))
        frame = Frame(Entity{S,D,I}, m)
        for j in 1 : m
            state = read(io, mime, S)
            def = read(io, mime, D)
            id = read(io, mime, I)
            push!(frame, Entity(state,def,id))
        end
        frames[i] = frame
    end

    return frames
end

max_stars_repo_name: tawheeler/Records

max_issues_repo_name: tawheeler/Records

max_forks_repo_name: tawheeler/Records


text:
addprocs( max(0, Sys.CPU_CORES-nprocs()) )

@everywhere module SetAndForget
  import Julz

  is_worker = ( dirname(@__FILE__) == "" )

  cd("$(string($(dirname(@__FILE__))))/..") do
    if is_worker
      include("config/bootload.jl")
    else
      include("../config/bootload.jl")
    end
  end

  #Define elementary operations on remote data
  import Base: length, +, *
  length(r1::Future)=length(Base.fetch(r1))
  +(r1::Future,r2::Future)=@spawnat r2.where Base.fetch(r1)+Base.fetch(r2)
  *(r1::Future,r2::Future)=@spawnat r2.where Base.fetch(r1)*Base.fetch(r2)

  function main()
    load_input("input.jl", true)

    gc_enable(false)

    r , s = @setup_procs
    n_run = Int( max(1, 5 - floor(n/1000)) )

    t_ser = @time_function n_run        prefix_serial!(  s, *)

    t_par = @time_function n_run @sync( prefix_parallel!(r, *) )

    open("output.jl", "w") do cur_file
      write(cur_file, "output = $(string(t_ser / t_par))" )
    end

    gc_enable(true)
  end
end


text:
"""
    ChainedVector(arrays::Vector{<:AbstractVector})

Create a `ChainedVector` of a `Vector` of homogenously-typed `AbstractVector`.
The "chain" of input vectors will be treated as a single, long vector.
A full set of typical mutable operations are supported (e.g. `push!`, `append!`, etc.).

As implementation details, mutable operations on single elements (e.g. `setindex!`, `push!`)
operate in-place or mutate an existing chained array, while `append!`/`prepend!` are optimized
to "chain" the incoming array to the existing chained arrays.
"""
struct ChainedVector{T, A <: AbstractVector{T}} <: AbstractVector{T}
    arrays::Vector{A}
    inds::Vector{Int}
end

function ChainedVector(arrays::Vector{A}) where {A <: AbstractVector{T}} where {T}
    n = length(arrays)
    inds = Vector{Int}(undef, n)
    x = 0
    @inbounds for i = 1:n
        # note that arrays[i] can have zero length
        x += length(arrays[i])
        inds[i] = x
    end
    return ChainedVector{T, A}(arrays, inds)
end

Base.IndexStyle(::Type{<:ChainedVector}) = Base.IndexLinear()
Base.size(x::ChainedVector) = (length(x.inds) == 0 ? 0 : x.inds[end],)

@inline function index(A::ChainedVector, i::Integer)
    chunk = searchsortedfirst(A.inds, i)
    return chunk, i - (chunk == 1 ? 0 : @inbounds A.inds[chunk - 1])
end

Base.@propagate_inbounds function Base.getindex(A::ChainedVector, i::Integer)
    @boundscheck checkbounds(A, i)
    chunk, ix = index(A, i)
    @inbounds x = A.arrays[chunk][ix]
    return x
end

Base.@propagate_inbounds function Base.isassigned(A::ChainedVector, i::Integer)
    @boundscheck checkbounds(A, i)
    chunk, ix = index(A, i)
    return @inbounds isassigned(A.arrays[chunk], ix)
end

# Base.@propagate_inbounds function Base.getindex(A::ChainedVector, chunk::Integer, i::Integer)
#     @boundscheck checkbounds(A.arrays, chunk)
#     @boundscheck checkbounds(A.arrays[chunk], i)
#     @inbounds x = A.arrays[chunk][i]
#     return x
# end

Base.@propagate_inbounds function Base.setindex!(A::ChainedVector, v, i::Integer)
    @boundscheck checkbounds(A, i)
    chunk, ix = index(A, i)
    @inbounds A.arrays[chunk][ix] = v
    return v
end

# efficient iteration
@inline function Base.iterate(A::ChainedVector)
    length(A) == 0 && return nothing
    i = 1
    chunk = 1
    chunk_i = 1
    chunk_len = A.inds[1]
    while i > chunk_len
        chunk += 1
        @inbounds chunk_len = A.inds[chunk]
    end
    x = A.arrays[chunk][1]
    # find next valid index
    i += 1
    if i > chunk_len
        while true
            chunk += 1
            chunk > length(A.inds) && break
            @inbounds chunk_len = A.inds[chunk]
            i <= chunk_len && break
        end
    else
        chunk_i += 1
    end
    return x, (i, chunk, chunk_i, chunk_len, length(A))
end

@inline function Base.iterate(A::ChainedVector, (i, chunk, chunk_i, chunk_len, len))
    i > len && return nothing
    @inbounds x = A.arrays[chunk][chunk_i]
    i += 1
    if i > chunk_len
        chunk_i = 1
        while true
            chunk += 1
            chunk > length(A.inds) && break
            @inbounds chunk_len = A.inds[chunk]
            i <= chunk_len && break
        end
    else
        chunk_i += 1
    end
    return x, (i, chunk, chunk_i, chunk_len, len)
end

# other AbstractArray functions
function Base.empty!(A::ChainedVector)
    empty!(A.arrays)
    empty!(A.inds)
    return A
end

function Base.copy(A::ChainedVector{T, AT}) where {T, AT}
    B = similar(AT, length(A))
    off = 1
    for arr in A.arrays
        n = length(arr)
        copyto!(B, off, arr, 1, n)
        off += n
    end
    return B
end

function Base.resize!(A::ChainedVector{T, AT}, len) where {T, AT}
    len >= 0 || throw(ArgumentError("`len` must be >= 0 when resizing ChainedVector"))
    len′ = length(A)
    if len′ < len
        # growing
        push!(A.arrays, similar(AT, len - len′))
        push!(A.inds, len)
    else
        # shrinking
        chunk = searchsortedfirst(A.inds, len)
        # get rid of any excess chunks
        resize!(A.arrays, chunk)
        resize!(A.inds, chunk)
        # resize individual chunk
        resize!(A.arrays[chunk], A.inds[chunk] - len)
        A.inds[chunk] -= A.inds[chunk] - len
    end
    return A
end

function Base.push!(A::ChainedVector{T, AT}, val) where {T, AT}
    if length(A.arrays) == 0
        push!(A.arrays, similar(AT, 0))
        push!(A.inds, 0)
    end
    @inbounds push!(A.arrays[end], val)
    @inbounds A.inds[end] += 1
    return A
end

function Base.pushfirst!(A::ChainedVector{T, AT}, val) where {T, AT}
    if length(A.arrays) == 0
        push!(A.arrays, similar(AT, 0))
        push!(A.inds, 0)
    end
    @inbounds pushfirst!(A.arrays[1], val)
    for i = 1:length(A.inds)
        @inbounds A.inds[i] += 1
    end
    return A
end

Base.@propagate_inbounds function Base.deleteat!(A::ChainedVector, i::Integer)
    @boundscheck checkbounds(A, i)
    chunk, ix = index(A, i)
    deleteat!(A.arrays[chunk], ix)
    for j = chunk:length(A.inds)
        @inbounds A.inds[j] -= 1
    end
    # check if we should remove an empty chunk
    if length(A.arrays[chunk]) == 0
        deleteat!(A.arrays, chunk)
        deleteat!(A.inds, chunk)
    end
    return A
end

Base.@propagate_inbounds function Base.deleteat!(A::ChainedVector, inds)
    y = iterate(inds)
    y === nothing && return A
    i, s = y
    chunk = 1
    N = length(A.inds)
    prevind = 0
    ind = chunk > N ? 0 : A.inds[chunk]
    todelete = Int[]
    while y !== nothing
        # find chunk where deleting starts
        while ind < i
            chunk += 1
            chunk > N && throw(BoundsError(A, i))
            prevind = ind
            @inbounds ind = A.inds[chunk]
        end
        # gather all indices for this chunk
        while i <= ind
            push!(todelete, i - prevind)
            y = iterate(inds, s)
            y === nothing && break
            i, s = y
        end
        # delete indices from this chunk
        @inbounds deleteat!(A.arrays[chunk], todelete)
        empty!(todelete)
    end
    # reset inds
    x = 0
    @inbounds for j = 1:N
        # note that A.arrays[j] can have zero length
        len = length(A.arrays[j])
        if len == 0
            push!(todelete, j)
        else
            x += len
            A.inds[j] = x
        end
    end
    for j in Iterators.reverse(todelete)
        deleteat!(A.arrays, j)
        deleteat!(A.inds, j)
    end
    return A
end

Base.@propagate_inbounds function Base.deleteat!(A::ChainedVector, inds::AbstractVector{Bool})
    length(inds) == length(A) || throw(BoundsError(A, inds))
    for i = length(A):-1:1
        if inds[i]
            deleteat!(A, i)
        end
    end
    return A
end

function Base.pop!(A::ChainedVector)
    if isempty(A)
        throw(ArgumentError("array must be non-empty"))
    end
    item = A[end]
    deleteat!(A, length(A))
    return item
end

function Base.popfirst!(A::ChainedVector)
    if isempty(A)
        throw(ArgumentError("array must be non-empty"))
    end
    item = A[1]
    deleteat!(A, 1)
    return item
end

Base.@propagate_inbounds function Base.insert!(A::ChainedVector{T, AT}, i::Integer, item) where {T, AT <: AbstractVector{T}}
    if i == 1 && length(A.arrays) == 0
        push!(A.arrays, similar(AT, 0))
        push!(A.inds, 0)
        chunk, ix = 1, 1
    else
        @boundscheck checkbounds(A, i)
        chunk, ix = index(A, i)
    end
    insert!(A.arrays[chunk], ix, item)
    for j = chunk:length(A.inds)
        @inbounds A.inds[j] += 1
    end
    return A
end

function Base.vcat(A::ChainedVector{T, AT}, arrays::ChainedVector{T, AT}...) where {T, AT <: AbstractVector{T}}
    newarrays = vcat(A.arrays, map(x->x.arrays, arrays)...)
    n = length(A.inds)
    inds = Vector{Int}(undef, n + sum(x->length(x.inds), arrays))
    copyto!(inds, 1, A.inds, 1, n)
    m = n + 1
    for x in arrays
        for y in x.arrays
            @inbounds inds[m] = ((m - 1) == 0 ? 0 : inds[m - 1]) + length(y)
            m += 1
        end
    end
    return ChainedVector{T, AT}(newarrays, inds)
end

function Base.append!(A::ChainedVector{T, AT}, B::AT) where {T, AT <: AbstractVector{T}}
    push!(A.arrays, B)
    push!(A.inds, A.inds[end] + length(B))
    return A
end

function Base.append!(A::ChainedVector{T, AT}, B::ChainedVector{T, AT}) where {T, AT <: AbstractVector{T}}
    append!(A.arrays, B.arrays)
    n = length(A.inds)
    m = length(B.inds)
    resize!(A.inds, n + m)
    for i = 1:m
        @inbounds A.inds[n + i] = ((n + i - 1) == 0 ? 0 : A.inds[n + i - 1]) + length(B.arrays[i])
    end
    return A
end

function Base.append!(A::ChainedVector{T}, B) where {T}
    for x in B
        push!(A, x)
    end
    return A
end

function Base.prepend!(A::ChainedVector{T, AT}, B::AT) where {T, AT <: AbstractVector{T}}
    pushfirst!(A.arrays, B)
    n = length(B)
    pushfirst!(A.inds, n)
    for i = 2:length(A.inds)
        @inbounds A.inds[i] += n
    end
    return A
end

function Base.prepend!(A::ChainedVector{T, AT}, B::ChainedVector{T, AT}) where {T, AT <: AbstractVector{T}}
    prepend!(A.arrays, B.arrays)
    n = length(A.inds)
    m = length(B.inds)
    M = length(B)
    prepend!(A.inds, B.inds)
    for i = 1:n
        @inbounds A.inds[m + i] += M
    end
    return A
end

function Base.prepend!(A::ChainedVector{T}, B) where {T}
    for (i, x) in enumerate(B)
        insert!(A, i, x)
    end
    return A
end


text:
# heatmap sampler (experimental)

@info "IncrementalInference.jl is loading tools related to Interpolations.jl."

using .Interpolations

# only export on Requires.jl
export HeatmapGridDensity, PackedHeatmapGridDensity
export LevelSetGridNormal, PackedLevelSetGridNormal

export sampleHeatmap

##

getManifold(hgd::HeatmapGridDensity) = getManifold(hgd.densityFnc)
getManifold(lsg::LevelSetGridNormal) = getManifold(lsg.heatmap)

"""
    $SIGNATURES

Get the grid positions at the specified height (within the provided spreads)

DevNotes
- TODO Should this be consolidated with AliasingScalarSampler? See IIF #1341
"""
function sampleHeatmap( roi::AbstractMatrix{<:Real},
                        x_grid::AbstractVector{<:Real},
                        y_grid::AbstractVector{<:Real},
                        thres::Real = 0  )
  #

  # mask the region of interest above the sampling threshold value
  mask = thres .<= roi

  idx2d = findall(mask)  # 2D indices
  pos = (v->[x_grid[v[1]],y_grid[v[2]]]).(idx2d)
  weights = (v->roi[v[1],v[2]]).(idx2d)
  weights ./= sum(weights)

  pos, weights
end

# TODO make n-dimensional, and later on-manifold
# TODO better standardize for heatmaps on manifolds w MKD
function fitKDE(support,
                weights,
                x_grid::AbstractVector{<:Real},
                y_grid::AbstractVector{<:Real};
                bw_factor::Real=0.7  )
  #
  # 1. set the bandwidth
  x_spacing = Statistics.mean(diff(x_grid))
  y_spacing = Statistics.mean(diff(y_grid))
  kernel_ = bw_factor*0.5*(x_spacing + y_spacing) # 70% of the average spacing
  kernel_bw = [kernel_; kernel_]                  # same bw in x and y
  # fit KDE
  kde!(support, kernel_bw, weights)
end

# Helper function to construct HGD
function HeatmapGridDensity(data::AbstractMatrix{<:Real},
                            domain::Tuple{<:AbstractVector{<:Real},<:AbstractVector{<:Real}},
                            hist_callback::Union{<:Function, Nothing}=nothing,
                            bw_factor::Real=0.7;  # kde spread between domain points
                            N::Int=10000  )
  #
  pos, weights_ = sampleHeatmap(data, domain..., 0)
  # recast to the appropriate shape
  @cast support_[i,j] := pos[j][i]

  # constuct a pre-density from which to draw intermediate samples
  # TODO remove extraneous collect()
  density_ = fitKDE(collect(support_), weights_, domain...; bw_factor=bw_factor)
  pts_preIS, = sample(density_, N)

  @cast vec_preIS[j][i] := pts_preIS[i,j]

  # weight the intermediate samples according to interpolation of raw data
  hm = Interpolations.LinearInterpolation( domain, data ) # interpolated heatmap
  d_scalar = Vector{Float64}( undef, length(vec_preIS) )

  # interpolate d_scalar for intermediate test points
  for (i,u) in enumerate(vec_preIS)
    if maximum(domain[1]) < abs(u[1]) || maximum(domain[2]) < abs(u[2])
      d_scalar[i] = 0.0
      continue
    end
    d_scalar[i] = hm(u...)
  end

  #
  weights = exp.(-d_scalar) # unscaled Gaussian
  weights ./= sum(weights)  # normalized

  # final samplable density object
  # TODO better standardize for heatmaps on manifolds
  bw = getBW(density_)[:,1]
  @cast pts[i,j] := vec_preIS[j][i]
  bel = kde!(collect(pts), bw, weights)
  density = ManifoldKernelDensity(TranslationGroup(Ndim(bel)), bel)

  # return `<:SamplableBelief` object
  HeatmapGridDensity(data, domain, hist_callback, bw_factor, density)
end

# legacy construct helper
function LevelSetGridNormal(data::AbstractMatrix{<:Real},
                            domain::Tuple{<:AbstractVector{<:Real},<:AbstractVector{<:Real}},
                            level::Real,
                            sigma::Real;
                            sigma_scale::Real=3,
                            hist_callback::Union{<:Function, Nothing}=nothing,
                            bw_factor::Real=0.7,  # kde spread between domain points
                            N::Int=10000  )
  #

  # select the support from raw data
  roi = data.-level
  # make Gaussian
  roi .^= 2
  roi .*= 0.5/(sigma^2)
  roi .-= sigma_scale^2
  roi .*= -1
  # truncate sigma_scale*sigma below zero
  #   h = heatmap;  z = measurement
  #   l = 1/2 (h-z/σ)^2
  #   masked_roi = 0 .< κ^2 - l

  hgd = HeatmapGridDensity(data, domain, hist_callback, bw_factor; N=N)

  LevelSetGridNormal(level, sigma, float(sigma_scale), hgd)
end

#


text:
function mcsweep(model::IsingModel)
  ndata::Int64 = length(model.data)
  randinds = rand(model.rng,1:ndata,ndata)
  for ind in randinds # loop entries in the data
    ediff::Float64 = localediff(ind,model)
    if ediff < 0
      flipbit!(ind,model)
    else if rand(model.rng) < exp(-ediff/model.T)
      flipbit!(ind,model)
    end
  end
end

function mcsweeps(nsweeps::Int64,model::IsingModel;cbevery::Int64=1,cb::Function=defaultcb,cbargs=nothing)
  ndata::Int64 = length(model.data)
  for i in 1:nsweeps # loop over sweeps
    if i%cbevery == 0 && !(cb == nothing)
      cb(model,cbargs)
    end
    mcsweep(model)
  end
end


text:
module TableViewExample

using Blink, CSV, DataFrames, WebIO, CSSUtil, Interact
using TableView.TableViewAdvanced

println("Starting [$(splitext(basename(@__FILE__))[1])]")

#=======================================
load the data for this example
=======================================#
current_file_folder = dirname(@__FILE__)
csv = CSV.File(joinpath(current_file_folder, "Data", "Data01.csv"))
df = DataFrame(csv)
# convert all the column to type Any to make it easier to process any type of imput and replace all missing values with spaces
for col in names(df)
    df[!,Symbol(col)] = Vector{Any}(df[!,Symbol(col)])
    replace!(df[!, col], missing => "");
end

#=======================================
grid event handlers
=======================================#
function OnCellValueChanged(; msg, table::Table, kwargs...)
    win = table_user_data(table)
    data = table_data(table)    # assuming a data frame

    # update the back-end data structure
    row = msg["row"]
    col = msg["col"]
    @info "Cell Change: $row, $col -> [$(msg["newValue"])]"     # debug message
    data[row, col] = msg["newValue"]
end

#=======================================
define table options
=======================================#
table_options = Dict{Symbol, Any}(
    :columnDefs => Dict(
        :ObjectID => Dict(
            :headerName => "Object ID",
            :headerTooltip => "A Unique name to identify this row",
        ),
        :Comment => Dict(
            :headerName => "Comment",
            :cellEditor => "agLargeTextCellEditor",
        ),
        :ObjectType => Dict(
            :headerName => "Object Type",
            :cellEditor => "agSelectCellEditor",
            :cellEditorParams => Dict(
                values => ["SOURCE", "Spherical Lens", "STOP", "DETECTOR"],
            ),
        ),
        :Position => Dict(
            :headerName => "Position(x,y,z,θx,θy,θz)",
        ),
        :RelativeTo => Dict(
            :headerName => "Relative To",
        ),
        :Material => Dict(
            :headerName => "Material",
        ),
        :Parameters => Dict(
            :headerName => "Parameters (radium, thinkness, conic, etc)",
        ),

    ),
    :defaultColDef => Dict(
        :editable => true,
        :sortable => true,
        :resizable => true,
    ),
    :events => Dict(
        :CellValueChanged => OnCellValueChanged,
    )
)

table = TableViewAdvanced.Table(df, table_options)

#=======================================
Create Blink Window
=======================================#
window_defaults = Blink.@d(
    :title => "Grid Example 01",
    :width=>1200,
    :height=>800,
)
win = Window(window_defaults)

# set the window as the table's "user data"  so we can use it inside grid events to run JS code
set_table_user_data(table, win)

#=======================================
create some of the ui widgets
=======================================#
toggle_dev_tools_button = button("Toggle Dev Tools")
on(toggle_dev_tools_button) do val
    Blink.AtomShell.tools(win)
end

show_in_repl_button = button("Show Data in REPL")
on(show_in_repl_button) do val
    @info df
end

checkbox_test = checkbox(label="""Hide the "Position" Column""")
on(checkbox_test) do val
    global table
    win = table_user_data(table)

    @js_ win begin
        x = [
            Dict(:field => "Position", :hide => $(val)),
        ]
        TVJS_UpdateColumns($(table_scope_id(table)), x)
    end
end

ui_widgets = hbox(
    toggle_dev_tools_button,
    Node(:span, attributes=Dict(:style => "width: 10px;")),
    show_in_repl_button,
    checkbox_test,
)

#=======================================
create some of the ui widgets
=======================================#

ui = Node(:dom,
    JS_loader_node(joinpath(dirname(@__FILE__), "examples.css")),
    Node(:h1, "Grid Example 02"),
    Node(:h2, "Some of the features we show in this example:."),
    Node(:ul,
        Node(:li, "Large text editor for editing the Comment field"),
        Node(:li, "ComboBox style editor for the ObjectType field"),
        Node(:li, "Hide/Show column toggle that show how to update columns properties in run-time."),
        Node(:li, "CellValueChanged - a grid event that shows how to update the back-end data structure (a DataFRame in this case)"),
    ),
    Node(:hr),
    ui_widgets,
    Node(:hr),
    Node(:div, table_scope(table)),
)

body!(win, ui, async=false)

println("Ending [$(splitext(basename(@__FILE__))[1])]")
end # module TableViewExample