text:
# Formulation K2: (if regul==:classic, adds additional regularization parmeters -ρ (top left) and δ (bottom right))
# [-Q - D     A' ] [x] = rhs
# [ A         0  ] [y]
export K2LDLDenseParams

"""
Type to use the K2 formulation with a LDLᵀ factorization.
The outer constructor

    sp = K2LDLDenseParams(; ρ0 = sqrt(eps()) * 1e5, δ0 = sqrt(eps()) * 1e5)

creates a [`RipQP.SolverParams`](@ref) that should be used to create a [`RipQP.InputConfig`](@ref).
"""
mutable struct K2LDLDenseParams <: AugmentedParams
  uplo::Symbol
  fact_alg::Symbol
  ρ0::Float64
  δ0::Float64
end

function K2LDLDenseParams(;
  fact_alg::Symbol = :bunchkaufman,
  ρ0::Float64 = sqrt(eps()) * 1e5,
  δ0::Float64 = sqrt(eps()) * 1e5,
)
  uplo = :L # mandatory for LDL fact
  return K2LDLDenseParams(uplo, fact_alg, ρ0, δ0)
end

mutable struct PreallocatedDataK2LDLDense{T <: Real, S, M <: AbstractMatrix{T}} <:
               PreallocatedDataAugmentedLDL{T, S}
  D::S # temporary top-left diagonal
  regu::Regularization{T}
  K::M # augmented matrix
  fact_alg::Symbol
end

# outer constructor
function PreallocatedData(
  sp::K2LDLDenseParams,
  fd::QM_FloatData{T},
  id::QM_IntData,
  itd::IterData{T},
  pt::Point{T},
  iconf::InputConfig{Tconf},
) where {T <: Real, Tconf <: Real}

  # init Regularization values
  D = similar(fd.c, id.nvar)
  if iconf.mode == :mono
    regu = Regularization(T(sp.ρ0), T(sp.δ0), 1e-5 * sqrt(eps(T)), 1e0 * sqrt(eps(T)), :classic)
    D .= -T(1.0e0) / 2
  else
    regu =
      Regularization(T(sp.ρ0), T(sp.δ0), T(sqrt(eps(T)) * 1e0), T(sqrt(eps(T)) * 1e0), :classic)
    D .= -T(1.0e-2)
  end

  K = Symmetric(zeros(T, id.nvar + id.ncon, id.nvar + id.ncon), :L)
  K.data[1:(id.nvar), 1:(id.nvar)] .= .-fd.Q.data .+ Diagonal(D)
  K.data[(id.nvar + 1):(id.nvar + id.ncon), 1:(id.nvar)] .= fd.A
  K.data[view(diagind(K), (id.nvar + 1):(id.nvar + id.ncon))] .= regu.δ

  if sp.fact_alg == :bunchkaufman
    bunchkaufman!(K)
  elseif sp.fact_alg == :ldl
    ldl_dense!(K)
  end

  return PreallocatedDataK2LDLDense(
    D,
    regu,
    K, #K
    sp.fact_alg,
  )
end

# function used to solve problems
# solver LDLFactorization
function solver!(
  dd::AbstractVector{T},
  pad::PreallocatedDataK2LDLDense{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
  step::Symbol,
) where {T <: Real}
  ldiv_dense!(pad.K, dd)
  return 0
end

function update_pad!(
  pad::PreallocatedDataK2LDLDense{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
) where {T <: Real}
  if cnts.k != 0
    update_regu!(pad.regu)
  end

  pad.D .= -pad.regu.ρ
  pad.D[id.ilow] .-= pt.s_l ./ itd.x_m_lvar
  pad.D[id.iupp] .-= pt.s_u ./ itd.uvar_m_x
  pad.K.data[1:(id.nvar), 1:(id.nvar)] .= .-fd.Q.data .+ Diagonal(pad.D)
  pad.K.data[(id.nvar + 1):(id.nvar + id.ncon), 1:(id.nvar)] .= fd.A
  pad.K.data[(id.nvar + 1):(id.nvar + id.ncon), (id.nvar + 1):(id.nvar + id.ncon)] .= zero(T)
  pad.K.data[view(diagind(pad.K), (id.nvar + 1):(id.nvar + id.ncon))] .= pad.regu.δ

  if pad.fact_alg == :bunchkaufman
    bunchkaufman!(pad.K)
  elseif pad.fact_alg == :ldl
    ldl_dense!(pad.K)
  end

  return 0
end


text:
module BrokemiaHelperCursedTeleportTrigger

using ..Ahorn, Maple

@mapdef Trigger "brokemiahelper/cursedTeleportTrigger" CursedTeleportTrigger(x::Integer, y::Integer, width::Integer=8, height::Integer=8, endRoom::String="")

const placements = Ahorn.PlacementDict(
    "CursedTeleportTrigger" => Ahorn.EntityPlacement(
        CursedTeleportTrigger,
        "rectangle",
        Dict{String, Any}()
    )
)

end

ext: jl

max_stars_repo_name: Brokemia/BrokemiaHelper

max_issues_repo_name: Brokemia/BrokemiaHelper

max_forks_repo_name: Brokemia/BrokemiaHelper


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))

using Dates, Printf, ProgressMeter, Random

include("common.jl")

largefile_coln = 20
space_per_row = 64 * 20 # 64 bit (for a Float64) times 20 columns
space_per_row = space_per_row / 8 / 1024 / 1024 / 1024 # Convert to GB
largefile_rown = convert(Int, round(10 / space_per_row)) # We want to use 10 GB for the results

if isfile(joinpath(@__DIR__), "local_writeconfig.jl")
    include("local_writeconfig.jl")
end

function our_copy(source, dest)
    if Sys.iswindows()
        readlines(`cmd /c copy $source $dest`)
    else
        readlines(`cp $source $dest`)
    end
end

function write_file(col_types, rown, withna, filename; make_warmup_copy=true, folder_path=nothing)
    coln = length(col_types)
    if folder_path===nothing
        folder_path = ourpath(rown, coln, withna)
    end
    file_path = joinpath(folder_path, filename)
    warmup_file_path = joinpath(folder_path, "warmup_" * filename)
    mkpath(folder_path)
    open(file_path, "w") do f
        # Write header
        for c in 1:coln
            c>1 && print(f, ",")
            print(f, "col", c)
        end
        println(f)

        # Write data
        @showprogress for r in 1:rown
            for c in 1:coln
                c>1 && print(f, ",")
                if col_types[c] == :float64
                    if !withna || rand(Bool)
                        print(f, rand(Float64))
                    end
                elseif col_types[c] == :float64short
                    if !withna || rand(Bool)
                        print(f, @sprintf("%.4f", rand(Float64)))
                    end
                elseif col_types[c] == :int64
                    if !withna || rand(Bool)
                        print(f, rand(Int64))
                    end
                elseif col_types[c] == :datetime
                    if !withna || rand(Bool)
                        print(f, DateTime(rand(1950:2000), rand(1:12), rand(1:28), rand(1:23), rand(1:59), rand(1:50)))
                    end
                elseif col_types[c] == :string
                    if !withna || rand(Bool)
                        print(f, '"', randstring(20), '"')
                    end
                elseif col_types[c] == :stringcat
                    if !withna || rand(Bool)
                        print(f, '"', "Categorical string ", rand(1:5), '"')
                    end
                elseif col_types[c] == :stringescaped
                    if !withna || rand(Bool)
                        print(f, '"', randstring(20), "\"\"", randstring(10), '"')
                    end
                end
            end
            println(f)
        end
    end
    make_warmup_copy && our_copy(file_path, warmup_file_path)
end

uniform_types = [:float64, :float64short, :int64, :datetime, :string, :stringcat, :stringescaped]

for n in ns, c in cs, withna in [true,false]
    println("Writing rows=$n, columns=$c, withna=$withna")
    for typ in uniform_types
        println("    $typ")
        write_file(fill(typ, c), n, withna, string("uniform_", lowercase(string(typ)), ".csv"))
    end
    println("    mixed")
    write_file(repeat(uniform_types, div(c, length(uniform_types)) + 1)[1:c], n, withna, "mixed.csv")
end

write_file(fill(:float64, 20), largefile_rown, true, "float64.csv", make_warmup_copy = false, folder_path = joinpath(@__DIR__, "..", "data", "large"))

max_stars_repo_name: davidanthoff/csv-comparison

max_issues_repo_name: davidanthoff/csv-comparison

max_forks_repo_name: davidanthoff/csv-comparison


text:
@testset "Utils" begin
    s = @inferred(LightGraphs.sample!([1:10;], 3))
    @test length(s) == 3
    for e in s
        @test 1 <= e <= 10
    end

    s = @inferred(LightGraphs.sample!([1:10;], 6, exclude=[1, 2]))
    @test length(s) == 6
    for e in s
        @test 3 <= e <= 10
    end

    s = @inferred(LightGraphs.sample(1:10, 6, exclude=[1, 2]))
    @test length(s) == 6
    for e in s
        @test 3 <= e <= 10
    end
end


text:
function real_to_parametric(p::Point2D, poly::QuadraticPolygon{N, 2, T}) where {N, T}
    return real_to_parametric(p, poly, 30)
end

function real_to_parametric(p::Point2D{T}, poly::QuadraticPolygon{N, 2, T},
                            max_iters::Int64) where {N, T}
    # Convert from real coordinates to the polygon's local parametric coordinates using
    # Newton-Raphson.
    # If a conversion doesn't exist, the minimizer is returned.
    # Initial guess at polygon centroid
    if N === 6 # Triangle
        rs = SVector{2,T}(1//3, 1//3)
    else # Quadrilateral
        rs = SVector{2,T}(1//2, 1//2)
    end
    for i ∈ 1:max_iters
        Δrs = inv(𝗝(poly, rs[1], rs[2]))*(p - poly(rs[1], rs[2]))
        if Δrs ⋅ Δrs < T((1e-8)^2)
            break
        end
        rs += Δrs
    end
    return Point2D{T}(rs[1], rs[2])
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# v0.6

function bestshuffle(str::String)::Tuple{String,Int}
    s = Vector{Char}(str)

    # Count the supply of characters.
    cnt = Dict{Char,Int}(c => 0 for c in s)
    for c in s; cnt[c] += 1 end

    # Allocate the result
    r = similar(s)
    for (i, x) in enumerate(s)
        # Find the best character to replace x.
        best = x
        rankb = -2
        for (c, rankc) in cnt
            # Prefer characters with more supply.
            # (Save characters with less supply.)
            # Avoid identical characters.
            if c == x; rankc = -1 end
            if rankc > rankb
                best = c
                rankb = rankc
            end
        end

        # Add character to list. Remove it from supply.
        r[i] = best
        cnt[best] -= 1
        if cnt[best] == 0; delete!(cnt, best) end
    end

    # If the final letter became stuck (as "ababcd" became "bacabd",
    # and the final "d" became stuck), then fix it.
    i = length(s)
    if r[i] == s[i]
        for j in 1:i
            if r[i] != s[j] && r[j] != s[i]
                r[i], r[j] = r[j], r[i]
                break
            end
        end
    end

    score = sum(x == y for (x, y) in zip(r, s))
    return r, score
end

for word in ("abracadabra", "seesaw", "elk", "grrrrrr", "up", "a")
    shuffled, score = bestshuffle(word)
    println("$word: $shuffled ($score)")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# helper module with misc SharedArray utilities
module joSAutils
    using Distributed
    using SharedArrays
    using JOLI: jo_convert, joAbstractLinearOperator, joPAsetup

    # A::joAbstractLinearOperator
    function jo_x_mv!(A::joAbstractLinearOperator,din::Tuple{Vararg{UnitRange{INT}}},dout::Tuple{Vararg{UnitRange{INT}}},
            in::SharedArray{ADDT,2},out::SharedArray{ARDT,2}) where {ADDT,ARDT,INT<:Integer}

            out[dout...]=jo_convert(ARDT,F*in[din...])
        return nothing
    end
    function jo_x_mv!(A::joAbstractLinearOperator,din::joPAsetup,dout::joPAsetup,
            in::SharedArray{ADDT,2},out::SharedArray{ARDT,2}) where {ADDT,ARDT}

        @assert out.pids==in.pids "SharedArrays pids in($(in.pids)) and out($(out.pids)) do not match"
        P::Int=length(out.pids)
        @sync begin
            for w=1:P
                p=out.pids[w]
                @async remotecall_wait(jo_x_mv!,p,A,din.idxs[w],dout.idxs[w],in,out)
            end
        end
        return nothing
    end
    # F::Function
    function jo_x_mv!(F::Function,din::Tuple{Vararg{UnitRange{INT}}},dout::Tuple{Vararg{UnitRange{INT}}},
            in::SharedArray{ADDT,2},out::SharedArray{ARDT,2}) where {ADDT,ARDT,INT<:Integer}

            out[dout...]=jo_convert(ARDT,F(in[din...]))
        return nothing
    end
    function jo_x_mv!(F::Function,din::joPAsetup,dout::joPAsetup,
            in::SharedArray{ADDT,2},out::SharedArray{ARDT,2}) where {ADDT,ARDT}

        @assert out.pids==in.pids "SharedArrays pids in($(in.pids)) and out($(out.pids)) do not match"
        P::Int=length(out.pids)
        @sync begin
            for w=1:P
                p=out.pids[w]
                @async remotecall_wait(jo_x_mv!,p,F,din.idxs[w],dout.idxs[w],in,out)
            end
        end
        return nothing
    end
end
using .joSAutils

export salloc
"""
    julia> salloc(d; [DT])

Allocates a SharedArray according to given distributor

# Signature

    salloc(d::joPAsetup;DT::DataType=d.DT)

# Arguments

- `d`: see help for joPAsetup
- `DT`: keyword argument to overwrite the type in joPAsetup

# Examples

- `salloc(d)`: allocate an array
- `salloc(d,DT=Float32)`: allocate array and overwite d.DT with Float32

"""
function salloc(d::joPAsetup;DT::DataType=d.DT)
    #bug S = SharedArray{DT}(d.dims; pids=d.procs)
    S = SharedArray{DT}(deepcopy(d.dims); pids=deepcopy(d.procs))
    #S = SharedArray{DT}(d.dims)
    return S
end

export szeros
function szeros(d::joPAsetup;DT::DataType=d.DT)
    fill=S->S[SharedArrays.localindices(S)]=zeros(DT,length(SharedArrays.localindices(S)))
    #bug S = SharedArray{DT}(d.dims; init=fill, pids=d.procs)
    S = SharedArray{DT}(deepcopy(d.dims); init=fill, pids=deepcopy(d.procs))
    #S = SharedArray{DT}(d.dims; init=fill)
    return S
end

export sones
function sones(d::joPAsetup;DT::DataType=d.DT)
    fill=S->S[SharedArrays.localindices(S)]=ones(DT,length(SharedArrays.localindices(S)))
    #bug S = SharedArray{DT}(d.dims; init=fill, pids=d.procs)
    S = SharedArray{DT}(deepcopy(d.dims); init=fill, pids=deepcopy(d.procs))
    #S = SharedArray{DT}(d.dims; init=fill)
    return S
end

export srand
function srand(d::joPAsetup;DT::DataType=d.DT)
    fill=S->S[SharedArrays.localindices(S)]=rand(DT,length(SharedArrays.localindices(S)))
    #bug S = SharedArray{DT}(d.dims; init=fill, pids=d.procs)
    S = SharedArray{DT}(deepcopy(d.dims); init=fill, pids=deepcopy(d.procs))
    #S = SharedArray{DT}(d.dims; init=fill)
    return S
end

export srandn
function srandn(d::joPAsetup;DT::DataType=d.DT)
    fill=S->S[SharedArrays.localindices(S)]=randn(DT,length(SharedArrays.localindices(S)))
    #bug S = SharedArray{DT}(d.dims; init=fill, pids=d.procs)
    S = SharedArray{DT}(deepcopy(d.dims); init=fill, pids=deepcopy(d.procs))
    #S = SharedArray{DT}(d.dims; init=fill)
    return S
end

export scatter
"""
    julia> scatter(A,d)

Scatters SharedArray according to given joPAsetup.

# Signature

    scatter(a::AbstractArray,d::joPAsetup)

# Arguments

- `A`: array to ditribute
- `d`: see help for joPAsetup

# Notes

- the type in joPAsetup is ignored here
- one of the dimensions must be large enough to hold at least one element on each worker

# Examples

- `scatter(A,d)`: scatter A using given distributor
- `scatter(A,joPAsetup(size(A)...))`: scatter A using default distributor settings

"""
function scatter(A::AbstractArray,d::joPAsetup)
    size(A)==d.dims || throw(joPAsetupException("joPAsetup: array size does not match dims of joPAsetup"))
    #bug SA = SharedArray{eltype(A)}(d.dims; pids=d.procs)
    SA = SharedArray{eltype(A)}(deepcopy(d.dims); pids=deepcopy(d.procs))
    #SA = SharedArray{eltype(A)}(d.dims)
    SA[:] = A[:]
    return SA
end


text:
using Distributions
using StatsBase

function Binom(k, N, p)
    return pdf(Binomial(N, p), k)
end

function Binom_sum(N, p)
    if iseven(N)
        N -= 1
    end
    return sum([Binom(k, N, p) for k = N:-1:N/2])
end

function compute_analytical_accuracy(; N, p_GC_sig, p_GC_bkg, weight = [1, 1])
    p_α_sig = 1 - p_GC_sig
    p_α_bkg = 1 - p_GC_bkg
    return mean([Binom_sum(N, p_α_sig), Binom_sum(N, p_α_bkg)], weights(weight))
end

function compute_analytical_accuracies(;
    N_half,
    p_GC_sig = 0.4,
    p_GC_bkg = 0.5,
    weight = [1, 1],
)
    bases_included = 1:N_half
    analytical_accuracies = [
        compute_analytical_accuracy(
            N = i * 2,
            p_GC_sig = p_GC_sig,
            p_GC_bkg = p_GC_bkg,
            weight = weight,
        ) for i in bases_included
    ]
    return analytical_accuracies
end

#%%

using CairoMakie
using ColorSchemes

N_half = 38
x = 1:N_half

analytical_accuracies = compute_analytical_accuracies(;N_half=N_half)

colormap = [x for x in ColorSchemes.Set1_9.colors]

#%%

f = Figure()
ax = Axis(
    f[1, 1],
    title = "Analytical accuracy",
    xlabel = "# bases included (centered)",
    ylabel = "Accuracy",
    # limits = (0.5, half_seq_length + 0.5, 0.634, 0.701),
    # limits = (0.5, half_seq_length + 0.5, ylimits...),
    xticks = 1:2:N_half,
)

scatterlines!(
    ax,
    x,
    analytical_accuracies,
    color = colormap[3],
    markercolor = colormap[3],
    label = "Analytical",
)

axislegend(position = :rb)
f

#%%

compute_analytical_accuracy(N=1, p_GC_sig = 0.4, p_GC_bkg = 0.5)
compute_analytical_accuracy(N=2, p_GC_sig = 0.4, p_GC_bkg = 0.5)
compute_analytical_accuracy(N=3, p_GC_sig = 0.4, p_GC_bkg = 0.5)
compute_analytical_accuracy(N=4, p_GC_sig = 0.4, p_GC_bkg = 0.5)
compute_analytical_accuracy(N=5, p_GC_sig = 0.4, p_GC_bkg = 0.5)
compute_analytical_accuracy(N=6, p_GC_sig = 0.4, p_GC_bkg = 0.5)

max_stars_repo_name: ChristianMichelsen/reference-free

max_issues_repo_name: ChristianMichelsen/reference-free

max_forks_repo_name: ChristianMichelsen/reference-free


text:
module OverflowContexts
__precompile__(false)

include("macros.jl")
include("base_ext.jl")

export @default_checked, @default_unchecked, @checked, @unchecked,
    unchecked_neg, unchecked_add, unchecked_sub, unchecked_mul, unchecked_negsub, unchecked_pow, unchecked_abs,
    checked_neg, checked_add, checked_sub, checked_mul, checked_pow, checked_negsub, checked_abs,
    SignedBitInteger

end # module


text:
# This file was generated, do not modify it. # hide
r = range(xgb, :num_round, lower=50, upper=500)
curve = learning_curve(xgbm, range=r, resolution=50,
                        measure=L1HingeLoss())

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@testset "Plotting" begin

    fname = joinpath(dirname(@__FILE__), "../data", "test_Hz19.5-testing.bdf")
    s = read_SSR(fname)
    s = rereference(s, "Cz")
    s = highpass_filter(s)
    s = extract_epochs(s)
    s = create_sweeps(s, epochsPerSweep = 2)
    s = ftest(s)

    @testset "Spectrum" begin
        p = plot_spectrum(s, "20Hz_SWN_70dB_R", targetFreq = 3.0)
        #= display(p) =#

        p = plot_spectrum(vec(s.data[:, 1]), Int(samplingrate(s)), dBPlot = false)
        #= display(p) =#

        p = plot_spectrum(s, 3, targetFreq = 40.0390625)
        #= display(p) =#
    end

    @testset "Filter reponse" begin

        p = plot_filter_response(s.processing["filter1"], Int(samplingrate(s)))
        #= display(p) =#
    end

    s = trim_channel(s, 8192*3)

    @testset "Multi channel time series" begin
        plot1 = plot_timeseries(s)
        #= display(plot1) =#

        plot2 = plot_timeseries(s, channels=["40Hz_SWN_70dB_R", "Cz"])
        #= display(plot2) =#

        s = rereference(s, "car")
        plot3 = plot_timeseries(s, channels=["40Hz_SWN_70dB_R", "Cz"])
        #= display(plot3) =#
    end

    @testset "Single channel time series" begin
        plot4 = plot_timeseries(s, channels=["40Hz_SWN_70dB_R"])
        #= display(plot4) =#

        plot5 = plot_timeseries(s, channels="Cz")
        #= display(plot5) =#

        keep_channel!(s, ["40Hz_SWN_70dB_R"])
        plot6 = plot_timeseries(s)
        #= display(plot6) =#
    end
end


text:
## --- Sorting and counting array elements

    """
    ```julia
    n = count_unique!(A)
    ```
    Sort the array `A` in-place (if not already sorted), move unique elements to
    the front, and return the number of unique elements found.

    `A[1:count_unique!(A)]` should return an array equivalent to `unique(A)`
    """
    function count_unique!(A)
        issorted(A) || sort!(A)
        n = 1
        last = A[1]
        @inbounds for i=2:length(A)
            if A[i] != last
                n += 1
                last = A[n] = A[i]
            end
        end
        return n
    end
    export count_unique!

## --- Convert between bin centers and bin edges

    """
    ```julia
    cntr(edges::AbstractArray)
    ```
    Given an array of bin edges, return a corresponding vector of bin centers
    """
    function cntr(edges::AbstractArray)
        centers = (edges[1:end-1] + edges[2:end]) ./ 2
        return centers
    end
    export cntr

## --- Searching arrays

    """
    ```julia
    findmatches(source, target)
    ```
    Return the index of the first value in `target` (if any) that is equal to
    a given value in `source` for each value in `source`; else 0.
    """
    function findmatches(source, target)
        index = Array{Int64}(undef, size(source))
        return findmatches!(index, source, target)
    end
    function findmatches!(index, source, target)
        # Loop through source and find first match for each (if any)
        @inbounds for i = 1:length(index)
            for j = 1:length(target)
                index[i] = 0
                if isequal(source[i], target[j])
                    index[i] = j
                    break
                end
            end
        end
        return index
    end
    export findmatches, findmatches!

    """
    ```julia
    findclosest(source, target)
    ```
    Return the index of the numerically closest value in the indexable collection
    `target` for each value in `source`.
    If muliple values are equally close, the first one is used
    """
    function findclosest(source, target)
        index = Array{Int64}(undef, size(source))
        return findclosest!(index, source, target)
    end
    function findclosest!(index, source, target)
        # Find closest (numerical) match in target for each value in source
        @inbounds for i = 1:length(source)
            d = abs(target[1] - source[i])
            index[i] = 1
            for j = 2:length(target)
                d_prop = abs(target[j] - source[i])
                if d_prop < d
                    d = d_prop
                    index[i] = j
                end
            end
        end
        return index
    end
    export findclosest, findclosest!

    """
    ```julia
    findclosestbelow(source, target)
    ```
    Return the index of the nearest value of the indexable collection `target`
    that is less than (i.e., "below") each value in `source`.
    If no such target values exist in `target`, returns an index of 0.
    """
    function findclosestbelow(source, target)
        index = Array{Int64}(undef, size(source))
        return findclosestbelow!(index, source, target)
    end
    function findclosestbelow!(index, source, target)
        @inbounds for i = 1:length(source)
            index[i] = d = j = 0
            while j < length(target)
                j += 1
                if target[j] < source[i]
                    d = source[i] - target[j]
                    index[i] = j
                    break
                end
            end
            while j < length(target)
                j += 1
                if target[j] < source[i]
                    d_prop = source[i] - target[j]
                    if d_prop < d
                        d = d_prop
                        index[i] = j
                    end
                end
            end
        end
        return index
    end
    export findclosestbelow, findclosestbelow!

    """
    ```julia
    findclosestabove(source, target)
    ```
    Return the index of the nearest value of the indexable collection `target`
    that is greater than (i.e., "above") each value in `source`.
    If no such values exist in `target`, returns an index of 0.
    """
    function findclosestabove(source, target)
        index = Array{Int64}(undef, size(source))
        return findclosestabove!(index,source,target)
    end
    function findclosestabove!(index, source, target)
        @inbounds for i = 1:length(source)
            index[i] = d = j = 0
            while j < length(target)
                j += 1
                if target[j] > source[i]
                    d = target[j] - source[i]
                    index[i] = j
                    break
                end
            end
            while j < length(target)
                j += 1
                if target[j] > source[i]
                    d_prop = target[j] - source[i]
                    if d_prop < d
                        d = d_prop
                        index[i] = j
                    end
                end
            end
        end
        return index
    end
    export findclosestabove, findclosestabove!

    """
    ```julia
    findnth(t::AbstractArray{Bool}, n::Integer)
    ```
    Return the index of the `n`th true value of `t`, else length(`t`)
    """
    function findnth(t::AbstractArray{Bool}, n::Integer)
        N = 0
        @inbounds for i ∈ eachindex(t)
            if t[i]
                N += 1
            end
            if N == n
                return i
            end
        end
        return length(t)
    end
    export findnth

## --- String matching

    # if ~ @isdefined contains
    #     """
    #     ```julia
    #     contains(haystack, needle)
    #     ```
    #
    #     Converts both `haystack` and `needle` to strings (if not already strings)
    #     and checks whether `string(haystack)` contains `string(needle)`.
    #     """
    #     contains(haystack::AbstractString, needle::Union{AbstractString,Regex,AbstractChar}) = occursin(needle, haystack)
    #     contains(haystack, needle) = occursin(string(needle), string(haystack))
    #     export contains
    # end

    """
    ```julia
    containsi(haystack, needle)
    ```

    Converts both `haystack` and `needle` to strings and checks whether
    `string(haystack)` contains `string(needle)`, ignoring case.
    """
    containsi(haystack::AbstractString, needle::Union{AbstractString,AbstractChar}) = occursin(lowercase(needle), lowercase(haystack))
    containsi(haystack, needle) = occursin(lowercase(string(needle)), lowercase(string(haystack)))
    export containsi

## --- Drawing a pseudorandom array from a numerically specified distribution

    """
    ```julia
    x = draw_from_distribution(dist::AbstractArray{<:AbstractFloat}, n::Integer)
    ```
    Draw `n` random floating point numbers from a continuous probability distribution
    specified by a vector `dist` defining the PDF curve thereof.
    """
    function draw_from_distribution(dist::AbstractArray{<:AbstractFloat}, n::Integer)
        x = Array{eltype(dist)}(undef, n)
        draw_from_distribution!(x, dist)
        return x
    end
    export draw_from_distribution

    """
    ```julia
    draw_from_distribution!(dist::AbstractArray{<:AbstractFloat}, x::Array{<:AbstractFloat})
    ```
    Fill an existing variable `x` with random floating point numbers drawn from
    a continuous probability distribution specified by a vector `dist`
    defining the PDF curve thereof.
    """
    function draw_from_distribution!(x::Array{<:AbstractFloat}, dist::AbstractArray{<:AbstractFloat})
        # Fill the array x with random numbers from the distribution 'dist'
        dist_ymax = maximum(dist)
        dist_xmax = prevfloat(length(dist) - 1.0)

        @inbounds for i ∈ eachindex(x)
            while true
                # Pick random x value
                rx = rand(eltype(x)) * dist_xmax
                # Interpolate corresponding distribution value
                f = floor(Int,rx)
                y = dist[f+2]*(rx-f) + dist[f+1]*(1-(rx-f))
                # See if x value is accepted
                ry = rand(Float64) * dist_ymax
                if (y > ry)
                    x[i] = rx / dist_xmax
                    break
                end
            end
        end
    end
    export draw_from_distribution!

## --- Numerically integrate a 1-d distribution

    """
    ```julia
    trapezoidalquadrature(edges, values)
    ```
    Add up the area under a curve with y positions specified by a vector of `values`
    and x positions specfied by a vector of `edges` using trapezoidal integration.
    Bins need not be evenly spaced, though it helps (integration will be faster
    if `edges` are specified as an AbstractRange).
    """
    function trapezoidalquadrature(edges::AbstractRange, values::AbstractArray)
        result = zero(eltype(values))
        @turbo for i = (firstindex(edges)+1):lastindex(edges)
            result += values[i-1]+values[i]
        end
        dx = (edges[end]-edges[1])/(length(edges) - 1)
        return result * dx / 2
    end
    function trapezoidalquadrature(edges::AbstractArray, values::AbstractArray)
        result = zero(promote_type(eltype(edges), eltype(values)))
        @turbo for i = (firstindex(edges)+1):lastindex(edges)
            result += (values[i-1] + values[i]) * (edges[i] - edges[i-1])
        end
        return result / 2
    end
    export trapezoidalquadrature
    trapz = trapezoidalquadrature
    export trapz

    """
    ```julia
    midpointquadrature(bincenters, values)
    ```
    Add up the area under a curve with y positions specified by a vector of `values`
    and x positions specfied by a vector of `bincenters` using midpoint integration.
    """
    function midpointquadrature(bincenters::AbstractRange, values::AbstractArray)
        sum(values) * (last(bincenters)-first(bincenters)) / (length(bincenters) - 1)
    end
    export midpointquadrature

## --- End of File

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
module TestLoggedFunctions

using Test
import ThreadPools: LoggedStaticPool
using ThreadPools

include("util.jl")

@testset "log functions" begin

    @testset "read/write logs" begin
        N = 2 * Threads.nthreads()
        pool = LoggedStaticPool()
        tforeach(pool, x->sleep(0.01*x), 1:N)
        close(pool)
        dumplog("_tmp.log", pool)
        log2 = readlog("_tmp.log")
        @test pool.log == log2
        rm("_tmp.log")
    end

    @testset "showactivity" begin
        io = IOBuffer()
        showactivity(io, "$(@__DIR__)/testlog.txt", 0.1, nthreads=4)
        @test replace(String(take!(io)), r"\s+\n"=>"\n") == """0.000   -   -   -   -
0.100   4   1   3   2
0.200   4   5   3   2
0.300   4   5   3   6
0.400   4   5   7   6
0.500   8   5   7   6
0.600   8   5   7   6
0.700   8   -   7   6
0.800   8   -   7   6
0.900   8   -   7   -
1.000   8   -   7   -
1.100   8   -   -   -
1.200   8   -   -   -
1.300   -   -   -   -
1.400   -   -   -   -
1.500   -   -   -   -
"""
    end

    @testset "showstats" begin
        io = IOBuffer()
        showstats(io, "$(@__DIR__)/testlog.txt")
        @test String(take!(io)) == """

        Total duration: 1.212 s
        Number of jobs: 8
        Average job duration: 0.457 s
        Minimum job duration: 0.11 s
        Maximum job duration: 0.805 s

        Thread 1: Duration 1.211 s, Gap time 0.0 s
        Thread 2: Duration 0.616 s, Gap time 0.0 s
        Thread 3: Duration 1.024 s, Gap time 0.0 s
        Thread 4: Duration 0.805 s, Gap time 0.0 s
    """
    end

end

end # module


text:
# AbstractRGB

abstract type AbstractStorageRGB{T, fr, fg, fb, en} <: AbstractRGB{T} end

const AbstractRGB16{en} = Union{
    AbstractStorageRGB{N0f8, 5, 6, 5, en},
    AbstractStorageRGB{N0f8, 5, 5, 5, en},
    AbstractStorageRGB{N0f8, 4, 4, 4, en}}

"""
    RGB565LE

A 16-bit RGB type with 5-bit red, 6-bit green and 5-bit blue components.
Unlike [`RGB565BE`](@ref), the blue component is included in the first byte as
follows.
```
  7  6  5  4  3  2  1  0     7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
|G2 G1 G0|B4 B3 B2 B1 B0|  |R4 R3 R2 R1 R0|G5 G4 G3|
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
       first byte                 second byte
```
As the blue component is before the red component, `RGB565LE` is sometimes
called "BGR565". Note that "BGR565" may also refer to a variant of "RGB565"
swapped by component, not by byte.
"""
struct RGB565LE <: AbstractStorageRGB{N0f8, 5, 6, 5, :le}
    color::UInt16
    RGB565LE(color::UInt16, ::Val{:no}) = new(color)
end

"""
    RGB565BE

A 16-bit RGB type with 5-bit red, 6-bit green and 5-bit blue components.
Unlike [`RGB565LE`](@ref), the red component is included in the first byte as
follows.
```
  7  6  5  4  3  2  1  0     7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
|R4 R3 R2 R1 R0|G5 G4 G3|  |G2 G1 G0|B4 B3 B2 B1 B0|
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
       first byte                 second byte
```
"""
struct RGB565BE <: AbstractStorageRGB{N0f8, 5, 6, 5, :be}
    color::UInt16
    RGB565BE(color::UInt16, ::Val{:no}) = new(color)
end

"""
    RGB555LE

A 16-bit RGB type with 5-bit red, green and blue components.
Unlike [`RGB555BE`](@ref), the blue component is included in the first byte as
follows.
```
  7  6  5  4  3  2  1  0     7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
|G2 G1 G0|B4 B3 B2 B1 B0|  |X |R4 R3 R2 R1 R0|G4 G3|
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
       first byte                 second byte
```
As the blue component is before the red component, `RGB555LE` is sometimes
called "BGR555". Note that "BGR555" may also refer to a variant of "RGB555"
swapped by component, not by byte.
"""
struct RGB555LE <: AbstractStorageRGB{N0f8, 5, 5, 5, :le}
    color::UInt16
    RGB555LE(color::UInt16, ::Val{:no}) = new(color)
end

"""
    RGB555BE

A 16-bit RGB type with 5-bit red, green and blue components.
Unlike [`RGB555LE`](@ref), the red component is included in the first byte as
follows.
```
  7  6  5  4  3  2  1  0     7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
|X |R4 R3 R2 R1 R0|G4 G3|  |G2 G1 G0|B4 B3 B2 B1 B0|
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
       first byte                 second byte
```
"""
struct RGB555BE <: AbstractStorageRGB{N0f8, 5, 5, 5, :be}
    color::UInt16
    RGB555BE(color::UInt16, ::Val{:no}) = new(color)
end

"""
    XRGB444LE

A 16-bit RGB type with 4-bit red, green and blue components.
Unlike [`XRGB444BE`](@ref), the blue component is included in the first byte as
follows.
```
  7  6  5  4  3  2  1  0     7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
|G3 G2 G1 G0|B3 B2 B1 B0|  |X  X  X  X |R3 R2 R1 R0|
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
       first byte                 second byte
```
"""
struct XRGB444LE <: AbstractStorageRGB{N0f8, 4, 4, 4, :le}
    color::UInt16
    XRGB444LE(color::UInt16, ::Val{:no}) = new(color)
end

"""
    XRGB444BE

A 16-bit RGB type with 4-bit red, green and blue components.
Unlike [`XRGB444LE`](@ref), the red component is included in the first byte as
follows.
```
  7  6  5  4  3  2  1  0     7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
|X  X  X  X |R3 R2 R1 R0|  |G3 G2 G1 G0|B3 B2 B1 B0|
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
       first byte                 second byte
```
"""
struct XRGB444BE <: AbstractStorageRGB{N0f8, 4, 4, 4, :be}
    color::UInt16
    XRGB444BE(color::UInt16, ::Val{:no}) = new(color)
end

"""
    RGBX444LE

A 16-bit RGB type with 4-bit red, green and blue components.
Unlike [`RGBX444BE`](@ref), the blue component is included in the first byte as
follows.
```
  7  6  5  4  3  2  1  0     7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
|B3 B2 B1 B0|X  X  X  X |  |R3 R2 R1 R0|G3 G2 G1 G0|
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
       first byte                 second byte
```
"""
struct RGBX444LE <: AbstractStorageRGB{N0f8, 4, 4, 4, :le}
    color::UInt16
    RGBX444LE(color::UInt16, ::Val{:no}) = new(color)
end

"""
    RGBX444BE

A 16-bit RGB type with 4-bit red, green and blue components.
Unlike [`RGBX444LE`](@ref), the red component is included in the first byte as
follows.
```
  7  6  5  4  3  2  1  0     7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
|R3 R2 R1 R0|G3 G2 G1 G0|  |B3 B2 B1 B0|X  X  X  X |
+--+--+--+--+--+--+--+--+  +--+--+--+--+--+--+--+--+
       first byte                 second byte
```
"""
struct RGBX444BE <: AbstractStorageRGB{N0f8, 4, 4, 4, :be}
    color::UInt16
    RGBX444BE(color::UInt16, ::Val{:no}) = new(color)
end

"""
    RGB332

An 8-bit RGB type with 3-bit red, 3-bit green and 2-bit blue components.
```
  7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+
|R2 R1 R0|G2 G1 G0|B1 B0|
+--+--+--+--+--+--+--+--+
```
"""
struct RGB332 <: AbstractStorageRGB{N0f8, 3, 3, 2, :no}
    color::UInt8
    RGB332(color::UInt8, ::Val{:no}) = new(color)
end

"""
    RGB242

An 8-bit RGB type with 2-bit red, 4-bit green and 2-bit blue components.

This type has a single green component, which is different from the "RGGB" (with
two green components) used in the Bayer pattern.
```
  7  6  5  4  3  2  1  0
+--+--+--+--+--+--+--+--+
|R1 R0|G3 G2 G1 G0|B1 B0|
+--+--+--+--+--+--+--+--+
```
"""
struct RGB242 <: AbstractStorageRGB{N0f8, 2, 4, 2, :no}
    color::UInt8
    RGB242(color::UInt8, ::Val{:no}) = new(color)
end

# TransparentRGB

struct ARGB4444LE <: AlphaColor{XRGB444LE, N0f8, 4}
    color::UInt16
    ARGB4444LE(color::UInt16, ::Val{:no}) = new(color)
end

struct ARGB4444BE <: AlphaColor{XRGB444BE, N0f8, 4}
    color::UInt16
    ARGB4444BE(color::UInt16, ::Val{:no}) = new(color)
end

struct ARGB1555LE <: AlphaColor{RGB555LE, N0f8, 4}
    color::UInt16
    ARGB1555LE(color::UInt16, ::Val{:no}) = new(color)
end

struct ARGB1555BE <: AlphaColor{RGB555BE, N0f8, 4}
    color::UInt16
    ARGB1555BE(color::UInt16, ::Val{:no}) = new(color)
end

# Type Aliases

const RGB565 = BE ? RGB565BE : RGB565LE
const RGB555 = BE ? RGB555BE : RGB555LE
const XRGB444 = BE ? XRGB444BE : XRGB444LE
const RGBX444 = BE ? RGBX444BE : RGBX444LE
const ARGB4444 = BE ? ARGB4444BE : ARGB4444LE
const ARGB1555 = BE ? ARGB1555BE : ARGB1555LE

# Constructors

(::Type{C})(color::UInt16, ::Val{:le}) where {C <: AbstractRGB16{:le}} =
    C(htol(color), Val(:no))
(::Type{C})(color::UInt16, ::Val{:be}) where {C <: AbstractRGB16{:be}} =
    C(hton(color), Val(:no))

# Here we quantize the input RGB values, but there will be another quantization
# error in the conversion to `N0f8`. It is possible to reduce the error by using
# unequally spaced steps, but for simplicity, we will apply the equally spaced
# steps here.
function (::Type{C})(r::Real, g::Real, b::Real) where {C <: Union{AbstractRGB16,
                                                                  RGB332, RGB242}}
    # For simplicity, the 0.5eps margins are not accepted.
    (min(r, g, b) >= 0) & (max(r, g, b) <= 1) || throw_rgb_error(C, r, g, b)
    T = rawtype(C)
    ri = unsafe_trunc(T, round(Int32, float(r) * saturatedr(C)))
    gi = unsafe_trunc(T, round(Int32, float(g) * saturatedg(C)))
    bi = unsafe_trunc(T, round(Int32, float(b) * saturatedb(C)))
    C(ri << nshiftr(C) + gi << nshiftg(C) + bi << nshiftb(C), Val(endian(C)))
end

function (::Type{C})(r::Integer, g::Integer, b::Integer) where {C <: AbstractRGB16}
    (r | g | b) >>> 0x1 == 0 || throw_rgb_error(C, r, g, b)
    ri = ifelse(r % Bool, maskr(C), UInt16(0))
    gi = ifelse(g % Bool, maskg(C), UInt16(0))
    bi = ifelse(b % Bool, maskb(C), UInt16(0))
    C(ri + gi + bi, Val(endian(C)))
end

function (::Type{C})(r::N0f8, g::N0f8, b::N0f8) where {C <: Union{RGB565LE, RGB565BE}}
    r8, g8, b8 = reinterpret(r), reinterpret(g), reinterpret(b)
    ri = (widemul(r8, 0xF9) + 0x400) & maskr(C)
    gi = (widemul(g8, 0xFD) + 0x200) >> nbitsfr(C)
    bi = (widemul(b8, 0xF9) + 0x400) >> nshiftr(C)
    C(ri + (gi & maskg(C)) + bi, Val(endian(C)))
end

function (::Type{C})(r::N0f8, g::N0f8, b::N0f8) where {C <: Union{RGB555LE, RGB555BE}}
    r8, g8, b8 = reinterpret(r), reinterpret(g), reinterpret(b)
    ri = (widemul(r8, 0xF9) + 0x400) >> 0x1
    gi = (widemul(g8, 0xF9) + 0x400) >> (nbitsfr(C) + 1)
    bi = (widemul(b8, 0xF9) + 0x400) >> (nshiftr(C) + 1)
    C((ri &  maskr(C)) + (gi & maskg(C)) + bi, Val(endian(C)))
end

function (::Type{C})(r::N0f8, g::N0f8, b::N0f8) where {C <: Union{XRGB444LE, XRGB444BE,
                                                                  RGBX444LE, RGBX444BE}}
    r8, g8, b8 = reinterpret(r), reinterpret(g), reinterpret(b)
    ri = (UInt16(r8) & saturatedr(C)) << nshiftr(C)
    gi = (UInt16(g8) & saturatedg(C)) << nshiftg(C)
    bi = (UInt16(b8) & saturatedb(C)) << nshiftb(C)
    C(ri + gi + bi, Val(endian(C)))
end

@noinline function throw_rgb_error(@nospecialize(C),
                                   @nospecialize(r), @nospecialize(g), @nospecialize(b))
    throw(ArgumentError(
        "Each RGB value must be in the interval [0,1]. $C cannot represent ($r, $g, $b)."))
end


text:
using VEM
using Documenter

DocMeta.setdocmeta!(VEM, :DocTestSetup, :(using VEM); recursive=true)

makedocs(;
    modules=[VEM],
    authors="Carlos Fernando Baptista",
    repo="https://github.com/EvolvingDynamics/VEM.jl/blob/{commit}{path}#{line}",
    sitename="VEM.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://EvolvingDynamics.github.io/VEM.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/EvolvingDynamics/VEM.jl",
    devbranch="master",
)


text:
import MKL_jll

struct VMLAccuracy
    mode::UInt
end

const VML_LA = VMLAccuracy(0x00000001)
const VML_HA = VMLAccuracy(0x00000002)
const VML_EP = VMLAccuracy(0x00000003)

Base.show(io::IO, m::VMLAccuracy) = print(io, m == VML_LA ? "VML_LA" :
                                              m == VML_HA ? "VML_HA" : "VML_EP")

vml_get_mode() = ccall((:vmlGetMode, MKL_jll.libmkl_rt), Cuint, ())
vml_set_mode(mode::Integer) = (ccall((:vmlSetMode, MKL_jll.libmkl_rt), Cuint, (UInt,), mode); nothing)

vml_set_accuracy(m::VMLAccuracy) = vml_set_mode((vml_get_mode() & ~0x03) | m.mode)
vml_get_accuracy() = VMLAccuracy(vml_get_mode() & 0x3)

function vml_check_error()
    vml_error = ccall((:vmlClearErrStatus, MKL_jll.libmkl_rt), Cint, ())
    if vml_error != 0
        if vml_error == 1
            throw(DomainError(-1, "This function does not support arguments outside its domain"))
        elseif vml_error == 2 || vml_error == 3 || vml_error == 4
            # Singularity, overflow, or underflow
            # I don't think Base throws on these
        elseif vml_error == 1000
            warn("IntelVectorMath does not support $(vml_get_accuracy); lower accuracy used instead")
        else
            error("an unexpected error occurred in IntelVectorMath ($vml_error)")
        end
    end
end

function vml_prefix(t::DataType)
    if t == Float32
        return "vs"
    elseif t == Float64
        return "vd"
    elseif t == Complex{Float32}
        return "vc"
    elseif t == Complex{Float64}
        return "vz"
    end
    error("unknown type $t")
end

function def_unary_op(tin, tout, jlname, jlname!, mklname;
        vmltype = tin)
    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(vmltype))$mklname"))
    exports = Symbol[]
    (@isdefined jlname) || push!(exports, jlname)
    (@isdefined jlname!) || push!(exports, jlname!)
    @eval begin
        function ($jlname!)(out::Array{$tout}, A::Array{$tin})
            size(out) == size(A) || throw(DimensionMismatch())
            ccall(($mklfn, MKL_jll.libmkl_rt), Nothing, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, out)
            vml_check_error()
            return out
        end
        $(if tin == tout
            quote
                function $(jlname!)(A::Array{$tin})
                    ccall(($mklfn, MKL_jll.libmkl_rt), Nothing, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, A)
                    vml_check_error()
                    return A
                end
            end
        end)
        function ($jlname)(A::Array{$tin})
            out = similar(A, $tout)
            ccall(($mklfn, MKL_jll.libmkl_rt), Nothing, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, out)
            vml_check_error()
            return out
        end
        $(isempty(exports) ? nothing : Expr(:export, exports...))
    end
end

function def_binary_op(tin, tout, jlname, jlname!, mklname, broadcast)
    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(tin))$mklname"))
    exports = Symbol[]
    (@isdefined jlname) || push!(exports, jlname)
    (@isdefined jlname!) || push!(exports, jlname!)
    @eval begin
        $(isempty(exports) ? nothing : Expr(:export, exports...))
        function ($jlname!)(out::Array{$tout}, A::Array{$tin}, B::Array{$tin})
            size(out) == size(A) == size(B) || throw(DimensionMismatch("Input arrays and output array need to have the same size"))
            ccall(($mklfn, MKL_jll.libmkl_rt), Nothing, (Int, Ptr{$tin}, Ptr{$tin}, Ptr{$tout}), length(A), A, B, out)
            vml_check_error()
            return out
        end
        function ($jlname)(A::Array{$tout}, B::Array{$tin})
            size(A) == size(B) || throw(DimensionMismatch("Input arrays need to have the same size"))
            out = similar(A)
            ccall(($mklfn, MKL_jll.libmkl_rt), Nothing, (Int, Ptr{$tin}, Ptr{$tin}, Ptr{$tout}), length(A), A, B, out)
            vml_check_error()
            return out
        end
    end
end


text:
using GalerkinSparseGrids
using Test
using HCubature
using StaticArrays
using SparseArrays

#---------------------------------------
# Testing regular hier DG reconstruction
#---------------------------------------

@info "Testing full discontinuous Galerkin (DG) reconstruction 1-D... "
@testset "hier_DG.jl" begin
    @testset "Regular hier DG reconstruction" begin
        for k in 1:5
            for l in 1:6
                dict = coeffs_DG(1, k, l, x->sin(4*x[1]); scheme="full")
                err = x->(reconstruct_DG(dict, [x])-sin(4*x))^2
                @test hquadrature(err, 0, 1; atol=1.0e-10)[1] < 1/(1<<(l+k-1))
            end
        end

        for k in 1:5
            for l in 1:6
                dict = coeffs_DG(1, k, l, x->sin(4*x[1]); scheme="sparse")
                err = x->(reconstruct_DG(dict, [x])-sin(4*x))^2
                @test hquadrature(err, 0, 1; atol=1.0e-10)[1] < 1/(1<<(l+k-1))
            end
        end
    end

    # Sparsification should be trivial in 1D
    @info "Testing that both reconstructions are equivalent... "
    @testset "Equivalent Reconstructions" begin
        for k in 1:5
            for l in 1:6
                full_coeffs = coeffs_DG(1, k, l, x->sin(4*x[1]); scheme="full")
                sparse_coeffs = coeffs_DG(1, k, l, x->sin(4*x[1]); scheme="full")
                diff_full = x->(reconstruct_DG(full_coeffs, [x])-sin(4*x))^2
                diff_sparse = x->(reconstruct_DG(sparse_coeffs, [x])-sin(4*x))^2
                diff_both = x->(reconstruct_DG(full_coeffs, [x])-reconstruct_DG(sparse_coeffs, [x]))^2

                @test abs(hquadrature(diff_full,0,1; atol=1.0e-9)[1]-
                          hquadrature(diff_sparse,0,1; atol=1.0e-9)[1]) < 1.0e-14
                @test hquadrature(diff_both,0,1; atol=1.0e-9)[1] < 1.0e-14
            end
        end
    end

# Testing sparse DG reconstruction in multidimensional space:

    @info "Testing sparse DG basis reconstruction 2-D..."
    @testset "Sparse DG reconstruction in multidimensional space" begin
        for k in 1:5
            for l in 1:6
                dict = coeffs_DG(2, k, l, x->sin(4*x[1]+x[2]))
                err = x->(reconstruct_DG(dict, [x[1],x[2]])-sin(4*x[1]+x[2]))^2
                @test hcubature(err, [0,0], [1,1]; atol=1.0e-10, maxevals=500)[1] < 1/(1<<(l+k-2))
            end
        end
    end
end


text:
macro ereplace(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    error("""\n
        The syntax is:
            ereplace <expression> [if] <expression>
        or
            bysort <varlist> (<varlist>): ereplace <expression> [if] <expression>
    """)
end

# short form. This is a copy of the above for `erep`.
macro erep(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@ereplace($t, $by, $sort, $arguments, $filter, $use, $options)
        end
    )
end

macro ereplace(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr, Nothing},
    use::Nothing,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@replace_byvec!($t, $by, $sort, $arguments, $filter, $options )
        end
    )
end

# implementation
macro replace_byvec!(t::Symbol,
    varlist_by::Union{Vector{Symbol}, Nothing},
    varlist_sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr,Nothing},
    options::Union{Dict{String,Any}, Nothing})

    # assert that `arguments` is an assignment
    (arguments.head == :(.=)) && error("`erep` expects a vector-wise assignment operation, e.g. `:x = :y + :z`. Do not broadcast the assignment operator.")
    (arguments.head == :(=)) || error("`erep` expects an assignment operation, e.g. :x = :y + :z")

    # construct the QuoteNode
    if isexpr(arguments.args[1])
        # this needs to take the form :x::Type
        ( arguments.args[1].head == Symbol("::") ) && error("`ereplace` does not allow explicit type declarations`")
        assigned_var_qn = arguments.args[1].args[1]::QuoteNode
    elseif isa(arguments.args[1], QuoteNode)
        assigned_var_qn = arguments.args[1]::QuoteNode
    end

    # make sure transformation is an Expr
    if isexpr(arguments.args[2])
        transformation = arguments.args[2]
    elseif isa(arguments.args[2], QuoteNode)
        transformation = Expr(:call, identity, arguments.args[2])
    else
        transformation = Expr(:call, identity, arguments.args[2])
    end

    return esc(
        quote
            # check variable is present
            ($(assigned_var_qn) ∈ propertynames($t)) || error("Variable $($(assigned_var_qn)) not present in DataFrame.")
            # do the assignment
            Douglass.@replace_byvec_allobs!($t, $varlist_by, $varlist_sort, $assigned_var_qn, $transformation, $filter)
        end
    )
end

macro replace_byvec_allobs!(t::Symbol,
    varlist_by::Vector{Symbol},
    varlist_sort::Union{Vector{Symbol}, Nothing},
    assigned_var_qn::QuoteNode,
    transformation,
    filter::Expr)

    # go through the expression tree and pick up QuoteNodes.
    qn_vec = Vector{QuoteNode}()
    transformation_converted = deepcopy(transformation)
    filter_converted = deepcopy(filter)

    # replace all QuoteNode's by args[x]
    replace_QuoteNodes!(transformation_converted, :args, qn_vec)
    # and the same for the filter condition
    replace_QuoteNodes!(filter_converted, :args, qn_vec)
    # add $assigned_var_qn if not already present
    idx_assigned_var_qn = 0
    idx = findfirst(isequal(assigned_var_qn), qn_vec)
    if !isnothing(idx)
        # already exists in qn_vec
        idx_assigned_var_qn = idx
    else
        push!(qn_vec, assigned_var_qn)
        idx_assigned_var_qn = length(qn_vec)
    end
    # add an index subscript to all arg's: arg[#] -> arg[#][index]
    # so that we can condition on a filter
    push_index!(transformation_converted, :indices, :args)

    # a vector of the Symbol's for each QuoteNode
    symbol_vec = [qn_vec[i].value for i=1:length(qn_vec)]

    return esc(
        quote
            # sort, if we need to, (first by-variables, then sort-variables)
            if !isnothing($(varlist_sort)) && !isempty($(varlist_sort))
                sort!($t, vcat($varlist_by, $varlist_sort))
            end

            # define a barrier function, which will be compiled and type stable
            # see Ch. 11 of Bogumil Kaminski's DataFrames tutorial
            local function my_f(args...)
                # indices are the indices of the input columns to be used.
                indices = BitArray($filter_converted)
                # the following is a wrapper for args[idx_assigned_var_qn][indices] = $transformation_converted
                # but if $transformation_converted expands into a scalar, it's broadcast
                return Douglass.assign_helper_rep!(args[$idx_assigned_var_qn],indices,$transformation_converted)
            end
            transform!(groupby($t, $varlist_by, sort=false, skipmissing = true ), $(symbol_vec) => my_f => $(assigned_var_qn))
        end

    )
end
# `by` but no `if`
macro replace_byvec_allobs!(t::Symbol,
    varlist_by::Vector{Symbol},
    varlist_sort::Union{Vector{Symbol}, Nothing},
    assigned_var_qn::QuoteNode,
    transformation,
    filter::Nothing)

    # go through the expression tree and pick up QuoteNodes.
    qn_vec = Vector{QuoteNode}()
    transformation_converted = deepcopy(transformation)
    filter_converted = deepcopy(filter)

    # replace all QuoteNode's by args[x]
    replace_QuoteNodes!(transformation_converted, :args, qn_vec)

    # a vector of the Symbol's for each QuoteNode
    symbol_vec = [qn_vec[i].value for i=1:length(qn_vec)]

    return esc(
        quote
            # sort, if we need to, (first by-variables, then sort-variables)
            if !isnothing($(varlist_sort)) && !isempty($(varlist_sort))
                sort!($t, vcat($varlist_by, $varlist_sort))
            end

            # define a barrier function, which will be compiled and type stable
            # see Ch. 11 of Bogumil Kaminski's DataFrames tutorial
            local function my_f(args...)
                return $transformation_converted
            end
            transform!(groupby($t, $varlist_by, sort=false, skipmissing = true ), $(symbol_vec) => my_f => $(assigned_var_qn))
        end

    )
end

# version without groups but with if
macro replace_byvec_allobs!(t::Symbol,
    varlist_by::Nothing,
    varlist_sort::Union{Vector{Symbol}, Nothing},
    assigned_var_qn::QuoteNode,
    transformation,
    filter::Expr)

    # go through the expression tree and pick up QuoteNodes.
    qn_vec = Vector{QuoteNode}()
    transformation_converted = deepcopy(transformation)
    filter_converted = deepcopy(filter)

    # replace all QuoteNode's by args[x]
    replace_QuoteNodes!(transformation_converted, :args, qn_vec)
    # and the same for the filter condition
    replace_QuoteNodes!(filter_converted, :args, qn_vec)
    # add $assigned_var_qn if not already present
    idx_assigned_var_qn = 0
    idx = findfirst(isequal(assigned_var_qn), qn_vec)
    if !isnothing(idx)
        # already exists in qn_vec
        idx_assigned_var_qn = idx
    else
        push!(qn_vec, assigned_var_qn)
        idx_assigned_var_qn = length(qn_vec)
    end
    # add an index subscript to all arg's: arg[#] -> arg[#][index]
    # so that we can condition on a filter
    push_index!(transformation_converted, :indices, :args)

    # a vector of the Symbol's for each QuoteNode
    symbol_vec = [qn_vec[i].value for i=1:length(qn_vec)]

    return esc(
        quote
            # define a barrier function, which will be compiled and type stable
            # see Ch. 11 of Bogumil Kaminski's DataFrames tutorial
            local function my_f(args...)
                # indices are the indices of the input columns to be used.
                indices = BitArray($filter_converted)
                # the following is a wrapper for a[indices] = $transformation_converted
                # but if $transformation_converted expands into a scalar, it's broadcast
                return Douglass.assign_helper_rep!(args[$idx_assigned_var_qn],indices,$transformation_converted)
            end
            transform!($t, $(symbol_vec) => my_f => $(assigned_var_qn))
        end
    )
end
# version without `by` and `if`
macro replace_byvec_allobs!(t::Symbol,
    varlist_by::Nothing,
    varlist_sort::Union{Vector{Symbol}, Nothing},
    assigned_var_qn::QuoteNode,
    transformation,
    filter::Nothing)

    # go through the expression tree and pick up QuoteNodes.
    qn_vec = Vector{QuoteNode}()
    transformation_converted = deepcopy(transformation)

    # replace all QuoteNode's by args[x]
    replace_QuoteNodes!(transformation_converted, :args, qn_vec)

    # a vector of the Symbol's for each QuoteNode
    symbol_vec = [qn_vec[i].value for i=1:length(qn_vec)]

    return esc(
        quote
            # define a barrier function, which will be compiled and type stable
            # see Ch. 11 of Bogumil Kaminski's DataFrames tutorial
            local function my_f(args...)
                return $transformation_converted
            end
            transform!($t, $(symbol_vec) => my_f => $(assigned_var_qn))
        end
    )
end


text:
function test()
    # some testing
    stepc = 5500
    epsilon = 0.4
    gamma = 0.9
    w = map(KnetArray, initialize_weights(actionc_accel*actionc_steer))
    state, env = initialize_simple_road();
    tree = build_tree(w, env, state, stepc, epsilon,gamma);
    backprop_weighted_q!(tree, gamma)
    @test (stepc+1)==length(tree.nodes)==length(tree.sensors)==length(tree.children_qs)==length(tree.states)==length(tree.accumulated_rewards)==length(tree.dones)
    @test is_root(tree.nodes[1])
    @test is_leaf(tree.nodes[end])
    visitedc = zeros(Int32, length(tree.nodes))
    for nix in length(tree.nodes):-1:1
        node = tree.nodes[nix]
        for ch in values(node.children)
            visitedc[nix] += visitedc[ch.id]
        end
        visitedc[nix] += 1
    end
    @test visitedc[1] == stepc+1
    for (ix,n) in enumerate(tree.nodes)
        tree.dones[n.id] && @test length(n.children)==0 || all_children_done(tree, n)
        if !is_leaf(n)
            @test visitedc[n.id] == 1+sum([visitedc[c.id] for c in values(n.children)])
        end
        if is_leaf(n)
            @test visitedc[n.id] == 1
            tree.dones[n.id] && @test abs(sum(tree.children_qs[n.id])) < 7.0e-5
            parent = get(n.parent)
            mn = Float32(mean(tree.children_qs[n.id]))
            reward = tree.accumulated_rewards[n.id] - tree.accumulated_rewards[parent.id]
            q = reward + gamma*mn
            @test abs(tree.children_qs[parent.id][n.action_ix] - q) < 7.0e-5
        elseif !is_root(n)
            parent = get(n.parent)
            q::Float32 = 0.0
            for (i,c) in n.children
                q += visitedc[c.id] * tree.children_qs[n.id][i]
            end
            q /= visitedc[n.id]-1
            q *= gamma
            reward = tree.accumulated_rewards[n.id] - tree.accumulated_rewards[parent.id]
            q += reward
            @test abs(tree.children_qs[parent.id][n.action_ix] - q) < 7.0e-5
        end
    end
    println("[V] all tests passed!")
end


text:
module SciPyOptimize

export fmin_powell, fmin_de

using PyCall

function __init__()
    py"""
    import warnings
    import numpy as np
    from scipy.optimize import minimize, differential_evolution
    from typing import Callable

    warnings.filterwarnings('ignore')

    def modified_powell(
            objective: Callable,
            n_gene: int,
            population: np.ndarray,
            ip: np.ndarray,
            maxiter: int,
    ) -> np.ndarray:
        lower = np.min(population[:, :n_gene], axis=0)
        upper = np.max(population[:, :n_gene], axis=0)
        direc = np.identity(n_gene) * 0.3 * (upper - lower)
        res = minimize(
            objective,
            population[ip[0], :n_gene],
            method='Powell',
            bounds=tuple(zip(lower, upper)),
            callback=lambda xk: True
            if objective(xk) < objective(population[ip[0], :n_gene])
            else False,
            options={
                #'disp': True,
                'xtol': 0.1,
                'ftol': 1.0,
                'maxiter': maxiter,
                'maxfev': 1000 * n_gene,
                'direc': direc,
            }
        )
        obj_val = objective(res.x)
        if obj_val < objective(population[ip[0], :n_gene]):
            population[ip[0], :n_gene] = res.x
            population[ip[0], -1] = obj_val
        return population

    def de_best2bin(
            objective: Callable,
            n_gene: int,
            population: np.ndarray,
            ip: np.ndarray,
            maxiter: int,
    ) -> np.ndarray:
        lower = np.min(population[:, :n_gene], axis=0)
        upper = np.max(population[:, :n_gene], axis=0)
        res = differential_evolution(
            objective,
            tuple(zip(lower, upper)),
            strategy='best2bin',
            mutation=(0.0, 0.3),
            recombination=0.9,
            maxiter=maxiter,
            popsize=1,
            polish=False,
            init=population[ip, :n_gene],
        )
        obj_val = objective(res.x)
        if obj_val < objective(population[ip[0], :n_gene]):
            population[ip[0], :n_gene] = res.x
            population[ip[0], -1] = obj_val
        return population
    """
end

function fmin_powell(
        objective::Function,
        n_gene::Int,
        population::Matrix{Float64},
        ip::Vector{Int},
        maxiter::Int
)::Matrix{Float64}
    return py"modified_powell"(objective, n_gene, population, ip .- 1, maxiter)
end

function fmin_de(
        objective::Function,
        n_gene::Int,
        population::Matrix{Float64},
        ip::Vector{Int},
        maxiter::Int
)::Matrix{Float64}
    return py"de_best2bin"(objective, n_gene, population, ip .- 1, maxiter)
end

end # module


text:
function huber_loss(ŷ, y)
    δ = 0.5
    error = abs.(ŷ - y)
    cond = error .< δ
    squared_loss = 0.5 * error.^2
    linear_loss = δ .* (error .- δ^2 .* 0.5)

    is_squared = Int.(cond)
    is_linear = Int.(.!cond)

    # why this line doesn't work?
    # Int.(cond) .* squared_loss + Int.(.!cond) .* (10 .* linear_loss)
    #=
    ERROR: LoadError: MethodError: no method matching !(::ForwardDiff.Dual{Nothing,Bool,1})
    Closest candidates are:
        !(!Matched::Missing) at missing.jl:83
        !(!Matched::Bool) at bool.jl:35
        !(!Matched::Function) at operators.jl:853
    =#
    is_squared .* squared_loss + is_linear .* (10 .* linear_loss)
end

function huber_loss_mean(ŷ, y)
    mean(huber_loss(ŷ, y))
end

function mse_rnn(x, y)
    l = Flux.mse(m(x), y)
    Flux.reset!(m)
    return l
end


text:
# TMI solutions were originally saved in mat format. Here convert to NetCDF using
using Revise, MAT, LinearAlgebra, SparseArrays, TMI, DrWatson, NCDatasets, Test

# choose one version
TMIversion = "modern_90x45x33_GH10_GH12"
TMIversion = "modern_180x90x33_GH11_GH12"
TMIversion = "modern_90x45x33_unpub12"
TMIversion = "modern_90x45x33_G14_v2"
TMIversion = "modern_90x45x33_G14"
TMIversion = "LGM_90x45x33_G14"
TMIversion = "LGM_90x45x33_G14A"
TMIversion = "LGM_90x45x33_GPLS1"
TMIversion = "LGM_90x45x33_GPLS2"
TMIversion = "LGM_90x45x33_OG18"

# original NetCDF version
#A, Alu, γ, TMIfile = config(TMIversion)

@time A, Alu, γ, TMIfile, L, B = config_from_mat(TMIversion);
#@time A2, Alu2, γ2, TMIfile2, L2, B2 = config_from_mat(TMIversion2);

config2nc(TMIversion,A,γ,L,B)

filenetcdf = datadir("TMI_"*TMIversion*".nc")
NCDataset(filenetcdf)

# read from NetCDF?
@time A2, Alu2, γ2, TMIfile2, L2, B2 = config_from_nc(TMIversion);

# Do NetCDF and mat-files agree?
@test A2 == A
@test L2 == L
@test B2 == B


text:
# represents a service - i.e. is a service running on a given day?
using Dates

struct Service
    monday::Bool
    tuesday::Bool
    wednesday::Bool
    thursday::Bool
    friday::Bool
    saturday::Bool
    sunday::Bool
    start_date::Date
    end_date::Date
    added_dates::Array{Date}
    removed_dates::Array{Date}
end

function is_service_running(service::Service, date::Date)::Bool
    # first check calendar_dates
    if any(service.added_dates .== date)
        return true
    elseif any(service.removed_dates .== date)
        return false
    else
        dow = Dates.dayofweek(date)
        # first check calendar
        return (
            date >= service.start_date &&
            date <= service.end_date && (
                (dow == 1 && service.monday) ||
                (dow == 2 && service.tuesday) ||
                (dow == 3 && service.wednesday) ||
                (dow == 4 && service.thursday) ||
                (dow == 5 && service.friday) ||
                (dow == 6 && service.saturday) ||
                (dow == 7 && service.sunday)
            )
        )
    end
end


text:
module WK3

using MAT
using Interpolations
using Distributions

# conversions
include("constants.jl")

# 5th order Runge-Kutta solver w/ adaptive stepsize (from Numerical Recipes)
include("solverparams.jl")
include("modelparams.jl")
include("rkck.jl")
include("rkqs.jl")
include("odeint.jl")

# type definitions
include("buildall.jl")

# odes for 3-element Windkessel w/ ventricular pressure source & diode valve
include("wk3odes.jl")

# ventricular elastance function
include("elastancefn.jl")

# data assimilation
include("patdatainterp.jl")
include("builderrors.jl")
include("paramwalk.jl")

# Fisher information matrix
include("fdjac.jl")

export odeint
export rkqs
export rkck
export wk3odes
export elastancefn
export SolverParams
export ModelParams
export buildall
export patdatainterp
export builderrors
export paramwalk!
export fdjac

end


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Translation action" begin
    M = Euclidean(2, 3)
    G = TranslationGroup(2, 3)
    A = TranslationAction(Euclidean(2, 3), G)

    @test repr(A) == "TranslationAction($(repr(M)), $(repr(G)), LeftAction())"
    @test repr(switch_direction(A)) ==
          "TranslationAction($(repr(M)), $(repr(G)), RightAction())"

    types_a = [Matrix{Float64}]

    types_m = [Matrix{Float64}]

    @test g_manifold(A) == M
    @test base_group(A) == G
    @test base_manifold(G) == M

    for (T_A, T_M) in zip(types_a, types_m)
        a_pts =
            convert.(
                T_A,
                [
                    [0.0 1.0 2.0; 2.0 4.0 1.0],
                    [-1.0 0.0 2.0; -2.0 3.0 3.0],
                    [1.0 1.0 2.0; 3.0 2.0 1.0],
                ],
            )
        m_pts =
            convert.(
                T_M,
                [
                    [0.0 1.0 2.0; 2.0 4.0 1.0],
                    [-1.0 0.0 2.0; -2.0 3.0 3.0],
                    [1.0 1.0 2.0; 3.0 2.0 1.0],
                ],
            )
        v_pts =
            convert.(
                T_M,
                [
                    [0.0 1.0 2.0; 2.0 4.0 1.0],
                    [-1.0 0.0 2.0; -2.0 3.0 3.0],
                    [1.0 1.0 2.0; 3.0 2.0 1.0],
                ],
            )

        atol = if eltype(T_M) == Float32
            1e-7
        else
            1e-15
        end
        @test apply_diff(A, a_pts[1], m_pts[1], v_pts[1]) === v_pts[1]
        @test inverse_apply_diff(A, a_pts[1], m_pts[1], v_pts[1]) === v_pts[1]
        test_action(
            A,
            a_pts,
            m_pts,
            v_pts;
            test_optimal_alignment=false,
            test_diff=true,
            atol=atol,
        )
    end
end


text:
module CIM

include("dummy.jl")

include("core.jl")

end


text:
include(joinpath(@__DIR__, "LibASICamera_types.jl"))
include(joinpath(@__DIR__, "LibASICamera_ccalls.jl"))

"""
    ASICamera

The struct which contains information about the camera.
Has fields .info and .control_caps.
"""
struct ASICamera
    info::ASI_CAMERA_INFO
    control_caps::Vector{ASI_CONTROL_CAPS}
end

"""
    allocate_buffer(width::Integer, height::Integer, img_type::ASI_IMG_TYPE)

Allocates an image buffer for the camera to write to.

# Args:
    width: Image width
    height: Image height
    img_type: One of
        -ASI_IMG_RAW8
        -ASI_IMG_Y8
        -ASI_IMG_RAW16
        -ASI_IMG_RAW24

# Returns:
    A zero-initialized array of the appropriate shape.

# Throws:
    ASIWrapperError if an unsupported image type is given.
"""
function allocate_buffer(width::Integer, height::Integer, img_type::ASI_IMG_TYPE)
    if img_type==ASI_IMG_RAW8 || img_type==ASI_IMG_Y8
        return zeros(UInt8, width, height)
    elseif img_type==ASI_IMG_RAW16
        return zeros(UInt16, width, height)
    elseif img_type==ASI_IMG_RAW24
        return zeros(UInt8, width, height, 3)
    else
        throw(ASIWrapperError("Image type $img_type not implemented."))
    end
end

# This allows to use an ellipsis: allocate_buffer(get_roi_format(cam)...)
allocate_buffer(width::Integer, height::Integer, unused, img_type::ASI_IMG_TYPE) =
    allocate_buffer(width, height, img_type)

"""
    get_num_connected_cameras()

This function returns the count of connected cameras and should be called first.
"""
get_num_connected_cameras() = ASIGetNumOfConnectedCameras()

"""
    get_camera_property(id::Integer)

Fetches the camera properties for a given ID.

# Args:
    id: Camera id

# Returns:
    ASI_CAMERA_INFO object

# Throws:
    ASIError in case of failure
"""
function get_camera_property(id::Integer)
    camera_info = _ASI_CAMERA_INFO()
    err = ASIGetCameraProperty(camera_info, id)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
    return ASI_CAMERA_INFO(camera_info) # convert to human-readable type
end

get_camera_property(cam::ASICamera) = get_camera_property(cam.info.CameraID)

"""
    open_camera(id::Integer)

Opens the ASI camera connection.

Open the camera before any interaction with the camera,
this will not affect the camera which is capturing.
Then you must call init_camera() to perform any actions.

# Args:
    id: Camera ID

# Throws:
	ASIError
"""
function open_camera(id::Integer)
    err = ASIOpenCamera(id)
    if err != ASI_SUCCESS
        print("[ERROR] Could not connect to camera. Make sure you can access the camera without being root by calling
        \"sudo install asi.rules /lib/udev/rules.d\" from the lib subdir of the SDK
        and then relogging / rebooting.")
        throw(ASIError(err))
    end
end

open_camera(cam::ASICamera) = open_camera(cam.info.CameraID)

"""
    init_camera(id::Integer)

Initialize the camera. Needs to be called before capturing any data.

# Args:
    id: Camera id

# Throws:
    ASIError
"""
function init_camera(id::Integer)
    err = ASIInitCamera(id)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

init_camera(cam::ASICamera) = init_camera(cam.info.CameraID)

"""
    close_camera(id::Integer)

Closes the ASI camera connection.

# Args:
    id: Camera id

# Throws:
	ASIError
"""
function close_camera(id::Integer)
    err = ASICloseCamera(id)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

close_camera(cam::ASICamera) = close_camera(cam.info.CameraID)

"""
    get_control_caps(id::Integer)

Returns the control properties available for this camera.

The camera needs to be open.

# Args:
    id: Camera id

# Returns:
    Vector of ASI_CONTROL_CAPS structs.

# Throws:
    ASIError
"""
function get_control_caps(id::Integer)
    num_controls = Ref{Cint}(0)

    err = ASIGetNumOfControls(id, num_controls)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end

    control_caps = Vector{ASI_CONTROL_CAPS}(undef, num_controls.x)

    for i in 1:num_controls.x
        cap = _ASI_CONTROL_CAPS()
        err = ASIGetControlCaps(id, i-1, cap)
        if err != ASI_SUCCESS
            throw(ASIError(err))
        end
        control_caps[i] = ASI_CONTROL_CAPS(cap) # convert to human-readable type
    end

    return control_caps
end

get_control_caps(cam::ASICamera) = get_control_caps(cam.info.CameraID)

"""
    get_connected_devices()

Returns a list of the connected devices.
"""
function get_connected_devices()
    num_cameras = get_num_connected_cameras()
    if num_cameras == 0
        throw(ASIWrapperError("No cameras found."))
    end

    cameras = Vector{ASICamera}(undef, num_cameras)
    for i in 0:num_cameras-1
        open_camera(i)
        init_camera(i)
        info = get_camera_property(i)
        control_caps = get_control_caps(i)
        cameras[i+1] = ASICamera(info, control_caps)
    end
    return cameras
end

"""
    get_control_value(id::Integer, control_type::ASI_CONTROL_TYPE)

Fetches the current setting of the control value, e.g. exposure or gain.

# Args:
    id: Camera id
    control_type: The control type to fetch, e.g. exposure or gain.

# Returns:
    A tuple (value, is_auto)

# Throws:
    ASIError
"""
function get_control_value(id::Integer, control_type::ASI_CONTROL_TYPE)
    value = Ref{Clong}(0)
    auto  = Ref{ASI_BOOL}(ASI_FALSE)
    err = ASIGetControlValue(id, control_type, value, auto)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
    return (value[], Bool(auto[]))
end

get_control_value(cam::ASICamera, control_type::ASI_CONTROL_TYPE) = get_control_value(cam.info.CameraID, control_type)

get_temperature(id) = get_control_value(id, ASI_TEMPERATURE)[1] * 0.1

"""
    set_control_value(id::Integer, control_type::ASI_CONTROL_TYPE, value, auto::Bool=false)

Sets a control (e.g. exposure) to the given value. Automatically sets the
minimum or maximum if the given value is out of bounds.

# Args:
    id: Camera id
    control_type: The control type to set, e.g. exposure or gain.
    value: The value to which the control is set.
    auto: Whether or not the control should be automatically set.
        Check if this is supported for the given control beforehand.

# Throws:
    ASIError
"""
function set_control_value(id::Integer, control_type::ASI_CONTROL_TYPE, value, auto::Bool=false)
    err = ASISetControlValue(id, control_type, value, auto)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

set_control_value(cam::ASICamera, control_type::ASI_CONTROL_TYPE, value, auto::Bool=false) = set_control_value(cam.info.CameraID, control_type, value, auto)

set_gain(id, gain, auto=false)              = set_control_value(id, ASI_GAIN, gain, auto)
set_exposure(id, exposure_μs, auto=false)   = set_control_value(id, ASI_EXPOSURE, exposure_μs, auto)
set_gamma(id, gamma, auto=false)            = set_control_value(id, ASI_GAMMA, gamma, auto)
set_bandwidth(id, bandwidth, auto=false)    = set_control_value(id, ASI_BANDWIDTHOVERLOAD, bandwidth, auto)
set_flip(id, flip)                              = set_control_value(id, ASI_FLIP, flip)
set_autoexp_max_gain(id, val)                   = set_control_value(id, ASI_AUTO_MAX_GAIN, val)
set_autoexp_max_exp(id, exposure_ms)            = set_control_value(id, ASI_AUTO_MAX_EXP, exposure_ms)
set_autoexp_target_brightness(id, brightness)   = set_control_value(id, ASI_AUTO_MAX_BRIGHTNESS, brightness)
set_highspeed_mode(id, active)                  = set_control_value(id, ASI_HIGH_SPEED_MODE, active)
# ...

"""
    get_status(id::Integer)

Returns the status of all camera parameters.
"""
function get_status(id::Integer)
    control_caps = get_control_caps(id)
    ret = []
    for t in control_caps
        push!(ret, (t.Name, get_control_value(id, t.ControlType)[1]))
    end
    return ret
end

get_status(cam::ASICamera) = get_status(cam.info.CameraID)

"""
    set_roi_format(id::Integer, width, height, binning, img_type::ASI_IMG_TYPE)

Sets the region of interest (roi). Do so before capturing.

The width and height are the values *after* binning, i.e. you need to set the
width to 640 and the height to 480 if you want to run at 640x480 @ BIN2.

ASI120's data size must be a multiple of 1024 which means width*height%1024==0.

# Args:
    id: Camera id
    width: ROI width
    height: ROI height
    binning: The binning mode; 2 means to read out 2x2 pixels together. Check
        which binning values are supported in the ASI_CAMERA_INFO struct of the
        camera struct or by calling get_camera_property(id).

# Throws:
    ASIError
"""
function set_roi_format(id::Integer, width, height, binning, img_type::ASI_IMG_TYPE)
    if width%8 != 0
        throw("Width must be a multiple of 8.")
    end
    if height%2 != 0
        throw("Height must be a multiple of 2.")
    end
    if (width*height)%1024 != 0 && get_camera_property(id).Name in ["ZWO ASI 120MM", "ZWO ASI 120MC"]
        throw("Width times height must be a multiple of 2.")
    end
    err = ASISetROIFormat(id, width, height, binning, img_type)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

set_roi_format(cam::ASICamera, width, height, binning, img_type::ASI_IMG_TYPE) = set_roi_format(cam.info.CameraID, width, height, binning, img_type)

"""
    get_roi_format(id::Integer)

Fetches the current region of interest settings.
"""
function get_roi_format(id::Integer)
    width = Ref{Cint}(0)
    height = Ref{Cint}(0)
    binning = Ref{Cint}(false)
    img_type = Ref{ASI_IMG_TYPE}(ASI_IMG_RAW8)
    err = ASIGetROIFormat(id, width, height, binning, img_type)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
    return (width.x, height.x, binning.x, img_type.x)
end

get_roi_format(cam::ASICamera) = get_roi_format(cam.info.CameraID)

"""
    set_roi_start(id::Integer, startx, starty)

Sets the position of the top-left corner of the region of interest.

You can call this while the camera is streaming to move the ROI. By default,
the ROI will be centered. In binned mode, the start values are relative to the
binned sensor size.

# Throws:
    ASIError
"""
function set_roi_start(id::Integer, startx, starty)
    err = ASISetStartPos(id, startx, starty)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

set_roi_start(cam::ASICamera, startx, starty) = set_roi_start(cam.info.CameraID, startx, starty)

"""
    get_roi_start(id::Integer)

Returns the region of interest start position (start_x, start_y).
"""
function get_roi_start(id::Integer)
    startx = Ref{Cint}(0)
    starty = Ref{Cint}(0)
    err = ASIGGetStartPos(id, startx, starty)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
    return (startx[], starty[])
end

get_roi_start(cam::ASICamera) = get_roi_start(cam.info.CameraID)

"""
    get_dropped_frames(id::Integer)

Returns the number of dropped frames.

Frames are dropped when the USB is bandwidth is low or the
harddisk write speed is slow. The count is reset to 0 after capturing stops.
"""
function get_dropped_frames(id::Integer)
    dropped_frames = Ref{Cint}(0)
    err = ASIGetDroppedFrames(id, dropped_frames)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

get_dropped_frames(cam::ASICamera) = get_dropped_frames(cam.info.CameraID)

"""
"""
function enable_dark_subtract(id::Integer, path)
    err = ASIEnableDarkSubtract(id, path)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

enable_dark_subtract(cam::ASICamera) = enable_dark_subtract(cam.info.CameraID)

"""
"""
function disable_dark_subtract(id::Integer)
    err = ASIDisableDarkSubtract(id)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

disable_dark_subtract(cam::ASICamera) = disable_dark_subtract(cam.info.CameraID)

"""
    start_video(id::Integer)

Start video capture.

# Throws:
    ASIError
"""
function start_video(id::Integer)
    err = ASIStartVideoCapture(id)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

start_video(cam::ASICamera) = start_video(cam.info.CameraID)

"""
    stop_video(id::Integer)

Stops a running video capture.

# Throws:
    ASIError
"""
function stop_video(id::Integer)
    err = ASIStopVideoCapture(id)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

stop_video(cam::ASICamera) = stop_video(cam.info.CameraID)

"""
    get_video_data!(id::Integer, buffer, timeout_ms)

Writes a video frame to the given buffer. Make sure the buffer is large
enough to fit the frame. Call this as fast as possible in a loop and check
whether the return value equals ASI_SUCCESS.

# Args:
    id: Camera id
    buffer: A buffer to write the video frame to.
    timeout_ms: Time to wait for a frame. Recommendation: 2 * exposure_μs + 500 ms <- inconsistent units?!

# Returns:
    An ASI_ERROR_CODE, which should be ASI_SUCCESS.
"""
function get_video_data!(id::Integer, buffer, timeout_ms=500)
    return ASIGetVideoData(id, buffer, sizeof(buffer), Int32(round(timeout_ms)))
end

get_video_data!(cam::ASICamera, buffer, timeout_ms) = get_video_data!(cam.info.CameraID, buffer, timeout_ms)

"""
    pulse_guide_on(id::Integer, direction::ASI_GUIDE_DIRECTION)

Activates the pulse guide in the given direction.

# Args:
    id: Camera id
    direction: Guiding direction; call 'instances(ASI_GUIDE_DIRECTION)'

# Throws:
    ASIError
"""
function pulse_guide_on(id::Integer, direction::ASI_GUIDE_DIRECTION)
    err = ASIPulseGuideOn(id, direction)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

pulse_guide_on(cam::ASICamera) = pulse_guide_on(cam.info.CameraID)

"""
    pulse_guide_off(id::Integer, direction::ASI_GUIDE_DIRECTION)

Deactivates the pulse guide in the given direction.

# Args:
    id: Camera id
    direction: Guiding direction; call 'instances(ASI_GUIDE_DIRECTION)' for options.

# Throws:
    ASIError
"""
function pulse_guide_off(id::Integer, direction::ASI_GUIDE_DIRECTION)
    err = ASIPulseGuideOff(id, direction)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

pulse_guide_off(cam::ASICamera) = pulse_guide_off(cam.info.CameraID)

"""
    start_exposure(id::Integer, is_dark=false)

Starts an exposure. All relevant parameters (exposure time, gain) have to be
set beforehand by calling set_control_value(...) or e.g. set_gain(...).
"""
function start_exposure(id::Integer, is_dark=false)
    err = ASIStartExposure(id, is_dark)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

start_exposure(cam::ASICamera, is_dark=false) = start_exposure(cam.info.CameraID, is_dark)

"""
    stop_exposure(id::Integer)

Stops an ongoing exposure.
"""
function stop_exposure(id::Integer)
    err = ASIStopExposure(id)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

stop_exposure(cam::ASICamera) = stop_exposure(cam.info.CameraID)

"""
    get_exp_status(id::Integer)

Returns the status of an ongoing exposure.
See 'instances(ASI_EXP_STATUS)'.

# Throws:
    ASIError
"""
function get_exp_status(id::Integer)
    status = Ref{ASI_EXPOSURE_STATUS}(ASI_EXP_IDLE)
    err = ASIGetExpStatus(id, status)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
    return status[]
end

get_exp_status(cam::ASICamera) = get_exp_status(cam.info.CameraID)

"""
    get_data_after_exp!(id::Integer, buffer)

Fetches the data after a successful exposure and writes it into buffer.
"""
function get_data_after_exp!(id::Integer, buffer)
    err = ASIGetDataAfterExp(id, buffer, sizeof(buffer))
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

get_data_after_exp!(cam::ASICamera, buffer) = get_data_after_exp!(cam.info.CameraID, buffer)

"""
    get_id(id::Integer)

Returns the camera id stored in flash, only available for USB3 cameras.
"""
function get_id(id::Integer)
    camera_id = ASI_ID()
    err = ASIGetID(id, camera_id)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

get_id(cam::ASICamera) = get_id(cam.info.CameraID)

"""
    capture_still(id::Integer)

Captures a still image. You have to set gain, exposure etc. beforehand using
set_control_value(...).

# Returns:
    An array containing the image.

# Throws:
    ASIWrapperError if the image format is not supported by the wrapper, or
    ASIError in other unfortunate cases
"""
function capture_still(id::Integer)
    exposure_μs, _ = get_control_value(id, ASI_EXPOSURE)
    width, height, binning, img_type = get_roi_format(id)

    buffer = allocate_buffer(width, height, img_type)

    start_exposure(id, false)
    sleep(0.05)
    # option 2: wait almost the entire exposure time instead of polling
    # sleep(exposure_μs/1000*0.9)

    while get_exp_status(id) == ASI_EXP_WORKING
        sleep(0.01)
    end

    if get_exp_status(id) == ASI_EXP_SUCCESS
        get_data_after_exp!(id, buffer)
        return buffer[end:-1:1, :]
    else
        throw("Exposure failed")
    end
end

capture_still(cam::ASICamera) = capture_still(cam.info.CameraID)

"""
    get_gain_offset(id::Integer)

Get the presets for offset and gain values at different "sweet spots".

# Args:
    id: Camera id

# Returns:
    A dictionary containing:
        1. The offset at highest dynamic range
        2. The offset at unity gain
        3. The gain with lowest readout noise
        4. The offset with lowest readout noise

# Throws:
    ASIError
"""
function get_gain_offset(id::Integer)
    offset_highest_dr = Ref{Cint}(0)  # Offset at highest dynamic range
    offset_unity_gain = Ref{Cint}(0)  # Offset at unity gain
    gain_lowest_rn    = Ref{Cint}(0)  # Gain at lowest readout noise
    offset_lowest_rn  = Ref{Cint}(0)  # Offset at lowest readout noise

    err = ASIGetGainOffset(id, offset_highest_dr, offset_unity_gain, gain_lowest_rn, offset_lowest_rn)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end

    return Dict(["offset_highest_dr" => offset_highest_dr[],
                 "offset_unity_gain" => offset_unity_gain[],
                 "gain_lowest_rn" => gain_lowest_rn[],
                 "offset_lowest_rn" => offset_lowest_rn[]])
end

get_gain_offset(cam::ASICamera) = get_gain_offset(cam.info.CameraID)

"""
    get_sdk_version()

Returns the SDK version.
"""
function get_sdk_version()
    return unsafe_string(ASIGetSDKVersion())
end

"""
    get_supported_modes(id::Integer)

Returns the supported camera modes, only need to call when the
IsTriggerCam in the CameraInfo is true.

# Throws:
    ASIError
"""
function get_supported_modes(id::Integer)
    supported_modes = ASI_SUPPORTED_MODE()

    err = ASIGetCameraSupportMode(id, supported_modes)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end

    return [m for m in supported_modes.SupportedCameraModes if m > ASI_MODE_END]
end

get_supported_modes(cam::ASICamera) = get_supported_modes(cam.info.CameraID)

"""
    get_camera_mode(id::Integer)

Get the current camera mode, only needed to call when the IsTriggerCam
in the CameraInfo is true.
"""
function get_camera_mode(id::Integer)
    mode = Ref{ASI_CAMERA_MODE}(ASI_MODE_END)

    err = ASIGetCameraMode(id, mode)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end

    return mode[]
end

get_camera_mode(cam::ASICamera) = get_camera_mode(cam.info.CameraID)

"""
    set_camera_mode(id::Integer, mode::ASI_CAMERA_MODE)

Set the camera mode, only needed to call when the IsTriggerCam
in the CameraInfo is true.
"""
function set_camera_mode(id::Integer, mode::ASI_CAMERA_MODE)
    err = ASISetCameraMode(id, mode)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

set_camera_mode(cam::ASICamera) = set_camera_mode(cam.info.CameraID)

"""
    send_soft_trigger(id::Integer, start::ASI_BOOL)

From original docs: Send a softTrigger. For edge trigger, it only need to set
true which means send a rising trigger to start exposure. For level trigger,
it need to set true first means start exposure, and set false means stop
exposure. Only needed to call when the IsTriggerCam in the CameraInfo is true.
"""
function send_soft_trigger(id::Integer, start::ASI_BOOL)
    err = ASISendSoftTrigger(id, start)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

send_soft_trigger(cam::ASICamera) = send_soft_trigger(cam.info.CameraID)

"""
    get_serial_number(id::Integer)

Returns the camera serial number.
"""
function get_serial_number(id::Integer)
    sn = ASI_SN()

    err = ASIGetSerialNumber(id, sn)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end

    return unsafe_string(sn.id)
end

get_serial_number(cam::ASICamera) = get_serial_number(cam.info.CameraID)

"""
    set_trigger_output_config(id::Integer, pin::ASI_TRIG_OUTPUT_PIN, high::ASI_BOOL, delay, duration)

Configure the output pin (A or B) of Trigger port. If duration <= 0,
this output pin will be closed. Only need to call when the IsTriggerCam
in the CameraInfo is true.

# Args:
    id: Camera id
    pin: Select the pin for output
    high: If true, the selected pin will output a high level as a signal
					when it is effective.
    delay: The delay between the camera receiving a trigger signal and the
            output of the valid level. From 0 μs - 2,000,000,000 μs.
    duration: The duration of the valid level output. Same range as delay.

"""
function set_trigger_output_config(id::Integer, pin::ASI_TRIG_OUTPUT_PIN, high::ASI_BOOL, delay_μs, duration_μs)
    err = ASISetTriggerOutputIOConf(id, pin, high, delay_μs, duration_μs)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end
end

set_trigger_output_config(cam::ASICamera) = set_trigger_output_config(cam.info.CameraID)

"""
    get_trigger_output_config(id::Integer)

Get the output pin configuration, only needed to call when the IsTriggerCam
in the CameraInfo is true.
"""
function get_trigger_output_config(id::Integer)
    pin = ASI_TRIG_OUTPUT_NONE
    high = ASI_FALSE
    delay = Ref{Clong}(0)
    duration = Ref{Clong}(0)

    err = ASIGetTriggerOutputIOConf(id, pin, high, delay, duration)
    if err != ASI_SUCCESS
        throw(ASIError(err))
    end

    return Dict(["pin" => pin,
                 "high" => high,
                 "delay" => delay,
                 "duration" => duration])
end

get_trigger_output_config(cam::ASICamera) = get_trigger_output_config(cam.info.CameraID)

# """
#     capture_video(id)
#
# Starts capturing a video in a while loop, displaying it in a Makie scene.
#
# Args:
#     id: Camera object or int
#
# Throws:
#     An error when the camera returns an image type which is unsupported by this package.
# """
# function capture_video(id)
#     # Camera stuff setup
#     exposure_μs, _ = get_control_value(id, ASI_EXPOSURE)
#     width, height, binning, img_type = get_roi_format(id)
#
#     buffer = allocate_buffer(width, height, img_type)
#
#     # Makie scene setup
#     scene = Scene() #resolution = (size(buffer,2),size(buffer,1)))
#     im = image!(scene, buffer, show_axis = false, scale_plot = false, colorrange=(0,255))[end]
#     display(scene)
#
#     function video_loop()
#         while isopen(scene)
#             get_video_data!(cam, buffer, 0.2*exposure_μs+500)
#             im[1] = buffer[end:-1:1, :]  # flip x axis for some reason...
#             yield()
#         end
#     end
#
#     start_video(id)
#
#     try
#         video_loop()
#     catch err
#         rethrow(err)
#     finally
#         stop_video(id)
#     end
#
#     return nothing
# end
#
# capture_video(cam::ASICamera, callback::Function) = capture_video(cam.info.CameraID, callback)


text:
module MiniGNN

#  to import back in python https://stackoverflow.com/questions/33811864/how-do-i-within-python-retrieve-a-sparse-matrix-stored-by-julia-in-a-jld-file
#using GeometricFlux
using SparseArrays
using LinearAlgebra: issymmetric
using GraphSignals: adjacency_matrix

using Flux
using Functors: @functor, functor, fmap
using Flux: glorot_uniform, glorot_normal
using Flux: MaxPool, relu, Dropout

include("graph_layers.jl")
export make_rescaled_mat, adddims, process_data, create_batch

function make_rescaled_mat(adj_mat::AbstractMatrix{R}) where R<: Real
    deg = vec(sum(adj_mat,dims=1))
    denom = 1 ./ .√(1 .+ deg)
    denom .* permutedims(denom) .*adj_mat
end

function make_rescaled_mat(g::SparseMatrixCSC{R}) where R<:Real
    make_rescaled_mat(Matrix{Float32}(adjacency_matrix(g)))
end

function adddims(x::AbstractArray,k::Integer)
    @assert k<= ndims(x)+1
    add_last = reshape(x, (size(x)...,1))
    perm = Int[i for i = 1:ndims(add_last)]
    perm[k] = ndims(add_last)
    for i = k+1:ndims(add_last)
        perm[i] -= 1
    end
    permutedims(add_last,perm)
end

function process_data(data::AbstractArray)
    d = permutedims(data)
    return adddims(d,1)
end

function create_batch(bs::Int,train_loader)
    batch = []
    labels = []
    for (b,l) in train_loader
        batch = b[:,1:bs]
        labels = l[:,1:bs]
        break
    end
    batch, labels
end

end # module


text:
__precompile__()

module Celeste

# import this before StaticArrays has a chance to load, to
# workaround a bug
import ForwardDiff

include("aliasscopes.jl")
include("config.jl")

# submodules
include("Log.jl")

include("SensitiveFloats.jl")
include("BivariateNormals.jl")

include("Coordinates.jl")

include("Model.jl")
include("Infer.jl")
include("Transform.jl")
include("PSF.jl")
include("SDSSIO.jl")
include("SEP.jl")

include("MCMC.jl")
include("StochasticVI.jl")
include("DeterministicVI.jl")
include("ParallelRun.jl")

include("Synthetic.jl")
include("AccuracyBenchmark.jl")
include("GalsimBenchmark.jl")

include("ArgumentParse.jl")

include("settings.jl")

import .ParallelRun: BoundingBox, OptimizedSource
import .SDSSIO: RunCamcolField, CatalogEntry
export BoundingBox, OptimizedSource, RunCamcolField, CatalogEntry

end # module


text:
function set_startvals!(
   start_vals,
   va, vm, pg, qg, p, q, f
   )
   # Set starting values if given
   for (subnet_idx, subnet) in start_vals["sn"]
      subnet_idx_int = parse(Int64,subnet_idx)
      if subnet_idx_int in keys(vm)
         for (bus_i,bus) in subnet["bus"]
            bus_i_int = bus["bus_i"]
            # println("bus_i_int: $bus_i_int")
            # println("keys(vm[subnet_idx_int]): $(axes(vm[subnet_idx_int])[1])")
            if (bus_i_int in axes(vm[subnet_idx_int])[1]) && ("vm" in keys(bus))
               JuMP.set_start_value(vm[subnet_idx_int][bus_i_int], bus["vm"])
               JuMP.set_start_value(va[subnet_idx_int][bus_i_int], bus["va"])
               println("set bus $bus_i_int vm and va starting values")
            end
         end
         for (gen_i,gen) in subnet["gen"]
            gen_i_int = gen["index"]
            if (gen_i_int in axes(pg[subnet_idx_int])[1]) && ("pg" in keys(gen))
               JuMP.set_start_value(pg[subnet_idx_int][gen_i_int], gen["pg"])
               JuMP.set_start_value(qg[subnet_idx_int][gen_i_int], gen["qg"])
               println("set gen $gen_i_int pg and qg starting values")
            end
         end
         for (branch_i,branch) in subnet["branch"]
            f_idx = (branch["index"], branch["f_bus"], branch["t_bus"])
            t_idx = (branch["index"], branch["t_bus"], branch["f_bus"])
            if (f_idx in axes(p[subnet_idx_int])[1]) && ("pf" in keys(branch))
               JuMP.set_start_value(p[subnet_idx_int][f_idx], branch["pf"])
               JuMP.set_start_value(p[subnet_idx_int][t_idx], branch["pt"])
               JuMP.set_start_value(q[subnet_idx_int][f_idx], branch["qf"])
               JuMP.set_start_value(q[subnet_idx_int][t_idx], branch["qt"])
               println("set branch $(branch["index"]) p and q starting values")
            end
         end
         if "frequency" in keys(subnet)
            if ref[:sn][subnet_idx_int][:variable_f]
               JuMP.set_start_value(f[subnet_idx], subnet["frequency"])
               println("set subnet $subnet_idx_int frequency starting value")
            end
         end
      end
   end

end


text:
#Ensemble Kalman Sampler: specific structures and function definitions

"""
    Sampler{FT<:AbstractFloat,IT<:Int} <: Process

An ensemble Kalman Sampler process.
"""
struct Sampler{FT <: AbstractFloat} <: Process
    ""
    prior_mean::Vector{FT}
    ""
    prior_cov::Array{FT, 2}
end

function update_ensemble!(ekp::EnsembleKalmanProcess{FT, IT, Sampler{FT}}, g_in::Array{FT, 2}) where {FT, IT}

    #catch works when g_in non-square
    if !(size(g_in)[2] == ekp.N_ens)
        throw(DimensionMismatch("ensemble size in EnsembleKalmanProcess and g_in do not match, try transposing or check ensemble size"))
    end

    # u: N_ens × N_par
    # g: N_ens × N_obs
    u_old = get_u_final(ekp)
    u_old = permutedims(u_old, (2, 1))
    u = u_old
    g = permutedims(g_in, (2, 1))

    # u_mean: N_par × 1
    u_mean = mean(u', dims = 2)
    # g_mean: N_obs × 1
    g_mean = mean(g', dims = 2)
    # g_cov: N_obs × N_obs
    g_cov = cov(g, corrected = false)
    # u_cov: N_par × N_par
    u_cov = cov(u, corrected = false)

    # Building tmp matrices for EKS update:
    E = g' .- g_mean
    R = g' .- ekp.obs_mean
    # D: N_ens × N_ens
    D = (1 / ekp.N_ens) * (E' * (ekp.obs_noise_cov \ R))

    Δt = 1 / (norm(D) + 1e-8)

    noise = MvNormal(u_cov)

    implicit =
        (1 * Matrix(I, size(u)[2], size(u)[2]) + Δt * (ekp.process.prior_cov' \ u_cov')') \
        (u' .- Δt * (u' .- u_mean) * D .+ Δt * u_cov * (ekp.process.prior_cov \ ekp.process.prior_mean))

    u = implicit' + sqrt(2 * Δt) * rand(noise, ekp.N_ens)'

    # store new parameters (and model outputs)
    push!(ekp.u, DataContainer(u, data_are_columns = false))
    push!(ekp.g, DataContainer(g, data_are_columns = false))
    # u_old is N_ens × N_par, g is N_ens × N_obs,
    # but stored in data container with N_ens as the 2nd dim

    compute_error!(ekp)

end


text:
export Iris

"""
Fisher's classic iris dataset.

Measurements from 3 different species of iris: setosa, versicolor and
virginica.  There are 50 examples of each species.

There are 4 measurements for each example: sepal length, sepal width, petal
length and petal width.  The measurements are in centimeters.

The module retrieves the data from the [UCI Machine Learning Repository](https://archive.ics.uci.edu/ml/datasets/iris).

NOTE: no pre-defined train-test split for this dataset, `features` and `labels` return the whole dataset.

# Interface

- [`Iris.features`](@ref)
- [`Iris.labels`](@ref)

# Utilities

- [`Iris.download`](@ref)
"""
module Iris

using DataDeps
using ..MLDatasets: bytes_to_type, datafile, download_dep, download_docstring
using DelimitedFiles

export features, labels, download

const DEPNAME = "Iris"
const LINK = "https://archive.ics.uci.edu/ml/machine-learning-databases/iris/"
const DOCS = "https://archive.ics.uci.edu/ml/datasets/Iris"
const DATA = "iris.data"

"""
    download([dir]; [i_accept_the_terms_of_use])

Trigger the (interactive) download of the full dataset into
"`dir`". If no `dir` is provided the dataset will be
downloaded into "~/.julia/datadeps/$DEPNAME".

This function will display an interactive dialog unless
either the keyword parameter `i_accept_the_terms_of_use` or
the environment variable `DATADEPS_ALWAY_ACCEPT` is set to
`true`. Note that using the data responsibly and respecting
copyright/terms-of-use remains your responsibility.
"""
download(args...; kw...) = download_dep(DEPNAME, args...; kw...)

function __init__()
    register(DataDep(
        DEPNAME,
        """
        Dataset: The Iris dataset
        Website: $DOCS
        """,
        LINK .* [DATA],
        "6f608b71a7317216319b4d27b4d9bc84e6abd734eda7872b71a458569e2656c0"  # if checksum omitted, will be generated by DataDeps
    ))
end

"""
    labels(; dir = nothing)

Return a string vector of length 150 containing observations' labels.
"""
function labels(; dir = nothing)
    path = datafile(DEPNAME, DATA, dir)
    iris = readdlm(path, ',')
    Vector{String}(iris[:, end])
end

"""
    features(; dir = nothing)

Return a 4x150 matrix containing the 4-dimensional features of each observation.
"""
function features(; dir = nothing)
    path = datafile(DEPNAME, DATA, dir)
    iris = readdlm(path, ',')
    Matrix{Float64}(iris[:, 1:4])' |> collect
end

end # module