text:
module PDMats

    using LinearAlgebra, SparseArrays, SuiteSparse

    import Base: +, *, \, /, ==, convert, inv, Matrix

    export
        # Types
        AbstractPDMat,
        PDMat,
        PDSparseMat,
        PDiagMat,
        ScalMat,

        # Functions
        dim,
        full,
        whiten,
        whiten!,
        unwhiten,
        unwhiten!,
        pdadd,
        pdadd!,
        add_scal,
        add_scal!,
        quad,
        quad!,
        invquad,
        invquad!,
        X_A_Xt,
        Xt_A_X,
        X_invA_Xt,
        Xt_invA_X,
        test_pdmat

    """
    The base type for positive definite matrices.
    """
    abstract type AbstractPDMat{T<:Real} end

    const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD)

    # source files

    include("chol.jl")
    include("utils.jl")

    include("pdmat.jl")

    if HAVE_CHOLMOD
        include("pdsparsemat.jl")
    end

    include("pdiagmat.jl")
    include("scalmat.jl")

    include("generics.jl")
    include("addition.jl")

    include("testutils.jl")

    include("deprecates.jl")

end # module


text:
using DataFrames: DataFrame
using StructuredLaTeXTables
using Test

# Example 1: Handbook of Writing for the Mathematical Sciences
# source: https://nhigham.com/2019/11/19/better-latex-tables-with-booktabs/

trigmv = "\\trigmv"
trig_expmv = "\\trigexpmv"
trig_block = "\\trigblock"
expleja = "\\expleja"
tols = "\$\\tol=\\tols\$"
told = "\$\\tol=\\told\$"

df1 = DataFrame(
    method=repeat([trigmv, trig_expmv, trig_block, expleja], 2),
    mv=[11034, 21952, 15883, 11180, 15846, 31516, 32023, 17348],
    error=[1.3e-7, 1.3e-7, 5.2e-8, 8.0e-9, 2.7e-11, 2.7e-11, 1.1e-11, 1.5e-11],
    time=[3.9, 6.2, 7.1, 4.3, 5.6, 8.8, 1.4e1, 6.6],
    tol=repeat([tols, told], inner=4),
)

observed1 = make_latex_tabular(df1,
    data=[:mv => raw"$mv$", :error => raw"Rel.~err", :time => raw"Time"],
    metadata=[:method => "", :tol => "not used here"],
    emphasis=:tol, # need emphasis in metadata
    swap_headers=true,
)

# NOTE: entry '14.0' should be '1.4e1' to match original.
# NOTE: Add trailing '\n'?
expected1 = replace(
raw"\begin{tabular}{lcccccc}
\toprule
& \multicolumn{3}{c}{$\tol=\tols$} & \multicolumn{3}{c}{$\tol=\told$} \\
\cmidrule(lr){2-4}\cmidrule(lr){5-7}
& $mv$ & Rel.~err & Time & $mv$ & Rel.~err & Time \\
\midrule
\trigmv & 11034 & 1.3e-7 & 3.9 & 15846 & 2.7e-11 & 5.6 \\
\trigexpmv & 21952 & 1.3e-7 & 6.2 & 31516 & 2.7e-11 & 8.8 \\
\trigblock & 15883 & 5.2e-8 & 7.1 & 32023 & 1.1e-11 & 14.0 \\
\expleja & 11180 & 8.0e-9 & 4.3 & 17348 & 1.5e-11 & 6.6 \\
\bottomrule
\end{tabular}",
"\r\n" => "\n")

@test observed1 == expected1

ext: jl


text:
@deprecate fitDiscrete fitdiscrete
@deprecate phyloNetworklm phylolm
@deprecate sigma2_estim sigma2_phylo
@deprecate mu_estim mu_phylo


text:
# 2D Field tracing on a regular grid.

"""
    bilin_reg(x, y, Q00, Q01, Q10, Q11)

Bilinear interpolation for x1,y1=(0,0) and x2,y2=(1,1)
Q's are surrounding points such that Q00 = F[0,0], Q10 = F[1,0], etc.
"""
function bilin_reg(x, y, Q00, Q10, Q01, Q11)
   fout =
      Q00*(1.0-x)*(1.0-y) +
      Q10* x *    (1.0-y) +
      Q01* y *    (1.0-x) +
      Q11* x * y
end

"""
    grid_interp(x, y, field, ix, iy)

Interpolate a value at (x,y) in a field. `ix` and `iy` are indexes for x,y
locations (0-based).
"""
grid_interp(x, y, field::Array, ix, iy) =
   bilin_reg(x-ix, y-iy,
   field[ix+1, iy+1],
   field[ix+2, iy+1],
   field[ix+1, iy+2],
   field[ix+2, iy+2])

"""
    DoBreak(iloc, jloc, iSize, jSize)

Check to see if we should break out of an integration.
"""
function DoBreak(iloc, jloc, iSize, jSize)
   ibreak = false
   if iloc ≥ iSize-1 || jloc ≥ jSize-1; ibreak = true end
   if iloc < 0 || jloc < 0; ibreak = true end
   ibreak
end

"Create unit vectors of field."
function normalize_field(ux, uy, dx, dy)
   fx, fy = similar(ux), similar(uy)
   dxInv, dyInv = 1/dx, 1/dy
   @inbounds for i = 1:length(ux)
      uInv = 1.0 / hypot(ux[i]*dxInv, uy[i]*dyInv)
      fx[i] = ux[i] * dxInv * uInv
      fy[i] = uy[i] * dyInv * uInv
   end
   fx, fy
end

"""
    euler(maxstep, ds, startx, starty, xGrid, yGrid, ux, uy)

Fast 2D tracing using Euler's method. It takes at most `maxstep` with step size `ds` tracing
the vector field given by `ux,uy` starting from `(startx,starty)` in the Cartesian grid
specified by ranges `xGrid` and `yGrid`. Step size is in normalized coordinates within the
range [0, 1]. Return footprints' coordinates in (`x`, `y`).
"""
function euler(maxstep, ds, startx, starty, xGrid, yGrid, ux, uy)

   @assert size(ux) == size(uy) "field array sizes must be equal!"

   x = Vector{eltype(startx)}(undef, maxstep)
   y = Vector{eltype(starty)}(undef, maxstep)

   iSize, jSize = size(xGrid,1), size(yGrid,1)

   # Get starting points in normalized/array coordinates
   dx = xGrid[2] - xGrid[1]
   dy = yGrid[2] - yGrid[1]
   x[1] = (startx - xGrid[1]) / dx
   y[1] = (starty - yGrid[1]) / dy

   # Create unit vectors from full vector field
   f1, f2 = normalize_field(ux, uy, dx, dy)

   nstep = 0
   # Perform tracing using Euler's method
   for n = 1:maxstep-1
      # Find surrounding points
      ix = floor(Int, x[n])
      iy = floor(Int, y[n])

      # Break if we leave the domain
      if DoBreak(ix, iy, iSize, jSize)
         nstep = n; break
      end

      # Interpolate unit vectors to current location
      fx = grid_interp(x[n], y[n], f1, ix, iy)
      fy = grid_interp(x[n], y[n], f2, ix, iy)

      if isnan(fx) || isnan(fy) || isinf(fx) || isinf(fy)
         nstep = n
         break
      end

      @muladd x[n+1] = x[n] + ds * fx
      @muladd y[n+1] = y[n] + ds * fy

      nstep = n
   end

   # Convert traced points to original coordinate system.
   @inbounds for i = 1:nstep
      @muladd x[i] = x[i]*dx + xGrid[1]
      @muladd y[i] = y[i]*dy + yGrid[1]
   end
   x[1:nstep], y[1:nstep]
end

"""
    rk4(maxstep, ds, startx, starty, xGrid, yGrid, ux, uy)

Fast and reasonably accurate 2D tracing with 4th order Runge-Kutta method and constant step
size `ds`. See also [`euler`](@ref).
"""
function rk4(maxstep, ds, startx, starty, xGrid, yGrid, ux, uy)

   @assert size(ux) == size(uy) "field array sizes must be equal!"

   x = Vector{eltype(startx)}(undef, maxstep)
   y = Vector{eltype(starty)}(undef, maxstep)

   iSize, jSize = size(xGrid,1), size(yGrid,1)

   # Get starting points in normalized/array coordinates
   dx = xGrid[2] - xGrid[1]
   dy = yGrid[2] - yGrid[1]
   x[1] = (startx - xGrid[1]) / dx
   y[1] = (starty - yGrid[1]) / dy

   # Create unit vectors from full vector field
   fx, fy = normalize_field(ux, uy, dx, dy)

   nstep = 0
   # Perform tracing using RK4
   for n = 1:maxstep-1
      # SUBSTEP #1
      ix = floor(Int, x[n])
      iy = floor(Int, y[n])
      if DoBreak(ix, iy, iSize, jSize); nstep = n; break end

      f1x = grid_interp(x[n], y[n], fx, ix, iy)
      f1y = grid_interp(x[n], y[n], fy, ix, iy)
      if isnan(f1x) || isnan(f1y) || isinf(f1x) || isinf(f1y)
         nstep = n; break
      end
      # SUBSTEP #2
      @muladd xpos = x[n] + f1x*ds*0.5
      @muladd ypos = y[n] + f1y*ds*0.5
      ix = floor(Int, xpos)
      iy = floor(Int, ypos)
      if DoBreak(ix, iy, iSize, jSize); nstep = n; break end

      f2x = grid_interp(xpos, ypos, fx, ix, iy)
      f2y = grid_interp(xpos, ypos, fy, ix, iy)

      if isnan(f2x) || isnan(f2y) || isinf(f2x) || isinf(f2y)
         nstep = n; break
      end
      # SUBSTEP #3
      @muladd xpos = x[n] + f2x*ds*0.5
      @muladd ypos = y[n] + f2y*ds*0.5
      ix = floor(Int, xpos)
      iy = floor(Int, ypos)
      if DoBreak(ix, iy, iSize, jSize); nstep = n; break end

      f3x = grid_interp(xpos, ypos, fx, ix, iy)
      f3y = grid_interp(xpos, ypos, fy, ix, iy)
      if isnan(f3x) || isnan(f3y) || isinf(f3x) || isinf(f3y)
         nstep = n; break
      end

      # SUBSTEP #4
      @muladd xpos = x[n] + f3x*ds
      @muladd ypos = y[n] + f3y*ds
      ix = floor(Int, xpos)
      iy = floor(Int, ypos)
      if DoBreak(ix, iy, iSize, jSize); nstep = n; break end

      f4x = grid_interp(xpos, ypos, fx, ix, iy)
      f4y = grid_interp(xpos, ypos, fy, ix, iy)
      if isnan(f4x) || isnan(f4y) || isinf(f4x) || isinf(f4y)
         nstep = n; break
      end

      # Peform the full step using all substeps
      @muladd x[n+1] = x[n] + ds/6.0 * (f1x + f2x*2.0 + f3x*2.0 + f4x)
      @muladd y[n+1] = y[n] + ds/6.0 * (f1y + f2y*2.0 + f3y*2.0 + f4y)

      nstep = n
   end

   # Convert traced points to original coordinate system.
   @inbounds for i = 1:nstep
      @muladd x[i] = x[i]*dx + xGrid[1]
      @muladd y[i] = y[i]*dy + yGrid[1]
   end
   x[1:nstep], y[1:nstep]
end

"""
	 trace2d_rk4(fieldx, fieldy, startx, starty, gridx, gridy;
		 maxstep=20000, ds=0.01, gridtype="ndgrid", direction="both")

Given a 2D vector field, trace a streamline from a given point to the edge of the vector
field. The field is integrated using Runge Kutta 4. Slower than Euler, but more accurate.
The higher accuracy allows for larger step sizes `ds`. Step size is in normalized
coordinates within the range [0,1]. See also [`trace2d_euler`](@ref).
"""
function trace2d_rk4(fieldx, fieldy, startx, starty, gridx, gridy;
   maxstep=20000, ds=0.01, gridtype="ndgrid", direction="both")

   @assert ndims(gridx) == 1 "Grid must be given in 1D range or vector!"

   if gridtype == "ndgrid"
      fx = fieldx
      fy = fieldy
   else # meshgrid
      fx = permutedims(fieldx)
      fy = permutedims(fieldy)
   end

   if direction == "forward"
      xt, yt = rk4(maxstep, ds, startx, starty, gridx, gridy, fx, fy)
   elseif direction == "backward"
      xt, yt = rk4(maxstep, ds, startx, starty, gridx, gridy, -fx, -fy)
   else
      x1, y1 = rk4(floor(Int,maxstep/2), ds, startx, starty, gridx, gridy,
         -fx, -fy)
      blen = length(x1)
      x2, y2 = rk4(maxstep-blen, ds, startx, starty, gridx, gridy, fx, fy)
      # concatenate with duplicates removed
      xt = vcat(reverse!(x1), x2[2:end])
      yt = vcat(reverse!(y1), y2[2:end])
   end

   xt, yt
end

"""
	 trace2d_euler(fieldx, fieldy, startx, starty, gridx, gridy;
		 maxstep=20000, ds=0.01, gridtype="ndgrid", direction="both")

Given a 2D vector field, trace a streamline from a given point to the edge of the vector
field. The field is integrated using Euler's method, which is faster but less accurate than
RK4. Only valid for regular grid with coordinates' range `gridx` and `gridy`. Step size is
in normalized coordinates within the range [0,1]. The field can be in both `meshgrid` or
`ndgrid` (default) format. Supporting `direction` of {"both","forward","backward"}.
"""
function trace2d_euler(fieldx, fieldy, startx, starty, gridx, gridy;
   maxstep=20000, ds=0.01, gridtype="ndgrid", direction="both")

   @assert ndims(gridx) == 1 "Grid must be given in 1D range or vector!"

   if gridtype == "ndgrid"
      fx = fieldx
      fy = fieldy
   else # meshgrid
      fx = permutedims(fieldx)
      fy = permutedims(fieldy)
   end

   if direction == "forward"
      xt, yt = euler(maxstep, ds, startx, starty, gridx, gridy, fx, fy)
   elseif direction == "backward"
      xt, yt = euler(maxstep, ds, startx, starty, gridx, gridy, -fx, -fy)
   else
      x1, y1 = euler(floor(Int,maxstep/2), ds, startx, starty, gridx, gridy,
         -fx, -fy)
      blen = length(x1)
      x2, y2 = euler(maxstep-blen, ds, startx, starty, gridx, gridy, fx, fy)
      # concatenate with duplicates removed
      xt = vcat(reverse!(x1), x2[2:end])
      yt = vcat(reverse!(y1), y2[2:end])
   end

   xt, yt
end

"""
    trace2d_euler(fieldx, fieldy, startx, starty, grid::CartesianGrid; kwargs...)
"""
function trace2d_euler(fieldx, fieldy, startx, starty, grid::CartesianGrid;
   kwargs...)

   gridmin = coordinates(minimum(grid))
   gridmax = coordinates(maximum(grid))
   Δx = spacing(grid)

   gridx = range(gridmin[1], gridmax[1], step=Δx[1])
   gridy = range(gridmin[2], gridmax[2], step=Δx[2])

   trace2d_euler(fieldx, fieldy, startx, starty, gridx, gridy; kwargs...)
end


text:
push!(LOAD_PATH,"../src/") # load module from local directory

using PolaronMobility, Documenter

makedocs(sitename="PolaronMobility.jl documentation")

deploydocs(repo="github.com/jarvist/PolaronMobility.jl.git",)


text:
module Smoothing

export binomial
"""
    binomial(data::Vector{any}, Nₚ::Int)

Takes a 1d data array and smooths the data with a 2Nₚ+1 binomial filter.
This is computed using the method outlined in :

*Binomial smoothing filter: A way to avoid some pitfalls of least‐squares
polynomial smoothing* in Review of Scientific Instruments 54, 1034 (1983);
https://doi.org/10.1063/1.1137498

or see

https://aip.scitation.org/doi/abs/10.1063/1.1137498

The method implements a 2Nₚ+1 binomial filter by Nₚ passes of the three point
binomial filter; and this is computationally implemented by two passess of the
{1,1}/2 smoothing implemented below. The filter preserves the values starting
and ending values of the initial data set.

# Examples (you can easily manually verify data1's smoothing is correct.)
```julia-repl
julia> data1 = [1.0, 2.0, 3.0, 4.0, 5.0]
julia> data2 = [3.0, 6.0, 1.0, 4.0, 0.0]
julia> binomial(data1, 1)

5-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0
 5.0

julia> binomial(data2, 20)

5-element Vector{Float64}:
 3.0
 2.3162835515104234
 1.5937387603335083
 0.8162830746732652
 0.0
```
"""
function binomial(data::Vector, Nₚ::Int)
    y = data
    for pass in 1:Nₚ
        x = zeros(length(data))

        for i in 1:length(y)
            if i==1
                x[i] = (y[i] + y[i+1])/2.0
            elseif i>1 && i<length(y)
                x[i] = (y[i] + y[i+1])/2.0
            end
        end

        y = zeros(length(data))

        for j in 1:length(x)
            if j==1
                y[j] = data[j]
            elseif j>1 && j<length(x)
                y[j] = (x[j-1] + x[j])/2.0
            else
                y[j] = data[j]
            end
        end

    end

    return y
end

end # module


text:
"""
    Berlage(; <keyword arguments>)

Create a Berlage wavelet.

# Arguments

**Keyword arguments**

* `dt::Real=0.002`: sampling interval in secs.
* `f0::Real=20.0`: central frequency in Hz.
* `m::Real=2`: exponential parameter of Berlage wavelet.
* `alpha::Real=180.0`: alpha parameter of Berlage wavelet in rad/secs.
* `phi0::Real`: phase rotation in radians.

# Example
```julia
julia> w = Berlage(); plot(w);
```

**Reference**

* Aldridge, David F., 1990, The berlage wavelet: GEOPHYSICS, 55, 1508--1511.
"""

function Berlage(; dt::Real=0.002, f0::Real=20.0, m::Real=2, alpha::Real=180.0,
                 phi0::Real=0.0)

    nw = floor(Int, 2.2/(f0*dt))
    t = dt*collect(0:1:nw-1)
    w = (t.^m).*exp(-alpha*t).*cos(2*pi*f0*t + phi0);
    w = w/maximum(w)

end


text:
function pluto_to_bookdown(pluto)
    #Read the Pluto notebook
    str = open(f->read(f, String), pluto)
    #Convert the string file in a list, each element of  the list contains a Pluto cell
    list = split(str, "# ╔═╡ ")
    popfirst!(list)

    #The last element of the list contains the order in which the cells should be rendered
    cell_order_str = pop!(list)
    cell_order = split(cell_order_str,"# ")[2:end]

    # We obtain a list of the cells order
    for i in 1:length(cell_order)
        x =findfirst(r"[a-zA-Z0-9]",cell_order[i])[1]
        cell_order[i] = cell_order[i][x:x+35]
    end
    #create a dictionary that match: cell Id => cell content
    cells_dict = Dict(list[1][1:36] => list[1][37:end])
    for cell in list
        cells_dict[cell[1:36]] = cell[37:end]
    end

    #Create de bookdown string that to store the bookdown content
    #All the bookdown will start the same way
    bookdown = "# {TITLE}\n\n```{julia, echo=FALSE}\nusing Markdown\nusing InteractiveUtils\n```\n\n"

    # Iterate through all the Pluto cells in order and handle every cell depending of their content
    for id in cell_order

        cell = cells_dict[id]

        # Handle md pluto cells
        if contains(cell,"md\"")

            if contains(cell,"md\"\"\"")
                start = findfirst("md\"\"\"",cell)[end] + 1
                last = findlast("\"\"\"",cell)[1] - 1
            elseif contains(cell,"md\"")
                start = findfirst("md\"",cell)[end] + 1
                last = findlast("\"",cell)[1] - 1
            end

            if contains(cell,"\";")
            bookdown = bookdown * "```{julia, echo=FALSE} \r\n"* cell[start:last] * "```\n"
            else
                bookdown = bookdown * cell[start:last] * "\n"
            end

        # Handle code cells
        else

            if contains(cell,"end;")
            bookdown = bookdown * "```{julia,results = FALSE} \n" * cell * "```\n"
            else
            bookdown = bookdown *  "```{julia} \r\n"* cell * "```\n"
            end

        end
    end

    #Create the Rmd file
    open(pluto[1:end-3]*".Rmd", "w") do f
        write(f, bookdown)
    end

    return pluto[1:end-3]*".Rmd" * " successfully created."
end

max_stars_repo_name: lambdaclass/pluto-to-bookdown

max_issues_repo_name: lambdaclass/pluto-to-bookdown

max_forks_repo_name: lambdaclass/pluto-to-bookdown


text:
mutable struct BufferWord{T} <: AbstractWord{T}
    storage::Vector{T}
    lidx::Int
    ridx::Int

    function BufferWord{T}(v::AbstractVector{<:Integer},
        freeatbeg=8, freeatend=8) where T

        storage = Vector{T}(undef, freeatbeg + length(v) + freeatend)
        bw = new{T}(storage, freeatbeg+1, freeatbeg + length(v))

        # placing content in the middle
        # storage[freeatbeg+1:end-freeatend] .= v
        @inbounds for i in 1:length(v)
            bw[i] = v[i]
        end

        return bw
    end

    function BufferWord{T}(freeatbeg, freeatend) where T
        l = freeatbeg + freeatend
        return new{T}(Vector{T}(undef, l), freeatbeg+1, freeatbeg)
    end
end

BufferWord(v::Union{<:Vector{<:Integer},<:AbstractVector{<:Integer}}) =
    BufferWord{UInt16}(v)

# helper functions for growing storage:

internal_length(bw::BufferWord) = length(bw.storage)

function _growatbeg!(bw::BufferWord, k::Integer)
    @assert k ≥ 0
    k == 0 && return bw
    resize!(bw.storage, internal_length(bw) + k)

    @inbounds for i in 0:length(bw)-1
        bw.storage[bw.ridx+k-i] = bw.storage[bw.ridx-i]
    end

    bw.ridx += k
    bw.lidx += k

    return bw
end

_growatend!(bw::BufferWord, k::Integer) =
    resize!(bw.storage, internal_length(bw) + k)

# AbstractWord Interface:

Base.@propagate_inbounds function Base.getindex(bw::BufferWord, n::Integer)
    @boundscheck checkbounds(bw, n)
    return @inbounds bw.storage[bw.lidx+n-1]
end

Base.@propagate_inbounds function Base.setindex!(bw::BufferWord, val, n::Integer)
    @boundscheck checkbounds(bw, n)
    return @inbounds bw.storage[bw.lidx+n-1] = val
end
Base.length(bw::BufferWord) = length(bw.lidx:bw.ridx)

function Base.push!(bw::BufferWord, k::Integer)
    if internal_length(bw) == bw.ridx
        _growatend!(bw, max(length(bw), 16))
    end
    bw.ridx += 1
    @inbounds bw[end] = k
    return bw
end

function Base.pushfirst!(bw::BufferWord, k::Integer)
    if bw.lidx ≤ firstindex(bw.storage)
        _growatbeg!(bw, max(length(bw), 16))
    end
    bw.lidx -= 1
    @inbounds bw[1] = k
    return bw
end

function Base.pop!(bw::BufferWord)
    @assert !isempty(bw)
    @inbounds val = bw[end]
    bw.ridx -= 1
    return val
end

function Base.popfirst!(bw::BufferWord)
    @assert !isempty(bw)
    @inbounds val = bw[1]
    bw.lidx +=1
    return val
end

function Base.prepend!(bw::BufferWord, w::AbstractVector)

    free_space = bw.lidx - 1
    lw = length(w)

    if (free_space - lw) ≤ 0
        _growatbeg!(bw, lw)
    end

    @inbounds for i in 1:lw
        bw.storage[bw.lidx-lw+i-1] = w[i]
    end
    bw.lidx -= lw

    return bw
end

function Base.append!(bw::BufferWord, w::AbstractVector)

    free_space = internal_length(bw) - bw.ridx
    lw = length(w)

    if (free_space - lw) ≤ 0
        _growatend!(bw, lw)
    end

    @inbounds for i in 1:lw
        bw.storage[bw.ridx+i] = w[i]
    end
    bw.ridx += lw

    return bw
end

function Base.resize!(bw::BufferWord, nl::Integer)
    l = length(bw)
    if nl > l
        free_space = internal_length(bw) - bw.ridx
        if nl-l > free_space
            _growatend!(bw, nl-l)
        end
    elseif nl != l
        if nl < 0
            throw(ArgumentError("new length must be ≥ 0"))
        end
    end
    bw.ridx += nl-l
    return bw
end

function Base.:*(bw::BufferWord{S}, bv::BufferWord{T}) where {S, T}
    res = BufferWord{promote_type(S,T)}(
        bw.lidx-1+length(bw) + length(bv),
        internal_length(bv) - bv.ridx
    )
    res.lidx = bw.lidx
    res.ridx = res.lidx+length(bw)+length(bv)-1
    @inbounds for i in 1:length(bw)
        res[i] = bw[i]
    end
    @inbounds for i in 1:length(bv)
        res[length(bw)+i] = bv[i]
    end
    return res
end

Base.similar(bw::BufferWord, ::Type{S}) where S =
    (l = internal_length(bw)÷2; BufferWord{S}(l, l))

function Base.empty!(bw::BufferWord)
    bw.lidx = 1
    bw.ridx = 0
    return bw
end

# performance methods overloaded from AbstractWord:
# one less allocation:
Base.one(::Type{BufferWord{T}}) where T = BufferWord{T}(8, 8)

# @view constructor
Base.view(bw::BufferWord, u::UnitRange{Int}) =
    KnuthBendix.SubWord(view(bw.storage, u.+(bw.lidx-1)))


text:
using JLD2, Distributions, Serialization
import GEAS:sim_pt_QTL, breeding_program, top_QTL

"""
    function the_simulation(jld, nqtl, rt, nrpt; rst = "")
---
## Description
- 3 comparisons
  - Normal genomic selection, no known QTL
  - Known 5 QTL as fixed effects
  - Edit the biggest of known 5 QTL in generation 2-6

## parameters
1. jld: the base
2. nqtl: e.g., 500
3. rt: correlation of QTL effects between 2 traits
4. nrpt: number of repeats
5. rst: result path, in `2021-06-30/`
"""
function the_simulation(jld, nqtl, rt, nrpt; rst = "")
    @load jld base
    if length(rst) > 1
        isdir(rst) || mkpath(rst)
    end

    w = ndigits(nrpt)
    for r in 1:nrpt
        d = MvNormal(zeros(2), [1 rt; rt 1])
        qtl = sim_pt_QTL(base, nqtl, d)

        suffix = lpad(string(r), w, '0')
        serialize(joinpath(rst, "qtl-$suffix.ser"), qtl)
        par = begin
            Parameters = Dict(
                :nSire=> 100,
                :nDam => 200,   # male:female = 1:2
                :nSib => 40,
                :nC7e => 30,    # number of sib to be challenged
                :nG8n => 11,    # number of generations
                :nQTL => [nqtl, nqtl],
                :h²   => [.5, .5],
                :p8e  => .5,    # death rate after challenge
                :e19e => 1.,    # edit_successsful_rate
                :w4t  => 1,     # weight on the binary trait EBV
                :nk3n => 5,     # number of known QTL on binary trait
                :dd   => 0.01, # A small value to be added to diagonals
                :log  => joinpath(rst, "qtl.log")
            )
            (; Parameters...)       # named tuple.  contents as above
        end

        # genomic selection with no known QTL
        prd, snp = breeding_program(base, par, qtl)
        serialize(joinpath(rst, "prd-a-$suffix.ser"), prd)
        serialize(joinpath(rst, "snp-a-$suffix.ser"), snp)

        # genomic selection with known QTL as fixed effects
        prd, snp = breeding_program(base, par, qtl, fixed=true)
        serialize(joinpath(rst, "prd-b-$suffix.ser"), prd)
        serialize(joinpath(rst, "snp-b-$suffix.ser"), snp)

        # genomic selection, edit biggest known QTL each generation
        prd, snp = breeding_program(base, par, qtl, edit=true)
        serialize(joinpath(rst, "prd-c-$suffix.ser"), prd)
        serialize(joinpath(rst, "snp-c-$suffix.ser"), snp)
    end
end

function simulation_2021_06_30()
    the_simulation("dat/run/base.jld2", 500, 0, 10, rst="dat/2021-06-30")
end

function simulation_2021_07_01()
    the_simulation("dat/run/base.jld2", 500, -.2, 10, rst="dat/2021-07-01")
end

function simulation_2021_07_02()
    the_simulation("dat/run/base.jld2", 500, .2, 10, rst="dat/2021-07-02")
end


text:
__precompile__()
module AprilTags

function __init__()
    depfile = joinpath(dirname(@__FILE__),"../deps/loadpath.jl")
    isfile(depfile) ? include(depfile) : error("AprilTags.jl not properly installed. Please run: Pkg.build(\"AprilTags\")")
end

using LinearAlgebra, Statistics
using Colors, ImageDraw, FixedPointNumbers
import Base.convert

export
#helpers
AprilTag,
AprilTagDetector,
freeDetector!,
getTagDetections,
homography_to_pose,
homographytopose,
threadcalldetect,
getAprilTagImage,
detectAndPose,
tagOrthogonalIteration,
# wrappers
apriltag_detector_create,
tag36h11_create,
tag36h11_destroy,

apriltag_detector_add_family,
apriltag_detector_detect,
apriltag_detections_destroy,
apriltag_detector_destroy,
threadcall_apriltag_detector_detect,

#drawing and plotting
drawTagBox!,
drawTagAxes!

include("wrapper.jl")
include("helpers.jl")
include("tagdraw.jl")

end # module


text:
using makePackage
using Test

@testset "makePackage.jl" begin
    # Write your tests here.
end


text:
"""
One dimensional chain.
"""
struct Chain <: AbstractCubicLattice
    sites::Int
    neighs::Matrix{Int} # row = right, left; col = siteidx

    # for generic checkerboard decomposition
    n_bonds::Int
    bonds::Matrix{Int} # src, trg, type
end

# constructors
"""
    Chain(nsites::Int)

Create a chain with `nsites`.
"""
function Chain(nsites::Int)
    neighs = build_neighbortable(Chain, nsites)

    # for generic checkerboard decomposition
    n_bonds = nsites
    bonds = zeros(n_bonds, 3)
    bondid = 1
    for src in 1:nsites
        nright = neighs[1, src]
        bonds[bondid,:] .= [src,nright,0]
        bondid += 1
    end

    return Chain(nsites, neighs, n_bonds, bonds)
end

function build_neighbortable(::Type{Chain}, nsites::Int)
    c = 1:nsites
    right = circshift(c,-1)
    left = circshift(c,1)
    return vcat(right[:]',left[:]')
end

@inline nsites(c::Chain) = c.sites
@inline neighbors_lookup_table(c::Chain) = c.neighs

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
import Base: fill!, rand, randn
import Random: seed!, rand!, randn!
import Adapt: adapt
export distribute

function fill!(a::MPIArray, x)
    forlocalpart!(y -> fill!(y, x), a)
    a
end

"""
    distribute(a::AbstractArray; root=0, A=nothing)
distributes a dense array.
"""
function distribute(a::AbstractArray{T,N}; root=0, A=nothing) where {T, N}
    if A == nothing
        TA = typeof(a)
        A = hasproperty(TA, :name) ? TA.name.wrapper : TA
    end
    reqs = MPI.Request[]
    size_a = MPI.bcast(size(a), root, MPI.COMM_WORLD)
    rslt = MPIArray{T,N,A}(undef, size_a...)
    if Rank() == root
        for r = 0:Size() - 1
            r == Rank() && continue
            sync()
            push!(reqs, Isend(@view(a[repeat([:], N-1)..., rslt.partitioning[r+1][end]]), r))
        end
        rslt.localarray .= @view(a[repeat([:], N-1)..., rslt.partitioning[Rank()+1][end]])
    else
        sync()
        push!(reqs, Irecv!(rslt.localarray, root))
    end
    for r in reqs
        MPI.Wait!(r)
    end
    rslt
end

"""
    distribute(a::AbstractSparseMatrix; root=0, A=Array)
distributes a sparse array into a dense MPIMatrix.
"""
function distribute(a::AbstractSparseMatrix; root=0, A=Array)
    T = eltype(a)
    reqs = MPI.Request[]
    size_a = MPI.bcast(size(a), root, MPI.COMM_WORLD)
    rslt = MPIMatrix{T,A}(undef, size_a[1], size_a[2])
    if Rank() == root
        for r = 0:Size() - 1
            r == Rank() && continue
            sync()
            push!(reqs, MPI.isend(a[:, rslt.partitioning[r+1][2]], r, 0, MPI.COMM_WORLD))
        end
        local_sparse = a[:, rslt.partitioning[Rank()+1][2]]
        rslt.localarray .= A(local_sparse)
    else
        sync()
        local_sparse = MPI.recv(root, 0, MPI.COMM_WORLD)[1]
        rslt.localarray .= A(local_sparse)
    end

    for r in reqs
        MPI.Wait!(r)
    end
    rslt
end

for (fname!, fname) in [(:rand!, :rand), (:randn!, :randn)]
    @eval begin
        function $fname!(a::MPIArray{T,N,A};seed=nothing, common_init=false, root=0) where {T,N,A}
            # for both local initialization and initialization from master, from Array{Float64} (for the sake of reproducibility)
            if seed ≠ nothing
                # Rank is added here in order to take different samples for each rank
                common_init ? seed!(seed) : seed!(seed + Rank())
            end
            if !common_init
                forlocalpart!(y -> $fname!(y), a)
                a
            else
                tmp = MPIArray{Float64, N, Array}(undef, size(a)...)
                if Rank() == root
                    src = $fname(size(a)...) # Intentionally creating Array{Float64}
                else
                    src = zeros(zeros(Int, N)...)
                end
                tmp = distribute(src; root=root)
                forlocalpart!(tmp, a) do x, y
                    y .= adapt(A{T}, x)
                end
                a
            end
        end
    end
end


text:
Base.summary(p::DiscreteDynamicProblem) =
    string(TYPE_COLOR, nameof(typeof(p)),
    NO_COLOR, " with ",
    TYPE_COLOR, length(p.tStateVectors),
    NO_COLOR, " state variables and ",
    TYPE_COLOR, length(p.tChoiceVectors),
    NO_COLOR, " choice variable(s) ",
    )

function Base.show(io::IO, p::DiscreteDynamicProblem)
  println(io,summary(p))
end

TreeViews.hastreeview(x::StructuralDDP.DiscreteDynamicProblem) = true
function TreeViews.treelabel(io::IO,x::StructuralDDP.DiscreteDynamicProblem,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Base.Text(Base.summary(x)))
end


text:
@testset "area.jl" begin

    angular_distance_x = 34.49
    angular_distance_y = 28.11
    angular_distance_z = 29.0

    @test area(angular_distance_z, angular_distance_y, angular_distance_x) ≈ 0.1239 atol = 0.01

    p1a = Point(0.0, 0.0)
    p2a = Point(0.0, 90.0)
    p3a = Point(89.99, 0.0)

    @test area(p1a, p2a, p3a) ≈ 1.5708 atol = 0.001
    @test area(p1a, Point(0.0, 179.99), p3a) ≈ 3.1416 atol = 0.001

    pi = Point(32.0, 47.0)
    poly = Polygon(pi, [p1a, p2a, p3a])

    @test area(poly) ≈ 1.5708 atol = 0.001
    @test area(Polygon(pi, [p1a, Point(0.0, 179.99), p3a])) ≈ 3.1416 atol = 0.001
    @test area(Polygon(pi, [p1a, p2a, p3a, Point(0.0, 270.001), Point(-89.99, 0.0)])) ≈ 4.7124 atol = 0.001
    @test area(Polygon(pi, [p2a, p3a, Point(0.0, 270.001), Point(-89.99, 0.0)])) ≈ 2*π atol = 0.001

    #concave polygon
    p1 = Point(64.0, 27.0)
    p2 = Point(47.0, 0.0)
    p3 = Point(42.0, 25.0)
    p4 = intermediate_point(p2,p3,0.33)
    p5 = Point(52.0, 16.0)
    p6 = intermediate_point(p2,p3,0.67)
    p7 = Point(55.0, -7.0)
    p8 = intermediate_point(p2,p1,0.33)
    p9 = Point(62.0, 4.0)
    p10 = intermediate_point(p2,p1,0.67)
    p11 = Point(56.0, 21.0)

    ## triangle areas
    at1 = area(p1,p2,p3)
    at2 = area(p4,p5,p6)
    at3 = area(p2,p7,p8)
    at4 = area(p8,p9,p10)

    ## area polygon
    expected_area_polygon = at1-at2+at3+at4

    polygon_complex_clockwise = Polygon(p11, [p1,p3,p6,p5,p4,p2,p7,p8,p9,p10])
    polygon_complex_counterclockwise = Polygon(p11, [p10,p9,p8,p7,p2,p4,p5,p6,p3,p1])

    @test area(polygon_complex_clockwise) ≈ expected_area_polygon atol = 0.001
    @test area(polygon_complex_counterclockwise) ≈ expected_area_polygon atol = 0.001
end


text:
# An exception type thrown when a method stub is called.
type NotImplementedException <: Exception end

# Exception types thrown when a Vertex is added to a graph that already belongs to a graph
type VertexAlreadyBelongsToAnotherGraphException <: Exception end
type VertexAlreadyBelongsToSpecifiedGraphException <: Exception end

# Exception types thrown when an Edge is added to a graph that already belongs to a graph
type EdgeAlreadyBelongsToAnotherGraphException <: Exception end
type EdgeAlreadyBelongsToSpecifiedGraphException <: Exception end

# An exception type thrown when the tail of an edge being added does not belong to the graph
type EdgeTailDoesNotBelongToGraphException <: Exception end

# An exception type thrown when the head of an edge being added does not belong to the graph
type EdgeHeadDoesNotBelongToGraphException <: Exception end

# Exceptions thrown on invalid query operation sequence, or invalid query state
type InvalidQueryOptionOutputTypeCombinationException <: Exception end
type InvalidInputQueryResultTypeException <: Exception end
type InvalidOutputQueryResultTypeException <: Exception end
type MissingPreviousQueryContextException <: Exception end
type MissingPreviousQueryContextResultException <: Exception end
type InvalidQueryOptionOutputTypeCombinationException <: Exception end
type InvalidInputOutputQueryOptionCombinationException <: Exception end
type CompositeQueryResultNotRealisedException <: Exception end
type InvalidStepCountForMergeException <: Exception end
type InvalidStepCountForResultCountException <: Exception end
type UnexpectedTypeOnMoveToVertices <: Exception end
type UnexpectedTypeOnMoveToEdges <: Exception end

type EdgeDoesNotBelongToGraphException <: Exception end
type VertexDoesNotBelongToGraphException <: Exception end


text:
s1 = "12.125"; v1 = 12.125
s2 = "-0.125"; v2 = -0.125
s3 = "5"     ; v3 = 5

@testset "parse $T" for T in (Double16, Double32, Double64)

@test T(s1) === T(v1)
@test T(s2) === T(v2)
@test T(s3) === T(v3)

end

@testset "T_str" begin

@test d64"12.125" === Double64(v1)
@test d32"12.125" === Double32(v1)
@test d16"12.125" === Double16(v1)

end


text:
function initialize(r, tuning::AbstractParticleSwarm)
    return initialize(tuning.rng, r, tuning.n_particles)
end

function retrieve!(state::ParticleSwarmState, tuning::AbstractParticleSwarm)
    return retrieve!(tuning.rng, state)
end

function pbest!(state::ParticleSwarmState, measurements, tuning::AbstractParticleSwarm)
    return pbest!(state, measurements, tuning.prob_shift)
end


text:
"""
    MainSolver <: AbstractSolver

Main solver. Handle the solving of a model, and optional multithreaded and/or distributed subsolvers.

# Arguments:
- `model::Model`: A formal description of the targeted problem
- `state::_State`: An internal state to store the info necessary to a solving run
- `options::Options`: User options for this solver
- `subs::Vector{_SubSolver}`: Optional subsolvers
"""
mutable struct MainSolver <: MetaSolver
    meta_local_id::Tuple{Int, Int}
    model::_Model
    options::Options
    pool::Pool
    rc_report::RemoteChannel
    rc_sol::RemoteChannel
    rc_stop::RemoteChannel
    remotes::Dict{Int, Future}
    state::State
    status::Symbol
    strategies::MetaStrategy
    subs::Vector{_SubSolver}
    time_stamps::TimeStamps
end

make_id(::Int, id, ::Val{:lead}) = (id, 0)

function solver(model = model();
    options = Options(),
    pool = pool(),
    strategies = MetaStrategy(model),
)
    mlid = (1, 0)
    rc_report = RemoteChannel(() -> Channel{Nothing}(length(workers())))
    rc_sol = RemoteChannel(() -> Channel{Pool}(length(workers())))
    rc_stop = RemoteChannel(() -> Channel{Nothing}(1))
    remotes = Dict{Int, Future}()
    subs = Vector{_SubSolver}()
    ts = TimeStamps(model)
    return MainSolver(mlid, model, options, pool, rc_report, rc_sol, rc_stop, remotes, state(), :not_called, strategies, subs, ts)
end

# Forwards from TimeStamps
@forward MainSolver.time_stamps add_time!, time_info, get_time

"""
    empty!(s::Solver)

"""
function Base.empty!(s::MainSolver)
    empty!(s.model)
    s.state = state()
    empty!(s.subs)
    # TODO: empty remote solvers
end

"""
    status(solver)
Return the status of a MainSolver.
"""
status(s::MainSolver) = s.status

function _init!(s::MainSolver)
    _init!(s, :global)
    _init!(s, :remote)
    _init!(s, :meta)
    _init!(s, :local)
end

function stop_while_loop(s::MainSolver, ::Atomic{Bool}, iter, start_time)
    remote_limit = isready(s.rc_stop) # Add ! when MainSolver is passive
    iter_limit = iter < get_option(s, "iteration")
    time_limit = time() - start_time < get_option(s, "time_limit")
    if !remote_limit
        s.status = :solution_limit
        return false
    end
    if !iter_limit
        s.status = :iteration_limit
        return false
    end
    if !time_limit
        s.status = :time_limit
        return false
    end
    return true
end

function remote_dispatch!(s::MainSolver)
    for (w, ls) in s.remotes
        remote_do(solve!, w, fetch(ls))
    end
end

function post_process(s::MainSolver)
    path = get_option(s, "info_path")
    sat = is_sat(s)
    if s.status == :not_called
        s.status = :solution_limit
    end
    if !isempty(path)
        info = Dict(
            :solution => has_solution(s) ? collect(best_values(s)) : nothing,
            :time => time_info(s),
            :type => sat ? "Satisfaction" : "Optimization",
        )
        !sat && has_solution(s) && push!(info, :value => best_value(s))
        write(path, JSON.json(info))
    end
end


text:
using ProjAssistant
@quickactivate

@time begin
    import Chemostat_EColi
    const ChE = Chemostat_EColi

    import Chemostat_Kayser2005
    const ChK = Chemostat_Kayser2005

    import Chemostat_Nanchen2006
    const ChN = Chemostat_Nanchen2006

    import Chemostat_Folsom2014
    const ChF = Chemostat_Folsom2014

    import Chemostat_InSilico
    const ChIn = Chemostat_InSilico

    using Base.Threads
    using Dates
end

## ------------------------------------------------------
# utils
function _info(msg; kwargs...)
    println()
    msg = rpad(string(msg, " "), 60, "-")
    time = now()
    @info(msg, time, kwargs...)
    println()
    println()
end

function run_script(Proj::Module, arg, args...; flags::Vector = [], nthrs = NTHREADS)
    src = scriptsdir(Proj, arg, args...)
    !isfile(src) && error("script file is missing, file: ", src)

    _info("Starting"; src, nthrs)
    cmd = Cmd(string.([JULIA_CMD, "-t$(nthrs)", "--startup-file=no", "--project", src, flags...]))
    run(cmd; wait = true)
    _info("Finished"; src, nthrs)

    return nothing
end

## ------------------------------------------------------
# ARGS
using ArgParse

set = ArgParseSettings()
@add_arg_table! set begin
    "--julia-cmd"
        help = "the juila command to use"
        default = Base.julia_exename()
        arg_type=String
    "--pdflatex-cmd"
        help = "the pdflatex command to use"
        default = "pdflatex"
        arg_type=String
    "--bibtex-cmd"
        help = "the bibtex command to use"
        default = "bibtex"
        arg_type=String
    "--nthrs", "-t"
        help = "the juila command to use"
        default = max(Sys.CPU_THREADS - 1, 1)
        arg_type=Int
end

if isinteractive()
    # Dev values
    JULIA_CMD = Base.julia_exename()
    NTHREADS = max(Sys.CPU_THREADS - 1, 1)
    PDFLATEX_CMD = "pdflatex"
    BIBTEX_CMD = "bibtex"
else
    parsed_args = parse_args(set)
    JULIA_CMD = parsed_args["julia-cmd"]
    NTHREADS = parsed_args["nthrs"]
    PDFLATEX_CMD = parsed_args["pdflatex-cmd"]
    BIBTEX_CMD = parsed_args["bibtex-cmd"]
end

PROJ_ROOT = dirname(projdir(ChE))

_info("globals"; PROJ_ROOT, NTHREADS, JULIA_CMD, PDFLATEX_CMD, BIBTEX_CMD)

## ------------------------------------------------------
# testing juila
try; run_script(ChE, "0.0_test.jl")
    catch err
    @error("Problem calling julia", JULIA_CMD)
    rethrow(err)
end

# check proj dir
let
    paths = readdir(PROJ_ROOT)
    for proj in [
            # "MaxEnt_EColi_paper", # TODO
            "Chemostat_EColi",
            "Chemostat_Kayser2005",
            "Chemostat_Nanchen2006",
            "Chemostat_Folsom2014",
        ]
        !(proj in paths) && error("proj ", proj, " is missing from the root dir")
    end
    _info("Root dir test passed")
end

## ------------------------------------------------------
# run simulation
let
    # sim2D
    _info("Run simulation")
    run_script(ChIn.Dynamic, "1.0_compute_Spaces.jl")
    run_script(ChIn.Dynamic, ["sim3D"], "2.0_run_sim3D.jl")
    run_script(ChIn.Dynamic, ["sim3D"], "3.0_run_MaxEnt.jl")
end

## ------------------------------------------------------
# run EColi packages
let
    _info("Run EColi packages")

    ## ------------------------------------------------------
    # Kayser2005
    run_script(ChK.KayserData, "1_converted_data.jl")
    run_script(ChK.iJR904, "1.0_prepare_gem.jl")
    run_script(ChK.iJR904, "2.0_maxent_ep_variants.jl")

    ## ------------------------------------------------------
    # Nanchen2005
    run_script(ChN.iJR904, "1.0_prepare_gem.jl")
    run_script(ChN.iJR904, "2.0_maxent_ep_variants.jl")

    ## ------------------------------------------------------
    # Folsom2014
    run_script(ChF.iJR904, "1.0_prepare_gem.jl")
    run_script(ChF.iJR904, "2.0_maxent_ep_variants.jl")

end

## ------------------------------------------------------
# collect
let
    run_script(ChE, ["experiments"], "0.1_LP.jl")
    run_script(ChE, ["experiments"], "0.2_collect_DAT.jl")
end

## ------------------------------------------------------
# make plots
let
    run_script(ChE, ["simulation"], "2.0_sim_plots.jl")
    run_script(ChE, ["experiments"], "1.0_exp_plots.jl")
end

# ## ------------------------------------------------------
# # run latex
# let
#     tex_proj = joinpath(PROJ_ROOT, "MaxEnt_EColi_paper")
#     teximgs_dir = joinpath(tex_proj, "images")
#     !isdir(teximgs_dir) && error("image dir not found, path: ", imgs_dir)

#     # cp to image dir
#     _info("Copying figures"; teximgs_dir)
#     for figi in 1:10
#         imgname = string("figure_", figi, ".png")
#         srcfile = plotsdir(ChE, imgname)
#         !isfile(srcfile) && error("image file is missing, file: ", srcfile)
#         destfile = joinpath(teximgs_dir, imgname)
#         cp(srcfile, destfile; force = true)
#         @info("Copied figure file", srcfile, destfile)
#     end

#     # try to run pdflatex
#     _info("Trying pdflatex"; tex_proj)
#     tex_mainfile = joinpath(tex_proj, "MaxEnt_EColi.tex")
#     !isfile(tex_mainfile) && error("main tex file not found, file: ", tex_mainfile)
#     pdflatex_cmd = Cmd(
#         [PDFLATEX_CMD, "-synctex=1", "-interaction=nonstopmode", "-file-line-error", "-recorder",
#         string("-output-directory=", tex_proj), tex_mainfile]
#     )
#     bibtex_cmd = Cmd(`$(BIBTEX_CMD)  "MaxEnt_EColi"`)

#     try;
#         cd(tex_proj)

#         # from vscode recipe
#         run(pdflatex_cmd; wait = true)
#         run(bibtex_cmd; wait = true)
#         run(pdflatex_cmd; wait = true)
#         run(bibtex_cmd; wait = true)
#         run(pdflatex_cmd; wait = true)

#         catch err
#             @info("pdflatex/bibtex failed, but the tex project is updated. You must compile it manually. ", tex_proj)
#             rethrow(err)
#         finally; cd(projdir(ChE))
#     end

# end

_info("The end")
exit()


text:
@testset "Scale" begin
    @test typeof(@inferred(Scale(1))) <: Scale <: Augmentor.AffineOperation
    @testset "constructor" begin
        @test_throws MethodError Scale()
        @test_throws MethodError Scale(())
        @test_throws MethodError Scale{0}(())
        @test_throws MethodError Scale(:a)
        @test_throws MethodError Scale([:a])
        @test_throws ArgumentError Scale([])
        @test_throws ArgumentError Scale{1}(([:a],))
        @test_throws ArgumentError Scale(1:3,1:2)
        @test_throws ArgumentError Scale(1:2,1:2,1)
        @test @inferred(Scale(1,2)) == @inferred(Scale{2}((1,2)))
        @test @inferred(Scale(1)) == @inferred(Scale{2}((1:1,1:1)))
        @test @inferred(Scale((1,))) == Scale{1}((1:1,))
        @test @inferred(Scale(1,2)) == Scale{2}((1:1,2:2))
        @test @inferred(Scale(1.)) == Scale{2}((1.0:1:1,1.0:1:1))
        @test @inferred(Scale(1,2,3.)) == Scale{3}((1.0:1,2.0:2,3.0:3))
        @test @inferred(Scale(1,2,3)) == Scale{3}((1:1,2:2,3:3))
        @test @inferred(Scale(1:1,2,3)) == Scale{3}((1:1,2:2,3:3))
        @test @inferred(Scale((1,2,3))) == Scale{3}((1:1,2:2,3:3))
        @test @inferred(Scale(1:2,2:3)) == Scale{2}((1:2,2:3))
        @test @inferred(Scale((1:2,2:3))) == Scale{2}((1:2,2:3))
        op = @inferred(Scale([1,2],2:3))
        @test op.factors == ([1.,2.],[2.,3.])
        @test str_show(op) == "Augmentor.Scale{2}(([1.0, 2.0], [2.0, 3.0]))"
        @test str_showconst(op) == "Scale([1.0, 2.0], [2.0, 3.0])"
        @test str_showcompact(op) == "Scale by I ∈ {1.0×2.0, 2.0×3.0}"
        op = @inferred(Scale([1,2],[2.,3]))
        @test op.factors == ([1.,2.],[2.,3.])
        op = @inferred(Scale([1,2]))
        @test op.factors == ([1.,2.],[1.,2.])
        op = @inferred(Scale(1))
        @test str_show(op) == "Augmentor.Scale{2}((1, 1))"
        @test str_showconst(op) == "Scale(1, 1)"
        @test str_showcompact(op) == "Scale by 1×1"
        op = @inferred(Scale(0.8,0.9))
        @test str_show(op) == "Augmentor.Scale{2}((0.8, 0.9))"
        @test str_showconst(op) == "Scale(0.8, 0.9)"
        @test str_showcompact(op) == "Scale by 0.8×0.9"
        op = @inferred(Scale(1:2,3:4))
        @test str_show(op) == "Augmentor.Scale{2}((1:2, 3:4))"
        @test str_showconst(op) == "Scale(1:2, 3:4)"
        @test str_showcompact(op) == "Scale by I ∈ {1×3, 2×4}"
        op = @inferred(Scale(1:3,3:5))
        @test str_show(op) == "Augmentor.Scale{2}((1:3, 3:5))"
        @test str_showconst(op) == "Scale(1:3, 3:5)"
        @test str_showcompact(op) == "Scale by I ∈ {1×3, 2×4, 3×5}"
        op = @inferred(Scale(1,2,3))
        @test str_show(op) == "Augmentor.Scale{3}((1, 2, 3))"
        @test str_showconst(op) == "Scale(1, 2, 3)"
        @test str_showcompact(op) == "Scale by 1×2×3"
    end
    @testset "randparam" begin
        @test @inferred(Augmentor.randparam(Scale(45), rect)) === (45.0, 45.0)
        @test @inferred(Augmentor.randparam(Scale(20.0, 15.0), rect)) === (20.0, 15.0)
        @test @inferred(Augmentor.randparam(Scale(20.0, 15.0), (rect, square2))) === (20.0, 15.0)
        @test @inferred(Augmentor.randparam(Scale(1:2,2:3), rect)) in [(1,2), (2,3)]
        @test @inferred(Augmentor.toaffinemap(Scale(2,3), rect, (2.,3.))) ≈ AffineMap([2. 0.; 0. 3.], [-1.5,-4.0])
        @test @inferred(Augmentor.toaffinemap(Scale([0.9,0.4],[0.8,0.3]), rect, (0.9,0.8))) ≈ AffineMap([.9 0.; 0. .8], [0.15,0.4])
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(Scale(10), nothing)
        @test Augmentor.supports_eager(Scale) === false
        # TODO: actual content tests (maybe test_reference)
        img_out1 = @inferred Augmentor.applyeager(Scale(1.5), square2)
        img_out2 = @inferred Augmentor.applyeager(Scale(0.2), square2)
        @test indices(img_out1) == (0:5, 0:5)
        @test indices(img_out2) == (2:3, 2:3)
        imgs = [
            (square2),
            (view(square2, :, :)),
            (Augmentor.prepareaffine(square2)),
            (OffsetArray(square2, -1, -1)),
            (view(square2, IdentityRange(1:4), IdentityRange(1:4))),
        ]
        @testset "fixed parameter" begin
            for img_in in imgs
                res = @inferred(Augmentor.applyeager(Scale(1.5), img_in))
                @test parent(res) ≈ parent(img_out1)
                @test typeof(res) == typeof(img_out1)
                res = @inferred(Augmentor.applyeager(Scale(0.2), img_in))
                @test parent(res) == parent(img_out2)
                @test typeof(res) == typeof(img_out2)
                # test same with tuple of images
                res1, res2 = @inferred(Augmentor.applyeager(Scale(1.5), (img_in, N0f8.(img_in))))
                @test parent(res1) ≈ parent(img_out1)
                @test parent(res2) == parent(img_out1)
                @test typeof(res1) == typeof(img_out1)
                @test typeof(res2) <: OffsetArray{N0f8}
                res1, res2 = @inferred(Augmentor.applyeager(Scale(0.2), (img_in, N0f8.(img_in))))
                @test parent(res1) ≈ parent(img_out2)
                @test parent(res2) == parent(img_out2)
                @test typeof(res1) == typeof(img_out2)
                @test typeof(res2) <: OffsetArray{N0f8}
            end
            # check that the affine map is computed for each image
            res1, res2 = @inferred(Augmentor.applyeager(Scale(1.5), (square, OffsetArray(square,-5,-5))))
            @test collect(res1) == collect(res2)
            @test indices(res1) != indices(res2)
        end
        @testset "random parameter" begin
            for img_in in imgs
                res1, res2 = @inferred(Augmentor.applyeager(Scale(.5:.1:2., .6:.1:2.1), (img_in, N0f8.(img_in))))
                # make sure same scales are used
                @test res1 == res2
                @test typeof(res1) == typeof(img_out1)
                @test typeof(res2) <: OffsetArray{N0f8}
            end
        end
    end
    ref = @inferred Augmentor.applyeager(Scale(2,3), square)
    @testset "affine" begin
        @test Augmentor.supports_affine(Scale) === true
        @test_throws MethodError Augmentor.applyaffine(Scale(90), nothing)
        wv = @inferred Augmentor.applyaffine(Scale(2,3), Augmentor.prepareaffine(square))
        @test wv == ref
        # TODO: better tests
        @test parent(wv).itp.coefs === square
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(Scale) === true
        @test_throws MethodError Augmentor.applyaffineview(Scale(90), nothing)
        wv = @inferred Augmentor.applyaffineview(Scale(2,3), Augmentor.prepareaffine(square))
        @test wv == ref
        # TODO: better tests
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(Scale) === true
        wv = @inferred Augmentor.applylazy(Scale(2,3), square)
        @test wv == ref
        # TODO: better tests
        @test parent(wv).itp.coefs === square
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "view" begin
        @test Augmentor.supports_view(Scale) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(Scale) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(Scale) === false
    end
end


text:
# using LinearAlgebra

"""
randomized SVD.
"""
function rsvd(A::AbstractArray{T}, k::Int=min(size(A)...), oversample::Int=10, power::Int=10, ortho::Bool=false) where T
    m, n = size(A)
    p = min(n,oversample*k)
    Y = A * randn(T, n,p)

    for i = 1:power
        Y = A * (ortho ? LinearAlgebra.qr(A' * Y).Q : A' * Y)
    end
    Q, R = LinearAlgebra.qr(Y)
    B = Q' * A

    U, s ,V = LinearAlgebra.svd(B)
    U = Q * U
    k = min(k,size(U, 2))
    U = U[:,1:k]
    V = V[:,1:k]
    s = s[1:k]
    return U,s,V
end


text:
## Cootes images (from ICAAM project)

immutable CootesDataset <: FaceDataset
end

const COOTES_DATA_DIR = joinpath(Pkg.dir(), "FaceDatasets", "data", "cootes", "data")
const COOTES_IMG_HEIGHT = 480

function load_shape_from_mat(path::AbstractString)
    return matread(path)["annotations"]
end

function load_shapes(::Type{CootesDataset}; count=-1)
    files = sort(filter(x -> endswith(x, ".mat"), readdir(COOTES_DATA_DIR)))
    paths = map(x->joinpath(COOTES_DATA_DIR, x), files)
    n_use = count > 0 ? count : length(paths)
    shapes = Array(Matrix{Float64}, n_use)
    for k=1:n_use
        shape_xy = load_shape_from_mat(paths[k])
        shape_ij = [COOTES_IMG_HEIGHT .- shape_xy[:, 2] shape_xy[:, 1]]
        shapes[k] = shape_ij
    end
    return shapes
end

function load_images(::Type{CootesDataset}; count=-1)
    files = sort(filter(x->endswith(x, ".bmp"), readdir(COOTES_DATA_DIR)))
    paths = map(x->joinpath(COOTES_DATA_DIR, x), files)
    n_use = count > 0 ? count : length(paths)
    imgs = Array(Array{Float64, 3}, n_use)
    for i=1:n_use
        img_rgb = load(paths[i])
        img_rgb = convert(Image{RGB{Float64}}, img_rgb)  # ensure 3-channel RGB
        imgs[i] = convert(Array, separate(img_rgb))
    end
    return imgs
end


text:
@everywhere begin
    using StochasticSearch
    function rosenbrock(x::Configuration, parameters::Dict{Symbol, Any})
        return (1.0 - x["i0"].value)^2 + 100.0 * (x["i1"].value - x["i0"].value^2)^2
    end
end

configuration = Configuration([FloatParameter(-2.0, 2.0, 0.0,"i0"),
                               FloatParameter(-2.0, 2.0, 0.0,"i1")],
                               "rosenbrock_config")

tuning_run = Run(cost               = rosenbrock,
                 starting_point     = configuration,
                 methods            = [[:simulated_annealing 1];
                                       [:iterative_first_improvement 1];
                                       [:randomized_first_improvement 1];
                                       [:iterative_greedy_construction 1];
                                       [:iterative_probabilistic_improvement 1];])

search_task = @task optimize(tuning_run)
result = consume(search_task)
print(result)
while result.is_final == false
    result = consume(search_task)
    print(result)
end


text:
import KnuthBendix: KnuthBendix.BufferWord

@testset "KnuthBendix.BufferWord internals" begin
    @test one(KnuthBendix.BufferWord{Int}) isa KnuthBendix.BufferWord{Int}
    @test isone(one(KnuthBendix.BufferWord{UInt16}))
    @test KnuthBendix.BufferWord([1,2,3]) isa KnuthBendix.BufferWord{UInt16}
    @test !isone(KnuthBendix.BufferWord([1,2,3]))

    let bw = one(KnuthBendix.BufferWord{UInt16});
        append!(bw, [1,2,3,4])
        il = KnuthBendix.internal_length(bw)
        W = Word(bw)
        KnuthBendix._growatend!(bw, 3)
        @test KnuthBendix.internal_length(bw) == il + 3
        @test bw == W
        KnuthBendix._growatbeg!(bw, 3)
        @test bw == W
        @test KnuthBendix.internal_length(bw) == il + 3 + 3

        @test pop!(bw) == 4
        @test bw == [1,2,3]
        @test popfirst!(bw) == 1
        @test bw == [2,3]

        @test KnuthBendix.internal_length(bw) == il + 3 + 3
    end

    let bw = KnuthBendix.BufferWord{UInt16}([1,2,3,4], 0, 0)
        @test KnuthBendix.internal_length(bw) == 4
        W = Word(bw)
        KnuthBendix._growatend!(bw, 4);
        @test bw == W
        @test KnuthBendix.internal_length(bw) == 8
        KnuthBendix._growatbeg!(bw, 4);
        @test bw == W
        @test KnuthBendix.internal_length(bw) == 12
        @test Word(bw) == W
    end
end

@testset "KnuthBendix.BufferWord push/append" begin

    let bw_orig = one(KnuthBendix.BufferWord{UInt16})

        bw = deepcopy(bw_orig)
        @test push!(bw,1) == [1]

        bw = deepcopy(bw_orig)
        @test pushfirst!(bw,1) == [1]

        bw = deepcopy(bw_orig)
        @test prepend!(bw,[1]) == [1]

        bw = deepcopy(bw_orig)
        @test append!(bw,[1]) == [1]
    end

    let bw_orig = KnuthBendix.BufferWord([2])

        bw = deepcopy(bw_orig)
        @test push!(bw,1) == [2,1]

        bw = deepcopy(bw_orig)
        @test pushfirst!(bw,1) == [1,2]

        bw = deepcopy(bw_orig)
        @test prepend!(bw,[1]) == [1, 2]

        bw = deepcopy(bw_orig)
        @test append!(bw,[1]) == [2,1]
    end

    let bw_orig = KnuthBendix.BufferWord{UInt16}([2], 0, 0)

        bw = deepcopy(bw_orig)
        @test push!(bw,1) == [2,1]

        bw = deepcopy(bw_orig)
        @test pushfirst!(bw,1) == [1,2]

        bw = deepcopy(bw_orig)
        @test prepend!(bw,[1]) == [1, 2]

        bw = deepcopy(bw_orig)
        @test append!(bw,[1]) == [2,1]
    end

    let bw_orig = KnuthBendix.BufferWord{UInt16}([2,2,2,2], 0, 0)

        bw = deepcopy(bw_orig)
        @test push!(bw,1) == [2,2,2,2,1]

        bw = deepcopy(bw_orig)
        @test pushfirst!(bw,1) == [1,2,2,2,2]

        bw = deepcopy(bw_orig)
        @test prepend!(bw,[1,1,1,1]) == [1,1,1,1,2,2,2,2]

        bw = deepcopy(bw_orig)
        @test append!(bw,[1,1,1,1]) == [2,2,2,2,1,1,1,1]
    end
end

abstract_word_conformance_test(KnuthBendix.BufferWord)


text:
using Test, LibCEED, LinearAlgebra, StaticArrays

showstr(x) = sprint(show, MIME("text/plain"), x)
summarystr(x) = sprint(summary, x)
getoutput(fname) =
    chomp(read(joinpath(@__DIR__, "output", string(CeedScalar), fname), String))

function checkoutput(str, fname)
    if str != getoutput(fname)
        write(fname, str)
        return false
    end
    return true
end

mutable struct CtxData
    io::IOBuffer
    x::Vector{Float64}
end

const run_dev_tests = !isrelease() || ("--run-dev-tests" in ARGS)

if run_dev_tests
    include("rundevtests.jl")
end

if !LibCEED.ceedversion_ge(LibCEED.minimum_libceed_version) && !run_dev_tests
    @warn "Skipping tests because of incompatible libCEED versions."
else
    @testset "LibCEED Release Tests" begin
        @testset "LibCEED" begin
            @test ceedversion() isa VersionNumber
            @test isrelease() isa Bool
            @test isfile(get_libceed_path())
        end

        @testset "Ceed" begin
            res = "/cpu/self/ref/serial"
            c = Ceed(res)
            @test isdeterministic(c)
            @test getresource(c) == res
            @test !iscuda(c)
            @test get_preferred_memtype(c) == MEM_HOST
            @test_throws LibCEED.CeedError create_interior_qfunction(c, "")
            @test showstr(c) == """
                Ceed
                  Ceed Resource: $res
                  Preferred MemType: host"""
        end

        @testset "Context" begin
            c = Ceed()
            data = zeros(CeedScalar, 3)
            ctx = Context(c, data)
            @test showstr(ctx) == """
                CeedQFunctionContext
                  Context Data Size: $(sizeof(data))"""
            @test_throws Exception set_data!(ctx, MEM_HOST, OWN_POINTER, data)
        end

        @testset "CeedVector" begin
            n = 10
            c = Ceed()
            v = CeedVector(c, n)
            @test size(v) == (n,)
            @test length(v) == n
            @test axes(v) == (1:n,)
            @test ndims(v) == 1
            @test ndims(CeedVector) == 1

            v[] = 0.0
            @test @witharray(a = v, all(a .== 0.0))

            v1 = rand(CeedScalar, n)
            v2 = CeedVector(c, v1)
            @test @witharray_read(a = v2, mtype = MEM_HOST, a == v1)
            @test Vector(v2) == v1
            v[] = v1
            for p ∈ [1, 2, Inf]
                @test norm(v, p) ≈ norm(v1, p)
            end
            @test_throws Exception norm(v, 3)
            @test witharray_read(sum, v) == sum(v1)
            reciprocal!(v)
            @test @witharray(a = v, mtype = MEM_HOST, all(a .== CeedScalar(1.0)./v1))

            witharray(x -> x .= 1.0, v)
            @test @witharray(a = v, all(a .== 1.0))

            @test summarystr(v) == "$n-element CeedVector"
            @test sprint(show, v) == @witharray_read(a = v, sprint(show, a))
            io = IOBuffer()
            summary(io, v)
            println(io, ":")
            @witharray_read(a = v, Base.print_array(io, a))
            s1 = String(take!(io))
            @test showstr(v) == s1

            setarray!(v, MEM_HOST, USE_POINTER, v1)
            syncarray!(v, MEM_HOST)
            @test @witharray_read(a = v, a == v1)
            p = takearray!(v, MEM_HOST)
            @test p == pointer(v1)

            m = rand(CeedScalar, 10, 10)
            vm = CeedVector(c, vec(m))
            @test @witharray_read(a = vm, size = size(m), a == m)

            @test CeedVectorActive()[] == LibCEED.C.CEED_VECTOR_ACTIVE[]
            @test CeedVectorNone()[] == LibCEED.C.CEED_VECTOR_NONE[]

            w1 = rand(CeedScalar, n)
            w2 = rand(CeedScalar, n)
            w3 = rand(CeedScalar, n)

            cv1 = CeedVector(c, w1)
            cv2 = CeedVector(c, w2)
            cv3 = CeedVector(c, w3)

            alpha = rand(CeedScalar)

            scale!(cv1, alpha)
            w1 .*= alpha
            @test @witharray_read(a = cv1, a == w1)

            pointwisemult!(cv1, cv2, cv3)
            w1 .= w2.*w3
            @test @witharray_read(a = cv1, a == w1)

            axpy!(alpha, cv2, cv1)
            axpy!(alpha, w2, w1)
            @test @witharray_read(a = cv1, a ≈ w1)
        end

        @testset "Basis" begin
            c = Ceed()
            dim = 3
            ncomp = 1
            p = 4
            q = 6
            b1 = create_tensor_h1_lagrange_basis(c, dim, ncomp, p, q, GAUSS_LOBATTO)

            @test checkoutput(showstr(b1), "b1.out")
            @test getdimension(b1) == 3
            @test gettopology(b1) == HEX
            @test getnumcomponents(b1) == ncomp
            @test getnumnodes(b1) == p^dim
            @test getnumnodes1d(b1) == p
            @test getnumqpts(b1) == q^dim
            @test getnumqpts1d(b1) == q

            q1d, w1d = lobatto_quadrature(3, AbscissaAndWeights)
            @test q1d ≈ CeedScalar[-1.0, 0.0, 1.0]
            @test w1d ≈ CeedScalar[1/3, 4/3, 1/3]

            q1d, w1d = gauss_quadrature(3)
            @test q1d ≈ CeedScalar[-sqrt(3/5), 0.0, sqrt(3/5)]
            @test w1d ≈ CeedScalar[5/9, 8/9, 5/9]

            b1d = CeedScalar[1.0 0.0; 0.5 0.5; 0.0 1.0]
            d1d = CeedScalar[-0.5 0.5; -0.5 0.5; -0.5 0.5]
            q1d = CeedScalar[-1.0, 0.0, 1.0]
            w1d = CeedScalar[1/3, 4/3, 1/3]
            q, p = size(b1d)
            d2d = zeros(CeedScalar, 2, q*q, p*p)
            d2d[1, :, :] = kron(b1d, d1d)
            d2d[2, :, :] = kron(d1d, b1d)

            dim2 = 2
            b2 = create_tensor_h1_basis(c, dim2, 1, p, q, b1d, d1d, q1d, w1d)
            @test getinterp(b2) == kron(b1d, b1d)
            @test getinterp1d(b2) == b1d
            @test getgrad(b2) == d2d
            @test getgrad1d(b2) == d1d
            @test checkoutput(showstr(b2), "b2.out")

            b3 = create_h1_basis(c, LINE, 1, p, q, b1d, reshape(d1d, 1, q, p), q1d, w1d)
            @test getqref(b3) == q1d
            @test getqweights(b3) == w1d
            @test checkoutput(showstr(b3), "b3.out")

            v = rand(CeedScalar, 2)
            vq = apply(b3, v)
            vd = apply(b3, v; emode=EVAL_GRAD)
            @test vq ≈ b1d*v
            @test vd ≈ d1d*v

            @test BasisCollocated()[] == LibCEED.C.CEED_BASIS_COLLOCATED[]
        end

        @testset "Request" begin
            @test RequestImmediate()[] == LibCEED.C.CEED_REQUEST_IMMEDIATE[]
            @test RequestOrdered()[] == LibCEED.C.CEED_REQUEST_ORDERED[]
        end

        @testset "Misc" begin
            for dim = 1:3
                D = CeedDim(dim)
                J = rand(CeedScalar, dim, dim)
                @test det(J, D) ≈ det(J)
                J = J + J' # make symmetric
                @test setvoigt(SMatrix{dim,dim}(J)) == setvoigt(J, D)
                @test getvoigt(setvoigt(J, D)) == J
                V = zeros(CeedScalar, dim*(dim + 1)÷2)
                setvoigt!(V, J, D)
                @test V == setvoigt(J, D)
                J2 = zeros(CeedScalar, dim, dim)
                getvoigt!(J2, V, D)
                @test J2 == J
            end
        end

        @testset "QFunction" begin
            c = Ceed()

            id = create_identity_qfunction(c, 1, EVAL_INTERP, EVAL_INTERP)
            Q = 10
            v = rand(CeedScalar, Q)
            v1 = CeedVector(c, v)
            v2 = CeedVector(c, Q)
            apply!(id, Q, [v1], [v2])
            @test @witharray(a = v2, a == v)

            @test showstr(create_interior_qfunction(c, "Poisson3DApply")) == """
                Gallery CeedQFunction Poisson3DApply
                  2 Input Fields:
                    Input Field [0]:
                      Name: "du"
                      Size: 3
                      EvalMode: "gradient"
                    Input Field [1]:
                      Name: "qdata"
                      Size: 6
                      EvalMode: "none"
                  1 Output Field:
                    Output Field [0]:
                      Name: "dv"
                      Size: 3
                      EvalMode: "gradient\""""

            @interior_qf id2 = (c, (a, :in, EVAL_INTERP), (b, :out, EVAL_INTERP), b .= a)
            v2[] = 0.0
            apply!(id2, Q, [v1], [v2])
            @test @witharray(a = v2, a == v)

            ctxdata = CtxData(IOBuffer(), rand(CeedScalar, 3))
            ctx = Context(c, ctxdata)
            dim = 3
            @interior_qf qf = (
                c,
                dim=dim,
                ctxdata::CtxData,
                (a, :in, EVAL_GRAD, dim),
                (b, :in, EVAL_NONE),
                (c, :out, EVAL_INTERP),
                begin
                    c[] = b*sum(a)
                    show(ctxdata.io, MIME("text/plain"), ctxdata.x)
                end,
            )
            set_context!(qf, ctx)
            in_sz, out_sz = LibCEED.get_field_sizes(qf)
            @test in_sz == [dim, 1]
            @test out_sz == [1]
            v1 = rand(CeedScalar, dim)
            v2 = rand(CeedScalar, 1)
            cv1 = CeedVector(c, v1)
            cv2 = CeedVector(c, v2)
            cv3 = CeedVector(c, 1)
            apply!(qf, 1, [cv1, cv2], [cv3])
            @test String(take!(ctxdata.io)) == showstr(ctxdata.x)
            @test @witharray_read(v3 = cv3, v3[1] == v2[1]*sum(v1))

            @test QFunctionNone()[] == LibCEED.C.CEED_QFUNCTION_NONE[]
        end

        @testset "Operator" begin
            c = Ceed()
            @interior_qf id = (
                c,
                (input, :in, EVAL_INTERP),
                (output, :out, EVAL_INTERP),
                begin
                    output[] = input
                end,
            )
            b = create_tensor_h1_lagrange_basis(c, 3, 1, 3, 3, GAUSS_LOBATTO)
            n = getnumnodes(b)
            offsets = Vector{CeedInt}(0:n-1)
            r = create_elem_restriction(c, 1, n, 1, 1, n, offsets)
            op = Operator(
                c;
                qf=id,
                fields=[
                    (:input, r, b, CeedVectorActive()),
                    (:output, r, b, CeedVectorActive()),
                ],
            )
            @test showstr(op) == """
                CeedOperator
                  2 Fields
                  1 Input Field:
                    Input Field [0]:
                      Name: "input"
                      Active vector
                  1 Output Field:
                    Output Field [0]:
                      Name: "output"
                      Active vector"""

            v = rand(CeedScalar, n)
            v1 = CeedVector(c, v)
            v2 = CeedVector(c, n)
            apply!(op, v1, v2)
            @test @witharray_read(a1 = v1, @witharray_read(a2 = v2, a1 == a2))
            apply_add!(op, v1, v2)
            @test @witharray_read(a1 = v1, @witharray_read(a2 = v2, a1 + a1 == a2))

            diag_vector = create_lvector(r)
            LibCEED.assemble_diagonal!(op, diag_vector)
            @test @witharray_read(a = diag_vector, a == ones(n))
            # TODO: change this test after bug-fix in libCEED
            diag_vector[] = 0.0
            LibCEED.assemble_add_diagonal!(op, diag_vector)
            @test @witharray(a = diag_vector, a == fill(1.0, n))

            comp_op = create_composite_operator(c, [op])
            apply!(comp_op, v1, v2)
            @test @witharray_read(a1 = v1, @witharray_read(a2 = v2, a1 == a2))
        end

        @testset "ElemRestriction" begin
            c = Ceed()
            n = 10
            offsets = Vector{CeedInt}([0:n-1; n-1:2*n-2])
            lsize = 2*n - 1
            r = create_elem_restriction(c, 2, n, 1, lsize, lsize, offsets)
            @test getcompstride(r) == lsize
            @test getnumelements(r) == 2
            @test getelementsize(r) == n
            @test getlvectorsize(r) == lsize
            @test getnumcomponents(r) == 1
            @test length(create_lvector(r)) == lsize
            @test length(create_evector(r)) == 2*n
            lv, ev = create_vectors(r)
            @test length(lv) == lsize
            @test length(ev) == 2*n
            mult = getmultiplicity(r)
            mult2 = ones(lsize)
            mult2[n] = 2
            @test mult == mult2
            rand_lv = rand(CeedScalar, lsize)
            rand_ev = [rand_lv[1:n]; rand_lv[n:end]]
            @test apply(r, rand_lv) == rand_ev
            @test apply(r, rand_ev; tmode=TRANSPOSE) == rand_lv.*mult
            @test showstr(r) == string(
                "CeedElemRestriction from (19, 1) to 2 elements ",
                "with 10 nodes each and component stride 19",
            )

            strides = CeedInt[1, n, n]
            rs = create_elem_restriction_strided(c, 1, n, 1, n, strides)
            @test showstr(rs) == string(
                "CeedElemRestriction from (10, 1) to 1 elements ",
                "with 10 nodes each and strides [1, $n, $n]",
            )

            @test ElemRestrictionNone()[] == LibCEED.C.CEED_ELEMRESTRICTION_NONE[]
        end
    end
end

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
using Test

using LinearAlgebra

include("../../lib/KTBC.jl")
import .KTBC.CreateKTBC

# 10
KTBC₃ = CreateKTBC(3)
K₃ = convert(Array{Int}, KTBC₃[1])

res = 4 * K₃

invRes = inv(res)

u(x) = x - x^2

@test u.([
  1/4
  1/2
  3/4
]) == round.(invRes * (1 / 2) * [
  1
  1
  1
], digits=4)

# 11
K = [
   12 -3
  -3   12
]

inv(K) * [
  1/3
  2/3
]

f(x) = - (x^3)/6 + x/6

[
  f(1/3)
  f(2/3)
]


text:
export SpHTypes
"""
    const SpHTypes = Union{Float64,Complex{Float64}}

The types supported by FastSphericalHarmonics (the same types as for
FastTransforms).
"""
const SpHTypes = Union{Float64,Complex{Float64}}

export sph_points
"""
    Θ, Φ = sph_points(N::Integer)
    Θ::Vector{Float64}
    Φ::Vector{Float64}

Calculate the locations of points on the sphere when using `N` points
in the θ (latitudinal) direction.

It is `length(Θ) = N` and `length(Φ) = M` where `M = 2N-1`.
"""
function sph_points(N::Int)
    @assert N > 0
    M = 2 * N - 1
    return π / N * (0.5:(N - 0.5)), 2π / M * (0:(M - 1))
end
sph_points(N::Integer) = sph_points(Int(N))

export sph_lmax
"""
    lmax = sph_lmax(N::Integer)
    lmax::Int

Calculate the maximum `l` mode that can be represented with `N`
points. It is `lmax = N - 1`.
"""
function sph_lmax(N::Int)
    @assert N > 0
    return N - 1
end
sph_lmax(N::Integer) = sph_lmax(Int(N))


text:
#!/usr/bin/env julia

using Gtk
using GtkBuilderAid

example_app = GtkApplication("io.github.matt5sean3.GtkBuilderAid.first", 0)

ft = @GtkFunctionTable begin

@guarded function click_ok(
    widget,
    qsdata)
  println("OK clicked!")
  return nothing
end

@guarded function close_window(
    widget,
    window)
  destroy(window)
  return nothing
end

end

start_application(ft, "io.github.matt5sean3.GtkBuilderAid.first", "main_window", "resources/main.ui", nothing)


text:
"""
    orient_strands(seqs::Array{String}, phreds::Union{Array{Vector{Phred},1},Void},
                   names::Union{Array{String,1},Void}, ref::String; k::Int=6)

Orients sequences (with phreds and names, which may be `nothing`) relative to a reference sequence.
`k` is kmer size for computing kmer vectors.
"""
function orient_strands(seqs::Array{String}, phreds::Union{Array{Vector{Phred},1},Void}, names::Union{Array{String,1},Void}, ref::String; k::Int=6)
    kmers = [kmer_count(s, k) for s in seqs]
    ref_kmer = kmer_count(ref, k)
    rev_ref_kmer = kmer_count(reverse_complement(ref), k)
    #return kmer_orient_strands(seqs, ref, kmers, ref_kmer, rev_ref_kmer)

    # if this method crashes use the other return statement, it shouldnt work but it might
    return orient_strands(seqs, phreds, names, ref, kmers, ref_kmer, rev_ref_kmer)
end

"""
    orient_strands(seqs::Array{String}, ref::String; k::Int=6)

Orients sequences relative to a reference sequence.
`k` is kmer size for computing kmer vectors.
"""
orient_strands(seqs::Array{String}, ref::String; k::Int=6) = orient_strands(seqs, nothing, nothing, ref, k=k)[1]

"""
    orient_strands(seqs::Array{String,1}, phreds::Union{Array{Vector{Phred},1},Void}, names::Union{Array{String,1},Void},
                        ref::String, kmers::Array{Array{T,1},1}, ref_kmer::Array{T,1}, rev_ref_kmer::Array{T,1}) where {T <: Real}

Orients sequences with given kmer vectors relative to a reference sequence with given kmer vector.
"""
function orient_strands(seqs::Array{String,1}, phreds::Union{Array{Vector{Phred},1},Void}, names::Union{Array{String,1},Void},
                        ref::String, kmers::Array{Array{T,1},1}, ref_kmer::Array{T,1}, rev_ref_kmer::Array{T,1}) where {T <: Real}
    oriented = String[]
    newphreds = phreds==nothing ? nothing : Vector{Phred}[]
    newnames = names==nothing ? nothing : String[]

    for i in 1:length(seqs)
        vec = kmers[i]
        s = seqs[i]
        d_fwd = sqeuclidean(vec, ref_kmer)
        d_rev = sqeuclidean(vec, rev_ref_kmer)
        if d_fwd < d_rev
            push!(oriented, s)
            if phreds != nothing && length(phreds) > 0
                push!(newphreds, phreds[i])
            end
            if names != nothing && length(names) > 0
                push!(newnames, names[i])
            end
        else
            push!(oriented, reverse_complement(s))
            if phreds != nothing && length(phreds) > 0
                push!(newphreds, reverse(phreds[i]))
            end
            if names != nothing && length(names) > 0
                push!(newnames, names[i])
            end
        end
    end
    return oriented, newphreds, newnames
end

"""
    orient_to_refs_file(seqs::Array{String,1}, phreds::Union{Array{Vector{Phred},1},Void},
                        names::Union{Array{String,1},Void}, refspath::String)

Orients each sequence in `seqs` to nearest reference sequence in panel of references (`refspath`).
Distance determined by amino acid similarity (kmer vector dot prod).
"""
function orient_to_refs_file(seqs::Array{String,1}, phreds::Union{Array{Vector{Phred},1},Void}, names::Union{Array{String,1},Void}, refspath::String)
    if length(refspath) < 5
        return seqs
    end
    if refspath[end-5:end] == ".fasta"
        refs, _ = read_fasta_with_names_in_other_order(refspath, seqtype=String)
    elseif refspath[end-5:end] == ".fastq"
        refs, _, _ = read_fastq(refspath, seqtype=String)
    else
        error("Unsupported file type: $(infile[end-5:end])")
    end
    return orient_to_refs(seqs, phreds, names, refs)
end

"""
    orient_to_refs(seqs::Array{String,1}, phreds::Union{Array{Vector{Phred},1},Void},
                   names::Union{Array{String,1},Void}, refs::Array{String,1}; k::Int = 6)

Orients each sequence in `seqs` to nearest reference sequence in panel of references (`refs`).
Distance determined by amino acid similarity (kmer vector dot prod).
"""
function orient_to_refs(seqs::Array{String,1}, phreds::Union{Array{Vector{Phred},1},Void}, names::Union{Array{String,1},Void}, refs::Array{String,1}; k::Int = 6)
    if length(refs) == 0
        println("Not orienting: no reference sequences")
        return seqs
    end
    ref_kmers = [sparse_aa_kmer_count(rf, k) for rf in refs]
    rev_ref_kmers = [sparse_aa_kmer_count(reverse_complement(rf), k) for rf in refs]
    oriented = seqs[:]
    phredsexist = phreds != nothing && length(phreds) == length(seqs)
    namesexist = names != nothing && length(names) == length(seqs)
    function orientate(sq::String, ph::Vector{Phred}, name::String)
        best_match = 0
        best_match_i = 0
        rev = false
        vec = sparse_aa_kmer_count(sq, k)
        # compare to every reference, update min dists
        for i in 1:length(refs)
            match = dot(vec, ref_kmers[i])
            if match > best_match
                best_match = match
                best_match_i = i
                rev = false
            end
            match = dot(vec, rev_ref_kmers[i])
            if match > best_match
                best_match = match
                best_match_i = i
                rev = true
            end
        end
        if rev
            return reverse_complement(sq), reverse(ph), name
        elseif best_match_i > 0
            return sq, ph, name
        else
            # this probably happens because it does not share any kmers with any references
            # fix by shortening k?
            # make flag to throw these away?
            return sq, ph, name
        end
    end
    oriented = pmap(i->orientate(oriented[i],
                                 phredsexist ? phreds[i] : Phred[],
                                 namesexist ? names[i] : ""),
                    [1:length(oriented)...])
    return [o[1] for o in oriented], phredsexist ? [o[2] for o in oriented] : nothing, namesexist ? [o[3] for o in oriented] : nothing
end

#=
"""Orient sequences without separate references. Coarse clusters and takes
refined consensus of largest cluster as reference"""
function orient_strands(seqs::Array{String, 1}; radius = 0.3, k = 4)
=#


text:
# Circular ordering of edges around vertices
function edgeSlopeOrdering(VE,V,EV)
	VE_sorted = Array{Array{Int64,1},1}()
	forwardBackward = Tuple{Array{Int64,1},Array{Int64,1},Array{Int64,1}}[]
	for (v,ve) in enumerate(VE)
		ve_angle = []
		if ve != []
			for edge in ve
				w = pop!(setdiff(Set(EV[edge]), Set(v)))
				(x,y) = V[:,w] - V[:,v]
				angle = atan2(y,x)
				push!(ve_angle, 180*angle/pi)
			end
		end
		pairs = sort(collect(zip(ve_angle,ve)))
		sortedEdges = [pair[2] for pair in pairs]
		push!(VE_sorted, sortedEdges)
		theNexts,thePrevs = nextPrev(sortedEdges)
		push!(forwardBackward, (sortedEdges,theNexts,thePrevs))

	end
	forwardBackward
end

function boundary_1(EV)
	we = cellComplex(EV);
	rows = rowvals(we);
	for e in 1:size(we,2)
		i,j = collect(nzrange(we,e))
		h,k = rows[i], rows[j]
		if h < k
			we[h,e] = -1
		else
			we[k,e] = -1
		end
	end;
	we
end

# Compute the codomain of functions (next,prev) for the edge indices
function nextPrev(edges::Array{Int64,1})
	theNexts,thePrevs = Array{Int64,1}(), Array{Int64,1}()
	n = length(edges)
	for k in 1:n
		push!(theNexts, k==n ? edges[1] : edges[k+1] )
		push!(thePrevs, k==1 ? edges[n] : edges[k-1] )
	end
	theNexts,thePrevs
end

function rowColIncidence(BND_1,n)
	vs,es = findn(BND_1)
	VE = [Int64[] for k=1:n]
	for k in 1:length(vs)
		push!(VE[vs[k]], es[k])
	end
	return VE
end

function edgeSlopeOrdering(VE,V,EV)
 	VE_sorted = Array{Array{Int64,1},1}()
 	forwardBackward = Tuple{Array{Int64,1},Dict{Int64,Int64},Dict{Int64,Int64}}[]
 	for (v,ve) in enumerate(VE)
 		ve_angle = []
 		if ve != []
 			for edge in ve
 				w = pop!(setdiff(Set(EV[:,edge]), Set(v)))
 				(x,y) = V[:,w] - V[:,v]
 				angle = atan2(y,x)
 				push!(ve_angle, 180*angle/pi)
 			end
 		end
 		pairs = sort(collect(zip(ve_angle,ve)))
 		sortedEdges = [pair[2] for pair in pairs]
 		push!(VE_sorted, sortedEdges)
 		theNexts,thePrevs = nextPrev(sortedEdges)
 		next = Dict(zip(sortedEdges,theNexts))
 		prev = Dict(zip(sortedEdges,thePrevs))
 		push!(forwardBackward, (sortedEdges,next,prev))
 	end
 	forwardBackward
end

function checkOrient(BND_1,m,Cell_0,pivot,next)
	sign = Cell_0[1,pivot]
	C_1 = sparse([next],[1],[1],m,1)
	C_0 = (BND_1 * C_1)'
	chain_0 = findnz(C_0)
	_,cells_0,vals_0 = chain_0
	k = findfirst((x->x==pivot),cells_0)
	if vals_0[k]*sign==sign
		return -sign
	else
		return sign
	end
end

function cellTracking(BND_1,m,n,forwardBackward, V,EV, thecell,thesign)
	chain = [thecell]
	signs = [thesign]
	signedChain = [thecell*thesign]
	C_1 = sparse(chain,[1],signs,n,1)
	C_0 = (BND_1 * C_1)'
	chain_0 = findnz(C_0)
	_,cells_0,vals_0 = chain_0
	while cells_0 != []
		for k in 1:length(cells_0)
			pivot = cells_0[k]
			Cell_0 = sparse([1],[pivot],[vals_0[k]],1,m)
			C_1 = (Cell_0 * BND_1)
			chain_1 = findnz(C_1)
			_,cells_1,vals_1 = chain_1
			hinge = pop!(intersect(Set(cells_1), Set(chain)))
			if Cell_0[1,pivot] == -1
				next = forwardBackward[pivot][3][hinge]
			else
				next = forwardBackward[pivot][2][hinge] end
			thesign = checkOrient(BND_1,n,Cell_0,pivot,next)
			signedCell = Int64(thesign*next)
			if !(signedCell in signedChain)
				signedChain = push!(signedChain,signedCell) end
		end
		chain = [abs(cell) for cell in signedChain]
		signs = [sign(cell) for cell in signedChain]
		C_1 = sparse(chain,[1 for cell in chain],signs,n,1)
		C_0 = (BND_1 * C_1)'
		chain_0 = findnz(C_0)
		_,cells_0,vals_0 = chain_0
	end

	return signedChain
end

function boundary(V::Array{Float64,2},EV::Array{Int64,2})
	BND_1 = boundary_1(EV)
	m = size(BND_1,1)
	n = size(BND_1,2)
	VE = rowColIncidence(BND_1,n)
	forwardBackward = edgeSlopeOrdering(VE,V,EV)

	out = Array{Int64,1}[]
	store1 = [100 for k=1:size(EV,2)]
	store2 = [100 for k=1:size(EV,2)]
	cell = 1
	thesign = 1
	while cell != 0
		#@show cell

		signedChain = cellTracking(BND_1,m,n,forwardBackward,V,EV, cell,thesign)
		push!(out,signedChain)
		# store the cells in signedChain
		for c in signedChain
			if store1[abs(c)] == 100
				store1[abs(c)] = sign(c)
			elseif store2[abs(c)] == 100
				store2[abs(c)] = sign(c)
			end
		end
		# get the next seed
		pairs = collect(zip(store1,store2))
		cell = findfirst(c->99<=sum(c)<=101, pairs)
		if cell != 0
			thesign = -(sum(pairs[cell])-100)
		else
			cell = findfirst(c->sum(c)==200, pairs)
			if cell==0 break end
		end
	end
	pairs = collect(zip(store1,store2))
	out
end

function boundaryOp(EV::Array{Int64,2}, chains::Array{Array{Int64,1},1})
	I,J,V = Int64[],Int64[],Int64[]
	for (k,chain) in enumerate(chains)
		append!(I, abs(chain))
		append!(J, [k for row in chain])
		append!(V, sign(chain))
	end
	return sparse(I,J,V,size(EV,2),length(chains))
end


text:
# https://en.wikipedia.org/wiki/Letter_frequency
function text_frequencies()
    (
        0.082,
        0.015,
        0.028,
        0.043,
        0.130,
        0.022,
        0.020,
        0.061,
        0.070,
        0.0015,
        0.0077,
        0.040,
        0.025,
        0.067,
        0.075,
        0.019,
        0.00095,
        0.060,
        0.063,
        0.091,
        0.028,
        0.0098,
        0.024,
        0.0015,
        0.020,
        0.00074,
    )
end

# https://en.wikipedia.org/wiki/Letter_frequency
function dictionary_frequencies()
    (
        0.078,
        0.02,
        0.04,
        0.038,
        0.11,
        0.014,
        0.03,
        0.023,
        0.082,
        0.0021,
        0.025,
        0.053,
        0.027,
        0.072,
        0.061,
        0.028,
        0.0024,
        0.073,
        0.087,
        0.067,
        0.033,
        0.01,
        0.0091,
        0.0027,
        0.016,
        0.0044,
    )
end

"""
Maximize text frequency of unique letters.
"""
function text_frequency_score(candidate::Candidate)
    letter_flags = LetterFlags(candidate)
    unnormalized_score = dot(letter_flags.flags, text_frequencies())
    # unnormalized_score / (5 * maximum(text_frequencies()))
    unnormalized_score / 0.65
end

"""
Maximize dictionary frequency of unique letters.
"""
function dictionary_frequency_score(candidate::Candidate)
    letter_flags = LetterFlags(candidate)
    unnormalized_score = dot(letter_flags.flags, dictionary_frequencies())
    # unnormalized_score / (5 * maximum(dictionary_frequencies()))
    unnormalized_score / 0.55
end

"""
Which letters have never been seen in previous trials?
"""
function never_seen(trials::Vector{Trial})
    previously_seen = fill(false, MVector{26, Bool})

    for trial in trials
        previously_seen .|= LetterFlags(trial.candidate).flags
    end

    .!previously_seen
end

"""
Maximize number of unique letters that were not already guessed.
"""
function novelty_score(candidate::Candidate, never_seen)
    sum(LetterFlags(candidate).flags .& never_seen) / 5
end

function score_candidates(
    candidates::Vector{Candidate},
    trials::Vector{Trial},
    weights::Vector,
)
    never_seen_chars = never_seen(trials)

    score_functions = (
        c -> text_frequency_score(c),
        c -> dictionary_frequency_score(c),
        c -> novelty_score(c, never_seen_chars),
    )

    @assert length(score_functions) == length(weights)

    scores = fill(0.0, length(candidates))

    for (i, f) in enumerate(score_functions)
        for (j, candidate) in enumerate(candidates)
            scores[j] += weights[i] * f(candidate)
        end
    end

    scores
end