text:
module TestExamples3DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_3d_dgsem")

@testset "MHD" begin
  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.017590099293094203, 0.017695875823827714, 0.017695875823827686, 0.017698038279620777, 0.07495006099352074, 0.010391801950005755, 0.010391801950005759, 0.010393502246627087, 2.524766553484067e-16],
      linf = [0.28173002819718196, 0.3297583616136297, 0.32975836161363004, 0.356862935505337, 1.2893514981209626, 0.10950981489747313, 0.10950981489747136, 0.11517234329681891, 2.0816911067714202e-15])
  end

  @trixi_testset "elixir_mhd_ec.jl with initial_condition=initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [4.270231310667203e-16, 2.4381208042014784e-15, 5.345107673575357e-15, 3.00313882171883e-15, 1.7772703118758417e-14, 1.0340110783830874e-15, 1.1779095371939702e-15, 9.961878521814573e-16, 8.1201730630719145e-16],
      linf = [2.4424906541753444e-15, 2.881028748902281e-14, 2.4646951146678475e-14, 2.3092638912203256e-14, 2.3447910280083306e-13, 1.7763568394002505e-14, 1.0436096431476471e-14, 2.042810365310288e-14, 7.057203733035201e-15],
      atol = 1000*eps(),
      initial_condition=initial_condition_constant)
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.0032217291057246157, 0.009511644936958913, 0.004217358459420256, 0.011591709179125335, 0.009456218722393708, 0.00916500047763897, 0.005069863732625444, 0.011503011541926135, 0.003988175543749985],
      linf = [0.01188593784273051, 0.03638015998373141, 0.01568200398945724, 0.04666974730787579, 0.031235294705421968, 0.03316343064943483, 0.011539436992528018, 0.04896687646520839, 0.018714054039927555])
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with flux_derigs_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.003755235939722358, 0.009062519246840721, 0.004096299856228109, 0.011429935838448906, 0.006897420817511043, 0.00900886245212482, 0.004926537542780259, 0.01153285554590683, 0.0037842060148666886],
      linf = [0.012982853115883541, 0.0320228076558316, 0.011575276754611022, 0.04425778643430531, 0.02478109022285846, 0.03198699034954189, 0.009761077061886558, 0.04433669321441455, 0.01618905441148782],
      volume_flux = (flux_derigs_etal, flux_nonconservative_powell))
  end

  @trixi_testset "elixir_mhd_alfven_wave_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"),
      l2   = [0.001879021634926363, 0.007032724521848316, 0.0032793932234187325, 0.009056594733320348, 0.007514150120617965, 0.007328739509868727, 0.00309794018112387, 0.009026356949274878, 0.0035732583778049776],
      linf = [0.013734346970999622, 0.06173467158736011, 0.02183946452704291, 0.06258216169457917, 0.03672304497348122, 0.055120532123884625, 0.018202716205672487, 0.06133688282205586, 0.019888161885935608],
      tspan = (0.0, 0.25))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with Orszag-Tang setup + flux_hll" begin
    # OBS! This setup does not make much sense and is only used to exercise all components of the
    # flux_hll implementation
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.004391143689111404, 0.04144737547475548, 0.041501307637678286, 0.04150353006408862, 0.03693135855995625, 0.021125605214031118, 0.03295607553556973, 0.03296235755245784, 7.16035229384135e-6],
      linf = [0.017894703320895378, 0.08486850681397005, 0.0891044523165206, 0.08492024792056754, 0.10448301878352373, 0.05381260695579509, 0.0884774018719996, 0.07784546966765199, 7.71609149516089e-5],
      initial_condition = initial_condition_orszag_tang,
      surface_flux = (flux_hll, flux_nonconservative_powell),
      volume_flux  = (flux_central, flux_nonconservative_powell),
      coordinates_min = (0, 0, 0),
      coordinates_max = (1, 1, 1),
      initial_refinement_level=3,
      cfl = 1.1,
      tspan = (0.0, 0.06))
  end
end

end # module


text:
# an ordering to be used in sorting DictArrays or Labeledarrays.
immutable AbstractArrayLT{N,O,F} <: Base.Ordering
  ords::O #NTuple{N,Base.Ordering}
  fields::F #NTuple{N,Vector}
end

AbstractArrayLT(arr::Union{DictArray,LabeledArray}, axis::Integer, field_names...; kwargs...) = begin
  kv = Dict(kwargs)
  ords = map(field_names) do field
    lt = get(kv, symbol(field, "_lt"), isless)
    by = get(kv, symbol(field, "_by"), identity)
    rev = get(kv, symbol(field, "_rev"), false)
    order = get(kv, symbol(field, "_ord"), Base.Forward)
    Base.Order.ord(lt, by, rev, order)
  end
  ndimsarr = ndims(arr)
  fields = map(field_names) do name
    #collect(selectfield(arr, name)[get(kv, symbol(name, "_coords"), ntuple(d->d==axis ? Colon() : 1, ndimsarr))...])
    selectfield(arr, name)[get(kv, symbol(name, "_coords"), ntuple(d->d==axis ? Colon() : 1, ndimsarr))...].a
  end
  AbstractArrayLT{length(field_names),typeof(ords),typeof(fields)}(ords, fields)
end

@generated Base.Sort.lt{N}(ltobj::AbstractArrayLT{N}, x::Integer, y::Integer) = begin
  comparison(i) = quote
    atix = getindexvalue(ltobj.fields[$i], x)
    atiy = getindexvalue(ltobj.fields[$i], y)
    ord = ltobj.ords[$i]
    if base_lt_helper(ord, atix, atiy)
      return true
    elseif base_lt_helper(ord, atiy, atix)
      return false
    end
  end
  exprs = map(n->comparison(n), 1:N)
  Expr(:block, exprs..., :(return false))
end

# @inline macro actually enhanced the performance.
@inline base_lt_helper{O,X,Y}(ord::O, x::X, y::Y) = Base.lt(ord, x, y)
@inline base_lt_helper{O,F<:AbstractFloat}(ord::O, x::F, y::F) =
  (isnan(x) && !isnan(y)) || (!isnan(x) && !isnan(y) && Base.lt(ord, x, y))
@inline base_lt_helper{F<:AbstractFloat}(ord::Base.Order.ReverseOrdering, x::F, y::F) =
  (!isnan(x) && isnan(y)) || (!isnan(x) && !isnan(y) && Base.lt(ord, x, y))
@inline base_lt_helper{O,X,Y}(ord::O, x::Nullable{X}, y::Nullable{Y}) =
  (x.isnull && !y.isnull) || (!x.isnull && !y.isnull && Base.lt(ord, x.value, y.value))
@inline base_lt_helper{X,Y}(ord::Base.Order.ReverseOrdering, x::Nullable{X}, y::Nullable{Y}) =
  (!x.isnull && y.isnull) || (!x.isnull && !y.isnull && Base.lt(ord, x.value, y.value))

sortpermbase(arr::Union{DictArray, LabeledArray}, axis::Integer, algorithm, order::Base.Ordering) = begin
  sorted_indices = collect(1:size(arr, axis))
  sort!(sorted_indices, algorithm, order)
  permuted_coords = ntuple(ndims(arr)) do d
    if d == axis
      sorted_indices
    else
      Colon()
    end
  end
  permuted_coords
end
Base.sortperm(arr::Union{DictArray,LabeledArray}, axis::Integer, fields...; alg=Base.Sort.defalg(arr), kwargs...) = begin
  ordering = AbstractArrayLT(arr, axis, fields...;kwargs...)
  sortpermbase(arr, axis, alg, ordering)
end
sortbase(arr::Union{DictArray,LabeledArray}, axis::Integer, algorithm, order::Base.Ordering) =
  arr[sortpermbase(arr, axis, algorithm, order)...]

"""

`sort(arr, axis, fields... [; alg=..., ...])`

Sort a `DictArray` or `LabeledArray` along some axis.

##### Arguments

* `arr` : either a `DictArray` or a `LabeledArray`.
* `axis` : an axis direction integer to denote which direction to sort along. If omitted, axis=1.
* `fields...` : the names of fields to determine the order. The preceding ones have precedence over the later ones. Note only the components [1,...,1,:,1,...1], where : is placed at the axis position, will be used out of each field. If omitted, all fields will be used in their order for `DictArray` and the axis along the `axis` direction for `LabeledArray`.
* optionally, `alg=algorithm` determines the sorting algorithm. `fieldname_lt=ltfunc` sets the less-than function for the field fieldname, and similarly for `by`/`rev`/`ord`.

##### Examples

```julia
julia> t = larr(a=[3 3 2;7 5 3], b=[:b :a :c;:d :e :f], axis1=["X","Y"])
2 x 3 LabeledArray

  |1   |2   |3
--+----+----+----
  |a b |a b |a b
--+----+----+----
X |3 b |3 a |2 c
Y |7 d |5 e |3 f

julia> sort(t, 1, :a)
2 x 3 LabeledArray

  |1   |2   |3
--+----+----+----
  |a b |a b |a b
--+----+----+----
X |3 b |3 a |2 c
Y |7 d |5 e |3 f

julia> sort(t, 2, :a)
2 x 3 LabeledArray

  |1   |2   |3
--+----+----+----
  |a b |a b |a b
--+----+----+----
X |2 c |3 b |3 a
Y |3 f |7 d |5 e

julia> sort(t, 2, :a, :b)
2 x 3 LabeledArray

  |1   |2   |3
--+----+----+----
  |a b |a b |a b
--+----+----+----
X |2 c |3 a |3 b
Y |3 f |5 e |7 d

julia> sort(t, 2, :a, :b, a_rev=true)
2 x 3 LabeledArray

  |1   |2   |3
--+----+----+----
  |a b |a b |a b
--+----+----+----
X |3 a |3 b |2 c
Y |5 e |7 d |3 f
```

"""
Base.sort(arr::Union{DictArray,LabeledArray}, axis::Integer, fields...; alg=Base.Sort.defalg(arr), kwargs...) = begin
  ordering = AbstractArrayLT(arr, axis, fields...;kwargs...)
  sortbase(arr, axis, alg, ordering)
end
Base.sort(arr::DictArray, axis::Integer; alg=Base.Sort.defalg(arr), kwargs...) = sort(arr, axis, keys(arr)...; alg=alg, kwargs...)
Base.sort(arr::LabeledArray, axis::Integer; alg=Base.Sort.defalg(arr), kwargs...) = begin
  coords = sortperm(arr, axis; alg=alg, kwargs...)
  arr[coords...]
end
Base.sort(arr::Union{DictArray,LabeledArray}; alg=Base.Sort.defalg(arr), kwargs...) = sort(arr, 1; alg=alg, kwargs...)
Base.sortperm(arr::DictArray, axis::Integer; alg=Base.Sort.defalg(arr), kwargs...) = sortperm(arr, axis, keys(arr)...; alg=alg, kwargs...)
Base.sortperm(arr::LabeledArray, axis::Integer; alg=Base.Sort.defalg(arr), kwargs...) = begin
  axis_processed = sortperm_inner_convert_to_dictarray_if_necessary(pickaxis(arr,axis))
  ordering = AbstractArrayLT(axis_processed, 1, keys(axis_processed)...;kwargs...)
  axisorder = sortpermbase(axis_processed, 1, alg, ordering)[1]
  coords = ntuple(ndims(arr)) do d
    if d == axis
      axisorder
    else
      Colon()
    end
  end
  coords
end
Base.sortperm(arr::Union{DictArray,LabeledArray}; alg=Base.Sort.defalg(arr), kwargs...) = sortperm(arr, 1; alg=alg, kwargs...)
sortperm_inner_convert_to_dictarray_if_necessary(arr::DictArray) = arr
sortperm_inner_convert_to_dictarray_if_necessary(arr::AbstractArray) = create_dictarray_nocheck(create_ldict_nocheck(:dummy=>arr))

ext: jl


text:
using LinearAlgebra
A = [0 1/2 1/4; 1 1 2; 2 1 1.5]
F = lu(A) # like "A^{-1}"
##
println("P = ")
display(F.P)
println("L = ")
display(F.L)
println("U = ")
display(F.U)
##
@show norm(F.P*A - F.L*F.U)


text:
"""
    CoordinateBoundaryConditions(left, right)

A set of two `BoundaryCondition`s to be applied along a coordinate x, y, or z.

The `left` boundary condition is applied on the negative or lower side of the coordinate
while the `right` boundary condition is applied on the positive or higher side.
"""
mutable struct CoordinateBoundaryConditions{L, R}
     left :: L
    right :: R
end

#####
##### Some aliases to make life easier.
#####

const CBC = CoordinateBoundaryConditions

PeriodicBCs() = CBC(PeriodicBC(), PeriodicBC())

# Here we overload setproperty! and getproperty to permit users to call
# the 'left' and 'right' bcs in the z-direction 'bottom' and 'top'
# and the 'left' and 'right' bcs in the y-direction 'south' and 'north'.
Base.setproperty!(cbc::CBC, side::Symbol, bc) = setbc!(cbc, Val(side), bc)

setbc!(cbc::CBC, ::Val{S}, bc) where S = setfield!(cbc, S, bc)
setbc!(cbc::CBC, ::Val{:bottom}, bc) = setfield!(cbc, :left,  bc)
setbc!(cbc::CBC, ::Val{:top},    bc) = setfield!(cbc, :right, bc)
setbc!(cbc::CBC, ::Val{:south},  bc) = setfield!(cbc, :left,  bc)
setbc!(cbc::CBC, ::Val{:north},  bc) = setfield!(cbc, :right, bc)

Base.getproperty(cbc::CBC, side::Symbol) = getbc(cbc, Val(side))

getbc(cbc::CBC, ::Val{S}) where S = getfield(cbc, S)
getbc(cbc::CBC, ::Val{:bottom}) = getfield(cbc, :left)
getbc(cbc::CBC, ::Val{:top})    = getfield(cbc, :right)
getbc(cbc::CBC, ::Val{:south})  = getfield(cbc, :left)
getbc(cbc::CBC, ::Val{:north})  = getfield(cbc, :right)


text:
using BenchmarkTools
using Test

using Groups
using Groups.New

function wl_ball(F; radius::Integer)
    g, state = iterate(F)
    while length(word(g)) <= radius
        res = iterate(F, state)
        isnothing(res) && break
        g, state = res
    end
    elts = collect(state.seen)
    elts = resize!(elts, length(elts)-1)
    return elts
end

@testset "Benchmarks" begin
    N = 4

    @testset "iteration: FreeGroup" begin
        FN = FreeGroup(N)
        R = 8

        let G = FN
            S = unique([gens(G); inv.(gens(G))])

            sizes1 = last(Groups.wlmetric_ball(S, radius=R, threading=false))
            sizes2 = last(Groups.wlmetric_ball(S, radius=R, threading=true))

            l = length(wl_ball(G, radius=R))

            @test sizes1 == sizes2
            @test last(sizes1) == l

            @info "Ball of radius $R in $(parent(first(S)))" sizes=sizes1
            @info "serial"
            @time Groups.wlmetric_ball(S, radius=R, threading=false)
            @info "threaded"
            @time Groups.wlmetric_ball(S, radius=R, threading=true)
            @info "iteration"
            @time wl_ball(G, radius=R)
        end
    end

    @testset "iteration: SAut(F_n)" begin
        R = 4
        FN = FreeGroup(N)
        SAutFN = New.SpecialAutomorphismGroup(FN)

        let G = SAutFN
            S = unique([gens(G); inv.(gens(G))])

            sizes1 = last(Groups.wlmetric_ball(S, radius=R, threading=false))
            sizes2 = last(Groups.wlmetric_ball(S, radius=R, threading=true))

            l = length(wl_ball(G, radius=R))

            @test sizes1 == sizes2
            @test last(sizes1) == l

            @info "Ball of radius $R in $(parent(first(S)))" sizes=sizes1
            @info "serial"
            @time Groups.wlmetric_ball(S, radius=R, threading=false)
            @info "threaded"
            @time Groups.wlmetric_ball(S, radius=R, threading=true)
            @info "iteration"
            @time wl_ball(G, radius=R)
        end
    end
end


text:
module TestAqua

import Aqua
import Baselet
using Test

Aqua.test_all(Baselet)

@testset "Compare Project.toml and test/Project.toml" begin
    Aqua.test_project_extras(Baselet)
end

end  # module


text:
#==============================================================================#
# Recipe for plotting field grids
#==============================================================================#

# TODO: Add plotting recipes for the system and fields
# Plot in terms of scaled cell coordinates, not grid index


text:
using Documenter, RainFARM

push!(LOAD_PATH,"../src/")
makedocs(sitename="RainFARM documentation")

deploydocs(
    repo = "github.com/jhardenberg/RainFARM.jl.git",
)


text:
# TODO: concreate typing

abstract type WeaveChunk end
abstract type Inline end
abstract type WeaveFormat end

mutable struct WeaveDoc
    source::AbstractString
    basename::AbstractString
    path::AbstractString
    chunks::Vector{WeaveChunk}
    cwd::AbstractString
    format::Any
    doctype::String
    header_script::String
    header::Dict
    chunk_defaults::Dict{Symbol,Any}
end

struct ChunkOutput
    code::String
    stdout::String
    rich_output::String
    figures::Vector{String}
end

mutable struct CodeChunk <: WeaveChunk
    content::String
    number::Int
    start_line::Int
    optionstring::String
    options::Dict{Symbol,Any}
    output::AbstractString
    rich_output::AbstractString
    figures::Vector{String}
    result::Vector{ChunkOutput}
end

function CodeChunk(content, number, start_line, optionstring, options)
    return CodeChunk(
        string(rstrip(content), '\n'), # normalize end of chunk)
        number,
        start_line,
        optionstring,
        options,
        "",
        "",
        AbstractString[],
        ChunkOutput[]
    )
end

mutable struct DocChunk <: WeaveChunk
    content::Vector{Inline}
    number::Int
    start_line::Int
end

struct InlineText <: Inline
    content::String
    number::Int
end

mutable struct InlineCode <: Inline
    content::String
    number::Int
    ctype::Symbol
    output::String
    rich_output::String
    figures::Vector{String}
end
InlineCode(content, number, ctype) = InlineCode(content, number, ctype, "", "", String[])


text:
module AcmeFlat

"""
    mutable struct Product

A product from Acme's catalog

Fields:
   `productId`: The unique identifier for a product
   `productName`: Name of the product
   `price`: The price of the product
   `tags`: Tags for the product
"""
mutable struct Product
    productId::Int
    productName::String
    price::Float64
    tags::Array{String}
end

end # end module


text:
# ---
# title: 765. Couples Holding Hands
# id: problem765
# author: Indigo
# date: 2021-02-14
# difficulty: Hard
# categories: Greedy, Union Find, Graph
# link: <https://leetcode.com/problems/couples-holding-hands/description/>
# hidden: true
# ---
#
# N couples sit in 2N seats arranged in a row and want to hold hands. We want to
# know the minimum number of swaps so that every couple is sitting side by side.
# A _swap_ consists of choosing **any** two people, then they stand up and
# switch seats.
#
# The people and seats are represented by an integer from `0` to `2N-1`, the
# couples are numbered in order, the first couple being `(0, 1)`, the second
# couple being `(2, 3)`, and so on with the last couple being `(2N-2, 2N-1)`.
#
# The couples' initial seating is given by `row[i]` being the value of the
# person who is initially sitting in the i-th seat.
#
# **Example 1:**
#
#
#
#     Input: row = [0, 2, 1, 3]
#     Output: 1
#     Explanation: We only need to swap the second (row[1]) and third (row[2]) person.
#
#
# **Example 2:**
#
#
#
#     Input: row = [3, 2, 0, 1]
#     Output: 0
#     Explanation: All couples are already seated side by side.
#
#
# **Note:**
#
#   1. `len(row)` is even and in the range of `[4, 60]`.
#   2. `row` is guaranteed to be a permutation of `0...len(row)-1`.
#
#
## @lc code=start
using LeetCode

function min_swaps_couples(row::Vector{Int})
    len, ps, res = length(row), length(row) ÷ 2, 0
    grp = [Set{Int}() for i in 1:ps]
    for i in 1:2:len
        l, r = row[i] ÷ 2 + 1, row[i + 1] ÷ 2 + 1
        (l != r) && (push!(grp[l], r); push!(grp[r], l))
    end
    visited = fill(false, ps)
    q = Queue{Int}()
    for i in 1:ps
        visited[i] && continue
        enqueue!(q, i)
        visited[i] = true
        cnt = 0
        while !isempty(q)
            frt = dequeue!(q)
            cnt += 1
            for neighbor in grp[frt]
                visited[neighbor] && continue
                enqueue!(q, neighbor)
                visited[neighbor] = true
            end
        end
        res += cnt - 1
    end
    res
end
## @lc code=end


text:
using MAT
using Regression
using DecisionTree

type PredictionCenter
	TEST_MODE_MIXTURE_TPS #default value:0
	TEST_MODE_DATASET #default value:1
	NUM_TPS_SAMPLES  #default value:100
	testVar
	taskName
	workloadName
	lockType
	learnLock
	whichTransactionToPlot
	ioConf
	lockConf

	trainConfig
	testMode
	testConfig
	testMixture
	testMinTPS
	testMaxTPS
	testWorkloadRatio

	testSampleTPS
	testSampleTransactionCount

	throttleLatencyType
	throttleTargetLatency
	throttleTargetTransactionIndex
	throttlePenalty
	throttleIndividualTransactions #default value: false;
end
PredictionCenter()=PredictionCenter(0,0,100,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,false)

type PredictionConfig
	io_conf
	lock_conf
	transactionType
	groupingStrategy
	datasetList #default value: {};
	initialized #default value: false;
	mv # this can be grouped or ungrouped depending on group parameters.
	mvUngrouped # this is always ungrouped.
	configSummary

	transactionCount # C
	averageCpuUsage # P
	diskWrite # IO
	transactionLatency # L
	transactionLatencyPercentile
	lockWaitTime # W, NumOfWaitDueToLocks
	currentLockWait # LocksBeingWaitedFor
	TPS # TPS
	rowsChanged # RowsChanged
	pagesFlushed # PagesFlushed
	transactionMixture # Mixture
	logicalReads
	physicalReads
	physicalReadsMB
	networkSendKB
	networkRecvKB
	logWriteMB # LogIOw

	transactionCountUngrouped # C
	averageCpuUsageUngrouped # P
	diskWriteUngrouped # IO
	transactionLatencyUngrouped # L
	transactionLatencyPercentileUngrouped
	lockWaitTimeUngrouped # W, NumOfWaitDueToLocks
	currentLockWaitUngrouped # LocksBeingWaitedFor
	TPSUngrouped # TPS
	rowsChangedUngrouped # RowsChanged
	pagesFlushedUngrouped # PagesFlushed
	transactionMixtureUngrouped # Mixture
	logicalReadsUngrouped
	physicalReadsUngrouped
	physicalReadsMBUngrouped
	networkSendKBUngrouped
	networkRecvKBUngrouped
	logWriteMBUngrouped # LogIOw
end
PredictionConfig()=PredictionConfig(0,0,0,0,{},false,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)

type typeOfPredictions
	R
	T_total
	M_total
	Vp
	V
	W
	Pcon
	totalWaits
	TimeSpentWaiting
	LocksBeingHeld
end

type typeOfFlushRate
	io_conf
	workloadName
end

type typeOfConcurrencyLevel
	lock_conf
	workloadName
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
## Relational Disjunctive Domain A1 ∨ A2 ∨ A3 ∨ ... ∨ An
## =================================================================

# This is a simple non-relational domain of disjunctive domains (paramaterised by T)
immutable Relational{T} <: Domain{T}
  values::Set{T}
end

# This domain is relational, duh
isrelational(::Union{Relational, Type{Relational}}) = true

## Set Operations
## ==============
issubset(x::SimpleDisjunctive, y::SimpleDisjunctive) = issubset(x.values, y.values)
isintersect(x::SimpleDisjunctive, y::SimpleDisjunctive) = intersect(x.values, y.values)
domaineq(x::SimpleDisjunctive, y::SimpleDisjunctive) = x.values == y.values
⊔{T}(x::SimpleDisjunctive{T},y::T) = push!(x.values, y)

# Apply f to every abstract element in disjunction
function setmap{T}(f::Function, s::Set{T})
  out = Set{T}()
  for x in s
    push!(out,f(x))
  end
  out
end


text:
module Benchmarker

export cachescrub, gcscrub

using Statistics
using DelimitedFiles

include("Results.jl")
include("CSV.jl")

gcscrub() = (GC.gc(); GC.gc(); GC.gc(); GC.gc())

A = rand(parse(Int32, ENV["LAMP_L3_CACHE_SIZE"]))
cachescrub() = (A .+= 0.0001)

function measure(iters, f, args...)
  timings = Array{Float64}(undef, iters)

  local elapsed_time::Float64 = 0.0
  for i=1:iters
    elapsed_time, result = f(args...)
    timings[i] = elapsed_time
  end
  return Results(iters, timings)
end
end


text:
using SimpleARules
using Test
# include("../src/arules.jl")

transactions = [["bla","blubb"], ["bla"]]

@testset "some tests" begin
    unique_items_expected = Dict("bla"=>1, "blubb"=>2)
    unique_items = SimpleARules.get_unique_items(transactions)
    @test unique_items == unique_items_expected

    sorted_integer_transactions_expected = [[1,2], [1]]
    sorted_integer_transactions = SimpleARules.get_sorted_integer_transactions(transactions, unique_items_expected)
    @test sorted_integer_transactions == sorted_integer_transactions_expected

    support_dict_expected_1 = Dict([1]=>2, [2]=>1)
    support_dict = SimpleARules.get_support_dict(sorted_integer_transactions_expected, 1, 1)
    @test support_dict == support_dict_expected_1

    support_dict_expected_2 = Dict([1,2]=>1)
    support_dict = SimpleARules.get_support_dict(sorted_integer_transactions_expected, 2, 1)
    @test support_dict == support_dict_expected_2

    support_dictionaries_expected = [support_dict_expected_1, support_dict_expected_2]
    support_dictionaries = SimpleARules.get_support_dictionaries(sorted_integer_transactions_expected, 2, 1)
    @test support_dictionaries == support_dictionaries_expected

    frequent_dict_expected = Dict(["bla"]=>2 ,["blubb"] => 1, ["bla","blubb"]=>1)
    frequent_df = SimpleARules.frequent(transactions, max_length=2, min_support=1)
    @test Dict(zip(frequent_df.set,frequent_df.support)) == frequent_dict_expected

    rules_expected = [(["bla"], ["blubb"], 2, 1, 0.5), (["blubb"], ["bla"], 1, 1, 1.)]
    rules = SimpleARules.get_rules(support_dictionaries_expected, unique_items_expected)
    @test rules == rules_expected

    min_confidence=0.7
    rules_df = SimpleARules.get_rules_df(support_dictionaries_expected, min_confidence, unique_items_expected)
    @test minimum(rules_df.confidence) >= min_confidence

    apriori_df = SimpleARules.apriori(transactions, max_length=2, min_support=0, min_confidence=0)
    @test length(apriori_df[:, :subset]) == 2
    @test Vector(apriori_df[1, :]) == [["bla"], ["blubb"], 2, 1, 0.5]
    @test Vector(apriori_df[2, :]) == [["blubb"], ["bla"], 1, 1, 1]
end


text:
export interactive_orbitdiagram, scaleod

"""
    interactive_orbitdiagram(ds::DiscreteDynamicalSystem, p_index, pmin, pmax, i0 = 1;
                             u0 = get_state(ds), parname = "p", title = "")

Open an interactive application for exploring orbit diagrams (ODs) of discrete
dynamical systems. Requires `DynamicalSystems`.

Keywords control the name of the parameter, the initial state (used for _any_ parameter)
or whether to add a title above the orbit diagram.

## Interaction
The application is separated in the "OD plot" (left) and the "control panel" (right).
On the OD plot you can interactively click
and drag with the left mouse button to select a region in the OD. This region is then
**re-computed** at a higher resolution.

The options at the control panel are straight-forward, with
* `n` amount of steps recorded for the orbit diagram (not all are in the zoomed region!)
* `t` transient steps before starting to record steps
* `d` density of x-axis (the parameter axis)
* `α` alpha value for the plotted points.

Notice that at each update `n*t*d` steps are taken.
You have to press `update` after changing these parameters.
Press `reset` to bring the OD in the original
state (and variable). Pressing `back` will go back through the history of your exploration
History is stored when the "update" button is pressed or a region is zoomed in.

You can even decide which variable to get the OD for
by choosing one of the variables from the wheel!
Because the y-axis limits can't be known when changing variable, they reset to the size
of the selected variable.

## Accessing the data
What is plotted on the application window is a _true_ orbit diagram, not a plotting
shorthand. This means that all data are obtainable and usable directly.
Internally we always scale the orbit diagram to [0,1]² (to allow `Float64` precision
even though plotting is `Float32`-based). This however means that it is
necessary to transform the data in real scale. This is done through the function
[`scaleod`](@ref) which accepts the 5 arguments returned from the current function:
```julia
oddata = interactive_orbitdiagram(...)
ps, us = scaleod(oddata)
```
"""
function interactive_orbitdiagram(ds, p_index, p_min, p_max, i0 = 1;
    u0 = DynamicalSystems.get_state(ds), parname = "p", title = "")

    scene, layout = layoutscene(resolution = (1400, 600), backgroundcolor = DEFAULT_BG)
    display(scene)
    odax = layout[1,1] = Axis(scene)
    for z in (:xpanlock, :ypanlock, :xzoomlock, :yzoomlock)
        setproperty!(odax, z, true)
    end
    controllayout = layout[1, 2] = GridLayout(height = Auto(false))
    colsize!(layout, 1, Relative(3/5))
    display(scene)

    controllayout.tellheight[] = false
    odax.tellheight = true

    nslider, Tslider, dslider, n, Ttr, d, α, i,
    ▢update, ▢back, ▢reset, ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊ =
    add_controls!(controllayout, scene, DynamicalSystems.dimension(ds), parname, i0)

    # Initial Orbit diagram data
    integ = DynamicalSystems.integrator(ds, u0)
    p₋, p₊ = p_min, p_max
    odinit, xmin, xmax = minimal_normalized_od(integ, i[], p_index, p₋, p₊, d[], n[], Ttr[], u0)
    od_obs = Observable(odinit)
    dinit = d[]; ninit = n[]; Ttrinit = Ttr[]

    # History stores the variable index and true diagram limits
    history = [(i[], p₋, p₊, xmin, xmax, n[], Ttr[], d[])]
    update_controls!(history[end], i, n, Ttr, d,  ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊)

    color = lift(a -> RGBAf0(0,0,0,a), α)
    scatter!(odax, od_obs; marker = MARKER, markersize = 1px, color = color, strokewidth = 0.0)

    xlims!(odax, 0, 1)
    ylims!(odax, 0, 1)
    odax.xticks = (0:0.25:1, ["$(parname)₋", " ", " ", " ", "$(parname)₊"])
    odax.yticks = (0:0.25:1, ["u₋", " ", " ", " ", "u₊"])
    odax.ylabel = "u"*subscript(i[])
    on(i) do o; odax.ylabel = "u"*subscript(o); end
    odax.xlabel = parname
    MakieLayout.deactivate_interaction!(odax, :rectanglezoom)
    rect = select_rectangle(odax.scene)

    # # Uppon interactively selecting a rectangle, with value `r` (in [0,1]²)
    on(rect) do r
        sp₋, sxmin = r.origin
        sp₊, sxmax = r.origin + r.widths
        # Convert p,x to true values
        j, pp₋, pp₊, pxmin, pxmax = history[end]
        pdif = pp₊ - pp₋; xdif = pxmax - pxmin
        p₋ = sp₋*pdif + pp₋
        p₊ = sp₊*pdif + pp₋
        xmin = sxmin*xdif + pxmin
        xmax = sxmax*xdif + pxmin

        od_obs[] = minimal_normalized_od(
            integ, j,  p_index, p₋, p₊,
            d[], n[], Ttr[], u0, xmin, xmax
        )
        # update history and controls
        push!(history, (j, p₋, p₊, xmin, xmax, n[], Ttr[], d[]))
        update_controls!(history[end], i, n, Ttr, d,  ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊)
    end

    # Upon hitting the update button
    on(▢update) do clicks
        j, p₋, p₊, xmin, xmax, m, T, dens = history[end]
        # Check if there was any change:
        if !(⬜p₋[] == p₋ && ⬜p₊[] == p₊ && ⬜u₋[] == xmin && ⬜u₊[] == xmax &&
            j == i[] && m == n[] && T == Ttr[] && dens == d[])

            p₋, p₊, xmin, xmax = ⬜p₋[], ⬜p₊[], ⬜u₋[], ⬜u₊[]
            newj, m, T, dens = i[], n[], Ttr[], d[]

            if newj ≠ j # a new variable is selected, so x limits must be recomputed
                od_obs[], xmin, xmax = minimal_normalized_od(
                    integ, newj, p_index, p₋, p₊, d[], n[], Ttr[], u0
                )
            else
                od_obs[] = minimal_normalized_od(
                    integ, newj, p_index, p₋, p₊, d[], n[], Ttr[], u0, xmin, xmax
                )
            end
            # Update history and controls
            push!(history, (newj, p₋, p₊, xmin, xmax, m, T, dens))
            update_controls!(history[end], i, n, Ttr, d,  ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊)
        end
    end

    # Upon hitting the "reset" button
    on(▢reset) do clicks
        if length(history) > 1
            deleteat!(history, 2:length(history))
            od_obs[] = odinit
            update_controls!(history[end], i, n, Ttr, d,  ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊)
        end
    end

    # Upon hitting the "back" button
    on(▢back) do clicks
        if length(history) > 1
            pop!(history)
            j, p₋, p₊, xmin, xmax, m, T, dens = history[end]
            od_obs[] = minimal_normalized_od(
                integ, j, p_index, p₋, p₊, dens, m, T, u0, xmin, xmax
            )
            update_controls!(history[end], i, n, Ttr, d,  ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊)
        end
    end

    # for the following two buttons the slider values must be updated
    onany(▢reset, ▢back) do val1, val2
        j, p₋, p₊, xmin, xmax, m, T, dens = history[end]
        set_close_to!(nslider, m)
        set_close_to!(Tslider, T)
        set_close_to!(dslider, dens)
    end

    if title ≠ ""
        layout[0, 1] = Label(scene, title, textsize = 30)
    end
    MakieLayout.trim!(layout)

    return od_obs, ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊
end

"""
    minimal_normalized_od(integ, i, p_index, p₋, p₊, d, n, Ttr, u0) → od, xmin, xmax
    minimal_normalized_od(integ, i, p_index, p₋, p₊, d, n, Ttr, u0, xmin, xmax) → od

Compute and return a minimal and normalized orbit diagram (OD).

All points are stored in a single vector of `Point2f0` to ensure fastest possible
plotting. In addition all numbers are scaled to [0, 1]. This allows us to have
64-bit precision while display is only 32-bit!

The version with `xmin, xmax` only keeps points with limits between the
real `xmin, xmax` (in the normal units of the dynamical system).
"""
function minimal_normalized_od(integ, i, p_index, p₋, p₊,
                              d::Int, n::Int, Ttr::Int, u0)

    pvalues = range(p₋, stop = p₊, length = d)
    pdif = p₊ - p₋
    od = Vector{Point2f0}() # make this pre-allocated
    xmin = eltype(integ.u)(Inf); xmax = eltype(integ.u)(-Inf)
    @inbounds for (j, p) in enumerate(pvalues)
        pp = (p - p₋)/pdif # p to plot, in [0, 1]
        DynamicalSystems.reinit!(integ, u0)
        integ.p[p_index] = p
        DynamicalSystems.step!(integ, Ttr)
        for z in 1:n
            DynamicalSystems.step!(integ)
            x = integ.u[i]
            push!(od, Point2f0(pp, integ.u[i]))
            # update limits
            if x < xmin
                xmin = x
            elseif x > xmax
                xmax = x
            end
        end
    end
    # normalize x values to [0, 1]
    xdif = xmax - xmin
    @inbounds for j in eachindex(od)
        x = od[j][2]; p = od[j][1]
        od[j] = Point2f0(p, (x - xmin)/xdif)
    end
    return od, xmin, xmax
end

function minimal_normalized_od(integ, i, p_index, p₋, p₊,
                              d::Int, n::Int, Ttr::Int, u0, xmin, xmax)

    pvalues = range(p₋, stop = p₊, length = d)
    pdif = p₊ - p₋; xdif = xmax - xmin
    od = Vector{Point2f0}()
    @inbounds for p in pvalues
        pp = (p - p₋)/pdif # p to plot, in [0, 1]
        DynamicalSystems.reinit!(integ, u0)
        integ.p[p_index] = p
        DynamicalSystems.step!(integ, Ttr)
        for z in 1:n
            DynamicalSystems.step!(integ)
            x = integ.u[i]
            if xmin ≤ x ≤ xmax
                push!(od, Point2f0(pp, (integ.u[i] - xmin)/xdif))
            end
        end
    end
    return od
end

function  update_controls!(h, i, n, Ttr, d, ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊)
    j, p₋, p₊, xmin, xmax, m, T, dens = h
    i[] = j; n[] = m; Ttr[] = T; d[] = dens
    ⬜p₋[] = p₋; ⬜p₊[] = p₊
    ⬜u₋[] = xmin; ⬜u₊[] = xmax
    return
end

"""
    scaleod(oddata) -> ps, us
Given the return values of [`interactive_orbitdiagram`](@ref), produce
orbit diagram data scaled correctly in data units. Return the data as a vector of
parameter values and a vector of corresponding variable values.
"""
scaleod(r::Tuple) = scaleod(r...)
function scaleod(od, p₋, p₊, u₋, u₊)
    oddata = od[]; L = length(oddata);
    T = promote_type(typeof(u₋[]), Float32)
    ps = zeros(T, L); us = copy(ps)
    udif = u₊[] - u₋[]; um = u₋[]
    pdif = p₊[] - p₋[]; pm = p₋[]
    @inbounds for i ∈ 1:length(oddata)
        p, u = oddata[i]
        ps[i] = pm + pdif*p; us[i] = um + udif*u
    end
    return ps, us
end

# TODO: Allow initial state to be a function of paramter (define function `get_u(f, p)`)

function observable_slider!(layout, i, j, scene, Label, r;
    wl = 40, wr = nothing, startvalue = r[1])
    slider = Slider(scene, range = r, startvalue = startvalue)
    text_prev = Label(scene, "$Label =", halign = :right)
    text_after = Label(scene, lift(a -> "$(string(a))", slider.value),
    halign = :left)
    layout[i, j] = hbox!(text_prev, slider, text_after)
    return slider
end

function od_sliders!(scene, controllayout)
    sliders = []
    for (i, (l, vals)) in enumerate(zip(("n =", "t =", "d =", "α ="),
                         (1000:100:100000, 1000:1000:100000, 100:100:10000, 0.001:0.001:1)))

        startvalue = l[1] == 'α' ? 0.1 : l[1] == 'd' ? 1000 : vals[1]
        sll = labelslider!(scene, l, vals; sliderkw = Dict(:startvalue => startvalue))
        push!(sliders, sll.slider)
        controllayout[i, :] = sll.layout
    end
    return sliders
end

function add_controls!(controllayout, scene, D, parname, i0)
    # Sliders
    # nslider = observable_slider!(controllayout, 1, :, scene, "n", 1000:100:100000)
    # Tslider = observable_slider!(controllayout, 2, :, scene, "t", 1000:1000:100000)
    # dslider = observable_slider!(controllayout, 3, :, scene, "d", 100:100:10000; startvalue = 1000)
    # αslider = observable_slider!(controllayout, 4, :, scene, "α", 0.001:0.001:1; startvalue = 0.1)
    nslider, Tslider, dslider, αslider = od_sliders!(scene, controllayout)
    # Buttons (incl. variable chooser)
    ▢update = Button(scene, label = "update")
    ▢back = Button(scene, label = "← back")
    ▢reset = Button(scene, label = "reset")
    imenu = Menu(scene, options = [string(j) for j in 1:D], width = 60)
    imenu.i_selected = i0
    controllayout[5, 1] = hbox!(
        ▢update, ▢back, ▢reset,
        Label(scene, "variable:"), imenu, width = Auto(false)
    )
    # Limit boxes. Unfortunately can't be made observables yet...
    ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊ = Observable.((0.0, 1.0, 0.0, 1.0))
    tsize = 16
    text_p₋ = Label(scene, lift(o -> "$(parname)₋ = $(o)", ⬜p₋),
        halign = :left, width = Auto(false), textsize = tsize)
    text_p₊ = Label(scene, lift(o -> "$(parname)₊ = $(o)", ⬜p₊),
        halign = :left, width = Auto(false), textsize = tsize)
    text_u₋ = Label(scene, lift(o -> "u₋ = $(o)", ⬜u₋),
        halign = :left, width = Auto(false), textsize = tsize)
    text_u₊ = Label(scene, lift(o -> "u₊ = $(o)", ⬜u₊),
        halign = :left, width = Auto(false), textsize = tsize)
    controllayout[6, 1] = grid!([text_p₋ text_p₊ ; text_u₋ text_u₊])
    ⬜p₋[], ⬜p₊[], ⬜u₋[], ⬜u₊[] = rand(4)
    return nslider, Tslider, dslider,
           nslider.value, Tslider.value, dslider.value, αslider.value,
           imenu.i_selected, ▢update.clicks, ▢back.clicks, ▢reset.clicks,
           ⬜p₋, ⬜p₊, ⬜u₋, ⬜u₊
end


text:
# simple program to test affinity propagation

using Test
using Distances
using Clustering
using LinearAlgebra
using Random
using Statistics
include("test_helpers.jl")

@testset "affinityprop() (affinity propagation)" begin
    @testset "Argument checks" begin
        Random.seed!(34568)
        @test_throws ArgumentError affinityprop(randn(2, 3))
        @test_throws ArgumentError affinityprop(randn(1, 1))
        @test_throws ArgumentError affinityprop(randn(2, 2), tol=0.0)
        @test_throws ArgumentError affinityprop(randn(2, 2), damp=-0.1)
        @test_throws ArgumentError affinityprop(randn(2, 2), damp=1.0)
        @test affinityprop(randn(2, 2), damp=0.5, tol=0.5) isa AffinityPropResult
        for disp in keys(Clustering.DisplayLevels)
            @test affinityprop(randn(2, 2), tol=0.1, display=disp) isa AffinityPropResult
        end
    end

    Random.seed!(34568)

    d = 10
    n = 500
    x = rand(d, n)
    S = -pairwise(Euclidean(), x, x, dims=2)

    # set diagonal value to median value
    S = S - diagm(0 => diag(S)) + median(S)*I

    R = affinityprop(S)

    @test isa(R, AffinityPropResult)
    k = length(R.exemplars)
    @test 0 < k < n
    @test length(R.assignments) == n
    @test all(R.assignments .>= 1) && all(R.assignments .<= k)
    @test all(R.assignments[R.exemplars] .== collect(1:k))

    @test length(R.counts) == k
    @test sum(R.counts) == n
    for i = 1:k
        @test R.counts[i] == count(==(i), R.assignments)
    end

    @testset "Support for arrays other than Matrix{T}" begin
        @testset "$(typeof(M))" for M in equivalent_matrices(S)
            R2 = affinityprop(M)
            @test R2.assignments == R.assignments
        end
    end

    #= compare with python result
    the reference assignments were computed using python sklearn:
    ```julia
    using PyCall

    @pyimport sklearn.cluster as cl
    af = cl.AffinityPropagation(affinity="precomputed")[:fit]( S )

    ref_assignments = af[:labels_] .+ 1
    ```
    =#
    ref_assignments = [11, 47, 23, 1, 3, 8, 54, 2, 6, 29, 11, 43, 2, 20, 13, 40,
                       3, 48, 3, 32, 19, 59, 4, 28, 59, 24, 5, 36, 6, 23, 49, 19,
                       42, 14, 40, 29, 20, 49, 53, 48, 51, 55, 57, 32, 32, 10, 48,
                       12, 13, 7, 1, 57, 52, 7, 3, 47, 43, 10, 53, 40, 52, 8, 42,
                       15, 23, 47, 32, 9, 10, 42, 52, 43, 49, 32, 49, 12, 5, 53,
                       54, 58, 14, 21, 7, 16, 14, 17, 12, 15, 38, 15, 33, 11, 27,
                       42, 56, 14, 42, 41, 34, 26, 55, 46, 47, 9, 25, 23, 19, 47,
                       34, 38, 53, 10, 9, 9, 43, 59, 43, 51, 12, 39, 29, 12, 13,
                       9, 52, 1, 13, 36, 11, 21, 58, 38, 2, 31, 18, 19, 16, 14, 6,
                       50, 15, 16, 21, 26, 17, 21, 43, 25, 18, 39, 57, 19, 41, 4,
                       31, 20, 15, 16, 34, 10, 28, 45, 21, 13, 13, 21, 43, 2, 46,
                       52, 12, 26, 21, 14, 52, 21, 36, 6, 22, 25, 45, 25, 12, 1, 3,
                       25, 31, 55, 28, 32, 9, 32, 44, 33, 6, 17, 2, 58, 32, 38, 20,
                       3, 12, 51, 23, 51, 1, 32, 12, 35, 13, 39, 48, 40, 7, 57, 4,
                       38, 51, 57, 24, 20, 30, 25, 28, 55, 32, 2, 28, 26, 27, 23,
                       46, 46, 9, 12, 49, 52, 17, 24, 44, 23, 54, 46, 9, 28, 32,
                       20, 27, 22, 57, 60, 61, 59, 56, 41, 8, 26, 20, 12, 36, 26,
                       23, 48, 17, 29, 19, 41, 20, 4, 29, 55, 43, 11, 24, 52, 42,
                       30, 37, 32, 11, 59, 6, 53, 59, 33, 52, 11, 31, 3, 52, 36,
                       34, 40, 23, 48, 50, 16, 32, 52, 55, 47, 56, 25, 13, 12, 33,
                       35, 20, 57, 61, 47, 40, 29, 31, 34, 37, 19, 2, 59, 57, 35,
                       44, 50, 30, 36, 27, 25, 50, 5, 38, 61, 8, 60, 23, 60, 37,
                       9, 6, 9, 25, 33, 23, 27, 55, 25, 7, 38, 55, 12, 44, 29, 7,
                       14, 34, 56, 39, 43, 14, 27, 17, 29, 51, 38, 40, 41, 42,
                       16, 42, 14, 19, 44, 43, 34, 60, 44, 20, 3, 45, 61, 57, 33,
                       46, 44, 40, 22, 8, 12, 43, 46, 20, 26, 43, 7, 61, 47, 59,
                       38, 48, 1, 49, 48, 35, 22, 50, 44, 43, 1, 53, 50, 52, 53,
                       41, 49, 6, 51, 51, 40, 5, 12, 20, 58, 57, 16, 1, 57, 58,
                       34, 17, 52, 25, 44, 20, 41, 50, 58, 53, 47, 58, 14, 40, 39,
                       10, 53, 1, 34, 54, 53, 49, 60, 31, 54, 35, 39, 23, 10, 55,
                       12, 56, 57, 48, 57, 15, 35, 23, 13, 10, 22, 4, 25, 17, 47,
                       5, 36, 51, 6, 44, 26, 27, 40, 58, 59, 15, 55, 19, 21, 7,
                       32, 45, 13, 41, 47, 30, 36, 1, 60, 43, 48, 61, 51]
    @test randindex(R.assignments, ref_assignments)[2] == 1.0
end


text:
export EuclidDistance, PackedEuclidDistance

"""
$(TYPEDEF)

Default linear offset between two scalar variables.
"""
struct EuclidDistance{T <: SamplableBelief} <: AbstractRelativeMinimize
  Z::T
end

EuclidDistance(::UniformScaling=LinearAlgebra.I) = EuclidDistance(Normal())

# consider a different group?
getManifold(::InstanceType{EuclidDistance}) = TranslationGroup(1)
getDimension(::InstanceType{<:EuclidDistance}) = 1

function getSample(cf::CalcFactor{<:EuclidDistance})
  rand(cf.factor.Z, 1)
end

# new and simplified interface for both nonparametric and parametric
function (s::CalcFactor{<:EuclidDistance})(z, x1, x2)
  # @info "distance?" z x1 x2
  return z .- norm(x2 .- x1)
end

Base.convert(::Type{<:MB.AbstractManifold}, ::InstanceType{EuclidDistance}) = Manifolds.TranslationGroup(1)

"""
$(TYPEDEF)
Serialization type for `EuclidDistance` binary factor.
"""
mutable struct PackedEuclidDistance <: PackedInferenceType
  _type::String
  Z::String
end

function convert(::Type{PackedEuclidDistance}, d::EuclidDistance)
  PackedEuclidDistance( "/application/JuliaLang/PackedSamplableBelief",
                        convert(PackedSamplableBelief, d.Z) )
end

function convert(::Type{<:EuclidDistance}, d::PackedEuclidDistance)
  EuclidDistance(convert(SamplableBelief, d.Z))
end

#


text:
#Implements a Matched Filter for BEAMFORMING

mutable struct MatchedFilter
    filter::Array
    chirpF::Array{Complex{Float64}}
    nFFT::Int
    nPhones::Int
    MatchedFilter() = new()
end
function (mf::MatchedFilter)(dataIn::Array,dataOut::Array)
    if size(dataOut)!=(mf.nFFT,mf.nPhones) error("Data Out has the wrong dims") end
    dataOut[1:size(dataIn,1),:] = dataIn
    fft!(dataOut,1)
    dataOut ./= (abs.(dataOut) .+ 1e-9) # OPTIONAL PHAT transform
    dataOut .*= mf.filter
    ifft!(dataOut,1)
    nothing
end

function prepMF(chirpIn::Array,n::Int,m::Int)
    if n<length(chirpIn)
        error("FFT range smaller than chirp")
    end
    thisMF = MatchedFilter()
    thisMF.nFFT=n
    thisMF.nPhones=m
    thisMF.chirpF = zeros(Complex{Float64},n)
    thisMF.chirpF[1:length(chirpIn)] = chirpIn
    conj!(fft!(thisMF.chirpF))
    thisMF.filter = repeat(thisMF.chirpF,1,m) # repmat
    return thisMF
end

#


text:
using GAIO
using StaticArrays
using HostCPUFeatures
using Test

@testset "exported functionality" begin
    f(x) = x .^ 2
    test_points = [
        (-1.0, -1.0), (-1.0, 1.0), (1.0, -1.0), (1.0, 1.0),
        (-1.0,  0.0), (1.0, 0.0),  (0.0, -1.0), (0.0, 1.0)
    ]
    center = SVector(0.0, 0.0)
    radius = SVector(1.0, 1.0)
    domain = Box(center, radius)
    g = PointDiscretizedMap(f, domain, test_points)
    @testset "basics" begin
        @test typeof(g) <: SampledBoxMap
        partition = BoxPartition(domain, (32,32))
        p1 = (0.0, 0.0)
        p2 = (0.5, 0.0)
        p3 = (0.0, -0.5)
        boxset = partition[(p1, p2, p3)]
        # re-implement a straight forward box map to have a ground truth
        # terrible implementation in any real scenario
        mapped_points = []
        for box in boxset
            push!(mapped_points, f(box.center .- box.radius))
            push!(mapped_points, f(box.center .+ box.radius))
            x = f((box.center[1] .- box.radius[1], box.center[2] .+ box.radius[2]))
            push!(mapped_points, x)
            y = f((box.center[1] .+ box.radius[1], box.center[2] .- box.radius[2]))
            push!(mapped_points, y)
        end
        image = partition[mapped_points]
        full_boxset = partition[:]
        mapped1 = g(boxset)

        @test !isempty(boxset) && !isempty(image)
        # easiest way right now to check for equality
        @test length(union(image, mapped1)) == length(intersect(image, mapped1))
    end
    g = PointDiscretizedMap(f, domain, test_points, Val(:cpu))
    @testset "basics with :cpu" begin
        @test typeof(g) <: SampledBoxMap
        partition = BoxPartition(domain, (32,32))
        p1 = SVector(0.0, 0.0)
        p2 = SVector(0.5, 0.0)
        p3 = SVector(0.0, -0.5)
        boxset = partition[(p1, p2, p3)]
        # re-implement a straight forward box map to have a ground truth
        # terrible implementation in any real scenario
        mapped_points = []
        for box in boxset
            push!(mapped_points, f(box.center .- box.radius))
            push!(mapped_points, f(box.center .+ box.radius))
            x = f(SVector(box.center[1] .- box.radius[1], box.center[2] .+ box.radius[2]))
            push!(mapped_points, x)
            y = f(SVector(box.center[1] .+ box.radius[1], box.center[2] .- box.radius[2]))
            push!(mapped_points, y)
        end
        image = partition[mapped_points]
        full_boxset = partition[:]
        mapped1 = g(boxset)

        @test !isempty(boxset) && !isempty(image)
        # easiest way right now to check for equality
        @test length(union(image, mapped1)) == length(intersect(image, mapped1))
    end
    g = PointDiscretizedMap(f, domain, test_points, Val(:cpu))
    @testset "basics with SIMD" begin
        @test typeof(g) <: SampledBoxMap
        partition = BoxPartition(domain, (32,32))
        p1 = SVector(0.0, 0.0)
        p2 = SVector(0.5, 0.0)
        p3 = SVector(0.0, -0.5)
        boxset = partition[(p1, p2, p3)]
        # re-implement a straight forward box map to have a ground truth
        # terrible implementation in any real scenario
        mapped_points = []
        for box in boxset
            push!(mapped_points, f(box.center .- box.radius))
            push!(mapped_points, f(box.center .+ box.radius))
            x = f(SVector(box.center[1] - box.radius[1], box.center[2] + box.radius[2]))
            push!(mapped_points, x)
            y = f(SVector(box.center[1] + box.radius[1], box.center[2] - box.radius[2]))
            push!(mapped_points, y)
        end
        image = partition[mapped_points]
        full_boxset = partition[:]
        mapped1 = g(boxset)

        @test !isempty(boxset) && !isempty(image)
        # easiest way right now to check for equality
        @test length(union(image, mapped1)) == length(intersect(image, mapped1))
    end
    @testset "points in boxmaps" begin
        x = (-2.0, 3.0)
        y = (4.0, 1.0)
        @test_throws MethodError g(x)
        @test_throws MethodError g(y)
    end
end


text:
@testset "TriangleIterator" begin
    using IntrinsicPersistence.equilaterals

    n = 100
    nrepeats = 10

    @testset "ordertoedge" begin
        using IntrinsicPersistence.ordertoedge

        for dim in 2:nrepeats+1
            pts = rand(dim, n)
            M   = pairwise(Euclidean(), pts)
            ord = ordertoedge(M)

            @test issorted(M[i, j] for (i, j) in ord)
            @test all(i > j for (i, j) in ord)
        end
    end

    @testset "type stability and iterator traits" begin
        using IntrinsicPersistence.Triangle

        for T in [Float64, Float32, Int64, Int32]
            if T <: Integer
                pts = rand(T, 3, n) .÷ T(16)
            else
                pts = rand(T, 3, n)
            end
            M = pairwise(Cityblock(), pts) # preserves Ints.

            @test_inferred equilaterals(M)
            Δs = equilaterals(M)
            @test_inferred start(Δs)
            st = start(Δs)
            @test_inferred done(Δs, st)
            @test_inferred next(Δs, st)
            @test_inferred first(Δs)

            @test Base.iteratoreltype(Δs) == Base.HasEltype()
            @test eltype(Δs) == Triangle{T}
            @test typeof(first(Δs)) == eltype(Δs)
            @test Base.iteratorsize(Δs) == Base.SizeUnknown()
        end
    end

    @testset "unique, sorted and equilateral?" begin
        using IntrinsicPersistence.sampledensity

        for dim in 2:nrepeats+1
            pts = rand(dim, n)
            M   = pairwise(Euclidean(), pts)
            Δs  = collect(equilaterals(M))
            s   = sampledensity(M)

            # Sensible length.
            @test length(Δs) < binomial(n, 3)

            # Unique.
            @test Δs == unique(Δs)

            # Sorted.
            for i in 1:length(Δs)-1
                @test Δs[i].r ≤ Δs[i+1].r
            end

            # Equilateral.
            for Δ in Δs
                i, j, k = Δ.vertices
                a = M[i, j]
                b = M[j, k]
                c = M[k, i]

                @test isapprox(a, b, atol = 2s) &&
                      isapprox(b, c, atol = 2s) &&
                      isapprox(c, a, atol = 2s)
                @test max(a, b, c) == Δ.r
            end
        end
    end

    @testset "counts" begin
        h = √3/2
        #   *
        #  * *
        # * * *
        pts = [0 1 2 0.5 1.5 1;
               0 0 0 h   h   2h]
        M  = pairwise(Euclidean(), pts)
        # General position:
        M .+= 0.001 * Symmetric(rand(size(M)))
        Δs = collect(equilaterals(M, tol = 0.1))

        @test Δs == unique(Δs)
        @test length(Δs) == 5

        #   * * *
        #  * * *
        # * * *
        h   = √3/2
        pts = [0 1 2 0.5 1.5 2.5 1  2  3;
               0 0 0 h   h   h   2h 2h 2h]
        M  = pairwise(Euclidean(), pts)
        # General position:
        M .+= 0.001 * Symmetric(rand(size(M)))
        Δs = collect(equilaterals(M, tol = 0.1))

        @test Δs == unique(Δs)
        @test length(Δs) == 12

        #  * *
        # *   *
        #  * *
        pts = [1 2 2.5 2  1  0.5;
               0 0 h   2h 2h h]
        M  = pairwise(Euclidean(), pts)
        # General position:
        M .+= 0.001 * Symmetric(rand(size(M)))
        Δs = collect(equilaterals(M, tol = 0.1))
        @test Δs == unique(Δs)
        @test length(Δs) == 2

        # All points the same distance apart with noise.
        M = ones(n, n) + 0.001 * Symmetric(rand(n, n))
        for i in 1:n
            M[i, i] = 0
        end
        Δs = collect(equilaterals(M, tol = 0.1))

        @test Δs == unique(Δs)
        @test length(Δs) == binomial(n, 3)

        # All points the same distance apart without noise --
        # each triangle is returned three times.
        M = ones(n, n)
        for i in 1:n
            M[i, i] = 0
        end
        Δs = collect(equilaterals(M, tol = 0))

        @test length(Δs) == binomial(n, 3) * 3
    end
end


text:
import Base.length
import Base.push!

using Random

"""
Replay buffer implementation based on https://github.com/openai/baselines/blob/master/baselines/deepq/replay_buffer.py
"""
mutable struct ReplayBuffer
    size
    storage
    next_idx
end

"""
Constructor
"""
ReplayBuffer(size) = ReplayBuffer(size, [], 1)

"""
Length of the buffer
"""
length(buf::ReplayBuffer) = length(buf.storage)

"""
Tells whether there is enough data in the buffer or not
"""
can_sample(buf::ReplayBuffer, batch_size) = batch_size <= length(buf)

"""
Add new transition to the buffer
"""
function push!(buf::ReplayBuffer, obs_t, action, reward, obs_tp1, done)
    data = (obs_t, action, reward, obs_tp1, done)

    if buf.next_idx > length(buf)
        push!(buf.storage, data)
    else
        buf.storage[buf.next_idx] = data
    end

    buf.next_idx = max(1, (buf.next_idx + 1) % (buf.size+1))
end

function encode_sample(buf::ReplayBuffer, idxes; stack=1)
    bs = length(idxes)
    inpsize = size(buf.storage[1][1])[1:(end-1)] # get the dimensions of the input
    stacksize = nothing
    if length(inpsize) < 3
        stacksize = (inpsize[1]*stack, inpsize[2:end]...)
    else
        stacksize = (inpsize[1], inpsize[2], inpsize[3]*4, inpsize[4:end]...)
    end

    obses_t = zeros(Float32, stacksize..., bs)
    actions = zeros(Int, bs)
    rewards = zeros(Float32, bs)
    obses_tp1 = zeros(Float32, stacksize..., bs)
    dones = zeros(Float32, 1, bs)

    indx = 0
    for ind in idxes
        indx += 1
        for i=-(stack-1):0
            curr = ind + i
            curr = curr < 1 ? 1 : curr
            data = buf.storage[curr] #use available data instead of filling the history with 0

            obs_t, action, reward, obs_tp1, done = data
            if length(inpsize) < 3
                obses_t[map(t->(t*(i-1+stack)+1):(t*(i-1+stack)+t), size(obs_t)[1:(end-1)])..., indx] = obs_t #stack frames
                obses_tp1[map(t->(t*(i-1+stack)+1):(t*(i-1+stack)+t), size(obs_t)[1:(end-1)])..., indx] = obs_tp1 #stack frames
            else
                obses_t[map(t->1:t, size(obs_t)[1:(end-2)])..., i+stack, indx] = obs_t #stack frames
                obses_tp1[map(t->1:t, size(obs_t)[1:(end-2)])..., i+stack, indx] = obs_tp1 #stack frames
            end
            if i == 0
                actions[indx] = action
                rewards[indx] = reward
                dones[1, indx] = done ? 0.0 : 1.0
            end
        end
    end
    obses_t, actions, rewards, obses_tp1, dones
end

function encode_recent(buf::ReplayBuffer, obs_t; stack=1)
    if length(buf) != 0
        recent_idx = max(1, (buf.next_idx - 1) % (buf.size+1))
        obses_t, _, _, _,_ = encode_sample(buf, [recent_idx]; stack=stack)
        inpsize = size(obs_t)[1:(end-1)] # get the dimensions of the input
        i = 0
        if length(inpsize) < 3
            obses_t[map(t->(t*(i-1+stack)+1):(t*(i-1+stack)+t), size(obs_t)[1:(end-1)])..., 1] = obs_t #stack frames
        else
            obses_t[map(t->1:t, size(obs_t)[1:(end-2)])..., i+stack, 1] = obs_t #stack frames
        end
    else
        bs = 1
        inpsize = size(obs_t)[1:(end-1)] # get the dimensions of the input
        stacksize = nothing
        if length(inpsize) < 3
            stacksize = (inpsize[1]*stack, inpsize[2:end]...)
        else
            stacksize = (inpsize[1], inpsize[2], inpsize[3]*4, inpsize[4:end]...)
        end

        obses_t = zeros(Float32, stacksize..., bs)
        for i=-(stack-1):0
            if length(inpsize) < 3
                obses_t[map(t->(t*(i-1+stack)+1):(t*(i-1+stack)+t), size(obs_t)[1:(end-1)])..., 1] = obs_t #stack frames
            else
                obses_t[map(t->1:t, size(obs_t)[1:(end-2)])..., i+stack, 1] = obs_t #stack frames
            end
        end
    end
    return obses_t
end

function sample_batch(buf::ReplayBuffer, batchsize; stack=1)
    idxes = randperm(length(buf))[1:batchsize]
    return encode_sample(buf, idxes; stack=stack)
end


text:
using Test
import Pluto: Notebook, ServerSession, ClientSession, Cell, load_notebook, load_notebook_nobackup, save_notebook, WorkspaceManager, cutename, numbered_until_new, readwrite, without_pluto_file_extension
import Random
import Pkg
import UUIDs: UUID

# We define some notebooks explicitly, and not as a .jl notebook file, to avoid circular reasoning 🤔
function basic_notebook()
    Notebook([
        Cell("100*a + b"),
        Cell("a = 1"),
        Cell("💩 = :💩"), # ends with 4-byte character
        Cell("b = let\n\tx = a + a\n\tx*x\nend"),
        Cell("html\"
Hoi!
\n
My name is kiki
\""),
        # test included Markdown import
        Cell("""md"# Cześć!
        My name is **baba** and I like \$\\LaTeX\$ _support!_

        \$\$\\begin{align}
        \\varphi &= \\sum_{i=1}^{\\infty} \\frac{\\left(\\sin{x_i}^2 + \\cos{x_i}^2\\right)}{i^2} \\\\
        b &= \\frac{1}{2}\\,\\log \\exp{\\varphi}
        \\end{align}\$\$

        ### The spectacle before us was indeed sublime.
        Apparently we had reached a great height in the atmosphere, for the sky was a dead black, and the stars had ceased to twinkle. By the same illusion which lifts the horizon of the sea to the level of the spectator on a hillside, the sable cloud beneath was dished out, and the car seemed to float in the middle of an immense dark sphere, whose upper half was strewn with silver. Looking down into the dark gulf below, I could see a ruddy light streaming through a rift in the clouds."
        """),
        # test included InteractiveUtils import
        Cell("subtypes(Number)"),
    ]) |> init_packages!
end

function shuffled_notebook()
    Notebook([
        Cell("z = y"),
        Cell("v = u"),
        Cell("y = x"),
        Cell("x = w"),
        Cell("using Dates"),
        Cell("t = 1"),
        Cell("w = v"),
        Cell("u = t"),
    ]) |> init_packages!
end

function shuffled_with_imports_notebook()
    Notebook([
        Cell("c = uuid1()"),
        Cell("a = (b, today())"),
        Cell("y = 2"),
        Cell("using UUIDs"),
        Cell("y"),
        Cell("x = y"),
        Cell("b = base64encode"),
        Cell("""
        begin
            using Dates
            using Base64
        end"""),
        Cell("BasicREPL"),
        Cell("""
        begin
            x
            using REPL
        end"""),
    ]) |> init_packages!
end

function init_packages!(nb::Notebook)
    nb.topology = Pluto.updated_topology(nb.topology, nb, nb.cells)
    Pluto.sync_nbpkg_core(nb)
    return nb
end

function bad_code_notebook()
    Notebook([
        Cell("z = y"),
        Cell("y = z"),
        Cell(""";lka;fd;jasdf;;;\n\n\n\n\nasdfasdf

        [[["""),
        Cell("using Aasdfdsf"),
    ]) |> init_packages!
end

function bonds_notebook()
    Notebook([
        Cell("y = x"),
        Cell("@bind x html\"<input type='range'>\""),
        Cell("@assert y === missing"),
        Cell("""struct Wow
            x
        end"""),
        Cell("Base.get(w::Wow) = w.x"),
        Cell("Base.show(io::IO, ::MIME\"text/html\", w::Wow) = nothing"),
        Cell("w = Wow(10)"),
        Cell("@bind z w"),
        Cell("@assert z == 10"),
    ]) |> init_packages!
end

function project_notebook()
    Notebook([
        Cell("using Dates"),
        Cell("using Example"),
    ]) |> init_packages!
end

@testset "Notebook Files" begin
    nbs = [String(nameof(f)) => f() for f in [basic_notebook, shuffled_notebook, shuffled_with_imports_notebook, bad_code_notebook, bonds_notebook, project_notebook]]

    @testset "Sample notebooks " begin
        # Also adds them to the `nbs` list
        for file in ["Basic.jl", "Tower of Hanoi.jl", "Interactivity.jl"]
            path = normpath(Pluto.project_relative_path("sample", file))

            @testset "$(file)" begin
                nb = @test_nowarn load_notebook_nobackup(path)
                @test length(nb.cells) > 0
                push!(nbs, "sample " * file => nb)
            end
        end
    end

    🍭 = ServerSession()
    for (name, nb) in nbs
        nb.path = tempname() * "é🧡💛.jl"

        client = ClientSession(Symbol("client", rand(UInt16)), nothing)
        client.connected_notebook = nb

        🍭.connected_clients[client.id] = client
    end

    @testset "I/O basic" begin
        @testset "$(name)" for (name, nb) in nbs
            save_notebook(nb)
            # @info "File" name Text(read(nb.path,String))
            result = load_notebook_nobackup(nb.path)
            @test_notebook_inputs_equal(nb, result)
        end
    end

    @testset "I/O overloaded" begin
        @testset "$(name)" for (name, nb) in nbs
            let
                tasks = []
                for i in 1:16
                    push!(tasks, @async save_notebook(nb))
                    if i <= 8
                        sleep(0.01)
                    end
                end
                wait.(tasks)
                result = load_notebook_nobackup(nb.path)
                @test_notebook_inputs_equal(nb, result)
            end
        end
    end

    @testset "Bijection test" begin
        @testset "$(name)" for (name, nb) in nbs
            new_path = tempname()
            @assert !isfile(new_path)
            readwrite(nb.path, new_path)

            # load_notebook also does parsing and analysis - this is needed to save the notebook with cells in their correct order
            # load_notebook is how they are normally loaded, load_notebook_nobackup
            new_nb = load_notebook(new_path)

            before_contents = read(new_path, String)

            after_path = tempname()
            write(after_path, before_contents)

            after = load_notebook(after_path)
            after_contents = read(after_path, String)

            if name != String(nameof(bad_code_notebook))
                @test Text(before_contents) == Text(after_contents)
            end
        end
    end

    @testset "Recover from bad cell order" begin
        contents = """
        ### A Pluto.jl notebook ###
        # v0.17.3

        using Markdown
        using InteractiveUtils

        # ╔═╡ cdd40e28-61be-11ec-28fd-111111111111
        x = 1

        # ╔═╡ cdd40e28-61be-11ec-28fd-222222222222
        y = 2

        # ╔═╡ cdd40e28-61be-11ec-28fd-333333333333
        z = 3

        # ╔═╡ Cell order:
        # ╠═cdd40e28-61be-11ec-28fd-111111111111
        # ╠═cdd40e28-61be-11ec-28fd-333333333333
        # ╠═cdd40e28-61be-11ec-28fd-444444444444
        """

        path = tempname()
        write(path, contents)

        nb = load_notebook(path)

        @test nb.cell_order == UUID.([
            "cdd40e28-61be-11ec-28fd-111111111111",
            "cdd40e28-61be-11ec-28fd-333333333333",
            "cdd40e28-61be-11ec-28fd-222222222222",
        ])
        @test keys(nb.cells_dict) == Set(nb.cell_order)
    end

    # Some notebooks are designed to error (inside/outside Pluto)
    expect_error = [String(nameof(bad_code_notebook)), String(nameof(project_notebook)), "sample Interactivity.jl"]

    @testset "Runnable without Pluto" begin
        @testset "$(name)" for (name, nb) in nbs
            new_path = tempname()
            @assert !isfile(new_path)
            cp(nb.path, new_path)

            new_nb = load_notebook(new_path)

            # println(read(new_nb.path, String))

            if name ∉ expect_error
                @test jl_is_runnable(new_nb.path; only_undefvar=false)
            end
        end
    end

    @testset "Runnable with Pluto" begin
        @testset "$(name)" for (name, nb) in nbs
            if name ∉ expect_error
                @test nb_is_runnable(🍭, nb)
                WorkspaceManager.unmake_workspace((🍭, nb))
            end
        end
    end

    @testset "Line endings" begin
        @testset "$(name)" for (name, nb) in nbs
            file_contents = sprint(save_notebook, nb)

            let
                result = sread(load_notebook_nobackup, file_contents, nb.path)
                @test_notebook_inputs_equal(nb, result)
            end

            let
                file_contents_windowsed = replace(file_contents, "\n" => "\r\n")
                result_windowsed = sread(load_notebook_nobackup, file_contents_windowsed, nb.path)
                @test_notebook_inputs_equal(nb, result_windowsed)
            end
        end
    end

    @testset "Backups" begin
        @testset "$(name)" for (name, nb) in nbs
            save_notebook(nb)

            new_dir = mktempdir()
            new_path = joinpath(new_dir, "nb.jl")
            cp(nb.path, new_path)
            @test_nowarn load_notebook(new_path)
            @test num_backups_in(new_dir) == 0

            # Delete last line
            # cp(nb.path, new_path, force=true)
            # to_write = readlines(new_path)[1:end - 1]
            # write(new_path, join(to_write, '\n'))
            # @test_logs (:warn, r"Backup saved to") load_notebook(new_path)
            # @test num_backups_in(new_dir) == 1

            # Duplicate first line
            cp(nb.path, new_path, force=true)
            to_write = let
                old_lines = readlines(new_path)
                [old_lines[1], old_lines...]
            end
            write(new_path, join(to_write, '\n'))

            @test_nowarn load_notebook(new_path)
            @test num_backups_in(new_dir) == 0
            @test readdir(new_dir) == ["nb.jl"]

            # Extra stuff in preamble
            cp(nb.path, new_path, force=true)
            to_write = let
                old_content = read(new_path, String)
                replace(old_content, "using Markdown" => "using Markdown\n1 + 1")
            end
            write(new_path, to_write)

            @test_nowarn load_notebook(new_path)
            @test num_backups_in(new_dir) == 0
            @test readdir(new_dir) == ["nb.jl"]

            # Extra stuff at the end of the file
            cp(nb.path, new_path, force=true)
            to_write = let
                old_lines = readlines(new_path)
                [old_lines..., "", "1 + 1", Pluto._cell_id_delimiter * "heyyyy", "# coolio"]
            end

            write(new_path, join(to_write, '\n'))

            @test_logs (:warn, r"Backup saved to") load_notebook(new_path)
            @test num_backups_in(new_dir) == 1

            @test readdir(new_dir) == ["nb backup 1.jl", "nb.jl"]

            # AGAIN

            write(new_path, join(to_write, '\n'))

            @test_logs (:warn, r"Backup saved to") load_notebook(new_path)
            @test num_backups_in(new_dir) == 2

            @test Set(readdir(new_dir)) == Set(["nb backup 2.jl", "nb backup 1.jl", "nb.jl"] )
        end
    end

    @testset "Import & export HTML" begin
        nb = basic_notebook()
        export_html = Pluto.generate_html(nb)

        embedded_jl = Pluto.embedded_notebookfile(export_html)
        jl_path = tempname()
        write(jl_path, embedded_jl)

        result = load_notebook_nobackup(jl_path)
        @test_notebook_inputs_equal(nb, result, false)

        filename = "howdy.jl"

        export_html = Pluto.generate_html(nb; notebookfile_js=filename)
        @test occursin(filename, export_html)
        @test_throws ArgumentError Pluto.embedded_notebookfile(export_html)
    end

    @testset "Utilities" begin
        @testset "Cute file names" begin
            trash = mktempdir()
            for i in 1:200
                numbered_until_new(joinpath(trash, cutename()); create_file=true)
            end

            @test all(!isfile, [numbered_until_new(joinpath(trash, cutename()); create_file=false) for _ in 1:200])
        end
    end

    @test without_pluto_file_extension("julia编程.jl") == "julia编程"
    @test without_pluto_file_extension("julia编程.jl🐼") == "julia编程.jl🐼"
    @test without_pluto_file_extension("asdf.pluto.jl") == "asdf"
    @test without_pluto_file_extension("asdf.jl") == "asdf"

    # TODO: test bad dirs, filenames, permissions
end


text:
import Random
rng0 = Random.MersenneTwister(1234);

const Mat = Matrix{Float64}
const Vec = Vector{Float64}

struct RNN
    Wx :: Mat  # Transform input to hidden
    Wh :: Mat  # Transform hidden to hidden
    b  :: Vec  # Bias hidden to hidden
end

# RNN model: Given input sequence xs, compute output sequence
function rnn(weights :: RNN, xs :: Array{Vec}) :: Array{Vec}
    h = rand(rng, length(weights.b))
    [
      h = tanh.(weights.Wh * h) + tanh.(weights.Wx * x) + weights.b
         for x in xs
    ]
end

# RNN model: Given input sequence xs, compute output sequence
function rnn2(weights :: RNN, xs :: Array{Vec})
    h = rand(rng, length(weights.b))
    show(h)
    out = []
    for x in xs
        h = tanh.(weights.Wh * h) + tanh.(weights.Wx * x) + weights.b
        push!(out, h)
    end
    out
end

weights = RNN(rand(5,3), rand(5,5), rand(5))

xs = [rand(rng, 3) for i in 1:5];

rng = copy(rng0);
r1 = rnn(weights, xs)

rng = copy(rng0);
r2 = rnn2(weights, xs)

max_stars_repo_name: tomjaguarpaw/knossos-ksc

max_issues_repo_name: tomjaguarpaw/knossos-ksc

max_forks_repo_name: tomjaguarpaw/knossos-ksc


text:
# Deprecated functions

### Deprecated at version v0.2.14

@Base.deprecate sum_fdiff sumfdiff
@Base.deprecate max_fdiff maxfdiff
@Base.deprecate min_fdiff minfdiff
@Base.deprecate wsum_fdiff wsumfdiff

@Base.deprecate sum_fdiff! sumfdiff!
@Base.deprecate max_fdiff! maxfdiff!
@Base.deprecate min_fdiff! minfdiff!
@Base.deprecate wsum_fdiff! wsumfdiff!

# @Base.deprecate asum sumabs  # since asum exists in BLAS, we keep both asum & sumabs
@Base.deprecate amax maxabs
@Base.deprecate amin minabs
@Base.deprecate sqsum sumsq

@Base.deprecate asum! sumabs!
@Base.deprecate amax! maxabs!
@Base.deprecate amin! minabs!
@Base.deprecate sqsum! sumsq!

@Base.deprecate adiffsum sumabsdiff
@Base.deprecate adiffmax maxabsdiff
@Base.deprecate adiffmin minabsdiff
@Base.deprecate sqdiffsum sumsqdiff

@Base.deprecate adiffsum! sumabsdiff!
@Base.deprecate adiffmax! maxabsdiff!
@Base.deprecate adiffmin! minabsdiff!
@Base.deprecate sqdiffsum! sumsqdiff!

@Base.deprecate sum_xlogx sumxlogx
@Base.deprecate sum_xlogy sumxlogy

@Base.deprecate sum_xlogx! sumxlogx!
@Base.deprecate sum_xlogy! sumxlogy!

@Base.deprecate wasum wsumabs
@Base.deprecate wadiffsum wsumabsdiff
@Base.deprecate wsqsum wsumsq
@Base.deprecate wsqdiffsum wsumsqdiff

@Base.deprecate wasum! wsumabs!
@Base.deprecate wadiffsum! wsumabsdiff!
@Base.deprecate wsqsum! wsumsq!
@Base.deprecate wsqdiffsum! wsumsqdiff!

@Base.deprecate vdiffnorm vnormdiff
@Base.deprecate vdiffnorm! vnormdiff!

@Base.deprecate unsafe_view view


text: function f() :: Int64 "foo" end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using .CuArrays: CuArrays, @cufunc, CuMatOrAdj, CuOrAdj, CUBLAS

### Base.:\

import Base: \

# Tracker

using Tracker: Tracker, TrackedArray, @grad, data

Base.inv(x::CuMatOrAdj{<:AbstractFloat}) = CUBLAS.matinv_batched([x])[2][1]

A::CuMatOrAdj   \ B::TrackedArray = Tracker.track(\, A, B)
A::TrackedArray \ B::CuOrAdj      = Tracker.track(\, A, B)

@grad function (A::Union{CuMatOrAdj,TrackedArray} \ B::Union{CuOrAdj,TrackedArray})
    return data(A) \ data(B), function (Δ)
        AtransposedivΔ = transpose(A) \ Δ
        ∇A = transpose(A \ B * transpose(-AtransposedivΔ))
        return (∇A,  AtransposedivΔ)
    end
end

# Zygote

using Zygote: @adjoint

@adjoint Base.:\(A::CuMatOrAdj, B::CuOrAdj) = A \ B, function (Δ)
    AtransposedivΔ = transpose(A) \ Δ
    ∇A = transpose(A \ B * transpose(-AtransposedivΔ))
    return (∇A,  AtransposedivΔ)
end

### StatsFuns

import StatsFuns: logit, logistic, log1pexp, logexpm1

@cufunc logit(x) = log(x / (one(x) - x))
@cufunc logistic(x) = inv(exp(-x) + one(x))

# Helper

function seed!(s::Int)
    Random.seed!(s)
    CuArrays.CURAND.seed!(s)
end


text:
#methods from https://doi.org/10.1016/j.cam.2009.11.010
@muladd function perform_step!(integrator,cache::WangLi3SMil_AConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  #stage 1
  k  = integrator.f(uprev,p,t)
  u  = uprev + dt*k

  #stage 2
  k  = integrator.g(u,p,t)
  tmp = u + k*integrator.sqdt
  k₁ = integrator.g(tmp,p,t)
  k₁ = (k₁ - k)/(integrator.sqdt)
  u  = u - (dt/2)*k₁

  #stage 3
  k  = integrator.g(u,p,t)
  tmp = u + k*integrator.sqdt
  k₁ = integrator.g(tmp,p,t)
  k₁ = (k₁ - k)/(integrator.sqdt)
  u  = u .+ k .* W.dW .+ 0.5 .* (W.dW .^ 2) .* k₁
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_ACache,f=integrator.f)
  @unpack tmp,k,k₁ = cache
  @unpack t,dt,uprev,u,W,p = integrator

  #stage 1
  integrator.f(k,uprev,p,t)
  @.. u  = uprev + dt*k

  #stage 2
  integrator.g(k,u,p,t)
  @.. tmp = u + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁ = (k₁ - k)/(integrator.sqdt)
  @.. u  = u - (dt/2)*k₁

  #stage 3
  integrator.g(k,u,p,t)
  @.. tmp = u + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁ = (k₁ - k)/(integrator.sqdt)
  @.. u  = u + k*W.dW + (W.dW^2/2)*k₁
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_BConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  #stage 1
  k   = integrator.g(uprev,p,t)
  tmp = uprev + k*integrator.sqdt
  k₁  = integrator.g(tmp,p,t)
  k₁  = (k₁ - k)/integrator.sqdt
  u  = uprev .+ k .* W.dW .+ 0.5 .* (W.dW .^ 2) .* k₁

  #stage 2
  k  = integrator.f(u,p,t)
  u  = u + dt*k

  #stage 3
  k   = integrator.g(u,p,t)
  tmp = u + k*integrator.sqdt
  k₁  = integrator.g(tmp,p,t)
  k₁  = (k₁ - k)/integrator.sqdt
  u   = u - (dt/2)*k₁
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_BCache,f=integrator.f)
  @unpack tmp,k,k₁ = cache
  @unpack t,dt,uprev,u,W,p = integrator

    #stage 1
  integrator.g(k,uprev,p,t)
  @.. tmp = uprev + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁  = (k₁ - k)/integrator.sqdt
  @.. u   = uprev + W.dW*k + (W.dW^2/2)*k₁

  #stage 2
  integrator.f(k,u,p,t)
  @.. u  = u + dt*k

  #stage 3
  integrator.g(k,u,p,t)
  @.. tmp = u + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁  = (k₁ - k)/integrator.sqdt
  @.. u   = u - (dt/2)*k₁

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_CConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  #stage 1
  k  = integrator.f(uprev,p,t)
  u  = uprev + dt*k

  #stage 2
  k  = integrator.g(u,p,t)
  tmp = u + k*integrator.sqdt
  k₁ = integrator.g(tmp,p,t)
  k₁ = (k₁ - k)/(integrator.sqdt)
  u  = u .+ k .* W.dW .+ 0.5 .* (W.dW .^ 2) .* k₁

  #stage 3
  k  = integrator.g(u,p,t)
  tmp = u + k*integrator.sqdt
  k₁ = integrator.g(tmp,p,t)
  k₁ = (k₁ - k)/(integrator.sqdt)
  u  = u - (dt/2)*k₁
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_CCache,f=integrator.f)
  @unpack tmp,k,k₁ = cache
  @unpack t,dt,uprev,u,W,p = integrator

  #stage 1
  integrator.f(k,uprev,p,t)
  @.. u  = uprev + dt*k

  #stage 2
  integrator.g(k,u,p,t)
  @.. tmp = u + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁ = (k₁ - k)/(integrator.sqdt)
  @.. u  = u + k*W.dW + (W.dW^2/2)*k₁

  #stage 3
  integrator.g(k,u,p,t)
  @.. tmp = u + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁ = (k₁ - k)/(integrator.sqdt)
  @.. u  = u - (dt/2)*k₁
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_DConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  #stage 1
  k  = integrator.g(uprev,p,t)
  tmp = uprev + k*integrator.sqdt
  k₁ = integrator.g(tmp,p,t)
  k₁ = (k₁ - k)/(integrator.sqdt)
  u  = uprev - (dt/2)*k₁

  #stage 2
  k  = integrator.f(u,p,t)
  u  = u + dt*k

  #stage 3
  k  = integrator.g(u,p,t)
  tmp = u + k*integrator.sqdt
  k₁ = integrator.g(tmp,p,t)
  k₁ = (k₁ - k)/(integrator.sqdt)
  u  = u .+ k .* W.dW .+ 0.5 .* (W.dW .^ 2) .* k₁
 integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_DCache,f=integrator.f)
  @unpack tmp,k,k₁ = cache
  @unpack t,dt,uprev,u,W,p = integrator

  #stage 1
  integrator.g(k,uprev,p,t)
  @.. tmp = uprev + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁ = (k₁ - k)/(integrator.sqdt)
  @.. u  = uprev - (dt/2)*k₁

  #stage 2
  integrator.f(k,u,p,t)
  @.. u  = u + dt*k

  #stage 3
  integrator.g(k,u,p,t)
  @.. tmp = u + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁ = (k₁ - k)/(integrator.sqdt)
  @.. u  = u + k*W.dW + (W.dW^2/2)*k₁
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_EConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  #stage 1
  k  = integrator.g(uprev,p,t)
  tmp = uprev + k*integrator.sqdt
  k₁ = integrator.g(tmp,p,t)
  k₁ = (k₁ - k)/(integrator.sqdt)
  u  = uprev - (dt/2)*k₁

  #stage 2
  k  = integrator.g(u,p,t)
  tmp = u + k*integrator.sqdt
  k₁ = integrator.g(tmp,p,t)
  k₁ = (k₁ - k)/(integrator.sqdt)
  u  = u .+ k .* W.dW .+ 0.5 .* (W.dW .^ 2) .* k₁

  #stage 3
  k  = integrator.f(u,p,t)
  u  = u + dt*k
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_ECache,f=integrator.f)
  @unpack tmp,k,k₁ = cache
  @unpack t,dt,uprev,u,W,p = integrator

  #stage 1
  integrator.g(k,uprev,p,t)
  @.. tmp = uprev + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁ = (k₁ - k)/(integrator.sqdt)
  @.. u  = uprev - (dt/2)*k₁

  #stage 2
  integrator.g(k,u,p,t)
  @.. tmp = u + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁ = (k₁ - k)/(integrator.sqdt)
  @.. u  = u + k*W.dW + (W.dW^2/2)*k₁

  #stage 3
  integrator.f(k,u,p,t)
  @.. u  = u + dt*k
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_FConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  #stage 1
  k   = integrator.g(uprev,p,t)
  tmp = uprev + k*integrator.sqdt
  k₁  = integrator.g(tmp,p,t)
  k₁  = (k₁ - k)/integrator.sqdt
  u  = uprev .+ k .* W.dW .+ 0.5 .* (W.dW .^ 2) .* k₁

  #stage 2
  k   = integrator.g(u,p,t)
  tmp = u + k*integrator.sqdt
  k₁  = integrator.g(tmp,p,t)
  k₁  = (k₁ - k)/integrator.sqdt
  u   = u - (dt/2)*k₁

  #stage 3
  k  = integrator.f(u,p,t)
  u  = u + dt*k
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::WangLi3SMil_FCache,f=integrator.f)
  @unpack tmp,k,k₁ = cache
  @unpack t,dt,uprev,u,W,p = integrator

    #stage 1
  integrator.g(k,uprev,p,t)
  @.. tmp = uprev + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁  = (k₁ - k)/integrator.sqdt
  @.. u   = uprev + W.dW*k + (W.dW^2/2)*k₁

  #stage 2
  integrator.g(k,u,p,t)
  @.. tmp = u + k*integrator.sqdt
  integrator.g(k₁,tmp,p,t)
  @.. k₁  = (k₁ - k)/integrator.sqdt
  @.. u   = u - (dt/2)*k₁

  #stage 3
  integrator.f(k,u,p,t)
  @.. u  = u + dt*k
  integrator.u = u
end


text:
cfg = get_config("test.yaml")
cfg["innovation_max"] = 5
cfg["p_mut_weights"] = 1.0

ind = NEATInd(cfg)
clone = deepcopy(ind)

@testset "Mutate" begin
    mut_w = NEAT.mutate_weights(ind, cfg)
    @test length(mut_w.connections) == length(ind.connections)
    for ci in eachindex(mut_w.connections)
        @test mut_w.connections[ci].in_node == ind.connections[ci].in_node
        @test mut_w.connections[ci].out_node == ind.connections[ci].out_node
        @test mut_w.connections[ci].weight != ind.connections[ci].weight
    end
    # test that mutation didn't effect ind
    test_identical(ind, clone)

    # add a neuron
    mut_add_neur = NEAT.mutate_add_neuron(ind, cfg)
    @test length(mut_add_neur.connections) == length(ind.connections) + 1
    @test length(mut_add_neur.neurons) == length(ind.neurons) + 1
    @test length(clone.connections) == length(ind.connections)
    @test length(clone.neurons) == length(ind.neurons)
    test_identical(ind, clone)

    # add a connection to original network, no change expected
    mut_add_con = NEAT.mutate_add_connection(ind, cfg)
    @test length(mut_add_con.connections) == length(mut_add_neur.connections)
    test_identical(ind, mut_add_con)

    # add connection to network with additional neuron, expect new connection
    clone2 = deepcopy(mut_add_neur)
    mut_add_con = NEAT.mutate_add_connection(mut_add_neur, cfg)
    @test length(mut_add_con.connections) == length(mut_add_neur.connections) + 1
    test_identical(mut_add_neur, clone2)

    # toggle enabled on n_times connections (not really used in practice)
    mut_en = NEAT.mutate_enabled(ind; n_times=length(ind.connections))
    @test length(mut_en.connections) == length(ind.connections)
    for ci in eachindex(mut_en.connections)
        @test mut_en.connections[ci].enabled != ind.connections[ci].enabled
    end
    test_identical(ind, clone)

    # mutate
    mut = mutate(ind, cfg)
    @test mut.n_in == ind.n_in
    @test mut.n_out == ind.n_out
    test_identical(ind, clone)
end


text:
# ---
# title: 202. Happy Number
# id: problem202
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/happy-number/description/>
# hidden: true
# ---
#
# Write an algorithm to determine if a number `n` is "happy".
#
# A happy number is a number defined by the following process: Starting with any
# positive integer, replace the number by the sum of the squares of its digits,
# and repeat the process until the number equals 1 (where it will stay), or it
# **loops endlessly in a cycle** which does not include 1. Those numbers for
# which this process **ends in 1** are happy numbers.
#
# Return True if `n` is a happy number, and False if not.
#
# **Example:  **
#
#
#
#     Input: 19
#     Output: true
#     Explanation: 12 + 92 = 82
#     82 + 22 = 68
#     62 + 82 = 100
#     12 + 02 + 02 = 1
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Lexer
using Error

export lex

function lex(source)
  if source[1] != '('
    try
      ret =  parse(source)
      if source[1] == '['
        if source[end] == ']'
          return eval(ret)
        else
          throw(LispError("Incomplete array Epression"))
        end
      end

      return ret
    catch
      throw(LispError("Incomplete Expression: $source"))
    end
	end
  if source[end] != ')'
		throw(LispError("Incomplete expression: $source"))
	end
  source = source[2:end - 1]
	expr = split_exps(source)
	out = Any[]
  for x in expr
    i = lex(x)
    push!(out, i)
  end
  return out
end

##
## Below are a few useful utility functions. These should come in handy when
## implementing `parse`. We don't want to spend the day implementing parenthesis
## counting, after all.
##

# Remove from a string anything in between a ; and a linebreak
function remove_comments(source)
  replace(source, r";.*\n", "\n")
end

## Given a string and the index of an opening parenthesis, determines
## the index of the matching closing paren.
function find_matching_paren(source, char_start, char_end, start=1, )
  @assert source[start] == char_start
  pos = start
  open_brackets = 1
  while open_brackets > 0
    pos += 1
    #if length(source) == pos
    if length(source) < pos
      throw(LispError("Incomplete expression: $(source[start:end])"))
    end
    if source[pos] == char_start
      open_brackets += 1
    end
    if source[pos] == char_end
      open_brackets -= 1
    end
  end
  return pos
end

## Splits a source string into subexpressions
## that can be parsed individually.

## Example:

##     > split_exps("foo bar (baz 123)")
##     ["foo", "bar", "(baz 123)"]
function split_exps(source)
  rest = strip(source)
  exps = AbstractString[]
  while !isempty(rest)
    exp, rest = first_expression(rest)
    push!(exps, exp)
  end
  return exps
end

## Split string into (exp, rest) where exp is the
## first expression in the string and rest is the
## rest of the string after this expression.
function first_expression(source)
  source = strip(source)
  if source[1] == '('
    last = find_matching_paren(source, '(', ')')
    return source[1:last], source[last + 1:end]
  elseif source[1] == '['
    last = find_matching_paren(source, '[', ']')
    return source[1:last], source[last + 1:end]
  else
    m = match(r"^[^\s)'\]]+", source)
    if typeof(m) == Void
      throw(ParseError("Syntax Error"))
    end
    _end = m.offset + length(m.match) - 1
    atom = source[1:_end]
    return atom, source[_end+1:end]
  end
end
end # module

max_stars_repo_name: scottcrossen/JuliaDevelopment

max_issues_repo_name: scottcrossen/JuliaDevelopment

max_forks_repo_name: scottcrossen/JuliaDevelopment


text:
@testset "smart_authorization.jl" begin
    smart_result = Foo()
    @test !has_fhir_patient_id(smart_result)
    @test !has_fhir_encounter_id(smart_result)
end


text:
# Utility file I/O functions. Not fit for general purpose use.
function load_dipha(filename)
    open(filename, "r") do f
        magic = read(f, Int)
        @assert magic == 8067171840
        type = read(f, Int)
        @assert type == 1
        len = read(f, Int)
        dim = read(f, Int)
        if dim == 2
            m = read(f, Int)
            n = read(f, Int)
            result = Array{Float64}(undef, (m, n))
            read!(f, result)

            return result
        elseif dim == 3
            m = read(f, Int)
            n = read(f, Int)
            o = read(f, Int)
            result = Array{Float64}(undef, (m, n, o))
            read!(f, result)

            return result
        else
            error("not implemented for dim=$dim")
        end
    end
end

function save_dipha(filename, data)
    open(filename, "w") do f
        write(f, 8067171840)
        write(f, 1)
        write(f, length(data))
        write(f, length(size(data)))
        for s in size(data)
            write(f, s)
        end
        for i in data
            write(f, Float64(i))
        end
    end
end

function load_points(filename)
    # ripser uses Float32
    table = CSV.read(filename, header=0, type=Float32, delim=' ')
    nrow, dim = size(table)
    result = Vector{NTuple{dim, Float32}}(undef, nrow)
    for i in 1:nrow
        result[i] = Tuple(table[i, :])
    end
    return result
end

function load_dist(filename)
    return Matrix(CSV.read(filename, header=0, type=Float32, delim=' '))
end

function save_sparse(filename, data)
    I, J, V = findnz(sparse(LowerTriangular(data)))
    return CSV.write(joinpath(filename), (I=I.-1, J=J.-1, V=V), writeheader=false, delim=' ')
end

function load_sparse(filename)
    table = CSV.read(filename, header=0, type=Float32, delim=' ')
    I = table[:, 1] .+ 1
    J = table[:, 2] .+ 1
    V = table[:, 3] .+ 1
    sparse([I; J], [J; I], [V; V])
end

function load_data(filename)
    ext = splitext(filename)[2]
    if ext == ".pts" || ext == ".alpha"
        return load_points(filename)
    elseif ext == ".dist"
        return load_dist(filename)
    elseif ext == ".spdist"
        return load_sparse(filename)
    elseif ext == ".dipha"
        return load_dipha(filename)
    else
        error("unsupported filetype \"$(ext)\"")
    end
end


text:
# SoftMax Arm selection algorithm

"""
    SoftMax Arm Selection
"""
struct SoftMax <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    τ::Float64              # temparature parameter
    play_count::Vector{Int64}    # Number of times each arm is played
    avgValue::Vector{Float64}   # Average value of each arm

    function SoftMax( noOfArms::Integer, τ::AbstractFloat )
        new( noOfArms,
             0,
             0,
             τ,
             zeros(Int64,noOfArms),
             zeros(Float64,noOfArms)
        )
    end
end

function get_arm_index( agent::SoftMax )
    pVec    = exp.(agent.avgValue/agent.τ) / sum(exp.(agent.avgValue/agent.τ))
    agent.lastPlayedArm = rand( Categorical(pVec) )
    return agent.lastPlayedArm
end

function update_reward!( agent::SoftMax, r::Real )
    agent.avgValue[agent.lastPlayedArm] = (agent.play_count[agent.lastPlayedArm]*agent.avgValue[agent.lastPlayedArm] + r) / (agent.play_count[agent.lastPlayedArm]+1)
    agent.play_count[agent.lastPlayedArm] += 1
    agent.noOfSteps += 1

    nothing
end

function reset!( agent::SoftMax )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0

    agent.play_count    = zeros( Int64, agent.noOfArms )
    agent.avgValue      = zeros( Float64, agent.noOfArms )

    nothing
end

function info_str( agent::SoftMax, latex::Bool )
    if latex
        return @sprintf( "SoftMax (\$\\tau = %4.3f\$)", agent.τ )
    else
        return @sprintf( "SoftMax (τ = %4.3f)", agent.τ )
    end
end


text:
using POMDPModels
using Test
using POMDPTesting

pomdp = MiniHallway()

@test has_consistent_distributions(pomdp)

@test simulate(RolloutSimulator(max_steps=100), pomdp, RandomPolicy(pomdp)) >= 0.0

@test stateindex(pomdp, 5) == 5
@test actionindex(pomdp, 1) == 1
@test obsindex(pomdp, 1) == 1

@test reward(pomdp, 11, 1) == 1