text:
function switcheroo(input::AbstractString)::String
    translation = Dict('a'=>'b','b'=>'a','c'=>'c')
    string( map(c->translation[c], input) )
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
@testset "Sigma points" begin
    @testset "Weighting parameters" begin
        num_states = 5

        weight_params = WanMerweWeightingParameters(0.5, 2, 0)
        @test @inferred(KalmanFilters.lambda(weight_params, num_states)) == -3.75
        @test @inferred(KalmanFilters.calc_mean_weights(weight_params, num_states)) == (-3, 0.4)
        @test @inferred(KalmanFilters.calc_cov_weights(weight_params, num_states)) == (-0.25, 0.4)
        @test @inferred(KalmanFilters.calc_cholesky_weight(weight_params, num_states)) == 1.25

        weight_params = MeanSetWeightingParameters(0.5)
        @test @inferred(KalmanFilters.calc_mean_weights(weight_params, num_states)) == (0.5, 0.05)
        @test @inferred(KalmanFilters.calc_cov_weights(weight_params, num_states)) == (0.5, 0.05)
        @test @inferred(KalmanFilters.calc_cholesky_weight(weight_params, num_states)) == 10

        weight_params = GaussSetWeightingParameters(3)
        @test all(@inferred(KalmanFilters.calc_mean_weights(weight_params, num_states)) .≈ (-2/3, 1/6))
        @test all(@inferred(KalmanFilters.calc_cov_weights(weight_params, num_states)) .≈ (-2/3, 1/6))
        @test @inferred(KalmanFilters.calc_cholesky_weight(weight_params, num_states)) == 3

        weight_params = ScaledSetWeightingParameters(0.5, 2, 1)
        @test @inferred(KalmanFilters.calc_mean_weights(weight_params, num_states)) == (-19, 2)
        @test @inferred(KalmanFilters.calc_cov_weights(weight_params, num_states)) == (-16.25, 2)
        @test @inferred(KalmanFilters.calc_cholesky_weight(weight_params, num_states)) == 0.25
    end

    @testset "Create sigma points" begin
        weight_params = ScaledSetWeightingParameters(0.5, 2, 1)
        A = [4.0 0.0; 0.0 4.0]
        x = randn(2)
        χ = @inferred KalmanFilters.calc_sigma_points(x, A, weight_params)
        @test χ.x0 == x
        @test χ.P_chol == cholesky(A .* 0.25).L

        P_chol_temp = zeros(2,2)
        χ = @inferred KalmanFilters.calc_sigma_points!(P_chol_temp, x, A, weight_params)
        @test χ.x0 == x
        @test χ.P_chol == cholesky(A .* 0.25).L
    end

    @testset "Weighted means" begin
        weight_params = ScaledSetWeightingParameters(0.5, 2, 1)
        χ = KalmanFilters.TransformedSigmaPoints(ones(5), hcat(ones(5,5) .* 4, ones(5,5) .* 2), weight_params)
        x = @inferred KalmanFilters.mean(χ)
        @test x == ones(5) .* -19 .+ ones(5) .* 20 .* 2 .+ ones(5) .* 10 .* 2

        χ = KalmanFilters.TransformedSigmaPoints(ones(5), hcat(ones(5,5) .* 4, ones(5,5) .* 2), weight_params)
        x = zeros(5)
        x = @inferred KalmanFilters.mean!(x, χ)
        @test x == ones(5) .* -19 .+ ones(5) .* 20 .* 2 .+ ones(5) .* 10 .* 2
    end

    @testset "substract mean" begin
        weight_params = ScaledSetWeightingParameters(0.5, 2, 1)
        χ = KalmanFilters.TransformedSigmaPoints(ones(5), hcat(ones(5,5) .* 4, ones(5,5) .* 2), weight_params)
        x = ones(5) .* 2
        unbiased_χ = @inferred KalmanFilters.substract_mean(χ, x)
        @test unbiased_χ == [-ones(5) hcat(ones(5,5) .* 2, zeros(5,5))]
        @test unbiased_χ.weight_params == weight_params

        other_weight_params = MeanSetWeightingParameters(0.0)
        unbiased_χ_temp = KalmanFilters.TransformedSigmaPoints(randn(5), randn(5,10), other_weight_params)
        unbiased_χ = @inferred KalmanFilters.substract_mean!(unbiased_χ_temp, χ, x)
        @test unbiased_χ == [-ones(5) hcat(ones(5,5) .* 2, zeros(5,5))]
        @test unbiased_χ.weight_params == weight_params
    end

    @testset "Covariance" begin
        weight_params = ScaledSetWeightingParameters(0.5, 2, 1)
        𝓨 = KalmanFilters.TransformedSigmaPoints(ones(5), hcat(ones(5,5) .* 4, ones(5,5) .* 2), weight_params)
        noise = Diagonal(ones(5))
        P = @inferred KalmanFilters.cov(𝓨, noise)
        @test P == -16.25 .* ones(5) * ones(5)' .+ 2 .* ones(5,5) * ones(5,5)' .* 16 .+
            2 .* ones(5,5) * ones(5,5)' .* 4 .+ noise

        P = zeros(5,5)
        P = @inferred KalmanFilters.cov!(P, 𝓨, noise)
        @test P == -16.25 .* ones(5) * ones(5)' .+ 2 .* ones(5,5) * ones(5,5)' .* 16 .+
            2 .* ones(5,5) * ones(5,5)' .* 4 .+ noise

        P_chol = LowerTriangular(ones(5,5))
        x = ones(5) .* 4
        χ = KalmanFilters.SigmaPoints(x, P_chol, weight_params)
        P = @inferred KalmanFilters.cov(χ, 𝓨)
        @test P == 2 .* LowerTriangular(ones(5,5)) * ones(5,5)' .* 4 .+
            2 .* LowerTriangular(-ones(5,5)) * ones(5,5)' .* 2

        P = zeros(5,5)
        P = @inferred KalmanFilters.cov!(P, χ, 𝓨)
        @test P == 2 .* LowerTriangular(ones(5,5)) * ones(5,5)' .* 4 .+
            2 .* LowerTriangular(-ones(5,5)) * ones(5,5)' .* 2
    end

    @testset "Transform sigma points" begin
        P = [2.0 0.0; 0.0 2.0]
        x = [1.0, 1.0]
        weight_params = ScaledSetWeightingParameters(0.5, 2, 1)
        χ = KalmanFilters.SigmaPoints(x, LowerTriangular(P), weight_params)
        F(x) = x .* 2
        𝓨 = @inferred KalmanFilters.transform(F, χ)
        @test 𝓨 == [ones(2) .* 2 [6 2; 2 6] [-2 2; 2 -2]]
        @test 𝓨.weight_params == weight_params

        other_weight_params = MeanSetWeightingParameters(0.5)
        𝓨_temp = KalmanFilters.TransformedSigmaPoints(zeros(2), zeros(2,4), other_weight_params)
        F!(y, x) = y .= x .* 2
        xi_temp = zeros(length(x))
        𝓨 = @inferred KalmanFilters.transform!(𝓨_temp, xi_temp, F!, χ)
        @test 𝓨 == [ones(2) .* 2 [6 2; 2 6] [-2 2; 2 -2]]
        @test 𝓨.weight_params == weight_params
    end
end

ext: jl


text:
module DPClustering

# packages
using DataFrames
using ProgressMeter
using KernelDensity
using Plots
using Distributions
using Colors

import Base.show

export dpclustering,
DPresults,
plotresults,
show

include("types.jl")
include("plotting.jl")
include("util.jl")
include("clustering.jl")

end # module


text:
struct SampleBased <: ConversionTrait end

function convert_arguments(::SampleBased, args::NTuple{N,AbstractVector{<:Number}}) where {N}
    return args
end

function convert_arguments(P::SampleBased, positions::Vararg{AbstractVector})
    return convert_arguments(P, positions)
end

function convert_arguments(::SampleBased, positions::NTuple{N,AbstractVector}) where {N}
    x = first(positions)
    if any(n-> length(x) != length(n), positions)
        error("all vector need to be same length. Found: $(length.(positions))")
    end
    labels = categoric_labels.(positions)
    xyrange = categoric_range.(labels)
    newpos = map(positions, labels) do pos,lab
        el32convert(categoric_position.(pos, Ref(lab)))
    end
    PlotSpec(newpos...; tickranges = xyrange, ticklabels = labels)
end

max_stars_repo_name: piever/StatMakie

max_issues_repo_name: piever/StatMakie

max_forks_repo_name: piever/StatMakie


text:
t = [0, 0, 0, 1, 1, 1]
x = [-10, -5, -1, 3, 4, 8]

#=
    Model: y = 1 / (1 + e^(-w' x))

    Cost Function =
=#


text:
module GoldenSectionSearch

export golden_section_search

Base.@irrational INV_PHI    0.61803398874989484820  (sqrt(big(5)) - 1) / 2
Base.@irrational INV_PHI_2  0.38196601125010515180  (3 - sqrt(big(5))) / 2

function golden_section_search(f, a, b, n, args...)
    h = b - a
    c = a + INV_PHI_2 * h
    d = a + INV_PHI * h
    yc = f(c, args...)
    yd = f(d, args...)
    for i = 1 : n
        if yc < yd
            b = d
            d = c
            yd = yc
            h *= INV_PHI
            c = a + INV_PHI_2 * h
            yc = f(c, args...)
        else
            a = c
            c = d
            yc = yd
            h *= INV_PHI
            d = a + INV_PHI * h
            yd = f(d, args...)
        end
    end
    if yc < yd
        c, yc
    else
        d, yd
    end
end

end # module GoldenSectionSearch


text:
module Solvers

export
    BatchedTridiagonalSolver, solve_batched_tridiagonal_system!,
    PressureSolver, solve_for_pressure!

using CUDA
using KernelAbstractions
using KernelAbstractions.Extras.LoopInfo: @unroll

using Oceananigans.Architectures: device, @hascuda, CPU, GPU, array_type
using Oceananigans.Utils
using Oceananigans.Grids
using Oceananigans.Grids: unpack_grid

include("solver_utils.jl")
include("batched_tridiagonal_solver.jl")

include("discrete_eigenvalues.jl")
include("plan_transforms.jl")
include("pressure_solver.jl")
include("triply_periodic_pressure_solver.jl")
include("horizontally_periodic_pressure_solver.jl")
include("channel_pressure_solver.jl")
include("box_pressure_solver.jl")
include("index_permutations.jl")
include("solve_for_pressure.jl")

end


text:
using Pkg
#Pkg.activate(dirname(@__FILE__))
Pkg.activate(@__DIR__)
using Fire
using PhysicsTutorials

@main function generate(filename::String)
    pns = _splitpath(abspath(filename))
    length(pns) < 3 && error("Invalid path to source file: $filename")

    category, tutorial_name = pns[end-2:end-1]
    fname, ext = splitext(pns[end])

    fname==tutorial_name || error("tutorial name $tutorial_name is not consistent with filename $fname, expected format: `tutorials/quantum_computing/<name>/<name>.ipynb`")

    if ext == ".ipynb"
        source = PhysicsTutorials.NotebookSource()
    elseif ext ==".jmd"
        source = PhysicsTutorials.WeaveSource()
    elseif ext ==".jl"
        source = PhysicsTutorials.LiterateSource()
    else
        error("Expected file extension: `*.jl`, `*.jmd` or `*.ipynb`, but got $ext.")
    end
    PhysicsTutorials.convert_tutorial(category,tutorial_name,source)
end

# compatibility patch
const path_dir_splitter = r"^(.*?)([/\\]+)([^/\\]*)$"

_splitdir_nodrive(path::String) = _splitdir_nodrive("", path)
function _splitdir_nodrive(a::String, b::String)
    m = match(path_dir_splitter,b)
    m === nothing && return (a,b)
    a = string(a, isempty(m.captures[1]) ? m.captures[2][1] : m.captures[1])
    a, String(m.captures[3])
end

function _splitpath(p::String)
    drive, p = splitdrive(p)
    out = String[]
    isempty(p) && (pushfirst!(out,p))  # "" means the current directory.
    while !isempty(p)
        dir, base = _splitdir_nodrive(p)
        dir == p && (pushfirst!(out, dir); break)  # Reached root node.
        if !isempty(base)  # Skip trailing '/' in basename
            pushfirst!(out, base)
        end
        p = dir
    end
    if !isempty(drive)  # Tack the drive back on to the first element.
        out[1] = drive*out[1]  # Note that length(out) is always >= 1.
    end
    return out
end

max_issues_repo_name: JuliaPhysics/PhysicsTutorials

max_forks_repo_name: JuliaPhysics/PhysicsTutorials


text:
using Mocking

complete_multipart_upload_1_patch = @patch function complete_multipart_upload(a...; b...)
    return Dict{String, Any}("location" => "/111122223333/vaults/my-vault/archives/NkbByEejwEggmBz2fTHgJrg0XBoDfjP4q6iu87-TjhqG6eGoOY9Z8i1_AUyUsuhPAdTqLHy8pTl5nfCFJmDl2yEZONi5L26Omw12vcs01MNGntHEQL8MBfGlqrEXAMPLEArchiveId", "checksum" => "9628195fcdbcbbe76cdde456d4646fa7de5f219fb39823836d81f0cc0e18aa67", "archiveId" => "NkbByEejwEggmBz2fTHgJrg0XBoDfjP4q6iu87-TjhqG6eGoOY9Z8i1_AUyUsuhPAdTqLHy8pTl5nfCFJmDl2yEZONi5L26Omw12vcs01MNGntHEQL8MBfGlqrEXAMPLEArchiveId")
end

create_vault_1_patch = @patch function create_vault(a...; b...)
    return Dict{String, Any}("location" => "/111122223333/vaults/my-vault")
end

describe_job_1_patch = @patch function describe_job(a...; b...)
    return Dict{String, Any}("InventoryRetrievalParameters" => Dict{String, Any}("Format" => "JSON"), "Completed" => false, "StatusCode" => "InProgress", "CreationDate" => "2015-07-17T20:23:41.616Z", "VaultARN" => "arn:aws:glacier:us-west-2:0123456789012:vaults/my-vault", "Action" => "InventoryRetrieval", "JobId" => "zbxcm3Z_3z5UkoroF7SuZKrxgGoDc3RloGduS7Eg-RO47Yc6FxsdGBgf_Q2DK5Ejh18CnTS5XW4_XqlNHS61dsO4CnMW")
end

describe_vault_1_patch = @patch function describe_vault(a...; b...)
    return Dict{String, Any}("VaultName" => "my-vault", "SizeInBytes" => 0, "NumberOfArchives" => 0, "CreationDate" => "2016-09-23T19:27:18.665Z", "VaultARN" => "arn:aws:glacier:us-west-2:111122223333:vaults/my-vault")
end

get_data_retrieval_policy_1_patch = @patch function get_data_retrieval_policy(a...; b...)
    return Dict{String, Any}("Policy" => Dict{String, Any}("Rules" => Any[Dict{String, Any}("Strategy" => "BytesPerHour", "BytesPerHour" => 10737418240)]))
end

get_job_output_1_patch = @patch function get_job_output(a...; b...)
    return Dict{String, Any}("body" => "inventory-data", "status" => 200, "acceptRanges" => "bytes", "contentType" => "application/json")
end

get_vault_access_policy_1_patch = @patch function get_vault_access_policy(a...; b...)
    return Dict{String, Any}("policy" => Dict{String, Any}("Policy" => "{\"Version\":\"2012-10-17\",\"Statement\":[{\"Sid\":\"Define-owner-access-rights\",\"Effect\":\"Allow\",\"Principal\":{\"AWS\":\"arn:aws:iam::999999999999:root\"},\"Action\":\"glacier:DeleteArchive\",\"Resource\":\"arn:aws:glacier:us-west-2:999999999999:vaults/examplevault\"}]}"))
end

get_vault_lock_1_patch = @patch function get_vault_lock(a...; b...)
    return Dict{String, Any}("ExpirationDate" => "exampledate", "CreationDate" => "exampledate", "Policy" => "{\"Version\":\"2012-10-17\",\"Statement\":[{\"Sid\":\"Define-vault-lock\",\"Effect\":\"Deny\",\"Principal\":{\"AWS\":\"arn:aws:iam::999999999999:root\"},\"Action\":\"glacier:DeleteArchive\",\"Resource\":\"arn:aws:glacier:us-west-2:999999999999:vaults/examplevault\",\"Condition\":{\"NumericLessThanEquals\":{\"glacier:ArchiveAgeinDays\":\"365\"}}}]}", "State" => "InProgress")
end

get_vault_notifications_1_patch = @patch function get_vault_notifications(a...; b...)
    return Dict{String, Any}("vaultNotificationConfig" => Dict{String, Any}("SNSTopic" => "arn:aws:sns:us-west-2:0123456789012:my-vault", "Events" => Any["InventoryRetrievalCompleted", "ArchiveRetrievalCompleted"]))
end

initiate_job_1_patch = @patch function initiate_job(a...; b...)
    return Dict{String, Any}("jobId" => " HkF9p6o7yjhFx-K3CGl6fuSm6VzW9T7esGQfco8nUXVYwS0jlb5gq1JZ55yHgt5vP54ZShjoQzQVVh7vEXAMPLEjobID", "location" => "/111122223333/vaults/examplevault/jobs/HkF9p6o7yjhFx-K3CGl6fuSm6VzW9T7esGQfco8nUXVYwS0jlb5gq1JZ55yHgt5vP54ZShjoQzQVVh7vEXAMPLEjobID")
end

initiate_multipart_upload_1_patch = @patch function initiate_multipart_upload(a...; b...)
    return Dict{String, Any}("location" => "/111122223333/vaults/my-vault/multipart-uploads/19gaRezEXAMPLES6Ry5YYdqthHOC_kGRCT03L9yetr220UmPtBYKk-OssZtLqyFu7sY1_lR7vgFuJV6NtcV5zpsJ", "uploadId" => "19gaRezEXAMPLES6Ry5YYdqthHOC_kGRCT03L9yetr220UmPtBYKk-OssZtLqyFu7sY1_lR7vgFuJV6NtcV5zpsJ")
end

initiate_vault_lock_1_patch = @patch function initiate_vault_lock(a...; b...)
    return Dict{String, Any}("lockId" => "AE863rKkWZU53SLW5be4DUcW")
end

list_jobs_1_patch = @patch function list_jobs(a...; b...)
    return Dict{String, Any}("JobList" => Any[Dict{String, Any}("RetrievalByteRange" => "0-3145727", "Action" => "ArchiveRetrieval", "SNSTopic" => "arn:aws:sns:us-west-2:0123456789012:my-vault", "ArchiveSizeInBytes" => 3145728, "VaultARN" => "arn:aws:glacier:us-west-2:0123456789012:vaults/my-vault", "StatusCode" => "InProgress", "CreationDate" => "2015-07-17T21:16:13.840Z", "JobDescription" => "Retrieve archive on 2015-07-17", "JobId" => "l7IL5-EkXyEY9Ws95fClzIbk2O5uLYaFdAYOi-azsX_Z8V6NH4yERHzars8wTKYQMX6nBDI9cMNHzyZJO59-8N9aHWav", "ArchiveId" => "kKB7ymWJVpPSwhGP6ycSOAekp9ZYe_--zM_mw6k76ZFGEIWQX-ybtRDvc2VkPSDtfKmQrj0IRQLSGsNuDp-AJVlu2ccmDSyDUmZwKbwbpAdGATGDiB3hHO0bjbGehXTcApVud_wyDw", "Completed" => false, "ArchiveSHA256TreeHash" => "9628195fcdbcbbe76cdde932d4646fa7de5f219fb39823836d81f0cc0e18aa67", "SHA256TreeHash" => "9628195fcdbcbbe76cdde932d4646fa7de5f219fb39823836d81f0cc0e18aa67"), Dict{String, Any}("InventoryRetrievalParameters" => Dict{String, Any}("Format" => "JSON"), "Completed" => false, "StatusCode" => "InProgress", "CreationDate" => "2015-07-17T20:23:41.616Z", "VaultARN" => "arn:aws:glacier:us-west-2:0123456789012:vaults/my-vault", "Action" => "InventoryRetrieval", "JobId" => "zbxcm3Z_3z5UkoroF7SuZKrxgGoDc3RloGduS7Eg-RO47Yc6FxsdGBgf_Q2DK5Ejh18CnTS5XW4_XqlNHS61dsO4CnMW")])
end

list_multipart_uploads_1_patch = @patch function list_multipart_uploads(a...; b...)
    return Dict{String, Any}("Marker" => "null", "UploadsList" => Any[Dict{String, Any}("PartSizeInBytes" => 4194304, "MultipartUploadId" => "xsQdFIRsfJr20CW2AbZBKpRZAFTZSJIMtL2hYf8mvp8dM0m4RUzlaqoEye6g3h3ecqB_zqwB7zLDMeSWhwo65re4C4Ev", "ArchiveDescription" => "archive 1", "CreationDate" => "2012-03-19T23:20:59.130Z", "VaultARN" => "arn:aws:glacier:us-west-2:012345678901:vaults/examplevault"), Dict{String, Any}("PartSizeInBytes" => 4194304, "MultipartUploadId" => "nPyGOnyFcx67qqX7E-0tSGiRi88hHMOwOxR-_jNyM6RjVMFfV29lFqZ3rNsSaWBugg6OP92pRtufeHdQH7ClIpSF6uJc", "ArchiveDescription" => "archive 2", "CreationDate" => "2012-04-01T15:00:00.000Z", "VaultARN" => "arn:aws:glacier:us-west-2:012345678901:vaults/examplevault"), Dict{String, Any}("PartSizeInBytes" => 4194304, "MultipartUploadId" => "qt-RBst_7yO8gVIonIBsAxr2t-db0pE4s8MNeGjKjGdNpuU-cdSAcqG62guwV9r5jh5mLyFPzFEitTpNE7iQfHiu1XoV", "ArchiveDescription" => "archive 3", "CreationDate" => "2012-03-20T17:03:43.221Z", "VaultARN" => "arn:aws:glacier:us-west-2:012345678901:vaults/examplevault")])
end

list_parts_1_patch = @patch function list_parts(a...; b...)
    return Dict{String, Any}("PartSizeInBytes" => 4194304, "Marker" => "null", "MultipartUploadId" => "OW2fM5iVylEpFEMM9_HpKowRapC3vn5sSL39_396UW9zLFUWVrnRHaPjUJddQ5OxSHVXjYtrN47NBZ-khxOjyEXAMPLE", "ArchiveDescription" => "archive description", "CreationDate" => "2012-03-20T17:03:43.221Z", "VaultARN" => "arn:aws:glacier:us-west-2:012345678901:vaults/demo1-vault", "Parts" => Any[Dict{String, Any}("RangeInBytes" => "0-4194303", "SHA256TreeHash" => "01d34dabf7be316472c93b1ef80721f5d4"), Dict{String, Any}("RangeInBytes" => "4194304-8388607", "SHA256TreeHash" => "0195875365afda349fc21c84c099987164")])
end

list_provisioned_capacity_1_patch = @patch function list_provisioned_capacity(a...; b...)
    return Dict{String, Any}("ProvisionedCapacityList" => Any[Dict{String, Any}("ExpirationDate" => "2016-12-12T00:00:00.000Z", "StartDate" => "2016-11-11T20:11:51.095Z", "CapacityId" => "zSaq7NzHFQDANTfQkDen4V7z"), Dict{String, Any}("ExpirationDate" => "2017-01-15T00:00:00.000Z", "StartDate" => "2016-12-13T20:11:51.095Z", "CapacityId" => "yXaq7NzHFQNADTfQkDen4V7z")])
end

list_tags_for_vault_1_patch = @patch function list_tags_for_vault(a...; b...)
    return Dict{String, Any}("Tags" => Dict{String, Any}("id" => "1234", "date" => "july2015"))
end

list_vaults_1_patch = @patch function list_vaults(a...; b...)
    return Dict{String, Any}("VaultList" => Any[Dict{String, Any}("VaultName" => "my-vault", "SizeInBytes" => 3178496, "NumberOfArchives" => 1, "CreationDate" => "2015-04-06T21:23:45.708Z", "VaultARN" => "arn:aws:glacier:us-west-2:0123456789012:vaults/my-vault", "LastInventoryDate" => "2015-04-07T00:26:19.028Z")])
end

purchase_provisioned_capacity_1_patch = @patch function purchase_provisioned_capacity(a...; b...)
    return Dict{String, Any}("capacityId" => "zSaq7NzHFQDANTfQkDen4V7z")
end

upload_archive_1_patch = @patch function upload_archive(a...; b...)
    return Dict{String, Any}("location" => "/0123456789012/vaults/my-vault/archives/kKB7ymWJVpPSwhGP6ycSOAekp9ZYe_--zM_mw6k76ZFGEIWQX-ybtRDvc2VkPSDtfKmQrj0IRQLSGsNuDp-AJVlu2ccmDSyDUmZwKbwbpAdGATGDiB3hHO0bjbGehXTcApVud_wyDw", "checksum" => "969fb39823836d81f0cc028195fcdbcbbe76cdde932d4646fa7de5f21e18aa67", "archiveId" => "kKB7ymWJVpPSwhGP6ycSOAekp9ZYe_--zM_mw6k76ZFGEIWQX-ybtRDvc2VkPSDtfKmQrj0IRQLSGsNuDp-AJVlu2ccmDSyDUmZwKbwbpAdGATGDiB3hHO0bjbGehXTcApVud_wyDw")
end

upload_multipart_part_1_patch = @patch function upload_multipart_part(a...; b...)
    return Dict{String, Any}("checksum" => "c06f7cd4baacb087002a99a5f48bf953")
end


text:
using ImgTools
using Test

@testset "ImgTools.jl" begin

    # single plot
    p = ImgTools.Plots.plot(rand(100); label = "")
    figfile = ImgTools.sfig(p,
        [tempdir()], "test", (;A = 1), "png"
    )
    @test isfile(figfile)

    # miltiple plots
    ps = map(1:10) do _
        p = ImgTools.Plots.plot(rand(100); label = "")
    end
    figfile = ImgTools.sfig(ps,
        [tempdir()], "test", (;A = 1), "png"
    )
    @test isfile(figfile)

    # gifs
    figfile = ImgTools.sgif(ps,
        [tempdir()], "test", (;A = 1), "gif"
    )
    @test isfile(figfile)

    rm(figfile)

end


text:
using DSP

import Base: getindex, size, length, full, issparse

@compat abstract type Frequency end

nfreqs(f::Frequency) = error("Not implemented")
freqs(f::Frequency) = error("Not implemented")

# Geometrically spaced frequency
# fₖ = min * 2^(1/bins)ᵏ
immutable GeometricFrequency <: Frequency
    min::Real
    max::Real
    bins::Real # the number of bins per octave

    function GeometricFrequency(min, max, bins=24)
        min <= max || throw(ArgumentError("max must be larger than min"))
        bins > 0 || throw(ArgumentError("bins must be positive"))
        new(min, max, bins)
    end
end

function GeometricFrequency(; min::Real=55, max::Real=5000, bins::Real=24)
    GeometricFrequency(min, max, bins)
end

nbins_per_octave(freq::GeometricFrequency) = freq.bins

# Q-factor
q(bins, qrate=1.0) = 1.0 / (2 ^ (1.0/bins) - 1) * qrate
q(f::GeometricFrequency, qrate=1.0) = q(f.bins, qrate)

function nfreqs(f::GeometricFrequency)
    convert(Int, round(nbins_per_octave(f)  * (log2(f.max / f.min))))
end

function freqs(f::GeometricFrequency)
    f.min * 2 .^ ((0:nfreqs(f)-1) / nbins_per_octave(f))
end

# Kernel property
immutable KernelProperty
    fs::Real
    freq::GeometricFrequency
    win::Function

    function KernelProperty(fs, freq, win)
        applicable(win, 512) || error("Invalid window function")
        new(fs, freq, win)
    end
end

# Kernel matrices
@compat abstract type KernelMatrix{T} <: AbstractMatrix{T} end

# Spectral kernel (freqency-domain kernel)
immutable SpectralKernelMatrix{T<:Complex} <: KernelMatrix{T}
    data::AbstractMatrix{T}
    property::KernelProperty
end

# Temporal kernel
immutable TemporalKernelMatrix{T<:Complex} <: KernelMatrix{T}
    data::AbstractMatrix{T}
    property::KernelProperty
end

property(k::KernelMatrix) = k.property
rawdata(k::KernelMatrix) = k.data

size(k::KernelMatrix) = size(k.data)
length(k::KernelMatrix) = length(k.data)

getindex(k::KernelMatrix, i::Real) = getindex(k.data, i)
getindex(k::KernelMatrix, i::Real...) = getindex(k.data, i...)
getindex(k::KernelMatrix, i::Real, j::Real) = getindex(k.data, i, j)
getindex(k::KernelMatrix, i::Range, j::Real) = getindex(k.data, i, j)
getindex(k::KernelMatrix, i::Real, j::Range) = getindex(k.data, i, j)

issparse(k::KernelMatrix) = issparse(k.data)
full(k::KernelMatrix) = full(k.data)

function tempkernel(T::Type,
                   fs::Real,
                   freq::GeometricFrequency=GeometricFrequency(55, fs/2),
                   win::Function=hamming)
    prop = KernelProperty(fs, freq, win)
    data = _tempkernel(T, prop)
    TemporalKernelMatrix(data, prop)
end

function _tempkernel(T::Type, prop::KernelProperty)
    _tempkernel(T, prop.fs, prop.freq, prop.win)
end

function _tempkernel(T::Type, fs::Real, freq::GeometricFrequency, win::Function)
    Q = q(freq)
    f = freqs(freq)
    winsizes = map(Int, map(trunc, (fs ./ f * Q)))
    fftlen = nextpow2(winsizes[1])

    K = zeros(Complex{T}, fftlen, length(winsizes))

    for k = 1:length(winsizes)
        Nk = winsizes[k]
        kernel = win(Nk) .* exp.(-im*2π*Q/Nk .* (1:Nk)) / Nk
        s = (fftlen - Nk) >> 1
        copy!(K, fftlen*(k-1) + s, kernel, 1, Nk)
    end

    K
end

function speckernel(T::Type,
                    fs::Real,
                    freq::GeometricFrequency=GeometricFrequency(55, fs/2),
                    win::Function=hamming,
                    threshold=0.005)
    prop = KernelProperty(fs, freq, win)
    data = _speckernel(T, prop, threshold)
    SpectralKernelMatrix(data, prop)
end

function _speckernel(T::Type, prop::KernelProperty, threshold=0.005)
    _speckernel(T, prop.fs, prop.freq, prop.win, threshold)
end

# Compute sparse kernel matrix in frequency-domain
function _speckernel(T::Type,
                     fs::Real,
                     freq::GeometricFrequency,
                     win::Function,
                     threshold)
    K = _tempkernel(T, fs, freq, win)
    conj!(K)

    fftlen = size(K, 1)

    # to frequency domain
    fft!(K, 1)

    # make it sparse
    K[abs.(K) .< threshold] = 0.0
    Kˢ = sparse(K)

    # take complex conjugate
    conj!(Kˢ)

    # normalize by fftlen
    Kˢ ./= fftlen

    Kˢ
end

function sym!(symfftout, fftout)
    copy!(symfftout, 1, fftout, 1, length(fftout))
    @inbounds for i = 1:length(fftout)-1
        symfftout[end-i+1] = conj(fftout[i+1])
    end
end

# J. C. Brown and M. S. Puckette, "An efficient algorithm for the calculation
# of a constant Q transform," J. Acoust. Soc. Amer., vol. 92, no. 5,
# pp. 2698–2701, 1992.
function cqt{T}(x::Vector{T},
                fs::Real,
                freq::GeometricFrequency=GeometricFrequency(55, fs/2),
                hopsize::Int=convert(Int, round(fs*0.005)),
                win::Function=hamming,
                K::AbstractMatrix = _speckernel(T, fs, freq, win, 0.005)
                )
    Q = q(freq)
    freqaxis = freqs(freq)

    winsizes = map(Int, map(trunc, (fs ./ freqaxis * Q)))
    nframes = div(length(x), hopsize)

    fftlen = nextpow2(winsizes[1])
    size(K) == (fftlen, length(winsizes)) || error("inconsistent kernel size")

    # Create padded signal
    xpadd = zeros(T, length(x) + fftlen)
    copy!(xpadd, fftlen>>1+1, x, 1, length(x))

    # FFT workspace
    fftin = Vector{T}(fftlen)
    fftout = Vector{Complex{T}}(fftlen>>1+1)
    fplan = plan_rfft(fftin)
    symfftout = Vector{Complex{T}}(fftlen)

    # Constant-Q spectrogram
    X = Array{Complex{T},2}(length(freqaxis), nframes)

    # tmp used in loop
    freqbins = Vector{Complex{T}}(length(freqaxis))

    for n = 1:nframes
        s = hopsize * (n - 1) + 1
        # copy to input buffer
        copy!(fftin, 1, xpadd, s, fftlen)
        # FFT
        A_mul_B!(fftout, fplan, fftin)
        # get full fft bins (rest half are complex conjugate)
        sym!(symfftout, fftout)
        # multiply in frequency-domain
        At_mul_B!(freqbins, K, symfftout)
        # copy to output buffer
        copy!(X, length(freqaxis)*(n-1) + 1, freqbins, 1, length(freqaxis))
    end

    timeaxis = [0:hopsize:nframes;]/fs
    X, timeaxis, freqaxis
end

# CQT with frequency-domain kernel matrix
function cqt(x::Vector,
             fs::Real,
             K::SpectralKernelMatrix;
             hopsize::Int = convert(Int, round(fs * 0.005))
             )
    prop = property(K)
    fs == prop.fs || error("Inconsistent kernel")
    cqt(x, prop.fs, prop.freq, hopsize, prop.win, K.data)
end

# J. C. Brown. "Calculation of a constant Q spectral transform,"
# Journal of the Acoustical Society of America,, 89(1):
# 425–434, 1991.
function time_domain_cqt{T}(x::Vector{T},
                            fs::Real,
                            freq::GeometricFrequency=GeometricFrequency(55, fs/2),
                            hopsize::Int=convert(Int, round(fs*0.005)),
                            win::Function=hamming,
                            K::AbstractMatrix = _tempkernel(T, fs, freq, win)
                            )
    Q = q(freq)
    freqaxis = freqs(freq)

    winsizes = map(Int, map(trunc, (fs ./ freqaxis * Q)))
    nframes = div(length(x), hopsize)

    fftlen = nextpow2(winsizes[1])
    size(K) == (fftlen, length(winsizes)) || error("inconsistent kernel size")

    # Create padded signal
    xpadd = zeros(T, length(x) + fftlen)
    copy!(xpadd, fftlen>>1+1, x, 1, length(x))

    # Constant-q spectrogram
    X = Array{Complex{T},2}(length(freqaxis), nframes)

    # buffer used in roop
    block = Vector{T}(fftlen)
    freqbins = Vector{Complex{T}}(size(X, 1))

    for n = 1:nframes
        s = hopsize * (n - 1) + 1
        # copy to input buffer
        copy!(block, 1, xpadd, s, fftlen)
        # multiply in time-domain
        At_mul_B!(freqbins, K, block)
        # copy to output buffer
        copy!(X, length(freqaxis)*(n-1) + 1, freqbins, 1, length(freqaxis))
    end

    timeaxis = [0:hopsize:nframes;]/fs
    X, timeaxis, freqaxis
end

# CQT with a time-domain kernel matrix
function cqt(x::Vector,
             fs::Real,
             K::TemporalKernelMatrix;
             hopsize::Int = convert(Int, round(fs * 0.005))
             )
    prop = property(K)
    fs == prop.fs || error("Inconsistent kernel")
    time_domain_cqt(x, fs, prop.freq, hopsize, prop.win, rawdata(K))
end


text:
struct PendulumParams{T} <:Any where {T<:Real}
    m::T 	# mass link
    l::T 	# length link
	mu::T    # friction constant
	g::T   	# gravity force (not directional)
	# fmag::T # max force magnitude

    PendulumParams() = new{Float64}(1., 1., 0.1, 9.8)
    PendulumParams(T::Type) = new{T}(1., 1., 0.1, 9.8)
	PendulumParams(T::Type, m, l, mu, g) = new{T}(m, l, mu, g)
end

function random_pendulum_params()
	m = 1.0
	l = 1.0
	g = 9.8
	mu = 0.1
	m = rand(Uniform(0.9*m, 1.1*m))  # mass of arm
	l = rand(Uniform(0.9*l, 1.1*l))    # length of arm
	mu = rand(Uniform(0.9*mu, 1.1*mu))  # friction coeff
	g = rand(Uniform(8., 10.))  # gravity
	params = PendulumParams(Float64, m, l, mu, g)
	return params
end

function pendulum_finitetime(;randomize=false, maxT=20.0, dt=0.01, droptime=true, stochastic_start=true)
	if randomize
		params = PendulumParams()
	else
		params = random_pendulum_params()
	end
	return create_pendulum_finitetime(params, maxT=maxT, dt=dt, droptime=droptime, stochastic_start=stochastic_start)
end

function create_pendulum_finitetime(params; maxT=20.0, dt=0.01, droptime=true, stochastic_start=true)
	X = zeros((2,2))
	X[1,:] .= [-π, π]           # theta1 range
	X[2,:] .= [-3. * π, 3. * π] # theta1Dot range
	S = ([0.0, maxT], X)

	A = [-1.0, 1.0]
	sim_steps = 10
	if droptime
		function get_outcome1(s,a,params,dt,sim_steps,maxT,stochastic_start)
			t, x, r, γ = pendulum_step!(s,a, params, dt, sim_steps, maxT, stochastic_start)
			s = (t,x)
			return s,x,r,γ
		end
        p = (s,a)->get_outcome1(s,a,params,dt,sim_steps,maxT,stochastic_start)

        function pdd0_obs()
            t,x = pendulum_sample_initial(stochastic_start)
            return (t,x), x
        end
        d0 = pdd0_obs
	else
		X = S
		function get_outcome2(s,a,params,dt,sim_steps, maxT,stochastic_start)
			t, x, r, γ = pendulum_step!(s,a, params, dt, sim_steps, maxT, stochastic_start)
			s = (t,x)
			return s,s,r,γ
		end
        p = (s,a)->get_outcome2(s,a,params,dt,sim_steps,maxT,stochastic_start)
        function pdd0()
            s = pendlum_sample_initial(stochastic_start)
            return s,s
        end
		d0 = pdd0
    end

	meta = Dict{Symbol,Any}()
    meta[:minreward] = -2.0
    meta[:maxreward] = 2.0
    meta[:minreturn] = -2 * maxT
    meta[:maxreturn] = 2 * maxT
    meta[:stochastic] = false
    meta[:minhorizon] = ceil(Int, maxT / dt)
    meta[:maxhorizon] = ceil(Int, maxT / dt)
	meta[:discounted] = false
	meta[:episodes] = 1000

	render = (state,clearplot=false)->pendulumplot(state, params)

	m = SequentialProblem(S,X,A,p,d0,meta,render)
end

function pendulum_sim!(state, u, params::PendulumParams, dt, simSteps::Int)
	theta, omega = state
	m, l, mu, g = params.m, params.l, params.mu, params.g

	thetaDot = 0.
	omegaDot = 0.
	subDt = dt / simSteps

	for i in 1:simSteps
		thetaDot = omega
		omegaDot = (-mu * omega - m*g*l*sin(theta) + u) / (m*l^2)
		theta += subDt*thetaDot
		omega += subDt*omegaDot

	end

	theta = mod(theta + π, 2. * π) - π
	omega = clamp(omega, -3. * π, 3. * π)

	state .= theta, omega
	return nothing
end

function pendulum_step!(state, action, params::PendulumParams, dt, simSteps, maxT, stochastic_start)
	u = clamp(action, -1., 1.) * 10.0
	t,x = state
	pendulum_sim!(x, u, params, dt, simSteps)
	t += dt

	reward = -cos(x[1]) - (x[2]^2 / 100.)
	γ = 1.0
	done = t ≥ maxT

	if done
        γ = 0.0
        t, x = pendulum_sample_initial(stochastic_start)
	end

	return t, x, reward, γ
end

function pendulum_sample_initial(stochastic_start=true)
	state = zeros(2)
	if stochastic_start
		state[1] = rand(Uniform(-π, π))
	end

	return 0.0, state
end

@userplot PendulumPlot
@recipe function f(ap::PendulumPlot)
    state, params = ap.args
    t,x = state
    theta = x[1]
	l = params.l
	maxlen = (l)*1.2
	x1,y1 = 0., 0.
	x2 = x1 - sin(theta)*l
	y2 = y1 - cos(theta)*l

	legend := false
	xlims := (-maxlen, maxlen)
	ylims := (-maxlen, maxlen)
	grid := false
	ticks := nothing
	foreground_color := :white
	aspect_ratio := 1.

	# target line
	@series begin
		linecolor := :red
		linewidth := 5

		[0., 0.], [0., maxlen]
	end

	# arms
	@series begin
		linecolor := :black
		linewidth := 10

		[x1, x2], [y1, y2]
	end

end


text:
using GeometricAlgebra
using BenchmarkTools

@basis "+++"

const mvec = Multivector(randn(2^N))
const b = 1.0v12
const kvec = fill(1.1, Bivector)
const store = zero(typeof(mvec))

@benchmark Ref($mvec * $b)[]
@benchmark Ref($mvec * $kvec)[]
@benchmark Ref($mvec * $mvec)[]

@benchmark Ref(geom!($store, $mvec, $b))[]
@benchmark Ref(geom!($store, $kvec, $kvec))[]
@benchmark Ref(geom!($store, $mvec, $kvec))[]
@benchmark Ref(geom!($store, $mvec, $mvec))[]

@code_warntype mvec * b
@code_warntype mvec * kvec

@benchmark reverse($mvec)


text:
# sparse linear algebra functions that perform operations between sparse matrices and dense
# vectors

using LinearAlgebra: AbstractTriangular

export sv2!, sv2, sv

"""
    mv!(transa::SparseChar, alpha::BlasFloat, A::CuSparseMatrix, X::CuVector, beta::BlasFloat, Y::CuVector, index::SparseChar)

Performs `Y = alpha * op(A) *X + beta * Y`, where `op` can be nothing (`transa = N`),
tranpose (`transa = T`) or conjugate transpose (`transa = C`). `X` is a sparse vector, and
`Y` is dense.
"""
mv!(transa::SparseChar, alpha::BlasFloat, A::CuSparseMatrix, X::CuVector, beta::BlasFloat, Y::CuVector, index::SparseChar)
for (fname,elty) in ((:cusparseSbsrmv, :Float32),
                     (:cusparseDbsrmv, :Float64),
                     (:cusparseCbsrmv, :ComplexF32),
                     (:cusparseZbsrmv, :ComplexF64))
    @eval begin
        function mv!(transa::SparseChar,
                     alpha::Number,
                     A::CuSparseMatrixBSR{$elty},
                     X::CuVector{$elty},
                     beta::Number,
                     Y::CuVector{$elty},
                     index::SparseChar)
            desc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            m,n = A.dims
            mb = div(m,A.blockDim)
            nb = div(n,A.blockDim)
            if transa == 'N'
                chkmvdims(X,n,Y,m)
            end
            if transa == 'T' || transa == 'C'
                chkmvdims(X,m,Y,n)
            end
            $fname(handle(), A.dir, transa, mb, nb,
                   nnz(A), alpha, desc, nonzeros(A), A.rowPtr,
                   A.colVal, A.blockDim, X, beta, Y)
            Y
        end
    end
end

"""
    sv2!(transa::SparseChar, uplo::SparseChar, alpha::BlasFloat, A::CuSparseMatrixBSR, X::CuVector, index::SparseChar)

Performs `X = alpha * op(A) \\ X `, where `op` can be nothing (`transa = N`), tranpose
(`transa = T`) or conjugate transpose (`transa = C`). `X` is a dense vector, and `uplo`
tells `sv2!` which triangle of the block sparse matrix `A` to reference.
"""
sv2!(transa::SparseChar, uplo::SparseChar, alpha::BlasFloat, A::CuSparseMatrixBSR, X::CuVector, index::SparseChar)
# bsrsv2
for (bname,aname,sname,elty) in ((:cusparseSbsrsv2_bufferSize, :cusparseSbsrsv2_analysis, :cusparseSbsrsv2_solve, :Float32),
                                 (:cusparseDbsrsv2_bufferSize, :cusparseDbsrsv2_analysis, :cusparseDbsrsv2_solve, :Float64),
                                 (:cusparseCbsrsv2_bufferSize, :cusparseCbsrsv2_analysis, :cusparseCbsrsv2_solve, :ComplexF32),
                                 (:cusparseZbsrsv2_bufferSize, :cusparseZbsrsv2_analysis, :cusparseZbsrsv2_solve, :ComplexF64))
    @eval begin
        function sv2!(transa::SparseChar,
                      uplo::SparseChar,
                      alpha::Number,
                      A::CuSparseMatrixBSR{$elty},
                      X::CuVector{$elty},
                      index::SparseChar)
            desc   = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, uplo, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            m,n      = A.dims
            if m != n
                throw(DimensionMismatch("A must be square, but has dimensions ($m,$n)!"))
            end
            mb = div(m,A.blockDim)
            mX = length(X)
            if mX != m
                throw(DimensionMismatch("X must have length $m, but has length $mX"))
            end
            info = bsrsv2Info_t[0]
            cusparseCreateBsrsv2Info(info)
            @workspace size=@argout(
                    $bname(handle(), A.dir, transa, mb, nnz(A),
                           desc, nonzeros(A), A.rowPtr, A.colVal, A.blockDim,
                           info[1], out(Ref{Cint}(1)))
                )[] buffer->begin
                    $aname(handle(), A.dir, transa, mb, nnz(A),
                           desc, nonzeros(A), A.rowPtr, A.colVal, A.blockDim,
                           info[1], CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                    posit = Ref{Cint}(1)
                    cusparseXbsrsv2_zeroPivot(handle(), info[1], posit)
                    if posit[] >= 0
                        error("Structural/numerical zero in A at ($(posit[]),$(posit[])))")
                    end
                    $sname(handle(), A.dir, transa, mb, nnz(A),
                           alpha, desc, nonzeros(A), A.rowPtr, A.colVal,
                           A.blockDim, info[1], X, X,
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                end
            cusparseDestroyBsrsv2Info(info[1])
            X
        end
    end
end

for elty in (:Float32, :Float64, :ComplexF32, :ComplexF64)
    @eval begin
        function sv2(transa::SparseChar,
                     uplo::SparseChar,
                     alpha::Number,
                     A::CuSparseMatrix{$elty},
                     X::CuVector{$elty},
                     index::SparseChar)
            sv2!(transa,uplo,alpha,A,copy(X),index)
        end
        function sv2(transa::SparseChar,
                     uplo::SparseChar,
                     A::CuSparseMatrix{$elty},
                     X::CuVector{$elty},
                     index::SparseChar)
            sv2!(transa,uplo,one($elty),A,copy(X),index)
        end
        function sv2(transa::SparseChar,
                     alpha::Number,
                     A::AbstractTriangular,
                     X::CuVector{$elty},
                     index::SparseChar)
            uplo = 'U'
            if istril(A)
                uplo = 'L'
            end
            sv2!(transa,uplo,alpha,A.data,copy(X),index)
        end
        function sv2(transa::SparseChar,
                     A::AbstractTriangular,
                     X::CuVector{$elty},
                     index::SparseChar)
            uplo = 'U'
            if istril(A)
                uplo = 'L'
            end
            sv2!(transa,uplo,one($elty),A.data,copy(X),index)
        end
    end
end

# csrsv2
for (bname,aname,sname,elty) in ((:cusparseScsrsv2_bufferSize, :cusparseScsrsv2_analysis, :cusparseScsrsv2_solve, :Float32),
                                 (:cusparseDcsrsv2_bufferSize, :cusparseDcsrsv2_analysis, :cusparseDcsrsv2_solve, :Float64),
                                 (:cusparseCcsrsv2_bufferSize, :cusparseCcsrsv2_analysis, :cusparseCcsrsv2_solve, :ComplexF32),
                                 (:cusparseZcsrsv2_bufferSize, :cusparseZcsrsv2_analysis, :cusparseZcsrsv2_solve, :ComplexF64))
    @eval begin
        function sv2!(transa::SparseChar,
                      uplo::SparseChar,
                      alpha::Number,
                      A::CuSparseMatrixCSR{$elty},
                      X::CuVector{$elty},
                      index::SparseChar)
            desc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, uplo, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            m,n = A.dims
            if m != n
                throw(DimensionMismatch("A must be square, but has dimensions ($m,$n)!"))
            end
            mX = length(X)
            if mX != m
                throw(DimensionMismatch("First dimension of A, $m, and of X, $(size(X)[1]) must match"))
            end
            info = csrsv2Info_t[0]
            cusparseCreateCsrsv2Info(info)
            @workspace size=@argout(
                    $bname(handle(), transa, m, nnz(A),
                           desc, nonzeros(A), A.rowPtr, A.colVal, info[1],
                           out(Ref{Cint}(1)))
                )[] buffer->begin
                    $aname(handle(), transa, m, nnz(A),
                           desc, nonzeros(A), A.rowPtr, A.colVal, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                    posit = Ref{Cint}(1)
                    cusparseXcsrsv2_zeroPivot(handle(), info[1], posit)
                    if posit[] >= 0
                        error("Structural/numerical zero in A at ($(posit[]),$(posit[])))")
                    end
                    $sname(handle(), transa, m,
                           nnz(A), alpha, desc, nonzeros(A), A.rowPtr,
                           A.colVal, info[1], X, X,
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                end
            cusparseDestroyCsrsv2Info(info[1])
            X
        end
    end
end

# cscsv2
for (bname,aname,sname,elty) in ((:cusparseScsrsv2_bufferSize, :cusparseScsrsv2_analysis, :cusparseScsrsv2_solve, :Float32),
                                 (:cusparseDcsrsv2_bufferSize, :cusparseDcsrsv2_analysis, :cusparseDcsrsv2_solve, :Float64),
                                 (:cusparseCcsrsv2_bufferSize, :cusparseCcsrsv2_analysis, :cusparseCcsrsv2_solve, :ComplexF32),
                                 (:cusparseZcsrsv2_bufferSize, :cusparseZcsrsv2_analysis, :cusparseZcsrsv2_solve, :ComplexF64))
    @eval begin
        function sv2!(transa::SparseChar,
                      uplo::SparseChar,
                      alpha::Number,
                      A::CuSparseMatrixCSC{$elty},
                      X::CuVector{$elty},
                      index::SparseChar)
            ctransa = 'N'
            cuplo = 'U'
            if transa == 'N'
                ctransa = 'T'
            end
            if uplo == 'U'
                cuplo = 'L'
            end
            desc   = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, cuplo, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            n,m      = A.dims
            if m != n
                throw(DimensionMismatch("A must be square, but has dimensions ($m,$n)!"))
            end
            mX = length(X)
            if mX != m
                throw(DimensionMismatch("First dimension of A, $m, and of X, $(size(X)[1]) must match"))
            end
            info = csrsv2Info_t[0]
            cusparseCreateCsrsv2Info(info)
            @workspace size=@argout(
                    $bname(handle(), ctransa, m, nnz(A),
                           desc, nonzeros(A), A.colPtr, rowvals(A), info[1],
                           out(Ref{Cint}(1)))
                )[] buffer->begin
                    $aname(handle(), ctransa, m, nnz(A),
                           desc, nonzeros(A), A.colPtr, rowvals(A), info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                    posit = Ref{Cint}(1)
                    cusparseXcsrsv2_zeroPivot(handle(), info[1], posit)
                    if posit[] >= 0
                        error("Structural/numerical zero in A at ($(posit[]),$(posit[])))")
                    end
                    $sname(handle(), ctransa, m,
                           nnz(A), alpha, desc, nonzeros(A), A.colPtr,
                           rowvals(A), info[1], X, X,
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                end
            cusparseDestroyCsrsv2Info(info[1])
            X
        end
    end
end


text:
Diagonal(I,n)## Acrobot
#TODO test
import TrajectoryOptimization
traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
urdf_folder = joinpath(traj_folder, "dynamics/urdf")
urdf_doublependulum = joinpath(urdf_folder, "doublependulum.urdf")

model = Model(urdf_doublependulum,[0.;1.]) # underactuated, only control for second "elbow" joint
n,m = model.n, model.m

# initial and goal states
x0 = [0.;0.;0.;0.]
xf = [pi;0.;0.;0.]

# costs
Q = 1e-4*Diagonal(I,n)
Qf = 250.0*Diagonal(I,n)
R = 1e-4*Diagonal(I,m)

# simulation
tf = 5.0
dt = 0.01

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)

acrobot = [model, obj_uncon]


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

"""
    AveragingCallback(semi::SemidiscretizationHyperbolic, tspan; output_directory="out",
                      filename="averaging.h5")

!!! warning "Experimental code"
    This callback is experimental and may change in any future release.

A callback that averages the flow field described by `semi` which must be a semidiscretization of
the compressible Euler equations in two dimensions. The callback records the mean velocity,
mean speed of sound, mean density, and mean vorticity for each node over the time interval given by
`tspan` and stores the results in an HDF5 file `filename` in the directory `output_directory`. Note
that this callback does not support adaptive mesh refinement ([`AMRCallback`](@ref)).
"""
struct AveragingCallback{TSpan, MeanValues, Cache}
  tspan::TSpan
  mean_values::MeanValues
  cache::Cache
  output_directory::String
  filename::String
end

function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:AveragingCallback})
  @nospecialize cb # reduce precompilation time
  averaging_callback = cb.affect!
  @unpack tspan = averaging_callback

  print(io, "AveragingCallback(tspan=", tspan, ")")
end

function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:AveragingCallback})
  @nospecialize cb # reduce precompilation time

  if get(io, :compact, false)
    show(io, cb)
  else
    averaging_callback = cb.affect!

    setup = [
             "Start time" => first(averaging_callback.tspan),
             "Final time" => last(averaging_callback.tspan)
            ]
    summary_box(io, "AveragingCallback", setup)
  end
end

function AveragingCallback(semi::SemidiscretizationHyperbolic{<:Any, <:CompressibleEulerEquations2D},
                           tspan; output_directory="out", filename="averaging.h5")
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
  mean_values = initialize_mean_values(mesh, equations, solver, cache)
  cache = create_cache(AveragingCallback, mesh, equations, solver, cache)

  averaging_callback = AveragingCallback(tspan, mean_values, cache, output_directory, filename)
  condition = (u, t, integrator) -> first(tspan) <= t <= last(tspan)

  return DiscreteCallback(condition, averaging_callback, save_positions=(false,false),
                          initialize=initialize!)
end

function initialize!(cb::DiscreteCallback{Condition,Affect!}, u_ode, t, integrator) where {Condition, Affect!<:AveragingCallback}
  averaging_callback = cb.affect!
  semi = integrator.p
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
  u = wrap_array(u_ode, mesh, equations, solver, cache)

  @trixi_timeit timer() "averaging" initialize_cache!(averaging_callback.cache, u,
                                                      mesh, equations, solver, cache)

  # avoid re-evaluating possible FSAL stages
  u_modified!(integrator, false)
  return nothing
end

# This function is called during time integration and updates the mean values according to the
# trapezoidal rule
function (averaging_callback::AveragingCallback)(integrator)
  @unpack mean_values = averaging_callback

  u_ode = integrator.u
  u_prev_ode = integrator.uprev
  semi = integrator.p
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
  u = wrap_array(u_ode, mesh, equations, solver, cache)
  u_prev = wrap_array(u_prev_ode, mesh, equations, solver, cache)

  dt = integrator.t - integrator.tprev
  tspan = averaging_callback.tspan

  integration_constant = 0.5 * dt / (tspan[2] - tspan[1]) # .5 due to trapezoidal rule

  @trixi_timeit timer() "averaging" calc_mean_values!(mean_values, averaging_callback.cache,
                                                      u, u_prev, integration_constant,
                                                      mesh, equations, solver, cache)

  # Store mean values in a file if this is the last time step
  if isfinished(integrator)
    save_averaging_file(averaging_callback, semi)
  end

  # avoid re-evaluating possible FSAL stages
  u_modified!(integrator, false)

  return nothing
end

function save_averaging_file(averaging_callback, semi::AbstractSemidiscretization)
  # Create output directory if it doesn't exist
  mkpath(averaging_callback.output_directory)

  save_averaging_file(averaging_callback, mesh_equations_solver_cache(semi)...)
end

function load_averaging_file(averaging_file, semi::AbstractSemidiscretization)
  load_averaging_file(averaging_file, mesh_equations_solver_cache(semi)...)
end

include("averaging_dg.jl")
include("averaging_dg2d.jl")

end # @muladd


text:
include("NfRel/NfRel.jl")
include("NfRel/NfRelNS.jl")
include("NfRel/NfRelOrd.jl")
include("NfRel/OrdElem.jl")
include("NfRel/Ideal.jl")
include("NfRel/FracIdeal.jl")
include("../Map/NfRelOrd.jl")
include("NfRel/ResidueRing.jl")
include("NfRel/NEQ.jl")


text: include("cpp.jl")


text:
################################################# FILE DESCRIPTION #########################################################

# This file an implementation of edge iteration.

################################################# IMPORT/EXPORT ############################################################

export EdgeIter

"""
An abstraction for alloc-free, fast edge iteration
"""
type EdgeIter <: EdgeList
   m::Int
   us::Vector{Int}
   vs::Vector{Int}
end

################################################# CONSTRUCTORS #############################################################

""" An iterator containing all edges in an adjacency matrix """
function EdgeIter(x::SparseMatrixCSC{Int})
   vs, us = findn(x)
   EdgeIter(nnz(x), us, vs)
end

""" An iterator containing all edges starting from a vertex """
function EdgeIter(x::SparseMatrixCSC{Int}, v::Int)
   adj = fadj(x, v)
   m = length(adj)
   EdgeIter(m, fill(v, m), adj)
end

""" An iterator containing all edges ending at a vertex """
function EdgeIter(x::SparseMatrixCSC{Int}, ::Colon, v::VertexID)
   EdgeIter(x.', v).'
end

""" An iterator containing all edges starting at a list of vertices """
function EdgeIter(x::SparseMatrixCSC{Int}, vlist::VertexList)
   Nerows = sum([outdegree(x, v) for v in vlist])
   us = Vector{Int}(Nerows)
   vs = Vector{Int}(Nerows)

   Cerows = 1
   for v in vlist
      p1 = x.colptr[v]
      p2 = x.colptr[v+1]
      sz = p2 - p1
      us[Cerows : (Cerows + sz - 1)] = v
      copy!(vs, Cerows, x.rowval, p1, sz)
      Cerows += sz
   end

   EdgeIter(Nerows, us, vs)
end

""" An iterator containing all edges ending at a list of vertices """
function EdgeIter(x::SparseMatrixCSC{Int}, ::Colon, vlist::VertexList)
   EdgeIter(x.', vlist).'
end

""" Split an input edge list into an iterator """
function EdgeIter(es::EdgeList)
   m = length(es)
   us = sizehint!(Vector{Int}(), m)
   vs = sizehint!(Vector{Int}(), m)

   for e in es
      push!(us, e.first)
      push!(vs, e.second)
   end

   EdgeIter(m, us, vs)
end

################################################# BASICS ###################################################################

(==)(x::EdgeIter, y::EdgeIter) = x.m == y.m && x.us == y.us && x.vs == y.vs

Base.length(x::EdgeIter) = x.m
Base.size(x::EdgeIter) = (x.m,)

Base.copy(x::EdgeIter) = EdgeIter(x.m, copy(x.us), copy(x.vs))
Base.deepcopy(x::EdgeIter) = EdgeIter(x.m, deepcopy(x.us), deepcopy(x.vs))

Base.issorted(x::EdgeIter) = true
Base.eltype(x::EdgeIter) = EdgeID

Base.transpose(x::EdgeIter) = EdgeIter(x.m, copy(x.vs), copy(x.us))

################################################# ITERATION ################################################################

Base.start(x::EdgeIter) = 1
Base.endof(x::EdgeIter) = length(x)
Base.eachindex(x::EdgeIter) = start(x) : endof(x)

Base.done(x::EdgeIter, i) = i > x.m

Base.next(x::EdgeIter, i) = (getindex(x, i), i+1)

function Base.collect(x::EdgeIter)
   n  = length(x)
   es = sizehint!(Vector{EdgeID}(), n)

   for i in 1 : n
      @inbounds u = x.us[i]
      @inbounds v = x.vs[i]
      push!(es, EdgeID(u, v))
   end

   return es
end

################################################# GETINDEX ##################################################################

""" Get the ith edge in the iterator """
Base.getindex(x::EdgeIter, i::Int) = EdgeID(x.us[i], x.vs[i])

""" Get a new iterator containing a subset of the edges """
function Base.getindex(x::EdgeIter, indxs::AbstractVector{Int})
   EdgeIter(length(indxs), x.us[indxs], x.vs[indxs])
end

""" Get a copy of the iterator """
Base.getindex(x::EdgeIter, ::Colon) = copy(x)

################################################# SHOW ######################################################################

Base.showarray(io::IO, x::EdgeIter) = show(io, x)

function Base.show(io::IO, x::EdgeIter)
   write(io, "Edge Iterator with $(x.m) values")
end

################################################# CONCATENATION #############################################################

""" Concatenate two iterators """
function Base.vcat(eit1::EdgeIter, eit2::EdgeIter)
   EdgeIter(eit1.m + eit2.m, vcat(eit1.us, eit2.us), vcat(eit1.vs, eit2.vs))
end


text:
using Flux
using Flux: flatten
using ExplainableAI: flatten_model, has_output_softmax, check_output_softmax, activation
using ExplainableAI: stabilize_denom, batch_dim_view, drop_batch_index
using Random

pseudorand(dims...) = rand(MersenneTwister(123), Float32, dims...)

# Test `activation`
@test activation(Dense(5, 2, gelu)) == gelu
@test activation(Conv((5, 5), 3 => 2, softplus)) == softplus
@test activation(BatchNorm(5, selu)) == selu
@test isnothing(activation(flatten))

# flatten_model
@test flatten_model(Chain(Chain(Chain(abs)), sqrt, Chain(relu))) == Chain(abs, sqrt, relu)
@test flatten_model(Chain(abs, sqrt, relu)) == Chain(abs, sqrt, relu)
@test_deprecated flatten_chain(Chain(identity))

# has_output_softmax
@test has_output_softmax(Chain(abs, sqrt, relu, softmax)) == true
@test has_output_softmax(Chain(abs, sqrt, relu, tanh)) == false
@test has_output_softmax(Chain(Chain(abs), sqrt, Chain(Chain(softmax)))) == true
@test has_output_softmax(Chain(Chain(abs), Chain(Chain(softmax)), sqrt)) == false
@test has_output_softmax(Chain(Dense(5, 5, softmax), Dense(5, 5, softmax))) == true
@test has_output_softmax(Chain(Dense(5, 5, softmax), Dense(5, 5, relu))) == false
@test has_output_softmax(Chain(Dense(5, 5, softmax), Chain(Dense(5, 5, softmax)))) == true
@test has_output_softmax(Chain(Dense(5, 5, softmax), Chain(Dense(5, 5, relu)))) == false

# check_output_softmax
@test_throws ArgumentError check_output_softmax(Chain(abs, sqrt, relu, softmax))

# strip_softmax
d_softmax = Dense(2, 2, softmax; init=pseudorand)
d_softmax2 = Dense(2, 2, softmax; init=pseudorand)
d_relu = Dense(2, 2, relu; init=pseudorand)
d_identity = Dense(2, 2; init=pseudorand)
# flatten to remove softmax
m = strip_softmax(Chain(Chain(abs), sqrt, Chain(Chain(softmax))))
@test m == Chain(abs, sqrt)
m1 = strip_softmax(Chain(d_relu, Chain(d_softmax)))
m2 = Chain(d_relu, d_identity)
x = rand(Float32, 2, 10)
@test typeof(m1) == typeof(m2)
@test m1(x) == m2(x)
# don't do anything if there is no softmax at the end
@test strip_softmax(Chain(Chain(abs), Chain(Chain(softmax)), sqrt)) ==
    Chain(Chain(abs), Chain(Chain(softmax)), sqrt)
@test strip_softmax(Chain(d_softmax, Chain(d_relu))) == Chain(d_softmax, Chain(d_relu))

# stabilize_denom
A = [1.0 0.0 1.0e-25; -1.0 -0.0 -1.0e-25]
S = @inferred stabilize_denom(A, 1e-3)
@test S ≈ [1.001 1e-3 1e-3; -1.001 1e-3 -1e-3]
S = @inferred stabilize_denom(Float32.(A), 1e-2)
@test S ≈ [1.01 1.0f-2 1.0f-2; -1.01 1.0f-2 -1.0f-2]

# batch_dim_view
A = [1 2; 3 4]
V = @inferred batch_dim_view(A)
@test size(V) == (2, 2, 1)
@test V[:, :, 1] == A

# drop_batch_index
I1 = CartesianIndex(5, 3, 2)
I2 = @inferred drop_batch_index(I1)
@test I2 == CartesianIndex(5, 3)
I1 = CartesianIndex(5, 3, 2, 6)
I2 = @inferred drop_batch_index(I1)
@test I2 == CartesianIndex(5, 3, 2)


text:
"""
Enumerated type for IO method used
"""
@enum IOMethods IO_METHOD_READ = 0 IO_METHOD_MMAP = 1 IO_METHOD_USERPTR = 3

function set_io_method(method::Int32 = Int32(0))
    ccall((:set_io_method,:libv4lcapture), Nothing, (Int32,), method)
end

"""
    function set_io_method(method::IOMethods)
Set the IO method to one of the following enumerated types:
```
    Video4Linux.IO_METHOD_READ
    Video4Linux.IO_METHOD_MMAP
    Video4Linux.IO_METHOD_USERPTR
```
"""
function set_io_method(method::IOMethods)
    ccall((:set_io_method,:libv4lcapture), Nothing, (Int32,), Int32(method))
end

"""
    open_device(deviceSting)

Open v4l2 video device [`/dev/video0`], returns device hanle [`Int32`]
"""
function open_device(devicename::String = "/dev/video0")
## open device
    fid = ccall((:open_device,:libv4lcapture), Int32, (Cstring,), devicename)
end

"""
    init_device(devicHandle)

Initialize already opened device. Returns 0 on success.
"""
function init_device(fid::Int32)
## init_device(fd, force_format);
    fid == -1 && error("Bad file descriptor")
    ccall((:init_device, :libv4lcapture),   Int32, ( Int32, Int32), fid, 0 )
end

"""
    start_capturing(devicHandle)

Start capturing on opened and initialized device. Returns 0 on success.
"""
function start_capturing(fid::Int32)
## start_capturing(fd);
    fid == -1 && error("Bad file descriptor")
    ccall((:start_capturing, :libv4lcapture),  Int32, ( Int32,), fid )
end

"""
    mainloop(devicHandle, frameCount)

Run main capute loop [frameCount] times. Returns 0 on success.
"""
function mainloop(fid::Int32, frame_count::Int32)
## mainloop(fd, frame_count);
    fid == -1 && error("Bad file descriptor")
    ccall((:mainloop, :libv4lcapture),  Int32, ( Int32, Int32), fid, frame_count)
end

mainloop(fid::Int32, frame_count::Int64) = mainloop(fid, convert(Int32,frame_count))

"""
    stop_capturing(devicHandle)

Stop capturing on opened device. Returns 0 on success.
"""
function stop_capturing(fid::Int32)
## stop_capturing(fd);
    fid == -1 && error("Bad file descriptor")
    ccall((:stop_capturing, :libv4lcapture), Int32, ( Int32,), fid )
end

"""
    uninit_device(devicHandle)

Uninitialize already opened device. Returns 0 on success.
"""
function uninit_device(fid::Int32)
## uninit_device();
    fid == -1 && error("Bad file descriptor")
    ccall((:uninit_device, :libv4lcapture),  Int32, () )
end

"""
    close_device(devicHandle)

Close already opened device. Returns 0 on success.
"""
function close_device(fid::Int32)
## close device
    ccall((:close_device,:libv4lcapture), Int32, (Int32,), fid)
end

##
struct RawFrame
    data::Vector{UInt8}
    width::Int
    height::Int
    pixelformat::Symbol
end

# struct to save the buffer data # and pixel format
struct Buffer8
    start::Ptr{UInt8}
    length::Csize_t
    width::UInt32
    height::UInt32
    pixelformat::UInt32
end

function get_global_bufferU8()
    # buffers = cglobal((:buffers, :libv4lcapture), Ptr{Buffer8})
    buffers = cglobal((:buffers, :libv4lcapture), Ptr{Buffer8})
    buf1 = unsafe_load(buffers,1)
    buf = unsafe_load(buf1,1)
end

"""
    copy_buffer_bytes(numBytes)

Copy [numBytes] from buffer to new buffer [`Array{UInt8,1}`]
"""
function copy_buffer_bytes(numbytes::Int64)

    buf = get_global_bufferU8()
    if numbytes > buf.length
        error("Trying to read more bytes than buffer size:")
    end
    im = zeros(UInt8, numbytes)
    unsafe_copyto!(pointer(im), buf.start, numbytes)
    return im
end

"""
    copy_buffer_bytes!(outbuffer, numBytes)

Copy [numBytes] from buffer to existing buffer outbuffer[`Array{UInt8,1}`]
"""
function copy_buffer_bytes!(im::Vector{UInt8}, numbytes::Int64)

    buf = get_global_bufferU8()
    if numbytes > buf.length
        error("Trying to read more bytes than buffer size:")
    end
    if length(im) != numbytes
        resize!(im, numbytes)
        warn("Output buffer size does not match number of bytes to be read, resizing.")
    end
    unsafe_copyto!(pointer(im), buf.start, numbytes)
    return nothing
end

"""
    copy_buffer_frame()

"""
function copy_buffer_frame()
    buf = get_global_bufferU8()
    newframe = RawFrame(zeros(UInt8, buf.length), Int(buf.width), Int(buf.height), Symbol(String(reinterpret(UInt8, [buf.pixelformat]))))
    unsafe_copyto!(pointer(newframe.data), buf.start, buf.length)
    return newframe
end

function get_v4l2_format(fid::Int32, pix::v4l2_pix_format)
    # get_v4l2_format(int fd, struct v4l2_pix_format *pix)
    fid == -1 && error("Bad file descriptor")
    ccall((:get_v4l2_format,:libv4lcapture), Int32, (Int32, Ref{v4l2_pix_format}), fid, pix)
end

function set_v4l2_format(fid::Int32, pix::v4l2_pix_format)
    # get_v4l2_format(int fd, struct v4l2_pix_format *pix)
    fid == -1 && error("Bad file descriptor")
    ccall((:set_v4l2_format,:libv4lcapture), Int32, (Int32, Ref{v4l2_pix_format}), fid, pix)
end


text:
include("blossomv.jl")
include("b-matching.jl")


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(pairwise),BrayCurtis,Adjoint{Float64, Matrix{Float64}},Adjoint{Float64, Matrix{Float64}}})   # time: 0.241804
end


text:
############################################################################################
# For the equivalent parts in OrdinaryDiffEq.jl, see:
# https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/alg_utils.jl
############################################################################################

OrdinaryDiffEq.alg_autodiff(alg::AbstractEK) = true
OrdinaryDiffEq.get_chunksize(alg::AbstractEK) = Val(0)
OrdinaryDiffEq.isfsal(alg::AbstractEK) = false

############################################
# Step size control
OrdinaryDiffEq.isadaptive(alg::AbstractEK) = true
OrdinaryDiffEq.alg_order(alg::AbstractEK) = alg.order+1
# OrdinaryDiffEq.alg_adaptive_order(alg::AbstractEK) =

# PI control is the default!
OrdinaryDiffEq.isstandard(alg::AbstractEK) = false # proportional
OrdinaryDiffEq.ispredictive(alg::AbstractEK) = false # not sure, maybe Gustafsson acceleration?

# OrdinaryDiffEq.qmin_default(alg::AbstractEK) =
# OrdinaryDiffEq.qmax_default(alg::AbstractEK) =
OrdinaryDiffEq.beta2_default(alg::AbstractEK) = 2//(5(alg.order+1))
OrdinaryDiffEq.beta1_default(alg::AbstractEK, beta2) = 7//(10(alg.order+1))
# OrdinaryDiffEq.gamma_default(alg::AbstractEK) =

# OrdinaryDiffEq.uses_uprev(alg::, adaptive::Bool) = adaptive


text:
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__)
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
CI && Pkg.instantiate()

using DynamicalBilliards
using Documenter, PyPlot, DocumenterTools

# %%
using DocumenterTools: Themes
# download the themes
for file in ("juliadynamics-lightdefs.scss", "juliadynamics-darkdefs.scss", "juliadynamics-style.scss")
    download("https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/$file", joinpath(@__DIR__, file))
end
# create the themes
for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "juliadynamics-style.scss"), String)
    theme = read(joinpath(@__DIR__, "juliadynamics-$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "juliadynamics-$(w).scss"), header*"\n"*theme)
end
# compile the themes
Themes.compile(joinpath(@__DIR__, "juliadynamics-light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "juliadynamics-dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

ioff()

makedocs(
modules=[DynamicalBilliards],
doctest=false,
depth = 1,
sitename= "DynamicalBilliards.jl",
root = @__DIR__,
format = Documenter.HTML(
    prettyurls = CI,
    assets = [
        "assets/logo.ico",
        asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css),
        ],
    ),
pages = [
    "Introduction" => "index.md",
    "High Level API" => "basic/high_level.md",
    "Visualizing & Animating" => "visualizing.md",
    "Phase Spaces" => "basic/phasespaces.md",
    "Ray-Splitting" => "ray-splitting.md",
    "Lyapunov Exponents" => "lyapunovs.md",
    "MushroomTools" => "mushroomtools.md",
    "Physics" => "physics.md",
    "Internals" => "basic/low_level.md",
    "Tutorials" => [
        "Defining a Billiard" => "tutorials/billiard_table.md",
        "Defining your own Obstacles" => "tutorials/own_obstacle.md",
        "Examples" => "tutorials/examples.md",
    ]
],
)

close("all")
ion()

if CI
    deploydocs(
        repo = "github.com/JuliaDynamics/DynamicalBilliards.jl.git",
        target = "build",
        push_preview = true
    )
end


text:
# Sequential assignment

# no longer used +++++

# Inputs are
# - `nc`: number of colleges
# - `totalCapacity` of all colleges; in multiples of  total `typeMass`
function make_test_entry_sequential(J, nc, totalCapacity)
    # objId = ObjectId(:entryOneStep);
    typeMass_jV = ones(J);
    typeMass = sum(typeMass_jV);
    capacityV = collect(range(1.0, 1.5, length = nc));
    capacityV = capacityV ./ sum(capacityV) .* totalCapacity .* typeMass;
    switches = EntrySequentialSwitches{Float64}(nTypes = J, nColleges = nc,
        typeMass_jV = typeMass_jV, capacityV = capacityV);
    return switches
    # return init_entry_decision(objId, switches)
end

function init_entry_decision(objId :: ObjectId,
    switches :: EntrySequentialSwitches{F1},
    st :: SymbolTable) where F1

    pEntryPref = init_entry_prefscale(switches, st);
    pvec = ParamVector(objId, [pEntryPref]);
    return EntrySequential(objId, pvec, ModelParams.value(pEntryPref), switches)
end

## ------------  Entry decisions

# Just a wrapper for consistent interface.
function entry_decisions(entryS :: EntrySequential{F1},
    admissionS :: AbstractAdmissionsRule{I1, F1},
    admProbFct :: AF1,
    vWork_jV :: Vector{F1}, vCollege_jcM :: Matrix{F1},
    endowPctV :: Vector{F1},
    rank_jV :: Vector{I2};
    prefShocks :: Bool = true) where
    {AF1 <: AbstractAdmProbFct{<: Real}, I1, I2 <: Integer, F1}

    return entry_sequential(entryS,  admissionS,  admProbFct,
        vWork_jV, vCollege_jcM, endowPctV,  rank_jV;
        prefShocks = prefShocks);
end

# A method that can be called directly. Avoids potential method ambiguities or dispatch errors.
function entry_sequential(entryS :: EntrySequential{F1},
    admissionS :: AbstractAdmissionsRule{I1, F1},
    vWork_jV :: Vector{F1}, vCollege_jcM :: Matrix{F1},
    endowPctV :: Vector{F1},
    rank_jV :: Vector{I2};
    prefShocks :: Bool = true) where {I1, I2 <: Integer, F1}

    nc = n_colleges(admissionS);
    nTypes = length(vWork_jV);
    enrollV = zeros(F1, nc);
    fullV = fill(false, nc);

    entryProb_jcM = zeros(F1, nTypes, nc);
    eVal_jV = zeros(F1, nTypes);

    # Loop over students in order of ranking
    for j in rank_jV
        entryProb_jcM[j, :], eVal_jV[j] = entry_decisions_one_student(
            entryS, admissionS, vWork_jV[j], vCollege_jcM[j, :],
            endowPctV[j], fullV);

        # Record enrollment
        enrollV .+= college_enrollment(entryS, entryProb_jcM[j,:], j);
        fullV = (enrollV .>= capacities(entryS));
    end

    er = EntryResults(entryS.switches, entryProb_jcM, eVal_jV, enrollV);
    @assert validate_er(er; validateFracLocal = prefShocks);
    return er
    # return entryProb_jcM, eVal_jV
end

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


text:
const CONTEXT_CACHE = Dict{Module, Context}()
const FLAGS = CXTranslationUnit_DetailedPreprocessingRecord | CXTranslationUnit_SkipFunctionBodies | CXTranslationUnit_PrecompiledPreamble #| CXTranslationUnit_SingleFileParse

Base.String(ctx::Context) = String(typeof(ctx))
header(ctx::Context) = header(typeof(ctx))
language(ctx::Context) = language(typeof(ctx))
language(::Type{Context{lang}}) where {lang} = lang

getcontext(mod::Module) = get(CONTEXT_CACHE, mod, nothing)

getjl(ctx::Context, args...; kwargs...) = getjl(typeof(ctx), args...; kwargs...)
getjl(::Type{C}, str::String; kwargs...) where {C<:Context} = getjl(C, Symbol(str); kwargs...)
getjl(::Type{C}, sym::Symbol) where {C<:Context} = esc(Symbol(replace(String(sym), r"^(struct|union|enum)\s+" => "")))

getname(ctx::Context, args...; kwargs...) = getname(typeof(ctx), args...; kwargs...)
getname(::Type{C}, str::String; kwargs...) where {C<:Context} = getname(C, Symbol(str); kwargs...)
function getname(::Type{C}, sym::Symbol; isbitstype::Bool = false) where {C<:Context}
	sym = Symbol("$(language(C))\"$(sym)\"")
	return isbitstype ? Symbol("($(sym))") : sym
end

gettype(ctx::Context, args...; kwargs...) = gettype(typeof(ctx), args...; kwargs...)
gettype(::Type{C}, str::String; kwargs...) where {C<:Context} = gettype(C, Symbol(str); kwargs...)
gettype(::Type{C}, sym::Symbol; kwargs...) where {C<:Context} = esc(getname(C, sym; kwargs...))

function getlib(ctx::Context, sym)
	for (handle, lib) in ctx.libs
		dlsym(handle, sym, throw_error = false) == C_NULL || return lib isa Symbol ? QuoteNode(lib) : lib
	end
	@warn "Failed to find `$(sym)` in any library: $(join(map(String, filter(x -> x isa Symbol, map(last, collect(ctx.libs)))), ", "))"
	return Nothing
end

function getexprs(ctx::Context, syms, blocks...)
	expr = quote end
	for block in blocks
		append!(expr.args, block.args)
	end

	blk = getblock(ctx)
	if blk.flags.jlsyms
		for (sym, jlsym, docs) in syms
			isnothing(jlsym) || push!(expr.args, :(const $(jlsym) = $(sym)))
		end
	end

	if !blk.flags.priv
		push!(expr.args, :(export $(map(((sym, jlsym, docs),) -> sym, syms)...)))
		blk.flags.jlsyms && push!(expr.args, :(export $(filter(!isnothing, map(((sym, jlsym, docs),) -> jlsym, syms))...)))
	end

	if !blk.flags.nodocs
		for (sym, jlsym, docs) in syms
			isnothing(docs) && continue
			push!(expr.args, :(@doc $(docs) $(sym.args[1])))
			isnothing(jlsym) && continue
			blk.flags.jlsyms && push!(expr.args, :(@doc $(docs) $(jlsym.args[1])))
		end
	end

	push!(expr.args, first(first(syms)))
	return expr
end

function getexprs_tu(ctx::Context, cursor::CXCursor)
	exprs = []

	for child in children(cursor)
		range = getlocation(child)
		isnothing(range) && continue  # ignore built-ins and such

		first(range).file == header(ctx) || haskey(ctx.hdrs, first(range).file) || continue
		first(range).file != header(ctx) || first(range).line > ctx.line || continue

		append!(exprs, getexprs(ctx, child))
	end

	return exprs
end

function getexprs_include(ctx::Context, cursor::CXCursor)
	exprs = []

	file = clang_getIncludedFile(cursor)
	file = _string(clang_getFileName, file)

	if !isnothing(file) && !haskey(ctx.hdrs, file)
		crng = getlocation(cursor)
		basedir = get(ctx.hdrs, first(crng).file, "")
		isimplicit = !isempty(basedir) && startswith(file, basedir)
		isexplicit = first(crng).file == header(ctx)

		if isexplicit || isimplicit
			push!(ctx.hdrs, file => isexplicit && getblock(ctx).flags.implic ? "$(dirname(file))/" : basedir)
			push!(exprs, :(include_dependency($(file))))
		end
	end

	return exprs
end

getblock(ctx::Context) = isempty(ctx.blocks) ? nothing : last(ctx.blocks)
getblock(ctx::Context, ind::Integer) = ind in eachindex(ctx.blocks) ? ctx.blocks[ind] : nothing
function getblock(ctx::Context, loc::CodeLocation)
	loc.file == header(ctx) || return nothing

	ind = searchsortedfirst(ctx.blocks, loc.line, lt = (a::CodeBlock, b::Integer) -> a.lines.stop < b)
	return getblock(ctx, ind)
end

function getlines(ctx::Context, blk::CodeBlock)
	lines = []

	io = IOBuffer(ctx.src.data[1:ctx.src.size])
	for (ind, line) in enumerate(eachline(io))
		ind in blk.lines && push!(lines, line)
	end

	return lines
end

function Base.pop!(ctx::Context)
	blk = getblock(ctx)
	if !isnothing(blk)
		pop!(ctx.blocks)
		io = IOBuffer(take!(ctx.src))
		ctx.src = IOBuffer()
		for (ind, line) in enumerate(eachline(io))
			ind in blk.lines || println(ctx.src, line)
		end
	end

	return ctx
end

# WARNING:  risky hack to trick Julia dlopen into ccalling the C dlopen with NULL (meaning to dlopen current process rather than a library)
struct _NullCString <: AbstractString end
Base.cconvert(::Type{Cstring}, ::_NullCString) = Cstring(C_NULL)

function configure!(ctx::Context, args::Vector{String})
	# find packaged include
	includedir = joinpath(dirname(dirname(libclang_path())), "lib", "clang")
	for entry in readdir(includedir)
		m = match(r"^(\d+\.\d+\.\d+)$", entry)
		if !isnothing(m)
			includedir = joinpath(includedir, m.captures[1], "include")
			break
		end
	end

	args = vcat(args, [
		"-Wno-empty-translation-unit",
		"-isystem", includedir,
	])

	libs = []
	libpaths = []
	for ind in eachindex(args)
		if args[ind] == "-L"
			if ind < length(args)
				push!(libpaths, args[ind+1])
				args[ind] = ""
				args[ind+1] = ""
			end
		elseif args[ind] == "-l"
			if ind < length(args)
				push!(libs, args[ind+1])
				args[ind] = ""
				args[ind+1] = ""
		end
		elseif startswith(args[ind], "-L")
			push!(libpaths, args[ind][3:end])
			args[ind] = ""
		elseif startswith(args[ind], "-l")
			push!(libs, args[ind][3:end])
			args[ind] = ""
		end
	end

	for lib in libs
		l = find_library(lib, libpaths)
		l = isempty(l) ? find_library("lib"*lib, libpaths) : l
		isempty(l) && error("Failed to find library $(lib) in $(join(libpaths, ", "))")

		lib = dlopen(l)
		lib == C_NULL && error("Failed to dlopen library $(l)")
		push!(ctx.libs, lib => Symbol(l))
	end
	lib = dlopen(_NullCString())
	lib == C_NULL && error("Failed to dlopen Julia process")
	push!(ctx.libs, lib => Nothing)

	ctx.ind = clang_createIndex(false, true)
	ctx.ind == C_NULL && error("Failed to create index")

	ptrs = map(pointer, args)

	unsaved = [CXUnsavedFile(
		Filename = pointer(header(ctx)),
		Contents = pointer(ctx.src.data),
		Length   = ctx.src.size,
	)]

	err = clang_parseTranslationUnit2(ctx.ind, header(ctx), ptrs, length(ptrs), unsaved, length(unsaved), FLAGS, ctx.tu)
	err == CXError_Success || error("Failed to parse translation unit $(err)")
end

function parse!(ctx::Context, mod::Module, loc::LineNumberNode, flags::NamedTuple, str::String)
	str = endswith(str, '\n') ? str : str*'\n'  # ensure new line here so countlines is consistent
	print(ctx.src, str)

	empty!(ctx.hdrs)

	blk = getblock(ctx)
	start = isnothing(blk) ? 1 : blk.lines.stop+1
	stop  = start + countlines(IOBuffer(str))-1
	push!(ctx.blocks, CodeBlock(
		mod,
		loc,
		start:stop,
		flags,
	))

	unsaved = [CXUnsavedFile(
		Filename = pointer(header(ctx)),
		Contents = pointer(ctx.src.data),
		Length   = ctx.src.size,
	)]

	err = clang_reparseTranslationUnit(ctx.tu[], length(unsaved), unsaved, FLAGS)
	err == CXError_Success || error("Failed to parse translation unit $(err)")
end

function check(ctx::Context)
	isfatal = false
	for ind in 1:clang_getNumDiagnostics(ctx.tu[])
		diag = nothing
		try
			diag = clang_getDiagnostic(ctx.tu[], ind-1)

			loc = clang_getDiagnosticLocation(diag)
			loc = getlocation(loc)
			isnothing(loc) && continue
			(loc.file == header(ctx) && loc.line > ctx.line) || continue

			blk = getblock(ctx, loc)
			isquiet = !isnothing(blk) && blk.flags.quiet
			if !isquiet
				msg = string(diag)
				msg = replace(msg, (header(ctx)*r"(?:\:\d+)*") => (!isnothing(blk) ? "$(blk.loc.file):$(blk.loc.line)" : @__FILE__))

				if !isnothing(blk)
					lines = getlines(ctx, blk)
					line = (loc.line - blk.lines.start)+1

					pre = length(lines[line]) - length(lstrip(lines[line]))
					here = lines[line][1:pre]*repeat(" ", loc.column-pre-1)*"^~~~~~ here"
					insert!(lines, line+1, here)

					msg *= "\n"*join(lines[max(1, line-2):min(length(lines), line+3)], '\n')
				end
			end

			sev = clang_getDiagnosticSeverity(diag)
			if sev in (CXDiagnostic_Error, CXDiagnostic_Fatal)
				isfatal = true
				isquiet || @error msg  (_module = !isnothing(blk) ? blk.mod : @__MODULE__)  (_file = !isnothing(blk) ? String(blk.loc.file) : @__FILE__)  (_line = !isnothing(blk) ? blk.loc.line : @__LINE__)
			elseif sev == CXDiagnostic_Warning
				isquiet || @warn  msg  (_module = !isnothing(blk) ? blk.mod : @__MODULE__)  (_file = !isnothing(blk) ? String(blk.loc.file) : @__FILE__)  (_line = !isnothing(blk) ? blk.loc.line : @__LINE__)
			else
				isquiet || @debug msg  (_module = !isnothing(blk) ? blk.mod : @__MODULE__)  (_file = !isnothing(blk) ? String(blk.loc.file) : @__FILE__)  (_line = !isnothing(blk) ? blk.loc.line : @__LINE__)
			end
		finally
			clang_disposeDiagnostic(diag)
		end
	end

	if isfatal
		pop!(ctx)
		error("Failed to create bindings, errors parsing $(String(ctx)) code")
	end
end

function advance!(ctx::Context)
	stop = ctx.line
	root = clang_getTranslationUnitCursor(ctx.tu[])
	for cursor in children(root)
		range = getlocation(cursor)
		isnothing(range) && continue

		first(range).file == header(ctx) || continue
		if first(range).line > ctx.line && last(range).line > stop
			stop = last(range).line
		end
	end
	ctx.line = stop
end

function clang_cmd(mod::Module, loc::LineNumberNode, lang::Symbol, str::String)
	cmd = Meta.parse("`$(escape_string(str))`")

	return quote
		CONTEXT_CACHE[@__MODULE__] = Context{$(QuoteNode(lang))}($(esc(cmd)).exec)
		nothing
	end
end

function clang_str(mod::Module, loc::LineNumberNode, lang::Symbol, str::String, opts::String = "")
	# TODO: allow interpolation in str?

	flags = (;
		implic = 'i' in opts,
		jlsyms = 'j' in opts,
		priv   = 'p' in opts,
		quiet  = 'q' in opts,
		ref    = 'r' in opts,
		skip   = 's' in opts,
		nodocs = 'u' in opts,
	)

	ref = getref(Context{lang}, str)
	isnothing(ref) || return ref
	flags.ref && error("Failed to parse $(String(Context{lang})) expression $(language(Context{lang}))$(repr(str))")

	haskey(CONTEXT_CACHE, mod) || error("Compiler context ($(language(Context{lang}))`...`) not created before declaring $(String(Context{lang})) bindings ($(language(Context{lang}))\"...\")")
	ctx = CONTEXT_CACHE[mod]
	language(ctx) === lang || error("Compiler context ($(language(Context{lang}))`...`) not created before declaring $(String(Context{lang})) bindings ($(language(Context{lang}))\"...\")")

	parse!(ctx, mod, loc, flags, str)
	check(ctx)
	root = clang_getTranslationUnitCursor(ctx.tu[])
	exprs = getexprs(ctx, root)
	advance!(ctx)

	return quote
		$(exprs...)
	end
end


text:
@testset "Knapsack" begin
    @testset "Knapsack: DpZeroOnePack!" begin
        dp = zeros(Int, 30)
        @test zero_one_pack!(20, [1, 3, 11], [2, 5, 30], dp) == 37
        dp = zeros(Int, 30)
        @test zero_one_pack!(10, [1, 3, 11], [20, 5, 80], dp) == 25
    end

    @testset "Knapsack: DpCompletePack!" begin
        dp = zeros(Int, 30)
        @test complete_pack!(20, [1, 2, 9], [1, 3, 20], dp) == 43
        dp = zeros(Int, 30)
        @test complete_pack!(10, [1, 3, 11], [20, 5, 80], dp) == 200
    end

    @testset "Knapsack: GaZeroOnePack!" begin
        dp = zeros(Int, 30)
        @test zero_one_pack!(20, [1, 3, 11], [2, 5, 30]) == 37
        dp = zeros(Int, 30)
        @test zero_one_pack!(10, [1, 3, 11], [20, 5, 80]) == 25
    end

    @testset "Knapsack: GaCompletePack!" begin
        dp = zeros(Int, 30)
        @test complete_pack!(20, [1, 2, 9], [1, 3, 20]) == 43
        dp = zeros(Int, 30)
        @test complete_pack!(10, [1, 3, 11], [20, 5, 80]) == 200
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
"""
    FoldsView(data, train_indices, val_indices, [obsdim])

Description
============

Create a vector-like representation of `data` where each
individual element is partition of `data` in the form of a tuple
of two data subsets (a training- and a validation subset).

The purpose of `FoldsView` is to apply a precomputed sequence of
subset assignment indices to some data container in a convenient
manner. By itself, `FoldsView` is agnostic to any particular
repartitioning strategy (such as k-folds). Instead, the
assignments, `train_indices` and `val_indices`, need to be
precomputed by such a strategy and then passed to `FoldsView`
with a concrete `data` container. The resulting object can then
be queried for its individual folds using `getindex`, or simply
iterated over.

Arguments
==========

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref) and
    [`nobs`](@ref) (see Details for more information).

- **`train_indices`** : Vector of integer vectors containing the
    sequence of training assignments. This means that each
    element is a vector of indices that describe each *training*
    subset. The length of this vector must match `val_indices`.

- **`val_indices`** : Vector of integer vectors containing the
    sequence of validation assignments. This means that each
    element is a vector of indices that describe each
    *validation* subset. The length of this vector must match
    `train_indices`.

- **`obsdim`** : Optional. If it makes sense for the type of
    `data`, `obsdim` can be used to specify which dimension of
    `data` denotes the observations. It can be specified in a
    type-stable manner as a positional argument (see
    `?LearnBase.ObsDim`), or more conveniently as a smart keyword
    argument.

Details
========

For `FoldsView` to work on some data structure, the desired type
`MyType` must implement the following interface:

- `LearnBase.getobs(data::MyType, idx, [obsdim::ObsDimension])` :
    Should return the observation(s) indexed by `idx`.
    In what form is up to the user.
    Note that `idx` can be of type `Int` or `AbstractVector`.

- `LearnBase.nobs(data::MyType, [obsdim::ObsDimension])` :
    Should return the total number of observations in `data`

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)
- Tom Breloff (Github: https://github.com/tbreloff)

Examples
=========

    # Load iris data for demonstration purposes
    X, y = MLDataUtils.load_iris()

    # Compute train- and validation-partition indices using kfolds
    train_idx, val_idx = kfolds(nobs(X), 10)

    # These two vectors contain the indices vector for each partitioning
    @assert typeof(train_idx) <: Vector{Vector{Int64}}
    @assert typeof(val_idx)   <: Vector{UnitRange{Int64}}
    @assert length(train_idx) == length(val_idx) == 10

    # Using the repartitioning as an oterator
    for (train_X, val_X) in FoldsView(X, train_idx, val_idx)
        @assert size(train_X) == (4, 135)
        @assert size(val_X) == (4, 15)
    end

    # Calling kfolds with the dataset will create
    # the FoldsView for you automatically.
    # Thus this code is equivalent to above
    for (train_X, val_X) in kfolds(X, 10)
        @assert size(train_X) == (4, 135)
        @assert size(val_X) == (4, 15)
    end

    # leavout is a shortcut for setting k = nobs(X)
    for (train_X, val_X) in leaveout(X)
        @assert size(val_X) == (4, 1)
    end

see also
=========

[`kfolds`](@ref), [`leaveout`](@ref), [`splitobs`](@ref),
[`DataSubset`](@ref)
"""
struct FoldsView{T,D,O,A1<:AbstractArray,A2<:AbstractArray} <: DataView{T,D}
    data::D
    train_indices::A1
    val_indices::A2
    obsdim::O

    function FoldsView{T,D,O,A1,A2}(
            data::D, train_indices::A1, val_indices::A2, obsdim::O) where {T,D,O,A1<:AbstractArray,A2<:AbstractArray}
        n = nobs(data, obsdim)
        (eltype(train_indices) <: AbstractArray{Int}) || throw(ArgumentError("The parameter \"train_indices\" must be an array of integer arrays"))
        (eltype(val_indices)  <: AbstractArray{Int}) || throw(ArgumentError("The parameter \"val_indices\" must be an array of integer arrays"))
        2 <= length(train_indices) <= n || throw(ArgumentError("The amount of train- and validation-indices must be within 2:$n respectively"))
        length(train_indices) == length(val_indices) || throw(DimensionMismatch("The amount of train- and validation-indices must match"))
        new{T,D,O,A1,A2}(data, train_indices, val_indices, obsdim)
    end
end

function FoldsView(data::D, train_indices::A1, val_indices::A2, obsdim::O) where {D,O,A1<:AbstractArray,A2<:AbstractArray}
    n = nobs(data, obsdim)
    # TODO: Move this back into the inner constructor after the
    #       "T = typeof(...)" line below is removed
    (1 <= minimum(minimum.(train_indices)) && maximum(maximum.(train_indices)) <= n) || throw(DimensionMismatch("All training indices must be within 1:$n"))
    (1 <= minimum(minimum.(val_indices))  && maximum(maximum.(val_indices))  <= n) || throw(DimensionMismatch("All validation indices must be within 1:$n"))
    # FIXME: In 0.6 it should be possible to compute just the return
    #        type without executing the function
    T = typeof((datasubset(data, train_indices[1], obsdim), datasubset(data, val_indices[1], obsdim)))
    FoldsView{T,D,O,A1,A2}(data, train_indices, val_indices, obsdim)
end

FoldsView(data, train_indices::AbstractArray, val_indices::AbstractArray; obsdim = default_obsdim(data)) =
    FoldsView(data, train_indices, val_indices, convert(LearnBase.ObsDimension,obsdim))

function FoldsView(data::T, train_indices::AbstractArray, val_indices::AbstractArray, obsdim) where T<:DataView
    @assert obsdim == data.obsdim
    @warn string("Trying to nest a ", T.name, " into an FoldsView, which is not supported. Returning FoldsView(parent(_)) instead")
    FoldsView(parent(data), train_indices, val_indices, obsdim)
end

# compare if both FoldsViews describe the same folds of the same data
# we don't care how the indices are stored, just that they match
# in their order and values
function Base.:(==)(fv1::FoldsView,fv2::FoldsView)
    fv1.data == fv2.data &&
        all(all(i1==i2 for (i1,i2) in zip(I1,I2)) for (I1,I2) in zip(fv1.train_indices,fv2.train_indices)) &&
        all(all(i1==i2 for (i1,i2) in zip(I1,I2)) for (I1,I2) in zip(fv1.val_indices,fv2.val_indices)) &&
        fv1.obsdim == fv2.obsdim
end

nobs(iter::FoldsView) = nobs(iter.data, iter.obsdim)
Base.parent(iter::FoldsView) = iter.data
Base.length(iter::FoldsView) = length(iter.train_indices)

function Base.getindex(iter::FoldsView, i::Int)
    (datasubset(iter.data, iter.train_indices[i], iter.obsdim),
     datasubset(iter.data, iter.val_indices[i], iter.obsdim))
end

function Base.getindex(iter::FoldsView, i::AbstractVector)
    FoldsView(iter.data, iter.train_indices[i], iter.val_indices[i], iter.obsdim)
end

# compatibility with nested functions
default_obsdim(iter::FoldsView) = iter.obsdim

function Base.summary(A::FoldsView)
    string(length(A), "-fold ", typeof(A).name,
           " of ", nobs(A), " observations")
end

_datastr(data) = summary(data)
_datastr(data::Tuple) = string("(", join(map(_datastr, data), ", "), ")")
function Base.show(io::IO, ::MIME"text/plain", A::FoldsView)
    println(io, summary(A), ":")
    println(io, "  data: ", _datastr(A.data))
    # TODO: consider the case where fold sizes differ
    ntrain = length(first(A.train_indices))
    println(io, "  training: ", ntrain, " observations/fold")
    nval = length(first(A.val_indices))
    println(io, "  validation: ", nval, " observations/fold")
    print(io, "  obsdim: ", obsdim_string(A.obsdim))
end

"""
    kfolds(n::Integer, [k = 5]) -> Tuple

Compute the train/validation assignments for `k` repartitions of
`n` observations, and return them in the form of two vectors. The
first vector contains the index-vectors for the training subsets,
and the second vector the index-vectors for the validation subsets
respectively. A general rule of thumb is to use either `k = 5` or
`k = 10`. The following code snippet generates the indices
assignments for `k = 5`

```julia
julia> train_idx, val_idx = kfolds(10, 5);
```

Each observation is assigned to the validation subset once (and
only once). Thus, a union over all validation index-vectors
reproduces the full range `1:n`. Note that there is no random
assignment of observations to subsets, which means that adjacent
observations are likely to be part of the same validation subset.

```julia
julia> train_idx
5-element Array{Array{Int64,1},1}:
 [3,4,5,6,7,8,9,10]
 [1,2,5,6,7,8,9,10]
 [1,2,3,4,7,8,9,10]
 [1,2,3,4,5,6,9,10]
 [1,2,3,4,5,6,7,8]

julia> val_idx
5-element Array{UnitRange{Int64},1}:
 1:2
 3:4
 5:6
 7:8
 9:10
```
"""
function kfolds(n::Integer, k::Integer = 5)
    2 <= k <= n || throw(ArgumentError("n must be positive and k must to be within 2:$(max(2,n))"))
    # Compute the size of each fold. This is important because
    # in general the number of total observations might not be
    # divideable by k. In such cases it is custom that the remaining
    # observations are divided among the folds. Thus some folds
    # have one more observation than others.
    sizes = fill(floor(Int, n/k), k)
    for i = 1:(n % k)
        sizes[i] = sizes[i] + 1
    end
    # Compute start offset for each fold
    offsets = cumsum(sizes) .- sizes .+ 1
    # Compute the validation indices using the offsets and sizes
    val_indices = map((o,s)->(o:o+s-1), offsets, sizes)
    # The train indices are then the indicies not in validation
    train_indices = map(idx->setdiff(1:n,idx), val_indices)
    # We return a tuple of arrays
    train_indices, val_indices
end

"""
    kfolds(data, [k = 5], [obsdim]) -> FoldsView

Repartition a `data` container `k` times using a `k` folds
strategy and return the sequence of folds as a lazy
[`FoldsView`](@ref). The resulting `FoldsView` can then be
indexed into or iterated over. Either way, only data subsets are
created. That means that no actual data is copied until
[`getobs`](@ref) is invoked.

Conceptually, a k-folds repartitioning strategy divides the given
`data` into `k` roughly equal-sized parts. Each part will serve
as validation set once, while the remaining parts are used for
training. This results in `k` different partitions of `data`.

In the case that the size of the dataset is not dividable by the
specified `k`, the remaining observations will be evenly
distributed among the parts.

```julia
for (x_train, x_val) in kfolds(X, k = 10)
    # code called 10 times
    # nobs(x_val) may differ up to ±1 over iterations
end
```

Multiple variables are supported (e.g. for labeled data)

```julia
for ((x_train, y_train), val) in kfolds((X, Y), k = 10)
    # ...
end
```

By default the folds are created using static splits. Use
[`shuffleobs`](@ref) to randomly assign observations to the
folds.

```julia
for (x_train, x_val) in kfolds(shuffleobs(X), k = 10)
    # ...
end
```

see [`FoldsView`](@ref) for more info, or [`leaveout`](@ref) for
a related function.
"""
function kfolds(data, k::Integer, obsdim)
    n = nobs(data, obsdim)
    train_indices, val_indices = kfolds(n, k)
    FoldsView(data, train_indices, val_indices, obsdim)
end

kfolds(data, k::Integer; obsdim = default_obsdim(data)) =
    kfolds(data, k, convert(LearnBase.ObsDimension,obsdim))

kfolds(data; k = 5, obsdim = default_obsdim(data)) =
    kfolds(data, k, convert(LearnBase.ObsDimension,obsdim))

kfolds(data, obsdim::Union{Tuple,ObsDimension}) =
    kfolds(data, 5, obsdim)

"""
    leaveout(n::Integer, [size = 1]) -> Tuple

Compute the train/validation assignments for `k ≈ n/size`
repartitions of `n` observations, and return them in the form of
two vectors. The first vector contains the index-vectors for the
training subsets, and the second vector the index-vectors for the
validation subsets respectively. Each validation subset will have
either `size` or `size+1` observations assigned to it. The
following code snippet generates the index-vectors for `size = 2`.

```julia
julia> train_idx, val_idx = leaveout(10, 2);
```

Each observation is assigned to the validation subset once (and
only once). Thus, a union over all validation index-vectors
reproduces the full range `1:n`. Note that there is no random
assignment of observations to subsets, which means that adjacent
observations are likely to be part of the same validation subset.

```julia
julia> train_idx
5-element Array{Array{Int64,1},1}:
 [3,4,5,6,7,8,9,10]
 [1,2,5,6,7,8,9,10]
 [1,2,3,4,7,8,9,10]
 [1,2,3,4,5,6,9,10]
 [1,2,3,4,5,6,7,8]

julia> val_idx
5-element Array{UnitRange{Int64},1}:
 1:2
 3:4
 5:6
 7:8
 9:10
```
"""
function leaveout(n::Integer, size::Integer = 1)
    1 <= size <= floor(n/2) || throw(ArgumentError("size must to be within 1:$(floor(Int,n/2))"))
    k = floor(Int, n / size)
    kfolds(n, k)
end

"""
    leaveout(data, [size = 1], [obsdim]) -> FoldsView

Repartition a `data` container using a k-fold strategy, where `k`
is chosen in such a way, that each validation subset of the
resulting folds contains roughly `size` observations. Defaults to
`size = 1`, which is also known as "leave-one-out" partitioning.

The resulting sequence of folds is returned as a lazy
[`FoldsView`](@ref), which can be index into or iterated over.
Either way, only data subsets are created. That means no actual
data is copied until [`getobs`](@ref) is invoked.

```julia
for (train, val) in leaveout(X, size = 2)
    # if nobs(X) is dividable by 2,
    # then nobs(val) will be 2 for each iteraton,
    # otherwise it may be 3 for the first few iterations.
end
```

see [`FoldsView`](@ref) for more info, or [`kfolds`](@ref) for a
related function.
"""
function leaveout(data, size, obsdim)
    n = nobs(data, obsdim)
    train_indices, val_indices = leaveout(n, size)
    FoldsView(data, train_indices, val_indices, obsdim)
end

leaveout(data, size::Integer; obsdim = default_obsdim(data)) =
    leaveout(data, size, convert(LearnBase.ObsDimension,obsdim))

leaveout(data; size = 1, obsdim = default_obsdim(data)) =
    leaveout(data, size, convert(LearnBase.ObsDimension,obsdim))

leaveout(data, obsdim::Union{Tuple,ObsDimension}) =
    leaveout(data, 1, obsdim)


text:
@testset "RCall-based plot Test" begin
  # testing for absence of errors, not for correctness

  #net = readTopology("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);")
  net = readTopology("(((Ag,(#H1:7.159::0.056,((Ak,(E:0.08,#H2:0.0::0.004):0.023):0.078,(M:0.0)#H2:::0.996):2.49):2.214):0.026,(((((Az:0.002,Ag2:0.023):2.11,As:2.027):1.697)#H1:0.0::0.944,Ap):0.187,Ar):0.723):5.943,(P,20):1.863,165);");
  @test_logs plot(net,:RCall);
  @test_logs (:warn, "At least one non-missing edge length: plotting any missing length as 1.0") plot(net,:RCall, useEdgeLength=true);
  @test_logs plot(net,:RCall, showTipLabel=false);
  @test_logs plot(net,:RCall, showNodeNumber=true, showIntNodeLabel=true);
  @test_logs plot(net,:RCall, tipOffset=1, showGamma=true);
  @test_logs plot(net,:RCall, showEdgeLength=true, showEdgeNumber=true);
  @test_logs plot(net,:RCall, edgeColor="tomato4",minorHybridEdgeColor="skyblue",
          majorHybridEdgeColor="tan");
  dat = DataFrame(node=[-5,-10,-1],bs=["90","95","100"],edge=[11,22,26]);
  @test_logs (:warn, "Some node numbers in the nodeLabel data frame are not found in the network:\n -1") plot(net,:RCall, nodeLabel=dat);
  @test_logs plot(net,:RCall, edgeLabel=dat[!,[:edge,:bs]]);

  @test_logs plot(net,:RCall, style=:majortree, arrowlen=0.1);
  @test_logs (:warn, "Style bogus is unknown. Defaulted to :fulltree.") plot(net,:RCall, style=:bogus);

  # coverage for the nomajorchild
  net2 = readTopology("((((B)#H1)#H2,((D,C,#H2)S1,(#H1:::.8,A)S2)S3)S4);")
  @test_logs plot(net2, :R, style=:majortree);

  # plot based on RCall and ape:
  tre = readTopology("(((((((1,2),3),4),5),(6,7)),(8,9)),10);");
  # fixit: plot(tre, :ape)
end

max_stars_repo_name: cecileane/PhyloPlots

max_issues_repo_name: cecileane/PhyloPlots

max_forks_repo_name: cecileane/PhyloPlots


text:
@testset "GroupedArray" begin
    N = 10
    a1 = collect(1:N)
    g1 = group(a1)
    @test g1 == GroupedArray(collect(UInt32, 1:N), N)
    @test size(g1) == (N,)
    @test length(g1) == N
    @test g1[1] == UInt32(1)
    @test g1[1:2] == [UInt32(1), UInt32(2)]
    @test g1[g1.<=2] == g1[[1,2]] == g1[1:2]

    a2 = [1,2]
    @test_throws DimensionMismatch group(a1, a2)

    a = rand(N)
    g = group(a)
    @test factorize!(g) == g
end