text:
#
# Denoise an SPD Example with Douglas Rachford
#
# Denoise an SPD Example with Parallel Douglas-Rachford to minimize the
#
# L2-TV functional with anisotropic TV
#
# where the example is the same data and cost as for SPD_Image_CP_vs_CPPA
#
# This example and its data is used for comparison in Section 6.2 of
#
# > R. Bergmann, R. Herzog, M. Silva Louzeiro, D. Tenbrinck, J. Vidal Núñez:
# > _Fenchel Duality Theory and a Primal-Dual Algorithm on Riemannian Manifolds_,
# > arXiv: [1908.02022](https://arxiv.org/abs/1908.02022)
#
using Manopt, Manifolds
using Images, CSV, DataFrames, LinearAlgebra, JLD2, Dates
#
# Settings
experiment_name = "SPD_Image_DR"
export_orig = true
export_result = true
export_table = true
use_debug = true
asy_render_detail = 2
results_folder = joinpath(@__DIR__, "Image_TV")
comparison_data = joinpath(
    @__DIR__, "..", "CyclicProximalPoint", "Image_TV", "SPD_Image_CPPA-cost.jld2"
)
!isdir(results_folder) && mkdir(results_folder)
#
# Parameters
η = 0.58
λ = 0.93
α = 6.0
#
# Manifold & Data
f = artificial_SPD_image2(32)
if export_orig
    fn = joinpath(results_folder, experiment_name * "-orig.asy")
    asymptote_export_SPD(fn; data=f, scale_axes=(7.5, 7.5, 7.5))
    render_asymptote(fn; render=asy_render_detail)
end
sC = StopAfterIteration(400)
try
    cppa_data = load(comparison_data)
    cost_threshold = cppa_data["cost_function_value"]
    cppa_iter = cppa_data["iterations"]
    global sC = StopWhenAny(StopAfterIteration(400), StopWhenCostLess(cost_threshold))
    @info "Comparison to CPPA (`SPDImage_CPPA.jl`) and its cost after $(cppa_iter) iterations of $cost_threshold."
catch e
    msg = sprint(showerror, e)
    @info "Comparison to CPPA (`CyclicProximalPoint/SPDImage_CPPA.jl`) not possible, data file $(comparison_data) either missing or corrupted.\n Error: $(msg).\n\n Starting with a default stopping criterion."
end
#
# Build Problem for L2-TV
pixelM = SymmetricPositiveDefinite(3);
M = PowerManifold(pixelM, NestedPowerRepresentation(), size(f)...)
d = length(size(f))
rep(d) = (d > 1) ? [ones(Int, d)..., d] : d
fidelity(x) = 1 / 2 * distance(M, x, f)^2
Λ(x) = forward_logs(M, x) # on T_xN
prior(x) = norm(norm.(Ref(pixelM), repeat(x, rep(d)...), Λ(x)), 1)
#
# Setup & Optimize
print("--- Douglas–Rachford with η: $(η) and λ: $(λ) ---\n")
cost(x) = 1 / α * fidelity(x[1]) + prior(x[1])
N = PowerManifold(M, NestedPowerRepresentation(), 5)
prox1 = (η, x) -> [prox_distance(M, η, f, x[1]), prox_parallel_TV(M, α * η, x[2:5])...]
prox2 = (η, x) -> fill(mean(M, x, GradientDescentEstimation(); stop_iter=4), 5)
x0 = fill(f, 5)
@time o = DouglasRachford(
    N,
    cost,
    [prox1, prox2],
    x0;
    λ=i -> η,
    α=i -> λ, # map from Paper notation of BPS16 to toolbox notation
    debug=use_debug ? [:Iteration, " | ", :Cost, "\n", 10, :Stop] : missing,
    record=export_table ? [:Iteration, :Cost] : missing,
    stopping_criterion=sC,
    return_options=true,
)
y = get_solver_result(o)[1]
export_table && (r = get_record(o))
#
# Result
numIter = length(r)
if export_result
    fn = joinpath(
        results_folder,
        experiment_name * "img-result-$(numIter)-α$(replace(string(α), "." => "-")).asy",
    )
    asymptote_export_SPD(fn; data=y, scale_axes=(7.5, 7.5, 7.5))
    render_asymptote(fn; render=asy_render_detail)
end
if export_table
    # scale cost back for saving such that its comparable with the other two results
    A = cat([ri[1] for ri in r], [ri[2] / α for ri in r]; dims=2)
    CSV.write(
        joinpath(results_folder, experiment_name * "-Cost.csv"),
        DataFrame(A);
        writeheader=false,
    )
end


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.9.1
#   kernelspec:
#     display_name: Julia 1.6.0
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# # Introduction
#
# In this notebook we will optimize disk spacings using the transformer algorithm. When optimizing for a certain frequency range $\nu_{min} - \nu_{max}$
# we minimize the cost function $C=-min(\beta ^2 (\nu) | \nu_{min} < \nu < \nu_{max})$. Essentially we take the worst frequency point of the boostfactor curve in this range and try to increase it (by decreasing the negative). This will naturally lead to a flat top hat shape.
#
# To reduce computation we only calculate the boostfactor in the range $\nu_{min} - \nu_{max}$ and keep the number of frequency points to a sensible minimum. With increasing disk number more points are needed to resolve the "wiggles" on top of the boostfactor.
#
# One also has to choose what boostfactor mode to optimize for. Fundamental mode, other modes, coupled power?
#
# Transformer calculates the mixing between modes which is computationally the most expensive part. In order to make optimization feasable we calculate the mixing for a selected number of points beforehand and then quickly interpolate the exact spacing during optimization.
#
# For gradient based methods we have to approximate the non-differentiable minimum() with a smooth function like the p-norm.
# Still, very close to the minimum the gradient might be unstable and gradient based methods might struggle.
#

# %%
push!(LOAD_PATH, pwd() * "/../src")
using BoostFractor
using Optim
using ForwardDiff #If we use gradient-based optimization
include("transformer_optim_utilities.jl") #Bunch of helper functions
using PyPlot
using Base.Threads
using Dates
using DelimitedFiles

# %% [markdown]
# # Setup

# %%
#Physical Parameters
n_disk = 20
epsilon = 24
diskR = 0.15
Mmax = 1
Lmax = 0
#We start with the most generic spacing possible: Equidistant spacing. Make sure the peak is roughly located in the desired
#frequency range.
init_spacing = 0.007118

# %%
#Frequency range to optimize. For 20 disks 8 points seems sufficient, 80 disks needs ~20-25
freq_center = 22.288e9
#freq_center = 20.e9
freq_width = 50e6
n_freq_optim = 8
#freq_min = freq_center-freq_width/2
#freq_max = freq_center+freq_width/2
frequencies_optim(center, width) = range(center - width / 2,stop=center + width / 2,
                                         length=n_freq_optim);
c_0 = 2.99792458e8
λ1 = c_0 / (freq_center * sqrt(epsilon))
λ2 = c_0 / freq_center
println("λ1: $λ1 λ2: $λ2, $(λ1/sqrt(epsilon))")
#init_spacing = 0.00808824
#println(init_spacing)

# %%

#For Plotting
df = 0.004*1e9
frequencies_plot = 22.1e9:df:22.6e9;
#frequencies_plot = 19.0e9:df:22.5e9

#Usual transformer setup
# 1D Coordinate System
X = [1e-9]
Y = [1e-9]
coords = SeedCoordinateSystem(X = X, Y = Y)
n_region = 2*n_disk+2
prop = propagator1D

# SetupBoundaries (note that this expects the mirror to be defined explicitly as a region)
eps = Array{Complex{Float64}}([i==1 ? 1e20 : i%2==0 ? 1 : epsilon for i=1:n_region])

distance = Array{Float64}([i==1 ? 0 : i%2==0 ? init_spacing : 1e-3 for i=1:n_region])
distance[end] = 0

#Uncomment to include thickness varation
init_thickness_variation = zeros(n_region,length(X),length(Y))
#xi = 0.035
#sigma = 2e-6;
#init_thickness_variation = surface_roughness(X,Y,n_region,mag=sigma,trunc=sigma*10, xi=xi, diskR=diskR)

sbdry_init = SeedSetupBoundaries(coords, diskno=n_disk, distance=distance, epsilon=eps,
                                 relative_surfaces=init_thickness_variation)

# Initialize modes
modes = SeedModes(coords, ThreeDim=false, Mmax=Mmax, Lmax=Lmax, diskR=diskR)

#  Mode-Vector defining beam shape to be reflected on the system
m_reflect = zeros(Mmax*(2*Lmax+1))
m_reflect[Lmax+1] = 1.0;

# # Plot Initial Boostfactor

#Calculate FULL frequency prop matrix of initial setup for plotting
prop_matrix_grid_plot = calc_propagation_matrices_grid(sbdry_init, coords, modes, 0,
                                                       frequencies_plot; diskR=diskR,
                                                       prop=prop);
prop_matrix_plot = [prop_matrix_grid_plot[r,f,1,1,1,:,:] for r=1:n_region, f=1:length(frequencies_plot)]
Eout_init = calc_boostfactor_modes(sbdry_init,coords,modes,frequencies_plot, prop_matrix_plot,
                                   diskR=diskR, prop=prop)

#Only plot fundamental mode
plot(frequencies_plot.*1e-9,abs2.(Eout_init)[1,Lmax+1,:])
xlabel("Frequency [GHz]")
ylabel("Power Boostfactor")

# %% [markdown]
# # Propagation Matrix Interpolation

# %%
#Calculate interpolation object for frequencies we want to optimize
#To avoid calculating the costly propagation matrix for each possible distance, we calculate it only at
#some points and interpolate the others. The spacing grid is relative distance changes from the provided sbdry struct.
#The grid steps shouldn't be larger than 50 um. Choose a sensible range from to safe computation (e.g. a small range if you
# are already close to the optimum, larger if far away)
spacing_grid = (-2500:50:2500)*1e-6;

prop_matrix_grid_sub(center) = calc_propagation_matrices_grid(sbdry_init, coords, modes,
                                    spacing_grid, frequencies_optim(center, freq_width),
                                    diskR=diskR, prop=prop);
#This is the interpolation object that is being called when interpolating
itp_sub(center) = construct_prop_matrix_interpolation(prop_matrix_grid_sub(center),
                                                      spacing_grid)

# %% [markdown]
# # Optimization

# %%
spacings = zeros(n_disk)
#We define a single argument cost function for the optimizer
#m_reflect defines the combination of modes we want to optimize. For now we are only interested in the fundamental mode
function cost_fun(fixed_disk, center)
    disk = copy(fixed_disk)
    itp = itp_sub(center)
    freq_o = frequencies_optim(center, freq_width)
    return x -> begin
        # if one disk is fixed, insert it into dist_shift as the optimizer then runs on one
        # dimension less but we still need the "full" list of disks
        if disk > 0
            # it's important to *copy* and not modify x here otherwise the optimizer gets confused
            x = vcat(x[1:disk - 1], spacings[disk], x[disk:length(x)])
        end

        calc_boostfactor_cost(x, itp, freq_o, sbdry_init, coords, modes, m_reflect,
                              diskR=diskR, prop=prop)
    end
end;

# %%
# With NelderMead. Slow but robust
#algorithm = NelderMead()
#options = Optim.Options(iterations=1000000)

#Quasi-Newton (gradient-based). Fast but likes to get stuck in local minima.
algorithm = BFGS()
options = Optim.Options(f_tol = 1e-6)

# %%
GC.gc()

#spacings = [-0.000247484690809835,
#             5.494161169410954e-5,
#             1.4082408393016912e-5,
#             3.474533260803849e-6,
#             1.4127821024835357e-6,
#            -7.062735634745444e-6,
#            -5.658759242525849e-5,
#            -0.0001736356257557033,
#             0.0003140458044571532,
#             0.000946244385671497,
#             6.889677168591535e-6,
#            -0.00015884778788364403,
#            -0.0001663066683797041,
#            -0.00035282832642484133,
#             0.0008533180045262031,
#            -0.002345992812444579,
#             0.0009458750290604845,
#            -0.0004042461326743376,
#            -0.00043332980022113657,
#             0.0004518582202935378]

function optimize_spacings(cost_function; starting_point=zeros(n_disk))
    spacings = Vector{Float64}()
    best_cost = Atomic{Float64}(1000.)
    lk = SpinLock()
    # Run initial optimization a few times and pick the best one
    @threads for i in 1:100
        # Add some random variation to start spacing.
        # Convergence very much depends on a good start point.
        x_0 = starting_point .+ 2 .* (rand(length(starting_point)).-0.5) .* 100e-6
        od = OnceDifferentiable(cost_function, x_0, autodiff=:forward)
        res = optimize(od, x_0, algorithm, options)
        cost = cost_function(Optim.minimizer(res))
        lock(lk)
        atomic_min!(best_cost, cost)
        if atomic_cas!(best_cost, cost, cost) === cost
            println("Best cost: $cost")
            spacings = Optim.minimizer(res)
        end
        unlock(lk)
    end
    spacings
end

# spacings = optimize_spacings(cost_fun(0, freq_center))
spacings = readdlm("optim_results/22.288GHz-2021-04-14T13:13.txt", ',')[:]

shift = 40 # Shift in MHz

# %%
results = Vector{Vector{Float64}}(undef, n_disk + 1)
for i in 0:n_disk
    println("Optimizing for fixed disk $i")
    n_moving_disks = n_disk - (i > 0 ? 1 : 0)
    x_0 = vcat(spacings[1:i-1], spacings[i+1:n_disk])
    results[i + 1] = optimize_spacings(cost_fun(i, freq_center + shift * 1e6), starting_point=x_0)
end

# %%
function calc_eout(distances; fixed_disk=0)
    dist = copy(distances)
    if fixed_disk > 0
        insert!(dist, fixed_disk, spacings[fixed_disk])
    end
    sbdry_optim = copy_setup_boundaries(sbdry_init, coords)
    sbdry_optim.distance[2:2:end-2] .+= dist

    #Calculate prop matrix grid at a dist shift of zero of optimized setup
    prop_matrix_grid_plot = calc_propagation_matrices_grid(sbdry_optim, coords,modes, 0,
                                                           frequencies_plot, diskR=diskR,
                                                           prop=prop)
    prop_matrix_plot = [prop_matrix_grid_plot[r,f,1,1,1,:,:] for r = 1:n_region,
                        f = 1:length(frequencies_plot)]
    calc_boostfactor_modes(sbdry_optim, coords, modes, frequencies_plot,
                                           prop_matrix_plot, diskR=diskR, prop=prop)
end
#Plot results
#Create optimized setup
GC.gc()
Eout_0 = calc_eout(spacings)

#results = readdlm("optim_results/22.288GHz+30MHz-2021-04-14T18:01.txt", ',')
#results = [filter(x -> x != "", results[i,:]) for i in 1:size(results, 1)]
Eout_optim = Array{Array{ComplexF64, 3}}(undef, n_disk + 1)
@time @threads for (i, res) in collect(enumerate(results))
    Eout_optim[i] = calc_eout(res, fixed_disk=i - 1)
end

# %%
function write_to_file(data, name)
    date = Dates.format(now(), "yyyy-mm-ddTHH:MM")
    open("optim_results/$name-$date.txt", "w") do io
        writedlm(io, data, ',')
    end
end

#write_to_file(spacings, "$(freq_center * 1e-9)GHz")
write_to_file(results, "$(freq_center * 1e-9)GHz+$(shift)MHz")

# %%
v = Dates.format(now(), "yyyy-mm-ddTHH:MM")

# %%
clf()
legend_text = Array{String}(undef, 0)
push!(legend_text, "freq_0")
plot(frequencies_plot .* 1e-9, abs2.(Eout_0[1, Lmax + 1, :]))
for i in 1:6
    out = Eout_optim[i]
    plot(frequencies_plot.*1e-9,abs2.(out[1, Lmax + 1, :]))
    push!(legend_text, "Fixed disk: $(i-1)")
end
xlabel("Frequency [GHz]")
ylabel("Power Boostfactor")
legend(legend_text)
savefig("+$(shift)MHz_1_$v", dpi=480)

# %%
clf()
legend_text = Array{String}(undef, 0)
push!(legend_text, "freq_0")
plot(frequencies_plot .* 1e-9, abs2.(Eout_0[1, Lmax + 1, :]))
for i in vcat(1, 7:11)
    out = Eout_optim[i]
    plot(frequencies_plot.*1e-9,abs2.(out)[1,Lmax+1,:])
    push!(legend_text, "Fixed disk: $(i-1)")
end
xlabel("Frequency [GHz]")
ylabel("Power Boostfactor")
legend(legend_text)
savefig("+$(shift)MHz_2_$v", dpi=480)

# %%
clf()
legend_text = Array{String}(undef, 0)
push!(legend_text, "freq_0")
plot(frequencies_plot .* 1e-9, abs2.(Eout_0[1, Lmax + 1, :]))
for i in vcat(1, 12:16)
    out = Eout_optim[i]
    plot(frequencies_plot.*1e-9,abs2.(out)[1,Lmax+1,:])
    push!(legend_text, "Fixed disk: $(i-1)")
end
xlabel("Frequency [GHz]")
ylabel("Power Boostfactor")
legend(legend_text)
savefig("+$(shift)MHz_3_$v", dpi=480)

# %%
clf()
legend_text = Array{String}(undef, 0)
push!(legend_text, "freq_0")
plot(frequencies_plot .* 1e-9, abs2.(Eout_0[1, Lmax + 1, :]))
for i in vcat(1, 17:21)
    out = Eout_optim[i]
    plot(frequencies_plot.*1e-9,abs2.(out)[1,Lmax+1,:])
    push!(legend_text, "Fixed disk: $(i-1)")
end
xlabel("Frequency [GHz]")
ylabel("Power Boostfactor")
legend(legend_text)
savefig("+$(shift)MHz_4_$v", dpi=480)

# %%

ext: jl

max_stars_repo_path: examples/Transformer, Fixed Disk.jl

max_issues_repo_path: examples/Transformer, Fixed Disk.jl

max_forks_repo_path: examples/Transformer, Fixed Disk.jl


text:
module ConstantArrays
export ConstantArray

struct ConstantArray{T, N} <: AbstractArray{T, N}
    _value::T
    _size::NTuple{N, Int64}
end

(::Type{ConstantArray})(value::Any, dims::Vararg{<: Integer}) =
    ConstantArray{typeof(value), length(dims)}(value,
                                               convert(NTuple{length(dims), Int64}, dims))

Base.size(q::ConstantArray) = q._size
Base.getindex(q::ConstantArray, ::Integer) = q._value
Base.getindex(q::ConstantArray, ::Vararg{<: Integer}) = q._value

Base.show(io::IO, q::ConstantArray) = show(io, "ConstantArray($(q.value), $(q.dims))")
end

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
"""
    MHAttention(nheads::Integer, qkv_layer, attn_drop, projection)

Multi-head self-attention layer.

# Arguments:

  - `nheads`: Number of heads
  - `qkv_layer`: layer to be used for getting the query, key and value
  - `attn_drop`: dropout rate after the self-attention layer
  - `projection`: projection layer to be used after self-attention
"""
struct MHAttention{P, Q, R}
    nheads::Int
    qkv_layer::P
    attn_drop::Q
    projection::R
end

"""
    MHAttention(planes::Integer, nheads::Integer = 8; qkv_bias::Bool = false, attn_drop = 0., proj_drop = 0.)

Multi-head self-attention layer.

# Arguments:

  - `planes`: number of input channels
  - `nheads`: number of heads
  - `qkv_bias`: whether to use bias in the layer to get the query, key and value
  - `attn_drop`: dropout rate after the self-attention layer
  - `proj_drop`: dropout rate after the projection layer
"""
function MHAttention(planes::Integer, nheads::Integer = 8; qkv_bias::Bool = false,
                     attn_drop = 0.0, proj_drop = 0.0)
    @assert planes % nheads==0 "planes should be divisible by nheads"
    qkv_layer = Dense(planes, planes * 3; bias = qkv_bias)
    attn_drop = Dropout(attn_drop)
    proj = Chain(Dense(planes, planes), Dropout(proj_drop))
    return MHAttention(nheads, qkv_layer, attn_drop, proj)
end

@functor MHAttention

function (m::MHAttention)(x::AbstractArray{T, 3}) where {T}
    nfeatures, seq_len, batch_size = size(x)
    x_reshaped = reshape(x, nfeatures, seq_len * batch_size)
    qkv = m.qkv_layer(x_reshaped)
    qkv_reshaped = reshape(qkv, nfeatures ÷ m.nheads, m.nheads, seq_len, 3 * batch_size)
    query, key, value = chunk(qkv_reshaped, 3; dims = 4)
    scale = convert(T, sqrt(size(query, 1) / m.nheads))
    key_reshaped = reshape(permutedims(key, (2, 1, 3, 4)), m.nheads, nfeatures ÷ m.nheads,
                           seq_len * batch_size)
    query_reshaped = reshape(permutedims(query, (1, 2, 3, 4)), nfeatures ÷ m.nheads, m.nheads, seq_len * batch_size)
    attention = m.attn_drop(softmax(batched_mul(query_reshaped, key_reshaped) .* scale))
    value_reshaped = reshape(permutedims(value, (1, 2, 3, 4)), nfeatures ÷ m.nheads, m.nheads, seq_len * batch_size)
    pre_projection = reshape(batched_mul(attention, value_reshaped),
                             (nfeatures, seq_len, batch_size))
    y = m.projection(reshape(pre_projection, size(pre_projection, 1), :))
    return reshape(y, :, seq_len, batch_size)
end


text:
using DocStringExtensions
using KernelAbstractions

using ..Mesh.Geometry
using ..VariableTemplates

import ..Mesh.Grids:
    _ξ1x1, _ξ2x1, _ξ3x1, _ξ1x2, _ξ2x2, _ξ3x2, _ξ1x3, _ξ2x3, _ξ3x3
import ..MPIStateArrays: array_device

"""
    VectorGradients{
        FT <: AbstractFloat,
        FTA2D <: AbstractArray{FT, 2},
        FTA3D <: AbstractArray{FT, 3},
    }

This data structure stores the spatial gradients of a velocity field.

# Fields

$(DocStringExtensions.FIELDS)

# Usage

    VectorGradients(data)

# Arguments for the inner constructor
 - `data`: 3-dimensional device array containing the spatial gradients
   (the second dimension must be 9)
"""
struct VectorGradients{
    FT <: AbstractFloat,
    FTA2D <: AbstractArray{FT, 2},
    FTA3D <: AbstractArray{FT, 3},
}
    "Device array storing the spatial gradient data"
    data::FTA3D
    "View of ∂u₁/∂x₁"
    ∂₁u₁::FTA2D
    "View of ∂u₁/∂x₂"
    ∂₂u₁::FTA2D
    "View of ∂u₁/∂x₃"
    ∂₃u₁::FTA2D
    "View of ∂u₂/∂x₁"
    ∂₁u₂::FTA2D
    "View of ∂u₂/∂x₂"
    ∂₂u₂::FTA2D
    "View of ∂u₂/∂x₃"
    ∂₃u₂::FTA2D
    "View of ∂u₃/∂x₁"
    ∂₁u₃::FTA2D
    "View of ∂u₃/∂x₂"
    ∂₂u₃::FTA2D
    "View of ∂u₃/∂x₃"
    ∂₃u₃::FTA2D

    function VectorGradients(
        data::AbstractArray{FT, 3},
    ) where {FT <: AbstractFloat}
        ∂₁u₁, ∂₂u₁, ∂₃u₁ =
            view(data, :, 1, :), view(data, :, 2, :), view(data, :, 3, :)
        ∂₁u₂, ∂₂u₂, ∂₃u₂ =
            view(data, :, 4, :), view(data, :, 5, :), view(data, :, 6, :)
        ∂₁u₃, ∂₂u₃, ∂₃u₃ =
            view(data, :, 7, :), view(data, :, 8, :), view(data, :, 9, :)

        return new{FT, typeof(∂₁u₁), typeof(data)}(
            data,
            ∂₁u₁,
            ∂₂u₁,
            ∂₃u₁,
            ∂₁u₂,
            ∂₂u₂,
            ∂₃u₂,
            ∂₁u₃,
            ∂₂u₃,
            ∂₃u₃,
        )
    end
end

"""
    VectorGradients(dg::DGModel, Q::MPIStateArray)

This constructor computes the spatial gradients of the velocity field.

# Arguments
 - `dg`: DGModel
 - `Q`: MPIStateArray containing the conservative state variables
"""
function VectorGradients(dg::DGModel, Q::MPIStateArray)
    bl = dg.balance_law
    FT = eltype(dg.grid)
    N = polynomialorder(dg.grid)
    Nq = N + 1
    npoints = Nq^3
    nrealelem = length(dg.grid.topology.realelems)

    g = similar(Q.realdata, npoints, nrealelem, 3, 3)
    data = similar(Q.realdata, npoints, 9, nrealelem)

    ind = [
        varsindex(vars_state_conservative(bl, FT), :ρ)
        varsindex(vars_state_conservative(bl, FT), :ρu)
    ]
    _ρ, _ρu, _ρv, _ρw = ind[1], ind[2], ind[3], ind[4]

    device = array_device(Q)
    workgroup = (Nq, Nq)
    ndrange = (nrealelem * Nq, Nq)

    kernel = vector_gradients_kernel!(device, workgroup)
    event = kernel(
        Q.realdata,
        dg.grid.D,
        dg.grid.vgeo,
        g,
        data,
        _ρ,
        _ρu,
        _ρv,
        _ρw,
        Val(Nq),
        ndrange = ndrange,
    )
    wait(event)

    return VectorGradients(data)
end

@kernel function vector_gradients_kernel!(
    sv::AbstractArray{FT},
    D::AbstractArray{FT, 2},
    vgeo::AbstractArray{FT},
    g::AbstractArray{FT, 4},
    vgrad_data::AbstractArray{FT, 3},
    _ρ::Int,
    _ρu::Int,
    _ρv::Int,
    _ρw::Int,
    ::Val{qm1},
) where {qm1, FT <: AbstractFloat}

    e = @index(Group, Linear)
    i, j = @index(Local, NTuple)

    s_D = @localmem FT (qm1, qm1)
    s_U = @localmem FT (qm1, qm1)
    s_V = @localmem FT (qm1, qm1)
    s_W = @localmem FT (qm1, qm1)

    s_D[i, j] = D[i, j]
    @synchronize
    # computing derivatives with respect to ξ1
    for t in 1:qm1, s in 1:qm1
        ijk = j + ((s - 1) + (t - 1) * qm1) * qm1
        s_U[i, j] = s_D[i, j] * (sv[ijk, _ρu, e] / sv[ijk, _ρ, e])
        s_V[i, j] = s_D[i, j] * (sv[ijk, _ρv, e] / sv[ijk, _ρ, e])
        s_W[i, j] = s_D[i, j] * (sv[ijk, _ρw, e] / sv[ijk, _ρ, e])
        @synchronize
        if j == 1
            for r in 2:qm1
                s_U[i, 1] += s_U[i, r]
                s_V[i, 1] += s_V[i, r]
                s_W[i, 1] += s_W[i, r]
            end
        end
        @synchronize
        if j == 1
            g[i + ((s - 1) + (t - 1) * qm1) * qm1, e, 1, 1] = s_U[i, 1] # ∂u₁∂ξ₁
            g[i + ((s - 1) + (t - 1) * qm1) * qm1, e, 2, 1] = s_V[i, 1] # ∂u₂∂ξ₁
            g[i + ((s - 1) + (t - 1) * qm1) * qm1, e, 3, 1] = s_W[i, 1] # ∂u₃∂ξ₁
        end
    end
    @synchronize
    # computing derivatives with respect to ξ2
    for t in 1:qm1, r in 1:qm1
        ijk = r + ((j - 1) + (t - 1) * qm1) * qm1
        s_U[i, j] = s_D[i, j] * (sv[ijk, _ρu, e] / sv[ijk, _ρ, e])
        s_V[i, j] = s_D[i, j] * (sv[ijk, _ρv, e] / sv[ijk, _ρ, e])
        s_W[i, j] = s_D[i, j] * (sv[ijk, _ρw, e] / sv[ijk, _ρ, e])
        @synchronize
        if j == 1
            for s in 2:qm1
                s_U[i, 1] += s_U[i, s]
                s_V[i, 1] += s_V[i, s]
                s_W[i, 1] += s_W[i, s]
            end
        end
        @synchronize
        if j == 1
            g[r + ((i - 1) + (t - 1) * qm1) * qm1, e, 1, 2] = s_U[i, 1] # ∂u₁∂ξ₂
            g[r + ((i - 1) + (t - 1) * qm1) * qm1, e, 2, 2] = s_V[i, 1] # ∂u₂∂ξ₂
            g[r + ((i - 1) + (t - 1) * qm1) * qm1, e, 3, 2] = s_W[i, 1] # ∂u₃∂ξ₂
        end
    end
    @synchronize
    # computing derivatives with respect to ξ3
    for s in 1:qm1, r in 1:qm1
        ijk = r + ((s - 1) + (j - 1) * qm1) * qm1
        s_U[i, j] = s_D[i, j] * (sv[ijk, _ρu, e] / sv[ijk, _ρ, e])
        s_V[i, j] = s_D[i, j] * (sv[ijk, _ρv, e] / sv[ijk, _ρ, e])
        s_W[i, j] = s_D[i, j] * (sv[ijk, _ρw, e] / sv[ijk, _ρ, e])
        @synchronize
        if j == 1
            for t in 2:qm1
                s_U[i, 1] += s_U[i, t]
                s_V[i, 1] += s_V[i, t]
                s_W[i, 1] += s_W[i, t]
            end
        end
        @synchronize
        if j == 1
            g[r + ((s - 1) + (i - 1) * qm1) * qm1, e, 1, 3] = s_U[i, 1] # ∂u₁∂ξ₃
            g[r + ((s - 1) + (i - 1) * qm1) * qm1, e, 2, 3] = s_V[i, 1] # ∂u₂∂ξ₃
            g[r + ((s - 1) + (i - 1) * qm1) * qm1, e, 3, 3] = s_W[i, 1] # ∂u₃∂ξ₃
        end
    end
    @synchronize

    ∂₁u₁, ∂₂u₁, ∂₃u₁ = 1, 2, 3
    ∂₁u₂, ∂₂u₂, ∂₃u₂ = 4, 5, 6
    ∂₁u₃, ∂₂u₃, ∂₃u₃ = 7, 8, 9

    for k in 1:qm1
        ijk = i + ((j - 1) + (k - 1) * qm1) * qm1

        ξ1x1 = vgeo[ijk, _ξ1x1, e]
        ξ1x2 = vgeo[ijk, _ξ1x2, e]
        ξ1x3 = vgeo[ijk, _ξ1x3, e]
        ξ2x1 = vgeo[ijk, _ξ2x1, e]
        ξ2x2 = vgeo[ijk, _ξ2x2, e]
        ξ2x3 = vgeo[ijk, _ξ2x3, e]
        ξ3x1 = vgeo[ijk, _ξ3x1, e]
        ξ3x2 = vgeo[ijk, _ξ3x2, e]
        ξ3x3 = vgeo[ijk, _ξ3x3, e]

        vgrad_data[ijk, ∂₁u₁, e] =
            g[ijk, e, 1, 1] * ξ1x1 +
            g[ijk, e, 1, 2] * ξ2x1 +
            g[ijk, e, 1, 3] * ξ3x1
        vgrad_data[ijk, ∂₂u₁, e] =
            g[ijk, e, 1, 1] * ξ1x2 +
            g[ijk, e, 1, 2] * ξ2x2 +
            g[ijk, e, 1, 3] * ξ3x2
        vgrad_data[ijk, ∂₃u₁, e] =
            g[ijk, e, 1, 1] * ξ1x3 +
            g[ijk, e, 1, 2] * ξ2x3 +
            g[ijk, e, 1, 3] * ξ3x3

        vgrad_data[ijk, ∂₁u₂, e] =
            g[ijk, e, 2, 1] * ξ1x1 +
            g[ijk, e, 2, 2] * ξ2x1 +
            g[ijk, e, 2, 3] * ξ3x1
        vgrad_data[ijk, ∂₂u₂, e] =
            g[ijk, e, 2, 1] * ξ1x2 +
            g[ijk, e, 2, 2] * ξ2x2 +
            g[ijk, e, 2, 3] * ξ3x2
        vgrad_data[ijk, ∂₃u₂, e] =
            g[ijk, e, 2, 1] * ξ1x3 +
            g[ijk, e, 2, 2] * ξ2x3 +
            g[ijk, e, 2, 3] * ξ3x3

        vgrad_data[ijk, ∂₁u₃, e] =
            g[ijk, e, 3, 1] * ξ1x1 +
            g[ijk, e, 3, 2] * ξ2x1 +
            g[ijk, e, 3, 3] * ξ3x1
        vgrad_data[ijk, ∂₂u₃, e] =
            g[ijk, e, 3, 1] * ξ1x2 +
            g[ijk, e, 3, 2] * ξ2x2 +
            g[ijk, e, 3, 3] * ξ3x2
        vgrad_data[ijk, ∂₃u₃, e] =
            g[ijk, e, 3, 1] * ξ1x3 +
            g[ijk, e, 3, 2] * ξ2x3 +
            g[ijk, e, 3, 3] * ξ3x3
    end
end

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

"""
    Vorticity{
        FT <: AbstractFloat,
        FTA2D <: AbstractArray{FT, 2},
        FTA3D <: AbstractArray{FT, 3},
    }

This data structure stores the vorticity of a velocity field.

# Fields

$(DocStringExtensions.FIELDS)

# Usage

    Vorticity(data)

# Arguments for the inner constructor
 - `data`: 3-dimensional device array containing the vorticity data
   (the second dimension must be 3)
"""
struct Vorticity{
    FT <: AbstractFloat,
    FTA2D <: AbstractArray{FT, 2},
    FTA3D <: AbstractArray{FT, 3},
}
    "Device array storing the vorticity data"
    data::FTA3D
    "View of x1 component of vorticity"
    Ω₁::FTA2D
    "View of x2 component of vorticity"
    Ω₂::FTA2D
    "View of x3 component of vorticity"
    Ω₃::FTA2D
    function Vorticity(data::AbstractArray{FT, 3}) where {FT <: AbstractFloat}
        Ω₁ = view(data, :, 1, :)
        Ω₂ = view(data, :, 2, :)
        Ω₃ = view(data, :, 3, :)
        return new{FT, typeof(Ω₁), typeof(data)}(data, Ω₁, Ω₂, Ω₃)
    end
end

"""
    Vorticity(
        dg::DGModel,
        vgrad::VectorGradients,
    )

This function computes the vorticity of the velocity field.

# Arguments
 - `dg`: DGModel
 - `vgrad`: vector gradients
"""
function Vorticity(dg::DGModel, vgrad::VectorGradients)
    bl = dg.balance_law
    FT = eltype(dg.grid)
    N = polynomialorder(dg.grid)
    Nq = N + 1
    npoints = Nq^3
    nrealelem = length(dg.grid.topology.realelems)

    data = similar(vgrad.data, npoints, 3, nrealelem)

    device = array_device(data)
    workgroup = (Nq, Nq)
    ndrange = (nrealelem * Nq, Nq)

    kernel = vorticity_kernel!(device, workgroup)
    event = kernel(vgrad.data, data, Val(Nq), ndrange = ndrange)
    wait(event)

    return Vorticity(data)
end

@kernel function vorticity_kernel!(
    vgrad_data::AbstractArray{FT, 3},
    vort_data::AbstractArray{FT, 3},
    ::Val{qm1},
) where {qm1, FT <: AbstractFloat}

    e = @index(Group, Linear)
    i, j = @index(Local, NTuple)

    ∂₁u₁, ∂₂u₁, ∂₃u₁ = 1, 2, 3
    ∂₁u₂, ∂₂u₂, ∂₃u₂ = 4, 5, 6
    ∂₁u₃, ∂₂u₃, ∂₃u₃ = 7, 8, 9
    Ω₁, Ω₂, Ω₃ = 1, 2, 3

    for k in 1:qm1
        ijk = i + ((j - 1) + (k - 1) * qm1) * qm1

        vort_data[ijk, Ω₁, e] =
            vgrad_data[ijk, ∂₂u₃, e] - vgrad_data[ijk, ∂₃u₂, e]
        vort_data[ijk, Ω₂, e] =
            vgrad_data[ijk, ∂₃u₁, e] - vgrad_data[ijk, ∂₁u₃, e]
        vort_data[ijk, Ω₃, e] =
            vgrad_data[ijk, ∂₁u₂, e] - vgrad_data[ijk, ∂₂u₁, e]
    end
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
#=
julia --color=yes make.jl
=#
using Documenter, LSMH

makedocs()


text:
function levelprune(w::AbstractArray, level)
    p = sortperm(vec(w); rev = true)
    n = ceil(Int, (1 - level) * length(p))

    return mask(w, p[1:n])
end

struct LevelPrune{T<:Real}
    level::T
end
(p::LevelPrune)(w) = levelprune(w, p.level)

thresholdprune(w::AbstractArray, threshold) = mask(w, abs(w) .>= threshold)

struct ThresholdPrune{T<:Real}
    threshold::T
end
(p::ThresholdPrune)(w) = thresholdprune(w, p.threshold)


text:
function frac_cov(t, d)
  dsave = d
  if d < -0.5;
   error("d must be > -0.5");
  end;
  if d > 1.5;
     error("d must be < 1.5");
  end;
  (d>0.50) && (d = d - 1)
  (d == 0) && return covmat = eye(t)
  (abs(d) == 0.5) && (d = 0.49999*sign(d))
  acv = zeros(t, 1)
  acv[1] = 1
  (d > 0) && (Γmd = gamma(d))
  (d < 0) && (Γmd = gamma(1+d)/d)
  Γ_md = gamma(1-d)
  fac = Γ_md/Γmd
  t1 = minimum([100, t-1])
  h = collect(1:t1)
  Γhd = gamma(h+d)
  Γhmd = gamma(h-d+1)
  acv[2:t1+1] = (Γhd./Γhmd)*fac
  if t1+1<t
    hvec = collect(t1+1:t-1)
    anum = lgamma(d+hvec)
    aden = lgamma(hvec-d+1)
    acv[t1+2:t] = fac*exp(anum-aden)
  end
  covmat = acv2cov(acv)
  var_scl=gamma(1-2*d)/(gamma(1-d))^2;
  covmat=covmat*var_scl;
  if dsave > .50
    sum_mat = tril(ones(t,t))
    covmat = sum_mat*covmat*sum_mat'
  end
   return covmat
end


text:
using Observables

abstract type AbstractWidget{T} <: Observables.AbstractObservable{T} end

Observables.observe(x::AbstractWidget) = x.value

struct Button{T} <: AbstractWidget{Bool}
    content::Observable{T}
    value::Observable{Bool}
    attributes::Dict{Symbol, Any}
end

function Button(content; kw...)
    return Button(
        Observable(content), Observable(false), Dict{Symbol, Any}(kw)
    )
end

function jsrender(button::Button)
    return DOM.input(
        type = "button",
        value = button.content,
        onclick = js"update_obs($(button.value), true);";
        button.attributes...
    )
end

struct TextField <: AbstractWidget{String}
    value::Observable{String}
    attributes::Dict{Symbol, Any}
end

function TextField(value::String; kw...)
    TextField(Observable(value), Dict{Symbol, Any}(kw))
end

function jsrender(tf::TextField)
    return DOM.input(
        type = "textfield",
        value = tf.value,
        onchange = js"update_obs($(tf.value),  this.value);";
        tf.attributes...
    )
end

struct NumberInput <: AbstractWidget{Float64}
    value::Observable{Float64}
    attributes::Dict{Symbol, Any}
end

function NumberInput(value::Float64; kw...)
    NumberInput(Observable(value), Dict{Symbol, Any}(kw))
end

function jsrender(ni::NumberInput)
    return DOM.input(
        type = "number",
        value = ni.value,
        onchange = js"update_obs($(ni.value), parseFloat(this.value));";
        ni.attributes...
    )
end

struct Slider{T <: AbstractRange, ET} <: AbstractWidget{T}
    range::Observable{T}
    value::Observable{ET}
    attributes::Dict{Symbol, Any}
end

to_node(x) = Observable(x)
to_node(x::Observable) = x

function Slider(range::T, value = first(range); kw...) where T <: AbstractRange
    Slider{T, eltype(range)}(
        to_node(range),
        to_node(value),
        Dict{Symbol, Any}(kw)
    )
end

function jsrender(slider::Slider)
    return DOM.input(
        type = "range",
        min = map(first, slider.range),
        max = map(last, slider.range),
        value = slider.value,
        step = map(step, slider.range),
        oninput = js"update_obs($(slider.value), parseFloat(value))";
        slider.attributes...
    )
end

@enum Orientation vertical horizontal

struct RangeSlider{T <: AbstractRange, ET <: AbstractArray} <: AbstractWidget{T}
    attributes::Dict{Symbol, Any}
    range::Observable{T}
    value::Observable{ET}
    connect::Observable{Bool}
    orientation::Observable{Orientation}
    tooltips::Observable{Bool}
    ticks::Observable{Dict{String, Any}}
end

function RangeSlider(range::T; value = [first(range)], kw...) where T <: AbstractRange
    RangeSlider{T, typeof(value)}(
        Dict{Symbol, Any}(kw),
        Observable(range),
        Observable(value),
        Observable(true),
        Observable(horizontal),
        Observable(false),
        Observable(Dict{String, Any}()),
    )
end

const noUiSlider = Dependency(
    :noUiSlider,
    [
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.js",
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.css"
    ]
)

function jsrender(session::Session, slider::RangeSlider)
    args = (slider.range, slider.connect, slider.orientation,
            slider.tooltips, slider.ticks)
    style = map(args...) do range, connect, orientation, tooltips, ticks
        value = slider.value[]
        return Dict(
            "range" => Dict(
                "min" => fill(range[1], length(value)),
                "max" => fill(range[end], length(value)),
            ),
            "step" => step(range),
            "start" => value,
            "connect" => connect,
            "tooltips" => tooltips,
            "direction" => "ltr",
            "orientation" => string(orientation),
            "pips" => isempty(ticks) ? nothing : ticks
        )
    end
    rangediv = DOM.div()
    create_slider = js"""function create_slider(style){
        var range = $(rangediv);
        console.log(style);
        range.noUiSlider.updateOptions(deserialize_js(style), true);
    }"""
    onload(session, rangediv, js"""function onload(range){
        var style = $(style[]);
        $(noUiSlider).create(range, style);
        range.noUiSlider.on('update', function (values, handle, unencoded, tap, positions){
            update_obs($(slider.value), [parseFloat(values[0]), parseFloat(values[1])]);
        });
    }""")
    onjs(session, style, create_slider)
    return rangediv
end


text:
ah1_file  = path*"/SampleFiles/AH/ah1.f"
ah1_fstr  = path*"/SampleFiles/AH/ah1.*"
ahc_file  = path*"/SampleFiles/AH/lhz.ah"
ah_resp   = path*"/SampleFiles/AH/BRV.TSG.DS.lE21.resp"
ah2_file  = path*"/SampleFiles/AH/ah2.f"
ah2_fstr  = path*"/SampleFiles/AH/ah2.*"

printstyled("  AH (Ad Hoc)\n", color=:light_green)

printstyled("    v1\n", color=:light_green)
redirect_stdout(out) do
  S = read_data("ah1", ah1_file, v=3)
  S = read_data("ah1", ah1_file, full=true)
  S = read_data("ah1", ah1_fstr, full=true, v=3)
  @test S.n == 4
  @test S.fs[1] == 4.0
  @test isapprox(S.gain[1], 64200.121094)
  @test isapprox(S.loc[1].lat, 35.599899)
  @test isapprox(S.loc[1].lon, -85.568802)
  @test isapprox(S.loc[1].el, 481.0)
  @test any([occursin("gdsn_tape",s) for s in S.notes[1]])
  @test any([occursin("demeaned",s) for s in S.notes[1]])
  @test length(S.resp[1].p) == 24
  @test length(S.resp[1].z) == 7
  @test u2d(S.t[1][1,2]*1.0e-6) == DateTime("1984-04-20T06:42:00.12")
  @test length(S.x[1]) == 720
  @test eltype(S.x[1]) == Float32
  @test isapprox(S.x[1][1:4], [-731.41247559, -724.41247559, -622.41247559, -470.4125061])

  C = read_data("ah1", ahc_file, v=3, full=true)[1]

  # Station
  @test isapprox(C.loc.lat, 36.5416984)
  @test isapprox(C.loc.lon, 138.2088928)
  @test isapprox(C.loc.el, 422.0)
  @test isapprox(C.gain, 1.178e8)
  @test length(C.resp.p) == 10
  @test isapprox(C.resp.p[1], -0.0123f0 + 0.0123f0im)
  @test isapprox(C.resp.p[2], -0.0123f0 - 0.0123f0im)
  @test length(C.resp.z) == 3

  # Data
  @test length(C.x) == 309
  @test eltype(C.x) == Float32
  @test C.fs == 1.0
  @test u2d(C.t[1,2]*1.0e-6) == DateTime("1990-05-12T04:49:54.49")

  # Event
  @test isapprox(C.misc["ev_lat"], 49.037)
  @test isapprox(C.misc["ev_lon"], 141.847)
  @test isapprox(C.misc["ev_dep"], 606.0)
  @test string(u2d(C.misc["ot"]*1.0e-6)) == "1990-05-12T04:50:08.7"
  @test startswith(C.misc["data_comment"], "Streckeisen STS-1V/VBB Seismometer")
  @test startswith(C.misc["event_comment"], "null")

  C = read_data("ah1", ah_resp, full=true)[1]
  @test isapprox(C.loc.lat, 53.058060)
  @test isapprox(C.loc.lon, 70.282799)
  @test isapprox(C.loc.el, 300.0)
  @test isapprox(C.gain, 0.05)
  @test isapprox(C.resp.a0, 40.009960)
  @test length(C.resp.p) == 7
  @test isapprox(C.resp.p[1], -0.1342653f0 + 0.1168836f0im)
  @test length(C.resp.z) == 4
  @test startswith(C.misc["data_comment"], "DS response in counts/nm")
  @test startswith(C.misc["event_comment"], "Calibration_for_hg_TSG")
  @test any([occursin("brv2ah: ahtedit",s) for s in C.notes])
  @test any([occursin("demeaned",s) for s in C.notes])
  @test any([occursin("modhead",s) for s in C.notes])
  @test any([occursin("ahtedit",s) for s in C.notes])

  printstyled("    v2\n", color=:light_green)
  S = read_data("ah2", ah2_file, v=3)
  S = read_data("ah2", ah2_file, v=3, full=true)
  S = read_data("ah2", ah2_fstr, v=3, full=true)
  @test S.n == 4
  @test S.fs[1] == 4.0
  @test isapprox(S.gain[1], 64200.121094)
  @test isapprox(S.loc[1].lat, 35.599899)
  @test isapprox(S.loc[1].lon, -85.568802)
  @test isapprox(S.loc[1].el, 481.0)
  @test any([occursin("gdsn_tape",s) for s in S.notes[1]])
  @test any([occursin("demeaned",s) for s in S.notes[1]])
  @test length(S.resp[1].p) == 24
  @test length(S.resp[1].z) == 7
  @test u2d(S.t[1][1,2]*1.0e-6) == DateTime("1984-04-20T06:42:00.12")
  @test length(S.x[1]) == 720
  @test eltype(S.x[1]) == Float32
  @test isapprox(S.x[1][1:4], [-731.41247559, -724.41247559, -622.41247559, -470.4125061])

end


text:
# g_memory_input_stream_new_from_data ()
# GInputStream *
# g_memory_input_stream_new_from_data (const void *data,
#                                      gssize len,
#                                      GDestroyNotify destroy);
# Creates a new GMemoryInputStream with data in memory of a given size.
# Parameters
# data
# input data.
# 	[array length=len][element-type guint8][transfer full]
# len
# length of the data, may be -1 if data is a nul-terminated string
# destroy
# function that is called to free data , or NULL.
# 	[nullable]
# Returns
# new GInputStream read from data of len bytes.

mutable struct GInputStream
	ptr::Ptr{Nothing}

    function GInputStream(ptr::Ptr{Nothing})
        self = new(ptr)
        finalizer(destroy, self)
        self
    end
end

function destroy(handle::GInputStream)
    if handle.ptr == C_NULL
        return
    end
    handle_free(handle)
    handle.ptr = C_NULL
    nothing
end

function glib_memory_input_stream_new_from_data(data::AbstractString)
	    ptr = ccall((:g_memory_input_stream_new_from_data, libgio), Ptr{Nothing},
                (Ptr{UInt8},UInt32,Ptr{Nothing}), string(data), sizeof(data),C_NULL)
    	GInputStream(ptr)
end


text:
# ---
# title: 36. Valid Sudoku
# id: problem36
# author: zhwang
# date: 2022-01-23
# difficulty: Medium
# categories: Hash Table
# link: <https://leetcode.com/problems/valid-sudoku/description/>
# hidden: true
# ---
#
# Determine if a `9 x 9` Sudoku board is valid. Only the filled cells need to be
# validated  **according to the following rules** :
#
#   1. Each row must contain the digits `1-9` without repetition.
#   2. Each column must contain the digits `1-9` without repetition.
#   3. Each of the nine `3 x 3` sub-boxes of the grid must contain the digits `1-9` without repetition.
#
# **Note:**
#
#   * A Sudoku board (partially filled) could be valid but is not necessarily solvable.
#   * Only the filled cells need to be validated according to the mentioned rules.
#
#
#
# **Example 1:**
#
# ![](https://upload.wikimedia.org/wikipedia/commons/thumb/f/ff/Sudoku-
# by-L2G-20050714.svg/250px-Sudoku-by-L2G-20050714.svg.png)
#
#
#
#     Input: board =
#     [["5","3",".",".","7",".",".",".","."]
#     ,["6",".",".","1","9","5",".",".","."]
#     ,[".","9","8",".",".",".",".","6","."]
#     ,["8",".",".",".","6",".",".",".","3"]
#     ,["4",".",".","8",".","3",".",".","1"]
#     ,["7",".",".",".","2",".",".",".","6"]
#     ,[".","6",".",".",".",".","2","8","."]
#     ,[".",".",".","4","1","9",".",".","5"]
#     ,[".",".",".",".","8",".",".","7","9"]]
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: board =
#     [["8","3",".",".","7",".",".",".","."]
#     ,["6",".",".","1","9","5",".",".","."]
#     ,[".","9","8",".",".",".",".","6","."]
#     ,["8",".",".",".","6",".",".",".","3"]
#     ,["4",".",".","8",".","3",".",".","1"]
#     ,["7",".",".",".","2",".",".",".","6"]
#     ,[".","6",".",".",".",".","2","8","."]
#     ,[".",".",".","4","1","9",".",".","5"]
#     ,[".",".",".",".","8",".",".","7","9"]]
#     Output: false
#     Explanation: Same as Example 1, except with the **5** in the top left corner being modified to **8**. Since there are two 8's in the top left 3x3 sub-box, it is invalid.
#
#
#
#
# **Constraints:**
#
#   * `board.length == 9`
#   * `board[i].length == 9`
#   * `board[i][j]` is a digit or `'.'`.
#
#
## @lc code=start
using LeetCode

function is_valid_sudoku(board::Vector{Vector{Char}})::Bool
    function check_line(i::Int)::Bool
        s = replace(join(board[i]), "." => "")
        return length(s) == length(Set{Char}(s))
    end
    function check_col(i::Int)::Bool
        s = replace(join(board[k][i] for k in 1:9), '.' => "")
        return length(s) == length(Set{Char}(s))
    end
    function check_box(i::Int, j::Int)::Bool
        s = replace(join(board[i + p1][j + p2] for p1 in 1:3 for p2 in 1:3), "." => "")
        return length(s) == length(Set{Char}(s))
    end
    return all(check_line(i) && check_col(i) for i in 1:9) &&
           all(check_box(3 * (i - 1), 3 * (j - 1)) for i in 1:3 for j in 1:3)
end

## @lc code=end


text:
using ProteinStructureReconstruction
using Test

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


text:
using Akagera
using Base.Test

@testset "Animator tests" begin
    @testset "Animator(interval::Int, frames::Vector{String}, animate_type::String)" begin
        interval = 1
        frames = ["T", "e", "s", "t", "!"]
        animate_type = "Test"

        @test try
            Animator(interval, frames, animate_type)
            true
        catch
            false
        end

        a = Animator(interval, frames, animate_type)

        @test a.flag == false
        @test a.interval == interval
        @test_throws UndefRefError a.msg
        @test a.frames == frames
        @test a.animate_type == animate_type
    end

    @testset "Animator(json::String; animate_type)" begin
        json = """
        {
            "interval": 2018,
            "frames": ["T", "e", "s", "t", "!"]
        }
        """

        @test try
            Animator(json, animate_type="linear")
            true
        catch
            false
        end

        a = Animator(json, animate_type="linear")

        @test a.frames == ["T", "e", "s", "t", "!"]
        @test a.interval == 2018

        @test_warn "animate_type is not specified. Default value `linear` will be set." Animator(json)

        json_without_interval = """
        {
            "frames": ["T", "e", "s", "t", "!"]
        }
        """

        @test_throws ArgumentError Animator(json_without_interval, animate_type="linear")

        json_without_frames = """
        {
            "interval": 2018
        }
        """

        @test_throws ArgumentError Animator(json_without_frames, animate_type="linear")
    end

    @testset "Animator(fc::AnimeContainer)" begin
        fc = arrow3

        @test try
            Animator(fc)
            true
        catch
            false
        end

        a = Animator(fc)

        @test a.flag == false
        @test a.interval == fc.interval
        @test_throws UndefRefError a.msg
        @test a.frames == fc.frames
        @test a.animate_type == fc.animate_type
    end

    @testset "init!()" begin
        interval = 1
        frames = ["T", "e", "s", "t", "!"]
        animate_type = "Test"

        a = Animator(interval, frames, animate_type)

        @test a.flag == false

        @test try
            init!(a)
            true
        catch
            false
        end

        @test a.flag == true
    end

    @testset "set_msg!()" begin
        interval = 1
        frames = ["T", "e", "s", "t", "!"]
        animate_type = "Test"

        a = Animator(interval, frames, animate_type)

        @test_throws UndefRefError a.msg

        @test try
            set_msg!(a, "TestMsg")
            true
        catch
            false
        end

        @test a.msg == "TestMsg"
    end

    @testset "start!()" begin
        @testset "render()" begin
            interval = 1
            frames = ["T", "e", "s", "t", "!"]
            animate_type = "Test"

            a = Animator(interval, frames, animate_type)

            set_msg!(a, "TestMsg")

            @test a.frames[1] * " " * a.msg == "T TestMsg"
        end

        a = Animator(dots)

        @test try
            start!(a)
            finish!(a)

            true
        catch
            false
        end

        linear = map(x->x, [1:5;])
        @test linear == [1, 2, 3, 4, 5]

        swing = map(x->
            begin
                if x <= 5
                    x
                else
                    10-x
                end
            end, [1:8;]
        )

        @test swing == [1, 2, 3, 4, 5, 4, 3, 2]
    end

    @testset "update_msg!()" begin
        interval = 1
        frames = ["T", "e", "s", "t", "!"]
        animate_type = "Test"

        a = Animator(interval, frames, animate_type)

        set_msg!(a, "TestMsg")

        @test a.msg == "TestMsg"

        @test try
            update_msg!(a, "TestUpdatingMsg")
            true
        catch
            false
        end

        @test a.msg == "TestUpdatingMsg"
    end

    @testset "finish!()" begin
        interval = 1
        frames = ["T", "e", "s", "t", "!"]
        animate_type = "Test"

        a = Animator(interval, frames, animate_type)

        init!(a)

        @test a.flag == true

        @test try
            finish!(a)
            true
        catch
            false
        end

        @test a.flag == false
    end
end


text:
using Format
using Compat.Test
using Compat.Printf
using Compat.Random

@testset "cformat"    begin include( "cformat.jl" ) end
@testset "fmtspec"    begin include( "fmtspec.jl" ) end
@testset "formatexpr" begin include( "formatexpr.jl" ) end
@testset "fmt"        begin include( "fmt.jl" ) end


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
mutable struct vertex
    layer::Int #index of the layer in the graph
    i::Int  #index of the variable
    w::Int  #accumulated weight
    zλ::Int #best profit possible for zλ
    z1::Int #best profit possible for z1
    z2::Int #best profit possible for z2
    has_parent_0::Bool
    has_parent_1::Bool
    parent_0::vertex #parent node, with an edge having a profit == 0
    parent_1::vertex #parent node, with an edge having a profit > 0
    vertex(l,i,w,zλ,z1,z2,hp0,hp1) = new(l,i,w,zλ,z1,z2,hp0,hp1)
    vertex(l,i,w,zλ,z1,z2,hp0,hp1,p) = new(l,i,w,zλ,z1,z2,hp0,hp1,p,p)
end

# Base.show(io::IO, v::vertex) = print(io, "v_",v.i,"_",v.w,"(",v.zλ,")(",v.layer,")")
Base.show(io::IO, v::vertex) = print(io, "v_",v.i,"_",v.w,"(",v.z1,",",v.z2,")(",v.layer,")")

#v_i_0
source(mono_pb::mono_problem) = vertex(1,mono_pb.variables[1],mono_pb.ω,mono_pb.min_profit_λ,mono_pb.min_profit_1,mono_pb.min_profit_2,false,false)

#inner degree of a vertex
inner_degree(v::vertex) = v.has_parent_0 + v.has_parent_1

#Creates a vertex v_i+1_w from a vertex v_i_w assuming we decided not to pick item i
function vertex_skip(v::vertex, mono_pb::mono_problem)
    varplus1 = v.layer == size(mono_pb) ? v.layer + 1 : mono_pb.variables[v.layer+1]
    return vertex(v.layer+1, varplus1, v.w, v.zλ, v.z1, v.z2, true, false, v)
end

#Creates a vertex v_i+1_w' from a vertex v_i_w assuming we decided to pick item i
function vertex_keep(v::vertex, mono_pb::mono_problem)
    #@assert v.mono_pb.variables[v.layer] == v.i
    var = v.i
    varplus1 = v.layer == size(mono_pb) ? v.layer + 1 : mono_pb.variables[v.layer+1]
    return vertex(v.layer+1,
            varplus1,
            v.w+mono_pb.w[var],
            v.zλ+mono_pb.p[var],
            v.z1+mono_pb.p1[var],
            v.z2+mono_pb.p2[var],
            false, true,
            v)
end

#merge two vertices
function merge!(a::vertex, b::vertex)
    if b.zλ >= a.zλ
        a.zλ, a.z1, a.z2 = b.zλ, b.z1, b.z2
    end
    a.has_parent_1 = true
    a.parent_1 = b.parent_1
end

#Comparison functions for searchsorted()
weight_lt(v::vertex, w::Int) = v.w < w
weight_lt(w::Int, v::vertex) = w < v.w

#Returns the unique parent of a vertex
function parent(v::vertex)
    return v.has_parent_0 ? v.parent_0 : v.parent_1
end

#Returns both parents of a vertex
function parents(v::vertex)
    return v.parent_0, v.parent_1
end

zλ(v::vertex) = v.zλ
z1(v::vertex) = v.z1
z2(v::vertex) = v.z2

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
import ..SamplingConstraints: SamplingConstraint
import ..UncertainValues: CertainValue

resample(x::Number) = x

resample(v::CertainValue) = v.value
resample(v::CertainValue, n::Int) = [v.value for i = 1:n]

resample(v::CertainValue, s::SamplingConstraint) = v.value
resample(v::CertainValue, s::SamplingConstraint, n::Int) = [v.value for i = 1:n]

constraints = [
    :(NoConstraint),
    :(TruncateLowerQuantile),
    :(TruncateUpperQuantile),
    :(TruncateQuantiles),
    :(TruncateMaximum),
    :(TruncateMinimum),
    :(TruncateRange),
    :(TruncateStd)
]

for constraint in constraints
    funcs = quote
        resample(x::CertainValue, constraint::$(constraint)) = x.value
        resample(x::CertainValue, constraint::$(constraint), n::Int) = [x.value for i = 1:n]
    end
    eval(funcs)
end


text:
# This file was generated, do not modify it. # hide
import MLJ: schema, std, mean, median, coerce, coerce!, scitype
using DataFrames
using UrlDownload
using PyPlot
ioff() # hide


text:
using CCBlade

import ForwardDiff, ReverseDiff
import SparseArrays, SparseDiffTools, FiniteDiff

import Statistics: median
using BenchmarkTools

# ------- setup problem for a number of flight states (run once) --------

function setup(nV)

    chord = 0.10
    D = 1.6
    RPM = 2100
    pitch = 1.0  # pitch distance in meters.

    turbine = false
    Rhub = 0.01
    Rtip = D/2
    B = 2  # number of blades

    R = D/2.0
    n = 11
    r = range(R/10, stop=9/10*R, length=n)
    theta = atan.(pitch./(2*pi*r))
    chord = chord*ones(n)

    function affunc(alpha, Re, M)

        cl = 6.2*alpha
        cd = 0.008 - 0.003*cl + 0.01*cl*cl

        return cl, cd
    end

    pitch = 0.0
    precone = 0.0

    if nV == 1
        Vinf = [30.0]
    else
        Vinf = range(29.0, 31.0, length=nV)
    end

    Omega = RPM * pi/30 * ones(nV)
    rho = 1.225 * ones(nV)

    x = [r; chord; theta; Rhub; Rtip; pitch; precone; Vinf; Omega; rho]
    y = zeros(2*nV)

    # ------- standard setup -----------

    # parameters that passthrough: airfoils, B, turbine, n, nV
    function ccbladewrapper!(y, x)

        r = x[1:n]
        chord = x[n+1:2*n]
        theta = x[2*n+1:3*n]
        Rhub = x[3*n+1]
        Rtip = x[3*n+2]
        pitch = x[3*n+3]
        precone = x[3*n+4]
        idx = 3*n+4
        Vinf = x[idx+1:idx+nV]
        Omega = x[idx+nV+1:idx+2*nV]
        rho = x[idx+2*nV+1:idx+3*nV]

        rotor = Rotor(Rhub, Rtip, B, turbine=turbine, precone=precone)
        sections = Section.(r, chord, theta, affunc)
        ops = simple_op.(Vinf', Omega', r, rho', pitch=pitch)

        outputs = solve.(Ref(rotor), sections, ops)

        for i = 1:nV
            T, Q = thrusttorque(rotor, sections, outputs[:, i])
            y[i] = T
            y[i+nV] = Q
        end

        return nothing
    end

    return x, y, ccbladewrapper!
end

# -------- compare forward AD and sparse forward AD ---------

nVvec = [1, 2, 4, 8, 16, 32, 64, 128, 256]
timead = zeros(9)
timeadsparse = zeros(9)
timeFDsparse = zeros(9)

for ii = 1:9
    nV = nVvec[ii]

    x, y, func = setup(nV)
    config = ForwardDiff.JacobianConfig(func, y, x)  #, ForwardDiff.Chunk{nChunk}())
    J = zeros(length(y), length(x))

    t = @benchmark ForwardDiff.jacobian!($J, $func, $y, $x, $config)
    timead[ii] = median(t).time*1e-9

    Jsparse = SparseArrays.sparse(J)
    colors = SparseDiffTools.matrix_colors(Jsparse)
    cache = SparseDiffTools.ForwardColorJacCache(func, x, dx=y, colorvec=colors, sparsity=Jsparse)  # nChunk,

    t = @benchmark SparseDiffTools.forwarddiff_color_jacobian!($Jsparse, $func, $x, $cache)
    timeadsparse[ii] = median(t).time*1e-9

    JsparseFD = copy(Jsparse)
    sparsecache = FiniteDiff.JacobianCache(x, Val{:forward}, colorvec=colors, sparsity=JsparseFD)

    t = @benchmark FiniteDiff.finite_difference_jacobian!($JsparseFD, $func, $x, $sparsecache)
    timeFDsparse[ii] = median(t).time*1e-9

    println(nV)
end

using PyPlot
rc("axes.spines", right=false, top=false)
rc("font", size=14.0)
rc("legend", frameon=false)

figure()
plot(nVvec, timead, "-o", color="#348ABD")
plot(nVvec, timeadsparse, "-o", color="#A60628")
xlabel("\\# inflow conditions")
ylabel("Jacobian time (s)")
legend(["AD", "AD w/ coloring"])
# savefig("adscaling.pdf")

# NOTE: the finite difference with coloring speed has improved since the paper was submitted.
figure()
plot(nVvec, timeFDsparse, "-o", color="#348ABD")
plot(nVvec, timeadsparse, "-o", color="#A60628")
xlabel("\\# inflow conditions")
ylabel("Jacobian time (s)")
legend(["FD w/ coloring", "AD w/ coloring"])
# savefig("adfd.pdf")

timead ./ timeadsparse
timeFDsparse ./ timeadsparse


text:
module K8s

include("client.jl")
include("informer.jl")
include("crd.jl")
include("core.jl")

export TaskListWatcher, TaskInformer, run, K8sClient, request, client

end # module

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
module FasterMinMax

export unsafe_min, unsafe_max, unsafe_minmax, unsafe_maxmin,
       minnum, maxnum, minmaxnum, maxminnum,
       minmag, maxmag, minmaxmag, maxminmag,

const FastFloat = Union{Float32, Float64}

function unsafe_min(x::T, y::T) where {T<:FastFloat}
    x_y = x - y
    (signbit(x_y) || isnan(x)) ? x : y
end

function unsafe_max(x::T, y::T) where {T<:FastFloat}
    y_x = y - x
    (signbit(y_x) || isnan(x)) ? x : y
end

function unsafe_minmax(x::T, y::T) where {T<:FastFloat}
    (unsafe_min(x,y), unsafe_max(x,y))
end

function unsafe_maxmin(x::T, y::T) where {T<:FastFloat}
    (signbit(y-x) || isnan(x)) ? (x,y) : (y,x)
end

function minnum(x::T, y::T) where {T<:FastFloat}
end

function maxnum(x::T, y::T) where {T<:FastFloat}
end

function minmaxnum(x::T, y::T) where {T<:FastFloat}
end

function maxminnum(x::T, y::T) where {T<:FastFloat}
end

function minmag(x::T, y::T) where {T<:FastFloat}
end

function maxmag(x::T, y::T) where {T<:FastFloat}
end

function minmaxmag(x::T, y::T) where {T<:FastFloat}
end

function maxminmag(x::T, y::T) where {T<:FastFloat}
end

end  # FasterMinMax


text:
"""
    delta(t::Real, lastspike, q)
    delta(t::Real, lastspike::AbstractArray{<:Real}, q::AbstractArray{<:Real})
    delta(t::Real, lastspike::CuVecOrMat{<:Real}, q::CuVecOrMat{<:Real})

Evaluate a Dirac-delta synapse.
Use `CuVector` instead of `Vector` for GPU support.

# Fields
- `t`: current time
- `lastspike`: last pre-synaptic spike time
- `q`: amplitude
"""
delta(t::Real, lastspike, q) = (t ≈ lastspike) * q
delta(t::Real, lastspike::AbstractArray{<:Real}, q::AbstractArray{<:Real}) = (t .≈ lastspike) .* q
delta(t::Real, lastspike::CuVecOrMat{<:Real}, q::CuVecOrMat{<:Real}) = (t .≈ lastspike) .* q


text:
function add_dr!(
    df;
    idlabel = :TrackID,
    xlabel = :POSITION_X,
    ylabel = :POSITION_Y
)
    df[!, :dx] .= prepend!(diff(df[:, xlabel]), NaN)
    df[!, :dy] .= prepend!(diff(df[:, ylabel]), NaN)
    df.dr2 = abs2.(df.dx) + abs2.(df.dy)
    df.dr = sqrt.(df.dr2)

    # df[df.FRAME2 .== 1, [:dx, :dy, :dr, :dr2]] .= NaN
    @inbounds for i = 2:nrow(df)
        if df[i, idlabel] != df[i-1, idlabel]
            df[i, [:dx, :dy, :dr, :dr2]] .= NaN
        end
    end
end

function dr2matrix(df, idlabel, datalabel, framelabel)
    N = length(Set(df[!, idlabel]))
    dr = Matrix{Union{Nothing,Float64}}(nothing, maximum(df[!, framelabel]), N)
    @inbounds for (i, n) = enumerate(sort!(collect(Set(df[!, idlabel]))))
        m = filter(!isnan, df[df[!, idlabel] .== n, datalabel])
        @simd for t = 1:length(m)
            dr[t, i] = m[t]
        end
    end
    dr
end

function add_s!(df, z; idlabel)
    df[!, :state] .= 0.0
    @inbounds for (i, n) = enumerate(sort!(collect(Set(df[!, idlabel]))))
        T = size(df[df[!, idlabel] .== n, :], 1)
        m = Vector{Any}(filter(!isnothing, z[1:T-1, i]))
        prepend!(m, NaN)
        df[df[:, idlabel] .== n, :state] .= m
    end
end

function data2matrix(
    df;
    idlabel = :TrackID,
    xlabel = :POSITION_X,
    ylabel = :POSITION_Y,
    framelabel = :FRAME2,
    δ::Int = 1
)
    N = length(Set(df[!, idlabel]))
    data = Matrix{Union{Nothing,Float64}}(undef, maximum(df[!, framelabel]) - δ, N)
    @inbounds for (i, n) = enumerate(sort!(collect(Set(df[!, idlabel]))))
        m = spt.extract(df, n, idlabel, [xlabel, ylabel])
        for Δt = 1:size(m, 1) - δ
            data[Δt, i] = displacement(m, Δt, δ)
        end
    end
    data
end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
using AtomicUnits
using Test

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


text:
# Linking of different PTX modules

export
    CuLink, CuLinkImage, add_data!, add_file!, complete

"""
    CuLink()

Creates a pending JIT linker invocation.
"""
mutable struct CuLink
    handle::CUlinkState
    ctx::CuContext

    options::Dict{CUjit_option,Any}
    optionKeys::Vector{CUjit_option}
    optionVals::Vector{Ptr{Cvoid}}

    function CuLink(options::Dict{CUjit_option,Any}=Dict{CUjit_option,Any}())
        handle_ref = Ref{CUlinkState}()

        options[JIT_ERROR_LOG_BUFFER] = Vector{UInt8}(undef, 1024*1024)
        @debug begin
            options[JIT_INFO_LOG_BUFFER] = Vector{UInt8}(undef, 1024*1024)
            options[JIT_LOG_VERBOSE] = true
            "JIT compiling code" # FIXME: remove this useless message
        end
        if Base.JLOptions().debug_level == 1
            # XXX: does not apply to the linker
            options[JIT_GENERATE_LINE_INFO] = true
        elseif Base.JLOptions().debug_level >= 2
            options[JIT_GENERATE_DEBUG_INFO] = true
        end
        optionKeys, optionVals = encode(options)

        cuLinkCreate_v2(length(optionKeys), optionKeys, optionVals, handle_ref)

        ctx = current_context()
        obj = new(handle_ref[], ctx, options, optionKeys, optionVals)
        finalizer(unsafe_destroy!, obj)
        return obj
    end
end

function unsafe_destroy!(link::CuLink)
    @finalize_in_ctx link.ctx cuLinkDestroy(link)
end

Base.unsafe_convert(::Type{CUlinkState}, link::CuLink) = link.handle

Base.:(==)(a::CuLink, b::CuLink) = a.handle == b.handle
Base.hash(link::CuLink, h::UInt) = hash(link.handle, h)

function Base.show(io::IO, link::CuLink)
    print(io, "CuLink(")
    @printf(io, "%p", link.handle)
    print(io, ", ", link.ctx, ")")
end

"""
    add_data!(link::CuLink, name::String, code::String)

Add PTX code to a pending link operation.
"""
function add_data!(link::CuLink, name::String, code::String)
    data = unsafe_wrap(Vector{UInt8}, code)

    # there shouldn't be any embedded NULLs
    checked_data = Base.unsafe_convert(Cstring, data)

    res = unsafe_cuLinkAddData_v2(link, JIT_INPUT_PTX, pointer(checked_data), length(data),
                                  name, 0, C_NULL, C_NULL)
    if res == ERROR_NO_BINARY_FOR_GPU ||
       res == ERROR_INVALID_IMAGE ||
       res == ERROR_INVALID_PTX
        throw(CuError(res, unsafe_string(pointer(link.options[JIT_ERROR_LOG_BUFFER]))))
    elseif res != SUCCESS
        throw_api_error(res)
    end
end

"""
    add_data!(link::CuLink, name::String, data::Vector{UInt8}, type::CUjitInputType)

Add object code to a pending link operation.
"""
function add_data!(link::CuLink, name::String, data::Vector{UInt8})
    res = unsafe_cuLinkAddData_v2(link, JIT_INPUT_OBJECT, pointer(data), length(data),
                                  name, 0, C_NULL, C_NULL)
    if res == ERROR_NO_BINARY_FOR_GPU ||
       res == ERROR_INVALID_IMAGE ||
       res == ERROR_INVALID_PTX
        throw(CuError(res, unsafe_string(pointer(link.options[JIT_ERROR_LOG_BUFFER]))))
    elseif res != SUCCESS
        throw_api_error(res)
    end

    return
end

"""
    add_file!(link::CuLink, path::String, typ::CUjitInputType)

Add data from a file to a link operation. The argument `typ` indicates the type of the
contained data.
"""
function add_file!(link::CuLink, path::String, typ::CUjitInputType)
    cuLinkAddFile_v2(link, typ, path, 0, C_NULL, C_NULL)

    return nothing
end

"""
The result of a linking operation.

This object keeps its parent linker object alive, as destroying a linker destroys linked
images too.
"""
struct CuLinkImage
    data::Array{UInt8}
    link::CuLink
end

"""
    complete(link::CuLink)

Complete a pending linker invocation, returning an output image.
"""
function complete(link::CuLink)
    cubin_ref = Ref{Ptr{Cvoid}}()
    size_ref = Ref{Csize_t}()

    res = unsafe_cuLinkComplete(link, cubin_ref, size_ref)
    if res == CUDA_ERROR_NO_BINARY_FOR_GPU || res == CUDA_ERROR_INVALID_IMAGE
        options = decode(link.optionKeys, link.optionVals)
        throw(CuError(res, options[JIT_ERROR_LOG_BUFFER]))
    elseif res != SUCCESS
        throw_api_error(res)
    end

    @debug begin
        options = decode(link.optionKeys, link.optionVals)
        if isempty(options[JIT_INFO_LOG_BUFFER])
            """JIT info log is empty"""
        else
            """JIT info log:
               $(options[JIT_INFO_LOG_BUFFER])"""
        end
    end

    data = unsafe_wrap(Array, convert(Ptr{UInt8}, cubin_ref[]), size_ref[])
    return CuLinkImage(data, link)
end

"""
    CuModule(img::CuLinkImage, ...)

Create a CUDA module from a completed linking operation. Options from `CuModule` apply.
"""
CuModule(img::CuLinkImage, args...) = CuModule(img.data, args...)
CuModule(img::CuLinkImage, options::Dict{CUjit_option,Any}) = CuModule(img.data, options)


text:
using ArgParse2
using DataFrames
using DataStructures
using BioSequences
using CodecZlib
using BioAlignments
using DataFramesMeta
using Statistics
using st16SseqJuliaTools
using XAM

function main()
    parser = ArgumentParser(prog = "analyse_alnself.jl",
                        description = "Detect chimeras",
                        epilog = "Removing potential chimeras and repair snps")
    add_argument!(parser, "--inbam", "-b", help = "input bam ",type = String)
    add_argument!(parser, "--ref", "-r", help = "reference fasta",type = String)
    add_argument!(parser, "--ideal", "-i", help = "ideal contigs - for debugging",type = String,default="")
    add_argument!(parser, "--out", "-o", help = "Fasta file to write out cleaned up sequences",default="chimerasfree.faSTA",type = String)
    add_argument!(parser, "--min-seqid", "-s", help = "Minimum sequence identity in match",type = Float64, default=0.95)
    add_argument!(parser, "--min-covid", "-c", help = "Minimum sequence length fraction in  match",type = Float64, default=0.96)
    #add_argument!(parser, "--output", "-o", help = "File with names",type = String)
    args = parse_args(parser)
    ref = fasta_to_dict(args.ref)

    #Limits for real match ... mayve should be options
    minl = 300
    minf_len = args.min_covid
    minf_id = args.min_seqid # Identity of first 250 fragments of contigs
    reader = open(BAM.Reader,args.inbam) #,index = args.inbam*".bai")
    args.ideal == "" ? known_good_contigs = [] : known_good_contigs = collect(keys(fasta_to_dict(args.ideal)))
    pairs = Array{Tuple{String,String},1}()
    chosenseqs = OrderedDict{String,BioSequences.LongSequence{BioSequences.DNAAlphabet{4}}}()
    ct_good = 0
    ct_aligned = 0
    ct_total = 0
    for record in reader
        ct_total += 1

        if BAM.ismapped(record) && BAM.isprimary(record)  #&& BAM.tempname(record)=="6_674;size=1"
            ct_aligned += 1
            rn = BAM.refname(record)
	        tempname = BAM.tempname(record)
            pos = BAM.position(record)
            seq = BAM.sequence(record)
            rpos = BAM.rightposition(record)
            aln = BAM.alignment(record)
            l_total = length(seq)
            l_aligned = abs(rpos-pos)
            mutations = call_variants(record)
            mismct = 0
            for m in mutations
                if m.class == 'M'
                    mismct += 1
                    #mismct += maximum([length(m.reference),length(m.variant)])
                end
                if m.class == 'D'
                    l_aligned -= length(m.reference)
                end
            end
            fl = round(l_aligned/l_total,digits=2)
            fid = round(1 - (length(mutations)/l_aligned),digits=2)
            #println(fl," ",fid," ")
            good = (tempname in known_good_contigs)
            if l_aligned >= minl && ( (fid >= minf_id && fl >= minf_len) || (fid >= 0.99 && fl >= 0.95) )
                ct_good += 1
                chosenseqs[tempname] = ref[tempname]
            end
        end
    end

    close(reader)
    println("$ct_total in total sequences, $ct_aligned  were found in alignments,  $ct_good sequences were  written to ",args.out)
    dict_to_fasta(chosenseqs,args.out)
end

main()


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
function repeatlength(A, i)
    k = 1
    @inbounds while i + k ≤ lastindex(A) && isequal(A[i+k], A[i])
        k += 1
    end
    k
end

"""
    RepVal(A)
Iterator generating the all (repeat length, value)-tuples in `A`.
"""
struct RepVal{T} A::T end
Base.IteratorSize(::Type{RepVal{T}}) where T = Base.SizeUnknown()
Base.eltype(x::RepVal) = Tuple{Int, eltype(x.A)} # (rep.len., value)

Base.iterate(x::RepVal) = iterate(x, firstindex(x.A))
function Base.iterate(x::RepVal, i::Int)
    i > lastindex(x.A) && return nothing
    k = repeatlength(x.A, i)
    (k, x.A[i]), i + k
end

maxrep_maxval(A) = maximum(RepVal(A))
negrep((k, v)) = (-k, v)
maxrep_minval(A) = negrep(minimum(negrep, RepVal(A)))

@doc RepVal

# %%
A = [1, 2, 3, 3, 3, 1, 1, 1, NaN, NaN, NaN, 2, 2, 3]
@show A
RepVal(A) |> collect

# %%
maxrep_maxval(A)

# %%
maxrep_minval(A)

# %%
@code_warntype iterate(RepVal(A), 1)

# %% [markdown]
# https://discourse.julialang.org/t/minimum-mode-problem-minimum-number-maximum-repetition/66749/2

# %%
function findn(x,y=zeros(Int,maximum(x)))
    m = 0
    fill!(y,0)
    @inbounds for i in x
        y[i] += 1
        m = max(y[i],m)
    end
    return findfirst(isequal(m),y)
end

# %%

# %%
using StatsBase, DataStructures

function findn_jling(x)
    cm = sort!(OrderedDict(countmap(x)); byvalue = true)
    last(cm.keys), last(cm.vals)
end

# %%

# %%
function minMode(c)
    minVal = -1
    i = 0

    i = i+1
    for row in eachrow(c)
    row = collect(row)
    temp = filter(x -> x!=0, row)
    count = counter(temp)
    sortedCollection = sort(collect(count), by=x->x[2], rev=true)
    minVal = sortedCollection[1][1]
    repetitions = sortedCollection[1][2]
    for (key, value) in sortedCollection
        if(value == repetitions)
            if(key < minVal)
                minVal = key
            end
        end
        end
    end

    return minVal
end

# %%
using StatsBase
swap_negval((val, rep)) = (rep, -val)
inv_swap_negval((rep, val)) = (-val, rep)
minmode(X) = inv_swap_negval(maximum(swap_negval, countmap(X)))k

# %%
X = rand(1:10000:10^8, 10^6)
X'

# %%
using BenchmarkTools

# %%
@btime maxrep_minval($X)

# %%
@btime findn($X)

# %%
@btime minMode($X)

# %%
@btime minmode($X)

# %%
@btime findn_jling($X)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using Test, LinearMaps, LinearAlgebra

@testset "wrapped maps" begin
    A = rand(10, 20)
    B = rand(ComplexF64, 10, 20)
    SA = A'A + I
    SB = B'B + I
    L = @inferred LinearMap{Float64}(A)
    @test summary(L) == "10×20 LinearMaps.WrappedMap{Float64}"
    @test occursin("10×20 LinearMaps.WrappedMap{Float64}", sprint((t, s) -> show(t, "text/plain", s), L))
    MA = @inferred LinearMap(SA)
    MB = @inferred LinearMap(SB)
    @test eltype(Matrix{Complex{Float32}}(LinearMap(A))) <: Complex
    @test size(L) == size(A)
    @test @inferred !issymmetric(L)
    @test @inferred issymmetric(MA)
    @test @inferred !issymmetric(MB)
    @test @inferred isposdef(MA)
    @test @inferred isposdef(MB)

    A = 2 * rand(ComplexF64, (20, 10)) .- 1
    v = rand(ComplexF64, 10)
    w = rand(ComplexF64, 20)
    u = rand(ComplexF64, 20, 1)
    V = rand(ComplexF64, 10, 3)
    W = rand(ComplexF64, 20, 3)
    Av = A * v
    AV = A * V
    M = @inferred LinearMap(A)
    N = @inferred LinearMap(M)

    @test @inferred M' * w == A' * w
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W
    @test @inferred transpose(M) * w == transpose(A) * w
    @test @inferred transpose(LinearMap(transpose(M))) * v == A * v
    @test @inferred LinearMap(M')' * v == A * v
    @test @inferred(transpose(transpose(M))) === M
    @test @inferred(adjoint(adjoint(M))) === M
    Mherm = @inferred LinearMap(A'A)
    @test @inferred ishermitian(Mherm)
    @test @inferred !issymmetric(Mherm)
    @test @inferred !issymmetric(transpose(Mherm))
    @test @inferred ishermitian(transpose(Mherm))
    @test @inferred ishermitian(Mherm')
    @test @inferred isposdef(Mherm)
    @test @inferred isposdef(transpose(Mherm))
    @test @inferred isposdef(adjoint(Mherm))
    @test @inferred !(transpose(M) == adjoint(M))
    @test @inferred !(adjoint(M) == transpose(M))
    @test @inferred transpose(M') * v ≈ transpose(A') * v
    @test @inferred transpose(LinearMap(M')) * v ≈ transpose(A') * v
    @test @inferred LinearMap(transpose(M))' * v ≈ transpose(A)' * v
    @test @inferred transpose(LinearMap(transpose(M))) * v ≈ Av
    @test @inferred adjoint(LinearMap(adjoint(M))) * v ≈ Av

    for w in (vec(u), u)
        @test @inferred mul!(copy(w), transpose(LinearMap(M')), v) ≈ transpose(A') * v
        @test @inferred mul!(copy(w), LinearMap(transpose(M))', v) ≈ transpose(A)' * v
        @test @inferred mul!(copy(w), transpose(LinearMap(transpose(M))), v) ≈ Av
        @test @inferred mul!(copy(w), adjoint(LinearMap(adjoint(M))), v) ≈ Av
    end
    @test @inferred mul!(copy(W), transpose(LinearMap(M')), V) ≈ transpose(A') * V
    @test @inferred mul!(copy(W), LinearMap(transpose(M))', V) ≈ transpose(A)' * V
    @test @inferred mul!(copy(W), transpose(LinearMap(transpose(M))), V) ≈ AV
    @test @inferred mul!(copy(W), adjoint(LinearMap(adjoint(M))), V) ≈ AV
    @test @inferred mul!(copy(V), transpose(M), W) ≈ transpose(A) * W
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W

    B = @inferred LinearMap(Symmetric(rand(10, 10)))
    @test transpose(B) == B
    @test B == transpose(B)

    B = @inferred LinearMap(Hermitian(rand(ComplexF64, 10, 10)))
    @test adjoint(B) == B
    @test B == B'

    N = 10
    Id = LinearMap(identity, N; issymmetric=true)
    A = rand(N, N)
    B = similar(A)
    @test mul!(copy(B), Id, A) == A
    @test mul!(B, A, Id) == B == A
    @test mul!(copy(B), A, LinearMap(Matrix(Id))) == A
    @test mul!(B, Id, A, true, true) ≈ 2A
    @test mul!(B, A, Id, true, true) == B == 3A
    @test mul!(copy(B), A, LinearMap(Matrix(Id)), true, true) == 4A
end


text:
using Test

argA = 0
argB = 1
function sigmoid(x::Float64)
    result = (argB - argA) / (1 + exp(-x)) + argA
    return result
end

@test sigmoid(0.0) ≈ 0.5
@test sigmoid(1.0) ≈ 0.7310585786300049
@test sigmoid(2.0) ≈ 0.8807970779778823
@test sigmoid(-1.0) ≈ 0.2689414213699951
@test sigmoid(-2.0) ≈ 0.11920292202211755


text:
type ClockModel{T} <: ClassicalSpinModel{T}
	spininfo::SpinInfo{T}
	factorW::Array{T,2}

	function ClockModel{T}(spininfo::SpinInfo{T})
		this = new{T}()
		this.spininfo = checkSpininfoForClock(spininfo)
		initialize!(this)
		return this
	end
end

ClockModel{T}(spininfo::SpinInfo{T}) = ClockModel{T}(spininfo)
ClockModel() = testClockModel()

#---
# functions for constructor
function checkSpininfoForClock{T}(spininfo::SpinInfo{T})
	if spininfo.modelname != "clock"
		println("warning:spininfo.modelname does not match to the type.
		  the modelname is changed to 'clock'")
		spininfo.modelname = "clock"
	end
	if isSymmetricFactorization == "sym"
		println("only 'asym'metric factorization for clock model.
		  change factorization to 'asym'.")
		spininfo.factorization = "asym"
	end
	return spininfo
end

function initialize!{T}(clockmodel::ClockModel{T})
	factorW = initializeFactorWasym(clockmodel::ClockModel)
	clockmodel.factorW = factorW
end

function getHamiltonian{T}(clockmodel::ClockModel{T})
	numberOfState = getStates(clockmodel)
	externalfield = getExternalfield(clockmodel)
	sigmah = 1
	Hamiltonian = zeros(T,numberOfState, numberOfState)
	for i = 1:numberOfState, j = 1:numberOfState
		Hamiltonian[i,j] += sigmaTerm(i,j,numberOfState,T) + externalfield * (sigmaTerm(i,sigmah,numberOfState,T) + sigmaTerm(j,sigmah,numberOfState,T))/4.0
	end
	return -Hamiltonian
end
function getAngle(n::Int, numberOfState::Int, T::DataType = Float64)
	return (2.0 * pi * (convert(T,(n-1))))/convert(T,numberOfState)
end
function sigmaTerm(i::Int, j::Int, numberOfState::Int, T::DataType = Float64)
	return cos(getAngle(i,numberOfState,T) - getAngle(j,numberOfState,T))
end

function getMeasureOperator{T}(spinmodel::ClockModel{T})
	NOS = getStates(spinmodel)
	measurementOperator = Array{T}(NOS)
	for i = 1:NOS
		measurementOperator[i] = sigmaTerm(i,1,NOS,T)
	end
	return measurementOperator
end

#---
# test to make a example SpinModel
function testClockModel(numberOfState::Int = 2)
	modelcode = "classical_clock_" * dec(numberOfState)
	factorization = "sym"
	externalfield = 1.0e-13
	temperature = 1.5
	spininfo = SpinInfo(modelcode, factorization, externalfield, temperature)
	testModel = buildSpinSystem(spininfo)
	return testModel
end


text:
using Echidna
using Test

@testset "Echidna" begin
    @testset "Test Problems" begin include("test_problems.jl"); end
    @testset "Core" begin include("test_core.jl"); end
    @testset "NSGA-II" begin include("test_nsgaii.jl"); end
end

# include("test_nsgaiii.jl")

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
#=
Implement the functions using the Unitful package to check units
=#

using Unitful

const uT = Unitful.K
const uP = Unitful.Pa
const uV = Unitful.m^3/Unitful.kg
const umV = Unitful.m^3/Unitful.mol
const uD = Unitful.kg/Unitful.m^3
const umD = Unitful.mol/Unitful.m^3
const perc = Unitful.cm/Unitful.m
const uH = Unitful.J/Unitful.kg
const umH = Unitful.J/Unitful.mol
const uS = Unitful.J/Unitful.kg/Unitful.K
const umS = Unitful.J/Unitful.mol/Unitful.K

dimless{T,U} = Quantity{T, Unitful.Dimensions{()}, U}
ndim = unit(dimless)

val(u, x::Quantity) = uconvert(u, x).val
val(u::Unitful.FreeUnits{(), Unitful.Dimensions{()}}, x) = uconvert(u, x)

function volume(::Type{DryAir}, Tk::Quantity, P::Quantity, u=uV)
    uconvert(u, volume(DryAir, val(uT, Tk), val(uP, P))*uV)
end

function volumem(::Type{DryAir}, Tk::Quantity, P::Quantity, u=umV)
    uconvert(u, volumem(DryAir, val(uT, Tk), val(uP, P))*umV)
end

function density(::Type{DryAir}, Tk::Quantity, P::Quantity, u=uD)
    uconvert(u, density(DryAir, val(uT, Tk), val(uP, P))*uD)
end

function enthalpy(::Type{DryAir}, Tk::Quantity, P::Quantity, u=uH)
    uconvert(u, enthalpy(DryAir, val(uT, Tk), val(uP, P))*uH)
end

function enthalpym(::Type{DryAir}, Tk::Quantity, P::Quantity, u=umH)
    uconvert(u, enthalpym(DryAir, val(uT, Tk), val(uP, P))*umH)
end

function entropy(::Type{DryAir}, Tk::Quantity, P::Quantity, u=uS)
    uconvert(u, entropy(DryAir, val(uT, Tk), val(uP, P))*uS)
end

function entropym(::Type{DryAir}, Tk::Quantity, P::Quantity, u=umS)
    uconvert(u, entropym(DryAir, val(uT, Tk), val(uP, P))*umS)
end

function compressfactor(::Type{DryAir}, Tk::Quantity, P::Quantity)
    compressfactor(DryAir, val(uT, Tk), val(uP, P))
end

function volume(::Type{Vapor}, Tk::Quantity, u=uV)
    uconvert(u, volume(Vapor, val(uT, Tk))*uV)
end

function volumem(::Type{Vapor}, Tk::Quantity, u=umV)
    uconvert(u, volumem(Vapor, val(uT, Tk))*umV)
end

function density(::Type{Vapor}, Tk::Quantity, u=uD)
    uconvert(u, density(Vapor, val(uT, Tk))*uD)
end

function enthalpy(::Type{Vapor}, Tk::Quantity, u=uH)
    uconvert(u, enthalpy(Vapor, val(uT, Tk))*uH)
end

function enthalpym(::Type{Vapor}, Tk::Quantity, u=umH)
    uconvert(u, enthalpym(Vapor, val(uT, Tk))*umH)
end

function entropy(::Type{Vapor}, Tk::Quantity, u=uS)
    uconvert(u, entropy(Vapor, val(uT, Tk))*uS)
end

function entropym(::Type{Vapor}, Tk::Quantity, u=umS)
    uconvert(u, entropym(Vapor, val(uT, Tk))*umS)
end

function compressfactor(::Type{Vapor}, Tk::Quantity)
    compressfactor(Vapor, val(uT, Tk))
end

function molarfrac(Tk::Quantity, ::Type{MolarFrac}, xv::Number, P::Quantity)
    xv
end
function molarfrac(Tk::Quantity, ::Type{HumRat}, w::Number, P::Quantity)
   w / (Mv/Ma + w)
end
function molarfrac(Tk::Quantity, ::Type{SpecHum}, q::Number, P::Quantity)
    q * Ma / (Mv + q*(Ma - Mv))
end

function molarfrac(Tk::Quantity, ::Type{RelHum}, rel, P::Quantity)
    molarfrac(val(uT, Tk), RelHum, rel, val(uP, P))
end

function molarfrac(Tk::Quantity, ::Type{DewPoint}, D::Quantity, P::Quantity)
    molarfrac(Tk, DewPoint, uconvert(uT, D).val, val(uP, P))
end

function molarfrac(Tk::Quantity, ::Type{WetBulb}, B::Quantity, P::Quantity)
    molarfrac(val(uT, Tk), WetBulb, uconvert(uT, B).val, val(uP, P))
end

function volume(::Type{MoistAir}, Tk::Quantity,
                ::Type{T}, y, P::Quantity, u=uV) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u, volume(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))*uV)
end

function volumem(::Type{MoistAir}, Tk::Quantity,
                 ::Type{T}, y, P::Quantity, u=umV) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u, volumem(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))*umV)
end

function density(::Type{MoistAir}, Tk::Quantity,
                 ::Type{T}, y, P::Quantity, u=uD) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u, density(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))*uD)
end

function enthalpy(::Type{MoistAir}, Tk::Quantity,
                  ::Type{T}, y, P::Quantity, u=uH) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u, enthalpy(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))*uH)
end

function enthalpym(::Type{MoistAir}, Tk::Quantity,
                   ::Type{T}, y, P::Quantity, u=umH) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u, enthalpym(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))*umH)
end

function entropy(::Type{MoistAir}, Tk::Quantity,
                  ::Type{T}, y, P::Quantity, u=uS) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u, entropy(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))*uS)
end

function entropym(::Type{MoistAir}, Tk::Quantity,
                  ::Type{T}, y, P::Quantity, u=umS) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u, entropym(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))*umS)
end

function compressfactor(::Type{MoistAir}, Tk::Quantity,
                        ::Type{T}, y, P::Quantity) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    compressfactor(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))
end

function dewpoint(::Type{MoistAir}, Tk::Quantity,
                  ::Type{T}, y, P::Quantity, u=u"°C") where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(Unitful.°C, dewpoint(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))*uT)
end

function relhum(::Type{MoistAir}, Tk::Quantity,
                ::Type{T}, y, P::Quantity, u=ndim) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u, relhum(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P)))
end

function wetbulb(::Type{MoistAir}, Tk::Quantity,
                 ::Type{T}, y, P::Quantity, u=Unitful.°C) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u"°C", wetbulb(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))*uT)
end

function humrat(::Type{MoistAir}, Tk::Quantity,
                ::Type{T}, y, P::Quantity, u=ndim) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    uconvert(u, humrat(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P)))
end

function molarfrac(::Type{MoistAir}, Tk::Quantity,
                   ::Type{T}, y, P::Quantity) where {T<:PsychroProperty}

    molarfrac(Tk, T, y, P)
end

function spechum(::Type{MoistAir}, Tk::Quantity,
                ::Type{T}, y, P::Quantity) where {T<:PsychroProperty}

    xv = molarfrac(Tk, T, y, P)
    spechum(MoistAir, val(uT, Tk), MolarFrac, xv, val(uP, P))
end


text:
## These functions define the possible iterations for the
## sorted containers.
## The prefix SDM is for SortedDict and SortedMultiDict
## The prefix SS is for SortedSet.  The prefix SA
## is for all sorted containers.
## The following two definitions now appear in tokens2.jl

# const SDMContainer = Union{SortedDict, SortedMultiDict}
# const SAContainer = Union{SDMContainer, SortedSet}

extractcontainer(s::SAContainer) = s
getrangeobj(s::SAContainer) = s

## This holds an object describing an exclude-last
## iteration.

abstract type AbstractExcludeLast{ContainerType <: SAContainer} end

struct SDMExcludeLast{ContainerType <: SDMContainer} <:
                              AbstractExcludeLast{ContainerType}
    m::ContainerType
    first::Int
    pastlast::Int
end

keytype(::SDMExcludeLast{T}) where {T <: SAContainer} = keytype(T)
keytype(::Type{SDMExcludeLast{T}}) where {T <: SAContainer} = keytype(T)
valtype(::SDMExcludeLast{T}) where {T <: SAContainer} = valtype(T)
valtype(::Type{SDMExcludeLast{T}}) where {T <: SAContainer} = valtype(T)
eltype(::SDMExcludeLast{T}) where {T <: SAContainer} = eltype(T)
eltype(::Type{SDMExcludeLast{T}}) where {T <: SAContainer} = eltype(T)

struct SSExcludeLast{ContainerType <: SortedSet} <:
                              AbstractExcludeLast{ContainerType}
    m::ContainerType
    first::Int
    pastlast::Int
end

eltype(::SSExcludeLast{T}) where {T <: SortedSet} = eltype(T)
eltype(::Type{SSExcludeLast{T}}) where {T <: SortedSet} = eltype(T)

extractcontainer(s::AbstractExcludeLast) = s.m
getrangeobj(s::AbstractExcludeLast) = s

## This holds an object describing an include-last
## iteration.

abstract type AbstractIncludeLast{ContainerType <: SAContainer} end

struct SDMIncludeLast{ContainerType <: SDMContainer} <:
                               AbstractIncludeLast{ContainerType}
    m::ContainerType
    first::Int
    last::Int
end

keytype(::SDMIncludeLast{T}) where {T <: SAContainer} = keytype(T)
keytype(::Type{SDMIncludeLast{T}}) where {T <: SAContainer} = keytype(T)
valtype(::SDMIncludeLast{T}) where {T <: SAContainer} = valtype(T)
valtype(::Type{SDMIncludeLast{T}}) where {T <: SAContainer} = valtype(T)
eltype(::SDMIncludeLast{T}) where {T <: SAContainer} = eltype(T)
eltype(::Type{SDMIncludeLast{T}}) where {T <: SAContainer} = eltype(T)

struct SSIncludeLast{ContainerType <: SortedSet} <:
                               AbstractIncludeLast{ContainerType}
    m::ContainerType
    first::Int
    last::Int
end

eltype(::SSIncludeLast{T}) where {T <: SortedSet} = eltype(T)
eltype(::Type{SSIncludeLast{T}}) where {T <: SortedSet} = eltype(T)

extractcontainer(s::AbstractIncludeLast) = s.m
getrangeobj(s::AbstractIncludeLast) = s

IteratorSize(::Type{T} where {T <: SAContainer}) = HasLength()
IteratorSize(::Type{T} where {T <: AbstractExcludeLast}) = SizeUnknown()
IteratorSize(::Type{T} where {T <: AbstractIncludeLast}) = SizeUnknown()

## The basic iterations are either over the whole sorted container, an
## exclude-last object or include-last object.

const SDMIterableTypesBase = Union{SDMContainer,
                                   SDMExcludeLast,
                                   SDMIncludeLast}

const SSIterableTypesBase = Union{SortedSet,
                                  SSExcludeLast,
                                  SSIncludeLast}

const SAIterableTypesBase = Union{SAContainer,
                                  AbstractExcludeLast,
                                  AbstractIncludeLast}

## The compound iterations are obtained by applying keys(..) or values(..)
## to the basic iterations of the SDM.. type.
## Furthermore, semitokens(..) can be applied
## to either a basic iteration or a keys/values iteration.

struct SDMKeyIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(::Type{SDMKeyIteration{T}}) where {T} = keytype(T)
eltype(s::SDMKeyIteration) = keytype(extractcontainer(s.base))
length(s::SDMKeyIteration{T} where T <: SDMContainer) = length(extractcontainer(s.base))

struct SDMValIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(::Type{SDMValIteration{T}}) where {T} = valtype(T)
eltype(s::SDMValIteration) = valtype(extractcontainer(s.base))
length(s::SDMValIteration{T} where T <: SDMContainer) = length(extractcontainer(s.base))

struct SDMSemiTokenIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(::Type{SDMSemiTokenIteration{T}}) where {T} =
    Tuple{IntSemiToken, keytype(T), valtype(T)}
eltype(s::SDMSemiTokenIteration) = Tuple{IntSemiToken,
                                         keytype(extractcontainer(s.base)),
                                         valtype(extractcontainer(s.base))}
length(s::SDMSemiTokenIteration{T} where T <: SDMContainer) = length(s.base)

struct SSSemiTokenIteration{T <: SSIterableTypesBase}
    base::T
end

eltype(::Type{SSSemiTokenIteration{T}}) where {T} =
    Tuple{IntSemiToken, eltype(T)}
eltype(s::SSSemiTokenIteration) = Tuple{IntSemiToken,
                                        eltype(extractcontainer(s.base))}
length(s::SSSemiTokenIteration{T} where T <: SortedSet) = length(s.base)

struct SDMSemiTokenKeyIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(::Type{SDMSemiTokenKeyIteration{T}}) where {T} =
    Tuple{IntSemiToken,
          keytype(T)}
eltype(s::SDMSemiTokenKeyIteration) = Tuple{IntSemiToken,
                                            keytype(extractcontainer(s.base))}
length(s::SDMSemiTokenKeyIteration{T} where T <: SDMContainer) = length(s.base)

struct SAOnlySemiTokensIteration{T <: SAIterableTypesBase}
    base::T
end

eltype(::Type{SAOnlySemiTokensIteration{T}} where {T}) = IntSemiToken
eltype(::SAOnlySemiTokensIteration) = IntSemiToken
length(s::SAOnlySemiTokensIteration{T} where T <: SAContainer) = length(s.base)

struct SDMSemiTokenValIteration{T <: SDMIterableTypesBase}
    base::T
end

eltype(::Type{SDMSemiTokenValIteration{T}}) where {T} =
    Tuple{IntSemiToken, valtype(T)}
eltype(s::SDMSemiTokenValIteration) = Tuple{IntSemiToken,
                                            valtype(extractcontainer(s.base))}
length(s::SDMSemiTokenValIteration{T} where T <: SDMContainer) = length(s.base)

const SACompoundIterable = Union{SDMKeyIteration,
                                 SDMValIteration,
                                 SDMSemiTokenIteration,
                                 SSSemiTokenIteration,
                                 SDMSemiTokenKeyIteration,
                                 SDMSemiTokenValIteration,
                                 SAOnlySemiTokensIteration}

extractcontainer(s::SACompoundIterable) = extractcontainer(s.base)
getrangeobj(s::SACompoundIterable) = getrangeobj(s.base)

const SAIterable = Union{SAIterableTypesBase, SACompoundIterable}

IteratorEltype(::Type{T} where {T <: SAIterable}) = HasEltype()
IteratorSize(::Type{SDMKeyIteration{T}}) where {T} = IteratorSize(T)
IteratorSize(::Type{SDMValIteration{T}}) where {T} = IteratorSize(T)
IteratorSize(::Type{SDMSemiTokenIteration{T}}) where {T} = IteratorSize(T)
IteratorSize(::Type{SSSemiTokenIteration{T}}) where {T} = IteratorSize(T)
IteratorSize(::Type{SDMSemiTokenKeyIteration{T}}) where {T} = IteratorSize(T)
IteratorSize(::Type{SDMSemiTokenValIteration{T}}) where {T} = IteratorSize(T)
IteratorSize(::Type{SAOnlySemiTokensIteration{T}}) where {T} = IteratorSize(T)

## All the loops maintain a state which is an object of the
## following type.

struct SAIterationState
    next::Int
    final::Int
end

exclusive(m::SDMContainer, ii::Tuple{IntSemiToken,IntSemiToken}) =
    SDMExcludeLast(m, ii[1].address, ii[2].address)
exclusive(m::SortedSet, ii::Tuple{IntSemiToken,IntSemiToken}) =
    SSExcludeLast(m, ii[1].address, ii[2].address)
exclusive(m::SAContainer, i1::IntSemiToken, i2::IntSemiToken) =
    exclusive(m, (i1, i2))

inclusive(m::SDMContainer, ii::Tuple{IntSemiToken,IntSemiToken}) =
    SDMIncludeLast(m, ii[1].address, ii[2].address)
inclusive(m::SortedSet, ii::Tuple{IntSemiToken,IntSemiToken}) =
    SSIncludeLast(m, ii[1].address, ii[2].address)
inclusive(m::SAContainer, i1::IntSemiToken, i2::IntSemiToken) =
    inclusive(m, (i1, i2))

# Next definition needed to break ambiguity with keys(AbstractDict) from Dict.jl

keys(ba::SortedDict) = SDMKeyIteration(ba)
keys(ba::SDMIterableTypesBase) = SDMKeyIteration(ba)

in(k, keyit::SDMKeyIteration{SortedDict{K,D,Ord}} where {K,D,Ord}) =
    haskey(extractcontainer(keyit.base), k)

in(k, keyit::SDMKeyIteration{SortedMultiDict{K,D,Ord}} where {K,D,Ord}) =
    haskey(extractcontainer(keyit.base), k)

# Next definition needed to break ambiguity with values(AbstractDict) from Dict.jl
values(ba::SortedDict) = SDMValIteration(ba)
values(ba::SDMIterableTypesBase) = SDMValIteration(ba)
semitokens(ba::SDMIterableTypesBase) = SDMSemiTokenIteration(ba)
semitokens(ba::SSIterableTypesBase) = SSSemiTokenIteration(ba)
semitokens(ki::SDMKeyIteration) = SDMSemiTokenKeyIteration(ki.base)
semitokens(vi::SDMValIteration) = SDMSemiTokenValIteration(vi.base)
onlysemitokens(ba::SAIterableTypesBase) = SAOnlySemiTokensIteration(ba)

function nexthelper(c::SAContainer, state::SAIterationState)
    sn = state.next
    (sn < 3 || !(sn in c.bt.useddatacells)) && throw(BoundsError())
    SAIterationState(nextloc0(c.bt, sn), state.final)
end

getitem(::SDMIterableTypesBase, dt, sn) = dt.k => dt.d
getitem(::SSIterableTypesBase, dt, sn) = dt.k
getitem(::SDMKeyIteration, dt, sn) = dt.k
getitem(::SDMValIteration, dt, sn) = dt.d
getitem(::SDMSemiTokenIteration, dt, sn) = (IntSemiToken(sn), dt.k, dt.d)
getitem(::SSSemiTokenIteration, dt, sn) = (IntSemiToken(sn), dt.k)
getitem(::SDMSemiTokenKeyIteration, dt, sn) = (IntSemiToken(sn), dt.k)
getitem(::SDMSemiTokenValIteration, dt, sn) = (IntSemiToken(sn), dt.d)
getitem(::SAOnlySemiTokensIteration, dt, sn) = IntSemiToken(sn)

function get_init_state(e::AbstractExcludeLast)
    (!(e.first in e.m.bt.useddatacells) || e.first == 1 ||
     !(e.pastlast in e.m.bt.useddatacells)) &&
     throw(BoundsError())
    if compareInd(e.m.bt, e.first, e.pastlast) < 0
        return SAIterationState(e.first, e.pastlast)
    else
        return SAIterationState(2, 2)
    end
end

function get_init_state(e::AbstractIncludeLast)
    (!(e.first in e.m.bt.useddatacells) || e.first == 1 ||
     !(e.last in e.m.bt.useddatacells) || e.last == 2) &&
     throw(BoundsError())
    if compareInd(e.m.bt, e.first, e.last) <= 0
        return SAIterationState(e.first, nextloc0(e.m.bt, e.last))
    else
        return SAIterationState(2, 2)
    end
end

get_init_state(m::SAContainer) = SAIterationState(beginloc(m.bt), 2)

function iterate(s::SAIterable, state = get_init_state(getrangeobj(s)))
    if state.next == state.final
        return nothing
    else
        c = extractcontainer(s)
        dt = isa(s, SAOnlySemiTokensIteration) ? nothing : c.bt.data[state.next]
        return (getitem(s, dt, state.next),
                nexthelper(c, state))
    end
end

eachindex(sd::SortedDict) = keys(sd)
eachindex(sdm::SortedMultiDict) = onlysemitokens(sdm)
eachindex(ss::SortedSet) = onlysemitokens(ss)
eachindex(sd::SDMExcludeLast{SortedDict{K,D,Ord}} where {K,D,Ord <: Ordering}) = keys(sd)
eachindex(smd::SDMExcludeLast{SortedMultiDict{K,D,Ord}} where {K,D,Ord <: Ordering}) =
     onlysemitokens(smd)
eachindex(ss::SSExcludeLast) = onlysemitokens(ss)
eachindex(sd::SDMIncludeLast{SortedDict{K,D,Ord}} where {K,D,Ord <: Ordering}) = keys(sd)
eachindex(smd::SDMIncludeLast{SortedMultiDict{K,D,Ord}} where {K,D,Ord <: Ordering}) =
     onlysemitokens(smd)
eachindex(ss::SSIncludeLast) = onlysemitokens(ss)

empty!(m::SAContainer) =  empty!(m.bt)
length(m::SAContainer) = length(m.bt.data) - length(m.bt.freedatainds) - 2
isempty(m::SAContainer) = length(m) == 0


text:
using Pkg; Pkg.activate(".")
using BKTrees
using Random
using StringDistances

lev(x::S,y::S) where S<:AbstractString = evaluate(Levenshtein(), x, y)
dictionary = [randstring(10) for _ in 1:10_000]
bkt = BKTree(lev, sort(dictionary))

target = randstring(10)

found = find(bkt, "bbb", 10, k=3)
@time found=find(bkt, target, 10, k=3)
@show target, found