text:
function writeStamp(dirname::String, t::Float64, surf::ThreeDSurfSimple, curfield::ThreeDFieldSimple)
    dirvec = readdir()
    if dirname in dirvec
        rm(dirname, recursive=true)
    end
    mkdir(dirname)
    cd(dirname)

    f = open("time", "w")
    Serialization.Serialization.serialize(f, ["#time \t"])
    DelimitedFiles.writedlm(f, [t]')
    close(f)

    for i = 1:surf.nspan
        f = open("kinem-$i", "w")
        Serialization.Serialization.serialize(f, ["alpha (deg) \t", "h/c \t", "u/uref \t"])
        DelimitedFiles.writedlm(f, [surf.s2d[i].kinem.alpha, surf.s2d[i].kinem.h, surf.s2d[i].kinem.u]')
        close(f)

        f = open("FourierCoeffs-$i", "w")
        Serialization.serialize(f, ["#Fourier coeffs (0-n) \t", "d/dt (Fourier coeffs) \n"])
        matfour = zeros(surf.s2d[i].naterm+1, 2)
        matfour[:,1] = [surf.s2d[i].a0[1]; surf.s2d[i].aterm[:]]
        matfour[1:4,2] = [surf.s2d[i].a0dot[1];surf.s2d[i].adot[:]]
        DelimitedFiles.writedlm(f, matfour)
        close(f)

        f = open("tev-$i", "w")
        Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
        tevmat = zeros(length(curfield.f2d[i].tev), 3)
        for iv = 1:length(curfield.f2d[i].tev)
            tevmat[iv,:] = [curfield.f2d[i].tev[iv].s curfield.f2d[i].tev[iv].x curfield.f2d[i].tev[iv].z]
        end
        DelimitedFiles.writedlm(f, tevmat)
        close(f)

        f = open("lev-$i", "w")
        Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
        levmat = zeros(length(curfield.f2d[i].lev), 3)
        for iv = 1:length(curfield.f2d[i].lev)
            levmat[iv,:] = [curfield.f2d[i].lev[iv].s curfield.f2d[i].lev[iv].x curfield.f2d[i].lev[iv].z]
        end
        DelimitedFiles.writedlm(f, levmat)
        close(f)

        f = open("boundv-$i", "w")
        Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
        bvmat = zeros(length(surf.s2d[i].bv), 3)
        for iv = 1:length(surf.s2d[i].bv)
            bvmat[iv,:] = [surf.s2d[i].bv[iv].s surf.s2d[i].bv[iv].x surf.s2d[i].bv[iv].z]
        end
        DelimitedFiles.writedlm(f, bvmat)
        close(f)
    end

    #Write spanwise properties of interest in a file
    f = open("spanwise-var", "w")
    Serialization.serialize(f, ["#y_le \t", "shed TEV str \t", "shed LEV str \t", "Gamma \t", "A0 2D \t", "A03D \t", "A0 tot \n"])
    mat = zeros(0,7)
    mat = [mat; [-surf.AR/2, 0. , 0. , 0. , 0. , 0. , 0. ]']
    for i = 1:surf.nspan
        yle = surf.yle[i]
        tev_s = curfield.f2d[i].tev[end].s
        if surf.s2d[i]. levflag == 1
            lev_s = curfield.f2d[i].lev[end].s
        else
            lev_s = 0.
        end
        gam = surf.s2d[i].uref*surf.s2d[i].c*pi*(surf.s2d[i].a0[1] + 0.5*surf.s2d[i].aterm[1])
        a02d = surf.s2d[i].a0[1] - surf.a03d[i]

        mat = [mat; [yle, tev_s, lev_s, gam, a02d, surf.a03d[i], surf.s2d[i].a0[1]]']
    end
    DelimitedFiles.writedlm(f, mat)
    close(f)

    cd("..")
end

function calc_forces(surf::ThreeDSurfSimple, field :: ThreeDFieldSimple, dt :: Float64)

    cl = zeros(surf.nspan)
    cd = zeros(surf.nspan)
    cm = zeros(surf.nspan)

    lev_s = zeros(surf.nspan)
    for i = 1:surf.nspan
        if surf.s2d[i].levflag[1] == 1
            lev_s[i] = field.f2d[i].lev[end].s
        end
    end

    for i = 1:surf.nspan
        # First term in eqn (2.30) Ramesh et al. in coefficient form
        wi = 0
        for n = 1:surf.nspan
            nn = 2*n -2
            wi += real(nn)*surf.bcoeff[n]*sin(nn*surf.psi[i])/sin(surf.psi[i])
        end

        cnc = 2*pi*(surf.s2d[i].kinem.u*cos(surf.s2d[i].kinem.alpha)/surf.s2d[i].uref + surf.s2d[i].kinem.hdot*sin(surf.s2d[i].kinem.alpha)/surf.s2d[i].uref + wi*sin(surf.s2d[i].kinem.alpha))*(surf.s2d[i].a0[1] + surf.s2d[i].aterm[1]/2.)

        # Second term in eqn (2.30) Ramesh et al. in coefficient form
        cnnc = 2*pi*(3*surf.s2d[i].c*surf.s2d[i].a0dot[1]/(4*surf.s2d[i].uref) + surf.s2d[i].c*surf.s2d[i].adot[1]/(4*surf.s2d[i].uref) + surf.s2d[i].c*surf.s2d[i].adot[2]/(8*surf.s2d[i].uref)) + (2*lev_s[i]/(dt*surf.s2d[i].uref^2))

        # Suction force given in eqn (2.31) Ramesh et al.
        cs = 2*pi*surf.s2d[i].a0[1]*surf.s2d[i].a0[1]

        #The components of normal force and moment from induced velocities are calulcated in dimensional units and nondimensionalized later
        nonl=0
        nonl_m=0
        for ib = 1:surf.s2d[i].ndiv-1
            nonl = nonl + (surf.s2d[i].uind[ib]*cos(surf.s2d[i].kinem.alpha) - surf.s2d[i].wind[ib]*sin(surf.s2d[i].kinem.alpha))*surf.s2d[i].bv[ib].s
            nonl_m = nonl_m + (surf.s2d[i].uind[ib]*cos(surf.s2d[i].kinem.alpha) - surf.s2d[i].wind[ib]*sin(surf.s2d[i].kinem.alpha))*surf.s2d[i].x[ib]*surf.s2d[i].bv[ib].s
        end
        nonl = nonl*2. /(surf.s2d[i].uref*surf.s2d[i].uref*surf.s2d[i].c)
        nonl_m = nonl_m*2. /(surf.s2d[i].uref*surf.s2d[i].uref*surf.s2d[i].c*surf.s2d[i].c)

        # Normal force coefficient
        cn = cnc + cnnc + nonl

        # Lift and drag coefficients
        cl[i] = cn*cos(surf.s2d[i].kinem.alpha) + cs*sin(surf.s2d[i].kinem.alpha)
        cd[i] = cn*sin(surf.s2d[i].kinem.alpha) - cs*cos(surf.s2d[i].kinem.alpha)

        #Pitching moment is clockwise or nose up positive
        cm[i] = cn*surf.s2d[i].pvt - 2*pi*((surf.s2d[i].kinem.u*cos(surf.s2d[i].kinem.alpha)/surf.s2d[i].uref + surf.s2d[i].kinem.hdot*sin(surf.s2d[i].kinem.alpha)/surf.s2d[i].uref)*(surf.s2d[i].a0[1]/4. + surf.s2d[i].aterm[1]/4. - surf.s2d[i].aterm[2]/8.) + (surf.s2d[i].c/surf.s2d[i].uref)*(7. *surf.s2d[i].a0dot[1]/16. + 3. *surf.s2d[i].adot[1]/16. + surf.s2d[i].adot[2]/16. - surf.s2d[i].adot[3]/64. )) - nonl_m
    end

    cl3d = 0
    cd3d = 0
    cm3d = 0

    for is = 1:surf.nspan
        if is == 1
            yl = -0.5*surf.AR*surf.cref
            yu = 0.5*(surf.yle[is] + surf.yle[is+1])
        elseif is == surf.nspan
            yl = 0.5*(surf.yle[is] + surf.yle[is-1])
            yu = 0.
        else
            yl = 0.5*(surf.yle[is] + surf.yle[is-1])
            yu = 0.5*(surf.yle[is] + surf.yle[is+1])
        end
        cl3d += cl[is]*(yu-yl)
        cd3d += cd[is]*(yu-yl)
        cm3d += cm[is]*(yu-yl)
    end
    cl3d = cl3d/(0.5*surf.AR*surf.cref)
    cd3d = cd3d/(0.5*surf.AR*surf.cref)
    cm3d = cm3d/(0.5*surf.AR*surf.cref)

    return cl3d, cd3d, cm3d, cl, cd, cm
end


text:
# abstract type LogdensityComputation end

# struct LogpdfComputation end

# struct BasemeasureComputation{B} <: LogdensityComputation
#     base::B
# end

# struct LogdensityThunk{C,M,X,L}
#     computation::C
#     measure::M
#     x::X
#     ℓ::L
# end

# function reduce_step(thunk :: LogdensityThunk, callback=Returns(nothing))

# reduce_step(LogdensityThunk{M,B,X,L}, callback=Returns(nothing)) -> Union{LogdensityThunk, L}

# function logdensity(μ, b, x; cb=Returns(nothing))
#     f(thunk) = reduce_step(thunk, cb)
#     fix(f, LogdensityThunk(μ, b, x, zero(Float64)))
# end

ext: jl


text:
struct TreeReader{T}
    io::T

    function TreeReader(input::String)
        io = isfile(input) ? open(input) : IOBuffer(input)
        return new{typeof(io)}(io)
    end
end

function next_tree(reader::TreeReader)
    depth = 0
    input, buf = reader.io, IOBuffer()
    while !eof(input)
        c = Char(read(reader.io, 1)[1])
        !isvalid(c) && continue
        c == '(' && (depth += 1)
        c == ')' && (depth -= 1)
        write(buf, c)
        c == ')' && depth == 0 && break
    end
    str = String(take!(buf))
    try
        tree = read_tree(str)
        readuntil(input, '(')
        return tree
    catch err
        @show err str
        throw(err)
    end
end

function Base.iterate(reader::TreeReader, state...)
    if eof(reader.io)
        close(reader.io)
        return nothing
    else
        return (next_tree(reader), position(reader.io))
    end
end

struct Treebank
    files::Vector{String}
end

"""
    Treebank(corpus)

A Treebank is an iterator over a corpus of trees in bracketed format.
"""
function Treebank(corpus::String)
    if isfile(corpus)
        return Treebank([corpus])
    elseif isdir(corpus)
        f = x -> joinpath(corpus, x)
        return Treebank(filter(isfile, f.(readdir(corpus))))
    else
        error("don't know how to read treebank '$corpus'")
    end
end

function Base.iterate(treebank::Treebank)
    iter = TreebankIterator(treebank, 1, TreeReader(treebank.files[1]))
    tree, state = iterate(iter)
    return (tree, (iter, state))
end
function Base.iterate(treebank::Treebank, state)
    (iter, st) = state
    next = iterate(iter, st)
    if next == nothing
        return nothing
    else
        tree, next_state = next
        return (tree, (iter, next_state))
    end
end

Base.IteratorEltype(::Type{<:Treebank}) = Base.HasEltype()
Base.eltype(::Type{<:Treebank})         = ConstituencyTree
Base.IteratorSize(::Type{<:Treebank})   = Base.SizeUnknown()

mutable struct TreebankIterator
    tb::Treebank
    i::Int
    r::TreeReader
end
function Base.iterate(t::TreebankIterator)
    return iterate(t.r)
end
function Base.iterate(t::TreebankIterator, state)
    next = iterate(t.r, state)
    if next === nothing
        close(t.r.io)
        if t.i < length(t.tb.files)
            t.i += 1
            t.r = TreeReader(t.tb.files[t.i])
            return iterate(t.r)
        end
    else
        return next
    end
end


text:
# E. Dynamic Parallelism

import CUDAdrv: CuDim3, CuStream_t
const cudaStream_t = CUDAdrv.CuStream_t

if VERSION >= v"1.2.0-DEV.512"
    @inline cudaLaunchDevice(buf::Ptr{Cvoid}, stream::CuStream) =
        ccall("extern cudaLaunchDeviceV2", llvmcall, cudaError_t,
              (Ptr{Cvoid}, cudaStream_t),
              buf, stream)
else
    import Base.Sys: WORD_SIZE
    # declare i32 @cudaLaunchDeviceV2(i8*, %struct.CUstream_st*)
    @eval @inline cudaLaunchDevice(buf::Ptr{Cvoid}, stream::CuStream) =
            Base.llvmcall(
                $("declare i32 @cudaLaunchDeviceV2(i8*, i8*)",
                  "%buf = inttoptr i$WORD_SIZE %0 to i8*
                   %stream = inttoptr i$WORD_SIZE %1 to i8*
                   %rv = call i32 @cudaLaunchDeviceV2(i8* %buf, i8* %stream)
                   ret i32 %rv"), cudaError_t,
                Tuple{Ptr{Cvoid}, cudaStream_t},
                buf, Base.unsafe_convert(cudaStream_t, stream))
end

# device-side counterpart of CUDAdrv.launch
@inline function launch(f::Ptr{Cvoid}, blocks::CuDim, threads::CuDim,
                        shmem::Int, stream::CuStream,
                        args...)
    blocks = CuDim3(blocks)
    threads = CuDim3(threads)

    buf = parameter_buffer(f, blocks, threads, shmem, args...)
    cudaLaunchDevice(buf, stream)

    return
end

if VERSION >= v"1.2.0-DEV.512"
    @inline cudaGetParameterBuffer(f::Ptr{Cvoid}, blocks::CuDim3, threads::CuDim3, shmem::Integer) =
        ccall("extern cudaGetParameterBufferV2", llvmcall, Ptr{Cvoid},
                            (Ptr{Cvoid}, CuDim3, CuDim3, Cuint),
                            f, blocks, threads, shmem)
else
    @inline cudaGetParameterBuffer(f::Ptr{Cvoid}, blocks::CuDim3, threads::CuDim3, shmem::Integer) =
        cudaGetParameterBuffer(f,
                               blocks.x, blocks.y, blocks.z,
                               threads.x, threads.y, threads.z,
                               convert(Cuint, shmem))
    # declare i8* @cudaGetParameterBufferV2(i8*, %struct.dim3, %struct.dim3, i32)
    @eval @inline cudaGetParameterBuffer(f::Ptr{Cvoid},
                                         blocks_x::Cuint, blocks_y::Cuint, blocks_z::Cuint,
                                         threads_x::Cuint, threads_y::Cuint, threads_z::Cuint,
                                         shmem::Cuint) =
            Base.llvmcall(
                $("declare i8* @cudaGetParameterBufferV2(i8*, {i32,i32,i32}, {i32,i32,i32}, i32)",
                  "%f = inttoptr i$WORD_SIZE %0 to i8*
                   %blocks.x = insertvalue { i32, i32, i32 } undef, i32 %1, 0
                   %blocks.y = insertvalue { i32, i32, i32 } %blocks.x, i32 %2, 1
                   %blocks.z = insertvalue { i32, i32, i32 } %blocks.y, i32 %3, 2
                   %threads.x = insertvalue { i32, i32, i32 } undef, i32 %4, 0
                   %threads.y = insertvalue { i32, i32, i32 } %threads.x, i32 %5, 1
                   %threads.z = insertvalue { i32, i32, i32 } %threads.y, i32 %6, 2
                   %rv = call i8* @cudaGetParameterBufferV2(i8* %f, {i32,i32,i32} %blocks.z, {i32,i32,i32} %threads.z, i32 %7)
                   %buf = ptrtoint i8* %rv to i$WORD_SIZE
                   ret i$WORD_SIZE %buf"), Ptr{Cvoid},
                Tuple{Ptr{Cvoid}, Cuint, Cuint, Cuint, Cuint, Cuint, Cuint, Cuint},
                f, blocks_x, blocks_y, blocks_z, threads_x, threads_y, threads_z, shmem)
end

@generated function parameter_buffer(f::Ptr{Cvoid}, blocks::CuDim3, threads::CuDim3,
                                     shmem::Int, args...)
    # allocate a buffer
    ex = quote
        buf = cudaGetParameterBuffer(f, blocks, threads, shmem)
    end

    # store the parameters
    #
    # > Each individual parameter placed in the parameter buffer is required to be aligned.
    # > That is, each parameter must be placed at the n-th byte in the parameter buffer,
    # > where n is the smallest multiple of the parameter size that is greater than the
    # > offset of the last byte taken by the preceding parameter. The maximum size of the
    # > parameter buffer is 4KB.
    offset = 0
    for i in 1:length(args)
        buf_index = Base.ceil(Int, offset / sizeof(args[i])) + 1
        offset = buf_index * sizeof(args[i])
        push!(ex.args, :(
            unsafe_store!(Base.unsafe_convert(Ptr{$(args[i])}, buf), args[$i], $buf_index)
        ))
    end

    push!(ex.args, :(return buf))

    return ex
end

if VERSION >= v"1.2.0-DEV.512"
    @inline synchronize() = ccall("extern cudaDeviceSynchronize", llvmcall, Cint, ())
else
    @eval @inline synchronize() = Base.llvmcall(
        $("declare i32 @cudaDeviceSynchronize()",
          "%rv = call i32 @cudaDeviceSynchronize()
           ret i32 %rv"), cudaError_t, Tuple{})
end

"""
    synchronize()

Wait for the device to finish. This is the device side version,
and should not be called from the host.

`synchronize` acts as a synchronization point for
child grids in the context of dynamic parallelism.
"""
synchronize


text:
include("../src/Groebner.jl")

if !isdefined(Main, :Groebner)
    import Groebner
end

import AbstractAlgebra
using BenchmarkTools
using Logging
global_logger(ConsoleLogger(stderr, Logging.Error))

# BenchmarkTools.DEFAULT_PARAMETERS.seconds = 6
BenchmarkTools.DEFAULT_PARAMETERS.samples = 4

function benchmark_system_my(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    # gb = Groebner.groebner(system)
    # println("length = $(length(gb))")
    # println("degree = $(maximum(AbstractAlgebra.total_degree, gb))")

    @btime Groebner.groebner($system, reduced=false)

end

function run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
    ("cyclic 7", reverse(Groebner.cyclicn(7, ground=ground))),
    ("cyclic 8", reverse(Groebner.cyclicn(8, ground=ground))),
    ("root 12", reverse(Groebner.rootn(12, ground=ground))),
    ("root 13", reverse(Groebner.rootn(13, ground=ground))),
    #("root 14", Groebner.rootn(14, ground=ground)),
        ("katsura 9", reverse(Groebner.katsuran(9, ground=ground))),
        ("noon 7"    ,Groebner.noonn(7, ground=ground)),
        ("noon 8"    ,Groebner.noonn(8, ground=ground)),
        ("eco 10"    ,Groebner.eco10(ground=ground))
    ]

    for (name, system) in systems
        println("$name")
        benchmark_system_my(system)
    end
end

println()
ground = AbstractAlgebra.GF(2^31 - 1)
run_f4_ff_degrevlex_benchmarks(ground)

#=
cyclic 7
  216.897 ms (96622 allocations: 70.60 MiB)
cyclic 8
  4.521 s (435192 allocations: 380.99 MiB)
root 12
  74.583 ms (179896 allocations: 32.40 MiB)
root 13
  242.372 ms (451340 allocations: 67.87 MiB)
katsura 9
  761.210 ms (123510 allocations: 59.33 MiB)
noon 7
  248.028 ms (335407 allocations: 86.89 MiB)
noon 8
  1.960 s (1501606 allocations: 444.08 MiB)
eco 10
  191.307 ms (119125 allocations: 41.66 MiB)
=#


text:
##
# Sparse BroadcastStyle
##
BroadcastStyle(::StructuredMatrixStyle{<:Diagonal}, ::BandedBlockBandedStyle) =
    BandedBlockBandedStyle()
BroadcastStyle(::BandedBlockBandedStyle, ::StructuredMatrixStyle{<:Diagonal}) =
    BandedBlockBandedStyle()

function blockbandwidths(P::PseudoBlockMatrix{<:Any,<:Diagonal})
    blockisequal(axes(P,1),axes(P,2)) || throw(DimensionMismatch())
    (0,0)
end

bandeddata(P::PseudoBlockMatrix) = bandeddata(P.blocks)
bandwidths(P::PseudoBlockMatrix) = bandwidths(P.blocks)

BroadcastStyle(::Type{<:SubArray{<:Any,2,<:PseudoBlockMatrix{<:Any,<:Diagonal},
                                <:Tuple{<:BlockSlice1,<:BlockSlice1}}}) = BandedStyle()

const BlockDiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Diagonal{VT}}

BlockDiagonal(A) = mortar(Diagonal(A))

function sizes_from_blocks(A::Diagonal, _)
    # for k = 1:length(A.du)
    #     size(A.du[k],1) == sz[1][k] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.du[k],2) == sz[2][k+1] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.dl[k],1) == sz[1][k+1] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    #     size(A.dl[k],2) == sz[2][k] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    # end
    (size.(A.diag, 1), size.(A.diag,2))
end

# Block Bi/Tridiagonal
const BlockTridiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Tridiagonal{VT}}
const BlockBidiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Bidiagonal{VT}}

BlockTridiagonal(A,B,C) = mortar(Tridiagonal(A,B,C))
BlockBidiagonal(A, B, uplo) = mortar(Bidiagonal(A,B,uplo))

function sizes_from_blocks(A::Tridiagonal, _)
    # for k = 1:length(A.du)
    #     size(A.du[k],1) == sz[1][k] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.du[k],2) == sz[2][k+1] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.dl[k],1) == sz[1][k+1] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    #     size(A.dl[k],2) == sz[2][k] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    # end
    (size.(A.d, 1), size.(A.d,2))
end

function sizes_from_blocks(A::Bidiagonal, _)
    # for k = 1:length(A.du)
    #     size(A.du[k],1) == sz[1][k] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.du[k],2) == sz[2][k+1] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.dl[k],1) == sz[1][k+1] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    #     size(A.dl[k],2) == sz[2][k] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    # end
    (size.(A.dv, 1), size.(A.dv,2))
end

blockbandwidths(A::BlockArray) = bandwidths(A.blocks)
isblockbanded(A::BlockArray) = isbanded(A.blocks)

@inline function getblock(block_arr::BlockBidiagonal{T,VT}, K::Int, J::Int) where {T,VT<:AbstractMatrix}
    @boundscheck blockcheckbounds(block_arr, K, J)
    l,u = blockbandwidths(block_arr)
    -l ≤ (J-K) ≤ u || return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
    block_arr.blocks[K,J]
end

@inline function getblock(block_arr::BlockTridiagonal{T,VT}, K::Int, J::Int) where {T,VT<:AbstractMatrix}
    @boundscheck blockcheckbounds(block_arr, K, J)
    abs(J-K) ≥ 2 && return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
    block_arr.blocks[K,J]
end

checksquareblocks(A) = blockisequal(axes(A)...) || throw(DimensionMismatch("blocks are not square: block dimensions are $(axes(A))"))

for op in (:-, :+)
    @eval begin
        function $op(A::BlockDiagonal, λ::UniformScaling)
            checksquareblocks(A)
            mortar(Diagonal(broadcast($op, A.blocks.diag, Ref(λ))))
        end
        function $op(λ::UniformScaling, A::BlockDiagonal)
            checksquareblocks(A)
            mortar(Diagonal(broadcast($op, Ref(λ), A.blocks.diag)))
        end

        function $op(A::BlockTridiagonal, λ::UniformScaling)
            checksquareblocks(A)
            mortar(Tridiagonal(A.blocks.dl, broadcast($op, A.blocks.d, Ref(λ)), A.blocks.du))
        end
        function $op(λ::UniformScaling, A::BlockTridiagonal)
            checksquareblocks(A)
            mortar(Tridiagonal(broadcast($op,A.blocks.dl), broadcast($op, Ref(λ), A.blocks.d), broadcast($op,A.blocks.du)))
        end
        function $op(A::BlockBidiagonal, λ::UniformScaling)
            checksquareblocks(A)
            mortar(Bidiagonal(broadcast($op, A.blocks.dv, Ref(λ)), A.blocks.ev, A.blocks.uplo))
        end
        function $op(λ::UniformScaling, A::BlockBidiagonal)
            checksquareblocks(A)
            mortar(Bidiagonal(broadcast($op, Ref(λ), A.blocks.dv), broadcast($op,A.blocks.ev), A.blocks.uplo))
        end
    end
end


text:
# transportMatrixTests.jl

# Age test, should be ≈ 1000 years
function computeIdealMeanAge(wet3d,T)
    l3d = copy(wet3d)
    l3d[:,:,2:end] .= 0
    iwet = find(wet3d)
  L = spdiagm(l3d[iwet])
    nwet = length(iwet)
  return (T + L) \ ones(nwet)
end
meanΓyrs = mean(computeIdealMeanAge(wet3d,T)) / (365 * spd)
@test 100 < meanΓyrs < 1000

# T 1 ≈ 0 test
Tcolmax = vec(maximum(abs.(T),1))
Trowmax = vec(maximum(abs.(T),2))
T1 = T * ones(nwet)
maximum(abs.(T1) ./ Tcolmax)
maximum(abs.(T1) ./ Trowmax)

# vᵀ T ≈ 0 test
vᵀT = T.' * v
maximum(abs.(vᵀT) ./ (v .* Tcolmax))
maximum(abs.(vᵀT) ./ (v .* Trowmax))

# S 1 ≈ 0 test
Scolmax = vec(maximum(abs.(S),1))
Srowmax = vec(maximum(abs.(S),2))
S1 = S * ones(nwet)
maximum(abs.(S1) ./ Scolmax)
maximum(abs.(S1) ./ Srowmax)
S13d = NaN * wet3d
S13d[iwet] .= S1
using Plots
gr()
plot(S13d[50,54,:],line = :scatter)

# vᵀ S ≈ 0 test
vᵀS = S.' * v
maximum(abs.(vᵀS) ./ (v .* Scolmax))
maximum(abs.(vᵀS) ./ (v .* Srowmax))

function isMassConserving(x,v)
  vᵀx = v.' * x
  return vᵀx


text:
@testset "compactify $Q" for Q in (Rational, FastRational)
    mid = Q(100, 111)
    tol = Q(1, 11)
    res = Q(5, 6)
    @test compactify(low=mid - tol, high=mid + tol) == res
    @test compactify(mid, tol) == res
    @test compactify(mid, eps()) == mid
    @test_throws ArgumentError compactify(mid, eps()/3)
    tm = typemax(Int64)>>12
    low = Q(tm, tm-1)
    high = Q(tm+1, tm-1)
    mid = Q(2tm+1, 2tm-2)
    tol = Q(1, 2tm-2)
    @test compactify(mid, tol) == Q(tm>>1+2, tm>>1+1) # rounding error effect?
    @test compactify(low=low, high=high) == Q(tm>>1+1, tm>>1) # better result
    #                                                                          >>>> !!FIXME!! throws MethodError
    # @test_throws ArgumentError compactify(mid, Q(0))
    @test compactify(low=high, high=high) == high
end


text:
@testset "Agent-based modelling" begin
    @enum Status susceptible infected recovered

    # Custom atom type
    mutable struct Person
        i::Int
        status::Status
        mass::Float64
        σ::Float64
        ϵ::Float64
    end

    Molly.mass(person::Person) = person.mass

    # Custom PairwiseInteraction
    struct SIRInteraction <: PairwiseInteraction
        nl_only::Bool
        dist_infection::Float64
        prob_infection::Float64
        prob_recovery::Float64
    end

    # Custom Logger
    struct SIRLogger
        n_steps::Int
        fracs_sir::Vector{Vector{Float64}}
    end

    # Custom force function
    function Molly.force(inter::SIRInteraction, dr, coord_i, coord_j, atom_i, atom_j, box_size)
        if (atom_i.status == infected && atom_j.status == susceptible) ||
                    (atom_i.status == susceptible && atom_j.status == infected)
            # Infect close people randomly
            r2 = sum(abs2, dr)
            if r2 < inter.dist_infection ^ 2 && rand() < inter.prob_infection
                atom_i.status = infected
                atom_j.status = infected
            end
        end
        # Workaround to obtain a self-interaction
        if atom_i.i == (atom_j.i + 1)
            # Recover randomly
            if atom_i.status == infected && rand() < inter.prob_recovery
                atom_i.status = recovered
            end
        end
        return zero(coord_i)
    end

    # Custom logging function
    function Molly.log_property!(logger::SIRLogger, s, neighbors, step_n; parallel=true)
        if step_n % logger.n_steps == 0
            counts_sir = [
                count(p -> p.status == susceptible, s.atoms),
                count(p -> p.status == infected   , s.atoms),
                count(p -> p.status == recovered  , s.atoms)
            ]
            push!(logger.fracs_sir, counts_sir ./ length(s))
        end
    end

    n_people = 500
    n_steps = 1_000
    box_size = SVector(10.0, 10.0)
    temp = 1.0
    n_starting = 2
    atoms = [Person(i, i <= n_starting ? infected : susceptible, 1.0, 0.1, 0.02) for i in 1:n_people]
    coords = place_atoms(n_people, box_size, 0.1)
    velocities = [velocity(1.0, temp; dims=2) for i in 1:n_people]
    pairwise_inters = (
        LennardJones=LennardJones(nl_only=true),
        SIR=SIRInteraction(false, 0.5, 0.06, 0.01),
    )
    neighbor_finder = DistanceNeighborFinder(
        nb_matrix=trues(n_people, n_people),
        n_steps=10,
        dist_cutoff=2.0,
    )
    simulator = VelocityVerlet(dt=0.02, coupling=AndersenThermostat(temp, 5.0))

    s = System(
        atoms=atoms,
        pairwise_inters=pairwise_inters,
        coords=coords,
        velocities=velocities,
        box_size=box_size,
        neighbor_finder=neighbor_finder,
        loggers=Dict(
            "coords" => CoordinateLogger(Float64, 10; dims=2),
            "SIR"    => SIRLogger(10, []),
        ),
        force_units=NoUnits,
        energy_units=NoUnits,
    )

    @time simulate!(s, simulator, n_steps; parallel=false)
end


text:
function mp_evaluate(equations, x)
    variables = MP.variables(equations)
    map(equations) do f
        f(variables => x)
    end
end

function mp_jacobian(equations, x)
    variables = MP.variables(equations)
    map(Iterators.product(equations, variables)) do arg
        (f, var) = arg
        MP.differentiate(f, var)(variables => x)
    end
end

function sp_evaluate(equations, x)
    F = SP.PolynomialSystem(equations)
    SP.evaluate(F, x)
end
function sp_jacobian(equations, x)
    F = SP.PolynomialSystem(equations)
    SP.jacobian(F, x)
end

function katsura(n)
  @polyvar x[0:n] # This creates variables x0, x1, ...

  return [
    (sum(x[abs(l)+1]*x[abs(m-l)+1] for l=-n:n if abs(m-l)<=n) -
    x[m+1] for m=0:n-1)...,
    x[1] + 2sum(x[i+1] for i=1:n) - 1
  ]
end

@testset "Katsura Evaluation tests" begin
    for T = [Float64, Complex{Float64}]
        for n=4:12
            nvars = n+1
            eqs = katsura(n)
            x = SVector{nvars}(rand(T, nvars))

            @test norm(sp_evaluate(eqs, x) - mp_evaluate(eqs, x)) < 1e-14
            @test norm(sp_evaluate(eqs, Vector(x)) - mp_evaluate(eqs, x)) < 1e-14

            @test norm(sp_jacobian(eqs, x) - mp_jacobian(eqs, x)) < 1e-14
            @test norm(sp_jacobian(eqs, Vector(x)) - mp_jacobian(eqs, x)) < 1e-14

            F = SP.PolynomialSystem(eqs)
            val, jac = SP.evaluate_and_jacobian(F, x)
            @test norm(val - mp_evaluate(eqs, x)) < 1e-14
            @test norm(F(x) - mp_evaluate(eqs, x)) < 1e-14
            @test norm(jac - mp_jacobian(eqs, x)) < 1e-14

            val, jac = SP.evaluate_and_jacobian(F, Vector(x))
            @test norm(val - mp_evaluate(eqs, x)) < 1e-14
            @test norm(F(x) - mp_evaluate(eqs, x)) < 1e-14
            @test norm(jac - mp_jacobian(eqs, x)) < 1e-14
        end
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

@testset "util" begin
    @testset "BAT._car_cdr" begin
        @test BAT._car_cdr_impl() == ()
        @test BAT._car_cdr_impl(4, 3.5, -1) == (4, (3.5, -1))
        @test BAT._car_cdr((3, -8, 2)) == (3, (-8, 2))
    end

    @testset "BAT._all_lteq" begin
        @test BAT._all_lteq(-0.4, -0.4:0.4:0.8 ,0.8)
        @test !BAT._all_lteq(-0.4, -0.4:0.4:0.8 ,0.7)
        @test !BAT._all_lteq(-0.3, -0.4:0.4:0.8 ,0.8)

        @test BAT._all_lteq([0.0, 0.5, -0.4], [0.3, 0.5, -0.4] , [0.3, 0.6, -0.3])
        @test !BAT._all_lteq([0.5, -0.4], [0.49, -0.4] , [0.6, -0.3])
        @test_throws DimensionMismatch BAT._all_lteq(
            [0.5, -0.4], [0.3, 0.5, -0.4] , [0.3, 0.6, -0.3])
    end
end


text:
@system Weather begin
    calendar(context) ~ ::Calendar(override)
    vp(context): vapor_pressure ~ ::VaporPressure

    s: store ~ provide(init=calendar.time, parameter)

    solrad: solar_radiation ~ drive(from=s, by=:SolRad, u"W/m^2")

    CO2 => 400 ~ preserve(u"μmol/mol", parameter)

    RH: relative_humidity ~ drive(from=s, by=:RH, u"percent")
    #RH => 0.6 ~ track # 0~1

    T_air: air_temperature ~ drive(from=s, by=:Tair, u"°C")
    #T_air => 25 ~ track # C

    Tk_air(T_air): absolute_air_temperature ~ track(u"K")

    wind: wind_speed ~ drive(from=s, by=:Wind, u"m/s")
    #wind => 2.0 ~ track # meters s-1

    #TODO: make P_air parameter?
    P_air: air_pressure => 100 ~ track(u"kPa")

    VPD(T_air, RH, D=vp.D): vapor_pressure_deficit => D(T_air, RH) ~ track(u"kPa")
    VPD_Δ(T_air, Δ=vp.Δ): vapor_pressure_saturation_slope_delta => Δ(T_air) ~ track(u"kPa/K")
    VPD_s(T_air, P_air, s=vp.s): vapor_pressure_saturation_slope => s(T_air, P_air) ~ track(u"K^-1")
end

#TODO: make @stub macro to automate this
@system WeatherStub begin
    weather ~ hold

    vp(x=weather.vp) => x ~ ::VaporPressure

    solrad(weather.solrad): solar_radiation ~ track(u"W/m^2")
    CO2(weather.CO2) ~ track(u"μmol/mol")
    RH(weather.RH): relative_humidity ~ track(u"percent")
    T_air(weather.T_air): air_temperature ~ track(u"°C")
    Tk_air(weather.Tk_air): absolute_air_temperature ~ track(u"K")
    wind(weather.wind): wind_speed ~ track(u"m/s")
    P_air(weather.P_air): air_pressure ~ track(u"kPa")

    VPD(weather.VPD): vapor_pressure_deficit ~ track(u"kPa")
    VPD_Δ(weather.VPD_Δ): vapor_pressure_saturation_slope_delta ~ track(u"kPa/K")
    VPD_s(weather.VPD_s): vapor_pressure_saturation_slope ~ track(u"K^-1")
end


text:
using AtBackslash: @\
using DataFrames
using Glob
using NamedTupleTools: @namedtuple
using Statistics
using ThreadsXBenchmarks
using UnPack: @unpack
using VegaLite

loadall(datadir) = (
    sort_target = copy(
        DataFrame,
        ThreadsXBenchmarks.Loading.results_by_benchname(
            "sort",
            readdir(glob"scaling_nthreads-*.json", datadir),
        ),
    ),
    sort_baseline = copy(
        DataFrame,
        ThreadsXBenchmarks.Loading.results_by_benchname(
            "sort",
            readdir(glob"scaling_nthreads_baseline.json", datadir),
        ),
    ),
    fold_target = copy(
        DataFrame,
        ThreadsXBenchmarks.Loading.results_by_benchname(
            !=("sort"),
            readdir(glob"scaling_nthreads-*.json", datadir),
        ),
    ),
    fold_baseline = copy(
        DataFrame,
        ThreadsXBenchmarks.Loading.results_by_benchname(
            !=("sort"),
            readdir(glob"scaling_nthreads_baseline.json", datadir),
        ),
    ),
)

function joinbaseline(target, baseline, on, agg)
    baseline_agg = by(baseline, on) do group
        @assert allunique(group.itime)
        # group = group[group.itime.>1, :]
        (baseline_time = agg(group.time),)
    end

    normalized = join(target, baseline_agg, on = on)
    normalized[!, :normalized_time] = normalized.time ./ normalized.baseline_time
    normalized[!, :speedup] = 1 ./ normalized.normalized_time

    return normalized
end

plotall(datadir::AbstractString; kw...) = plotall(loadall(datadir); kw...)

function plotall(data; agg = median)
    @unpack sort_target, sort_baseline, fold_target, fold_baseline = data

    sort_normalized =
        joinbaseline(sort_target, sort_baseline, [:distribution, :datasize], agg)
    fold_normalized = joinbaseline(fold_target, fold_baseline, [:benchname, :datasize], agg)

    fold_normalized = by(fold_normalized, [:benchname]) do group
        datasizes = sort!(unique(group.datasize))
        df = DataFrame(group)
        df[!, :datasize_label] = findfirst.((==).(group.datasize), Ref(datasizes))
        df
    end

    # Smaller `datasize` look better:
    sort_normalized = sort_normalized[sort_normalized.datasize.==1_000_000, :]

    # Larger `datasize` look better:
    fold_normalized = fold_normalized[fold_normalized.datasize_label.==2, :]

    sort_raw_plot = nothing
    sort_raw_plot = @vlplot(
        mark = {type = :point, tooltip = true},
        x = :nthreads,
        y = :speedup,
        color = :alg,
        row = :distribution,
        column = :datasize,
        data = sort_normalized,
    )

    sort_agg_plot = @vlplot(
        facet = {
            column = {field = :distribution, type = :nominal},
            # row = {field = :datasize, type = :nominal},
        },
        spec = {
            layer = [
                {
                    mark = {type = :point, tooltip = true},
                    encoding = {
                        x = {field = :nthreads},
                        y = {field = :speedup, axis = {title = "$agg(speedup)"}},
                        color = {field = :alg},
                    },
                },
                {
                    mark = {type = :rule, color = :black},
                    data = {values = [{ref = 1}]},
                    encoding = {y = {field = :ref, type = :quantitative}},
                },
            ],
        },
        data = by(
            sort_normalized,
            [:nthreads, :basesize, :alg, :distribution, :datasize];
            speedup = :speedup => agg,
        ),
    )

    fold_raw_plot = @vlplot(
        mark = {type = :point, tooltip = true},
        x = :nthreads,
        y = :speedup,
        color = :benchname,
        row = :benchname,
        column = :datasize,
        data = fold_normalized,
    )

    fold_agg_plot = @vlplot(
        layer = [
            {
                mark = {type = :point, tooltip = true},
                encoding = {
                    x = {field = :nthreads},
                    y = {field = :speedup, axis = {title = "$agg(speedup)"}},
                    color = {field = :benchname},
                },
            },
            {
                mark = {type = :rule, color = :black},
                data = {values = [{ref = 1}]},
                encoding = {y = {field = :ref, type = :quantitative}},
            },
        ],
        data = by(fold_normalized, [:nthreads, :benchname]; speedup = :speedup => agg),
    )

    return @namedtuple(data..., sort_raw_plot, sort_agg_plot, fold_raw_plot, fold_agg_plot)
end


text:
# LS-RTM of the 2D Marmousi model using stochastic gradient descent
# Author: pwitte.slim@gmail.com
# Date: December 2018
#
# Warning: The examples requires ~40 GB of memory per shot if used without optimal checkpointing.
#

using Statistics, Random, LinearAlgebra
using JUDI, SegyIO, HDF5, PyPlot

# Load migration velocity model
if ~isfile("$(JUDI.JUDI_DATA)/marmousi_model.h5")
    ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_model.h5")
end
n, d, o, m0 = read(h5open("$(JUDI.JUDI_DATA)/marmousi_model.h5", "r"), "n", "d", "o", "m0")

# Set up model structure
model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)

# Load data
if ~isfile("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
    ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_2D.segy")
end
block = segy_read("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
d_lin = judiVector(block)   # linearized observed data

# Set up wavelet
src_geometry = Geometry(block; key = "source", segy_depth_key = "SourceDepth")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.03)    # 30 Hz wavelet
q = judiVector(src_geometry, wavelet)

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

# Setup operators
opt = Options(optimal_checkpointing=true)  # ~40 GB of memory per source w/o checkpointing
M = judiModeling(model0, q.geometry, d_lin.geometry; options=opt)
J = judiJacobian(M, q)

# Right-hand preconditioners (model topmute)
Mr = judiTopmute(model0.n, 52, 10)

# Stochastic gradient
x = zeros(Float32, prod(model.n))
batchsize = 10
niter = 20
fval = zeros(Float32, niter)

# Main loop
for j = 1: niter
    println("Iteration: ", j)

    # Select batch and set up left-hand preconditioner
    i = randperm(d_lin.nsrc)[1: batchsize]
    Ml = judiMarineTopmute2D(30, d_lin[i].geometry)    # data topmute

    # Compute residual and gradient
    r = Ml*J[i]*Mr*x - Ml*d_lin[i]
    g = adjoint(Mr)*adjoint(J[i])*adjoint(Ml)*r

    # Step size and update variable
    fval[j] = .5f0*norm(r)^2
    t = norm(r)^2/norm(g)^2
    global x -= t*g
end

# Save final velocity model, function value and history
h5open("lsrtm_marmousi_sgd_result.h5", "w") do file
    write(file, "x", reshape(x, model0.n), "fval", fval)
end

# Plot final image
figure(); imshow(copy(adjoint(reshape(x, model0.n))), extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
title("LS-RTM with SGD"); xlabel("Lateral position [km]"); ylabel("Depth [km]")


text:
import Rotations
import StaticArrays: StaticVector, @SVector
import PoseComposition: Pose
import LinearAlgebra: I
import Base: convert

function Pose(orientation::Rotations.Rotation{3})::Pose
    return Pose(zeros(3), orientation)
end

const IDENTITY_ORIENTATION = Rotations.UnitQuaternion(1.0,0.0,0.0,0.0)

function Pose(pos::AbstractVector{<:Real})::Pose
    return Pose(pos, IDENTITY_ORIENTATION)
end

function Pose(x::Real, y::Real, z::Real)::Pose
    return Pose([x,y,z], IDENTITY_ORIENTATION)
end

const IDENTITY_POSE = Pose(zeros(3),IDENTITY_ORIENTATION)

rotation_pose(pose::Pose) = Pose(zeros(3), pose.orientation)
translation_pose(pose::Pose) = Pose(pose.pos, IDENTITY_ORIENTATION)

"""
    c_relative_to_a::Pose = get_c_relative_to_a(b_relative_to_a::Pose, c_relative_to_b::Pose)

Given the pose of b relative to a, and the pose of c relative to b, return the pose of c relative to a.
"""
function get_c_relative_to_a(b_relative_to_a::Pose, c_relative_to_b::Pose)
    pos = b_relative_to_a.pos + (b_relative_to_a.orientation * c_relative_to_b.pos)
    orientation = b_relative_to_a.orientation * c_relative_to_b.orientation
    return Pose(pos, orientation)
end

function pose_pow(pose::Pose, t::Int)
    if t == 0
        return IDENTITY_POSE
    end
    p = pose
    for i in 1:(t-1)
        p = get_c_relative_to_a(p, pose)
    end
    p
end

"""
    b_relative_to_a::Pose = inverse_pose(a_relative_to_b::Pose)

Given the pose of b relative to a, return the pose of a relative to b.
"""
function inverse_pose(a_relative_to_b::Pose)
    pos = a_relative_to_b.orientation' * (-a_relative_to_b.pos)
    orientation = inv(a_relative_to_b.orientation)
    b_relative_to_a = Pose(pos, orientation)
    return b_relative_to_a
end

"""
    points = get_points_in_frame_b(points_in_frame_a::Matrix{<:Real}, b_relative_to_a::Pose)

Given coordinates of points in one coordinate frame, get the coordinates of the points in another coordinate frame.

The points that are taken as input and returned as output are both 3 x n matrices.
"""
function get_points_in_frame_b(points_in_frame_a::Matrix{<:Real}, b_relative_to_a::Pose)
    if size(points_in_frame_a)[1] != 3
        error("expected an 3 x n matrix")
    end
    n = size(points_in_frame_a[2])
    return b_relative_to_a.orientation' * (points_in_frame_a .- b_relative_to_a.pos)
end

function move_points_to_frame_b(points_in_frame_a::Matrix{<:Real}, b_relative_to_a::Pose)
    if size(points_in_frame_a)[1] != 3
        error("expected an 3 x n matrix")
    end
    n = size(points_in_frame_a[2])
    return b_relative_to_a.orientation * points_in_frame_a .+ b_relative_to_a.pos
end

function convert(::Type{Matrix{Float32}}, pose::Pose)
    R = Matrix{Float32}(pose.orientation)
    mat = Matrix{Float32}(I, 4, 4)
    mat[1:3,1:3] = R
    mat[1:3,4] = pose.pos
    return mat
end

export Pose, close, isapprox, IDENTITY_ORIENTATION, IDENTITY_POSE, inverse_pose, get_points_in_frame_b,
        get_c_relative_to_a, move_points_to_frame_b, pose_pow, rotation_pose, translation_pose


text:
module TestSymSpell
include("preamble.jl")
using SymSpellChecker: DictWord, String32

@testset "dictionary creation" begin
    d = SymSpell()

    @test push!(d, "a", 2)
    @test d.words == [DictWord("a", 2)]
    @test isempty(d.below_threshold_words)
    @test d.deletes == Dict{String, Vector{UInt32}}("" => [1], "a" => [1])

    push!(d, "bc", 2)
    @test d.max_length == 2

    d = SymSpell(count_threshold = 2)
    push!(d, "a", 100)
    push!(d, "bc", 100)
    push!(d, "c", 1)
    @test length(d.words) == 2
    @test length(d.below_threshold_words) == 1

    push!(d, "c", 1)
    @test length(d.words) == 3
    @test length(d.below_threshold_words) == 0

    d = SymSpell(max_dictionary_edit_distance = 1)
    push!(d, "abc", 100)
    @test d.deletes == Dict{String, Vector{UInt32}}("bc" => [1], "ac" => [1], "ab" => [1], "abc" => [1])

    push!(d, "abd", 100)
    @test d.deletes["ab"] == UInt32[1, 2]

    d = SymSpell(max_dictionary_edit_distance = 1)
    push!(d, "world", 10)
    push!(d, "word", 5)
    @test d.words[d.words_idx["word"]] == DictWord("word", 5)
end

@testset "utf-8 dictionary" begin
    d = SymSpell()

    @test push!(d, "привет", 10)
end

@testset "dictionary load" begin
    d = SymSpell(count_threshold = 2)
    update!(d, joinpath(@__DIR__, "..", "assets", "test_dict.txt"))

    @test d.words[d.words_idx["key"]] == DictWord("key", 10)
    @test d.below_threshold_words["sad"] == 1
    @test length(d.words) == 3
end

end # module


text:
using LinearAlgebra, SparseArrays, Random, QPDAS
# Number of variables
n = 1000
# Number of halfspaces
mi = 50 # Inequalities
me = 50 # Equalities

Random.seed!(1)
# One point in polytope
x0 = randn(n)

# Create polytope containing x0
# Inequality
C = Matrix{Float64}(undef, mi, n)
d = randn(mi)

# Make sure x0 is in polytope by setting sign of inequality
for i = 1:mi
    v = randn(n)
    b = randn()
    if v'x0  <= b
        C[i,:] .= v
    else
        C[i,:] .= -v
    end
    d[i] = b
end

# Create equality
A = randn(me, n)
b = A*x0

###### Solve with only inequality

# Project from here
x = randn(n)

@time qp = QuadraticProgram(zeros(0,n), zeros(0), C, d, -x, I)
@time sol, val = solve!(qp)
# Test Feasibility
maximum(C*sol-d)

###### Solve with equality and inequality

# Project from here
x = randn(n)
qp = QuadraticProgram(A, b, C, d, -x, I)
solve!(qp)
# Test Feasibility
maximum(C*sol-d)

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
# This work is extending the Svyatoslav Pidgorny's work https://github.com/SP3269/Julia-Playground
# I juste structured it and encapsulated operations and added some functions

module BfInterpreter

    #using ResumableFunctions

    export getInstructionsDict, bfType, runnable!, execute, bfVM, floatToBfInstruction, BracketError, MermoryBfVM, genesToBfInstructions

    struct BracketError <: Exception
        message :: String
    end

    struct MermoryBfVM <: Exception
        message :: String
    end

    mutable struct bfVM
        cells :: Array{UInt8,1}         # default value : [0 for i=1:15000]
        nbCells :: Int
        cellno :: Int64                 # Current cell number :: initial value : 1
        n :: Int64                      # length of the bf code
        m_InstructionPointer :: Int64   # Pointer to the current position in the code, initial_value : 1
        nbExecLimit :: Int32            # default value : 200000
        nbExec :: Int32                 # initial value : 0
        storage :: UInt8
        instructionDict #:: AbstractDict
        ignoreJump :: Bool
        ignoreComment :: Bool
        output :: Array{Char,1}         # Data returned
        input #:: String
        indexInput :: Int
        l_input :: Int
        m_Stop :: Bool

        ## <summary>
        ## Holds the instruction pointer for the start of the loop. Used to bypass all inner-loops when searching for the end of the current loop.
        ## </summary>
        m_ExitLoopInstructionPointer :: Int # = 0

        m_CurrentCallStack :: Array

        #bfVM(n :: Int64) = new([0 for i=1:30000], 1, n, 1, 200000, 0)
    end

    function getInstructionsDict()
        instructionsSet = Dict()

        # Create the instruction set for Basic Brainfuck.
        instructionsSet['+'] = vm ->
            if (!vm.ignoreJump && !vm.ignoreComment)
                #@assert vm.cells[ vm.cellno] < 255
                if vm.cells[ vm.cellno] == 255
                    vm.cells[ vm.cellno] = 0 #vm.cells[ vm.cellno]+1
                    #throw(MermoryBfVM("Memory: 0 => 255 (cell $(vm.cellno), state $(vm.cells[vm.cellno]), cmd '+')"))
                else
                    vm.cells[ vm.cellno] +=  1
                end
                #vm.cells[ vm.cellno] = vm.cells[ vm.cellno]+1
            end
        instructionsSet['-'] = vm ->
            if (!vm.ignoreJump && !vm.ignoreComment)
                #@assert vm.cells[ vm.cellno] < 255
                if vm.cells[ vm.cellno] == 0
                    vm.cells[ vm.cellno] = 255 #vm.cells[ vm.cellno]+1
                    #throw(MermoryBfVM("Memory: 0 => 255 (cell $(vm.cellno), state $(vm.cells[vm.cellno]), cmd '+')"))
                else
                    vm.cells[ vm.cellno] -=  1
                end
                #vm.cells[ vm.cellno] = vm.cells[ vm.cellno]+1
            end
        instructionsSet['>'] = vm -> ((!vm.ignoreJump && !vm.ignoreComment) ?
            #=lcells = length(vm.cells)
            #@assert vm.cellno < lcells
            if vm.cellno == lcells
                throw(MermoryBfVM("Memory: 1 => 1 ($(vm.cells[1])) , $lcells => 1 ($(vm.cells[lcells]))"))
            end=#
            vm.cellno = vm.cellno + 1
            : nothing)
        instructionsSet['<'] = vm ->
            if !vm.ignoreJump && !vm.ignoreComment
                #=if vm.cellno == 1
                    lcells = vm.nbCells #length(vm.cells)
                    throw(MermoryBfVM("Memory: 1 => 1 ($(vm.cells[1])) , $lcells => 1 ($(vm.cells[lcells]))"))
                end=#
                #@show vm.cellno
                vm.cellno = vm.cellno - 1
            end
        instructionsSet['.'] = vm -> ((!vm.ignoreJump && !vm.ignoreComment) ?
            push!(vm.output, Char(vm.cells[vm.cellno]) )
            : nothing)
        #(@yield Char(vm.cells[vm.cellno]) ) : nothing
        instructionsSet[','] = vm ->
            if !vm.ignoreJump && !vm.ignoreComment
                if vm.indexInput <= vm.l_input # length(vm.input)
                    vm.cells[ vm.cellno] = UInt8(vm.input[vm.indexInput])
                    vm.indexInput += 1
                else
                    vm.cells[vm.cellno] = 0
                end
            end
        instructionsSet['['] = vm ->
            if !vm.ignoreComment
                if (!vm.ignoreJump && vm.cells[vm.cellno] == 0)
                    ## Jump forward to the matching ] and exit this loop (skip over all inner loops).
                    vm.ignoreJump = true;

                    ## Remember this instruction pointer, so when we get past all inner loops and finally pop this one off the stack, we know we're done.
                    vm.m_ExitLoopInstructionPointer = vm.m_InstructionPointer;
                end
                pushfirst!(vm.m_CurrentCallStack, vm.m_InstructionPointer)
                #this.m_CurrentCallStack.Push(vm.m_InstructionPointer);
            end
        instructionsSet[']'] = vm ->
            if !vm.ignoreComment
                temp = popfirst!(vm.m_CurrentCallStack) #vm.m_CurrentCallStack[end] #TODO

                if (!vm.ignoreJump)
                    vm.m_InstructionPointer = vm.cells[vm.cellno]!=0 ? temp -1 : vm.m_InstructionPointer
                elseif (temp == vm.m_ExitLoopInstructionPointer)
                    ## Continue executing after loop.
                    ## We've finally exited the loop.
                    vm.ignoreJump = false;
                    vm.m_ExitLoopInstructionPointer = 0;
                end
            end
        instructionsSet['$'] = vm ->  (!vm.ignoreJump && !vm.ignoreComment) ? (vm.storage = vm.cells[vm.cellno]) : nothing
        instructionsSet['!'] = vm ->  (!vm.ignoreJump && !vm.ignoreComment) ? (vm.cells[vm.cellno] = vm.storage) : nothing
        instructionsSet['@'] = vm -> (!vm.ignoreJump && !vm.ignoreComment) ? vm.m_Stop = true : nothing
        instructionsSet['%'] = vm -> (!vm.ignoreJump && !vm.ignoreComment) ? vm.cells[ vm.cellno] = trunc(UInt8, vm.cells[vm.cellno] % vm.storage) : nothing
        instructionsSet['/'] = vm -> (!vm.ignoreJump && !vm.ignoreComment) ? vm.cells[ vm.cellno] = trunc(UInt8, vm.cells[vm.cellno] / vm.storage) : nothing
        instructionsSet['*'] = vm -> (!vm.ignoreJump && !vm.ignoreComment) ? vm.cells[ vm.cellno] = trunc(UInt8, vm.cells[vm.cellno] * vm.storage) : nothing
        instructionsSet['#'] = vm -> (!vm.ignoreJump) ? (vm.ignoreComment = !vm.ignoreComment) : nothing

        instructionsSet
    end

    #=function execute(bfcode :: String, input :: Array{UInt8,1}, instructionSet)
        execute(bfcode, join(map((v -> Char(v)), input), ""), instructionSet)
    end =#

    function execute(bfcode :: String, instructionSet )
        execute(bfcode, "", instructionSet)
    end

    #=function execute(bfcode :: String, input :: Array{Char,1}, instructionSet )
        execute(bfcode, join(input, ""), instructionSet)
    end=#

    function execute(bfcode :: String, input , instructionSet )

        # initialize the vm
        vm = bfVM( [0 for i=1:2000], 2000, 1, length(bfcode),
                    0, 2000, -1 , 0, Dict(), false, false, [],
                     input, 1, length(input), false,  -1, [])

        #@show bfcode
        ignoreCharacter = vm.ignoreJump | vm.ignoreComment
        while vm.m_InstructionPointer< vm.n && vm.nbExec < vm.nbExecLimit &&  !vm.m_Stop
            vm.m_InstructionPointer = vm.m_InstructionPointer + 1
            vm.nbExec = vm.nbExec +1
            #global cellno
            #ignoreCharacter = ignoreJump | ignoreComment
            cmd = bfcode[ vm.m_InstructionPointer]
            try
                instructionSet[cmd](vm)
            catch y
                break
            end

            @debug "Executed $cmd at $m_InstructionPointer. Current cell # $(cellno)\t$(cells[cellno])"
        end
        vm.output, vm.nbExec
    end

    # Eliminating all extranious characters in the code
    function purifyCode(bfcode :: String)
        replace(bfcode, r"[^\+\-\<\>\.\,\[\]]" => s"")
    end
end


text:
import Base: isequal, convert

type KUparam{A,T,N}; arr; diff; init; initp; lr; gc; l1reg; l2reg; adagrad; ada; momentum; mom; nesterov; nes; average; avg; inc; KUparam()=new(); end

KUparam(w; o...)=init(setopt!(KUparam{atype(w),eltype(w),ndims(w)}(); arr=w, o...))
KUparam{A,T}(::Type{A}, ::Type{T}, d::Dims; o...)=KUparam(A(T,d); o...)
KUparam{A,T}(::Type{A}, ::Type{T}, d1::Int, d::Int...; o...)=KUparam(A,T,tuple(d1,d...); o...)
KUparam{T}(::Type{T}, d::Dims; o...)=KUparam((gpu()?CudaArray:Array),T,d; o...)
KUparam{T}(::Type{T}, d::Int...; o...)=KUparam(T, d; o...)
KUparam(d::Int...; o...)=KUparam(Float64, d; o...)
setopt!(p::KUparam; o...)=(for (n,v) in o; p.(n)=v; end; p)

# TODO: create an rgen type: makes things serializable
function init(p::KUparam, T::DataType=eltype(p), d::Dims=size(p.arr))
    (size(p.arr)==d && eltype(p.arr)==T) || (p.arr = similar(p.arr, T, d))
    # we want no init if params given in matrix
    if !isempty(p.arr) && nz(p,:init,nothing)
        initp = (!isdefined(p,:initp) ? () :
                 isa(p.initp,Tuple) ? p.initp :
                 (p.initp,))
        p.init(p.arr, initp...)
    end
    return p
end

function isequal(a::KUparam,b::KUparam)
    typeof(a) == typeof(b) || return false
    size(a) == size(b) || return false
    for n in fieldnames(a)
        if isdefined(a,n) && isdefined(b,n)
            isequal(a.(n), b.(n)) || return false
        elseif isdefined(a,n) || isdefined(b,n)
            return false
        end
    end
    return true
end

# BASIC ARRAY OPS:

for fname in (:eltype, :length, :ndims, :size, :strides, :pointer, :isempty, :vecnorm)
    @eval (Base.$fname)(a::KUparam)=$fname(a.arr)
end

for fname in (:size, :stride)
    @eval (Base.$fname)(a::KUparam,n)=$fname(a.arr,n)
end

atype{A}(::KUparam{A})=A
diff(a::KUparam)=a.diff
difnorm(a::KUparam)=(isdefined(a,:diff) ? vecnorm(a.diff) : 0)

# DEPRECATED:
# update!(::Nothing;o...)=nothing
# setopt!(::Nothing;o...)=nothing
# initdiff(w::KUparam; fill=nothing, o...)=(similar!(w, :diff, w.arr); fill!=nothing && fill!(w.diff,fill); w)

# We need to fix cpu/gpu copy so the type changes appropriately:

copy(x::KUparam)=deepcopy(x)

function cpucopy_internal{A<:CudaArray,T,N}(x::KUparam{A,T,N},d::ObjectIdDict)
    haskey(d,x) && return d[x]
    y = KUparam{Array,T,N}()
    for n in fieldnames(x)
        isdefined(x,n) || continue
        y.(n) = cpucopy_internal(x.(n),d)
    end
    d[x] = y
end

function gpucopy_internal{A<:Array,T,N}(x::KUparam{A,T,N},d::ObjectIdDict)
    haskey(d,x) && return d[x]
    y = KUparam{CudaArray,T,N}()
    for n in fieldnames(x)
        isdefined(x,n) || continue
        y.(n) = gpucopy_internal(x.(n),d)
    end
    d[x] = y
end

convert{A<:BaseArray}(::Type{A}, a::KUparam)=convert(A, a.arr)
convert{A<:BaseArray}(::Type{KUparam}, a::A)=KUparam(a)

# TODO: both weights and training parameters are called param; setopt! vs params is confusing
# Probably should rename this back to Param after moving others to DynamicArray.
# How does caffe deal with this problem?


text:
using Images, FileIO

img = load("data/bitmapOutputTest.jpg")
save("data/bitmapOutputTest.ppm", img)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function krest!(r)
    for side in [Nord,West,Sud,Ost]
        movemarkers!(r,side)
        moveback!(r,side)
    end
    putmarker!(r)
end
function movemarkers!(r,side)
    while isborder(r,side)==0
            move!(r,side)
            putmarker!(r)
    end
end
function moveback!(r,side)
    while ismarker(r)==1
    move!(r,HorizonSide((Int(side)+2)%4))
    end
end
krest!(r)


text:
# --------------------------------------------------
# mask functions

function apply_mask!(datacube::AbstractSpectralCube, maskcube::AbstractSpectralCube)
    datacube.data[isnan.(maskcube.data)] .= NaN;
    datacube.data[iszero.(maskcube.data)] .= NaN;
    return datacube
end

function apply_mask(datacube::AbstractSpectralCube, maskcube::AbstractSpectralCube)
    datacube2 = deepcopy(datacube)
    datacube2.data[isnan.(maskcube.data)] .= NaN;
    datacube2.data[iszero.(maskcube.data)] .= NaN;
    return datacube2
end

function apply_mask!(datacube::AbstractSpectralCube, maskmap::AbstractMap)
    mask2d = maskmap.data
    s = size(maskmap.data)
    mask3d = zeros(eltype(maskmap.data), (s[1], s[2], size(datacube.data, 3)))
    mask3d .= mask2d
    datacube.data[isnan.(mask3d)] .= NaN;
    datacube.data[iszero.(mask3d)] .= NaN;
    return datacube
end

function apply_mask!(datacube::AbstractSpectralCube, maskarr::Union{Array{T, 2} where T <: Real, BitArray})
    mask2d = maskarr
    s = size(mask2d)
    T = eltype(maskarr)
    mask3d = zeros(T, (s[1], s[2], size(datacube.data, 3)))
    mask3d .= mask2d
    datacube.data[isnan.(mask3d)] .= NaN;
    datacube.data[iszero.(mask3d)] .= NaN;
    return datacube
end


text:
module C3D

using VaxData

@enum Endian LE=1 BE=2

export readc3d,
       writetrc

export C3DFile

include("parameters.jl")
include("groups.jl")
include("header.jl")
include("validate.jl")

struct C3DFile
    name::String
    header::Header
    groups::Dict{Symbol, Group}
    point::Dict{String, Array{Union{Missing, Float32},2}}
    residual::Dict{String, Array{Union{Missing, Float32},1}}
    analog::Dict{String, Array{Float32,1}}
end

include("util.jl")

function C3DFile(name::String, header::Header, groups::Dict{Symbol,Group},
                 point::AbstractArray, residuals::AbstractArray, analog::AbstractArray;
                 missingpoints::Bool=true, strip_prefixes::Bool=false)
    fpoint = Dict{String,Array{Union{Missing, Float32},2}}()
    fresiduals = Dict{String,Array{Union{Missing, Float32},1}}()
    fanalog = Dict{String,Array{Float32,1}}()

    l = size(point, 1)
    allpoints = 1:l

    if strip_prefixes
        if haskey(groups, :SUBJECTS) && groups[:SUBJECTS][Int, :USES_PREFIXES] == 1
            rgx = Regex("("*join(groups[:SUBJECTS][Vector{String}, :LABEL_PREFIXES], '|')*
                        ")(?<label>\\w*)")
        else
            rgx = r":(?<label>\w*)"
        end
    end

    numpts = groups[:POINT][Int, :USED]
    if !iszero(numpts)
        for (idx, ptname) in enumerate(groups[:POINT][Vector{String}, :LABELS][1:numpts])
            if strip_prefixes
                m = match(rgx, ptname)
                if !isnothing(m) && !isnothing(m[:label])
                    ptname = m[:label]
                end
            end

            fpoint[ptname] = point[:,((idx-1)*3+1):((idx-1)*3+3)]
            fresiduals[ptname] = residuals[:,idx]
            if missingpoints
                invalidpoints = findall(x -> x === -1.0f0, fresiduals[ptname])
                calculatedpoints = findall(iszero, fresiduals[ptname])
                goodpoints = setdiff(allpoints, invalidpoints ∪ calculatedpoints)
                fpoint[ptname][invalidpoints, :] .= missing
                fresiduals[ptname][goodpoints] = calcresiduals(fresiduals[ptname], abs(groups[:POINT][Float32, :SCALE]))[goodpoints]
                fresiduals[ptname][invalidpoints] .= missing
                fresiduals[ptname][calculatedpoints] .= 0.0f0
            end
        end
    end

    if !iszero(groups[:ANALOG][Int, :USED])
        for (idx, name) in enumerate(groups[:ANALOG][Vector{String}, :LABELS][1:groups[:ANALOG][Int, :USED]])
            fanalog[name] = analog[:, idx]
        end
    end

    return C3DFile(name, header, groups, fpoint, fresiduals, fanalog)
end

function Base.show(io::IO, f::C3DFile)
    if get(io, :compact, true)
        print(io, "C3DFile(\"", f.name, "\")")
    else
        length = (f.groups[:POINT][Int, :FRAMES] == typemax(UInt16)) ?
            f.groups[:POINT][Int, :LONG_FRAMES] :
            f.groups[:POINT][Int, :FRAMES]

        print(io, "C3DFile(\"", f.name, "\", ",
              length, "sec, ",
              f.groups[:POINT][Int, :USED], " points, ",
              f.groups[:ANALOG][Int, :USED], " analog channels)")
    end
end

function calcresiduals(x::AbstractVector, scale)
    residuals = (reinterpret.(UInt32, x) .>> 16) .& 0xff .* scale
end

function readdata(io::IOStream, groups::Dict{Symbol,Group}, FEND::Endian, FType::Type{T}) where T <: Union{Float32,VaxFloatF}
    seek(io, (groups[:POINT][Int, :DATA_START]-1)*512)

    format = groups[:POINT][Float32, :SCALE] > 0 ? Int16 : FType

    # Read data in a transposed structure for better read/write speeds due to Julia being
    # column-order arrays
    numframes::Int = groups[:POINT][Int, :FRAMES]
    if numframes == typemax(UInt16) && haskey(groups, :TRIAL) && haskey(groups[:TRIAL][:params], :ACTUAL_END_FIELD)
        numframes = only(reinterpret(Int32, groups[:TRIAL][Vector{Int16}, :ACTUAL_END_FIELD]))
    end
    nummarkers = groups[:POINT][Int, :USED]
    hasmarkers = !iszero(nummarkers)
    if hasmarkers
        point = Array{Float32,2}(undef, nummarkers*3, numframes)
        residuals = Array{Float32,2}(undef, nummarkers, numframes)

        nb = nummarkers*4
        pointidxs = filter(x -> x % 4 != 0, 1:nb)
        residxs = filter(x -> x % 4 == 0, 1:nb)

        pointtmp = Vector{format}(undef, nb)
        pointview = view(pointtmp, pointidxs)
        resview = view(pointtmp, residxs)
    else
        point = Array{Float32,2}(undef, 0,0)
        residuals = Array{Float32,2}(undef, 0,0)
    end

    numchannels = groups[:ANALOG][Int, :USED]
    haschannels = !iszero(numchannels)
    if haschannels
        # Analog Samples Per Frame => ASPF
        aspf = convert(Int, groups[:ANALOG][Float32, :RATE]/groups[:POINT][Float32, :RATE])
        analog = Array{Float32,2}(undef, numchannels, aspf*numframes)

        analogtmp = Matrix{format}(undef, (numchannels,aspf))
    else
        analog = Array{Float32,2}(undef, 0,0)
    end

    @inbounds for i in 1:numframes
        if hasmarkers
            saferead!(io, pointtmp, FEND)
            point[:,i] = pointview
            residuals[:,i] = resview
        end
        if haschannels
            saferead!(io, analogtmp, FEND)
            analog[:,((i-1)*aspf+1):(i*aspf)] = analogtmp
        end
    end

    if hasmarkers && format == Int16
        # Multiply or divide by [:point][:scale]
        POINT_SCALE = groups[:POINT][Float32, :SCALE]
        point .*= abs(POINT_SCALE)
    end

    if haschannels
        if numchannels == 1
            ANALOG_OFFSET = groups[:ANALOG][Float32, :OFFSET]
            SCALE = (groups[:ANALOG][Float32, :GEN_SCALE] * groups[:ANALOG][Float32, :SCALE])
            analog .= (analog .- ANALOG_OFFSET) .* SCALE
        else
            VECANALOG_OFFSET = groups[:ANALOG][Vector{Int}, :OFFSET][1:numchannels]
            VECSCALE = groups[:ANALOG][Float32, :GEN_SCALE] .*
                            groups[:ANALOG][Vector{Float32}, :SCALE][1:numchannels]

            analog .= (analog .- VECANALOG_OFFSET) .* VECSCALE
        end
    end

    return (permutedims(point), permutedims(residuals), permutedims(analog))
end

function saferead(io::IOStream, ::Type{T}, FEND::Endian) where T
    if FEND == LE
        return ltoh(read(io, T))
    else
        return ntoh(read(io, T))
    end
end

function saferead!(io::IOStream, x::AbstractArray, FEND::Endian)
    if FEND == LE
        x .= ltoh.(read!(io, x))
    else
        x .= ntoh.(read!(io, x))
    end
    nothing
end

function saferead(io::IOStream, ::Type{T}, FEND::Endian, dims)::Array{T} where T
    if FEND == LE
        return ltoh.(read!(io, Array{T}(undef, dims)))
    else
        return ntoh.(read!(io, Array{T}(undef, dims)))
    end
end

function saferead(io::IOStream, ::Type{VaxFloatF}, FEND::Endian)::Float32
    if FEND == LE
        return convert(Float32, ltoh(read(io, VaxFloatF)))
    else
        return convert(Float32, ntoh(read(io, VaxFloatF)))
    end
end

function saferead(io::IOStream, ::Type{VaxFloatF}, FEND::Endian, dims::NTuple{N,Int})::Array{Float32,N} where N
    if FEND == LE
        return convert(Array{Float32,N}, ltoh.(read!(io, Array{VaxFloatF}(undef, dims))))
    else
        return convert(Array{Float32,N}, ntoh.(read!(io, Array{VaxFloatF}(undef, dims))))
    end
end

"""
    readc3d(fn)

Read the C3D file at `fn`.

# Keyword arguments
- `paramsonly::Bool = false`: Only reads the header and parameters
- `validateparams::Bool = true`: Validates parameters against C3D requirements
- `missingpoints::Bool = true`: Sets invalid points to `missing`
"""
function readc3d(fn::AbstractString; paramsonly=false, validate=true,
                 missingpoints=true, strip_prefixes=false)
    if !isfile(fn)
        error("File ", fn, " cannot be found")
    end

    io = open(fn, "r")

    groups, header, FEND, FType = _readparams(fn, io)

    if validate
        validatec3d(header, groups)
    end

    if paramsonly
        point = Dict{String,Array{Union{Missing, Float32},2}}()
        residual = Dict{String,Array{Union{Missing, Float32},1}}()
        analog = Dict{String,Array{Float32,1}}()
    else
        (point, residual, analog) = readdata(io, groups, FEND, FType)
    end

    close(io)

    res = C3DFile(fn, header, groups, point, residual, analog;
                  missingpoints, strip_prefixes)

    return res
end

function _readparams(fn::String, io::IOStream)
    params_ptr = read(io, UInt8)

    if read(io, UInt8) != 0x50
        error("File ", fn, " is not a valid C3D file")
    end

    # Jump to parameters block
    seek(io, (params_ptr - 1) * 512)

    # Skip 2 reserved bytes
    # TODO: store bytes for saving modified files
    skip(io, 2)

    paramblocks = read(io, UInt8)
    proctype = read(io, Int8) - 83

    FType = Float32

    # Deal with host big-endianness in the future
    if proctype == 1
        # little-endian
        FEND = LE
    elseif proctype == 2
        # DEC floats; little-endian
        FType = VaxFloatF
        FEND = LE
    elseif proctype == 3
        # big-endian
        FEND = BE
    else
        error("Malformed processor type. Expected 1, 2, or 3. Found ", proctype)
    end

    mark(io)
    header = readheader(io, FEND, FType)
    reset(io)
    unmark(io)

    gs = Array{Group,1}()
    ps = Array{Parameter,1}()
    moreparams = true
    fail = 0
    np = 0

    skip(io, 1)
    if read(io, Int8) < 0
        # Group
        skip(io, -2)
        push!(gs, readgroup(io, FEND, FType))
        np = gs[end].pos + gs[end].np + abs(gs[end].nl) + 2
        moreparams = gs[end].np != 0 ? true : false
    else
        # Parameter
        skip(io, -2)
        push!(ps, readparam(io, FEND, FType))
        np = ps[end].pos + ps[end].np + abs(ps[end].nl) + 2
        moreparams = ps[end].np != 0 ? true : false
    end

    while moreparams
        # Mark current position in file in case the pointer is incorrect
        mark(io)
        if fail === 0 && np != position(io)
                # @debug "Pointer mismatch at position $(position(io)) where pointer was $np"
                seek(io, np)
        elseif fail > 1 # this is the second failed attempt
            # @debug "Second failed parameter read attempt from $(position(io))"
            break
        end

        # Read the next two bytes to get the gid
        skip(io, 1)
        local gid = read(io, Int8)
        if gid < 0 # Group
            # Reset to the beginning of the group
            skip(io, -2)
            try
                push!(gs, readgroup(io, FEND, FType))
                np = gs[end].pos + gs[end].np + abs(gs[end].nl) + 2
                moreparams = gs[end].np != 0 ? true : false # break if the pointer is 0 (ie the parameters are finished)
                fail = 0 # reset fail counter following a successful read
            catch e
                # Last readgroup failed, possibly due to a bad pointer. Reset to the ending
                # location of the last successfully read parameter and try again. Count the failure.
                reset(io)
                # @debug "Read group failed, last parameter ended at $(position(io)), pointer at $np" fail
                fail += 1
            finally
                unmark(io) # Unmark the file regardless
            end
        elseif gid > 0 # Parameter
            # Reset to the beginning of the parameter
            skip(io, -2)
            try
                push!(ps, readparam(io, FEND, FType))
                np = ps[end].pos + ps[end].np + abs(ps[end].nl) + 2
                moreparams = ps[end].np != 0 ? true : false
                fail = 0
            catch e
                reset(io)
                # @debug "Read group failed, last parameter ended at $(position(io)), pointer at $np" fail
                fail += 1
            finally
                unmark(io)
            end
        else # Last parameter pointer is incorrect (assumption)
            # The group ID should never be zero, if it is, the most likely explanation is
            # that the pointer is incorrect (eg the pointer was not fixed when the previously
            # last parameter was deleted or moved)
            # @debug "Bad last position. Assuming parameter section is finished."
            break
        end
    end

    groups = Dict{Symbol,Group}()
    gids = Dict{Int8,Symbol}()

    for group in gs
        groups[group.symname] = group
        gids[abs(group.gid)] = group.symname
    end

    for param in ps
        groups[gids[param.gid]].params[param.symname] = param
    end

    return (groups, header, FEND, FType)
end

end # module


text:
###############################################################################
#
#   Error objects
#
###############################################################################

mutable struct ErrorConstrDimMismatch <: Exception
  expect_r::Int
  expect_c::Int
  get_r::Int
  get_c::Int
  get_l::Int

  function ErrorConstrDimMismatch(er::Int, ec::Int, gr::Int, gc::Int)
    e = new(er, ec, gr, gc, -1)
    return e
  end

  function ErrorConstrDimMismatch(er::Int, ec::Int, gl::Int)
    e = new(er, ec, -1, -1, gl)
    return e
  end

  function ErrorConstrDimMismatch(er::Int, ec::Int, a::Array{T, 2}) where {T}
    gr, gc = size(a)
    return ErrorConstrDimMismatch(er, ec, gr, gc)
  end

  function ErrorConstrDimMismatch(er::Int, ec::Int, a::Array{T, 1}) where {T}
    gl = length(a)
    return ErrorConstrDimMismatch(er, ec, gl)
  end
end

function Base.showerror(io::IO, e::ErrorConstrDimMismatch)
  if e.get_l == -1
    print(io, "Expected dimension $(e.expect_r) x $(e.expect_c), ")
    print(io, "got $(e.get_r) x $(e.get_c)")
  else
    print(io, "Expected an array of length $(e.expect_r * e.expect_c), ")
    print(io, "got $(e.get_l)")
  end
end

const error_dim_negative = ErrorException("Dimensions must be non-negative")


text:
function get_cirf(obj::Frames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getCIRF", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_ecliptic(obj::Frames, arg0::IERSConventions)
    return jcall(obj, "getEcliptic", Frame, (IERSConventions,), arg0)
end

function get_eme2000(obj::Frames)
    return jcall(obj, "getEME2000", FactoryManagedFrame, ())
end

function get_frame(obj::Frames, arg0::Predefined)
    return jcall(obj, "getFrame", Frame, (Predefined,), arg0)
end

function get_gcrf(obj::Frames)
    return jcall(obj, "getGCRF", Frame, ())
end

function get_gtod(obj::Frames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getGTOD", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_gtod(obj::Frames, arg0::jboolean)
    return jcall(obj, "getGTOD", FactoryManagedFrame, (jboolean,), arg0)
end

function get_icrf(obj::Frames)
    return jcall(obj, "getICRF", Frame, ())
end

function get_itrf(obj::Frames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getITRF", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_itrf(obj::Frames, arg0::ITRFVersion, arg1::IERSConventions, arg2::jboolean)
    return jcall(obj, "getITRF", VersionedITRF, (ITRFVersion, IERSConventions, jboolean), arg0, arg1, arg2)
end

function get_itrf_equinox(obj::Frames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getITRFEquinox", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_mod(obj::Frames, arg0::IERSConventions)
    return jcall(obj, "getMOD", FactoryManagedFrame, (IERSConventions,), arg0)
end

function get_mod(obj::Frames, arg0::jboolean)
    return jcall(obj, "getMOD", FactoryManagedFrame, (jboolean,), arg0)
end

function get_pz9011(obj::Frames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getPZ9011", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_teme(obj::Frames)
    return jcall(obj, "getTEME", FactoryManagedFrame, ())
end

function get_tirf(obj::Frames, arg0::IERSConventions)
    return jcall(obj, "getTIRF", FactoryManagedFrame, (IERSConventions,), arg0)
end

function get_tirf(obj::Frames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getTIRF", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_tod(obj::Frames, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getTOD", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_tod(obj::Frames, arg0::jboolean)
    return jcall(obj, "getTOD", FactoryManagedFrame, (jboolean,), arg0)
end

function get_veis1950(obj::Frames)
    return jcall(obj, "getVeis1950", FactoryManagedFrame, ())
end

function of(::Type{Frames}, arg0::TimeScales, arg1::CelestialBodies)
    return jcall(Frames, "of", Frames, (TimeScales, CelestialBodies), arg0, arg1)
end

function of(::Type{Frames}, arg0::TimeScales, arg1::Supplier)
    return jcall(Frames, "of", Frames, (TimeScales, Supplier), arg0, arg1)
end


text:
module Bezier

using StaticArrays
using LinearAlgebra

struct BezierCurve{N,D,T}
	p::SVector{N,SVector{D,T}}
	@inline BezierCurve(p::Vararg{Any,N}) where{N} = BezierCurve{N}(p...)
	@inline BezierCurve{N}(p::StaticVector{D}...) where{N,D} =
		new{N,D,promote_type(eltype.(p)...)}(p)
	@inline BezierCurve{N}(p::AbstractVector...) where{N} =
		BezierCurve{N}(SVector{length(first(p))}.(p)...)
end

@inline degree(::Type{<:BezierCurve{N}}) where{N} = N
@inline dim(::Type{<:BezierCurve{N,D}}) where{N,D} = D
@inline pttype(::Type{<:BezierCurve{N,D,T}}) where{N,D,T} = SVector{D,T}
@inline degree(b::B) where{B<:BezierCurve} = degree(B)
@inline dim(b::B) where{B<:BezierCurve} = dim(B)
@inline pttype(b::B) where{B<:BezierCurve} = pttype(B)
@inline controlpoint(b::BezierCurve, i::Integer) = b.p[i]

"""
    interpolate(b::BezierCurve, t::Real)

Computes the position of the point at parameter `t` on the curve `b`.
"""
@inline interpolate(b::BezierCurve, t::Real) = _interpolate(b.p, t)

# recursive definition for de Casteljau's algorithm:
@inline _interpolate(p::SVector{1}, _) = p[1]
@inline _interpolate(p::SVector{N,P}, t) where{N,P} =
	_interpolate(SVector{N-1,Base.return_types(*,(P,typeof(t)))[1]}(
		t*p[i+1]+(1-t)*p[i] for i in 1:N-1), t)

"""
    interpolate(b::BezierCurve, v::AbstractVector{<:Real})

Interpolates `b` at all values t ∈ `v`.
"""
@inline interpolate(b::BezierCurve, v::AbstractVector{<:Real}) =
	(interpolate(b, t) for t in v)

"""
    interpolate(b::BezierCurve, n::Integer)

Interpolates as `n` regularly-spaced points (in parameter space).
"""
@inline interpolate(b::BezierCurve, n::Integer) = interpolate(b, (0:n)/n)

end


text:
function firstnonconsecutive(arr::Array)
    # bisection search, more code but less runtime
    if sizeof(arr) < 2
        return nothing
    end
    a = firstindex(arr)
    b = lastindex(arr)
    offset = arr[a]-a
    if arr[b] == b+offset
        return nothing
    end
    while a < b
        c = (a+b)÷2
        if arr[c] == c+offset
            a = c+1
        else
            b = c
        end
    end
    arr[b]
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using Test
using ExcelReaders
using Mimi
using DataFrames
using CSV

include("../src/rice2010.jl")
using .Rice2010

m = getrice()
run(m)

parameter_filename = joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm")

f=openxl(parameter_filename)
regions = ["US", "EU", "Japan", "Russia", "Eurasia", "China", "India", "MidEast", "Africa", "LatAm", "OHI", "OthAsia"]

#Function to get true values form Rice2010 Excel
function Truth(range::AbstractString)
	true_vals=Array{Float64}(undef, 60, length(regions))
    for (i,r) = enumerate(regions)
        data=readxl(f,"$(r)!$(range)")
        for n=1:60
			true_vals[n,i] = data[n]
        end
    end
    return true_vals
end

# Test Precision
Precision = 1.0e-11

@testset "mimi-rice-2010" begin

#------------------------------------------------------------------------------
#   1. Run tests on the whole model
#------------------------------------------------------------------------------

@testset "mimi-rice-2010-model" begin

# TATM Test (temperature increase)
True_TATM = vec(readxl(f,"Global!B70:BI70"))
@test maximum(abs, m[:climatedynamics, :TATM] .- True_TATM) ≈ 0. atol=Precision

# MAT Test (carbon concentration atmosphere)
True_MAT = vec(readxl(f,"Global!B59:BI59"))
@test maximum(abs, m[:co2cycle, :MAT] .- True_MAT) ≈ 0. atol=Precision

# DAMFRAC Test (damages fraction)
True_DAMFRAC = Truth("B63:BI63")
@test maximum(abs, m[:damages, :DAMFRAC] .- True_DAMFRAC) ≈ 0. atol=Precision

# DAMAGES Test (damages $)
True_DAMAGES = Truth("B64:BI64")
@test maximum(abs, m[:damages, :DAMAGES] .- True_DAMAGES) ≈ 0. atol=Precision

# E Test (emissions)
True_E = vec(readxl(f,"Global!B55:BI55"))
@test maximum(abs, m[:emissions, :E] .- True_E) ≈ 0. atol=Precision

# YGROSS Test (gross output)
True_YGROSS = Truth("B61:BI61")
@test maximum(abs, m[:grosseconomy, :YGROSS] .- True_YGROSS) ≈ 0. atol=Precision

# CPC Test (per capita consumption)
True_CPC = Truth("B88:BI88")
@test maximum(abs, m[:neteconomy, :CPC] .- True_CPC) ≈ 0. atol=Precision

# FORC Test (radiative forcing)
True_FORC = vec(readxl(f,"Global!B71:BI71"))
@test maximum(abs, m[:radiativeforcing, :FORC] .- True_FORC) ≈ 0. atol=Precision

# TOTALSLR Test (total sea level rise)
True_TOTALSLR = vec(readxl(f,"SLR!B62:BI62"))
@test maximum(abs, m[:sealevelrise, :TOTALSLR] .- True_TOTALSLR) ≈ 0. atol=Precision

# SLRDAMAGES Test (damages from sea level rise)
True_SLRDAMAGES = Truth("B50:BI50")
@test maximum(abs, m[:sealeveldamages, :SLRDAMAGES] .- True_SLRDAMAGES) ≈ 0. atol=Precision

True_PERIODUTILITY = Truth("B94:BI94")
@test maximum(abs, m[:welfare, :CEMUTOTPER] .- True_PERIODUTILITY) ≈ 0. atol=Precision

True_UTILITY = readxl(f,"Global!B77:B77")
@test maximum(abs, m[:welfare, :UTILITY] .- True_UTILITY) ≈ 0. atol=Precision

end #mimi-rice-2010-model testset

#------------------------------------------------------------------------------
#   2. Run tests to make sure integration version (Mimi v0.5.0)
#   values match Mimi 0.4.0 values
#------------------------------------------------------------------------------

@testset "mimi-rice-2010-integration" begin

nullvalue = -999.999

for c in map(name, Mimi.compdefs(m)), v in Mimi.variable_names(m, c)

    #load data for comparison
    filepath = joinpath(@__DIR__, "..", "data", "validation_data_v040", "$c-$v.csv")
    results = m[c, v]

    if typeof(results) <: Number
        validation_results = CSV.read(filepath)[1,1]

    else
        validation_results = convert(Array, CSV.read(filepath))

        #remove NaNs
        results[ismissing.(results)] .= nullvalue
        results[isnan.(results)] .= nullvalue
        validation_results[isnan.(validation_results)] .= nullvalue

        #match dimensions
        if size(validation_results,1) == 1
            validation_results = validation_results'
        end
    end

    @test results ≈ validation_results atol = Precision

end #for loop

end #mimi-rice-2010-integration testset

end #mimi-rice-2010 testset

nothing


text:
using Weave, Compat
using Base.Test

function pljtest(source, resfile, doctype)
  weave("documents/$source", out_path = "documents/plotsjl/$resfile", doctype=doctype, plotlib=nothing)
  result = @compat readstring(open("documents/plotsjl/$resfile"))
  ref = @compat readstring(open("documents/plotsjl/$resfile.ref"))
  @test result == ref
end

pljtest("plotsjl_test.jmd", "plotsjl_test.md", "pandoc")
pljtest("plotsjl_test.jmd", "plotsjl_test.tex", "tex")


text:
include(joinpath(@__DIR__, "../PlotUtils/PlotUtils.jl"))
include(joinpath(@__DIR__, "../PlotUtils/imports.jl"))

using Distributed

rmprocs(filter(w -> w > 1, workers()))
addprocs(15)

@everywhere include(joinpath(@__DIR__, "../PlotUtils/PlotUtils.jl"))
@everywhere include(joinpath(@__DIR__, "../PlotUtils/imports.jl"))

if length(ARGS) == 0
  results_dir = "/scratchdata/alui2/cytof/results/repfam/test-sim-6-5"
else
  results_dir = ARGS[1]  # path to results directory
end

### Create metrics ###

# Name of output file
@everywhere OUTPUT_FILE = "output.bson"

println("Producing metrics ...")

# directories to loop through
results_dirs = ["$(results_dir)/seed_$(seed)/scale_$(scale)"
                for seed in 1:3
                for scale in (0, 1, 10)]

# Number of directories
num_dirs = length(results_dirs)

parsekey(T, key, path) = parse(T, match(Regex("(?<=$(key)_)\\d+"), path).match)

# Construct a data frame for R
function fill_R_df!(path, metrics, R_df)
  # Parse seed
  seed = parsekey(Int, "seed", path)

  # Parse scale
  scale = parsekey(Int, "scale", path)

  for metric in metrics
    Kmcmc, m = metric
    I = length(m[:Rmean])
    for i in 1:I
      append!(R_df, DataFrames.DataFrame(Kmcmc=Kmcmc, i=i,
                                         seed=seed, scale=scale,
                                         Rmean=m[:Rmean][i],
                                         R_975=m[:R_975][i],
                                         R_025=m[:R_025][i]))
    end
  end
end

# Create make metrics function which only takes directory as input
@everywhere make_metrics(d) = PlotUtils.make_metrics(d, OUTPUT_FILE, thresh=.01)

# Get all metrics in parallel
allmetrics = pmap(make_metrics, results_dirs)

# Put results in a data frame.
R_df = let
  _R_df = DataFrames.DataFrame()
  for (d, metrics) in zip(results_dirs, allmetrics)
    fill_R_df!(d, metrics, _R_df)
  end
  _R_df
end

# TODO: Plot R-metrics nicely
# for each seed, make the R by Kmcmc graphs
for R_df_seed in groupby(R_df, :seed)
  seed = unique(R_df_seed.seed)[1]
  unique_i = unique(R_df_seed.i)
  I = length(unique_i)
  plt.figure()
  for i in 1:I
    R_df_seed_i = R_df_seed[R_df_seed.i .== i, :]
    plt.subplot(I, 1, i)
    # Put graphs for different scales in same graph for same seed
    for R_df_seed_i_scale in groupby(R_df_seed_i, :scale)
      scale = unique(R_df_seed_i_scale.scale)[1]
      df = sort(R_df_seed_i_scale, :Kmcmc)
      plt.plot(df.Kmcmc, df.Rmean, label="scale: $(scale)", lw=2, marker="o")
      plt.fill_between(df.Kmcmc, df.R_025, df.R_975, alpha=.5)
    end
    plt.ylabel("R$(i)")
  end

  outdir = "$(results_dir)/metrics/seed_$(seed)"
  mkpath(outdir)

  # plt.legend(loc="lower right")
  plt.legend(loc="best")
  plt.xlabel("K")
  plt.savefig("$(outdir)/R.pdf", bbox_inches="tight")
  plt.close()
end

# Remove extra processors
rmprocs(filter(w -> w > 1, workers()))

println("DONE!")


text:
# File: fSkew.jl
# Author:
# Date: 29-June-2020

# SKEW
function m3 = fSkew(x,flag,f)
# --------------------------------------------------------------------------
# x:     returns
# flag:  0 = sample, 1 = population
# f:     reporting frequency
# m3:    skew (sample or population)
# --------------------------------------------------------------------------

    m3 = skewness(x,flag)./sqrt(f);
end


text:
using FDFD, GeometryPrimitives
using Test
using LinearAlgebra

ENV["PARDISOLICMESSAGE"]=0

function compare_solvers_dipole()
    grd = Grid(0.01, [15 15], [-3 3], [-3 3]);
    dev = Device(grd, 2*π*200e12);
    setup_src!(dev, Point(0, 0));

    ENV["FDFD_SOLVER"] = "julia"
    field_julia = solve(dev);
    ENV["FDFD_SOLVER"] = "pardiso"
    field_pardiso = solve(dev);
    ENV["FDFD_SOLVER"] = "julia"

    equality_real = real.(field_julia) ≈ real.(field_pardiso)
    equality_imag = imag.(field_julia) ≈ imag.(field_pardiso)
    return equality_real && equality_imag
end

function compare_solvers_wg()
    grd = Grid(0.02, [15, 15], [0.0, 10.0], [-1.0, 1.0]);
    dev = Device(grd, 2*π*200e12);
    setup_ϵᵣ!(dev, [Box([5.0, 0.0, 0.0], [Inf, 0.3, Inf], Matrix{Float64}(I, 3, 3), 12)])
    add_mode!(dev, Mode(TM, x̂, 3.5, Point(1.0, 0), 0.8));

    ENV["FDFD_SOLVER"] = "julia"
    field_julia = solve(dev);
    ENV["FDFD_SOLVER"] = "pardiso"
    field_pardiso = solve(dev);
    ENV["FDFD_SOLVER"] = "julia"

    equality_real = real.(field_julia) ≈ real.(field_pardiso)
    equality_imag = imag.(field_julia) ≈ imag.(field_pardiso)
    return equality_real && equality_imag
end

function compare_solvers_modwg()
    nsidebands = 1
    L = 5.0
    grd = Grid(0.01, [15 15], [0.0 L], [-1.0 1.0]);
    dev = ModulatedDevice(grd, 2π*1.939e14, 4.541e14, nsidebands)
    a = 0.2202
    q = 2.9263
    region_wg(x,y) = -a/2<=y<=a/2
    region_modulation(x,y) = 1<=x<=(L-1) && -a/2<=y<=0
    modulation_function(x,y) = exp.(1im*q*x)
    setup_ϵᵣ!(dev, region_wg, 12.25)
    setup_Δϵᵣ!(dev, region_modulation, modulation_function)
    add_mode!(dev, Mode(TM, x̂, 3.5, [0.2, 0], 4*a))

    ENV["FDFD_SOLVER"] = "julia"
    field_julia = solve(dev);
    ENV["FDFD_SOLVER"] = "pardiso"
    field_pardiso = solve(dev);
    ENV["FDFD_SOLVER"] = "julia"

    equality_real = all([real.(field_julia[i]) ≈ real.(field_pardiso[i]) for i = 1:2*nsidebands+1])
    equality_imag = all([imag.(field_julia[i]) ≈ imag.(field_pardiso[i]) for i = 1:2*nsidebands+1])
    return equality_real && equality_imag
end

@testset "FDFD" begin
    @testset "Solver" begin
        @test compare_solvers_dipole()
        @test compare_solvers_wg()
        @test compare_solvers_modwg()
    end
end


text:
import gmsh

gmsh.initialize(append!(["gmsh"], ARGS))

gmsh.model.add("t2")

lc = 1e-2
gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(.1, 0,  0, lc, 2)
gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
gmsh.model.geo.addPoint(0, .3, 0, lc, 4)
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(3, 2, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(0, [1, 2], 1)
gmsh.model.addPhysicalGroup(1, [1, 2], 2)
gmsh.model.addPhysicalGroup(2, [1], 6)
gmsh.model.setPhysicalName(2, 6, "My surface")

gmsh.model.geo.addPoint(0, .4, 0, lc, 5)
gmsh.model.geo.addLine(4, 5, 5)

gmsh.model.geo.translate([(0, 5)], -0.02, 0, 0)

gmsh.model.geo.rotate([(0, 5)], 0,0.3,0, 0,0,1, -pi/4)

ov = gmsh.model.geo.copy([(0, 3)])
gmsh.model.geo.translate(ov, 0, 0.05, 0)

gmsh.model.geo.addLine(3, ov[1][2], 7)
gmsh.model.geo.addLine(ov[1][2], 5, 8)
gmsh.model.geo.addCurveLoop([5,-8,-7,3], 10)
gmsh.model.geo.addPlaneSurface([10], 11)

ov = gmsh.model.geo.copy([(2, 1), (2, 11)])
gmsh.model.geo.translate(ov, 0.12, 0, 0)

println("New surfaces ", ov[1][2], " and ", ov[2][2])

gmsh.model.geo.addPoint(0., 0.3, 0.12, lc, 100)
gmsh.model.geo.addPoint(0.1, 0.3, 0.12, lc, 101)
gmsh.model.geo.addPoint(0.1, 0.35, 0.12, lc, 102)

gmsh.model.geo.synchronize()
xyz = gmsh.model.getValue(0, 5, [])
gmsh.model.geo.addPoint(xyz[1], xyz[2], 0.12, lc, 103)

gmsh.model.geo.addLine(4, 100, 110)
gmsh.model.geo.addLine(3, 101, 111)
gmsh.model.geo.addLine(6, 102, 112)
gmsh.model.geo.addLine(5, 103, 113)
gmsh.model.geo.addLine(103, 100, 114)
gmsh.model.geo.addLine(100, 101, 115)
gmsh.model.geo.addLine(101, 102, 116)
gmsh.model.geo.addLine(102, 103, 117)

gmsh.model.geo.addCurveLoop([115, -111, 3, 110], 118)
gmsh.model.geo.addPlaneSurface([118], 119)
gmsh.model.geo.addCurveLoop([111, 116, -112, -7], 120)
gmsh.model.geo.addPlaneSurface([120], 121)
gmsh.model.geo.addCurveLoop([112, 117, -113, -8], 122)
gmsh.model.geo.addPlaneSurface([122], 123)
gmsh.model.geo.addCurveLoop([114, -110, 5, 113], 124)
gmsh.model.geo.addPlaneSurface([124], 125)
gmsh.model.geo.addCurveLoop([115, 116, 117, 114], 126)
gmsh.model.geo.addPlaneSurface([126], 127)

gmsh.model.geo.addSurfaceLoop([127, 119, 121, 123, 125, 11], 128)
gmsh.model.geo.addVolume([128], 129)

ov2 = gmsh.model.geo.extrude([ov[2]], 0, 0, 0.12)

gmsh.model.geo.mesh.setSize([(0,103), (0,105), (0,109), (0,102), (0,28),
                             (0, 24), (0,6), (0,5)], lc * 3)

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(3, [129,130], 1)
gmsh.model.setPhysicalName(3, 1, "The volume")

gmsh.model.mesh.generate(3)

gmsh.write("t2.msh")

gmsh.finalize()

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# inverse
ilab = 1:10
olab = [randstring(1) for i =1:10]
fst = WFST(ilab,olab)
add_arc!(fst, 1, 1, ilab[1] , olab[1] , rand())
add_arc!(fst, 1, 2, ilab[2] , olab[2] , rand())
add_arc!(fst, 2, 3, ilab[3] , olab[3] , rand())
add_arc!(fst, 2, 2, ilab[4] , olab[4] , rand())
add_arc!(fst, 3, 3, ilab[5] , olab[5] , rand())
add_arc!(fst, 3, 4, ilab[6] , olab[6] , rand())
add_arc!(fst, 4, 4, ilab[7] , olab[7] , rand())
add_arc!(fst, 4, 5, ilab[8] , olab[8] , rand())
add_arc!(fst, 5, 5, ilab[9] , olab[9] , rand())
add_arc!(fst, 5, 6, ilab[10], olab[10], rand())
initial!(fst,1)
final!(fst,6,rand())
ifst = inv(fst)
@test all([a in ilab for a in get_oalphabet(ifst)])
@test all([a in olab for a in get_ialphabet(ifst)])
@test isfinal(ifst,6)
@test isinitial(ifst,1)
@test get_weight(fst,6) == get_weight(ifst,6)


text:
"""
`CSV.readline(io::IO, q='"', e='\\', buf::IOBuffer=IOBuffer())` => `String`
`CSV.readline(source::CSV.Source)` => `String`

read a single line from `io` (any `IO` type) or a `CSV.Source` as a string, accounting for potentially embedded newlines in quoted fields (e.g. value1, value2, \"value3 with \n embedded newlines\"). Can optionally provide a `buf::IOBuffer` type for buffer reuse
"""
function readline end

function readline(io::IO,q::UInt8,e::UInt8,buf::IOBuffer=IOBuffer())
    while !eof(io)
        b = unsafe_read(io, UInt8)
        Base.write(buf, b)
        if b == q
            while !eof(io)
                b = unsafe_read(io, UInt8)
                Base.write(buf, b)
                if b == e
                    b = unsafe_read(io, UInt8)
                    Base.write(buf, b)
                elseif b == q
                    break
                end
            end
        elseif b == NEWLINE
            break
        elseif b == RETURN
            !eof(io) && unsafe_peek(io) == NEWLINE && Base.write(buf,unsafe_read(io, UInt8))
            break
        end
    end
    return takebuf_string(buf)
end
readline(io::IO, q='"', e='\\', buf::IOBuffer=IOBuffer()) = readline(io, UInt8(q), UInt8(e), buf)
readline(source::CSV.Source) = readline(source.data, source.options.quotechar, source.options.escapechar)

"""
`CSV.readsplitline(io, d=',', q='"', e='\\', buf::IOBuffer=IOBuffer())` => `Vector{String}`
`CSV.readsplitline(source::CSV.Source)` => `Vector{String}`

read a single line from `io` (any `IO` type) as a `Vector{String}` with elements being delimited fields (separated by a delimiter `d`). Can optionally provide a `buf::IOBuffer` type for buffer reuse
"""
function readsplitline end

function readsplitline(io::IO,d::UInt8,q::UInt8,e::UInt8,buf::IOBuffer=IOBuffer())
    vals = String[]
    while !eof(io)
        b = unsafe_read(io, UInt8)
        if b == q
            while !eof(io)
                b = unsafe_read(io, UInt8)
                if b == e
                    Base.write(buf, b)
                    b = unsafe_read(io, UInt8)
                    Base.write(buf, b)
                elseif b == q
                    break
                else
                    Base.write(buf, b)
                end
            end
        elseif b == d
            push!(vals,takebuf_string(buf))
        elseif b == NEWLINE
            break
        elseif b == RETURN
            !eof(io) && unsafe_peek(io) == NEWLINE && unsafe_read(io, UInt8)
            break
        else
            Base.write(buf, b)
        end
    end
    return push!(vals,takebuf_string(buf))
end
readsplitline(io::IO,d=',',q='"',e='\\',buf::IOBuffer=IOBuffer()) = readsplitline(io,UInt8(d),UInt8(q),UInt8(e),buf)
readsplitline(source::CSV.Source) = readsplitline(source.data, source.options.delim, source.options.quotechar, source.options.escapechar)

"""
`CSV.countlines(io::IO, quotechar, escapechar)` => `Int`
`CSV.countlines(source::CSV.Source)` => `Int`

count the number of lines in a file, accounting for potentially embedded newlines in quoted fields
"""
function countlines(io::IO,q::UInt8,e::UInt8)
    nl = 1
    b = 0x00
    while !eof(io)
        b = unsafe_read(io, UInt8)
        if b == q
            while !eof(io)
                b = unsafe_read(io, UInt8)
                if b == e
                    b = unsafe_read(io, UInt8)
                elseif b == q
                    break
                end
            end
        elseif b == NEWLINE
            nl += 1
        elseif b == RETURN
            nl += 1
            !eof(io) && unsafe_peek(io) == NEWLINE && unsafe_read(io, UInt8)
        end
    end
    return nl - (b == NEWLINE || b == RETURN)
end
countlines(io::IO, q, e) = countlines(io, UInt8(q), UInt8(e))
countlines(source::CSV.Source) = countlines(source.data, source.options.quotechar, source.options.escapechar)

function skipto!(f::IO,cur,dest,q,e)
    cur >= dest && return
    for _ = 1:(dest-cur)
        CSV.readline(f,q,e)
    end
    return
end

# NullField is used only during the type detection process
immutable NullField end

# try to infer the type of the value in `val`. The precedence of type checking is `Int` => `Float64` => `Date` => `DateTime` => `String`
function detecttype(val::AbstractString,format,null)
    (val == "" || val == null) && return NullField
    try
        v = parsefield(IOBuffer(replace(val, Char(COMMA), "")), Int)
        !isnull(v) && return Int
    end
    try
        v = CSV.parsefield(IOBuffer(replace(val, Char(COMMA), "")), Float64)
        !isnull(v) && return Float64
    end
    if format != EMPTY_DATEFORMAT
        try # it might be nice to throw an error when a format is specifically given but doesn't parse
            Date(IOBuffer(val),format)
            return Date
        end
        try
            DateTime(IOBuffer(val),format)
            return DateTime
        end
    else
        try
            v = CSV.parsefield(IOBuffer(val), Date)
            !isnull(v) && return Date
        end
        try
            v = parsefield(IOBuffer(val), DateTime)
            !isnull(v) && return DateTime
        end
    end
    return WeakRefString{UInt8}
end


text:
# file test/runtest.jl

include("diagramtest.jl")


text:
using FFTW;

struct Person
    name::String
    age::Int64
end

rfft(A [, dims])


text:
using Test, QSimulator

using LinearAlgebra: I
@testset "Primitives" begin
    r = Resonator("test", 3, QSimulator.ResonatorSpec(5.5))

    # test some trivial cases of the ladder operators
    @test raising(r) == [0 0 0;
                         1 0 0;
                         0 √2 0]
    @test lowering(r) == [0 1 0;
                          0 0 √2;
                          0 0 0]

    r = Resonator("test", 7, QSimulator.ResonatorSpec(5.5))

    # test harmonic oscillator operator identity a†a = N
    @test raising(r) * lowering(r) ≈ number(r)

    # test X = a + a†
    @test X(r) ≈ raising(r) + lowering(r)

    # test phase rotation of 0.25 on X → Y
    @test X(r, 0.25) ≈ Y(r)

    # test commutator relation [a, a†] = 1; because of truncation it doesn't hold for last row/columns
    commutator = lowering(r)*raising(r)- raising(r)*lowering(r)
    @test commutator[1:end-1,1:end-1] ≈ Matrix{eltype(commutator)}(I, dimension(r)-1, dimension(r)-1)

    # test trivial case of the `flip_flop` operator
    a = Resonator("a", 2, QSimulator.ResonatorSpec(5.5))
    b = Resonator("b", 2, QSimulator.ResonatorSpec(5.5))
    @test flip_flop(a,b) == [0 0 0 0;
                             0 0 1 0;
                             0 1 0 0;
                             0 0 0 0]

end


text:
using Documenter, OMETIFF

makedocs(
    authors="Tamas Nagy and contributors",
    sitename="OMETIFF.jl Documentation",
    pages = [
        "Home" => "index.md",
        "Library" => Any[
            "Internals" => joinpath("lib", "internals.md")
        ],
    ],
)

deploydocs(
    repo = "github.com/tlnagy/OMETIFF.jl.git",
    push_preview=true
)