text:
mutable struct Data{
    S <: AbstractString,
    A <: AbstractArray{T, 3} where T,
    B <: AbstractMatrix, C <: AbstractMatrix,
    D <: AbstractMatrix, E <: AbstractMatrix,
    F<: AbstractMatrix, G <: AbstractMatrix,
    H <: AbstractMatrix, I <: AbstractMatrix,
    J <: AbstractMatrix, K <: AbstractMatrix
    }
    特征名::Dict{S, Int}
    特征::A
    涨幅::B
    时间戳::C
    代码::D
    最新价::E
    买1价::F
    卖1价::G
    手续费率::H
    涨停::I
    跌停::J
    交易池::K
end

Base.:(==)(x::Data, y::Data) = all(s -> isequal(getfield(x, s), getfield(y, s)), fieldnames(Data))

afieldnames(t) = Symbol[fieldname(t, n) for n in 1:fieldcount(t) if fieldtype(t, n) <: AbstractArray]

featnames(data) = first.(sort(collect(data.特征名), by = last))

mapdata(f, data) = Data(data.特征名, [f(getfield(data, s)) for s in afieldnames(Data)]...)

for f in [:getindex, :view]
    @eval Base.$f(data::Data, is...) = mapdata(data) do x
        ndims(x) == 2 ? $f(x, is...) : $f(x, :, is...)
    end

    @eval function Base.$f(data::Data, ts::StringRange)
        ti = str2unix(ts.start)
        tf = str2unix(ts.stop)
        mask = ti .<= data.时间戳 .<= tf
        ns = findall(dropdims(any, mask, dims = 2))
        ts = findall(dropdims(any, mask, dims = 1))
        $f(data, arr2rng(ns), arr2rng(ts))
    end

    @eval function Base.$f(data::Data, code::String)
        code = eltype(data.代码)(code)
        mask = data.代码 .== code
        ns = findall(vec(any(mask, dims = 2)))
        ts = findall(vec(any(mask, dims = 1)))
        $f(data, arr2rng(ns), arr2rng(ts))
    end
end

Base.maybeview(data::Data, is...) = view(data, is...)

Base.reshape(data::Data, is...) = mapdata(data) do x
    ndims(x) == 2 ? reshape(x, is...) : reshape(x, size(x, 1), is...)
end

Base.vec(data::Data) = reshape(data, 1, length(data))

nfeats(data) = size(data.特征, 1)

ncodes(data) = size(data.涨幅, 1)

nticks(data) = size(data.涨幅, 2)

Base.size(data::Data) = size(data.涨幅)

Base.size(data::Data, dim) = size(data.涨幅, dim)

Base.length(data::Data) = prod(size(data))

Base.lastindex(data::Data, n) = size(data, n)

function ndays(data)
    Δts = Float32[]
    for n in 1:size(data, 1)
        ts = filter(!isna, data.时间戳[n, :])
        !isempty(ts) && return sortednunique(unix2date, ts)
    end
    return 0
end

nticksperday(data) = nticks(data) ÷ ndays(data)

const _edgemap = Dict{String, Vector{Float32}}()

const _sourcemap = Dict{UInt, String}()

sourceof(x) = get(_sourcemap, objectid(x), nothing)

function _loaddata(src; fload = nothing, ti = nothing, tf = nothing, ka...)
    if endswith(src, ".bson")
        data = isnothing(fload) ? bsload(src, Data; ka...) : fload(src; ka...)
    else
        data = isnothing(fload) ? h5load(src, Data; ka...) : fload(src; ka...)
        if "bin_edges" ∈ h5open(names, src)
            bin_edges = h5read(src, "bin_edges")
            for (c, f) in data.特征名
                _edgemap[c] = bin_edges[:, f]
            end
        end
    end
    if isnothing(ti) && isnothing(tf)
        _sourcemap[objectid(data)] = src
        return data
    else
        ti = something(ti, "20000101")
        tf = something(tf, "20501231")
        data′ = @view data[:, ti:tf]
        _sourcemap[objectid(data′)] = src
        return data′
    end
end

function loaddata(srcs::AbstractArray, a...; dims = -1, ka...)
    datas = @showprogress 10 "loaddata..." map(srcs) do src
        _loaddata(src, a...; ka...)
    end
    length(srcs) == 1 && return datas[1]
    concat(filter(!isempty, datas); dims)
end

function loaddata(pattern, a...; ka...)
    if isfile(pattern)
        srcs = [pattern]
    elseif occursin("*", pattern)
        if startswith(pattern, "/")
            srcs = glob(pattern[2:end], "/")
        else
            srcs = glob(pattern)
        end
    else
        srcs = []
    end
    if isempty(srcs) && !startswith(pattern, "/")
        job = get(ENV, "JOB", expanduser("~/job"))
        srcs = glob(pattern, joinpath(job, "data"))
    end
    isempty(srcs) && error(pattern * " not found")
    loaddata(srcs, a...; ka...)
end

function savedata(dst, data)
    if endswith(dst, ".bson")
        bssave(dst, data[:, :])
    elseif endswith(dst, ".h5")
        h5save(dst, data)
    end
    return dst
end

function reloaddata(data)
    dst = @sprintf("data-%s.h5", randstring())
    savedata(dst, data)
    loaddata(dst, mode = "r+")
end

function initdata(dst, eltype_, (F, N, T), feature = nothing)
    isfile(dst) && rm(dst)
    feature = something(feature, string.(1:F))
    h5open(dst, "w") do fid
        g_create(fid, "nonarray")
        @showprogress 10 "initdata..." for s in afieldnames(Data)
            if s == :时间戳
                d_zeros(fid, string(s), Float64, (N, T))
            elseif s == :代码
                fid["代码"] = [MLString{8}(string(n)) for n in 1:N, t in 1:T]
            elseif s == :涨停 || s == :跌停 || s == :交易池
                d_zeros(fid, string(s), UInt8, (N, T))
            elseif s == :特征
                d_zeros(fid, string(s), eltype_, (F, N, T))
            else
                d_zeros(fid, string(s), Float32, (N, T))
            end
        end
        fid["交易池"][:, :] = 1
        if !isempty(feature)
            特征名 = Dict(reverse(p) for p in enumerate(feature))
            write_nonarray(fid, "特征名", 特征名)
        end
    end
    return dst
end

function Base.show(io::IO, data::Data)
    isempty(data) && return
    compact = get(io, :compact, false)
    @printf(io, "特征数: %d\t", nfeats(data))
    @printf(io, "品种数: %d\t", ncodes(data))
    @printf(io, "TK数: %d\n", nticks(data))
    @printf(io, "交易池比例: %.2g\t", mean(data.交易池))
    @printf(io, "涨停比例: %.2g\t", mean(data.涨停))
    @printf(io, "跌停比例: %.2g\n", mean(data.跌停))
    @printf(io, "日期范围: %s/%s\t", Dates.format.(datespan(data), dateformat"yymmdd")...)
    @printf(io, "价格范围: %.3g/%.3g\n", extrema(filter(!isna, data.最新价))...)
    @printf(io, "涨幅范围: %.2g/%.2g\n", extrema(filter(!isna, data.涨幅))...)
    compact && return
    header, stats = String[], Array{String}[]
    for f in ("涨幅", "手续费率")
        x = vec(getfield(data, Symbol(f)))
        s = StatsBase.summarystats(x)
        push!(header, f)
        push!(stats, split(string(s), '\n')[2:end-1])
    end
    for (f, i) in data.特征名
        x = filter(z -> !isnan(z) & !isinf(z), data.特征[i, :, :])
        if any(!iszero, x)
            s = StatsBase.summarystats(x)
        else
            s = StatsBase.SummaryStats(zeros(6)..., length(x), 0)
        end
        push!(header, string(f))
        push!(stats, split(string(s), '\n')[2:end-1])
    end
    for (header′, stats′) in zip(Iterators.partition(header, 4), Iterators.partition(stats, 4))
        print_header_stats(io, header′, stats′)
    end
end

Base.show(io::IO, ::MIME"text/plain", data::Data) =
    show(IOContext(io, :compact => true), data)

function print_header_stats(io, header, stats)
    print(io, '\n', '\n')
    for h in header
        print(io, h, ['\t' for i in 1:4]...)
    end
    for r in 1:length(stats[1])
        print(io, '\n')
        for c in 1:length(stats)
            print(io, stats[c][min(r, end)], '\t')
        end
    end
end

function period(data)
    Δts = Float32[]
    for n in 1:size(data, 1)
        ts = filter(!isna, data.时间戳[n, :])
        length(ts) < 2 && continue
        Δt = median(diff(ts))
        2 * length(ts) > size(data, 2) && return Δt
        push!(Δts, Δt)
    end
    isempty(Δts) ? 0f0 : median(Δts)
end

downsample(data::Data, freq::String; ka...) =
    downsample(data, round(Int, min(nticks(data), parsefreq(freq) / period(data))); ka...)

function downsample(data::Data, freq::Int; phase = 1, ka...)
    freq <= 1 && return data[:, :]
    downsample(data, phase:freq:nticks(data); ka...)
end

function downsample(data::Data, ts::AbstractArray{Int}; average = false)
    data′ = data[:, ts]
    F, N, T = size(data.特征)
    fill!(data′.涨幅, 0)
    涨幅 = Array(data.涨幅)
    for t′ in 1:(length(ts) - 1)
        t⁻ = t′ == 1 ? 0 : ts[t′ - 1]
        @inbounds for t in (t⁻ + 1):ts[t′], n in 1:N
            data′.涨幅[n, t′] += 涨幅[n, t]
        end
    end
    if length(ts) == 1
        copyto!(data′.涨幅, sum(涨幅, dims = 2))
    end
    !average && return data′
    fill!(data′.特征, 0)
    @showprogress 10 for t′ in 1:(length(ts) - 1)
        t⁻ = t′ == 1 ? 0 : ts[t′ - 1]
        Δt = ts[t′] - t⁻
        @inbounds for t in (t⁻ + 1):ts[t′], n in 1:N, f in 1:F
            data′.特征[f, n, t′] += data.特征[f, n, t] / Δt
        end
    end
    return data′
end

function getfeat(data::Data, c; denorm = true, errors = "raise")
    if c in keys(data.特征名)
        f = data.特征名[c]
        x = @view data.特征[f, :, :]
        if haskey(_edgemap, c) && denorm
            undiscretize(x, _edgemap[c])
        end
        return x
    elseif Symbol(c) in fieldnames(Data)
        getfield(data, Symbol(c))
    elseif errors == "ignore"
        zeros(Float32, size(data))
    else
        error("$c not in data or data.特征名")
    end
end

function getfeat(data::Data, r::Regex; ka...)
    cs = getfeats(data, r; ka...)
    !isempty(cs) && return cs[1]
    getfeat(data, string(r); ka...)
end

getfeats(data::Data, cs; ka...) = [getfeat(data, c; ka...) for c in cs]

getfeats(data::Data, r::Regex; ka...) = getfeats(data, filter(c -> occursin(r, c), featnames(data)); ka...)

function dropfeats(data, cs::AbstractArray{Int})
    isempty(cs) && return data
    dict = to_dict(data)
    fs = setdiff(1:nfeats(data), cs)
    dict[:特征名] = idxmap(featnames(data)[fs])
    dict[:特征] = view(data.特征, fs, :, :)
    to_struct(Data, dict)
end

dropfeats(data, cs::Array{<:AbstractString}) = dropfeats(data, [get(data.特征名, c, -1) for c in cs])
dropfeats(data, r::Regex) = dropfeats(data, filter(c -> occursin(r, c), featnames(data)))
keepfeats(data, cs::AbstractArray{Int}) = dropfeats(data, setdiff(1:nfeats(data), cs))
keepfeats(data, cs::Array{<:AbstractString}) = dropfeats(data, setdiff(featnames(data), cs))
keepfeats(data, r::Regex) = keepfeats(data, filter(c -> occursin(r, c), featnames(data)))

getcats(data) = filter(!isempty, unique(String.(first.(split.(featnames(data), ':')))))
keepcats(data, cats) = keepfeats(data, filter(c -> any(cat -> startswith(c, cat) || !occursin(':', c), splat(cats)), featnames(data)))
dropcats(data, cats) = keepcats(data, setdiff(getcats(data), splat(cats)))
dropcats(data, r::Regex) = dropcats(data, filter(c -> occursin(r, c), getcats(data)))
keepcats(data, r::Regex) = keepcats(data, filter(c -> occursin(r, c), getcats(data)))

macro uncol(ex)
    cs, rhs = ex.args
    d, cs = gensym(), isa(cs, Symbol) ? [cs] : cs.args
    kd = [:($c = $getfeat($d, $(string(c)))) for c in cs]
    Expr(:block, :($d = $rhs), kd..., d) |> esc
end

Base.split(data::Data, t::String) = @views data["20010101":t], data[t:"30000101"]

datespan(data::Data) = (firstdate(data), lastdate(data))

firstdate(data::Data) = unix2date(minimum(t -> ifelse(isna(t), Inf, t), data.时间戳))

lastdate(data::Data) = unix2date(maximum(t -> ifelse(isna(t), -Inf, t), data.时间戳))

pct_change(data::Data, h::String) = pct_change(data, ceil(Int, min(nticks(data), parsefreq(h) / period(data))))

pct_change(data::Data, h::Int) = pct_change(Array(data.涨幅), h)

function pct_change(v, h)
    h = ceil(Int, h)
    N, T = size(v)
    chg = zeros(Float32, N, T)
    r = zeros(Float32, N)
    for t in T:-1:max(2, T - h + 1), n in 1:N
        r[n] += v[n, t]
        chg[n, t - 1] = r[n]
    end
    for t in (T - h):-1:2
        for n in 1:N
            r[n] += v[n, t] - v[n, t + h]
            chg[n, t - 1] = r[n]
        end
    end
    return chg
end

function setcomm(data, c)
    dict = to_dict(data)
    dict[:手续费率] = fill(Float32(c), size(data))
    to_struct(Data, dict)
end

function setpool(data, pool)
    ns = findall(vec(any(isone, pool, dims = 2)))
    data′, pool′ = data[ns, :], pool[ns, :]
    data′.交易池 = isone.(pool′)
    return data′
end

setpool(data, c::Union{String, Regex}) = setpool(data, getfeat(data, c))

function Base.repeat(data::Data, n)
    n == 1 && return data
    fvs = []
    for s in fieldnames(Data)
        x = getfield(data, s)
        if isa(x, AbstractMatrix)
            bsizes = (fill(size(x, 1), n), [size(x, 2)])
            x = _BlockArray(repeat([x], n, 1), bsizes...)
        end
        push!(fvs, x)
    end
    return Data(fvs...)
end

epochsof(data::Data) = sort(unique(data.时间戳))

datetimesof(data::Data) = map(unix2datetime, epochsof(data))

datesof(data::Data) = map(unix2date, epochsof(data))

codesof(data::Data) = sort(unique(data.代码))

normal_mask(data::Data) = @. iszero(data.涨停) | iszero(data.跌停)

function isaligned(datas)
    isempty(datas) && return true
    codes = datas[1].代码[:, 1]
    for data in datas[2:end]
        if data.代码[:, 1] != codes
            return false
        end
    end
    return true
end

function concat(datas; dims = -1)
    fvs = []
    for s in fieldnames(Data)
        xs = getfield.(datas, s)
        x1 = getfield(datas[1], s)
        if isa(x1, AbstractArray)
            d1 = dims > 0 ? dims : (ndims(x1) + 1 + dims)
            bsizes = [d == d1 ? size.(xs, d) : [size(x1, d)] for d in 1:ndims(x1)]
            xsizes = [d == d1 ? (:) : 1 for d in 1:ndims(x1)]
            x = _BlockArray(reshape(xs, xsizes...), bsizes...)
        else
            x = x1
        end
        push!(fvs, x)
    end
    return Data(fvs...)
end

Base.isempty(data::Data) = length(data) == 0 || !any(!isna, data.时间戳)

Base.copy(data::Data) = to_struct(Data, to_dict(data))

pivot(datas::AbstractArray{<:Data}, a...; ka...) = pivot(concat(map(vec, datas)), a...; ka...)

pivot(h5::String, a...; ka...) = pivot(loaddata(h5), a...; ka...)

function pivot(data::Data, dst = "pivot.h5"; meta_only = false)
    epochs = epochsof(data)
    codes = codesof(data)
    F = nfeats(data)
    N, T = size(data)
    N′ = length(codes)
    T′ = length(epochs)
    if meta_only
        F, dict = 1, to_dict(data)
        dict[:特征] = zeros(eltype(data), 1,N, T)
        dict[:特征名] = Dict("dummy" => 1)
        data = to_struct(Data, dict)
    end
    index = fill((0, 0), N, T)
    @unpack 代码, 时间戳 = data
    codemap = Dict(reverse.(enumerate(codes)))
    epochmap = Dict(reverse.(enumerate(epochs)))
    @showprogress 10 "pivot.index" for t in 1:T, n in 1:N
        n′ = codemap[代码[n, t]]
        t′ = epochmap[时间戳[n, t]]
        index[n, t] = (n′, t′)
    end
    data′ = similar(data, (F, N′, T′), dst)
    fill!(data′.涨停, 1)
    fill!(data′.跌停, 1)
    fill!(data′.交易池, 0)
    @inbounds for s in afieldnames(Data)
        src = getfield(data, s)
        dest = getfield(data′, s)
        desc = string("pivot.", s)
        p = Progress(T, desc = desc)
        if s == :特征
            Threads.@threads for t in 1:T
                for n in 1:N
                    n′, t′ = index[n, t]
                    for f in 1:F
                        dest[f, n′, t′] = src[f, n, t]
                    end
                end
                next!(p)
            end
        elseif s == :时间戳
            Threads.@threads for t′ in 1:T′
                for n′ in 1:N′
                    dest[n′, t′] = epochs[t′]
                end
                next!(p)
            end
        elseif s == :代码
            Threads.@threads for t′ in 1:T′
                for n′ in 1:N′
                    dest[n′, t′] = codes[n′]
                end
            end
        else
            Threads.@threads for t in 1:T
                for n in 1:N
                    n′, t′ = index[n, t]
                    dest[n′, t′] = src[n, t]
                end
            end
            s == :最新价 && Threads.@threads for t′ in 2:T′
                for n′ in 1:N′
                    if iszero(dest[n′, t′])
                        dest[n′, t′] = dest[n′, t′ - 1]
                    end
                end
                next!(p)
            end
            s == :最新价 && Threads.@threads for t′ in (T′ - 1):-1:1
                for n′ in 1:N′
                    if iszero(dest[n′, t′])
                        dest[n′, t′] = dest[n′, t′ + 1]
                    end
                end
                next!(p)
            end
        end
    end
    @. data′.跌停 = ifelse(isone(data′.交易池), data′.跌停, 0f0)
    @. data′.涨停 = ifelse(isone(data′.交易池), data′.涨停, 0f0)
    Mmap.sync!(data′)
    return data′, index
end

function pivot(src::AbstractArray{T, N}, index) where {T, N}
    dims′ = ntuple(d -> maximum(I -> I[d], index), N)
    dest = zeros(T, dims′...)
    for n in 1:length(index)
        dest[index[n]...] = src[n]
    end
    return dest
end

Base.show(io::IO, ::MIME"text/plain", ::Array{<:Data}) =
    print(io, "Array of Data")

function rolldata(data, rolltrn, rolltst)
    ti, tf = map(unix2date, extrema(data.时间戳))
    Δtb = Day(parsefreq(rolltrn) ÷ 3600 ÷ 24)
    Δtf = Day(parsefreq(rolltst) ÷ 3600 ÷ 24)
    inds = rollindices(ti, tf, Δtb, Δtf)
    @views ((data[t1], data[t2]) for (t1, t2) in inds)
end

function to_df(data::Data; columns = nothing, meta_only = false, cnvtdate = false)
    df = DataFrame()
    @showprogress 10 "to_df..." for s in afieldnames(Data)
        x = vec(getfield(data, s))
        !isnothing(columns) && string(s) ∉ columns && continue
        if s == :代码
            df[s] = to_category(x)
        elseif s == :时间戳
            df[s] = x
            if cnvtdate
                df[s] = df[s].mul(1e9).astype("datetime64[ns]")
            end
        elseif s == :特征 && !meta_only
            df = pdhcat(df, pd.DataFrame(pymat(data.特征), columns = featnames(data)))
        else
            df[s] = x
        end
    end
    return df
end

function to_data(df, dst; ka...)
    try
        _to_data(df, dst; ka...)
    catch e
        rm(dst, force = true)
        throw(e)
    end
end

function _to_data(df, dst; ncode = 0)
    isempty(df) && return
    F = length(featcols(df))
    N = ncode == 0 ? df["代码"].nunique() : ncode
    T = length(df) ÷ N
    @assert N * T == length(df)
    dst = initdata(dst, Float32, (F, N, T), featcols(df))
    h5open(dst, "r+") do fid
        @showprogress 10 "to_data.meta..." for c in metacols(df)
            if c == "时间戳" && df[c].dtype.kind == "M"
                x = df[c].view("int64").div(1e9).values
            elseif c == "代码"
                n = cld(df["代码"].str.len().max(), 8) * 8
                x = from_category(df[c], MLString{n})
            else
                x = df[c].values
            end
            fid[c][:, :] = reshape(x, N, :)
        end
        fid["涨停"][:, end] = 0
        fid["跌停"][:, end] = 0
        Δt = 1024^3 ÷ (4 * F * N)
        rows = LinearIndices((1:N, 1:T))
        @showprogress 10 "to_data.fea..." for t in 1:Δt:T
            ti, tf = t, min(t + Δt - 1, T)
            ri, rf = rows[1, ti], rows[end, tf]
            slc = df.loc[(ri - 1):(rf - 2), featcols(df)]
            slc = slc.astype("float32").T.values
            x = reshape(slc, F, N, :)
            fid["特征"][:, :, ti:tf] = x
        end
    end
    @eval GC.gc()
    return dst
end

featcols(df) = sort!(setdiff(df.columns, string.(afieldnames(Data))))

metacols(df) = sort!(setdiff(df.columns, featcols(df)))

split_metafeat(df) = df[metacols(df)], df[featcols(df)]

Base.eltype(data::Data) = eltype(data.特征)

function Base.similar(data::Data, dims, dst = randstring())
    initdata(dst, eltype(data), dims, featnames(data))
    loaddata(dst, mode = "r+")
end

function Mmap.sync!(data::Data)
    for s in afieldnames(Data)
        Mmap.sync!(getfield(data, s))
    end
end

function findsnap(data::Data, hsnap)
    stamps = data.时间戳
    for t in 1:size(stamps, 2)
        zmax = 0.0
        @inbounds @simd for n in 1:size(stamps, 1)
            z = stamps[n, t]
            zmax = ifelse(isna(z), zmax, max(zmax, z))
        end
        unix2hour(zmax) > hsnap && return max(1, t - 1)
    end
    return size(stamps, 2)
end

function splitday(data::Data)
    unix = vec(maximum(fillnan, data.时间戳, dims = 1))
    is = findall(diff(unix .÷ 86400) .>= 1)
    rs = UnitRange.([1; is .+ 1], [is; length(unix)])
    return [view(data, :, r) for r in rs]
end


text:
import LinearAlgebra: det
import StaticArrays: SizedMatrix
import .DelaunayTriangulations: DelaunayTriangulation

""" Re-zeros the array `a`. """
function rezero!(a)
    @inbounds for i in eachindex(a)
        a[i] = 0.0
    end
end

""" Fill the elements of vector `from` into vector `into`. """
function fill_into!(into, from)
    @inbounds for i in eachindex(into)
        into[i] = from[i]
    end
end

""" Fill vector `into` at indices `inds` with the elements at `from[inds]`. """
function fill_at_inds!(into, from, inds)
    @inbounds for i in 1:length(inds)
        into[inds[i]] = from[i]
    end
end

function mdet_fm_in_x_view!(M, v, n, starts, stops)
    @inbounds for i = 0:(n - 1)
        M[starts[i+1]:stops[i+1]] = view(v, starts[i+1]:stops[i+1])
    end
    @fastmath det(M)
end

function iscontained!(signs, s_arr, sx, point, dim, temp_arr, starts, stops)
    # Redefine the temporary simplex. This is in-place, so we don't allocate
    # memory. We could also have re-initialised `signs`, but since we're never
    # comparing more than two consecutive signs, this is not necessary.
    rezero!(s_arr)
    rezero!(signs)
    fill_into!(s_arr, sx)

    #Replace first vertex with the point
    fill_at_inds!(s_arr, point, 1:dim)

    # Signed volume
    signs[1] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops))

    rezero!(s_arr)
    fill_into!(s_arr, sx) #reset

    for κ = 2:dim # Check remaining signs and stop if sign changes
        # Replace the ith vertex with the point we're cheking (leaving the
        # 1 appended to Vi intact.)
        idxs = ((dim + 1)*(κ - 1)+1):((dim + 1)*(κ - 1)+ 1 + dim - 1)
        fill_at_inds!(s_arr, point, idxs) # ith change

        signs[κ] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops))

        if !(signs[κ-1] == signs[κ])
            return false
        end

        rezero!(s_arr)
        fill_into!(s_arr, sx)
    end

    # Last the last vertex with the point in question
    idxs = ((dim + 1)*(dim)+1):((dim+1)^2-1)
    fill_at_inds!(s_arr, point, idxs)

    signs[end] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops))

    if !(signs[end-1] == signs[end])
       return false
    else
        return true
    end
end

function innerloop_optim!(inds::Vector{Int}, signs, s_arr, Sj, pt, dim::Int, M, i::Int, temp_arr, starts, stops)
    for j in 1:length(inds)
        if iscontained!(signs, s_arr, Sj[j], pt, dim, temp_arr, starts, stops)
            M[inds[j], i] += 1.0
        end
    end
end

"""
    reshape_simplices(pts::Vector{Vector{T}}, DT::DelaunayTriangulation) where T

Creates alternative representations of the simplices that allow for efficient
(mostly non-allocating) checking if a point lies inside a simplex.

This function adds some information needed when calculating the orientations
of the simplices *before* we start the computation (appends a row of ones on top of the
matrix representation of the simplex). Otherwise, the fallback is machinery in Simplices.jl,
which is much slower because the additional information must be added
`(n_simplices^2)*n_sample_pts` times (creating that many new matrices), which takes forever.

Reshaping the simplices and appending that information beforehand avoids excessive
memory allocation.
"""
function reshape_simplices(pts, DT::DelaunayTriangulation)
    n_simplices = length(DT.indices)
    n_vertices = length(DT[1])
    dim = n_vertices - 1

    # The [:, j, i] th entry of these two arrays holds the jth vertex of the
    # ith simplex, but instead of having just `dim` vertices, we append a `1`
    # to the end of the vectors. This allows for efficent (non-allocating)
    # computation within the `contains_point_lessalloc!` function. If we instead
    # would have appended the 1's inside that function, we would be performing
    # memory-allocating operations, which are very expensive. Doing this instead
    # gives orders of magnitude speed-ups for sufficiently large triangulations.
    S1 = Array{Float64}(undef, dim + 1, dim + 1, n_simplices)

    # Collect simplices in the form of (dim+1)^2-length column vectors. This
    # also helps with the
    NDIM, N_SIMPLICES = (dim+1)^2, n_simplices
    simplices = SizedMatrix{NDIM}{N_SIMPLICES}(zeros((NDIM, N_SIMPLICES)))

    @inbounds for i in 1:n_simplices
        for j in 1:n_vertices
            S1[:, j, i] = vcat(pts[DT[i][j]], 1.0)
        end

        simplices[:, i] = reshape(S1[:, :, i], (dim+1)^2)
    end

    return simplices
end

"""
idxs_potentially_intersecting_simplices(all_pts::Vector{Vector{T}},
        DT::DelaunayTriangulation, idx::Int) where T

Given a delaunay triangulation `DT` of `all_pts[1:(end - 1)]`, find the indices of
the simplices that potentially intersect with the image simplex with index `idx`.
"""
function idxs_potentially_intersecting_simplices(all_pts, DT::DelaunayTriangulation, idx::Int)
	# Vector that will hold the indices of the simplices potentially
	# intersecting with the image of simplex #idx
	inds_potential_simplices = Int[]

	n_simplices = length(DT)

	original_pts = all_pts[1:(end - 1)]
	forward_pts = all_pts[2:end]

	simplices = [SimplexTypes.MutableSimplex(original_pts[DT[i]]) for i = 1:n_simplices]
	image_simplices = [SimplexTypes.MutableSimplex(forward_pts[DT[i]]) for i = 1:n_simplices]

	cs = [SimplexTypes.centroid(simplices[i]) for i = 1:n_simplices]
	cs_im = [SimplexTypes.centroid(image_simplices[i]) for i = 1:n_simplices]

	rs = [SimplexTypes.radius(simplices[i]) for i = 1:n_simplices]
	rs_im = [SimplexTypes.radius(image_simplices[i]) for i = 1:n_simplices]

	@inbounds for i = 1:n_simplices
	    δ = transpose(cs_im[idx] .- cs[i]) * (cs_im[idx] .- cs[i]) .- ((rs_im[idx] + rs[i])^2)
	    if δ[1] < 0
	        push!(inds_potential_simplices, i)
	    end
	end

	return inds_potential_simplices
end

""" Get the simplices at the given indices. """
function get_simplices_at_inds!(simps, inds::Vector{Int}, simplices)
    for i in 1:length(inds)
        simps[i] = simplices[:, inds[i]]
    end
end

ext: jl


text:
using TableLogger
using Documenter

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

makedocs(;
    modules=[TableLogger],
    authors="Ole Kröger <o.kroeger@opensourc.es> and contributors",
    repo="https://github.com/Wikunia/TableLogger.jl/blob/{commit}{path}#{line}",
    sitename="TableLogger.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Wikunia.github.io/TableLogger.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Wikunia/TableLogger.jl",
    devbranch="main",
)


text:
module CutPursuit

using CxxWrap

const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LibCutPursuit not installed properly, run Pkg.build(\"CutPursuit\"), restart Julia and try again")
end
include(depsjl_path)
check_deps()

function __init__()
    @initcxx
end

# Wrap cut pursuit c++ code

@wrapmodule libcpjl

export cut_pursuit

function test()
    # obs =
    # source =
    # target =
    # edge_weight =

    cut_pursuit(rand(100,2), Int32.(rand(1:100, 100)), Int32.(rand(1:100,100)), ones(Float64,100))
end

end # module


text:
export mean_abs_error

function mean_abs_error(ys::Vector, zs::Vector, transformer=nothing)
    data = collect(Iterators.flatten(map(x -> vec(abs.(x[1]-x[2])), data)))
    score = sum(data) / length(data)
end

function mean_squared_error()
end

function mae(data::Vector, trans)
    data1 = collect(Iterators.flatten(map(x -> vec(x[1]), data)))
    data1 = reshape(data1, 1, length(data1))
    data1 = trans(data1)
    data2 = collect(Iterators.flatten(map(x -> vec(x[2]), data)))
    data2 = reshape(data2, 1, length(data2))
    data2 = trans(data2)
    data = abs.(data1 - data2)

    # data = collect(Iterators.flatten(map(x -> vec(abs.(x[1]-x[2])), data)))
    score = sum(data) / length(data)
    println("MAE:")
    display(score)
    println()
end


text:
function triple01()
#
#  test polynomial suggested by Goedecker
#
    p = triple(5,5,5)
    z = [[0.9;1;1.1]  [5;5;5]]

    p, PolyZeros(z)
end


text:
using PyCall
pushfirst!(PyVector(pyimport("sys")["path"]), "")
@pyimport src.route.pydomain as pd

"""
    haversine(lon1, lat1, lon2, lat2)

Calculate the haversine distance and bearing. Distance is in nm.
"""
function haversine(lon1::Float64, lat1::Float64, lon2::Float64, lat2::Float64)
    R = 6372.8  # Earth radius in kilometers

    dLat = deg2rad(lat2 - lat1)
    dLon = deg2rad(lon2 - lon1)
    lat1 = deg2rad(lat1)
    lat2 = deg2rad(lat2)
    lon1 = deg2rad(lon1)
    lon2 = deg2rad(lon2)
    a = sin(dLat/2)^2 + cos(lat1)*cos(lat2)*sin(dLon/2)^2
    c = 2*asin(sqrt(a))
    theta = atan(sin(dLon)*cos(lat2),
                 cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(dLon))
    theta = (rad2deg(theta) + 360) % 360
    return R*c*0.5399565, theta
end

"""
    co_ordinates(start_long, finish_long, start_lat, finish_lat,
                 n_ranks, n_nodes, dist)

Return the co-ordinates of each point across the discretized domain.
"""
function co_ordinates(start_long, finish_long, start_lat, finish_lat,
                      n_ranks::Int, n_nodes::Int, dist)
    x, y, land = pd.return_co_ords(start_long, finish_long, start_lat,
                                   finish_lat, n_ranks, n_nodes,
                                   dist)
    return x, y, land
end

function output()
    x, y, land = co_ordinates(-6.0, -14.0, 0.0, 0.0, 10, 10, 5000.0)
    print(x)
    print(y)
    print(land)
end

function min_angle(a::Float64, b::Float64)
    abs(mod(a - b + 180.0, 360.0) - 180.0)
end


text:
Base.@kwdef mutable struct ToBeOrdered
    time::Int
    type_id::Int
end

vector_to_be_ordered = Vector{ToBeOrdered}()

for i in 1:100
    push!(vector_to_be_ordered, ToBeOrdered(
        time = i,
        type_id = rand([1,2])
    ))
end

function reorder(data::Vector{ToBeOrdered})::Vector{ToBeOrdered}

    new_vector = Vector{ToBeOrdered}()

    for i in eachindex(data)
        if data[i].type_id == 2
            push!(new_vector, data[i])
            push!(new_vector, ToBeOrdered(time = data[i].time, type_id = 3))
        else
            push!(new_vector, data[i])
        end
    end
    return new_vector
end

test = reorder(vector_to_be_ordered)

using DataFrames
check = DataFrame(test)


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastLoop
using GeometricExamples

const Δt     = 2.
const ntime  = 25000
const nloop  = 200
const nsave  = 25
const nplot  = 10
const run_id = "poincare_1st_fast_vprk_dt2"

include("../guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_1st(Δt, nloop, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_1st.jl")


text:
module BridgeFinder
include("BridgeFinder_.jl")
export debridgePE
export debridgeSE
export splice_mates
export revcompFQ
export readlength_summary
export readlength_SE_summary
end

max_stars_repo_name: straightlab/chartools-dev

max_issues_repo_name: straightlab/chartools-dev

max_forks_repo_name: straightlab/chartools-dev


text:
using Documenter, TravelingSalesmanExact, Asciicast

asciinema_version = "2.6.1"
makedocs(;
    modules=[TravelingSalesmanExact],
    authors="Eric P. Hanson",
    repo="https://github.com/ericphanson/TravelingSalesmanExact.jl/blob/{commit}{path}#L{line}",
    sitename="TravelingSalesmanExact.jl",
    format=Documenter.HTML(;
        ansicolor=true,
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ericphanson.github.io/TravelingSalesmanExact.jl",
        assets=[asset("https://cdnjs.cloudflare.com/ajax/libs/asciinema-player/$(asciinema_version)/asciinema-player.min.js"),
                asset("https://cdnjs.cloudflare.com/ajax/libs/asciinema-player/$(asciinema_version)/asciinema-player.min.css")],
    ),
    pages=[
        "Home" => "index.md",
    ],
    strict=:doctest,
)

# Very hacky fix to load the asciinema JS before Documenter's require.js
# <https://github.com/JuliaDocs/Documenter.jl/issues/1433>
requires_regex = r"""<script src="https://cdnjs\.cloudflare\.com/ajax/libs/require\.js/[0-9]+\.[0-9]+\.[0-9]+/require\.min\.js" data-main="assets/documenter\.js"></script>"""
asciinema_script = """<script src="https://cdnjs.cloudflare.com/ajax/libs/asciinema-player/$(asciinema_version)/asciinema-player.min.js"></script>"""
index_path = joinpath(@__DIR__, "build", "index.html")
index = read(index_path, String)
requires_script = match(requires_regex, index).match # get the right version numbers for requires.js
index = replace(index, asciinema_script => "")
index = replace(index, requires_script => asciinema_script*requires_script)
write(index_path, index)

deploydocs(;
    repo="github.com/ericphanson/TravelingSalesmanExact.jl",
    push_preview=true,
    devbranch="main"
)


text:
using SHA

if isempty(ARGS)
    error("need file to test sha perf")
elseif !isfile(ARGS[1])
    error("file $(ARGS[1]) does not exist")
end

function do_tests(filepath)
    # test performance
    print("read:    ")
    @time begin
        const fh = open(filepath, "r")
        const bytes = readbytes(fh)
    end
    gc()

    print("SHA-1:   ")
    sha1(bytes)
    gc()
    @time sha1(bytes)

    print("SHA2-256: ")
    sha256(bytes)
    gc()
    @time sha256(bytes)

    print("SHA2-512: ")
    sha512(bytes)
    gc()
    @time sha512(bytes)

    print("SHA3-256: ")
    sha3_256(bytes)
    gc()
    @time sha3_256(bytes)

    print("SHA3-512: ")
    sha3_512(bytes)
    gc()
    @time sha3_512(bytes)
end

do_tests(ARGS[1])


text:
using AlphaZero

depth = 5
gspec = Examples.games["tictactoe"]
computer = MinMax.Player(depth=depth, amplify_rewards=true, τ=0.2)
interactive!(gspec, computer, Human())


text:
# This type allows to compute length once and for all
struct StringWithLength{T <: AbstractString} <: AbstractString
    s::T
    l::Int
end
string_with_length(s::AbstractString) = StringWithLength(s, length(s))
Base.length(s::StringWithLength) = s.l
Base.iterate(s::StringWithLength, i::Integer = firstindex(s.s)) = iterate(s.s, i)
Base.nextind(s::StringWithLength, i::Int, n::Int = 1) = nextind(s.s, i, n)
Base.ncodeunits(s::StringWithLength) = ncodeunits(s.s)
Base.isvalid(s::StringWithLength, i::Int) = isvalid(s.s, i)

function reorder(s1::AbstractString, s2::AbstractString)
    s1 = string_with_length(s1)
    s2 = string_with_length(s2)
    (length(s1) <= length(s2)) ? (s1, s2) : (s2, s1)
end

function reorder(s1, s2)
    (length(s1) <= length(s2)) ? (s1, s2) : (s2, s1)
end

function common_prefix(s1, s2)
    x1 = iterate(s1)
    x2 = iterate(s2)
    l = 0
    while (x1 !== nothing) & (x2 !== nothing)
        ch1, state1 = x1
        ch2, state2 = x2
        ch1 != ch2 && break
        l += 1
        x1 = iterate(s1, state1)
        x2 = iterate(s2, state2)
    end
    return l, x1, x2
end

function _take(s, n::Integer)
    Base.Iterators.take(s, n)
end
function _take(s::AbstractString, n::Integer)
   SubString(s, firstindex(s), nextind(s, 0, n))
end

function _drop(s, n::Integer)
    Base.Iterators.drop(s, n)
end
function _drop(s::AbstractString, n::Integer)
   SubString(s, nextind(s, 0, n + 1), lastindex(s))
end

function _slice(s, n1::Integer, n2::Integer)
    Base.Iterators.take(Base.Iterators.drop(s, n1), n2 - n1)
end
function _slice(s::AbstractString, n1::Integer, n2::Integer)
   SubString(s, nextind(s, 0, n1 + 1),  nextind(s, 0, n2))
end


text:
using Test, ImageFeatures, Images, TestImages, Distributions

@testset "Generating brief params" begin
    brief_params = BRIEF(size = 8, window = 3, seed = 123)
    @test brief_params.size == 8
    @test brief_params.window == 3
    @test brief_params.seed == 123
    @test brief_params.sampling_type == gaussian
    @test brief_params.sigma == 2 ^ 0.5
end

@testset "Sampling Patterns" begin
    # Check Bounds
    for si in [32, 128, 256], wi in [5, 15, 25, 49], se in [123, 546, 178]
        s1, s2 = ImageFeatures.center_sample(si, wi, se)
        @test check_samples(s1, s2, si, wi) == true
        @test s1 == zeros(CartesianIndex{2}, si)
        s1, s2 = random_coarse(si, wi, se)
        @test check_samples(s1, s2, si, wi) == true
        s1, s2 = random_uniform(si, wi, se)
        @test check_samples(s1, s2, si, wi) == true
        s1, s2 = gaussian(si, wi, se)
        @test check_samples(s1, s2, si, wi) == true
        s1, s2 = gaussian_local(si, wi, se)
        @test check_samples(s1, s2, si, wi) == true
    end
end

@testset "Descriptor Calculation" begin
    img_1 = Gray{Float64}[  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
                            0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  ]
    keypoints_1 = Keypoints(imcorner(img_1))

    img_2 = Gray{Float64}[ 0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0
                            0 0 1 1 1 0 0 0 0 0
                            0 0 1 1 1 1 0 0 0 0
                            0 0 1 1 1 1 1 0 0 0
                            0 0 1 1 1 1 1 1 0 0
                            0 0 1 1 1 1 1 1 0 0
                            0 0 1 1 1 1 1 1 0 0
                            0 0 0 0 0 0 0 0 0 0
                            0 0 0 0 0 0 0 0 0 0 ]
    keypoints_2 = Keypoints(imcorner(img_2))

    brief_params = BRIEF(size = 8, window = 3, seed = 123)

    desc_1, ret_keypoints_1 = create_descriptor(img_1, keypoints_1, brief_params)
    desc_2, ret_keypoints_2 = create_descriptor(img_2, keypoints_2, brief_params)
    matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2)
    expected_matches = [[CartesianIndex((5,4)),CartesianIndex((3,3))],
                        [CartesianIndex((16,4)),CartesianIndex((8,3))],
                        [CartesianIndex((16,15)),CartesianIndex((8,8))]]
    @test all(matches .== expected_matches)

    img_1 = Gray{Float64}[ 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 1 0 0 0
                           0 0 1 1 1 1 0 0 0
                           0 0 1 1 1 1 0 0 0
                           0 0 1 1 1 1 0 0 0
                           0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0 ]
    keypoints_1 = Keypoints(imcorner(img_1))

    img_2 = Gray{Float64}[ 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 1 1 1 0
                           0 0 1 1 1 1 1 1 0
                           0 0 1 1 1 1 1 1 0
                           0 0 1 1 1 1 1 1 0
                           0 0 1 1 1 1 1 1 0
                           0 0 1 1 1 1 1 1 0
                           0 0 0 0 0 0 0 0 0 ]
    keypoints_2 = Keypoints(imcorner(img_2))

    brief_params = BRIEF(size = 8, window = 3, seed = 123)

    desc_1, ret_keypoints_1 = create_descriptor(img_1, keypoints_1, brief_params)
    desc_2, ret_keypoints_2 = create_descriptor(img_2, keypoints_2, brief_params)
    matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2)
    expected_matches = [[CartesianIndex((2,3)),CartesianIndex((2,3))],
                        [CartesianIndex((5,3)),CartesianIndex((7,3))],
                        [CartesianIndex((2,6)),CartesianIndex((2,8))],
                        [CartesianIndex((5,6)),CartesianIndex((7,8))]]
    @test all(matches .== expected_matches)

    img_1 = Gray{Float64}[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
    keypoints_1 = Keypoints(imcorner(img_1))

    img_2 = Gray{Float64}[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0
                           0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
                           0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
    keypoints_2 = Keypoints(imcorner(img_2))

    brief_params = BRIEF(size = 128, window = 5, seed = 123)

    desc_1, ret_keypoints_1 = create_descriptor(img_1, keypoints_1, brief_params)
    desc_2, ret_keypoints_2 = create_descriptor(img_2, keypoints_2, brief_params)
    matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.2)
    expected_matches = [[CartesianIndex{2}((3,3)),CartesianIndex{2}((9,10))],
                        [CartesianIndex{2}((6,3)),CartesianIndex{2}((13,10))],
                        [CartesianIndex{2}((3,7)),CartesianIndex{2}((9,15))],
                        [CartesianIndex{2}((6,7)),CartesianIndex{2}((13,15))],
                        [CartesianIndex{2}((8,13)),CartesianIndex{2}((3,4))],
                        [CartesianIndex{2}((13,13)),CartesianIndex{2}((7,4))]]
    @test all(matches .== expected_matches)

end

@testset "Testing with Standard Images - Lighthouse (Translation (100, 200))" begin
    img = testimage("lighthouse")
    img_array_1 = convert(Array{Gray}, img)
    img_array_2 = _warp(img_array_1, 100, 200)

    keypoints_1 = Keypoints(fastcorners(img_array_1, 12, 0.4))
    keypoints_2 = Keypoints(fastcorners(img_array_2, 12, 0.4))

    brief_params = BRIEF(size = 256, window = 10, seed = 123)

    desc_1, ret_keypoints_1 = create_descriptor(img_array_1, keypoints_1, brief_params)
    desc_2, ret_keypoints_2 = create_descriptor(img_array_2, keypoints_2, brief_params)
    matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.1)
    @test sum(m[1] + CartesianIndex(100, 200) == m[2] for m in matches) + 1 == length(matches)
end

@testset "Testing with Standard Images - Lena (Translation (10, 20))" begin
    img = testimage("lena_gray_512")
    img_array_1 = convert(Array{Gray}, img)
    img_array_2 = _warp(img_array_1, 10, 20)

    keypoints_1 = Keypoints(fastcorners(img_array_1, 12, 0.4))
    keypoints_2 = Keypoints(fastcorners(img_array_2, 12, 0.4))

    brief_params = BRIEF(size = 256, window = 10, seed = 123)

    desc_1, ret_keypoints_1 = create_descriptor(img_array_1, keypoints_1, brief_params)
    desc_2, ret_keypoints_2 = create_descriptor(img_array_2, keypoints_2, brief_params)
    matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.1)
    @test all(m[1] + CartesianIndex(10, 20) == m[2] for m in matches)
end


text:
#=

# 8: Water collapse (explicit, symplectic, and reversible)

```@raw html
	<img src='../assets/fixpa.png' alt='missing' width="50%" height="50%" /><br>
```

Simulation of a water column collapsing under its own weight onto dry bottom.
Here we use a symplectic scheme and get a reversible simulation.
At the end of the simulation, the velocities are reverted and the simulation goes back to its initial conditions.
Despite the reversibility, Boltzmann entropy grows and attains its maximum value just before the velocities are reverted.
=#

module collapse_symplectic

using Printf
using SmoothedParticles
using Parameters
using Plots
using DataFrames # to store the csv file
using CSV# to store the csv file
include("utils/FixPA.jl")
include("utils/entropy.jl")
using .FixPA
using .entropy

#using ReadVTK  #not implemented
#using VTKDataIO

#=
Declare constant parameters
=#

##physical
const dr = 1.0e-2          #average particle distance (decrease to make finer simulation)
const h = 3.0*dr           #size of kernel support
const rho0 = 1000.   	   #fluid density
const m = rho0*dr^2        #particle mass
const g = -9.8*VECY  #gravitational acceleration
const mu = 0.0#8.4e-4          #dynamic viscosity of water

##geometrical
const water_column_width = 1.0
const water_column_height = 2.0
const box_height = 3.0
const box_width = 4.0
const wall_width = 2.5*dr

##artificial
const c = 50.0             #numerical speed of sound
const dr_wall = 0.95*dr
const E_wall = 10*norm(g)*water_column_height
const eps = 1e-16

##temporal
const dt = 0.1*h/c
const t_end = 1.0
const dt_frame = t_end/100

##particle types
const FLUID = 0.
const WALL = 1.

@with_kw mutable struct Particle <: AbstractParticle
	x::RealVector #position
	v::RealVector = VEC0 #velocity
	a::RealVector = VEC0 #acceleration
	P::Float64 = 0. #pressure
	rho::Float64 = 0. #density
    rho0::Float64 = 0.
	type::Float64 #particle_type
end

#=
Define geometry and make particles
=#

function make_system()
	grid = Grid(dr, :square)
	box = Rectangle(0., 0., box_width, box_height)
	fluid = Rectangle(0., 0., water_column_width, water_column_height)
	walls = BoundaryLayer(box, grid, wall_width)
	#walls = Specification(walls, x -> (x[2] < box_height))
	domain = Rectangle(-box_width, -box_width, 2*box_width, 3*box_height)
	sys = ParticleSystem(Particle, domain, h)
	generate_particles!(sys, grid, fluid, x -> Particle(x = x, type = FLUID))
	generate_particles!(sys, grid, walls, x -> Particle(x = x, type = WALL))
	return sys
end

#=
Define particle interactions
=#

@inbounds function find_rho!(p::Particle, q::Particle, r::Float64)
    if p.type == FLUID && q.type == FLUID
		p.rho += m*wendland2(h,r)
	end
end

@inbounds function find_rho0!(p::Particle, q::Particle, r::Float64)
    if p.type == FLUID && q.type == FLUID
		p.rho0 += m*wendland2(h,r)
	end
end

function find_pressure!(p::Particle)
	p.P = c^2*(p.rho - p.rho0)
end

@inbounds function internal_force!(p::Particle, q::Particle, r::Float64)
	if p.type == FLUID && q.type == FLUID
		ker = m*rDwendland2(h,r)
		p.a += -ker*(p.P/rho0^2 + q.P/rho0^2)*(p.x - q.x)
		#p.a += +2*ker*mu/rho0^2*(p.v - q.v)
	elseif p.type == FLUID && q.type == WALL && r < dr_wall
		s = dr_wall/(r + eps)
		p.a += -E_wall/(r + eps)^2*(s^2 - s^4)*(p.x - q.x)
	end
end

function reset_a!(p::Particle)
    p.a = zero(RealVector)
end

function reset_rho!(p::Particle)
    p.rho = 0.0
end

function move!(p::Particle)
	if p.type == FLUID
		p.x = rev_add(p.x, dt*p.v)
	end
end

function accelerate!(p::Particle)
	if p.type == FLUID
		p.v = rev_add(p.v, 0.5*dt*(p.a + g))
	end
end

function LJ_potential(p::Particle, q::Particle, r::Float64)::Float64
	if q.type == WALL && p.type == FLUID && r < dr_wall
		s = dr_wall/(r + eps)
		return m*E_wall*(0.5s^2 - 0.25s^4 -0.25)
	else
		return 0.0
	end
end

function energy_kinetic(sys::ParticleSystem)::Float64
	return sum(p -> 0.5*m*dot(p.v, p.v), sys.particles)
end

function energy(sys::ParticleSystem, p::Particle)::Float64
	kinetic = 0.5*m*dot(p.v, p.v)
	internal =  0.5*m*c^2*(p.rho - p.rho0)^2/rho0^2
	gravity_potential = -m*dot(g, p.x)
	wall_potential = SmoothedParticles.sum(sys, LJ_potential, p)
	return kinetic + internal + gravity_potential + wall_potential
end

#=
Put everything into a time loop
=#

function verlet_step!(sys::ParticleSystem)
    apply!(sys, accelerate!)
    apply!(sys, move!)
    create_cell_list!(sys)
    apply!(sys, reset_rho!)
    apply!(sys, find_rho!, self = true)
    apply!(sys, find_pressure!)
    apply!(sys, reset_a!)
    apply!(sys, internal_force!)
    apply!(sys, accelerate!)
end

function save_results!(out::SmoothedParticles.DataStorage, sys::ParticleSystem, k::Int64)
    if (k %  Int64(round(dt_frame/dt)) == 0)
        @printf("t = %.6e\n", k*dt)
        #energy
        E = sum(p -> energy(sys,p), sys.particles)
        @show E
        println("# of part. = ", length(sys.particles))
        println()
        save_frame!(out, sys, :v, :a, :P, :rho, :rho0)
    end
end

function main(;revert = true) #if revert=true, velocities are inverted at the end of the simulation and the simulation then goes backward
	sys = make_system()
	out = new_pvd_file("results/collapse_fixpa")
    #initialization
    create_cell_list!(sys)
    apply!(sys, find_rho0!, self = true)
    apply!(sys, find_rho!, self = true)
    apply!(sys, find_pressure!)
    apply!(sys, internal_force!)

	N_of_particles = length(sys.particles)
	@show(N_of_particles)
	@show(m)

	step_final = Int64(round(t_end/dt))
	times = Float64[] #time instants
	Ss = Float64[] # Entropy values
	Ekin = Float64[] # Kinetic energy values
	for k = 0 : step_final
        verlet_step!(sys)
        save_results!(out, sys, k)
    	if k % round(step_final/100) == 0 # store a number of entropy values
			distr = velocity_histogram(sys, N = 100)
			S = entropy_2D_MB(distr)
			push!(times, k*dt)
			push!(Ss, S)
			push!(Ekin, energy_kinetic(sys))
			@show(S)
        	println()
		end
	end

	# Plotting the velocity distribution in comparison with Maxwell-Boltzmann
	T = plot_velocity_distr(sys, m, "energy_distribution_middle.pdf")

	# Plotting the entropy in time
	Sred_eq_E = [(1+log(Ekin[k]/(m*length(sys.particles)))) for k in 1:length(Ss)]
	Sred_eq_T= (1+log(kB*T/m))*ones(Float64, length(Ss))
	p = plot(times, [Ss Sred_eq_T Sred_eq_E], label = ["entropy" "S_eq(T)" "S_eq(E)"],legend=:bottomright)
	savefig(p, "entropy_middle.pdf")
	df = DataFrame(time_steps = times, S_Boltzmann = Ss, S_eq_T = Sred_eq_T, S_eq_E = Sred_eq_E)
	CSV.write("entropy_middle.csv", df)

	if revert
		#revert velocities
		println("--------------------")
		println("Reverting velocities")
		println("--------------------")
		for p in sys.particles
			p.v = -p.v
		end
		Ss_rev = Float64[]
		for k = step_final:-1:0
			verlet_step!(sys)
			save_results!(out, sys, k)
			if k % round(step_final/100) == 0 # store a number of entropy values
				distr = velocity_histogram(sys, v_max = sqrt(2*norm(g)*water_column_height), N = 100)
				S = entropy_2D_MB(distr)
				push!(Ss_rev, S)
				@show(S)
				println()
			end
		end
		plot_velocity_distr(sys, m, "energy_distribution_final.pdf")

		# Plotting the entropy in time
		p = plot(times, [Ss Ss_rev Sred_eq_T Sred_eq_E], label = ["entropy forward" "entropy backward" "S_eq(T)" "S_eq(E)"], legend=:bottomright)
		savefig(p, "entropy_final.pdf")
		df = DataFrame(time_steps = times, S_Boltzmann = Ss, S_eq_T = Sred_eq_T, S_eq_E = Sred_eq_E)
		CSV.write("entropy_final.csv", df)
	end

	save_pvd_file(out)

end ## function main

end ## module


text:
# ## 2.4 グラフのプロット

## #src
#-
using CSV
using DataFrames
using Measures
using Plots

## #src
#-
data = CSV.read("../data/house-prices-advanced-regression-techniques/train.csv", DataFrame)

## #src
#-
scatter(data[!, :GrLivArea], data[!, :SalePrice], label="", margin = 12mm)
xlabel!("GrLivArea")
ylabel!("SalePrice")

## #src
#-
scatter(data[!, :MSSubClass], data[!, :SalePrice], label="", margin = 12mm)
xlabel!("MSSubClass")
ylabel!("SalePrice")


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --project=@. --color=yes --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

import GuillotineModels.Data.read_from_file
using GuillotineModels.PPG2KP.Heuristic
using TimerOutputs
using RandomNumbers.Xorshifts # for Xoroshiro128Plus (for heuristic rng)
import Profile

function main(args)
	timer = TimerOutput()
	# Parsing.
	qt_args = length(args)
	if qt_args < 1 || qt_args > 3
		println("usage: ./run_heuristic instance_filename [rng_seed] [max_iter]")
		println("For now the instance type is fixed to `Val{:Classic_G2KP}`.")
		exit()
	end
	filename = args[1]
	seed = (qt_args > 1 ? parse(Int, args[2]) : 1)
	max_iter = (qt_args > 2 ? parse(Int, args[3]) : 1000000)

	# Initializing.
	instance = read_from_file(Val(:Classic_G2KP), filename)
	N, L, W, l, w, p, d = length(instance.l), instance.L, instance.W,
		instance.l, instance.w, instance.p, instance.d

	#GC.enable(false)
	# Running.
	rng = Xoroshiro128Plus(seed)
	bkv2, sel2, pat2 = @timeit timer "mock_fast" begin
		bkv, sel, pat = fast_iterated_greedy(
			d, p, l, w, L, W, rng, max_iter
		)
		#GC.gc()
		bkv, sel, pat
	end
	rng = Xoroshiro128Plus(seed)
	bkv1, sel1, pat1 = @timeit timer "mock_normal" begin
		bkv, sel, pat = iterated_greedy(
			d, p, l, w, L, W, rng, max_iter
		)
		#GC.gc()
		bkv, sel, pat
	end
	if (bkv1, sel1, pat1) == (bkv2, sel2, pat2)
		println("The two mocks give the same values.")
		println("bkv = $(bkv1)")
		println("sel = $(sel1)")
		println("pat = $(pat1)")
	else
		println("The two mocks give distinct values. They follow:")
		println("Note: two pieces of the same width changing positions is normal.")
		println("iterated_greedy return:")
		@show bkv1
		@show sel1
		@show pat1
		println("fast_iterated_greedy return:")
		@show bkv2
		@show sel2
		@show pat2
	end

	Profile.clear_malloc_data()
	rng = Xoroshiro128Plus(seed)
	bkv1, sel1, pat1 = @timeit timer "normal" begin
		bkv, sel, pat = iterated_greedy(
			d, p, l, w, L, W, rng, max_iter
		)
		#GC.gc()
		bkv, sel, pat
	end
	rng = Xoroshiro128Plus(seed)
	bkv2, sel2, pat2 = @timeit timer "fast" begin
		bkv, sel, pat = fast_iterated_greedy(
			d, p, l, w, L, W, rng, max_iter
		)
		#GC.gc()
		bkv, sel, pat
	end
	if (bkv1, sel1, pat1) == (bkv2, sel2, pat2)
		println("The two runs give the same values. They follow:")
		println("bkv = $(bkv1)")
		println("sel = $(sel1)")
		println("pat = $(pat1)")
	else
		println("The two runs give distinct values. They follow:")
		println("Note: two pieces of the same width changing positions is normal.")
		println("iterated_greedy return:")
		@show bkv1
		@show sel1
		@show pat1
		println("fast_iterated_greedy return:")
		@show bkv2
		@show sel2
		@show pat2
	end
	TimerOutputs.print_timer(timer)
	println()
end

main(ARGS)


text:
"""
    optZeros = ds_optzeros(n, opt=1)

A helper function for the synthesizeNTF function. Returns the zeros which
minimize the in-band noise power of a delta-sigma modulator's NTF.
"""
function ds_optzeros(n, opt=1)
    if opt == 0
        optZeros = zeros(ceil(Int, n/2))
    else
        if n == 1
            optZeros = [0.0]
        elseif n == 2
            if opt == 1
                optZeros = [sqrt(1/3)]
            else
                optZeros = [0.0]
            end
        elseif n == 3
            optZeros = [sqrt(3/5), 0]
        elseif n == 4
            if opt == 1
                discr = sqrt(9/49 - 3/35)
                tmp = 3/7
                optZeros = sqrt.([tmp+discr, tmp-discr])
            else
                optZeros = [0, sqrt(5/7)]
            end
        elseif n == 5
            discr = sqrt(25/81 - 5/21)
            tmp = 5/9
            optZeros = sqrt.([tmp+discr, tmp-discr, 0])
        elseif n == 6
            if opt == 1
                optZeros = [0.23862059, 0.66120988, 0.9324696]
            else
                discr = sqrt(56)/33
                tmp = 7/11;
                optZeros = sqrt.([0, tmp+discr, tmp-discr])
            end
        elseif n == 7
            optZeros = [0, 0.40584371, 0.74153078, 0.94910785]
        elseif n == 8
            if opt == 1
                optZeros = [0.18343709, 0.52553345, 0.79666684, 0.96028993]
            else
                optZeros = [0, 0.50563161, 0.79017286, 0.95914731]
            end
        elseif n == 9
            optZeros = [0, 0.32425101, 0.61337056, 0.83603082, 0.9681602]
        elseif n == 10
            if opt == 1
                optZeros = [0.1834370913, 0.5255334458, 0.7966668433, 0.9602899327]
            else
                optZeros = [0, 0.41572267, 0.67208682, 0.86238894, 0.97342121]
            end
        elseif n == 11
            optZeros = [0, 0.26953955, 0.51909468, 0.73015137, 0.88706238, 0.97822864]
        elseif n == 12
            if opt == 1
                optZeros = [0.12523875, 0.36783403, 0.58731921, 0.7699033, 0.90411753, 0.9815607]
            else
                optZeros = [0, 0.35222363, 0.58006251, 0.76647993, 0.90281326, 0.98132047]
            end
        elseif n == 13
            optZeros = [0, 0.23045331, 0.44849063, 0.64234828, 0.8015776, 0.91759824, 0.98418306]
        elseif n == 14
            if opt == 1
                optZeros = [0.10806212, 0.31911586, 0.51525046, 0.68729392, 0.82720185, 0.92843513, 0.98628389]
            else
                optZeros = [0, 0.30524384, 0.50836649, 0.6836066, 0.82537239, 0.92772336, 0.98615167]
            end
        else
            error("Optimized zeros for n>14 are not available")
            return
        end
    end

    # Sort the zeros and replicate them.
    z = sort(optZeros)
    optZeros = zeros(n)
    m = 1
    if isodd(n)
        optZeros[1] = z[1]
        z = z[2:length(z)]
        m = m + 1
    end
    for i = 1:length(z)
        optZeros[m]   =  z[i]
        optZeros[m+1] = -z[i]
        m = m + 2
    end

    return optZeros
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
""" testing the input handling
"""

import Particula as P
using Unitful, Test

@testset "℞X" begin

    @testset "in_temperature" begin

        @test P.℞T(300) == 300 * u"K"
        @test P.℞T(300 * u"K") == 300 * u"K"
        @test_throws ArgumentError P.℞T(1e-5 * u"kg/m")

    end

    @testset "in_viscosity" begin

        @test P.℞µ(1e-5) == 1e-5 * u"kg/m/s"
        @test P.℞µ(1e-5 * u"Pa*s") == 1e-5 * u"kg/m/s"
        @test_throws ArgumentError P.℞µ(1e-5 * u"kg/m")

    end

end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    SatelliteToolbox orbit propagator API: propagate_to_epoch!
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export propagate_to_epoch!

"""
    function propagate_to_epoch!(orbp, JD::Number) where T
    function propagate_to_epoch!(orbp, JD::AbstractVector) where T

If `t` is a number, then propagate `orbp` until the epoch `JD` [Julian Day].
Otherwise, if `JD` is an array, then propagate the orbit in `orbp` using the
epochs defined in the vector `t` [Julian Day].

In both cases, the orbit propagator algorithm is the one related to the
structure `orbp`.

The structure `orbp` will contain the elements at the last propagation instant.

# Returns

* The mean Keplerian elements represented in inertial frame in each time instant
  (see `Orbit`) [SI units].
* The position vector represented in inertial frame in each time instant [m].
* The velocity vector represented in inertial frame in each time instant [m].

If `JD` is an array, then those values will be an array containing the
information related to each epoch in `JD`.

# Remarks

The inertial frame in which the output is represented depends on which frame it
was used to generate the orbit parameters. If the orbit parameters are obtained
from a TLE, then the inertial frame will be TEME. Notice, however, that the
perturbation theory requires an inertial frame with true equator.

"""
function propagate_to_epoch!(orbp::OrbitPropagatorJ2{T},
                             JD::Union{Number,AbstractVector}) where T
    propagate!(orbp, (JD .- orbp.j2d.epoch)*86400)
end

function propagate_to_epoch!(orbp::OrbitPropagatorSGP4{T},
                             JD::Union{Number,AbstractVector}) where T
    propagate!(orbp, (JD .- orbp.sgp4d.epoch)*86400)
end

function propagate_to_epoch!(orbp::OrbitPropagatorTwoBody{T},
                             JD::Union{Number,AbstractVector}) where T
    propagate!(orbp, (JD .- orbp.tbd.epoch)*86400)
end


text:
module ColorTypes

using Color, FixedPointNumbers, Base.Cartesian
import Color: Fractional, _convert

import Base: ==, abs, abs2, clamp, convert, div, isfinite, isinf,
    isnan, isless, length, one, promote_array_type, promote_rule, zero,
    trunc, floor, round, ceil, itrunc, ifloor, iceil, iround, bswap,
    mod, rem

export ARGB, BGR, RGB1, RGB4, BGRA, AbstractGray, Gray, GrayAlpha, Gray24, AGray32, YIQ, AlphaColor, ColorType, noeltype

typealias ColorType Union(ColorValue, AbstractAlphaColorValue)

# An alpha-channel-first memory layout
immutable AlphaColor{C<:ColorValue, T<:Fractional} <: AbstractAlphaColorValue{C,T}
    alpha::T
    c::C

    AlphaColor(x1::Real, x2::Real, x3::Real, alpha::Real = 1.0) = new(alpha, C(x1, x2, x3))
    AlphaColor(c::ColorValue, alpha::Real) = new(alpha, c)
end
AlphaColor{T<:Fractional}(c::ColorValue{T}, alpha::T = one(T)) = AlphaColor{typeof(c),T}(c, alpha)

typealias ARGB{T} AlphaColor{RGB{T}, T}

# Little-endian RGB (useful for BGRA & Cairo)
immutable BGR{T<:Fractional} <: AbstractRGB{T}
    b::T
    g::T
    r::T

    BGR(r::Real, g::Real, b::Real) = new(b, g, r)
end
BGR(r::Integer, g::Integer, b::Integer) = BGR{Float64}(r, g, b)
BGR(r::Fractional, g::Fractional, b::Fractional) = (T = promote_type(typeof(r), typeof(g), typeof(b)); BGR{T}(r, g, b))

typealias BGRA{T} AlphaColorValue{BGR{T}, T}

# Some readers return a byte for an alpha channel even if it's not meaningful
immutable RGB1{T<:Fractional} <: AbstractRGB{T}
    alphadummy::T
    r::T
    g::T
    b::T

    RGB1(r::Real, g::Real, b::Real) = new(one(T), r, g, b)
end
RGB1(r::Integer, g::Integer, b::Integer) = RGB1{Float64}(r, g, b)
RGB1(r::Fractional, g::Fractional, b::Fractional) = (T = promote_type(typeof(r), typeof(g), typeof(b)); RGB1{T}(r, g, b))

immutable RGB4{T<:Fractional} <: AbstractRGB{T}
    r::T
    g::T
    b::T
    alphadummy::T

    RGB4(r::Real, g::Real, b::Real) = new(r, g, b, one(T))
end
RGB4(r::Integer, g::Integer, b::Integer) = RGB4{Float64}(r, g, b)
RGB4(r::Fractional, g::Fractional, b::Fractional) = (T = promote_type(typeof(r), typeof(g), typeof(b)); RGB4{T}(r, g, b))

for CV in subtypes(AbstractRGB)
    CV == RGB && continue
    @eval begin
        convert{T}(::Type{$CV{T}}, c::$CV{T}) = c
        convert(::Type{$CV}, c::$CV) = c
        convert{T}(::Type{$CV{T}}, c::$CV) = $CV{T}(convert(T, c.r), convert(T, c.g), convert(T, c.b))
        convert{T<:Fractional}(::Type{$CV}, c::ColorValue{T}) = convert($CV{T}, c)
        convert(::Type{$CV}, c::AbstractRGB) = $CV(c.r, c.g, c.b)
        convert{T}(::Type{$CV{T}}, c::AbstractRGB) = $CV{T}(c.r, c.g, c.b)
    end
end
convert(::Type{RGB}, c::AbstractRGB) = RGB(c.r, c.g, c.b)
convert{T}(::Type{RGB{T}}, c::AbstractRGB) = RGB{T}(c.r, c.g, c.b)

# Sometimes you want to be explicit about grayscale. Also needed for GrayAlpha.
abstract AbstractGray{T} <: ColorValue{T}
immutable Gray{T<:Fractional} <: AbstractGray{T}
    val::T
end
convert{T}(::Type{Gray{T}}, x::Gray{T}) = x
convert{T,S}(::Type{Gray{T}}, x::Gray{S}) = Gray{T}(x.val)
convert{T<:Real}(::Type{T}, x::Gray) = convert(T, x.val)
convert{T}(::Type{Gray{T}}, x::Real) = Gray{T}(x)

convert{T}(::Type{Gray{T}}, x::AbstractRGB) = convert(Gray{T}, 0.299*x.r + 0.587*x.g + 0.114*x.b)  # Rec 601 luma conversion

zero{T}(::Type{Gray{T}}) = Gray{T}(zero(T))
 one{T}(::Type{Gray{T}}) = Gray{T}(one(T))

immutable Gray24 <: ColorValue{Uint8}
    color::Uint32
end
Gray24() = Gray24(0)
Gray24(val::Uint8) = (g = uint32(val); g<<16 | g<<8 | g)
Gray24(val::Ufixed8) = Gray24(reinterpret(val))

convert(::Type{Uint32}, g::Gray24) = g.color

typealias GrayAlpha{T} AlphaColorValue{Gray{T}, T}

immutable AGray32 <: AbstractAlphaColorValue{Gray24, Uint8}
    color::Uint32
end
AGray32() = AGray32(0)
AGray32(val::Uint8, alpha::Uint8) = (g = uint32(val); uint32(alpha)<<24 | g<<16 | g<<8 | g)
AGray32(val::Ufixed8, alpha::Ufixed8) = AGray32(reinterpret(val), reinterpret(alpha))

convert(::Type{Uint32}, g::AGray32) = g.color

convert(::Type{RGB}, x::Gray) = RGB(x.val, x.val, x.val)
convert{T}(::Type{RGB{T}}, x::Gray) = (g = convert(T, x.val); RGB{T}(g, g, g))

# YIQ (NTSC)
immutable YIQ{T<:FloatingPoint} <: ColorValue{T}
    y::T
    i::T
    q::T

    YIQ(y::Real, i::Real, q::Real) = new(y, i, q)
end
YIQ(y::FloatingPoint, i::FloatingPoint, q::FloatingPoint) = (T = promote_type(typeof(y), typeof(i), typeof(q)); YIQ{T}(y, i, q))

clamp{T}(c::YIQ{T}) = YIQ{T}(clamp(c.y, zero(T), one(T)),
                             clamp(c.i, convert(T,-0.5957), convert(T,0.5957)),
                             clamp(c.q, convert(T,-0.5226), convert(T,0.5226)))

function convert{T}(::Type{YIQ{T}}, c::AbstractRGB)
    rgb = clamp(c)
    YIQ{T}(0.299*rgb.r+0.587*rgb.g+0.114*rgb.b,
           0.595716*rgb.r-0.274453*rgb.g-0.321263*rgb.b,
           0.211456*rgb.r-0.522591*rgb.g+0.311135*rgb.b)
end
convert{T}(::Type{YIQ}, c::AbstractRGB{T}) = convert(YIQ{T}, c)

function _convert{T}(::Type{RGB{T}}, c::YIQ)
    cc = clamp(c)
    RGB{T}(cc.y+0.9563*cc.i+0.6210*cc.q,
           cc.y-0.2721*cc.i-0.6474*cc.q,
           cc.y-1.1070*cc.i+1.7046*cc.q)
end

# Y'CbCr
immutable YCbCr{T<:FloatingPoint} <: ColorValue{T}
    y::T
    cb::T
    cr::T

    YCbCr(y::Real, cb::Real, cr::Real) = new(y, cb, cr)
end
function YCbCr(y::FloatingPoint, cb::FloatingPoint, cr::FloatingPoint)
    T = promote_type(typeof(y), typeof(cb), typeof(cr))
    YCbCr{T}(y, cb, cr)
end

clamp{T}(c::YCbCr{T}) = YCbCr{T}(clamp(c.y, convert(T,16), convert(T,235)),
                                 clamp(c.cb, convert(T,16), convert(T,240)),
                                 clamp(c.cr, convert(T,16), convert(T,240)))

function convert{T}(::Type{YCbCr{T}}, c::AbstractRGB)
    rgb = clamp(c)
    YCbCr{T}(16+65.481*rgb.r+128.553*rgb.g+24.966*rgb.b,
             128-37.797*rgb.r-74.203*rgb.g+112*rgb.b,
             128+112*rgb.r-93.786*rgb.g-18.214*rgb.b)
end
convert{T}(::Type{YCbCr}, c::AbstractRGB{T}) = convert(YCbCr{T}, c)

function _convert{T}(::Type{RGB{T}}, c::YCbCr)
    cc = clamp(c)
    ny = cc.y - 16
    ncb = cc.cb - 128
    ncr = cc.cr - 128
    RGB{T}(0.004567ny - 1.39135e-7ncb + 0.0062586ncr,
           0.004567ny - 0.00153646ncb - 0.0031884ncr,
           0.004567ny + 0.00791058ncb - 2.79201e-7ncr)
end

# HSI
immutable HSI{T<:FloatingPoint} <: ColorValue{T}
    h::T
    s::T
    i::T

    HSI(h::Real, s::Real, i::Real) = new(h, s, i)
end
function HSI(h::FloatingPoint, s::FloatingPoint, i::FloatingPoint)
    T = promote_type(typeof(h), typeof(s), typeof(i))
    HSI{T}(h, s, i)
end

function convert{T}(::Type{HSI{T}}, c::AbstractRGB)
    rgb = clamp(c)
    α = (2rgb.r - rgb.g - rgb.b)/2
    β = 0.8660254*(rgb.g - rgb.b)
    h = atan2(β, α)
    i = (rgb.r + rgb.g + rgb.b)/3
    s = 1-min(rgb.r, rgb.g, rgb.b)/i
    s = ifelse(i > 0, s, zero(s))
    HSI{T}(h, s, i)
end
convert{T}(::Type{HSI}, c::AbstractRGB{T}) = convert(HSI{T}, c)

# TODO: HSI->RGB

## Generic algorithms

length(cv::ColorType) = div(sizeof(cv), sizeof(eltype(cv)))
# Because this can be called as `length(RGB)`, we might need to fill in a default element type.
# But the compiler chokes if we ask it to create RGB{Float64}{Float64}, even if that's inside
# the non-evaluated branch of a ternary expression, so we have to be sneaky about this.
length{CV<:ColorValue}(::Type{CV}) = _length(CV, eltype(CV))
_length{CV<:ColorValue}(::Type{CV}, ::Type{Any}) = length(CV{Float64})
_length{CV<:ColorValue}(::Type{CV}, ::DataType)  = div(sizeof(CV), sizeof(eltype(CV)))
length{CV,T}(::Type{AlphaColorValue{CV,T}}) = length(CV)+1
length{CV,T}(::Type{AlphaColor{CV,T}}) = length(CV)+1

for ACV in (ColorValue, AbstractRGB, AbstractGray, AbstractAlphaColorValue)
    for CV in subtypes(ACV)
        length(CV.parameters) == 1 || continue
        @eval noeltype{T}(::Type{$CV{T}}) = $CV
    end
end

for f in (:trunc, :floor, :round, :ceil, :itrunc, :ifloor, :iceil, :iround, :abs, :abs2, :isfinite, :isnan, :isinf, :bswap)
    @eval $f{T}(g::Gray{T}) = Gray{T}($f(g.val))
    @eval @vectorize_1arg Gray $f
end

for f in (:mod, :rem, :mod1)
    @eval $f(x::Gray, m::Gray) = Gray($f(x.val, m.val))
    @eval @vectorize_2arg Gray $f
end

# Return types for arithmetic operations
multype(a::Type,b::Type) = typeof(one(a)*one(b))
sumtype(a::Type,b::Type) = typeof(one(a)+one(b))
divtype(a::Type,b::Type) = typeof(one(a)/one(b))

# Math on ColorValues. These implementations encourage inlining and,
# for the case of Ufixed types, nearly halve the number of multiplications.
for CV in subtypes(AbstractRGB)
    @eval begin
        (*){R<:Real,T}(f::R, c::$CV{T}) = $CV{multype(R,T)}(f*c.r, f*c.g, f*c.b)
        (*){R<:Real,T}(f::R, c::AlphaColorValue{$CV{T},T}) =
            AlphaColorValue{$CV{multype(R,T)},multype(R,T)}(f*c.c.r, f*c.c.g, f*c.c.b, f*c.alpha)
        (*){R<:Real,T}(f::R, c::AlphaColor{$CV{T},T}) =
            AlphaColor{$CV{multype(R,T)},multype(R,T)}(f*c.c.r, f*c.c.g, f*c.c.b, f*c.alpha)
        function (*){R<:FloatingPoint,T<:Ufixed}(f::R, c::$CV{T})
            fs = f/reinterpret(one(T))
            $CV{multype(R,T)}(fs*reinterpret(c.r), fs*reinterpret(c.g), fs*reinterpret(c.b))
        end
        function (*){R<:Ufixed,T<:Ufixed}(f::R, c::$CV{T})
            fs = reinterpret(f)/widen(reinterpret(one(T)))^2
            $CV{multype(R,T)}(fs*reinterpret(c.r), fs*reinterpret(c.g), fs*reinterpret(c.b))
        end
        function (/){R<:FloatingPoint,T<:Ufixed}(c::$CV{T}, f::R)
            fs = one(R)/(f*reinterpret(one(T)))
            $CV{divtype(R,T)}(fs*reinterpret(c.r), fs*reinterpret(c.g), fs*reinterpret(c.b))
        end
        (+){S,T}(a::$CV{S}, b::$CV{T}) = $CV{sumtype(S,T)}(a.r+b.r, a.g+b.g, a.b+b.b)
        (-){S,T}(a::$CV{S}, b::$CV{T}) = $CV{sumtype(S,T)}(a.r-b.r, a.g-b.g, a.b-b.b)
        (+){S,T}(a::AlphaColorValue{$CV{S},S}, b::AlphaColorValue{$CV{T},T}) =
            AlphaColorValue{$CV{sumtype(S,T)},sumtype(S,T)}(a.c.r+b.c.r, a.c.g+b.c.g, a.c.b+b.c.b, a.alpha+b.alpha)
        (+){S,T}(a::AlphaColor{$CV{S},S}, b::AlphaColor{$CV{T},T}) =
            AlphaColor{$CV{sumtype(S,T)},sumtype(S,T)}(a.c.r+b.c.r, a.c.g+b.c.g, a.c.b+b.c.b, a.alpha+b.alpha)
        (-){S,T}(a::AlphaColorValue{$CV{S},S}, b::AlphaColorValue{$CV{T},T}) =
            AlphaColorValue{$CV{sumtype(S,T)},sumtype(S,T)}(a.c.r-b.c.r, a.c.g-b.c.g, a.c.b-b.c.b, a.alpha-b.alpha)
        (-){S,T}(a::AlphaColor{$CV{S},S}, b::AlphaColor{$CV{T},T}) =
            AlphaColor{$CV{sumtype(S,T)},sumtype(S,T)}(a.c.r-b.c.r, a.c.g-b.c.g, a.c.b-b.c.b, a.alpha-b.alpha)
        function (.+){T}(A::AbstractArray{$CV{T}}, b::AbstractRGB)
            bT = convert($CV{T}, b)
            out = similar(A)
            add!(out, A, bT)
        end
        function (.+){T}(A::Union(AbstractArray{AlphaColorValue{$CV{T},T}}, AbstractArray{AlphaColor{$CV{T},T}}), b::AbstractAlphaColorValue)
            bT = convert($CV{T}, b)
            out = similar(A)
            add!(out, A, bT)
        end
        (.+){T}(b::AbstractRGB, A::AbstractArray{$CV{T}}) = (.+)(A, b)
        (.+){T}(b::AbstractAlphaColorValue, A::Union(AbstractArray{AlphaColorValue{$CV{T},T}}, AbstractArray{AlphaColor{$CV{T},T}})) =
            (.+)(A, b)
        function (.-){T}(A::AbstractArray{$CV{T}}, b::AbstractRGB)
            bT = convert($CV{T}, b)
            out = similar(A)
            sub!(out, A, bT)
        end
        function (.-){T}(A::Union(AbstractArray{AlphaColorValue{$CV{T},T}}, AbstractArray{AlphaColor{$CV{T},T}}), b::AbstractAlphaColorValue)
            bT = convert($CV{T}, b)
            out = similar(A)
            sub!(out, A, bT)
        end
        function (.-){T}(b::AbstractRGB, A::AbstractArray{$CV{T}})
            bT = convert($CV{T}, b)
            out = similar(A)
            sub!(out, bT, A)
        end
        function (.-){T}(b::AbstractAlphaColorValue, A::Union(AbstractArray{AlphaColorValue{$CV{T},T}}, AbstractArray{AlphaColor{$CV{T},T}}))
            bT = convert($CV{T}, b)
            out = similar(A)
            sub!(out, bT, A)
        end
        isfinite{T<:Ufixed}(c::$CV{T}) = true
        isfinite{T<:Ufixed}(c::AbstractAlphaColorValue{$CV{T},T}) = true
        isfinite{T<:FloatingPoint}(c::$CV{T}) = isfinite(c.r) && isfinite(c.g) && isfinite(c.b)
        isfinite{T<:FloatingPoint}(c::AbstractAlphaColorValue{$CV{T},T}) = isfinite(c.c.r) && isfinite(c.c.g) && isfinite(c.c.b) && isfinite(c.alpha)
        isnan{T<:Ufixed}(c::$CV{T}) = false
        isnan{T<:Ufixed}(c::AbstractAlphaColorValue{$CV{T},T}) = false
        isnan{T<:FloatingPoint}(c::$CV{T}) = isnan(c.r) || isnan(c.g) || isnan(c.b)
        isnan{T<:FloatingPoint}(c::AbstractAlphaColorValue{$CV{T},T}) = isnan(c.c.r) || isnan(c.c.g) || isnan(c.c.b) || isnan(c.alpha)
        isinf{T<:Ufixed}(c::$CV{T}) = false
        isinf{T<:Ufixed}(c::AbstractAlphaColorValue{$CV{T},T}) = false
        isinf{T<:FloatingPoint}(c::$CV{T}) = isinf(c.r) || isinf(c.g) || isinf(c.b)
        isinf{T<:FloatingPoint}(c::AbstractAlphaColorValue{$CV{T},T}) = isinf(c.c.r) || isinf(c.c.g) || isinf(c.c.b) || isinf(c.alpha)
        abs(c::$CV) = abs(c.r)+abs(c.g)+abs(c.b) # should this have a different name?
        abs{T<:Ufixed}(c::$CV{T}) = float32(c.r)+float32(c.g)+float32(c.b) # should this have a different name?
        one{T<:ColorType}(::T) = one(T)
        one{T}(::Type{$CV{T}}) = $CV{T}(one(T),one(T),one(T))
        one{T}(::Type{AlphaColorValue{$CV{T},T}}) = AlphaColorValue{$CV{T},T}(one(T),one(T),one(T),one(T))
        one{T}(::Type{AlphaColor{$CV{T},T}}) = AlphaColorValue{$CV{T},T}(one(T),one(T),one(T),one(T))
        zero{T<:ColorType}(::T) = zero(T)
        zero{T}(::Type{$CV{T}}) = $CV{T}(zero(T),zero(T),zero(T))
        zero{T}(::Type{AlphaColorValue{$CV{T},T}}) = AlphaColorValue{$CV{T},T}(zero(T),zero(T),zero(T),zero(T))
        zero{T}(::Type{AlphaColor{$CV{T},T}}) = AlphaColorValue{$CV{T},T}(zero(T),zero(T),zero(T),zero(T))
    end
end
(*)(c::AbstractRGB, f::Real) = (*)(f, c)
(*){CV<:AbstractRGB}(c::AlphaColorValue{CV}, f::Real) = (*)(f, c)
(*){CV<:AbstractRGB}(c::AlphaColor{CV}, f::Real) = (*)(f, c)
(.*)(f::Real, c::AbstractRGB) = (*)(f, c)
(.*){CV<:AbstractRGB}(f::Real, c::AlphaColorValue{CV}) = (*)(f, c)
(.*){CV<:AbstractRGB}(f::Real, c::AlphaColor{CV}) = (*)(f, c)
(.*)(c::AbstractRGB, f::Real) = (*)(f, c)
(.*){CV<:AbstractRGB}(c::AlphaColorValue{CV}, f::Real) = (*)(f, c)
(.*){CV<:AbstractRGB}(c::AlphaColor{CV}, f::Real) = (*)(f, c)
(/)(c::AbstractRGB, f::Real) = (one(f)/f)*c
(/){CV<:AbstractRGB}(c::AlphaColorValue{CV}, f::Real) = (one(f)/f)*c
(/){CV<:AbstractRGB}(c::AlphaColor{CV}, f::Real) = (one(f)/f)*c
(/)(c::AbstractRGB, f::Integer) = (one(eltype(c))/f)*c
(/){CV<:AbstractRGB}(c::AlphaColorValue{CV}, f::Integer) = (one(eltype(c))/f)*c
(/){CV<:AbstractRGB}(c::AlphaColor{CV}, f::Integer) = (one(eltype(c))/f)*c
(./)(c::AbstractRGB, f::Real) = (/)(c, f)
(./){CV<:AbstractRGB}(c::AlphaColorValue{CV}, f::Real) = (/)(c, f)
(./){CV<:AbstractRGB}(c::AlphaColor{CV}, f::Real) = (/)(c, f)
(+){CV<:AbstractRGB}(A::AbstractArray{CV}, b::AbstractRGB) = (.+)(A, b)
(+){CV<:AbstractRGB,T,CV2<:AbstractRGB}(A::AbstractArray{AlphaColorValue{CV,T}}, b::AbstractAlphaColorValue{CV2}) = (.+)(A, b)
(+){CV<:AbstractRGB,T,CV2<:AbstractRGB}(A::AbstractArray{AlphaColor{CV,T}}, b::AbstractAlphaColorValue{CV2}) = (.+)(A, b)
(-){CV<:AbstractRGB}(A::AbstractArray{CV}, b::AbstractRGB) = (.-)(A, b)
(-){CV<:AbstractRGB,T,CV2<:AbstractRGB}(A::AbstractArray{AlphaColorValue{CV,T}}, b::AbstractAlphaColorValue{CV2}) = (.-)(A, b)
(-){CV<:AbstractRGB,T,CV2<:AbstractRGB}(A::AbstractArray{AlphaColor{CV,T}}, b::AbstractAlphaColorValue{CV2}) = (.-)(A, b)
(+){CV<:AbstractRGB}(b::AbstractRGB, A::AbstractArray{CV}) = (.+)(b, A)
(+){CV<:AbstractRGB,T,CV2<:AbstractRGB}(b::AbstractAlphaColorValue{CV2}, A::AbstractArray{AlphaColorValue{CV,T}}) = (.+)(A, b)
(+){CV<:AbstractRGB,T,CV2<:AbstractRGB}(b::AbstractAlphaColorValue{CV2}, A::AbstractArray{AlphaColor{CV,T}}) = (.+)(A, b)
(-){CV<:AbstractRGB}(b::AbstractRGB, A::AbstractArray{CV}) = (.-)(b, A)
(-){CV<:AbstractRGB,T,CV2<:AbstractRGB}(b::AbstractAlphaColorValue{CV2}, A::AbstractArray{AlphaColorValue{CV,T}}) = (.-)(A, b)
(-){CV<:AbstractRGB,T,CV2<:AbstractRGB}(b::AbstractAlphaColorValue{CV2}, A::AbstractArray{AlphaColor{CV,T}}) = (.-)(A, b)

# Math on Gray
for CV in subtypes(AbstractGray)
    @eval begin
        (*){R<:Real,T}(f::R, c::$CV{T}) = $CV{multype(R,T)}(f*c.val)
        (*)(c::$CV, f::Real) = (*)(f, c)
        (.*)(f::Real, c::$CV) = (*)(f, c)
        (.*)(c::$CV, f::Real) = (*)(f, c)
        (/)(c::$CV, f::Real) = (one(f)/f)*c
        (/)(c::$CV, f::Integer) = (one(eltype(c))/f)*c
        (./)(c::$CV, f::Real) = (/)(c, f)
        (+){S,T}(a::$CV{S}, b::$CV{T}) = $CV{sumtype(S,T)}(a.val+b.val)
        (-){S,T}(a::$CV{S}, b::$CV{T}) = $CV{sumtype(S,T)}(a.val-b.val)
        (*){S,T}(a::$CV{S}, b::$CV{T}) = $CV{multype(S,T)}(a.val*b.val)
        (+)(A::AbstractArray{$CV}, b::AbstractGray) = (.+)(A, b)
        (-)(A::AbstractArray{$CV}, b::AbstractGray) = (.-)(A, b)
        (+)(b::AbstractGray, A::AbstractArray{$CV}) = (.+)(b, A)
        (-)(b::AbstractGray, A::AbstractArray{$CV}) = (.-)(b, A)
        function (.+){T}(A::AbstractArray{$CV{T}}, b::AbstractGray)
            bT = convert($CV{T}, b)
            out = similar(A)
            add!(out, A, bT)
        end
        (.+){T}(b::AbstractGray, A::AbstractArray{$CV{T}}) = (.+)(A, b)
        function (.-){T}(A::AbstractArray{$CV{T}}, b::AbstractGray)
            bT = convert($CV{T}, b)
            out = similar(A)
            sub!(out, A, bT)
        end
        function (.-){T}(A::AbstractArray{$CV{T}}, b::AbstractGray)
            bT = convert($CV{T}, b)
            out = similar(A)
            sub!(out, A, bT)
        end
        function (.-){T}(b::AbstractGray, A::AbstractArray{$CV{T}})
            bT = convert($CV{T}, b)
            out = similar(A)
            sub!(out, bT, A)
        end
        (*)(A::AbstractArray, b::AbstractGray) = A .* b
        (*)(b::AbstractGray, A::AbstractArray) = b .* A
        (/)(A::AbstractArray, b::AbstractGray) = A ./ b
        function (.*){T,S}(A::AbstractArray{$CV{T}}, b::AbstractGray{S})
            Tout = multype(T,S)
            out = similar(A, $CV{Tout})
            mul!(out, A, b)
        end
        (.*){T,S}(b::AbstractGray{S}, A::AbstractArray{$CV{T}}) = A .* b
        function (./){T,S}(A::AbstractArray{$CV{T}}, b::AbstractGray{S})
            Tout = divtype(T,S)
            out = similar(A, $CV{Tout})
            div!(out, A, b)
        end
        isfinite{T<:Ufixed}(c::$CV{T}) = true
        isfinite{T<:FloatingPoint}(c::$CV{T}) = isfinite(c.val)
        isnan{T<:Ufixed}(c::$CV{T}) = false
        isnan{T<:FloatingPoint}(c::$CV{T}) = isnan(c.val)
        isinf{T<:Ufixed}(c::$CV{T}) = false
        isinf{T<:FloatingPoint}(c::$CV{T}) = isinf(c.val)
        abs(c::$CV) = abs(c.val) # should this have a different name?
        abs{T<:Ufixed}(c::$CV{T}) = float32(c.val) # should this have a different name?

        (<)(c::$CV, r::Real) = c.val < r
        (<)(r::Real, c::$CV) = r < c.val
        isless(c::$CV, r::Real) = c.val < r
        isless(r::Real, c::$CV) = r < c.val
        (<)(a::$CV, b::AbstractGray) = a.val < b.val
    end
end
(/)(a::AbstractGray, b::AbstractGray) = a.val/b.val
div(a::AbstractGray, b::AbstractGray) = div(a.val, b.val)
(+)(a::AbstractGray, b::Number) = a.val+b
(-)(a::AbstractGray, b::Number) = a.val-b
(+)(a::Number, b::AbstractGray) = a+b.val
(-)(a::Number, b::AbstractGray) = a-b.val
(.+)(a::AbstractGray, b::Number) = a.val+b
(.-)(a::AbstractGray, b::Number) = a.val-b
(.+)(a::Number, b::AbstractGray) = a+b.val
(.-)(a::Number, b::AbstractGray) = a-b.val

@ngenerate N typeof(out) function add!{T,N}(out, A::AbstractArray{T,N}, b)
    @inbounds begin
        @nloops N i A begin
            @nref(N, out, i) = @nref(N, A, i) + b
        end
    end
    out
end
# need a separate sub! because of unsigned types
@ngenerate N typeof(out) function sub!{T,N}(out, A::AbstractArray{T,N}, b::ColorType)  # TODO: change to b::T when julia #8045 fixed
    @inbounds begin
        @nloops N i A begin
            @nref(N, out, i) = @nref(N, A, i) - b
        end
    end
    out
end
@ngenerate N typeof(out) function sub!{T,N}(out, b::ColorType, A::AbstractArray{T,N})
    @inbounds begin
        @nloops N i A begin
            @nref(N, out, i) = b - @nref(N, A, i)
        end
    end
    out
end
@ngenerate N typeof(out) function mul!{T,N}(out, A::AbstractArray{T,N}, b)
    @inbounds begin
        @nloops N i A begin
            @nref(N, out, i) = @nref(N, A, i) * b
        end
    end
    out
end
@ngenerate N typeof(out) function div!{T,N}(out, A::AbstractArray{T,N}, b)
    @inbounds begin
        @nloops N i A begin
            @nref(N, out, i) = @nref(N, A, i) / b
        end
    end
    out
end

# To help type inference
for ACV in (ColorValue, AbstractRGB, AbstractGray)
    for CV in subtypes(ACV)
        (length(CV.parameters) == 1 && !(CV.abstract)) || continue
        @eval promote_array_type{T<:Real,S<:Real}(::Type{T}, ::Type{$CV{S}}) = $CV{promote_type(T, S)}
        @eval promote_rule{T<:Fractional,S<:Fractional}(::Type{$CV{T}}, ::Type{$CV{S}}) = $CV{promote_type(T, S)}
        for AC in subtypes(AbstractAlphaColorValue)
            (length(AC.parameters) == 2 && !(AC.abstract)) || continue
            @eval promote_array_type{T<:Real,S<:Real}(::Type{T}, ::Type{$AC{$CV{S},S}}) = (TS = promote_type(T, S); $AC{$CV{TS}, TS})
            @eval promote_rule{T<:Fractional,S<:Fractional}(::Type{$CV{T}}, ::Type{$CV{S}}) = $CV{promote_type(T, S)}
            @eval promote_rule{T<:Fractional,S<:Integer}(::Type{$CV{T}}, ::Type{S}) = $CV{promote_type(T, S)} # for Array{RGB}./Array{Int}
        end
    end
end

for (CV, CVstr, fields) in ((BGR,  "BGR",  (:(c.r),:(c.g),:(c.b))),
                            (RGB1, "RGB1", (:(c.r),:(c.g),:(c.b))),
                            (RGB4, "RGB4", (:(c.r),:(c.g),:(c.b))),
                            (ARGB, "ARGB", (:(c.c.r),:(c.c.g),:(c.c.b),:(c.alpha))),
                            (BGRA, "BGRA", (:(c.c.r),:(c.c.g),:(c.c.b),:(c.alpha))),
                            (Gray, "Gray", (:(c.val),)),
                            (GrayAlpha, "GrayAlpha", (:(c.c.val),:(c.alpha))))
    Color.makeshow(CV, CVstr, fields)
end

for T in (RGB24, ARGB32, Gray24, AGray32)
    @eval begin
        ==(x::Uint32, y::$T) = x == convert(Uint32, y)
        ==(x::$T, y::Uint32) = ==(y, x)
    end
end
=={T}(x::Gray{T}, y::Gray{T}) = x.val == y.val
=={T}(x::T, y::Gray{T}) = x == convert(T, y)
=={T}(x::Gray{T}, y::T) = ==(y, x)

end


text:
@testset "Clebsch-Gordan" begin
using PyCall, Test, SphericalHarmonics

sympy = pyimport("sympy")
spin = pyimport("sympy.physics.quantum.spin")

pycg(j1, m1, j2, m2, j3, m3, T=Float64) =
      spin.CG(j1, m1, j2, m2, j3, m3).doit().evalf().__float__()

@info("Testing cg1 implementation against sympy ... ")
for j1 = 0:2, j2=0:2, j3=0:4
   for m1 = -j1:j1, m2=-j2:j2, m3=-j3:j3
      @test cg1(j1,m1,j2,m2,j3,m3) ≈ pycg(j1,m1, j2,m2, j3,m3)
   end
end

@info("Check CG Coefficients for some higher quantum numbers...")
j1 = 8
j2 = 11
j3 = j1+j2
for m1 = -j1:j1, m2=-j2:j2, m3=-j3:j3
   @test cg1(j1,m1,j2,m2,j3,m3) ≈ pycg(j1,m1, j2,m2, j3,m3)
end

@info("Checking the SphH expansion in terms of CG coeffs")
# expansion coefficients of a product of two spherical harmonics in terms a
# single spherical harmonic
# see e.g. https://en.wikipedia.org/wiki/Clebsch–Gordan_coefficients
# this is the magic formula that we need
for ntest = 1:10
   # two random Ylm  ...
   l1, l2 = rand(1:10), rand(1:10)
   m1, m2 = rand(-l1:l1), rand(-l2:l2)
   # ... evaluated at random spherical coordinates
   θ = rand() * π
   φ = (rand()-0.5) * 2*π
   # evaluate all relevant Ylms (up to l1 + l2)
   Ylm = compute_y(l1+l2, cos(θ), φ)
   # evaluate the product p = Y_l1_m1 * Y_l2_m2
   p = Ylm[index_y(l1,  m1)] * Ylm[index_y(l2,m2)]
   # and its expansion in terms of CG coeffs
   p2 = 0.0
   M = m1 + m2  # all other coeffs are zero

   for L = abs(M):(l1+l2)
      p2 += sqrt( (2*l1+1)*(2*l2+1) / (4 * π * (2*L+1)) ) *
            cg1(l1,  0, l2,  0, L, 0) *
            cg1(l1, m1, l2, m2, L, M) *
            Ylm[index_y(L, M)]
   end
   print(@test p ≈ p2); print(" ")
end
println()

end # @testset


text:
using Documenter, PlotDocs, Plots, PlotThemes

# Set matplotlib gui backend
ENV["MPLBACKEND"] = "agg"

# Initialize backends
gr()
plotlyjs()
pyplot()
pgfplotsx()
unicodeplots()
inspectdr()
gaston()

plotthemes_path = dirname(dirname(pathof(PlotThemes)))

cp(
    joinpath(plotthemes_path, "README.md"),
    joinpath(@__DIR__, "src", "generated", "plotthemes.md"),
    force = true,
)

const PAGES = Any[
    "Home" => "index.md",
    "Getting Started" => [
        "Installation" => "install.md",
        "Basics" => "basics.md",
        "Tutorial" => "tutorial.md",
    ],
    "Manual" => [
        "Input Data" => "input_data.md",
        "Output" => "output.md",
        "Attributes" => [
            "Overview" => "attributes.md",
            "Series Attributes" => "generated/attributes_series.md",
            "Plot Attributes" => "generated/attributes_plot.md",
            "Subplot Attributes" => "generated/attributes_subplot.md",
            "Axis Attributes" => "generated/attributes_axis.md",
        ],
        "Layouts" => "layouts.md",
        "Recipes" => "recipes.md",
        "Colors" => "colors.md",
        "ColorSchemes" => "generated/colorschemes.md",
        "Animations" => "animations.md",
        "Themes" => "generated/plotthemes.md",
        "Backends" => "backends.md",
        "Supported Attributes" => "generated/supported.md",
    ],
    "Learning" => "learning.md",
    "Contributing" => "contributing.md", # TODO: testing
    "Ecosystem" => [
        "Overview" => "ecosystem.md",
        "GraphRecipes" => [
            "Introduction" => "graphrecipes/introduction.md",
            "Examples" => "graphrecipes/examples.md",
            "Attributes" => "generated/graph_attributes.md",
        ],
    ],
    "Advanced Topics" => ["Internals" => "pipeline.md"],
    "Examples" => [
        "GR" => "generated/gr.md",
        "PlotlyJS" => "generated/plotlyjs.md",
        "PyPlot" => "generated/pyplot.md",
        "PGFPlotsX" => "generated/pgfplotsx.md",
        "UnicodePlots" => "generated/unicodeplots.md",
        "InspectDR" => "generated/inspectdr.md",
        "Gaston" => "generated/gaston.md",
    ],
]

generate_attr_markdown()
generate_supported_markdown()
generate_graph_attr_markdown()
generate_colorschemes_markdown()
for be in (:gr, :plotlyjs, :pyplot, :pgfplotsx, :unicodeplots, :inspectdr, :gaston)
    generate_markdown(be)
end
ansicolor = get(ENV, "PLOTDOCS_ANSICOLOR", "true") == "true"
@show ansicolor
@time makedocs(
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/favicon.ico"],
        ansicolor = ansicolor,
    ),
    sitename = "Plots",
    authors = "Thomas Breloff",
    pages = PAGES,
)

deploydocs(
    repo = "github.com/JuliaPlots/PlotDocs.jl.git",
    push_preview = true,
    versions = ["stable" => "v^", "v#.#", "dev" => "dev", "latest" => "dev"],
    forcepush = true
)


text:
#! /usr/bin/julia

# Rosetta Code, Compound data type

type Point{T<:Real}
    x::T
    y::T
end

==(u::Point, v::Point) = (u.x == v.x) & (u.y == v.y)
-(u::Point) = Point(-u.x, -u.y)
+(u::Point, v::Point) = Point(u.x + v.x, u.y + v.y)
-(u::Point, v::Point) = Point(u.x - v.x, u.y - v.y)

a = Point(1, 2)
b = Point(3, 7)
c = Point(2, 4)

println("a = ", a)
println("b = ", b)
println("c = ", c)

println("a + b =     ", a+b)
println("-a + b =    ", -a+b)
println("a - b =     ", a-b)
println("a + b + c = ", a+b+c)
println("a == c      ", a == c)
println("a + a == c  ", a + a == c)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module QUBOConstraints

# SECTION - Usings and imports
using ConstraintDomains
using DataFrames
using Flux
using LinearAlgebra
using PrettyTables

import Flux.Optimise: update!
import Flux: params
import DataFrames: describe

# SECTION - exports
export binarize
export debinarize
export is_valid
export train

export QUBO_linear_sum

# SECTION - includes
include("common.jl")
include("base.jl")

include("handmade/linear_sum.jl")

include("encoding/domain_wall.jl")
include("encoding/one_hot.jl")
include("encoding/conversion.jl")

include("learn.jl")

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/variables.jl
# Defines single variable constraints supported by optimizer and how to store.
#############################################################################

##### Access variable information from MOI variable index
has_upper_bound(m::Optimizer, vi::MOI.VariableIndex) = m._input_problem._variable_info[vi.value].has_upper_bound
has_lower_bound(m::Optimizer, vi::MOI.VariableIndex) = m._input_problem._variable_info[vi.value].has_lower_bound
is_fixed(m::Optimizer, vi::MOI.VariableIndex) = m._input_problem._variable_info[vi.value].is_fixed
is_integer(m::Optimizer, i::Int64) = is_integer(m._input_problem._variable_info[i])

##### Add unconstrained variables
function MOI.add_variable(m::Optimizer)
    m._input_problem._variable_count += 1
    push!(m._input_problem._variable_info, VariableInfo())
    return VI(m._input_problem._variable_count)
end
MOI.add_variables(m::Optimizer, n::Int) = [MOI.add_variable(m) for i in 1:n]

##### Supports function and add_constraint for single variable functions
const INEQ_SETS = Union{LT, GT, ET}
MOI.supports_constraint(::Optimizer, ::Type{SV}, ::Type{S}) where {S <: INEQ_SETS} = true

#=
function MOI.add_constraint(m::Optimizer, v::SV, zo::ZO)
    vi = v.variable
    check_inbounds!(m, vi)
    has_upper_bound(m, vi) && error("Upper bound on variable $vi already exists.")
    has_lower_bound(m, vi) && error("Lower bound on variable $vi already exists.")
    is_fixed(m, vi) && error("Variable $vi is fixed. Cannot also set upper bound.")
    m._input_problem._variable_info[vi.value].lower_bound = 0.0
    m._input_problem._variable_info[vi.value].upper_bound = 1.0
    m._input_problem._variable_info[vi.value].has_lower_bound = true
    m._input_problem._variable_info[vi.value].has_upper_bound = true
    m._input_problem._variable_info[vi.value].is_integer = true
    return CI{SV, ZO}(vi.value)
end
=#

function MOI.add_constraint(m::Optimizer, v::SV, lt::LT)
    vi = v.variable
    check_inbounds!(m, vi)
    if isnan(lt.upper)
        error("Invalid upper bound value $(lt.upper).")
    end
    if has_upper_bound(m, vi)
        error("Upper bound on variable $vi already exists.")
    end
    if is_fixed(m, vi)
        error("Variable $vi is fixed. Cannot also set upper bound.")
    end
    m._input_problem._variable_info[vi.value].upper_bound = lt.upper
    m._input_problem._variable_info[vi.value].has_upper_bound = true
    return CI{SV, LT}(vi.value)
end

function MOI.add_constraint(m::Optimizer, v::SV, gt::GT)
    vi = v.variable
    check_inbounds!(m, vi)
    if isnan(gt.lower)
        error("Invalid lower bound value $(gt.lower).")
    end
    if has_lower_bound(m, vi)
        error("Lower bound on variable $vi already exists.")
    end
    if is_fixed(m, vi)
        error("Variable $vi is fixed. Cannot also set lower bound.")
    end
    m._input_problem._variable_info[vi.value].lower_bound = gt.lower
    m._input_problem._variable_info[vi.value].has_lower_bound = true
    return CI{SV, GT}(vi.value)
end

function MOI.add_constraint(m::Optimizer, v::SV, eq::ET)
    vi = v.variable
    check_inbounds!(m, vi)
    if isnan(eq.value)
        error("Invalid fixed value $(gt.lower).")
    end
    if has_lower_bound(m, vi)
        error("Variable $vi has a lower bound. Cannot be fixed.")
    end
    if has_upper_bound(m, vi)
        error("Variable $vi has an upper bound. Cannot be fixed.")
    end
    if is_fixed(m, vi)
        error("Variable $vi is already fixed.")
    end
    m._input_problem._variable_info[vi.value].lower_bound = eq.value
    m._input_problem._variable_info[vi.value].upper_bound = eq.value
    m._input_problem._variable_info[vi.value].has_lower_bound = true
    m._input_problem._variable_info[vi.value].has_upper_bound = true
    m._input_problem._variable_info[vi.value].is_fixed = true
    return CI{SV, ET}(vi.value)
end


text:
module TEST_TILE_MAP
using NRES
using Test

struct Dummy_type
    some_id
    some_more
    function Dummy_type(arg=nothing)
        if !isnothing(arg)
            some_id = arg
        else
            some_id = new(rand())
        end
        new(some_id, Nothing)
    end
end

@testset "Constructing Tile_map with dimensjonality ZERO" begin
    case = NRES.Tile_map{Dummy_type}()
    @test isa(case, NRES.Representation)
    " NRES.Tile_map with empty argument "

    @test isnothing(case._range)
    " NRES with empty argument becomes a range-less zero-dimentional NRES "

    @test dimention_of_NRES(case) == 0
    " Range-less NRES: a zero-dimentional Tile_map (a point-representation) "

    @test isa(NRES.map_to_SAT(case), Dummy_type)
    " map_to_SAT(case) returns the one SAT of the 0dim NRES "

    @test isa(case._all_SAT, Array)
    " Even zero-dimensional NRES maps have an array of SAT "

    @test length(case._all_SAT) == 1
    " A zero-dimentional tile-map has exactly on SAT in this array "

    @test isa(case._all_SAT[1], Dummy_type)
    " Point-NRES contains one element of the supplied type T "

    case = NRES.Tile_map{Dummy_type}(the_sat = Dummy_type())
    @test isa(case._all_SAT[1], Dummy_type)
    " Tile map can be constructed by argument "

    the_dummy = Dummy_type()
    case = NRES.Tile_map{Dummy_type}(the_sat = the_dummy)
    @test case._all_SAT[1] === the_dummy
    " Equivalence: the supplied SAT becomes the 0-dim NRES consitional "

    NRES.Tile_map{Dummy_type}(the_sat=Dummy_type()) # OK
    " Spesified SAT (a single conditional) implies dim=0. If both args are specified: ArgumentError "

    @test case._number_of_intervals == 1
    " monosatian NRES have a single SAT => N=1 "
end

@testset "Constructing Tile_map with dimensjonality ONE" begin
    @testset "Range" begin
        steps_N = 2
        @test_throws MethodError NRES.Tile_map{Dummy_type}(:anna_enn_nothing_eller_tuple, N=steps_N)
        @test_throws MethodError NRES.Tile_map{Dummy_type}((:ikkje_Real ,:endpoints), N=steps_N)
        " Range skal være av type Tuple{<:Real, <:Real} : viser endepunkta på lineær akse "

        some_bolean_sat = Dummy_type()
        @test_throws ArgumentError NRES.Tile_map{Dummy_type}((0,1), N=steps_N, the_sat=some_bolean_sat)
        " If the_sat is spesified, i.e. NRES with a single SAT, specifying a range trows ArgumentError "

        case = NRES.Tile_map{Dummy_type}( (0.,2.), N=steps_N ) #[(0.,1.)])
        @test case._range == (0,2)
        " First argument is range: can be a tuple (A,B) that signifies from A to B "

        @test_throws MethodError NRES.Tile_map{Dummy_type}((1,2,3), N=steps_N)
        " Range med meir enn to endepunkt kaster MethodError "

        @test case._number_of_intervals == 2
        " How many steps N? "

        @test dimention_of_NRES(case) == 1
        " dimention_of_NRES with range should be one "

        @test_throws ArgumentError NRES.map_to_SAT(case)
        " Map NRES to SAT in NRES with range but no supplied coord? Error! "

    end
    " Testset som ser på contrucksjon: Range. "

    @testset "Diskretisering" begin
        unit_range = (0.0, 1.0)

        @test_throws ArgumentError NRES.Tile_map{Dummy_type}( unit_range, N=missing )

        case_N1 = NRES.Tile_map{Dummy_type}( unit_range, N=1 )
        @test length(case_N1._all_SAT) == 1
        " With N=1, you get one SAT "

        case_N2 = NRES.Tile_map{Dummy_type}( unit_range, N=2 )
        @test length(case_N2._all_SAT) == 2
        " With N=2, you get two SAT "

        case_N3 = NRES.Tile_map{Dummy_type}( unit_range, N=3 )
        @test length(case_N3._all_SAT) == 3
        " With N=3, you get three SAT "

        case = NRES.Tile_map{Dummy_type}( unit_range, N=100 )
        @test length(case._all_SAT) == 100
        " With N=100, you get 100 SAT "

        for it ∈ case_N3._all_SAT
            isa(it, Dummy_type)
        end
        " All elements in vector is of type T "

        @test case_N1._number_of_intervals == 1
        @test case_N2._number_of_intervals == 2
        @test case_N3._number_of_intervals == 3
        @test case._number_of_intervals == 100
        " member variable that holds N "

        @test first(case_N1._all_SAT) == last(case_N1._all_SAT)
        @test first(case_N2._all_SAT) != last(case_N2._all_SAT)
        " N1: first = last. N2 first != last. "

        # TODO Gjør BoundsError bedre: gjør en test i map_to_SAT som sjekker bounds?
        @test NRES.map_to_SAT(case_N2, coordinate = first(unit_range)) == first(case_N2._all_SAT)
        @test NRES.map_to_SAT(case_N2, coordinate = last(unit_range))  == last(case_N2._all_SAT)
        " Ytterpunkta: I N2 gir first coordinate den første SAT, og last coordinate siste SAT "

        @test_throws ArgumentError NRES.map_to_SAT(case_N2, coordinate = first(unit_range)-0.001)
        @test_throws ArgumentError NRES.map_to_SAT(case_N2, coordinate = last(unit_range)+0.001)
        " Utafor ytterpunkta: Kastar ArgumentError "
    end
end

end # module TEST_TILE_MAP

max_stars_repo_name: leikanger/NRES

max_issues_repo_name: leikanger/NRES

max_forks_repo_name: leikanger/NRES