text:
"""
    skipmissingpairs(x::RealOrMissingVector, y::RealOrMissingVector)
Returns a pair `(a,b)`, filtered copies of `x` and `y`, in which elements `x[i]` and `y[i]`
are "skipped" (filtered out) if either `ismissing(x[i])` or `ismissing(y[i])`.
"""
function skipmissingpairs(x::RealOrMissingVector{T}, y::RealOrMissingVector{U}) where T where U

    length(x) == length(y) || error("Vectors must have same length")

    T2 = x isa Vector{Missing} ? Missing : T
    U2 = y isa Vector{Missing} ? Missing : U

    nout::Int = 0
    @inbounds for i = 1:length(x)
        if !(ismissing(x[i]) || ismissing(y[i]))
            nout += 1
        end
    end

    res1 = Vector{T2}(undef, nout)
    res2 = Vector{U2}(undef, nout)
    j::Int = 0

    @inbounds for i = 1:length(x)
        if !(ismissing(x[i]) || ismissing(y[i]))
            j += 1
            res1[j] = x[i]
            res2[j] = y[i]
        end
    end

    res1, res2
end

"""
    skipmissingpairs(X::RealOrMissingMatrix)
Returns `A`, a filtered copy of `X`, in which the row `X[i,:]` is "skipped" (filtered out)
if `any(ismissing,X[i,:])`.
"""
function skipmissingpairs(X::RealOrMissingMatrix{T}) where T
    T2 = X isa Matrix{Missing} ? Missing : T
    nr,nc = size(X)

    chooser = fill(true,nr)

    nrout = nr
    @inbounds for i = 1:nr
        for j = 1:nc
            if ismissing(X[i,j])
                chooser[i]=false
                nrout -= 1
                break
            end
        end
    end

    res = Matrix{T2}(undef,nrout,nc)
    @inbounds for j = 1:nc
        k = 0
        for i = 1:nr
            if chooser[i]
                k+=1
                res[k,j]=X[i,j]
            end
        end
    end
    res
end

"""
    skipmissingpairs(X::RealOrMissingMatrix,Y::RealOrMissingMatrix)
Returns a pair `(A,B)`, filtered copies of `X` and `Y`, in which the rows `X[i,:]` and
`Y[i,:]` are "skipped" (filtered out) if either `any(ismissing,X[i,:])`  or
`any(ismissing,Y[i,:])`.
"""
function skipmissingpairs(X::RealOrMissingMatrix{T},Y::RealOrMissingMatrix{U}) where T where U

    size(X,1)==size(Y,1) || error("arrays must have the same number of rows,"
                                * " but got $(size(X,1)) and $(size(Y,1))")

    T2 = X isa Matrix{Missing} ? Missing : T
    U2 = Y isa Matrix{Missing} ? Missing : U

    nr,ncx = size(X)
    ncy = size(Y,2)

    chooser = fill(true,nr)

    nrout = nr
    @inbounds for i = 1:nr
        for j = 1:ncx
            if ismissing(X[i,j])
                chooser[i]=false
                nrout -= 1
                break
            end
        end
        if chooser[i]
            for j = 1:ncy
                if ismissing(Y[i,j])
                    chooser[i]=false
                    nrout -= 1
                    break
                end
            end
        end
    end

    res1 = Matrix{T2}(undef,nrout,ncx)
    @inbounds for j = 1:ncx
        k = 0
        for i = 1:nr
            if chooser[i]
                k+=1
                res1[k,j]=X[i,j]
            end
        end
    end

    res2 = Matrix{U2}(undef,nrout,ncy)
    @inbounds for j = 1:ncy
        k = 0
        for i = 1:nr
            if chooser[i]
                k+=1
                res2[k,j]=Y[i,j]
            end
        end
    end

    res1, res2

end

"""
    skipmissingpairs(x::RealOrMissingVector,Y::RealOrMissingMatrix)
Returns a pair `(a,B)`, filtered copies of `x` and `Y`, in which the elements `x[i]` and
rows `Y[i,:]` are "skipped" (filtered out) if either `ismissing(x[i])` or
`any(ismissing,Y[i,:])`.
"""
function skipmissingpairs(x::RealOrMissingVector{T},Y::RealOrMissingMatrix{U}) where T where U
    length(x)==size(Y,1) || error("vector length must must match number of rows in matrix,"
                                * " but got $(length(x)) and $(size(Y,1))")

    T2 = x isa Vector{Missing} ? Missing : T
    U2 = Y isa Matrix{Missing} ? Missing : T
    nr = length(x)
    ncy = size(Y,2)

    chooser = fill(true,nr)

    nrout = nr
    @inbounds for i = 1:nr
        if ismissing(x[i])
            chooser[i]=false
            nrout -= 1
        end
        if chooser[i]
            for j = 1:ncy
                if ismissing(Y[i,j])
                    chooser[i]=false
                    nrout -= 1
                    break
                end
            end
        end
    end

    res1 = Vector{T2}(undef,nrout)
    k = 0
    @inbounds for i = 1:nr
        if chooser[i]
            k+=1
            res1[k]=x[i]
        end
    end

    res2 = Matrix{U2}(undef,nrout,ncy)
    @inbounds for j = 1:ncy
        k = 0
        for i = 1:nr
            if chooser[i]
                k+=1
                res2[k,j]=Y[i,j]
            end
        end
    end

    res1, res2

end

function skipmissingpairs(X::RealOrMissingMatrix,y::RealOrMissingVector)
    res2,res1 = skipmissingpairs(y,X)
    res1,res2
end


text:
using StemCellModels
using DataFrames
using ProgressMeter
using StatsBase
using Statistics
using Random
using StatsBase
using ArgParse
using CSV
using Distributions

s = ArgParseSettings()
@add_arg_table s begin
    "--resultsfile1"
        help = "File 1 to save results"
    "--resultsfile2"
        help = "File 2 to save results"
end

Random.seed!(123)

parsed_args = parse_args(ARGS, s)
print(parsed_args)
include("optim.jl")

function intervaldnds_data(Msel, Mneut, mmin, Nmax; mgap = 5.0, mup = 1.0, mud = 1.0)
    #calculate cumulative dn/ds across cohort, normalizing for mutation rates
    dndscumsum = Float64[]
    dnvec = Float64[]
    dsvec = Float64[]
    driverVAFt = filter(x -> x >= mmin, Msel)
    passengerVAFt = filter(x -> x >= mmin, Mneut)
    freqs = collect(mmin:mgap:Nmax)
    for f in freqs
        dn = length(filter(x -> x< f, driverVAFt))/mud
        ds = length(filter(x -> x< f, passengerVAFt))/mup
        push!(dndscumsum, dn/ds)
        push!(dnvec, dn)
        push!(dsvec, ds)
    end
    return dndscumsum[2:end], freqs[2:end], dnvec[2:end], dsvec[2:end]
end

#function to simulate population of cells with different paramters
function simulatepopulation(;Δ = 0.0, ρ = 100.0, Amin = 1 ./ ρ, gap = 0.02, tend = 20.0,
        Ncells = 1000, r = 0.5, λ = 1.0,
        sample = false, nsample = 10, Amax = 0.0)

    rangeN = 1:10^4
    Nits = 100
    Cntd = zeros(length(rangeN), Nits)
    Cntp = zeros(length(rangeN), Nits)
    Msel = Int64[]
    Mneut = Int64[]
    Δvec = Float64[]
    SMtog = StemCellModel(Ncells, Δmut = Δ, μp = 0.001, μd = 0.001, tend = tend, r = r, λ = λ)
    for i in 1:Nits
        Δsample = rand(Exponential(Δ))
        push!(Δvec, Δsample)
        SMtog = StemCellModel(Ncells, Δmut = Δsample, μp = 0.001, μd = 0.001, tend = tend, r = r, λ = λ)
        scst = runsimulation(SMtog, progress = false, finish = "time", onedriver = true, restart = true, maxcells = 10^6)
        append!(Msel, filter!(x -> x > 0, counts(sort(StemCellModels.cellsconvert(scst.stemcells)[2]), rangeN)))
        append!(Mneut, filter!(x -> x > 0, counts(sort(StemCellModels.cellsconvert(scst.stemcells)[1]), rangeN)))
    end

    Asel = Msel ./ ρ
    Aneut = Mneut ./ ρ
    Asel = filter(x -> x >= Amin, Asel)
    Aneut = filter(x -> x >= Amin, Aneut)

    if sample
        prevlength = length(Asel)
        Asel = Asel[StatsBase.sample(1:length(Asel), nsample, replace = false)]
        Aneut = Aneut[StatsBase.sample(1:length(Aneut),
                Int64(round(nsample * (length(Aneut) / prevlength))), replace = false)]
    end

    if Amax == 0.0
        Amax = maximum(Asel)
    end

    x1 = intervaldnds_data(Asel, Aneut, Amin, Amax; mgap = gap, mup = SMtog.μp, mud = SMtog.μd)

    DFinterval = DataFrame(dnds = x1[1], A = x1[2], dn = x1[3] * SMtog.μd / Nits,
    ds = x1[4] * SMtog.μd / Nits)

    return DFinterval, SMtog, [DFinterval[:dn][end] * Nits, DFinterval[:ds][end] * Nits], mean(Δvec)
end

println("Generating some example data with fits")

#create empty dataframe to store data
myDF = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Int64],
[:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda, :rsq, :deltatrue,
    :deltatruesample, :lambdartrue, :simindex], 0)

global j = 0
for i in 1:1000
    print(i)
    for Δ in [0.05, 0.1, 0.2]
        println(Δ)
        for rlam in [0.25]
            global j = j + 1
            println(rlam)
            rlambda = rlam
            DF1, SM, x1, Δsample = simulatepopulation(;Δ = Δ, tend = 30.0, Amin = 0.05, ρ = 100.0, λ = rlambda, Amax = 25.0)
            x = LLoptimizationresults(DF1[:dnds], DF1[:A]; t = 30.0, Amin = 0.05, ρ = 100.0)
            x.DF[:deltatrue] = Δ
            x.DF[:deltatruesample] = Δsample
            x.DF[:lambdartrue] = rlambda * SM.r
            x.DF[:simindex] = j
            append!(myDF, x.DF)
        end
    end
end

println("Writing data to file")
@rput myDF
R"""
res <- dplyr::distinct(myDF, deltafit, lambdarfit, deltatrue, deltatruesample, simindex, rsq)
readr::write_csv(res, $(parsed_args["resultsfile1"]))
""";

#function to simulate population of cells with different paramters
function simulatepopulation(;Δ = 0.0, ρ = 100.0, Amin = 1 ./ ρ, gap = 0.02, tend = 20.0,
        Ncells = 1000, r = 0.5, λ = 1.0,
        sample = false, nsample = 10, Amax = 0.0)

    rangeN = 1:10^4
    Nits = 100
    Cntd = zeros(length(rangeN), Nits)
    Cntp = zeros(length(rangeN), Nits)
    Msel = Int64[]
    Mneut = Int64[]
    Δvec = Float64[]
    SMtog = StemCellModel(Ncells, Δmut = Δ, μp = 0.001, μd = 0.001, tend = tend, r = r, λ = λ)
    for i in 1:Nits
        Δsample = rand(Beta(0.3, 0.3 * ((1 - Δ) / Δ)))
        Δsample = Δ
        push!(Δvec, Δsample)
        SMtog = StemCellModel(Ncells, Δmut = Δsample, μp = 0.001, μd = 0.001, tend = tend, r = r, λ = λ)
        scst = runsimulation(SMtog, progress = false, finish = "time", onedriver = true, restart = true, maxcells = 10^6)
        append!(Msel, filter!(x -> x > 0, counts(sort(StemCellModels.cellsconvert(scst.stemcells)[2]), rangeN)))
        append!(Mneut, filter!(x -> x > 0, counts(sort(StemCellModels.cellsconvert(scst.stemcells)[1]), rangeN)))
    end

    Asel = Msel ./ ρ
    Aneut = Mneut ./ ρ
    Asel = filter(x -> x >= Amin, Asel)
    Aneut = filter(x -> x >= Amin, Aneut)

    if sample
        prevlength = length(Asel)
        Asel = Asel[StatsBase.sample(1:length(Asel), nsample, replace = false)]
        Aneut = Aneut[StatsBase.sample(1:length(Aneut),
                Int64(round(nsample * (length(Aneut) / prevlength))), replace = false)]
    end

    if Amax == 0.0
        Amax = maximum(Asel)
    end

    x1 = intervaldnds_data(Asel, Aneut, Amin, Amax; mgap = gap, mup = SMtog.μp, mud = SMtog.μd)

    DFinterval = DataFrame(dnds = x1[1], A = x1[2], dn = x1[3] * SMtog.μd / Nits,
    ds = x1[4] * SMtog.μd / Nits)

    return DFinterval, SMtog, [DFinterval[:dn][end] * Nits, DFinterval[:ds][end] * Nits], mean(Δvec)
end

println("Generating some example data with fits")

#create empty dataframe to store data
myDF = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Int64],
[:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda, :rsq, :deltatrue,
    :deltatruesample, :lambdartrue, :simindex], 0)

global j = 0
for i in 1:1000
    print(i)
    for Δ in [0.05, 0.1, 0.2]
        println(Δ)
        for rlam in [0.25]
            global j = j + 1
            println(rlam)
            rlambda = rlam
            DF1, SM, x1, Δsample = simulatepopulation(;Δ = Δ, tend = 30.0, Amin = 0.05, ρ = 100.0, λ = rlambda, Amax = 25.0)
            x = LLoptimizationresults(DF1[:dnds], DF1[:A]; t = 30.0, Amin = 0.05, ρ = 100.0)
            x.DF[:deltatrue] = Δ
            x.DF[:deltatruesample] = Δsample
            x.DF[:lambdartrue] = rlambda * SM.r
            x.DF[:simindex] = j
            append!(myDF, x.DF)
        end
    end
end

println("Writing data to file")
@rput myDF
R"""
res <- dplyr::distinct(myDF, deltafit, lambdarfit, deltatrue, deltatruesample, simindex, rsq)
readr::write_csv(res, $(parsed_args["resultsfile2"]))
""";

ext: jl


text:
module NCCL

using Printf

using CuArrays

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream, CuStream, CuPtr

using CEnum

const libnccl = "libnccl"

# core library
include("libnccl_common.jl")
include("error.jl")
include("libnccl.jl")

include("base.jl")
include("communicator.jl")
include("group.jl")
include("collective.jl")

function __init__()
    precompiling = ccall(:jl_generating_output, Cint, ()) != 0
    silent = parse(Bool, get(ENV, "JULIA_CUDA_SILENT", "false")) || precompiling
    verbose = parse(Bool, get(ENV, "JULIA_CUDA_VERBOSE", "false"))

    # if any dependent GPU package failed, expect it to have logged an error and bail out
    if !CUDAdrv.functional()
        verbose && @warn "NCCL.jl did not initialize because CUDAdrv.jl failed to"
        return
    end

    try
        if version() < v"2.4"
            silent || @warn "NCCL.jl only supports NCCL 2.4 or higher (you are using $(version()))"
        end
    catch ex
        # don't actually fail to keep the package loadable
        if !silent
            if verbose
                @error "NCCL.jl failed to initialize" exception=(ex, catch_backtrace())
            else
                @info "NCCL.jl failed to initialize, GPU functionality unavailable (set JULIA_CUDA_SILENT or JULIA_CUDA_VERBOSE to silence or expand this message)"
            end
        end
    end
end

end


text:
cols_1_to_30 = [
    -2.5   0.5         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   5.26957  -2.45701   -0.812567   0.0        0.0         0.0         0.0        0.0         0.0         0.0
    1.0  -3.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   5.26957  -2.45701   -0.812567   0.0        0.0         0.0         0.0        0.0         0.0         0.0
    0.0   0.0        -2.5   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   2.0       0.0        0.0        0.0        0.0         0.0         0.0        0.0         0.0         0.728864
    0.0   0.0         0.0  -2.5  -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       2.0        0.0        0.0        0.0         0.0         0.0        0.0         0.0         1.05443
    0.0   0.0         0.0   0.0  -2.5   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        2.0        0.0        0.0         0.0         0.0        0.0         0.0         0.564472
    0.0   0.0         0.0   0.0   0.0  -2.5   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        2.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  -2.5  -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        2.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  -2.5   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         2.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  -2.5   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         2.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  -2.5  -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        2.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  -2.5   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         2.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  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.91546  -1.1331     0.217645   0.0        0.0         0.0         0.0        0.0         0.0         0.5
    0.0   0.0         0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   4.21772  -2.19645   -0.021269   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  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   2.25789  -0.845166   0.587279   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  -2.5   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.91546   -1.1331      0.217645    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  -2.5  -0.0   0.0   0.0   0.0   0.0       0.0        0.0        4.21772   -2.19645    -0.021269    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  -2.5   0.0   0.0   0.0   0.0       0.0        0.0        2.25789   -0.845166    0.587279    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  -2.5   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.91546   -1.1331      0.217645    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  -2.5  -0.0   0.0       0.0        0.0        0.0        0.0         0.0         4.21772   -2.19645    -0.021269    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  -2.5   0.0       0.0        0.0        0.0        0.0         0.0         2.25789   -0.845166    0.587279    0.0
    0.0   0.0         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  -3.12854   0.0        0.0        2.97348   -1.38642    -0.458509    0.0        0.0         0.0         1.45773
    0.0   0.0         0.0   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      -3.12854   -1.16909    6.05379   -2.82266    -0.933492    0.0        0.0         0.0         2.10886
    0.0   0.0         0.0   0.0   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       1.16909   -3.12854   -0.106837   0.0498141   0.0164742   0.0        0.0         0.0         1.12894
    0.0   0.0         0.0   0.0   0.0   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       -3.12854    0.0         0.0         2.97348   -1.38642    -0.458509    0.0
    0.0   0.0         0.0   0.0   0.0   0.0   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       -3.12854    -1.16909     6.05379   -2.82266    -0.933492    0.0
    0.0   0.0         0.0   0.0   0.0   0.0   0.0   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        1.16909    -3.12854    -0.106837   0.0498141   0.0164742   0.0
    0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   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        -3.12854    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.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       -3.12854    -1.16909     0.0
    0.0   0.0         0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   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        1.16909    -3.12854     0.0
    0.0   1.69282     0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.91546  -1.1331     0.217645   0.0        0.0         0.0         0.0        0.0         0.0        -4.69282
    0.0   3.44646     0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   4.21772  -2.19645   -0.021269   0.0        0.0         0.0         0.0        0.0         0.0         0.0
    0.0  -0.0608229   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   2.25789  -0.845166   0.587279   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.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.91546   -1.1331      0.217645    0.0        0.0         0.0         4.46021
    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.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        4.21772   -2.19645    -0.021269    0.0        0.0         0.0         9.08068
    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.0   0.0   0.0   0.0   0.0       0.0        0.0        2.25789   -0.845166    0.587279    0.0        0.0         0.0        -0.160255
    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.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.91546   -1.1331      0.217645    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.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         4.21772   -2.19645    -0.021269    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.0   0.0       0.0        0.0        0.0        0.0         0.0         2.25789   -0.845166    0.587279    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.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
]

cols_31_to_40 = [
    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.283275    0.0544112    0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
   -0.549114   -0.00531726   0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
   -0.211292    0.14682      0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
    0.0         0.0          0.728864  -0.283275    0.0544112    0.0        0.0        0.0          0.0   0.0
    0.0         0.0          1.05443   -0.549114   -0.00531726   0.0        0.0        0.0          0.0   0.0
    0.0         0.0          0.564472  -0.211292    0.14682      0.0        0.0        0.0          0.0   0.0
    0.0         0.0          0.0        0.0         0.0          0.728864  -0.283275   0.0544112    0.0   0.0
    0.0         0.0          0.0        0.0         0.0          1.05443   -0.549114  -0.00531726   0.0   0.0
    0.0         0.0          0.0        0.0         0.0          0.564472  -0.211292   0.14682      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.5         0.0          0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
    0.0         0.5          0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
    0.0         0.0          0.5        0.0         0.0          0.0        0.0        0.0          0.0   0.0
    0.0         0.0          0.0        0.5         0.0          0.0        0.0        0.0          0.0   0.0
    0.0         0.0          0.0        0.0         0.5          0.0        0.0        0.0          0.0   0.0
    0.0         0.0          0.0        0.0         0.0          0.5        0.0        0.0          0.0   0.0
    0.0         0.0          0.0        0.0         0.0          0.0        0.5        0.0          0.0   0.0
    0.0         0.0          0.0        0.0         0.0          0.0        0.0        0.5          0.0   0.0
   -0.56655     0.108822     0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
   -1.09823    -0.0106345    0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
   -0.422583    0.29364      0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
    0.0         0.0          1.45773   -0.56655     0.108822     0.0        0.0        0.0          0.0   0.0
    0.0         0.0          2.10886   -1.09823    -0.0106345    0.0        0.0        0.0          0.0   0.0
    0.0         0.0          1.12894   -0.422583    0.29364      0.0        0.0        0.0          0.0   0.0
    0.0         0.0          0.0        0.0         0.0          1.45773   -0.56655    0.108822     0.0   1.12854
    0.0         0.0          0.0        0.0         0.0          2.10886   -1.09823   -0.0106345    0.0   2.29764
    0.0         0.0          0.0        0.0         0.0          1.12894   -0.422583   0.29364      0.0  -0.0405486
    0.0         0.0          0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
   -4.69282    -1.75364      0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
    1.75364    -4.69282      0.0        0.0         0.0          0.0        0.0        0.0          0.0   0.0
   -2.07963    -0.687763    -4.69282    0.0         0.0          0.0        0.0        0.0          0.0   0.0
   -4.23399    -1.40024      0.0       -4.69282    -1.75364      0.0        0.0        0.0          0.0   0.0
    0.0747212   0.0247113    0.0        1.75364    -4.69282      0.0        0.0        0.0          0.0   0.0
    0.0         0.0          4.46021   -2.07963    -0.687763    -4.69282    0.0        0.0          0.0   0.0
    0.0         0.0          9.08068   -4.23399    -1.40024      0.0       -4.69282   -1.75364      0.0   0.0
    0.0         0.0         -0.160255   0.0747212   0.0247113    0.0        1.75364   -4.69282      0.0   0.0
    0.0         0.0          0.0        0.0         0.0          1.31739   -0.614252  -0.203142    -2.5   2.0
    0.0         0.0          0.0        0.0         0.0          2.63479   -1.2285    -0.406283     1.0  -2.0
]

B_data = [cols_1_to_30 cols_31_to_40]


text:
module TestModelSearch

using Test
using MLJModels
import MLJBase
using ScientificTypes

pca = info("PCA", pkg="MultivariateStats")
cnst = info("ConstantRegressor", pkg="MLJModels")
tree = info("DecisionTreeRegressor", pkg="DecisionTree")

@test_throws ArgumentError info("Julia")

@test info(ConstantRegressor) == cnst
@test info(Standardizer()) == info("Standardizer", pkg="MLJModels")

@testset "localmodels" begin
    @test cnst in localmodels(modl=TestModelSearch)
    @test !(tree in localmodels(modl=TestModelSearch))
    import MLJDecisionTreeInterface.DecisionTreeRegressor
    @test tree in localmodels(modl=TestModelSearch)
end

@testset "matching" begin
    X = ones(2,3)
    y = ones(2)
    w = ones(2)
    w_class = Dict(:male => 1.0, :female => 10.0)

    task = matching(X)
    @test !task.is_supervised
    @test ismissing(task.supports_weights)
    @test ismissing(task.supports_class_weights)
    @test task.input_scitype == scitype(X)
    @test ismissing(task.target_scitype)

    task = matching(X, y)
    @test task.is_supervised
    @test ismissing(task.supports_weights)
    @test ismissing(task.supports_class_weights)
    @test task.input_scitype == scitype(X)
    @test task.target_scitype == scitype(y)

    task = matching(X, y, w)
    @test task.is_supervised
    @test task.supports_weights
    @test !task.supports_class_weights
    @test task.input_scitype == scitype(X)
    @test task.target_scitype == scitype(y)

    task = matching(X, y, w_class)
    @test task.is_supervised
    @test !task.supports_weights
    @test task.supports_class_weights
    @test task.input_scitype == scitype(X)
    @test task.target_scitype == scitype(y)
end

@testset "models() and localmodels()" begin
    t(model) = model.is_pure_julia
    mods = models(t)
    @test pca in mods
    @test cnst in mods
    @test !(info("SVC") in mods)
    mods = localmodels(t, modl=TestModelSearch)
    @test cnst in mods
    @test !(pca in mods)
    u(model) = !(model.is_supervised)
    @test pca in models(u, t)
    @test !(cnst in models(u, t))
end

@testset "models(matching())" begin
    X = MLJBase.table(ones(2,3))
    y = MLJBase.coerce(["a", "b"], Multiclass)
    w = ones(2)
    ms = models(matching(X, y, w))

    # by hand:
    ms2 = models() do m
        m.is_supervised &&
            AbstractVector{Multiclass{2}} <: m.target_scitype &&
            Table(Continuous) <: m.input_scitype &&
            m.supports_weights
    end
    @test ms == ms2
end

@testset "models(needle)" begin
    @test pca in models("PCA")
    @test pca ∉ models("PCA′")
    @test pca in models(r"PCA")
    @test pca ∉ models(r"PCA′")
end

end
true


text:
export dotrap, dodebugtrap

@inline dotrap() = ccall("llvm.trap", llvmcall, Cvoid, ())
@inline dodebugtrap() = ccall("llvm.debugtrap", llvmcall, Cvoid, ())


text:
function energie_kgr(u, v, epsilon, k, dx, llambda, sigma)

    z = (u + conj(v)) / 2
    dtz = 1im / 2 / epsilon * ifft(sp.sqrt(1 + epsilon * k .^ 2) * fft(u - conj(v)))

    Q = -epsilon * dx * sp.real(1im * sp.sum(dtz * conj(z)))
    E = epsilon * dx * np.linalg.norm(dtz) .^ 2 \
        + dx * np.linalg.norm(ifft(k * fft(z))) .^ 2 \
        + dx * np.linalg.norm(
        z) .^ 2 / epsilon \
        - llambda / (sigma + 1) * dx * sp.sum(abs(z) .^ (2 * sigma + 2))

    return Q, E

end


text:
"""
    AbstractBinnedResampling

Resampling schemes where data are binned.
"""
abstract type AbstractBinnedResampling <: AbstractUncertainDataResampling end

"""
    AbstractBinnedUncertainValueResampling <: AbstractUncertainDataResampling

Resampling schemes where data are binned and the data in each bin
are represented as an uncertain value.
"""
abstract type AbstractBinnedUncertainValueResampling <: AbstractBinnedResampling end

"""
    AbstractBinnedSummarisedResampling <: AbstractUncertainDataResampling

Resampling schemes where data are binned and the data in each bin
are summarised to a single value (e.g. the mean).
"""
abstract type AbstractBinnedSummarisedResampling <: AbstractBinnedResampling end

export
AbstractBinnedResampling,
AbstractBinnedSummarisedResampling,
AbstractBinnedUncertainValueResampling


text:
# indicator of a halfspace

export IndHalfspace

"""
**Indicator of a halfspace**

    IndHalfspace(a, b)

For an array `a` and a scalar `b`, returns the indicator of set
```math
S = \\{x : \\langle a,x \\rangle \\leq b \\}.
```
"""
struct IndHalfspace{R <: Real, T <: AbstractArray{R}} <: ProximableFunction
  a::T
  b::R
  norm_a::R
  function IndHalfspace{R, T}(a::T, b::R) where {R <: Real, T <: AbstractArray{R}}
    norm_a = norm(a)
    if norm_a == 0 && b < 0
        error("function is improper")
    end
    new(a, b, norm_a)
  end
end

IndHalfspace(a::T, b::R) where {R <: Real, T <: AbstractArray{R}} = IndHalfspace{R, T}(a, b)

is_convex(f::IndHalfspace) = true
is_set(f::IndHalfspace) = true
is_cone(f::IndHalfspace) = f.b == 0 || f.b == Inf

function (f::IndHalfspace{R})(x::AbstractArray{R}) where R
  if dot(f.a, x) - f.b <= eps(R)*f.norm_a*(1 + abs(f.b))
    return zero(R)
  end
  return R(Inf)
end

function prox!(y::AbstractArray{R}, f::IndHalfspace{R}, x::AbstractArray{R}, gamma::R=one(R)) where R
  s = dot(f.a, x)
  if s > f.b
    y .= x .- ((s - f.b)/f.norm_a^2) .* f.a
  else
    copyto!(y, x)
  end
  return zero(R)
end

fun_name(f::IndHalfspace) = "indicator of a halfspace"
fun_dom(f::IndHalfspace) = "AbstractArray{Real}"
fun_expr(f::IndHalfspace) = "x ↦ 0 if <a,x> ⩽ b, +∞ otherwise"
fun_params(f::IndHalfspace) =
  string( "a = ", typeof(f.a), " of size ", size(f.a), ", ",
          "b = $(f.b)")

function prox_naive(f::IndHalfspace{R}, x::AbstractArray{R}, gamma::R=one(R)) where R
  s = dot(f.a, x) - f.b
  if s <= 0
    return x, 0.0
  end
  return x - (s/norm(f.a)^2)*f.a, 0.0
end


text:
using ControlSystems, OrdinaryDiffEq

"""
    GainSchedulingSimulator
# Fields:
    s::Simulator
    f = (x,p,t) -> x
    y = (x,t) -> y
    r = (x,t) -> reference
    e = (x,t) -> r(x,t) .- x # Error
    controllers::Vector{StateSpace}
    conditions::Vector{Function}
"""
struct GainSchedulingSimulator <: ControlSystems.AbstractSimulator
    s::Simulator
    f
    y
    r
    e
    controllers::Vector{StateSpace}
    conditions::Vector{T} where T <: Function
    feedforward::Vector{Float64}
end
"""
    GainSchedulingSimulator(P,r,controllers::AbstractVector{LTISystem},conditions::AbstractVector{Function (x,y,r)->Bool}; inputfun=(u,t)->u, feedforward=zeros())
`r` is a function (x,t)-> reference signal
`controllers` is a list of LTIsystems representing the controllers
`conditions` is a list of function that accepts state, output and reference and returns a Bool indicating which controller is active. At least one function should return true for every input. If more than one condition is true, the index of the first condition in the list is chosen.
´feedforward´ is a vector of constants that are added to the control signal calculated by the active controller, i.e., `u = C*x + D*e + feedforward[active_index]`
For a usage example, see notebook in example folder.
```
"""
function GainSchedulingSimulator(P,ri,controllers::AbstractVector{Tu},conditions::AbstractVector{Tc};
    inputfun=(u,t)->u,
    feedforward=zeros(length(controllers))) where Tu <: StateSpace where Tc <: Function
    s = Simulator(P)
    pinds = 1:P.nx # Indices of plant-state derivative
    r = (x,t) -> ri(x[pinds],t)
    y(x,t) = s.y(x[pinds],t)
    y(sol::ODESolution,t) = P.C*sol(t)[pinds,:]
    e = (x,t) -> r(x,t) .- y(x,t)
    f = function(der,x,p,t)
        xyr = (x[pinds],y(x,t),r(x,t))
        index = findfirst(c->c(xyr...), conditions)
        @assert (index != nothing) "No condition returned true"
        et = e(x,t)
        ind = P.nx+1 # First index of currently updated controller
        for c in controllers
            c.nx == 0 && continue # Gain only, no states
            inds = ind:(ind+c.nx-1) # Controller indices
            der[inds] = c.A*x[inds] + c.B*et # Controller dynamics, driven by error
            ind += c.nx
        end
        c = controllers[index] # Active controller
        cind = P.nx + (index == 1 ? 1 : 1+sum(i->controllers[i].nx, 1:index-1)) # Get index of first controller state
        der[pinds] = P.A*x[pinds] # System dynamics
        u = c.C*x[cind:(cind+c.nx-1)] .+ c.D*et .+ feedforward[index] # Form control signal
        der[pinds] .+= vec(inputfun(P.B*u,t)) # Add input from active controller to system dynamics
        der
    end
    GainSchedulingSimulator(s,f,y,r,e,controllers,conditions,feedforward)
end

function GainSchedulingSimulator(P,ri,controllers::AbstractVector{Tu},conditions::AbstractVector{Tc};
    inputfun=(u,t)->u,
    feedforward=zeros(length(controllers))) where Tu <: LTISystem where Tc <: Function
    GainSchedulingSimulator(P,ri,ss.(controllers),conditions; inputfun=inputfun, feedforward=feedforward)
end
# ============================================================================================

"""
    sol = solve(s::AbstractSimulator, x0, tspan,  args...; kwargs...)
Simulate the system represented by `s` from initial state `x0` over time span `tspan = (t0,tf)`.
`args` and `kwargs` are sent to the `solve` function from `OrdinaryDiffEq`, e.g., `solve(s, x0, tspan,  Tsit5(), reltol=1e-5)` solves the problem with solver [`Tsit5()`](http://docs.juliadiffeq.org/stable/solvers/ode_solve.html) and relative tolerance 1e-5.

See also `Simulator` `lsim`
"""

DiffEqBase.solve(s::GainSchedulingSimulator, x0, tspan, solver=Tsit5(), args...; kwargs...) = solve(ODEProblem(s.f,vcat(x0, zeros(sum(c->c.nx, s.controllers))),tspan), solver, args...; kwargs...)


text:
include("basic_dqn.jl")
include("dqn.jl")
include("prioritized_dqn.jl")
include("qr_dqn.jl")
include("rem_dqn.jl")
include("rainbow.jl")
include("iqn.jl")
include("common.jl")


text:
include("Parser.jl")
include("Lexer.jl")
include("Token.jl")
include("Transpile.jl")

max_stars_repo_name: GabrielTrettel/Portran

max_issues_repo_name: GabrielTrettel/Portran

max_forks_repo_name: GabrielTrettel/Portran


text:
include("./main.jl")
module macros
using macroexpansion

max_stars_repo_name: timo-cmd/DANCE

max_issues_repo_name: timo-cmd/Twister

max_forks_repo_name: timo-cmd/DANCE


text:
# using SuperResolution
using Distributions
include("fctSR.jl")
include("modelSR.jl")

mutable struct SimSR
    sample::Array{Array{Float64,2},1} #each array entry corresponds to a sample
    x::Array{Float64,2} #gathered x values of all simulations (for plot purpose)
    y::Array{Float64,2} #gathered y values
    function SimSR(model::MixtureModel,simN,simIT)  #simN..number of samples drawn // #simIT..repetitions of sampling
         obj=new()
         obj.sample=Array{Array{Float64,2},1}(undef,simIT)
         obj.x=Array{Float64,2}(undef,simN,simIT)
         obj.y=Array{Float64,2}(undef,simN,simIT)
         for i=1:simIT obj.sample[i]=ctrans(rand(model,simN)) end
         obj.x[:,1]=obj.sample[1][:,1]
         obj.y[:,1]=obj.sample[1][:,2]
         map(i->obj.x[:,1:i]=hcat(obj.x[:,1:(i-1)],obj.sample[i][:,1]),2:length(obj.sample))
         map(i->obj.y[:,1:i]=hcat(obj.y[:,1:(i-1)],obj.sample[i][:,2]),2:length(obj.sample))
         obj
    end
end

#   !!for a selected model K_SR!! (i.e., there are assumably K_SR clusters in this model),
#   compute the inner distances and filter members for each clus_ID
mutable struct SegDistanceSR
    K_SR::Int64
    N_SR::Int64
    val_bool::BitArray{1}
    val_bool2::Array{BitArray{1},1}
    singleDist_SR::Array{Array{Float64,2},1}
    #singleDist_SR::Array{Float64,2}
    loc_SR::Array{Array{Array,1},1}
    frame_SR::Array{Array{Int64,1},1}
    cov_SR::Array{Array{Array{Float64,2},1},1}
    D_SR::Array{Float64,1}
    W_K::Float64
    weight_SR::Array{Array{Float64,2},1}
    function SegDistanceSR(modelSR::SuperResolution.ModelSR,K_SR::Int64,modINPUT::Int64)
        obj=new()
        obj.K_SR=K_SR
        obj.val_bool=length.(unique.(modelSR.member)).==K_SR
        obj.N_SR=length(modelSR.member[obj.val_bool])
        obj.D_SR=zeros(K_SR)

        obj.val_bool2=map(i->modelSR.member[obj.val_bool][1].==modelSR.clus_ID[obj.val_bool][1][i],1:length(modelSR.clus_ID[obj.val_bool][1])) #Array{BitArray{1},1}
        obj.loc_SR=map(i->dim2adj(modelSR.loc[1][obj.val_bool2[i],:]),1:length(obj.val_bool2))   #dim2adj(Array{Float64,2} of x,y loc) ->  Array{Array,1}
        obj.frame_SR=map(i->modelSR.frame[obj.val_bool2[i]],1:length(obj.val_bool2))
        obj.cov_SR=map(i->modelSR.cov[1][obj.val_bool2[i]],1:length(obj.val_bool2))
        #obj.weight_SR=ones(K_SR)#map(i->modelSR.weight[1][obj.val_bool2[i],:],1:length(obj.val_bool2))
        # dim2KLDfull_UNWEIGHTED(MU::Array{Array,1},S::Array{Array{Float64,2},1},modusKLD=2,weightParaKLD=NaN)
        if modINPUT==1  # weighted KL clustering
            obj.singleDist_SR=map(i->modelSR.wkld[1][obj.val_bool2[i],:],1:length(obj.val_bool2))
            obj.weight_SR=map(i->modelSR.wkld[1][obj.val_bool2[i],:],1:length(obj.val_bool2))
            obj.D_SR=map(i->sum(obj.singleDist_SR[i]),1:K_SR)
            obj.W_K=sum(i->sum(obj.D_SR[i]/(2*sum(obj.val_bool2[i]))),1:obj.K_SR)
            obj
        # if modINPUT==1  # weighted KL clustering
        #     obj.singleDist_SR=map(i->dim2KLDfull_WEIGHTED(obj.loc_SR[i],obj.cov_SR[i],obj.frame_SR[i])[1],1:K_SR)
        #     obj.weight_SR=map(i->dim2KLDfull_WEIGHTED(obj.loc_SR[i],obj.cov_SR[i],obj.frame_SR[i])[2],1:K_SR)
        #     obj.D_SR=map(i->sum(obj.singleDist_SR[i]),1:K_SR)
        #     obj.W_K=sum(i->sum(obj.D_SR[i]/(2*sum(obj.val_bool2[i]))),1:obj.K_SR)
        #     obj
        elseif modINPUT==2  # UNweighted KL clustering
        #dim2KLDfull_WEIGHTED(MU::Array{Array,1},S::Array{Array{Float64,2},1},SEG_FrNm)
            # obj.singleDist_SR=map(i->dim2KLDfull_UNWEIGHTED(obj.loc_SR[i],obj.cov_SR[i]),1:K_SR)
            # obj.D_SR=map(i->sum(obj.singleDist_SR[i]),1:K_SR)
            # obj.W_K=sum(i->sum(obj.D_SR[i]/(2*sum(obj.val_bool2[i]))),1:obj.K_SR)
            # obj
            obj.singleDist_SR=map(i->modelSR.wkld[1][obj.val_bool2[i],:],1:length(obj.val_bool2))
            obj.weight_SR=map(i->modelSR.wkld[1][obj.val_bool2[i],:],1:length(obj.val_bool2))
            obj.D_SR=map(i->sum(obj.singleDist_SR[i]),1:K_SR)
            obj.W_K=sum(i->sum(obj.D_SR[i]/(2*sum(obj.val_bool2[i]))),1:obj.K_SR)
            obj
        end
    end
end
mutable struct DistanceSR
    Distance::Array{SuperResolution.SegDistanceSR,1}
    function DistanceSR(model_SR::SuperResolution.ModelSR,modINPUT::Int64)
        obj=new()
        obj.Distance=Array{SuperResolution.SegDistanceSR,1}(undef,model_SR.K[1])
        if modINPUT==1  # weighted KL clustering
            print("##################################################################\n")
            print("Computing full DistanceSR struct WEIGHTED KLD (modINPUT==1)\n")
            @time(
                for j=1:model_SR.K[1]
                    obj.Distance[j]=SegDistanceSR(model_SR,j,modINPUT)
                end
            )
        elseif modINPUT==2  # UNweighted KL clustering
            print("##################################################################\n")
            print("Computing full DistanceSR struct UNweighted KLD (modINPUT==2)\n")
            @time(
                for j=1:model_SR.K[1]
                    obj.Distance[j]=SegDistanceSR(model_SR,j,modINPUT)
                end
            )
        end
        obj
    end
end
## Evtl. ein AIC Verfahren zur Modellbestimmung implementieren
mutable struct clusWKLD
    modus::Int64
    wModel::SuperResolution.ModelSR
    wkld::Array{Array{Float64,2},1}
    weight::Array{Array{Float64,2},1}
    member::Array{Array{Int64,1},1}
    wkldFilteredCEIL::Array{BitArray{2},1}
    wkldFilteredFLOOR::Array{BitArray{2},1}

    wkldMean::Array{Float64,1}
    wkldMedian::Array{Float64,1}
    wkldMax::Array{Float64,1}
    wkldMin::Array{Float64,1}
    wkldMinDelta::Array{Float64,1}

    max_mean_idx::Array{Int64,1}
    max_median_idx::Array{Int64,1}
    max_max_idx::Array{Int64,1}
    max_minDelta_idx::Array{Int64,1}

    KL_Distance::SuperResolution.DistanceSR
    W_K::Array{Float64,1}
    W_K_star::Array{Float64,1}
    diff_gap::Array{Float64,1}
    diff_gap_max_idx::Array{Int64,1}
    diff_WK::Array{Float64,1}
    diff_WK_ratio::Array{Float64,1}
    diff_WK_max_idx::Array{Int64,1}
    Gap_n::Array{Float64,1}
    Gap::Array{Float64,1}
    s_k::Array{Float64,1}
    sd_k::Array{Float64,1}
    k_hat::Int64
    function clusWKLD(ModelWKLD_SR::SuperResolution.ModelSR,modINPUT=1)
        obj=new()
        obj.modus=modINPUT  # 1..original method, 2..AIC ?

        obj.wkld=ModelWKLD_SR.wkld
        obj.weight=ModelWKLD_SR.weight
        obj.member=ModelWKLD_SR.member

        tuP=0
        obj.wkldFilteredCEIL=map(i->ModelWKLD_SR.wkld[i].<10^4,2:(size(ModelWKLD_SR.wkld,1)-tuP))   # remove artif added big numbers
        obj.wkldFilteredFLOOR=map(i->ModelWKLD_SR.wkld[i].>1e-10,2:(size(ModelWKLD_SR.wkld,1)-tuP)) # remove artif added small numbers

        obj.wkldMean=map(i->mean(ModelWKLD_SR.wkld[i][obj.wkldFilteredCEIL[i-1]]),2:(size(ModelWKLD_SR.wkld,1)-tuP))
        obj.wkldMedian=map(i->median(ModelWKLD_SR.wkld[i][obj.wkldFilteredCEIL[i-1]]),2:(size(ModelWKLD_SR.wkld,1)-tuP))
        obj.wkldMax=map(i->maximum(ModelWKLD_SR.wkld[i][obj.wkldFilteredCEIL[i-1]]),2:(size(ModelWKLD_SR.wkld,1)-tuP))
        obj.wkldMin=map(i->minimum(ModelWKLD_SR.wkld[i][obj.wkldFilteredFLOOR[i-1]]),2:(size(ModelWKLD_SR.wkld,1)-tuP))

        obj.wkldMinDelta=abs.((obj.wkldMin[1:(end-1)]-obj.wkldMin[2:end]))#./obj.wkldMedian[2:(end)])

        obj.max_mean_idx        =   findall(i->i==maximum(obj.wkldMean),obj.wkldMean)
        obj.max_median_idx      =   findall(i->i==maximum(obj.wkldMedian),obj.wkldMedian)
        obj.max_max_idx         =   findall(i->i==maximum(obj.wkldMax),obj.wkldMax)
        obj.max_minDelta_idx    =   findall(i->i==maximum(obj.wkldMinDelta),obj.wkldMinDelta)
        obj
    end
end
export clusGAP
export SimSR
export SegDistanceSR
export DistanceSR
export clusWKLD

#unique(tmp5.Distance[end-1].singleDist_SR)

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
using Test
using Yao.EasyBuild
using Yao.EasyBuild: pair_ring, pair_square, merged_rotor, rotor, rotorset

@testset "pairs geometries" begin
    @test pair_ring(3) == [1=>2,2=>3,3=>1]
    ps = pair_square(2, 2; periodic=false)
    @test length(ps) == 4
    for item in [1=>2, 3=>4, 1=>3, 2=>4]
        @test item in ps
    end

    ps = pair_square(2, 2; periodic=true)
    @test length(ps) == 8
    for item in [1=>2, 3=>4, 2=>1, 4=>3, 1=>3, 2=>4, 3=>1, 4=>2]
        @test item in ps
    end
end

@testset "rotter, collect_blocks, num_gradient, opgrad" begin
    @test  merged_rotor(true, true) == Rx(0.0)
    @test  merged_rotor(false, false) == merged_rotor() == chain(Rz(0.0), Rx(0.0), Rz(0.0))
    @test  merged_rotor(false, true) == chain(Rz(0.0), Rx(0.0))
    @test  merged_rotor(true, false) == chain(Rx(0.0), Rz(0.0))
    @test collect_blocks(RotationGate, rotorset(:Merged, 5, true, false)) |> length == 10

    @test rotor(5, 2, true, true) isa ChainBlock
    @test rotor(5, 2, true, true) |> length == 1
    @test rotor(5, 2, true, true) |> nqubits == 5
    @test collect_blocks(PutBlock{2, <:Any, <:RotationGate}, rotorset(:Split, 5, true, false)) |> length == 10
end

@testset "entangler" begin
    c = variational_circuit(5; entangler=(n,i,j)->put(n,(i,j)=>ConstGate.CZ))
    @test nparameters(c) == 50
    c = variational_circuit(5; entangler=(n,i,j)->put(n,(i,j)=>rot(ConstGate.CNOT, 0.0)))
    @test nparameters(c) == 65
    ps = randn(nparameters(c))
    dispatch!(c, ps)
    @test parameters(c) ≈ ps
end


text:
using Plots
using Random
x,y = randn(100),randn(100)
f = scatter(x, y)
display(f)


text:
using DiffEqFlux, Flux, OrdinaryDiffEq, Test
#using Plots

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = (α - β*y)x
  du[2] = dy = (δ*x - γ)y
end
prob = ODEProblem(lotka_volterra,[1.0,1.0],(0.0,10.0))
const len = length(range(0.0,stop=10.0,step=0.1)) # 101

# Reverse-mode

p = param([2.2, 1.0, 2.0, 0.4])
params = Flux.Params([p])
function predict_rd()
  vec(diffeq_rd(p,prob,Tsit5(),saveat=0.1))
end
loss_rd() = sum(abs2,x-1 for x in predict_rd())
loss_rd()

grads = Tracker.gradient(loss_rd, params, nest=true)
grads[p]

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function ()
  display(loss_rd())
  #display(plot(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6)))
end

# Display the ODE with the current parameter values.
cb()

Flux.train!(loss_rd, params, data, opt, cb = cb)

# Forward-mode, R^n -> R^m layer

p = param([2.2, 1.0, 2.0, 0.4])
params = Flux.Params([p])
function predict_fd()
  diffeq_fd(p,vec,2*len,prob,Tsit5(),saveat=0.1) # 2 times for 2 output variables
end
loss_fd() = sum(abs2,x-1 for x in predict_fd())
loss_fd()

@test_broken begin
  grads = Tracker.gradient(loss_fd, params, nest=true)
  grads[p]
end

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function ()
  display(loss_fd())
  #display(plot(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6)))
end

# Display the ODE with the current parameter values.
cb()

Flux.train!(loss_fd, params, data, opt, cb = cb)

# Forward-mode, R^n -> R loss

p = param([2.2, 1.0, 2.0, 0.4])
params = Flux.Params([p])
loss_reduction(sol) = sum(abs2,x-1 for x in vec(sol))
function predict_fd2()
  diffeq_fd(p,loss_reduction,nothing,prob,Tsit5(),saveat=0.1)
end
loss_fd2() = predict_fd2()
loss_fd2()

grads = Tracker.gradient(loss_fd2, params, nest=true)
grads[p]

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function ()
  display(loss_fd2())
  #display(plot(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6)))
end

# Display the ODE with the current parameter values.
cb()

Flux.train!(loss_fd2, params, data, opt, cb = cb)

# Adjoint sensitivity
p = param([2.2, 1.0, 2.0, 0.4])
params = Flux.Params([p])
function predict_adjoint()
    diffeq_adjoint(p,prob,Tsit5())
end
loss_adjoint() = loss_reduction(predict_adjoint())
loss_adjoint()

grads = Tracker.gradient(loss_adjoint, params, nest=true)
grads[p]

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function ()
  display(loss_adjoint())
  #display(plot(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6)))
end

# Display the ODE with the current parameter values.
cb()

Flux.train!(loss_adjoint, params, data, opt, cb = cb)


text:
function adjust_rect(base_rect, adjusting_rect)

    base_x = base_rect[2].x
    base_y = base_rect[2].y
    base_w = base_rect[3].x - base_x
    base_h = base_rect[1].y - base_y

    adjusting_x = adjusting_rect[2].x
    adjusting_y = adjusting_rect[2].y
    adjusting_w = adjusting_rect[3].x - adjusting_x
    adjusting_h = adjusting_rect[1].y - adjusting_y

    rectangle = rect(Point(base_x + base_w * adjusting_x, base_y + base_h * adjusting_y), base_w * adjusting_w, base_h * adjusting_h, vertices = true)
    return rectangle
end


text:
struct DensityAnalysis
    options::Dict{Symbol, Any}
end
DensityAnalysis(; kwargs...) = DensityAnalysis(Dict{Symbol, Any}(kwargs))

# Work around lack of length 1 tuple method
# TODO: also add weights support here
_kde(data::NTuple{1, Any}; kwargs...) = kde(data...; kwargs...)
_kde(data::Tuple; kwargs...) = kde(data; kwargs...)

applydatalimits(f::Function, d) = map(f, d)
applydatalimits(i::Union{AbstractArray, Tuple}, _) = i

function _density(data...; datalimits=extrema, npoints=200, kwargs...)
    k = _kde(data; kwargs...)
    es = applydatalimits(datalimits, data)
    rgs = map(e -> range(e...; length=npoints), es)
    res = pdf(k, rgs...)
    return (rgs..., res)
end

function (d::DensityAnalysis)(le::Entry)
    options = copy(d.options)
    get!(options, :datalimits) do
        return map(v -> mapreduce(extrema, extend_extrema, v), le.positional)
    end
    entry = map(le) do p, n
        return _density(p...; pairs(n)..., options...), (;)
    end
    N = length(le.positional)
    labels = copy(le.labels)
    labels[N + 1] = "pdf"
    plottypes = [LinesFill, Heatmap, Volume]
    default_plottype = plottypes[N]
    plottype = Makie.plottype(le.plottype, default_plottype)
    return Entry(entry; plottype, labels)
end

"""
    density(; datalimits, npoints, kernel, bandwidth)

Fit a kernel density estimation of `data`.
Here, `datalimits` specifies the range for which the density should be calculated, `npoints`
is the number of points used by Makie to draw the line and `kernel` and `bandwidth` are
forwarded to `KernelDensity.kde`.
"""
density(; kwargs...) = transformation(DensityAnalysis(; kwargs...))


text:
decode_key(codec, k::String) = k

# Minimal decoder

decode(codec::MinimalCodec, x::BSONType) = x

# Default BSON decoder

decode_key(codec::BSONCodec, k::String) = codec.as_symbol ? Symbol(k) : k
decode(codec::BSONCodec, x::BSONType) = x
decode(codec::BSONCodec, x::Binary) = x.subtype == 0x00 ? x.bytes : x
decode(codec::BSONCodec, x::BSONArray) = [decode(codec, elem.value) for elem in x.elist]

function decode(codec::BSONCodec, x::Document)
    codec.dicttype(decode_key(codec, elem.key) => decode(codec, elem.value) for elem in x.elist)
end

# Canonical Extended JSON decoder

decode_key(codec::JSONCodec, k::String) = codec.as_symbol ? Symbol(k) : k
typekey(codec::JSONCodec, t) = decode_key(codec, "\$t")
decode(codec::JSONCodec, x::BSONType) = x
decode(codec::JSONCodec, x::Maxkey) = codec.dicttype(typekey(codec, "maxKey") => 1)
decode(codec::JSONCodec, x::Minkey) = codec.dicttype(typekey(codec, "minKey") => 1)
decode(codec::JSONCodec, x::Missing) = codec.dicttype(typekey(codec, "undefined") => true)
decode(codec::JSONCodec, x::ObjectId) = codec.dicttype(typekey(codec, "oid") => bytes2hex(Vector{UInt8}(x)))
decode(codec::JSONCodec, x::DateTime) = codec.dicttype(typekey(codec, "date" => decode(codec, value(x) - UNIXEPOCH)))
decode(codec::JSONCodec, x::BSONSymbol) = codec.dicttype(typekey(codec, "symbol") => x.value)
decode(codec::JSONCodec, x::Int32) = codec.dicttype(typekey(codec, "numberInt") => string(x))
decode(codec::JSONCodec, x::Int64) = codec.dicttype(typekey(codec, "numberLong") => string(x))
decode(codec::JSONCodec, x::Code) = codec.dicttype(typekey(codec, "code") => x.code)
decode(codec::JSONCodec, x::BSONArray) = [decode(codec, elem.value) for elem in x.elist]

function decode(codec::JSONCodec, x::Document)
    codec.dicttype(decode_key(codec, elem.key) => decode(codec, elem.value) for elem in x.elist)
end

function decode(codec::JSONCodec, x::DBPointer)
    codec.dicttype(typekey(codec, "dbPointer") =>
                   codec.dicttype(typekey(codec, "ref") => x.ref,
                                  typekey(codec, "id") => decode(codec, x.id)))
end

function decode(codec::JSONCodec, x::CodeWithScope)
    codec.dicttype(typekey(codec, "code") => x.code,
                   typekey(codec, "scope") => decode(codec, x.scope))
end

function decode(codec::JSONCodec, x::BSONRegex)
    codec.dicttype(typekey(codec, "regularExpression") =>
                   codec.dicttype(decode_key(codec, "pattern") => x.pattern,
                                  decode_key(codec, "options") => x.options))
end

function decode(codec::JSONCodec, x::Binary)
    codec.dicttype(typekey(codec, "binary") =>
                   codec.dicttype(decode_key(codec, "base64") =>
                                  base64encode(x.bytes),
                                  decode_key(codec, "subType") =>
                                  string(x.subtype, base=16, pad=2)))
end

function decode(codec::JSONCodec, x::Timestamp)
    codec.dicttype(typekey(codec, "timestamp") =>
                   codec.dicttype(decode_key(codec, "t") => UInt32(x.value >> 32),
                                  decode_key(codec, "i") => UInt32(x.value & typemax(UInt32))))
end

function decode(codec::JSONCodec, x::Float64)
    if isfinite(x)
        v = string(x)
    elseif isnan(x)
        v = "NaN"
    elseif x > 0
        v = "Infinity"
    else
        v = "-Infinity"
    end
    codec.dicttype(typekey(codec, "numberDouble") => v)
end


text:
function first_part(input)
    twice, thrice = 0, 0

    for line in String.(split(input))
        freqs = Array{Int8}(undef, 256)

        for i in line
            freqs[Int8(i)] += 1
        end

        if any(x -> x == 2, freqs)
            twice += 1
        end
        if any(x -> x == 3, freqs)
            thrice += 1
        end
    end

    println(twice * thrice)
end

function second_part(args)
    ids = collect(String.(split(input)))

    for i in ids, j in ids
        # println(i," ", j)
        common = find_common_letters(i, j)
        if !isempty(common)
            println(common)
            break
        end
    end

end

function find_common_letters(id1, id2)
    if id1 == id2 || length(id1) != length(id2)
        return []
    end

    found_once = false
    ids_tuple = zip(id1, id2)
    for (c1, c2) in ids_tuple
        if c1 != c2
            if found_once
                return []
            end
            found_once = true
        end
    end

    return collect(
        map((c, _) -> c,
        filter((c1, c2) -> c1 == c2,
        ids_tuple)))
end

file = open("../input/input.txt")
input = read(file, String)

first_part(input)
second_part(input)

close(file)

max_stars_repo_name: GrayJack/Advent-Of-Code

max_issues_repo_name: GrayJack/Advent-Of-Code

max_forks_repo_name: GrayJack/Advent-Of-Code


text:
# Utilities for testing
#
#       One can use the facilities provided here to simplify the testing of
#       the implementation of a subtype of AbstractPDMat
#

using PDMats, SuiteSparse, Test

const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD)
const PDMatType = HAVE_CHOLMOD ? Union{PDMat, PDSparseMat, PDiagMat} : Union{PDMat, PDiagMat}

## driver function
function test_pdmat(C::AbstractPDMat, Cmat::Matrix;
                    verbose::Int=2,             # the level to display intermediate steps
                    cmat_eq::Bool=false,        # require Cmat and Matrix(C) to be exactly equal
                    t_diag::Bool=true,          # whether to test diag method
                    t_cholesky::Bool=true,      # whether to test cholesky method
                    t_scale::Bool=true,         # whether to test scaling
                    t_add::Bool=true,           # whether to test pdadd
                    t_logdet::Bool=true,        # whether to test logdet method
                    t_eig::Bool=true,           # whether to test eigmax and eigmin
                    t_mul::Bool=true,           # whether to test multiplication
                    t_rdiv::Bool=true,          # whether to test right division (solve)
                    t_quad::Bool=true,          # whether to test quad & invquad
                    t_triprod::Bool=true,       # whether to test X_A_Xt, Xt_A_X, X_invA_Xt, and Xt_invA_X
                    t_whiten::Bool=true         # whether to test whiten and unwhiten
                    )

    d = size(Cmat, 1)
    verbose >= 1 && printstyled("Testing $(typeof(C)) with dim = $d\n", color=:blue)

    pdtest_basics(C, Cmat, d, verbose)
    pdtest_cmat(C, Cmat, cmat_eq, verbose)

    t_diag && pdtest_diag(C, Cmat, cmat_eq, verbose)
    isa(C, PDMatType) && t_cholesky && pdtest_cholesky(C, Cmat, cmat_eq, verbose)
    t_scale && pdtest_scale(C, Cmat, verbose)
    t_add && pdtest_add(C, Cmat, verbose)
    t_logdet && pdtest_logdet(C, Cmat, verbose)

    t_eig && pdtest_eig(C, Cmat, verbose)
    Imat = inv(Cmat)

    n = 5
    X = rand(eltype(C),d,n) .- convert(eltype(C),0.5)

    t_mul && pdtest_mul(C, Cmat, X, verbose)
    t_rdiv && pdtest_rdiv(C, Imat, X, verbose)
    t_quad && pdtest_quad(C, Cmat, Imat, X, verbose)
    t_triprod && pdtest_triprod(C, Cmat, Imat, X, verbose)

    t_whiten && pdtest_whiten(C, Cmat, verbose)

    verbose >= 2 && println()
end

## core testing functions

_pdt(vb::Int, s) = (vb >= 2 && printstyled("    .. testing $s\n", color=:green))

function pdtest_basics(C::AbstractPDMat, Cmat::Matrix, d::Int, verbose::Int)
    _pdt(verbose, "dim")
    @test dim(C) == d

    _pdt(verbose, "size")
    @test size(C) == (d, d)
    @test size(C, 1) == d
    @test size(C, 2) == d
    @test size(C, 3) == 1

    _pdt(verbose, "ndims")
    @test ndims(C) == 2

    _pdt(verbose, "length")
    @test length(C) == d * d

    _pdt(verbose, "eltype")
    @test eltype(C) == eltype(Cmat)
#    @test eltype(typeof(C)) == eltype(typeof(Cmat))

    _pdt(verbose, "index")
    @test all(C[i] == Cmat[i] for i in 1:(d^2))
    @test all(C[i, j] == Cmat[i, j] for j in 1:d, i in 1:d)

    _pdt(verbose, "isposdef")
    @test isposdef(C)

    _pdt(verbose, "ishermitian")
    @test ishermitian(C)
end

function pdtest_cmat(C::AbstractPDMat, Cmat::Matrix, cmat_eq::Bool, verbose::Int)
    _pdt(verbose, "full")
    if cmat_eq
        @test Matrix(C) == Cmat
    else
        @test Matrix(C) ≈ Cmat
    end
end

function pdtest_diag(C::AbstractPDMat, Cmat::Matrix, cmat_eq::Bool, verbose::Int)
    _pdt(verbose, "diag")
    if cmat_eq
        @test diag(C) == diag(Cmat)
    else
        @test diag(C) ≈ diag(Cmat)
    end
end

function pdtest_cholesky(C::Union{PDMat, PDiagMat}, Cmat::Matrix, cmat_eq::Bool, verbose::Int)
    _pdt(verbose, "cholesky")
    if cmat_eq
        @test cholesky(C).U == cholesky(Cmat).U
    else
        @test cholesky(C).U ≈ cholesky(Cmat).U
    end
end

if HAVE_CHOLMOD
    function pdtest_cholesky(C::PDSparseMat, Cmat::Matrix, cmat_eq::Bool, verbose::Int)
        _pdt(verbose, "cholesky")
        # We special case PDSparseMat because we can't perform equality checks on
        # `SuiteSparse.CHOLMOD.Factor`s and `SuiteSparse.CHOLMOD.FactorComponent`s
        @test diag(cholesky(C)) ≈ diag(cholesky(Cmat).U)
        # NOTE: `==` also doesn't work because `diag(cholesky(C))` will return `Vector{Float64}`
        # even if the inputs are `Float32`s.
    end
end

function pdtest_scale(C::AbstractPDMat, Cmat::Matrix, verbose::Int)
    _pdt(verbose, "scale")
    @test Matrix(C * convert(eltype(C),2)) ≈ Cmat * convert(eltype(C),2)
    @test Matrix(convert(eltype(C),2) * C) ≈ convert(eltype(C),2) * Cmat
end

function pdtest_add(C::AbstractPDMat, Cmat::Matrix, verbose::Int)
    M = rand(eltype(C),size(Cmat))
    _pdt(verbose, "add")
    @test C + M ≈ Cmat + M
    @test M + C ≈ M + Cmat

    _pdt(verbose, "add_scal")
    @test pdadd(M, C, convert(eltype(C),2)) ≈ M + Cmat * convert(eltype(C),2)

    _pdt(verbose, "add_scal!")
    R = M + Cmat * convert(eltype(C),2)
    Mr = pdadd!(M, C, convert(eltype(C),2))
    @test Mr === M
    @test Mr ≈ R
end

function pdtest_logdet(C::AbstractPDMat, Cmat::Matrix, verbose::Int)
    _pdt(verbose, "logdet")
    # default tolerance in isapprox is different on 0.4. rtol argument can be deleted
    # ≈ form used when 0.4 is no longer supported
    @test isapprox(logdet(C), logdet(Cmat), rtol=sqrt(max(eps(real(eltype(C))), eps(real(eltype(Cmat))))))
end

function pdtest_eig(C::AbstractPDMat, Cmat::Matrix, verbose::Int)
    _pdt(verbose, "eigmax")
    @test eigmax(C) ≈ eigmax(Cmat)

    _pdt(verbose, "eigmin")
    @test eigmin(C) ≈ eigmin(Cmat)
end

function pdtest_mul(C::AbstractPDMat, Cmat::Matrix, verbose::Int)
    n = 5
    X = rand(eltype(C),dim(C), n)

    _pdt(verbose, "multiply")
    @test C * X ≈ Cmat * X

    for i = 1:n
        xi = vec(copy(X[:,i]))
        @test C * xi ≈ Cmat * xi
    end
end

function pdtest_mul(C::AbstractPDMat, Cmat::Matrix, X::Matrix, verbose::Int)
    _pdt(verbose, "multiply")
    @test C * X ≈ Cmat * X

    y = similar(C * X, size(C, 1))
    ymat = similar(Cmat * X, size(Cmat, 1))
    for i = 1:size(X,2)
        xi = vec(copy(X[:,i]))
        @test C * xi ≈ Cmat * xi

        mul!(y, C, xi)
        mul!(ymat, Cmat, xi)
        @test y ≈ ymat
    end
end

function pdtest_rdiv(C::AbstractPDMat, Imat::Matrix, X::Matrix, verbose::Int)
    _pdt(verbose, "rdivide")
    @test C \ X ≈ Imat * X

    for i = 1:size(X,2)
        xi = vec(copy(X[:,i]))
        @test C \ xi ≈ Imat * xi
    end
end

function pdtest_quad(C::AbstractPDMat, Cmat::Matrix, Imat::Matrix, X::Matrix, verbose::Int)
    n = size(X, 2)

    _pdt(verbose, "quad")
    r_quad = zeros(eltype(C),n)
    for i = 1:n
        xi = vec(X[:,i])
        r_quad[i] = dot(xi, Cmat * xi)
        @test quad(C, xi) ≈ r_quad[i]
        @test quad(C, view(X,:,i)) ≈ r_quad[i]
    end
    @test quad(C, X) ≈ r_quad

    _pdt(verbose, "invquad")
    r_invquad = zeros(eltype(C),n)
    for i = 1:n
        xi = vec(X[:,i])
        r_invquad[i] = dot(xi, Imat * xi)
        @test invquad(C, xi) ≈ r_invquad[i]
        @test invquad(C, view(X,:,i)) ≈ r_invquad[i]
    end
    @test invquad(C, X) ≈ r_invquad
end

function pdtest_triprod(C::AbstractPDMat, Cmat::Matrix, Imat::Matrix, X::Matrix, verbose::Int)
    Xt = copy(transpose(X))

    _pdt(verbose, "X_A_Xt")
    # default tolerance in isapprox is different on 0.4. rtol argument can be deleted
    # ≈ form used when 0.4 is no longer supported
    lhs, rhs = X_A_Xt(C, Xt), Xt * Cmat * X
    @test isapprox(lhs, rhs, rtol=sqrt(max(eps(real(float(eltype(lhs)))), eps(real(float(eltype(rhs)))))))

    _pdt(verbose, "Xt_A_X")
    lhs, rhs = Xt_A_X(C, X), Xt * Cmat * X
    @test isapprox(lhs, rhs, rtol=sqrt(max(eps(real(float(eltype(lhs)))), eps(real(float(eltype(rhs)))))))

    _pdt(verbose, "X_invA_Xt")
    @test X_invA_Xt(C, Xt) ≈ Xt * Imat * X

    _pdt(verbose, "Xt_invA_X")
    @test Xt_invA_X(C, X) ≈ Xt * Imat * X
end

function pdtest_whiten(C::AbstractPDMat, Cmat::Matrix, verbose::Int)
    Y = PDMats.chol_lower(Cmat)
    Q = qr(convert(Array{eltype(C),2},randn(size(Cmat)))).Q
    Y = Y * Q'                    # generate a matrix Y such that Y * Y' = C
    @test Y * Y' ≈ Cmat
    d = dim(C)

    _pdt(verbose, "whiten")
    Z = whiten(C, Y)
    @test Z * Z' ≈ Matrix{eltype(C)}(I, d, d)
    for i = 1:d
        @test whiten(C, Y[:,i]) ≈ Z[:,i]
    end

    _pdt(verbose, "whiten!")
    Z2 = copy(Y)
    whiten!(C, Z2)
    @test Z ≈ Z2

    _pdt(verbose, "unwhiten")
    X = unwhiten(C, Z)
    @test X * X' ≈ Cmat
    for i = 1:d
        @test unwhiten(C, Z[:,i]) ≈ X[:,i]
    end

    _pdt(verbose, "unwhiten!")
    X2 = copy(Z)
    unwhiten!(C, X2)
    @test X ≈ X2

    _pdt(verbose, "whiten-unwhiten")
    @test unwhiten(C, whiten(C, Matrix{eltype(C)}(I, d, d))) ≈ Matrix{eltype(C)}(I, d, d)
    @test whiten(C, unwhiten(C, Matrix{eltype(C)}(I, d, d))) ≈ Matrix{eltype(C)}(I, d, d)
end


text:
"""
    weighted_iterator(d)

Return an iterator through pairs of the values and probabilities in distribution `d`.

This is designed to speed up value iteration. Distributions are encouraged to provide a custom optimized implementation if possible.
"""
weighted_iterator(d) = (x=>pdf(d, x) for x in iterator(d))


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module Benchmark
include("load.jl")
include("../shared/load.jl")
using ADPerfTest
using TestLoader
using GMMData
using BAData
using HandData
using LSTMData
using DataUtils

export run_benchmark

const measurable_time_not_achieved = -1

"""
    find_repeats_for_minimum_measurable_time!(context, minimum_measurable_time, func)

Finds the number of repeats that `func` needs to be asked to compute its objective
in order to to run for at least `minimum_measurable_time` seconds.
`func` must have a method `func(::typeof(context), repeats::Int)`. That is the method
that will be tested.

Returns a tuple `(repeats, min_sample, total_time)`, where
`repeats` is the found necessary number of repeats,
`min_sample` is the minimal encountered value of
`@elapsed func(context, repeats) / repeats`
`total_time` is the total time spent on running `func` in seconds.
"""
function find_repeats_for_minimum_measurable_time!(context::Any, minimum_measurable_time::Float64, func::Function)
    total_time = 0.0
    min_sample = floatmax(Float64)
    repeats = 1
    max_possible_power_of_two = typemax(Int) ÷ 2 + 1
    while true
        current_run_time = @elapsed func(context, repeats)
        if current_run_time > minimum_measurable_time
            current_sample = current_run_time / repeats
            min_sample = min(min_sample, current_sample)
            total_time += current_run_time
            break
        end
        # The next iteration will overflow a loop counter that's why we recognize that we cannot reach the minimum measurable time.
        if repeats == max_possible_power_of_two
            repeats = measurable_time_not_achieved
            break
        end
        repeats *= 2
    end
    repeats, min_sample, total_time
end

"Benchmarks func according to the methodology described in docs/Methodology.md"
function measure_shortest_time!(context::Any, minimum_measurable_time::Float64, nruns::Int, time_limit::Float64, func::Function)::Float64
    precompile(func, (typeof(context), Int))
    repeats, min_sample, total_time = find_repeats_for_minimum_measurable_time!(context, minimum_measurable_time, func)
    if repeats == measurable_time_not_achieved
        throw(ErrorException("It was not possible to reach the number of repeats sufficient to achieve the minimum measurable time."))
    end
    run = 1
    while run < nruns && total_time < time_limit
        current_run_time = @elapsed func(context, repeats)
        min_sample = min(min_sample, current_run_time / repeats)
        total_time += current_run_time
        run += 1
    end
    min_sample
end

"Performs the entire benchmark process according to the methodology described in docs/Methodology.md"
function run_benchmark(
    input::Input,
    input_name::AbstractString,
    module_name::AbstractString,
    output_prefix::AbstractString,
    module_display_name::AbstractString,
    minimum_measurable_time::Float64,
    nruns_f::Int,
    nruns_J::Int,
    time_limit::Float64
) where Input
    output = create_empty_output_for(input)
    test = get_test_for(input, output, module_name)
    # Callbacks defined in test exist only in the latest world
    Base.invokelatest() do
        test.prepare!(test.context, input)

        objective_time = measure_shortest_time!(test.context, minimum_measurable_time, nruns_f, time_limit, test.calculate_objective!)
        derivative_time = measure_shortest_time!(test.context, minimum_measurable_time, nruns_J, time_limit, test.calculate_jacobian!)

        test.output!(output, test.context)

        save_time_to_file(times_file_name(output_prefix, input_name, module_display_name), objective_time, derivative_time)
        save_output_to_file(output, output_prefix, input_name, module_display_name)
    end
end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
"""
    pctOnly(df::DataFrame, ColPct::String)

Creates a new dataframe without rows without precipitation.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColPct::String`: The name of the column of `df` that allows to know if it has rained or not during the time division used.
"""
function pctOnly(df::DataFrame, ColPct::String)
    filter(Symbol(ColPct) => x -> x > 0, df)
end

"""
    pctOnly!(df::DataFrame, ColPct::Symbol)

Deletes rows without precipitation.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColPct::String`: The name of the column of `df` that allows to know if it has rained or not during the time division used.
"""
function pctOnly!(df::DataFrame, ColPct::String)
    filter!(Symbol(ColPct) => x -> x > 0, df)
end

"""
    overflowOnly(df::DataFrame, ColOverflow::Strings)

Creates a new dataframe without rows without overflow.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColOverflow::String`: The name of the column of `df` containing the binary vector of overflows.
"""
function overflowOnly(df::DataFrame, ColOverflow::String)
    filter(Symbol(ColOverflow) => x -> x == 1, df)
end

"""
    overflowOnly!(df::DataFrame, ColOverflow::String)

Deletes rows without overflow.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColOverflow::String`: The name of the column of `df` containing the binary vector of overflows.
"""
function overflowOnly!(df::DataFrame, ColOverflow::String)
    filter!(Symbol(ColOverflow) => x -> x == 1, df)
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
export ElectronVolt, Torr, Atmosphere, Degree

const ElectronVolt = NonSIUnit{typeof(Joule),:eV}()
convert(::Type{SIQuantity},::typeof(ElectronVolt)) = 1.60217656535e-19Joule

const Torr = NonSIUnit{typeof(Pascal),:torr}()
convert(::Type{SIQuantity},::typeof(Torr)) = 133.322368Pascal

const Atmosphere = NonSIUnit{typeof(Pascal),:atm}()
convert(::Type{SIQuantity},::typeof(Atmosphere)) = 101325Pascal

const Degree = NonSIUnit{typeof(Radian),:deg}()
convert(::Type{SIQuantity},::typeof(Degree)) = π/180.*Radian

for (func,funcd) in ((:sin,:sind),
                     (:cos,:cosd),
                     (:tan,:tand),
                     (:cot,:cotd),
                     (:sec,:secd),
                     (:csc,:cscd))
    @eval $func{T}(θ::NonSIQuantity{T,$(typeof(Degree))}) = $funcd(θ.val)
end


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # Robust uncertainty sets

# Computes the Value at Risk for a data-driven uncertainty set; see "Data-Driven
# Robust Optimization" (Bertsimas 2013), section 6.1 for details. Closed-form
# expressions for the optimal value are available.

using JuMP
import SCS
import LinearAlgebra
import Test

function example_robust_uncertainty()
    R = 1
    d = 3
    𝛿 = 0.05
    ɛ = 0.05
    N = ceil((2 + 2 * log(2 / 𝛿))^2) + 1
    c = randn(d)
    μhat = rand(d)
    M = rand(d, d)
    Σhat = 1 / (d - 1) * (M - ones(d) * μhat')' * (M - ones(d) * μhat')
    Γ1(𝛿, N) = R / sqrt(N) * (2 + sqrt(2 * log(1 / 𝛿)))
    Γ2(𝛿, N) = 2 * R^2 / sqrt(N) * (2 + sqrt(2 * log(2 / 𝛿)))
    model = Model(SCS.Optimizer)
    set_silent(model)
    @variable(model, Σ[1:d, 1:d], PSD)
    @variable(model, u[1:d])
    @variable(model, μ[1:d])
    @constraint(model, [Γ1(𝛿 / 2, N); μ - μhat] in SecondOrderCone())
    @constraint(model, [Γ2(𝛿 / 2, N); vec(Σ - Σhat)] in SecondOrderCone())
    @SDconstraint(model, [((1 - ɛ) / ɛ) (u - μ)'; (u - μ) Σ] >= 0)
    @objective(model, Max, LinearAlgebra.dot(c, u))
    optimize!(model)
    I = Matrix(1.0 * LinearAlgebra.I, d, d)
    exact =
        LinearAlgebra.dot(μhat, c) +
        Γ1(𝛿 / 2, N) * LinearAlgebra.norm(c) +
        sqrt((1 - ɛ) / ɛ) * sqrt(LinearAlgebra.dot(c, (Σhat + Γ2(𝛿 / 2, N) * I) * c))
    Test.@test objective_value(model) ≈ exact atol = 1e-3
    return
end

example_robust_uncertainty()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is part of Kpax3. License is MIT.

RecipesBase.@userplot Kpax3PlotP

RecipesBase.@recipe function plotp(g::Kpax3PlotP;
                                   clusterorder=:auto,
                                   clusterlabel=:auto)
  if (length(g.args) != 2) ||
     (typeof(g.args[1]) != Vector{Int}) ||
     (typeof(g.args[2]) != Matrix{Float64})
      throw(KInputError("Two arguments required. Got: $(typeof(g.args))"))
  end

  R, P = g.args

  n = length(R)
  k = maximum(R)

  if (clusterorder == :auto) || (length(clusterorder) != k)
    clusterorder = collect(1:k)
  end

  if (clusterlabel == :auto) || (length(clusterlabel) != k)
    clusterlabel = string.(clusterorder)
  end

  (ord, mid, sep) = reorderunits(R, P, clusterorder)

  mid .-= 0.5
  sep .-= 0.5

  N = n^2

  x = zeros(Int, N)
  y = zeros(Int, N)
  z = fill("", N)

  idx = 1
  pij = 0.0
  for j = 1:n, i = 1:n
    x[idx] = j
    y[idx] = i

    pij = P[ord[i], ord[j]]
    z[idx] = if pij <= 0.2
      "#FFFFFF"
    elseif pij <= 0.4
      "#F0F0F0"
    elseif pij <= 0.6
      "#BDBDBD"
    elseif pij <= 0.8
      "#636363"
    else
      "#000000"
    end

    idx += 1
  end

  # try to minimize the number of rectangles. We will do this greedly by
  # dividing the image into maximal adjacent rectangles
  # Hopefully the SVG file won't be too big
  #
  # TODO: improve the algorithm
  processed = falses(n, n)

  len = zero(Int)
  xmin = zeros(Int, len)
  xmax = zeros(Int, len)
  ymin = zeros(Int, len)
  ymax = zeros(Int, len)
  col = fill("", len)

  # elements on the diagonal have always p = 1.0
  val = "#000000"

  j = 1
  while j <= n
    jmax = expandsquarediag(j, n, val, z)

    push!(xmin, j - 1)
    push!(xmax, jmax)
    push!(ymin, j - 1)
    push!(ymax, jmax)
    push!(col, val)

    len += 1

    processed[j:jmax, j:jmax] .= true
    j = jmax + 1
  end

  # scan the lower triangular matrix
  j = 1
  while j < n
    # start a new column from the first non processed element
    i = j + 1
    while i <= n
      if !processed[i, j]
        val = z[LinearIndices((n, n))[i, j]]

        (imax, jmax) = expandrect(i, j, n, n, val, z, processed)

        # lower triangular
        push!(xmin, j - 1)
        push!(xmax, jmax)
        push!(ymin, i - 1)
        push!(ymax, imax)
        push!(col, val)

        len += 1

        # upper triangular (transpose)
        push!(xmin, i - 1)
        push!(xmax, imax)
        push!(ymin, j - 1)
        push!(ymax, jmax)
        push!(col, val)

        len += 1

        processed[i:imax, j:jmax] .= true
        i = imax + 1
      else
        while i <= n && processed[i, j]
          i += 1
        end
      end
    end

    j += 1
  end

  xborder = zeros(Float64, 5, k)
  yborder = zeros(Float64, 5, k)
  for i = 1:k
    xborder[:, i] = [sep[i]; sep[i]; sep[i + 1]; sep[i + 1]; sep[i]]
    yborder[:, i] = [sep[i]; sep[i + 1]; sep[i + 1]; sep[i]; sep[i]]
  end

  mcol = ["#FFFFFF" "#F0F0F0" "#BDBDBD" "#636363" "#000000"]
  mlab = ["[0.0, 0.2]" "(0.2, 0.4]" "(0.4, 0.6]" "(0.6, 0.8]" "(0.8, 1.0]"]

  rect_x = zeros(Float64, 5 * len)
  rect_y = zeros(Float64, 5 * len)
  for i = 1:len
    idx = LinearIndices((5, len))[:, i]
    rect_x[idx] .= [xmin[i]; xmax[i]; xmax[i]; xmin[i]; NaN]
    rect_y[idx] .= [ymin[i]; ymin[i]; ymax[i]; ymax[i]; NaN]
  end

  # start by plotting cluster borders
  RecipesBase.@series begin
    seriestype := :path
    linecolor := :black
    label := ""
    legend := :none

    xborder, yborder
  end

  # add the rectangles and fill them with color
  RecipesBase.@series begin
    seriestype := :shape
    fillcolor := col
    linealpha := 0.0
    linewidth := 0.0
    linestyle := :dot
    label := ""

    rect_x, rect_y
  end

  # complete the plot by adding the legend
  RecipesBase.@series begin
    seriestype := :scatter
    background_color := :white
    grid := false
    markershape := :rect
    markercolor := mcol
    label := mlab
    xlims := (0, n)
    xticks := (mid, clusterlabel)
    ylims := (0, n)
    yticks := (mid, clusterlabel)
    yflip := true

    html_output_format --> :svg
    size --> (800, 600)
    window_title --> ""
    title --> ""
    legend --> :right
    xlabel --> "Samples by cluster"
    ylabel --> "Samples by cluster"

    fill(-2.0, 1, length(mcol)), fill(-2.0, 1, length(mcol))
  end
end

RecipesBase.@userplot Kpax3PlotC

RecipesBase.@recipe function plotc(g::Kpax3PlotC)
  if (length(g.args) != 3) ||
     (typeof(g.args[1]) != Vector{Int}) ||
     (typeof(g.args[2]) != Vector{Float64}) ||
     (typeof(g.args[3]) != Matrix{Float64})
      throw(KInputError("Three arguments required. Got: $(typeof(g.args))"))
  end

  site, freq, C = g.args

  m = site[end]
  M = length(site)

  noise = zeros(Float64, m)
  weak = zeros(Float64, m)

  key = 1
  for b in 1:M
    if site[b] != key
      weak[key] += noise[key]
      key += 1
    end

    noise[key] += freq[b] * C[1, b]
    weak[key] += freq[b] * C[2, b]
  end

  weak[key] += noise[key]

  mcol = ["#FFFFFF" "#808080" "#000000"]
  mlab = ["Noise" "Weak signal" "Strong signal"]

  # start with a complete black plot (strong signal)
  RecipesBase.@series begin
    seriestype := :path
    linecolor := :black
    fill := (0, :black)
    label := ""
    legend := :none

    1:m, ones(Float64, m)
  end

  # add the weak signal on top
  RecipesBase.@series begin
    seriestype := :path
    linecolor := "#808080"
    fill := (0, "#808080")
    label :=""
    legend := :none

    1:m, weak
  end

  # finally add the white noise
  RecipesBase.@series begin
    seriestype := :path
    linecolor := :white
    fill := (0, :white)
    label :=""
    legend := :none

    1:m, noise
  end

  # add the legend
  RecipesBase.@series begin
    seriestype := :scatter
    background_color := :white
    grid := false
    markershape := :rect
    markercolor := mcol
    label := mlab
    legend := :right
    xlims := (1, m)
    ylims := (0, 1)

    html_output_format --> :svg
    size --> (800, 600)
    window_title --> ""
    title --> ""
    legend --> :right
    xlabel --> "Site"
    ylabel --> "Posterior probability"

    fill(-2, 1, length(mcol)), fill(-2, 1, length(mcol))
  end
end

RecipesBase.@userplot Kpax3PlotD

RecipesBase.@recipe function plotd(g::Kpax3PlotD;
                                   clusterorder=:auto,
                                   clusterlabel=:auto)
  if length(g.args) != 2
      throw(KInputError("Two arguments required. Got: $(typeof(g.args))"))
  end

  x, state = g.args

  (m, n) = size(x.data)
  M = length(x.ref)

  if (clusterorder == :auto) || (length(clusterorder) != state.k)
    clusterorder = collect(1:state.k)
  end

  if (clusterlabel == :auto) || length(clusterlabel) != state.k
    clusterlabel = string.(clusterorder)
  end

  v = zeros(Float64, state.k)
  for i in 1:n
    v[state.R[i]] += 1.0
  end
  v /= n

  v = [0.0; v[clusterorder]]

  cumsum!(v, v)

  xax = collect(1:M)

  yax = zeros(Float64, state.k)
  for g in 1:state.k
    yax[g] = (v[g] + v[g + 1]) / 2
  end
  round.(yax, digits=3)

  # colors
  colset = [("A", "Ala", "#FF3232");
            ("R", "Arg", "#FF7F00");
            ("N", "Asn", "#CAB2D6");
            ("D", "Asp", "#FDBF6F");
            ("C", "Cys", "#33A02C");
            ("E", "Glu", "#1F4E78");
            ("Q", "Gln", "#FF3300");
            ("G", "Gly", "#660066");
            ("H", "His", "#B2DF8A");
            ("I", "Ile", "#CC6600");
            ("L", "Leu", "#375623");
            ("K", "Lys", "#A6CEE3");
            ("M", "Met", "#CC3399");
            ("F", "Phe", "#FB9A99");
            ("P", "Pro", "#C81414");
            ("S", "Ser", "#1F78B4");
            ("T", "Thr", "#6A3D9A");
            ("W", "Trp", "#FFFF99");
            ("Y", "Tyr", "#FF66CC");
            ("V", "Val", "#993300");
            ("U", "Sec", "#00FFFF");
            ("O", "Pyl", "#FFFF00");
            ("B", "Asx", "#E4B9A3");
            ("Z", "Glx", "#8F413C");
            ("J", "Xle", "#825E12");
            ("X", "Xaa", "#000000");
            ("-",   "-", "#E0E0E0");
            ("*",   "*", "#000000");
            ( "",    "", "#FFFFFF")]

  coltable = Dict(map(a -> (a[1], a[3]), colset))

  # indices
  idx = 1
  s = 0
  t = 1
  b = 1
  g = 1
  h = 1

  # temporary values
  c = 0x01
  w = 0
  flag = false

  z = fill("#FFFFFF", state.k * M)

  for j in 1:M
    if x.ref[j] == UInt8('.')
      s += 1
      c = 0x01

      while (b <= m) && (x.key[b] == s)
        if state.C[state.cl[1], b] > c
          c = state.C[state.cl[1], b]
        end
        b += 1
      end

      if c > 0x02
        for g in 1:state.k
          h = clusterorder[g]
          w = t
          flag = false

          while !flag && w < b
            flag = state.C[state.cl[h], w] == 0x04
            w += 1
          end

          if flag
            z[idx] = coltable[uppercase(string(Char(x.val[w - 1])))]
          end

          idx += 1
        end
      else
        idx += state.k
      end

      t = b
    else
      idx += state.k
    end
  end

  # try to minimize the number of rectangles. We will do this greedly by
  # dividing the image into maximal adjacent rectangles
  # Hopefully the SVG file won't be too big
  #
  # TODO: improve the algorithm
  processed = falses(state.k, M)

  len = zero(Int)
  xmin = zeros(Int, len)
  xmax = zeros(Int, len)
  ymin = zeros(Float64, len)
  ymax = zeros(Float64, len)
  col = fill("", len)

  j = 1
  while j <= M
    # start a new column from the first non processed element
    i = 1
    while i <= state.k
      if !processed[i, j]
        val = z[LinearIndices((state.k, M))[i, j]]

        (imax, jmax) = expandrect(i, j, state.k, M, val, z, processed)

        push!(xmin, j - 1)
        push!(xmax, jmax)
        push!(ymin, v[i])
        push!(ymax, v[imax + 1])
        push!(col, val)

        len += 1

        processed[i:imax, j:jmax] .= true
        i = imax + 1
      else
        while i <= state.k && processed[i, j]
          i += 1
        end
      end
    end

    j += 1
  end

  mcol = reshape(map(a -> a[3], colset[1:26]), 1, 26)
  mlab = reshape(map(a -> a[2], colset[1:26]), 1, 26)

  rect_x = zeros(Float64, 5 * len)
  rect_y = zeros(Float64, 5 * len)
  for i = 1:len
    idx = LinearIndices((5, len))[:, i]
    rect_x[idx] .= [xmin[i]; xmax[i]; xmax[i]; xmin[i]; NaN]
    rect_y[idx] .= [ymin[i]; ymin[i]; ymax[i]; ymax[i]; NaN]
  end

  # draw the rectangles
  RecipesBase.@series begin
    seriestype := :shape
    fillcolor := col
    linecolor := col
    label := ""
    legend := :none

    rect_x, rect_y
  end

  # add cluster separators
  RecipesBase.@series begin
    seriestype := :hline
    linecolor := :black
    label := ""
    legend := :none

    v
  end

  # complete the plot by adding the legend
  RecipesBase.@series begin
    seriestype := :scatter
    background_color := :white
    grid := false
    markershape := :rect
    markercolor := mcol
    label := mlab
    xlims := (1, M)
    ylims := (0, 1)
    yticks := (yax, clusterlabel)
    yflip := true

    html_output_format --> :svg
    size --> (800, 600)
    window_title --> ""
    title --> ""
    legend --> :right
    xlabel --> "Site"
    ylabel --> "Cluster"

    fill(-2, 1, length(mcol)), fill(-2, 1, length(mcol))
  end
end


text:
import Base.alignment

# A bit of a mess but does the job...
function _blockarray_print_matrix_row(io::IO,
        X::AbstractVecOrMat, A::Vector,
        i::Integer, cols::AbstractVector, sep::AbstractString)
    cumul = 0
    block = 1

    row_buf = IOBuffer()

    row_sum = cumulsizes(X,1)[2:end] .- 1
    if ndims(X) == 2
        col_sum = (cumulsizes(X,2)[2:end] .- 1)[1:end-1]
    end

    # Loop over row
    for k = 1:length(A)
        n_chars = 0
        j = cols[k]

        if isassigned(X,Int(i),Int(j)) # isassigned accepts only `Int` indices
            x = X[i,j]
            a = Base.alignment(io, x)
            sx = sprint(show, x; context=io, sizehint=0)
        else
            a = Base.undef_ref_alignment
            sx = Base.undef_ref_str
        end
        l = repeat(" ",A[k][1]-a[1]) # pad on left and right as needed
        r = repeat(" ", A[k][2]-a[2])
        prettysx = Base.replace_in_print_matrix(X,i,j,sx)
        # Print the element
        print(io, l, prettysx, r)

        # Jump forward
        n_chars += length(l) + length(prettysx) + length(r) + 2

        cumul += 1
        if ndims(X) == 2
            # Have accumulated enough for the block, should print a |
            if block < length(cumulsizes(X,2)) - 1 && cumul == blocksize(X, 2, block)
                block += 1
                cumul = 0
                print(io, "  │")
                n_chars += 3
            end
        end

        if k == 1
            n_chars -= 2
        end

        if i in row_sum
            print(row_buf, "─"^(n_chars-1))
            if ndims(X) == 2 && k in col_sum
                print(row_buf, "┼")
            else
                print(row_buf, "─")
            end
        end

        if k < length(A); print(io, sep); end
    end

    if i < size(X, 1)
        row_str = String(take!(row_buf))
        if length(row_str) > 0
            print(io, "\n ")
            print(io, row_str)
        end
    end
end

Base.print_matrix_row(io::IO,
        X::AbstractBlockVecOrMat, A::Vector,
        i::Integer, cols::AbstractVector, sep::AbstractString) =
        _blockarray_print_matrix_row(io, X, A, i, cols, sep)


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 4, n}) where n
    io = openInput(puzzle)
    passports = split(read(io, String), "\n\n", keepempty=false)
    # Map passports to an array of fields where each field is given by a type and value
    map(passport -> split.(split(passport, r"\n| ", keepempty=false), ':'), passports)
end

struct PassportField{field}
    value::AbstractString
end

function parsePassword(s)
    range, char, password = split(s, ' ')
    left, right = parse.(Int, split(range, '-'))

    return (left:right, char[1], password)
end

function passpordIsValid1(passpord)
    return length(filter(x -> x[1] != "cid", passpord)) == 7
end

function solve(::Puzzle{2020, 4, 1}, input)
    credentials = input
    return filter(passpordIsValid1, credentials) |> length
end

isValid(x, y) = false
function isValid(::Val{:byr}, input)
    byr = parse(Int, input)
    if (! (byr in 1920:2002))
        return false
    end
    return true
end

function passpordIsValid2(passpord)
    fp = filter(x -> x[1] != "cid", passpord)
    if (length(fp) != 7)
        return false
    end
    p = Dict(fp);

    if (! isValidByr(parse(Int, p["byr"])))
        return false
    end

    if (! isValidIyr(parse(Int, p["iyr"])))
        return false
    end

    if (! isValidEyr(parse(Int, p["eyr"])))
        return false
    end

    if (! isValidHeight(p["hgt"]))
        return false
    end

    if (! isValidHcl(p["hcl"]))
        return false
    end

    if (! isValidEcl(p["ecl"]))
        return false
    end

    if (! isValidPid(p["pid"]))
       return false
    end

    return true
end

isValidByr(byr) = 1920 <= byr <= 2002
isValidIyr(iyr) = 2010 <= iyr <= 2020
isValidEyr(eyr) = 2020 <= eyr <= 2030

function isValidHeight(hgt)
    m = match(r"^(\d.*)(cm|in)$", hgt)

    if (isnothing(m))
        return false
    end

    unit = m.captures[2]
    height = parse(Int, m[1])

    if (unit == "cm")
        return 150 <= height <= 193
    elseif (unit == "in")
        return 59 <= height <= 76
    end

    return false
end

isValidEcl(ecl) = ecl in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
isValidHcl(hcl) = occursin(r"\#[0-9a-f]{6}", hcl)

function isValidPid(pid)
    try
        pidi = parse(Int, pid)
        if (length(pid) != 9)
            return false
        end
    catch e
        return false
    end
    return true
end

function solve(::Puzzle{2020, 4, 2}, input)
    credentials = input
    return  filter(passpordIsValid2, credentials) |> length
end

@testset "Day 4" begin
    @test(solve(Puzzle(4, 1), input(Puzzle(4, 1))) == 242)
    @test(solve(Puzzle(4, 2), input(Puzzle(4, 2))) == 186)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Trilinos"
version = v"12.12.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/trilinos/Trilinos.git", "512a9e81183c609ab16366a9b09d70d37c6af8d4")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mkdir trilbuild
cd trilbuild
install_license ${WORKSPACE}/srcdir/Trilinos/LICENSE
SRCDIR="/workspace/srcdir/Trilinos"
FLAGS="-O3 -fPIC"
cmake -G "Unix Makefiles" -DCMAKE_C_COMPILER=gcc -DCMAKE_CXX_COMPILER=g++ -DCMAKE_Fortran_COMPILER=gfortran -DCMAKE_CXX_FLAGS="$FLAGS" -DCMAKE_C_FLAGS="$FLAGS" -DCMAKE_Fortran_FLAGS="$FLAGS" -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_MAKE_PROGRAM="make" -DTrilinos_ENABLE_NOX=ON -DNOX_ENABLE_LOCA=ON -DTrilinos_ENABLE_EpetraExt=ON   -DEpetraExt_BUILD_BTF=ON   -DEpetraExt_BUILD_EXPERIMENTAL=ON -DEpetraExt_BUILD_GRAPH_REORDERINGS=ON -DTrilinos_ENABLE_TrilinosCouplings=ON -DTrilinos_ENABLE_Ifpack=ON -DTrilinos_ENABLE_Isorropia=ON -DTrilinos_ENABLE_AztecOO=ON -DTrilinos_ENABLE_Belos=ON -DTrilinos_ENABLE_Teuchos=ON -DTeuchos_ENABLE_COMPLEX=ON -DTrilinos_ENABLE_Amesos=ON -DAmesos_ENABLE_KLU=ON -DTrilinos_ENABLE_Sacado=ON -DTrilinos_ENABLE_Kokkos=OFF -DTrilinos_ENABLE_ALL_OPTIONAL_PACKAGES=OFF -DTrilinos_ENABLE_CXX11=ON -DTPL_ENABLE_AMD=ON -DAMD_LIBRARY_DIRS="/${libdir}" -DAMD_LIBRARY_NAMES="libsuitesparseconfig.${dlext};libamd.${dlext};libklu.${dlext};libcolamd.${dlext};libbtf.${dlext}" -DAMD_INCLUDE_DIRS="${prefix}/include" -DTPL_ENABLE_UMFPACK=ON -DUMFPACK_LIBRARY_DIRS="/${libdir}" -DAMD_LIBRARY_NAMES="libumfpack.${dlext}" -DTPL_UMFPACK_INCLUDE_DIRS="${prefix}/include" -DTPL_ENABLE_BLAS=ON -DTPL_ENABLE_LAPACK=ON -DBLAS_LIBRARY_DIRS="${prefix}/lib" -DBLAS_LIBRARY_NAMES="libopenblas.${dlext}" -DLAPACK_LIBRARY_DIRS="${prefix}/lib" -DLAPACK_LIBRARY_NAMES="libopenblas.${dlext}" -DCMAKE_BUILD_TYPE=Release $SRCDIR
make -j${nprocs}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter(p -> (!Sys.iswindows(p) && !Sys.isapple(p)), supported_platforms())

platforms = expand_cxxstring_abis(platforms)
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libaztecoo", :libaztecoo),
    LibraryProduct("liblocalapack", :liblocalapack),
    LibraryProduct("libifpack", :libifpack),
    LibraryProduct("libloca", :libloca),
    LibraryProduct("libzoltan", :libzoltan),
    LibraryProduct("libisorropia", :libisorropia),
    LibraryProduct("libbelos", :libbelos),
    LibraryProduct("libteuchosparameterlist", :libteuchosparameterlist),
    LibraryProduct("libnoxlapack", :libnoxlapack),
    LibraryProduct("libnox", :libnox),
    LibraryProduct("libteuchoscore", :libteuchoscore),
    LibraryProduct("libsimpi", :libsimpi),
    LibraryProduct("libteuchosremainder", :libteuchosremainder),
    LibraryProduct("liblocaepetra", :liblocaepetra),
    LibraryProduct("libteuchosnumerics", :libteuchosnumerics),
    LibraryProduct("libteuchoscomm", :libteuchoscomm),
    LibraryProduct("libtrilinoscouplings", :libtrilinoscouplings),
    LibraryProduct("libepetraext", :libepetraext),
    LibraryProduct("libtriutils", :libtriutils),
    LibraryProduct("libbelosepetra", :libbelosepetra),
    LibraryProduct("libepetra", :libepetra),
    LibraryProduct("libtrilinosss", :libtrilinosss),
    LibraryProduct("libnoxepetra", :libnoxepetra),
    LibraryProduct("libsacado", :libsacado),
    LibraryProduct("libamesos", :libamesos)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="SuiteSparse_jll", uuid="bea87d4a-7f5b-5778-9afe-8cc45184846c"))
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7.1.0")

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
@inline function ExaTronKernels.dmid(n::Int, x::CuDeviceArray{Float64,1},
                      xl::CuDeviceArray{Float64,1}, xu::CuDeviceArray{Float64,1})
    tx = threadIdx().x

    if tx <= n
        @inbounds x[tx] = max(xl[tx], min(x[tx], xu[tx]))
    end
    CUDA.sync_threads()

    return
end


text:
let
    State = TaskGraphs.State
    Action = TaskGraphs.Action

    f = pctapf_problem_1
    cost_model = SumOfMakeSpans()
    solver = NBSSolver(assignment_model = TaskGraphsMILPSolver(GreedyAssignment()))
    pc_tapf = f(solver;cost_function=cost_model,verbose=false)

    base_search_env = pc_tapf.env
    prob = formulate_assignment_problem(solver.assignment_model,pc_tapf)
    sched, cost = solve_assignment_problem!(solver.assignment_model,prob,pc_tapf)
    @test validate(sched)
    search_env = construct_search_env(solver,deepcopy(sched),base_search_env)
    pc_mapf = PC_MAPF(search_env)
    node = initialize_root_node(search_env)
    let
        env = build_env(solver,pc_mapf,search_env,node,AgentID(1))

        s = State()
        a = Action()
        get_next_state(s,a)
        get_next_state(env,s,a)

        @test State(1,2) == State(1,2)
        @test State() == State()
        @test states_match(State(),State())
        @test states_match(env,State(),State())
        CRCBS.wait(s)
        CRCBS.wait(env,s)

        get_cost_model(env)
        get_heuristic_model(env)
        @test CRCBS.is_valid(env,State()) == false
        @test CRCBS.is_valid(env,Action()) == false
        @test CRCBS.is_valid(env,State(1,2))
        @test CRCBS.is_valid(env,Action(Edge(1,2),1))

        get_heuristic_cost(env,State(1,2))
        # @test is_goal(env,State(1,0)) == false
        # @test is_goal(env,State(2,3)) == false
        # @test is_goal(env,State(1,3))
        @test length(get_possible_actions(env,State(-1,-1))) == 0
    end
    let
        for v in vertices(search_env.schedule)
            schedule_node = get_node_from_vtx(search_env.schedule,v)
            if isa(schedule_node,COLLECT)
                env = build_env(solver,pc_mapf,search_env,node,VtxID(v))
                @test length(get_possible_actions(env,State(1,0))) == 1
                @test length(get_possible_actions(env,State(-1,-1))) == 0
                break
            end
        end
    end
end


text:
module FactorModels

using DataArrays
using DataFrames

import GLM.predict

export FactorModel, DynamicFactorModel, predict, calculate_factors, scree_plot,
    lag_vector, lagged_matrix, factor_model_DGP, normalize, pseudo_out_of_sample_forecasts, MSE, RMSE, benchmark_forecasts, benchmark_ar, diebold_mariano, matrix_to_table, print_matrix,
    targeted_predictors,
    wild_bootstrap, residual_bootstrap,
    LR_test, LM_test, LM_test_gls, Wald_test, quandt_andrews, quandt_andrews_critical_value, chow_test

include("FactorModel.jl")
include("DynamicFactorModel.jl")
include("utils.jl")
include("targeted_predictors.jl")
include("bootstrap.jl")
include("chowtest.jl")

end # module


text:
using Documenter, StableDQMC

makedocs(
    modules = [StableDQMC],
    doctest = false,
    sitename = "StableDQMC.jl",
    pages = [
        "Home" => "index.md",
        "Stable Inversion" => "inversions.md",
        "Matrix Decompositions" => "decompositions.md",
        "Exports" => ["exports/qr_udt.md", "exports/svd.md"],
    ],
    # assets = ["assets/custom.css", "assets/custom.js"]
)

deploydocs(
    repo = "github.com/crstnbr/StableDQMC.jl.git",
    push_preview = true,
    # target = "site",
)