text:
struct BIOM
    data::SparseMatrixCSC
    sample::NamedTuple
    observation::NamedTuple
end

function BIOM(data::SparseMatrixCSC,
        sample::Union{Nothing, NamedTuple}=nothing,
        observation::Union{Nothing, NamedTuple}=nothing)::BIOM
    if isnothing(sample)
        sample = (ids=Array([]), metadata=Dict(), groupmetadata=Dict())
    end
    if isnothing(observation)
        observation = (ids=Array([]), metadata=Dict(), groupmetadata=Dict())
    end

    required_keys = (:ids, :metadata, :groupmetadata)
    if  keys(sample) != required_keys
        throw(ArgumentError("sample should have keys: $(required_keys), but has keys: $(keys(sample))"))
    end
    if  keys(observation) != required_keys
        throw(ArgumentError("observation should have keys: $(required_keys), but has keys: $(keys(observation))"))
    end

    BIOM(data, sample, observation)
end

function BIOM(data::Matrix,
        sample::Union{Nothing, NamedTuple}=nothing,
        observation::Union{Nothing, NamedTuple}=nothing)::BIOM
    BIOM(sparse(data), sample, observation)
end

function BIOM(biom::Dict)
    #= parse data =#
    ptrs = biom["sample"]["matrix"]["indptr"]
    s = cat([fill(i-1, ptrs[i] - ptrs[i-1]) for i in 2:length(ptrs)]..., dims=1)
    o = biom["sample"]["matrix"]["indices"] .+ 1
    d = biom["sample"]["matrix"]["data"]
    data = sparse(s, o, d)
    #= read metadata =#
    si = biom["sample"]["ids"]
    sm = biom["sample"]["metadata"]
    sg = biom["sample"]["group-metadata"]
    oi = biom["observation"]["ids"]
    om = biom["observation"]["metadata"]
    og = biom["observation"]["group-metadata"]
    sample      = (ids=si, metadata=sm, groupmetadata=sg)
    observation = (ids=oi, metadata=om, groupmetadata=og)
    #= struct =#
    BIOM(data, sample, observation)
end

function BIOM(s::String)::BIOM
    biomdict = h5open(s, "r") do h5
        read(h5)
    end
    BIOM(biomdict)
end

Base.size(S::BIOM, i...) = Base.size(S.data, i...)

"""
    getindex(biom::BIOM, sample_ind, observation_ind) -> BIOM

return a subset of the data, and metadata in `BIOM` as specified by `sample_ind` and `observation_ind`.

Requires exactly two indices, but otherwise behaves the same as `getindex()`
"""
function Base.getindex(S::BIOM, i, j)
	function dict_subset(dic, sub)
		K =  keys(dic)
		subdict = Dict()
		for k in K
            # if size(dic(k)) == 0 || dic(k) === nothing
            #     subdict[k] = dic[k]
            if ndims(dic[k]) == 1
				subdict[k] = dic[k][sub]
			elseif ndims(dic[k]) == 2
				subdict[k] = dic[k][:, sub]
			else
				throw(error("$(ndims(dic[k])) dims in dict not implemented"))
			end
		end
		return subdict
	end
    data = SparseMatrixCSC(Base.getindex(S.data, i, j))
	sample = (
		ids = S.sample.ids[i],
		metadata = dict_subset(S.sample.metadata, i),
		groupmetadata = dict_subset(S.sample.groupmetadata, i)
	)
	observation = (
		ids = S.observation.ids[j],
		metadata = dict_subset(S.observation.metadata, j),
		groupmetadata = dict_subset(S.observation.groupmetadata, j)
	)

	return BIOM(data, sample, observation)
end

function Base.hash(S::BIOM)
    d = hash(S.data)
    s = hash(S.sample, d)
    o = hash(S.observation, s)
    return o
end

function Base.isequal(S::BIOM, T::BIOM)
    hash(S) == hash(T)
end

function (==)(S::BIOM, T::BIOM)
    Base.isequal(S, T)
end

Base.setindex!(S::BIOM, v, i, j) = Base.setindex!(S.data, v, i, j)


text:
module PendulumEnv
# Ported from: https://github.com/openai/gym/blob/996e5115621bf57b34b9b79941e629a36a709ea1/gym/envs/classic_control/pendulum.py
#              https://github.com/openai/gym/wiki/Pendulum-v0

using Reinforce: AbstractEnvironment
using LearnBase: IntervalSet
using RecipesBase
using Distributions
using Random: seed!

import Reinforce: reset!, actions, finished, step!, state

export
  Pendulum,
  reset!,
  step!,
  actions,
  finished,
  state

const max_speed = 8.0
const max_torque = 2.0

angle_normalize(x) = ((x+π) % (2π)) - π

mutable struct PendulumState{T<:AbstractFloat} <: AbstractVector{T}
  θ::T
  θvel::T
end

PendulumState() = PendulumState(0., 0.)

Base.size(::PendulumState) = (2,)

function Base.getindex(s::PendulumState, i::Int)
  (i > length(s)) && throw(BoundsError(s, i))
  ifelse(i == 1, s.θ, s.θvel)
end

function Base.setindex!(s::PendulumState, x, i::Int)
  (i > length(s)) && throw(BoundsError(s, i))
  setproperty!(s, ifelse(i == 1, :θ, :θvel), x)
end

mutable struct Pendulum{V<:AbstractVector} <: AbstractEnvironment
  state::V
  reward::Float64
  a::Float64 # last action for rendering
  steps::Int
  maxsteps::Int
end

Pendulum(maxsteps = 500) = Pendulum(PendulumState(),0., 0., 0, maxsteps)

function reset!(env::Pendulum)
  env.state = PendulumState(rand(Uniform(-π, π)), rand(Uniform(-1., 1.)))
  env.reward = 0.0
  env.a = 0.0
  env.steps = 0
  env
end

actions(env::Pendulum, s) = IntervalSet(-max_torque, max_torque)

function step!(env::Pendulum, s, a)
  θ, θvel = env.state
  g = 10.0
  m = 1.0
  l = 1.0
  dt = 0.05

  env.a = a
  a = clamp(a, -max_torque, max_torque)
  env.reward = -(angle_normalize(θ)^2 + 0.1θvel^2 + 0.001a^2)

  # update state
  newθvel = θvel + (-1.5g/l * sin(θ+π) + 3/(m*l^2)*a) * dt
  newθ = θ + newθvel * dt
  newθvel = clamp(newθvel, -max_speed, max_speed)
  env.state = PendulumState(newθ, newθvel)

  env.steps += 1
  env.reward, env.state
end

function state(env::Pendulum)
  θ, θvel = env.state
  Float64[cos(θ), sin(θ), θvel]
end

finished(env::Pendulum, s′) = env.steps >= env.maxsteps

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

@recipe function f(env::Pendulum)
  legend := false
  xlims := (-1,1)
  ylims := (-1,1)
  grid := false
  ticks := nothing

	# pole
  @series begin
    w = 0.2
    x = [-w,w,w,-w]
    y = [-.1,-.1,1,1]
    θ = env.state.θ
    fillcolor := :red
    seriestype := :shape
    x*cos(θ) - y*sin(θ), y*cos(θ) + x*sin(θ)
  end

  # center
  @series begin
    seriestype := :scatter
    markersize := 10
    markercolor := :black
    annotations := [(0, -0.2, "a: $(round(env.a, digits = 4))", :top)]
    [0],[0]
  end
end

end

ext: jl


text:
@testset "PorousModels" begin

    @testset "BrooksAndCorey" begin
        # HF135 membrane, Brooks and Corey model
        # Data from Buser (PhD thesis, 2016)
        # http://hdl.handle.net/1773/38064
        θr = 1.07e-6
        θs = 0.936
        k = 5.50e-13  # m**2
        α = 1/1.96  # 1/m
        n = 0.669

        bc = BrooksAndCorey(n=n, α=α, k=k, θr=θr, θs=θs)

        htest = -10.0

        @test @inferred hθ(bc, @inferred θh(bc, htest)) ≈ htest
        @test @inferred Ch(bc, htest) ≈ derivative(h -> θh(bc, h), htest)
        @test @inferred Cθ(bc, @inferred θh(bc, htest)) ≈ @inferred Ch(bc, htest)
        @test @inferred Kθ(bc, @inferred θh(bc, htest)) ≈ @inferred Kh(bc, htest)

        θtest = 0.5

        @test Dθ(bc, θtest) ≈ 1.8877259646141242e-06

        @test @inferred θh(bc, @inferred hθ(bc, θtest)) ≈ θtest
        @test @inferred Ch(bc, @inferred hθ(bc, θtest)) ≈ @inferred Cθ(bc, θtest)
        @test @inferred Kh(bc, @inferred hθ(bc, θtest)) ≈ @inferred Kθ(bc, θtest)
        @test @inferred Dθ(bc, θtest) ≈ @inferred Kθ(bc, θtest)/@inferred Cθ(bc, θtest)
    end

    @testset "VanGenuchten" begin
        # HF135 membrane, Van Genuchten model
        # Data from Buser (PhD thesis, 2016)
        # http://hdl.handle.net/1773/38064
        θr = 0.0473
        θs = 0.945
        k = 5.50e-13  # m**2
        α = 0.2555  # 1/m
        n = 2.3521

        vg = VanGenuchten(n=n, α=α, k=k, θr=θr, θs=θs)

        htest = -10.0

        @test @inferred Ch(vg, htest) ≈ 0.02896477570729651
        @test Kh(vg, htest) ≈ 9.430485870291618e-9

        @test @inferred hθ(vg, @inferred θh(vg, htest)) ≈ htest
        @test @inferred Ch(vg, htest) ≈ derivative(h -> θh(vg, h), htest)
        @test @inferred Cθ(vg, @inferred θh(vg, htest)) ≈ @inferred Ch(vg, htest)
        @test @inferred Kθ(vg, @inferred θh(vg, htest)) ≈ @inferred Kh(vg, htest)

        θtest = 0.5

        @test Dθ(vg, θtest) ≈ 1.769723354269708e-6

        @test @inferred θh(vg, @inferred hθ(vg, θtest)) ≈ θtest
        @test @inferred Ch(vg, @inferred hθ(vg, θtest)) ≈ @inferred Cθ(vg, θtest)
        @test @inferred Kh(vg, @inferred hθ(vg, θtest)) ≈ @inferred Kθ(vg, θtest)
        @test @inferred Dθ(vg, θtest) ≈ @inferred Kθ(vg, θtest)/@inferred Cθ(vg, θtest)
    end
end


text:
"""
    struct HFcBuffer{
            T,D,TB,Tβ,Tσ,Ta,Tmat,Tvec
            } <: DiffusionDefinition.AbstractBuffer{T}
        data::Vector{T}
        B::TB
        β::Tβ
        σ::Tσ
        a::Ta
        mat::Tmat
        vec::Tvec
    end

A buffer for temporary computations of in-place ODE solvers solving for H,F,c
system.
"""
struct HFcBuffer{
        TB,Tβ,Tσ,Ta,Tmat,Tvec
        } <: DiffusionDefinition.AbstractBuffer
    B::TB
    β::Tβ
    σ::Tσ
    a::Ta
    mat::Tmat
    vec::Tvec

    function HFcBuffer(
            B::TB, β::Tβ, σ::Tσ, a::Ta, mat::Tmat, vec::Tvec
        ) where {TB,Tβ,Tσ,Ta,Tmat,Tvec}
        new{TB,Tβ,Tσ,Ta,Tmat,Tvec}(B, β, σ, a, mat, vec)
    end
end

function HFcBuffer{T}(D::Integer) where T
    B = zeros(T, (D, D))
    β = zeros(T, D)
    σ = zeros(T, (D, D))
    a = zeros(T, (D, D))
    m = zeros(T, (D, D))
    v = zeros(T, D)

    HFcBuffer(B, β, σ, a, m, v)
end

#=

"""
    struct HFcContainer{T,D,TH,TF,Tc} <: AbstractBuffer{T}
        data::Vector{T}
        H::TH
        F::TF
        c::Tc
    end

A buffer containing data for in-place computations of H,F,c terms.
"""
struct HFcContainer{T,D,TH,TF,Tc} <: DiffusionDefinition.AbstractBuffer{T}
    data::Vector{T}
    H::TH
    F::TF
    c::Tc

    function HFcContainer{T}(D::Integer) where T
        data = zeros(T, D*D+D+1)
        H = reshape( view(data, 1:(D*D)), (D,D) )
        F = view(data, (D*D+1):(D*D+D))
        c = view(data, (D*D+D+1):(D*D+D+1))
        new{T,D,typeof(H),typeof(F),typeof(c)}(data, H, F, c)
    end

    function HFcContainer(data::AbstractArray{T,1}, D) where T
        H = reshape( view(data, 1:(D*D)), (D,D) )
        F = view(data, (D*D+1):(D*D+D))
        c = view(data, (D*D+D+1):(D*D+D+1))
        new{T,D,typeof(H),typeof(F),typeof(c)}(data, H, F, c)
    end
end

function Base.similar(c::HFcContainer{T,D}) where {T,D}
    HFcContainer(similar(c.data), D)
end

function Base.similar(c::HFcContainer{T,D}, ::Type{K}) where {T,D,K}
    HFcContainer(similar(c.data, K), D)
end

"""
    struct HFcBuffer{
            T,D,TB,Tβ,Tσ,Ta,Tmat,Tvec
            } <: DiffusionDefinition.AbstractBuffer{T}
        data::Vector{T}
        B::TB
        β::Tβ
        σ::Tσ
        a::Ta
        mat::Tmat
        vec::Tvec
    end

A buffer for temporary computations of in-place ODE solvers solving for H,F,c
system.
"""
struct HFcBuffer{
        T,D,TB,Tβ,Tσ,Ta,Tmat,Tvec
        } <: DiffusionDefinition.AbstractBuffer{T}
    data::Vector{T}
    B::TB
    β::Tβ
    σ::Tσ
    a::Ta
    mat::Tmat
    vec::Tvec

    function HFcBuffer{T}(D::Integer) where T
        data = zeros(T, 4*D*D+D)

        B = reshape( view(data, 1:(D*D)), (D,D) )
        β = view(data, (D*D+1):(D*D+D))
        σ = reshape( view(data, (D*D+D+1):(2*D*D+D)), (D,D) )
        a = reshape( view(data, (2*D*D+D+1):(3*D*D+D)), (D,D) )
        m = reshape( view(data, (3*D*D+D+1):(4*D*D+D)), (D,D) )
        v =  view(data, (3*D*D+D+1):(3*D*D+2*D))
        new{T,D,typeof(B),typeof(β),typeof(σ),typeof(a),typeof(m),typeof(v)}(
            data, B, β, σ, a, m, v
        )
    end
end

=#

#=
"""
    struct HFcContainer{T,D,TH,TF,Tc} <: AbstractBuffer{T}
        data::Vector{T}
        H::TH
        F::TF
        c::Tc
    end

A buffer containing data for in-place computations of H,F,c terms.
"""
struct HFcContainer{T,TD,TH,TF,Tc} <: DiffusionDefinition.AbstractBuffer{T}
    data::TD
    H::TH
    F::TF
    c::Tc

    function HFcContainer{T}(
            data::TD, H::TH, F::TF, c::Tc
        ) where {T,TD,TH,TF,Tc}
        new{T,TD,TH,TF,Tc}(data, H, F, c)
    end
end

function HFcContainer{T}(D::Integer) where T
    H = zeros(T, (D,D))
    F = zeros(T, D)
    c = zeros(T, 1)
    data = ArrayPartition(H, F, c)
    HFcContainer{T}(data, H, F, c)
end

function HFcContainer(data::ArrayPartition{T}) where T
    HFcContainer{T}(data, data.x[1], data.x[2], data.x[3])
end

function Base.similar(d::HFcContainer{T}) where T
    new_data = similar(d.data)
    HFcContainer{T}(new_data, new_data.x[1], new_data.x[2], new_data.x[3])
end
# the one below is likely wrong, but it's also not used atm, come back later
function Base.similar(d::HFcContainer{K}, ::Type{T}) where {K,T}
    new_data = similar(d.data, T)
    HFcContainer{T}(new_data, new_data.x[1], new_data.x[2], new_data.x[3])
end

#function Base.similar(c::HFcContainer{T,D}) where {T,D}
#    HFcContainer(similar(c.data), D)
#end

#function Base.similar(c::HFcContainer{T,D}, ::Type{K}) where {T,D,K}
#    HFcContainer(similar(c.data, K), D)
#end
=#


text:
include("load.jl")

using SSE
using BenchmarkTools

const N = 27
const lt = Chain(N)
const spins = rand_pispins(N, 200)
const model = TFIModel(1.0, 0.5, 0.01, lt)
const sseiter = PathIntegralIter(model, spins)
println("neighbor(lt, 3)")
display(@benchmark neighbors($lt, 3))
println("\nneighbor(lt[3])")
display(@benchmark neighbors($lt[3]))
println("\nspins[neighbors(lt[3])] .= 1")
display(@benchmark @inbounds $spins[neighbors(lt[3])] .= 1)
println("\nspins[neighbors(lt, 3)] .= 1")
display(@benchmark @inbounds $spins[neighbors(lt, 3)] .= 1)
println("\niterate(sseiter)")
display(@benchmark iterate($sseiter))


text:
@testset "Looks like in CoreModel, detailed test" begin
    cp = test_LP()
    @test isempty(filter(looks_like_exchange_reaction, reactions(cp)))

    cp = test_simpleLP()
    @test isempty(filter(looks_like_exchange_reaction, reactions(cp)))

    cp = CoreModel(
        [-1.0 -1 -2; 0 -1 0; 0 0 0],
        zeros(3),
        ones(3),
        ones(3),
        ones(3),
        ["EX_m1"; "r2"; "r3"],
        ["m1"; "m2"; "m3"],
    )
    @test filter(looks_like_exchange_reaction, reactions(cp)) == ["EX_m1"]

    cp = CoreModel(
        [-1.0 0 0; 0 0 -1; 0 -1 0],
        zeros(3),
        ones(3),
        ones(3),
        ones(3),
        ["EX_m1"; "Exch_m3"; "Ex_m2"],
        ["m1"; "m2"; "m3_e"],
    )
    @test filter(looks_like_exchange_reaction, reactions(cp)) == ["EX_m1", "Exch_m3", "Ex_m2"]
    @test filter(x -> looks_like_exchange_reaction(x; exchange_prefixes = ["Exch_"]), reactions(cp)) == ["Exch_m3"]

    # this is originally the "toyModel1.mat"
    cp = test_toyModel()

    @test filter(looks_like_exchange_reaction, reactions(cp)) == ["EX_m1(e)", "EX_m3(e)", "EX_biomass(e)"]
    @test filter(x -> looks_like_exchange_reaction(x; exclude_biomass=true), reactions(cp)) == ["EX_m1(e)", "EX_m3(e)"]
    @test filter(looks_like_exchange_metabolite, metabolites(cp)) == ["m1[e]", "m3[e]"]
    @test filter(looks_like_biomass_reaction, reactions(cp)) == ["EX_biomass(e)", "biomass1"]
    @test filter(x -> looks_like_biomass_reaction(x; exclude_exchanges=true), reactions(cp)) == ["biomass1"]
end

@testset "Looks like functions, basic" begin
    model = load_model(model_paths["e_coli_core.json"])
    @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
    @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
    @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1

    model = load_model(model_paths["e_coli_core.xml"])
    @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
    @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
    @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1

    model = load_model(model_paths["e_coli_core.mat"])
    @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
    @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
    @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1

    model = convert(StandardModel, model)
    @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
    @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
    @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1

    model = convert(CoreModelCoupled, model)
    @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
    @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
    @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1

end


text:
# open Demo/ViewController/ViewController.xcworkspace
# Run

using Swifter

vc = initial("http://localhost:8080")

@query vc.view.backgroundColor = UIColor.yellowColor()

sleep(0.5)

@query vc.label.text = "Hello Swift"

sleep(0.5)

@query vc.label.numberOfLines = 2
@query vc.label.text = "Hello Swift\nHello Julia"

sleep(0.5)

@query vc.label.font = UIFont(name: "Helvetica", size: 35)

sleep(0.5)

@query vc.label.backgroundColor = UIColor.cyanColor()

sleep(0.5)

@query vc.view.subviews[0].backgroundColor = UIColor.orangeColor()

sleep(0.5)

(Left, Center, Right) = (0, 1, 2)
@query vc.label.textAlignment = Center

sleep(0.5)

@query vc.label.frame = "{{20, 300}, {280, 200}}"

sleep(0.5)

max_stars_repo_name: wookay/Console

max_issues_repo_name: wookay/Console

max_forks_repo_name: wookay/Console


text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v01.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice_pwscf.jl")

include("structure_factor.jl")
include("gen_dr.jl")
include("gen_rho.jl")
include("calc_ewald.jl")

function test_main()
    Ns = [64, 64, 64]
    LatVecs = gen_lattice_sc(16.0)
    #
    pw = PWGrid( Ns, LatVecs )
    #
    Npoints = pw.Npoints
    Ω = pw.Ω
    r = pw.r
    G = pw.G
    G2 = pw.G2
    Ns = pw.Ns

    @printf("dr, dVol = %10.5e %10.5e\n", (Ω/Npoints)^(1/3), Ω/Npoints)

    # Atomic positions and nuclear charge
    Xpos = reshape( [0.0, 0.0, 0.0], (3,1) )

    println("Xpos =")
    println(Xpos)

    Z = 1.0
    # Calculate structure factor
    Sf = structure_factor( Xpos, G )

    println("sum(Sf) = ", sum(Sf))

    E_nn = calc_ewald( pw, Xpos, Sf, sigma=0.3 )
    @printf("E_nn = %18.10f\n", E_nn)
end

@time test_main()


text:
# -------- Benchmarks --------
struct Benchmark
    label::String
    words::Int64 # bytes = words * sizeof(Float64) * N
    flops::Int64
    write_alloc_factor::Float64 # bytes = words * write_alloc_factor * sizeof(Float64) * N
end

const BENCHMARKS = [
    Benchmark("Init", 1, 0, 2.0),
    Benchmark("Copy", 2, 0, 3/2),
    Benchmark("Update", 2, 1, 1.0),
    Benchmark("Triad", 3, 2, 4/3),
    Benchmark("Daxpy", 3, 2, 1.0),
    Benchmark("STriad", 4, 2, 5/4),
    Benchmark("SDaxpy", 4, 2, 1.0)
]

const NBENCH = length(BENCHMARKS)

# -------- Kernel (first `nthreads` threads) --------
function init_kernel(a, scalar; thread_indices, nthreads)
    @threads :static for tid in 1:nthreads
        @inbounds for i in thread_indices[tid]
            a[i] = scalar
        end
    end
    return nothing
end

function copy_kernel(a, b; thread_indices, nthreads)
    @threads :static for tid in 1:nthreads
        @inbounds for i in thread_indices[tid]
            a[i] = b[i]
        end
    end
    return nothing
end

function update_kernel(a, scalar; thread_indices, nthreads)
    @threads :static for tid in 1:nthreads
        @inbounds for i in thread_indices[tid]
            a[i] = a[i] * scalar
        end
    end
    return nothing
end

function triad_kernel(a, b, c, scalar; thread_indices, nthreads)
    @threads :static for tid in 1:nthreads
        @inbounds for i in thread_indices[tid]
            a[i] = b[i] + scalar * c[i]
        end
    end
    return nothing
end

function daxpy_kernel(a, b, scalar; thread_indices, nthreads)
    @threads :static for tid in 1:nthreads
        @inbounds for i in thread_indices[tid]
            a[i] = a[i] + scalar * b[i]
        end
    end
    return nothing
end

function striad_kernel(a, b, c, d; thread_indices, nthreads)
    @threads :static for tid in 1:nthreads
        @inbounds for i in thread_indices[tid]
            a[i] = b[i] + d[i] * c[i]
        end
    end
    return nothing
end

function sdaxpy_kernel(a, b, c; thread_indices, nthreads)
    @threads :static for tid in 1:nthreads
        @inbounds for i in thread_indices[tid]
            a[i] = a[i] + b[i] * c[i]
        end
    end
    return nothing
end

# -------- Kernel (all threads) --------
function init_kernel_allthreads(a, scalar)
    @threads :static for i in eachindex(a)
        @inbounds a[i] = scalar
    end
    return nothing
end

function copy_kernel_allthreads(a, b)
    @threads :static for i in eachindex(a)
        @inbounds a[i] = b[i]
    end
    return nothing
end

function update_kernel_allthreads(a, scalar)
    @threads :static for i in eachindex(a)
        @inbounds a[i] = a[i] * scalar
    end
    return nothing
end

function triad_kernel_allthreads(a, b, c, scalar)
    @threads :static for i in eachindex(a)
        @inbounds a[i] = b[i] + scalar * c[i]
    end
    return nothing
end

function daxpy_kernel_allthreads(a, b, scalar)
    @threads :static for i in eachindex(a)
        @inbounds a[i] = a[i] + scalar * b[i]
    end
    return nothing
end

function striad_kernel_allthreads(a, b, c, d)
    @threads :static for i in eachindex(a)
        @inbounds a[i] = b[i] + d[i] * c[i]
    end
    return nothing
end

function sdaxpy_kernel_allthreads(a, b, c)
    @threads :static for i in eachindex(a)
        @inbounds a[i] = a[i] + b[i] * c[i]
    end
    return nothing
end

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
using BenchmarkTools
using VectorDataUtils

#Non vectorised Lp norms and metrics

l2inner_alt(x,y) = sum(x.*y)
# sql2norm_alt(x) = sum(x.*x)
l2norm_alt(x) = sqrt(sum(x.*x))
sql2dist_alt(x,y) = sum((x-y).*(x-y))
l2dist_alt(x,y) = sqrt(sql2dist(x,y))
l1norm_alt(x) = sum(abs.(x))
l1dist_alt(x,y) = sum(abs.(x-y))

function sql2norm_alt(a::AbstractVector{T}) where {T}
    s = zero(T)
    @simd for i ∈ eachindex(a)
        s += a[i] * a[i]
    end
    return s
end

function return_alt(f1::Function)
    f2_sym = Symbol(String(Symbol(f1))*"_alt")
    f2 = getfield(Main,f2_sym)
    return f2
end

f1 = l1norm

ndim = 15
Ftype = Float32
x = randn(Ftype,ndim)
y = randn(Ftype,ndim)
ip = (x,)

f2 = return_alt(f1)
m1 = median(@benchmark f1(ip...))
m2 = median(@benchmark f2(ip...))
f1(ip...) ≈ f2(ip...)
judge(m1,m2)

function fmatch(f1::Function,ip)
    f2 = return_alt(f1)
    return f1(ip...) ≈ f2(ip...)
end

function fcomp(f1::Function,ip)
    f2 = return_alt(f1)
    b1 = @benchmarkable f1(ip...)
    b2 = @benchmarkable f2(ip...)
    tune!(b1)
    tune!(b2)
    t1 = run(b1)
    t2 = run(b2)
    m1 = median(t1)
    m2 = median(t2)
    return m1,m2
end

function bench_median(f::Function,ip)
    b1 = @benchmark f(ip...)
    return median(b1)
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 653cd9ec-01a9-11eb-2e64-3f7ae853eb87
begin
	using Yao, YaoPlots
	plot(chain(1, put(1=>X)))
end

# ╔═╡ 6c350302-00b5-11eb-0537-cd523683f91c
md"# *Qubits*
Like the *bits* of a classical computer, quantum computers have their own fundamental unit of data called, a *qubit*."

# ╔═╡ 0260bd12-00b6-11eb-19e9-81c90f08fe7b
md"Qubits can make use of some quantum mechanical properties like superposition. For the sake of understanding how qubits work, imagine the quantum computer as a box, full of qubits. Now these qubits are objects which have some mathematical and physical properties, and can store data and can be used to manipulate data to get some computation done.

Bits are represented by two states, 0 and 1. At any time, a qubit is in a *superposition* of two states, represented by $ \; a|\;0> +\;b|\;1> $. When we *measure* a qubit, its state *collapses* to, either $ |\;0> $ or $ |\;1> $. The chances\(probability\) of a qubit collapsing to the state, $ |\;0>\; $ is $ a^2 $ and to $ |\;1>\; $ is $ b^2 $. Hence, it must be satisfy $ |a^2| + |b^2| = 1 $.
a and b are also known as *probability amplitudes* of a given qubit."

# ╔═╡ df131fc4-00b9-11eb-09ed-9dcb647f6fda
md"**Note** : The notation of $ \;|> and <|$ are known as Dirac's notation, or the bra-ket notation.
For this tutorial, you need to understand that column vectors\(matrices of size n × 1\) are called kets and are represented by $ \;|> $ , and row vectors\(matrices of size 1 × n\) are called bras, and are represented by $ <| \;$ ."

# ╔═╡ 91768fe2-019c-11eb-3ac7-7747249e9dca
md"# Working with qubits
So what do the terms, superposition and measurement, mean...? Well, to keep it simple, superposition of two or more states just means their linear combination. The superposition of states $ |\;0> and\; |\;1> $ just means a linear combination of $ |\;0> and\; |\;1> $. If you still don't get it, no problem! Your intuition about superposition will build up by the time this tutorial reaches multiple qubits\(probably\).

When we use terms like measurement, we refer to *looking* at the state of the qubit. And yeah, we can't know the state of the qubit, because just when we look at it, it's state changes to $ |\;0> or\; |\;1>, $ and all the information about that qubit is lost. So the state of the qubit\(the values of a and b\) can't be determined."

# ╔═╡ c48bb1c0-019f-11eb-1210-e148231b166c
md"## Quantum Gates and Circuits
So consider that the state of a qubit is $ a|\;0> +\; b|\;1> $. We want to do something with this..... say, we want to change the state to $ a|\;1> +\; b|\;0> $... how do we manipulate the qubits or perform any operation on them?

We use quantum gates, the building blocks of quantum circuit, to manipulate qubits to get some task done.

Note: For the information that follows, please note that we can represent qubits as a column vector. Consider the qubit, $ a|\;0> +\; b|\;1> $, we can represent this as
$ \begin{bmatrix} a \\\ b \end{bmatrix} $."

# ╔═╡ c101c26c-01a7-11eb-03b2-75598cd2d3bc
md"### Single qubit gates
As the name suggests, these gates take one qubit for an input, and give a qubit with changed state for an output.
##### 1. The X gate
It changes the state of the qubit from $ a|\;0> +\; b|\;1> $ to $ a|\;1> +\; b|\;0> $. It *flips* the state of the qubit.
Mathematically, its represented by
$ \begin{bmatrix} 0 & 1 \\\ 1 & 0 \end{bmatrix} $, and applying this gate to a qubit is mathematically equivalent to multiplying the vector representing the qubit to the above matrix. It looks somewhat like this, when implemented in a circuit."

# ╔═╡ 0378bb12-01aa-11eb-2671-09822e968fc3
md"##### 2. The Y gate
It changes the state of the qubit from $ a|\;0> +\; b|\;1> $, to $ \;b|\;0> -\; a|\;1> $. It does a *bit flip* and a *phase flip* at the same time.
The Y gate is mathematically represented by the matrix,
$ i \begin{bmatrix} 0 & -1 \\\ 1 & 0 \end{bmatrix} $.
Mathematically, passing a qubit through this gate is equivalent to multiplying the state of the qubit, i.e.,
$ \begin{bmatrix} a \\\ b \end{bmatrix} $,
to the above matrix.
It's represented in a circuit by"

# ╔═╡ d0832d3a-0228-11eb-1bc0-95096891f407
plot(chain(1, put(1=>Y)))

# ╔═╡ 809f6516-0228-11eb-32b6-0389c49e9ec6
md"##### 3. The Z gate
It changes the state of the qubit from $ a|\;0> +\; b|\;1> $, to $ \;a|\;0> -\; b|\;1> $. It does a *sign flip* on the qubit.
The Z gate is mathematically represented by the matrix,
$ \begin{bmatrix} 1 & 0 \\\ 0 & -1 \end{bmatrix} $.
Mathematically, passing a qubit through this gate is equivalent to multiplying the state of the qubit, i.e.,
$ \begin{bmatrix} a \\\ b \end{bmatrix} $,
to the above matrix. It looks like this when implemented in a circuit."

# ╔═╡ 9bce5884-0229-11eb-2dc5-c75d4eaa20bf
plot(chain(1,put(1=>Z)))

# ╔═╡ 144beaee-022a-11eb-1c9d-d359ca75064f
md" **_Note_**: *The matrix representation of the above three gates are known as Pauli's matrices, represented by ${} \sigma_{x}, \;\sigma_y\;and\; \sigma_z$, for the X gate, the Y gate and the Z gate, respectively*"

# ╔═╡ 0ba83070-0230-11eb-30e9-d579362bda98
md"##### 4. The H gate
When a qubit is passed through H gate, the $ |\;0> $ changes to $ \frac{1}{\sqrt2}\(|\;0> +\; |\;1>\) $, and the $ |\;1> $ changes to $ \frac{1}{\sqrt2}\(|\;0> -\; |\;1>\) \\ $.

For an example, $ a|\;0> +\; b|\;1>, $ changes to, $ \frac{a}{\sqrt2}\(|\;0> +\; |\;1>\) +  \frac{b}{\sqrt2}\(|\;0> -\; |\;1>\)$. It can be simplified to give, $ \frac{a+b}{\sqrt2}\ | \;0> + $ $ \frac{a-b}{\sqrt2}\ | \;1> $.
Mathematically, its represented by the matrix $ \frac{1}{\sqrt2} \begin{bmatrix} 1 & 1 \\\ 1 & -1 \end{bmatrix} $ and passing a qubit through the H gate is mathematically equivalent to multiplying the vector representing the state of the qubit, to the above matrix. In a circuit, the H gate is represented by"

# ╔═╡ 0a36b65c-0233-11eb-212c-5d6b03794af1
plot(chain(1,put(1=>H)))

# ╔═╡ 70fcb098-0250-11eb-051b-ff526b9be775
md"**_Note_**: *The state $ \frac{1}{\sqrt2}\(|\;0> +\; |\;1>\) $ is often called $ |\;+> $, and the state $ \frac{1}{\sqrt2}\(|\;0> -\; |\;1>\) $ is often called $ |\;-> $.* "

# ╔═╡ 46c5e0d2-0233-11eb-387b-ad52dbad8553
md"### Multiqubit Gates
As the title suggests, these gates takes in and operate on more than one qubits."

# ╔═╡ 7edff156-0233-11eb-024b-090b7ae6481a
md"##### The CNOT gate
It takes two qubits as an input, and if the first qubit is a $ |\;1> $ the state of the second qubit is flipped to $ |\;0> $ if it was $ |\;1> $, and $ |\;0> $ if it was $ |\;1> $. If the state of the first qubit is $ |\;0> $, no change is made to the second qubit.
Mathematically, its represented by the matrix $ \begin{bmatrix} 1 & 0 & 0 & 0 \\\ 0 & 1 & 0 & 0 \\\ 0 & 0 & 0 & 1 \\\ 0 & 0 & 1 & 0 \end{bmatrix} $, and passing two qubits through it is equivalent to first collecting the vector of two two qubits into one and multiplying it to the above matrix.
In a circuit, the CNOT or the **CX** gate is represented by"

# ╔═╡ f21da5fc-024f-11eb-355b-09ce90c6f7f6
plot(chain(2,control(1,2=>X)))

# ╔═╡ a9ab153c-0256-11eb-093f-71fa576c4034
md" **_Note:_** *To represent a system of two qubits, ${}a_1 |\;00> +\; a_2 |\;01> +\;\; a_3 |\;10> +\; a_4 |\;11>$ in vector form, we can write them as, ${} \begin{bmatrix} a_1 \\a_2 \\a_3\\a_4 \end{bmatrix}$, where ${} |a_1^2| + |a_2^2| + |a_3^2| + |a_4^2| = 1.$ Also, we can use the CNOT or Control NOT or any other control gate(CY and CZ) to entangle two qubits.*"

# ╔═╡ Cell order:
# ╟─6c350302-00b5-11eb-0537-cd523683f91c
# ╟─0260bd12-00b6-11eb-19e9-81c90f08fe7b
# ╟─df131fc4-00b9-11eb-09ed-9dcb647f6fda
# ╟─91768fe2-019c-11eb-3ac7-7747249e9dca
# ╟─c48bb1c0-019f-11eb-1210-e148231b166c
# ╟─c101c26c-01a7-11eb-03b2-75598cd2d3bc
# ╟─653cd9ec-01a9-11eb-2e64-3f7ae853eb87
# ╟─0378bb12-01aa-11eb-2671-09822e968fc3
# ╟─d0832d3a-0228-11eb-1bc0-95096891f407
# ╟─809f6516-0228-11eb-32b6-0389c49e9ec6
# ╟─9bce5884-0229-11eb-2dc5-c75d4eaa20bf
# ╟─144beaee-022a-11eb-1c9d-d359ca75064f
# ╟─0ba83070-0230-11eb-30e9-d579362bda98
# ╟─0a36b65c-0233-11eb-212c-5d6b03794af1
# ╟─70fcb098-0250-11eb-051b-ff526b9be775
# ╟─46c5e0d2-0233-11eb-387b-ad52dbad8553
# ╟─7edff156-0233-11eb-024b-090b7ae6481a
# ╟─f21da5fc-024f-11eb-355b-09ce90c6f7f6
# ╟─a9ab153c-0256-11eb-093f-71fa576c4034

max_stars_repo_name: logankilpatrick/tutorials

max_issues_repo_name: logankilpatrick/tutorials

max_forks_repo_name: logankilpatrick/tutorials


text:
#JOB FILEPATHS
hmm_output = "/bench/PhD/NGS_binaries/BBM/survey_chains"
refined_path = "/bench/PhD/NGS_binaries/BBM/refined_chains"
sample_output = "/bench/PhD/NGS_binaries/BBM/refinement_samples"

survey_folders = "/bench/PhD/NGS_binaries/BBM/survey_folders"

#GENERAL SETUP
@info "Loading libraries..."
using BioBackgroundModels, DataFrames, Distributed, Serialization

const delta_thresh=1e-4 #stopping/convergence criterion (log probability difference btw subsequent EM iterates)
const max_iterates=15000

#LOAD SAMPLES
@info "Loading samples from $sample_output..."
sample_dfs = deserialize(sample_output)

#BUILD TRAINING AND TEST SETS FROM SAMPLES
training_sets, test_sets = split_obs_sets(sample_dfs)

#GENERATE Chain_ID Vector AND CHAINS SUBSET
@info "Loading reports from $survey_folders"
report_folders=deserialize(survey_folders)
job_ids=Vector{Chain_ID}()
for (partition,folder) in report_folders
    report=folder.partition_report
    for id in report.best_repset
        push!(job_ids,id)
    end
end

if isfile(refined_path) #if some results have already been collected, load them
    @info "Loading incomplete results..."
    refined_chains=deserialize(refined_path)
else
    @info "Initialising new HMM results file at $refined_path"
    chains = deserialize(hmm_output)
    refined_chains=Dict{Chain_ID,Vector{EM_step}}()
    for job_id in job_ids
        refined_chains[job_id]=[EM_step(1, chains[job_id][end].hmm, 0, 0, false)]
    end
end

#DISTRIBUTED CLUSTERs CONSTANTS
remote_machine = "10.0.0.3"
no_local_processes = 1
no_remote_processes = 1
load_dict=Dict{Int64,LoadConfig}()
local_config=LoadConfig(1:6,0:2)
remote_config=LoadConfig(1:6,0:2)

#SETUP DISTRIBUTED BAUM WELCH LEARNERS
@info "Spawning local cluster workers..."
worker_pool=addprocs(no_local_processes, topology=:master_worker)
for worker in worker_pool
    load_dict[worker]=local_config
end

remote_pool=addprocs([(remote_machine,no_remote_processes)], tunnel=true, topology=:master_worker)
for worker in remote_pool
    load_dict[worker]=remote_config
end

worker_pool=vcat(worker_pool, remote_pool)

@info "Loading worker libraries everywhere..."
@everywhere using BioBackgroundModels

em_jobset = setup_EM_jobs!(job_ids, training_sets; delta_thresh=delta_thresh, chains=refined_chains)
execute_EM_jobs!(worker_pool, em_jobset..., refined_path; delta_thresh=delta_thresh, max_iterates=max_iterates, load_dict=load_dict)


text:
using LyapunovExponents
using Plots
demo = solve!(LyapunovExponents.beer_95(); progress=1)
report(demo)
plt = plot(demo)


text:
# ---
# title: 961. N-Repeated Element in Size 2N Array
# id: problem961
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/n-repeated-element-in-size-2n-array/description/>
# hidden: true
# ---
#
# In a array `A` of size `2N`, there are `N+1` unique elements, and exactly one
# of these elements is repeated `N` times.
#
# Return the element repeated `N` times.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,3,3]
#     Output: 3
#
#
# **Example 2:**
#
#
#
#     Input: [2,1,2,5,3,2]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: [5,1,5,2,5,3,5,4]
#     Output: 5
#
#
#
#
# **Note:**
#
#   * `4 <= A.length <= 10000`
#   * `0 <= A[i] < 10000`
#   * `A.length` is even
#
#
## @lc code=start
using LeetCode

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


text:
__precompile__()

module Records

export
    Entity,
    Frame,
    EntityFrame,
    RecordFrame,
    RecordState,
    ListRecord,
    QueueRecord,
    EntityQueueRecord,

    ListRecordFrameIterator,
    ListRecordStateByIdIterator,

    get_statetype,
    get_deftype,
    get_idtype,

    capacity,
    nframes,
    nstates,
    nids,
    frame_inbounds,
    pastframe_inbounds,
    n_objects_in_frame,
    id2index,
    get_ids,
    nth_id,
    get_state,
    get_def,
    get_time,
    get_timestep,
    get_elapsed_time,
    get_subinterval,
    get_by_id,
    findfirst_stateindex_with_id,
    findfirst_frame_with_id,
    findlast_frame_with_id,
    get_first_available_id,
    push_back_records!,
    update!,
    allocate_frame,
    get_sparse_lookup

include("common.jl")
include("entities.jl")
include("frames.jl")
include("listrecords.jl")
include("queuerecords.jl")
include("conversions.jl")

end # module


text:
#Cylinder impedance BC
function Z(t , θ::Float64 , p::Param)

    Z0 = p.Z0
    ϵ  = p.ϵ
    Ω  = p.Ω

    return Z0*(1 + ϵ*cos(θ - Ω*t)^2)
end

function dZ(t,θ,p::Param)

    Z0 = p.Z0
    ϵ  = p.ϵ
    Ω  = p.Ω

    return 2 * Z0 * Ω * ϵ * cos(θ - Ω * t) * sin(θ - Ω * t)
end


text:
@testset "Persistency (naive coarse quantizer): save/load_ivfadc_index" begin
    ivfadc = build_index_random_data(coarse_quantizer=:naive)
    filepath, io = mktemp()

    # write to disk
    try
        save_ivfadc_index(filepath, ivfadc)

        # load from disk and compare
        ivfadc2 = load_ivfadc_index(filepath)
        @test typeof(ivfadc) == typeof(ivfadc2)

        # coarse quantizer
        @test ivfadc.coarse_quantizer.vectors == ivfadc2.coarse_quantizer.vectors
        # residual quantizer
        @test ivfadc.residual_quantizer.k== ivfadc2.residual_quantizer.k
        @test ivfadc.residual_quantizer.dims == ivfadc2.residual_quantizer.dims
        @test ivfadc.residual_quantizer.distance== ivfadc2.residual_quantizer.distance
        @test ivfadc.residual_quantizer.quantization == ivfadc2.residual_quantizer.quantization
        @test ivfadc.residual_quantizer.rot == ivfadc2.residual_quantizer.rot
        for i in 1:length(ivfadc.residual_quantizer.codebooks)
            @test ivfadc.residual_quantizer.codebooks[i].codes == ivfadc2.residual_quantizer.codebooks[i].codes
            @test ivfadc.residual_quantizer.codebooks[i].vectors == ivfadc2.residual_quantizer.codebooks[i].vectors
        end
        #index
        for i in 1:length(ivfadc.inverse_index)
            @test ivfadc.inverse_index[i].idxs == ivfadc2.inverse_index[i].idxs
            @test ivfadc.inverse_index[i].codes == ivfadc2.inverse_index[i].codes
        end
        rm(filepath, force=true, recursive=true)
    catch e
        rm(filepath, force=true, recursive=true)
        @test e
    end
end

@testset "Persistency (hnsw coarse quantizer)" begin
    ivfadc = build_index_random_data(coarse_quantizer=:hnsw)
    filepath, io = mktemp()
    # try to write to disk and test for failure
    try
        save_ivfadc_index(filepath, ivfadc)

        # load from disk and compare
        ivfadc2 = load_ivfadc_index(filepath)
        @test typeof(ivfadc) == typeof(ivfadc2)

        # brevity assignments
        hnsw_orig = ivfadc.coarse_quantizer.hnsw
        hnsw_loaded = ivfadc2.coarse_quantizer.hnsw
        # coarse quantizer
        @test hnsw_orig.lgraph.linklist == hnsw_loaded.lgraph.linklist
        @test hnsw_orig.lgraph.M == hnsw_loaded.lgraph.M
        @test hnsw_orig.lgraph.M0 == hnsw_loaded.lgraph.M0
        @test hnsw_orig.lgraph.m_L == hnsw_loaded.lgraph.m_L
        @test hnsw_orig.added == hnsw_loaded.added
        @test hnsw_orig.data == hnsw_loaded.data
        @test hnsw_orig.ep == hnsw_loaded.ep
        @test hnsw_orig.entry_level == hnsw_loaded.entry_level
        @test hnsw_orig.vlp.num_elements == hnsw_loaded.vlp.num_elements
        for i in 1:length(hnsw_orig.vlp.pool)
            @test hnsw_orig.vlp.pool[i].list == hnsw_loaded.vlp.pool[i].list
            @test hnsw_orig.vlp.pool[i].visited_value == hnsw_loaded.vlp.pool[i].visited_value
        end
        @test hnsw_orig.metric == hnsw_loaded.metric
        @test hnsw_orig.efConstruction == hnsw_loaded.efConstruction
        @test hnsw_orig.ef == hnsw_loaded.ef
        # residual quantizer
        @test ivfadc.residual_quantizer.k== ivfadc2.residual_quantizer.k
        @test ivfadc.residual_quantizer.dims == ivfadc2.residual_quantizer.dims
        @test ivfadc.residual_quantizer.distance== ivfadc2.residual_quantizer.distance
        @test ivfadc.residual_quantizer.quantization == ivfadc2.residual_quantizer.quantization
        @test ivfadc.residual_quantizer.rot == ivfadc2.residual_quantizer.rot
        for i in 1:length(ivfadc.residual_quantizer.codebooks)
            @test ivfadc.residual_quantizer.codebooks[i].codes == ivfadc2.residual_quantizer.codebooks[i].codes
            @test ivfadc.residual_quantizer.codebooks[i].vectors == ivfadc2.residual_quantizer.codebooks[i].vectors
        end
        #index
        for i in 1:length(ivfadc.inverse_index)
            @test ivfadc.inverse_index[i].idxs == ivfadc2.inverse_index[i].idxs
            @test ivfadc.inverse_index[i].codes == ivfadc2.inverse_index[i].codes
        end
        rm(filepath, force=true, recursive=true)
    catch e
        rm(filepath, force=true, recursive=true)
        @test e
    end
end


text:
function isvalidcircle(i,j,r,m,n,options)
	return (options.bounds_row[1] < i/m < options.bounds_row[2]) &&
		(options.bounds_column[1] < j/n < options.bounds_column[2]) &&
		(options.bounds_radius[1] < r/m < options.bounds_radius[2])
end

"""
	initvals(img,thresh=0.5)
Try to find good initial guesses for the cornea position in image `img`, which should be full color. The `thresh` parameter is the fraction of the max in the green channel intensity that a pixel should have to be considered part of the sclera. Returns a vector of zero to two tuples of (center_i,center_j,radius)
"""
function initvals(X::AbstractMatrix{T} where T <: AbstractRGB,purk=[];options=get_defaults())
	m,n = size(X)
	if isempty(purk)
		purk = findpurkinje(X;options)
	end

	# vertically, use the purkinje if possible
	i_c = length(purk) > 50 ? median(i[1] for i in purk) : m/2

	# horizontal guess 1: use the sclera's brightness
	G = green.(X)
	G[purk] .= 0  # ignore the purkinje

	idx = findall(G.>0.5*maximum(G))
	## get distribution peaks
	j_c = r_c = jleft = jright = wleft = wright = 0
	try
		(jleft,wleft),(jright,wright) = findpeaks(idx,(m,n),2)
		@debug "jleft,jright = $jleft,$jright"
		j_c = (jleft+jright)/2
		r_c = (jright-jleft)/2
	catch
	end

	## reject if the radius is unrealistic
	@debug "init1 = ($i_c,$j_c,$r_c)"
	guess = isvalidcircle(i_c,j_c,r_c,m,n,options) ? [(i_c,j_c,r_c)] : []

	# horizontal guess 2: use the purkinje location heuristically
	if length(purk) > 50
		j_p = median(i[2] for i in purk)
		γ = options.purkinje_location
		if wleft > wright
			r_c = (jright - j_p)/(1-γ)
			j_c = jright - r_c
		else
			r_c = (j_p-jleft)/(1+γ)
			j_c = jleft + r_c
		end
		@debug "init2 = ($i_c,$j_c,$r_c)"
		if isvalidcircle(i_c,j_c,r_c,m,n,options)
			push!(guess,(i_c,j_c,r_c))
		end
	end

	return guess
end

"""
	totalgrad(icen,jcen,r,Z[,trange])
Computes the sum of the radial derivative in an image along the specified arc of a circle. `Z` must be a callable function that returns numerical values for real arguments (i.e., allow interpolation). If given, `trange` is a 2-vector defining the range of angles to be used, measured ccw from "straight down" in the usual image visualization (i.e., vertically flipped).
"""
totalgrad(ic,jc,r,Z,θ) = sum(circlegrad(ic,jc,r,Z,θ))
function circlegrad(ic,jc,r,Z,θ)
	g = zeros(size(θ))
	p2c = (r,t) -> (ic+r*cos(t),jc+r*sin(t))
    for (k,t) in enumerate(θ)
		Zplus = Z(p2c(r+2,t)...)
		Zminus = Z(p2c(r-2,t)...)
		g[k] = Zplus - Zminus
	end
    return g
end

"""
	fitcircle(imgfun,i0,j0,r0[,t_range][;options=missing])
Find the circle that optimizes the total r-difference criterion, given initial values for the center and radius. Optionally specify the range of angles to be used (measured from straight down in the image) and options for the `optimize` solver.
"""
function fitcircle(imgfun,m,n,i0,j0,r0,method,θ=[-π,π];options=get_defaults())
	# don't let the center get too close to the edges, or the radius get too small or too close to the nearest edge
	function penalty(u)
		f = (x,lo,hi) -> log( max(1e-12,4(x-lo)*(hi-x)/(hi-lo)^2) )^4
		i,j,r = u
		b1 = options.bounds_row
		b2 = options.bounds_column
		b3 = options.bounds_radius
		return f(i/m,b1...) + f(j/n,b2...) + f(r/m,b3...)
	end

	if length(θ)==2  # given a range
		# space angles at 1 degree, regardless of the size of the range
		θ = filter( t->θ[1] ≤ t ≤ θ[2], 2π*(-180:180)/360)
	end

	objective = u -> -totalgrad(u...,imgfun,θ) + 0.1penalty(u)
	res = optimize(objective,[i0,j0,r0],method,Optim.Options(x_tol = 1e-2,f_tol = 1e-8))
	if !Optim.converged(res)
		@warn "Optimization did not converge"
		#show(res)
	end
	return res.minimizer,res.minimum
end

max_stars_repo_name: tobydriscoll/cornea-detection

max_issues_repo_name: tobydriscoll/cornea-detection

max_forks_repo_name: tobydriscoll/cornea-detection


text:
module SuffixTrees

export Node, SuffixTree, dotraversal, getlongestrepeatedsubstring

const oo = typemax(Int)

"""
    mutable struct Node
The suffix-tree's node.
Note that these are referenced not by pointer but by an index into an array of nodes.
"""
mutable struct Node
    children::Dict{Char, Int}
    start::Int
    ending::Int
    suffixlink::Int
    suffixindex::Int
end

Node() = Node(Dict(), 0, oo, 0, -1)
Node(start, ending) = Node(Dict(), start, ending, 0, -1)

""" Ukkonen Suffix-Tree """
mutable struct SuffixTree
    nodes::Vector{Node}
    text::Vector{Char}
    root::Int
    position::Int
    currentnode::Int
    needsuffixlink::Int
    remainder::Int
    activenode::Int
    activelength::Int
    activeedge::Int
end

""" length of edge of a node (the portion of the sequence it covers) """
edgelength(st, n::Node) = min(n.ending, st.position + 1) - n.start

""" make a node for the tree """
function newnode(st, start, ending)
    st.currentnode += 1
    st.nodes[st.currentnode] = Node(start, ending)
    return st.currentnode
end

""" constructor for Ukkonen Suffix-Tree """
function SuffixTree(str::String)
    nodes = [Node() for _ in 1:length(str) * 2]
    st = SuffixTree(nodes, [c for c in str], 1, 0, 0, 0, 0, 1, 1, 1)
    st.root = newnode(st, 0, 0)
    st.activenode = st.root
    for i in 1:length(st.text)
        extendsuffixtree(st, i)
    end
    setsuffixindexbyDFS(st, st.nodes[st.root], 0)
    return st
end

""" add a link to tree """
function addsuffixlink(st, nodenum::Int)
    if st.needsuffixlink > 0
        st.nodes[st.needsuffixlink].suffixlink = nodenum
    end
    st.needsuffixlink = nodenum
end

activeedge(st) = st.text[st.activeedge]

""" walk down the tree to its active length """
function walkdown!(st, currnode::Int)
    len = edgelength(st, st.nodes[currnode])
    st.activelength < len && return false
    st.activeedge += len
    st.activelength -= len
    st.activenode = currnode
    return true
end

""" extend tree by a character """
function extendsuffixtree(st, pos)
    st.position = pos
    st.needsuffixlink = 0
    st.remainder += 1
    while st.remainder > 0
        st.activelength == 0 && (st.activeedge = st.position)
        if !haskey(st.nodes[st.activenode].children, activeedge(st))
            nodenum = newnode(st, st.position, oo)
            st.nodes[st.activenode].children[activeedge(st)] = nodenum
            addsuffixlink(st, st.activenode)
        else
            next = st.nodes[st.activenode].children[activeedge(st)]
            walkdown!(st, next) && continue
            if st.text[st.nodes[next].start + st.activelength] == st.text[pos]
                addsuffixlink(st, st.activenode)
                st.activelength += 1
                break
            end
            splt = newnode(st, st.nodes[next].start, st.nodes[next].start + st.activelength)
            st.nodes[st.activenode].children[activeedge(st)] = splt
            nodenum = newnode(st, st.position, oo)
            st.nodes[splt].children[st.text[pos]] = nodenum
            st.nodes[next].start += st.activelength
            st.nodes[splt].children[st.text[st.nodes[next].start]] = next
            addsuffixlink(st, splt)
        end
        st.remainder -= 1
        if st.activenode == st.root && st.activelength > 0
            st.activelength -= 1
            st.activeedge = st.position - st.remainder + 1
        elseif st.activenode != st.root
            st.activenode = st.nodes[st.activenode].suffixlink
        end
    end
end

""" set the index of the leaves of the tree within the sequence """
function setsuffixindexbyDFS(st, node, labelheight, verbose=false)
    verbose && node.start > 0 && print(st.text[node.start:min(node.ending, length(st.text))])
    isleaf = true
    for child in map(v -> st.nodes[v], collect(values(node.children)))
        verbose && isleaf && node.start > 0 && println(" [", node.suffixindex, "]")
        isleaf = false
        setsuffixindexbyDFS(st, child, labelheight + edgelength(st, child))
    end
    if isleaf
        idx = length(st.text) - labelheight
        node.suffixindex = idx
        verbose && println(" [$idx]")
    end
end

""" traverse the suffix tree """
function dotraversal(st)
    maxheight, substringstartindices = 0, [0]
    function traversal(node::Node, labelheight)
        if node.suffixindex == -1
            for child in map(v -> st.nodes[v], collect(values(node.children)))
                traversal(child, labelheight + edgelength(st, child))
            end
        elseif maxheight < labelheight - edgelength(st, node)
            maxheight = labelheight - edgelength(st, node)
            substringstartindices = [node.suffixindex + 1]
        elseif maxheight == labelheight - edgelength(st, node)
            push!(substringstartindices, node.suffixindex + 1)
        end
    end
    traversal(st.nodes[st.root], 0)
    return maxheight, substringstartindices
end

""" find the longest repeated suffix of the tree """
function getlongestrepeatedsubstring(st::SuffixTree, label="", printresult=true)
    len, starts = dotraversal(st)
    substring = len == 0 ? "" :
        join(unique(map(x -> String(st.text[x:x+len-1]), starts)), " (or) ")
    if printresult
        print("  ", label == "" ? String(st.text) : label, ": ")
        println(len == 0 ? "No repeated substring." : substring)
    end
    return substring
end

end # module


text:
#: Utilities for norm functions

export ptdot, ptnorm1, ptnorm2, ptnormInf
export mixed_norm, gradient_norm, gradient_norm_batch
export proj_options

# Options for projection-related operations

struct ProjOptions
    maxiter::Union{Nothing,Int64}
    xrtol::Union{Nothing,Real}
end

proj_options(; maxiter::Union{Nothing,Int64}=nothing, xrtol::Union{Nothing,Real}=nothing) = ProjOptions(maxiter, xrtol)

# Mixed norm

struct MixedNorm{T,D,N1,N2}<:ProximableFunction{T,D}
    proj_opt::ProjOptions
end

function mixed_norm(D::Number, N1::Number, N2::Number; T::DataType=Float32, proj_opt::Union{Nothing,ProjOptions}=nothing)
    proj_opt === nothing && (proj_opt = proj_options())
    (D == 2) && return MixedNorm{T,3,N1,N2}(proj_opt)
    (D == 3) && return MixedNorm{T,4,N1,N2}(proj_opt)
    throw("Dimension not supported")
end

# Mixed norm (batch)

struct MixedNormBatch{T,D,N1,N2}<:ProximableFunction{T,D}
    proj_opt::ProjOptions
end

function mixed_norm_batch(D::Number, N1::Number, N2::Number; T::DataType=Float32, proj_opt::Union{Nothing,ProjOptions}=nothing)
    proj_opt === nothing && (proj_opt = proj_options())
    (D == 2) && return MixedNormBatch{T,4,N1,N2}(proj_opt)
    (D == 3) && return MixedNormBatch{T,5,N1,N2}(proj_opt)
    throw("Dimension not supported")
end

# L22 norm

function proxy!(p::AbstractArray{CT,N}, λ::T, ::MixedNorm{CT,N,2,2}, q::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}
    np = norm22(p)
    np <= λ ? (return q .= CT(0)) : (return q .= (1-λ/np)*p)
end

function project!(p::AbstractArray{CT,N}, ε::T, ::MixedNorm{CT,N,2,2}, q::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}
    np = norm22(p)
    np <= ε ? (return q .= p) : (return q .= ε*p/np)
end

(::MixedNorm{T,N,2,2})(p::AbstractArray{T,N}) where {T,N} = norm22(p)

# L21 norm

function proxy!(p::AbstractArray{CT,N}, λ::T, ::MixedNorm{CT,N,2,1}, q::AbstractArray{CT,N}; ptn::Union{AbstractArray{T,N},Nothing}=nothing) where {T<:Real,N,CT<:RealOrComplex{T}}
    ptn === nothing && (ptn = ptnorm2(p; η=eps(T)))
    return q .= (1 .-λ./ptn).*(ptn .>= λ).*p
end

function project!(p::AbstractArray{CT,N}, ε::T, g::MixedNorm{CT,N,2,1}, q::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}
    ptn = ptnorm2(p; η=eps(T))
    sum(ptn) <= ε && (return q .= p)
    λ = pareto_search_proj21(ptn, ε; xrtol=g.proj_opt.xrtol, maxiter=g.proj_opt.maxiter)
    return proxy!(p, λ, g, q; ptn=ptn)
end

(::MixedNorm{T,N,2,1})(p::AbstractArray{T,N}) where {T,N} = norm21(p)

# L21 norm (batch, 2-D)

function proxy!(p::AbstractArray{CT,4}, λ::T, ::MixedNormBatch{CT,4,2,1}, q::AbstractArray{CT,4}; ptn::Union{AbstractArray{CT,4},Nothing}=nothing) where {T<:Real,CT<:RealOrComplex{T}}
    ptn === nothing && (ptn = ptnorm2_batch(p; η=eps(T)))
    nx,ny,nc,nb = size(p)
    p = reshape(p, nx,ny,2,div(nc,2)*nb)
    q = reshape(q, nx,ny,2,div(nc,2)*nb)
    ptn = reshape(ptn, nx,ny,1,div(nc,2)*nb)
    q .= (CT(1).-λ./ptn).*(ptn .>= λ).*p
    return reshape(q, nx,ny,nc,nb)
end

(::MixedNormBatch{T,4,2,1})(p::AbstractArray{T,4}) where T = norm21_batch(p)

# L2Inf norm

function proxy!(p::AbstractArray{CT,N}, λ::T, ::MixedNorm{CT,N,2,Inf}, q::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}
    project!(p, λ, mixed_norm(N-1, 2, 1; T=CT), q)
    return q .= p.-q
end

function project!(p::AbstractArray{CT,N}, ε::T, ::MixedNorm{CT,N,2,Inf}, q::AbstractArray{CT,N}) where {T<:Real,N,CT<:RealOrComplex{T}}
    ptn = ptnorm2(p; η=eps(T))
    val = ptn .>= ε
    q .= p.*(ε*val./ptn+(!).(val))
    return q
end

(::MixedNorm{T,N,2,Inf})(p::AbstractArray{T,N}) where {T,N} = norm2Inf(p)

# Gradient mixed norm

function gradient_norm(N1::Number, N2::Number, n::NTuple{D,Int64}, h::NTuple{D,S}; T::DataType=Float32, weight::Union{Nothing,AbstractLinearOperator}=nothing, proj_opt::Union{Nothing,ProjOptions}=nothing) where {D,S<:Real}
    ∇ = gradient_operator(n, h; T=T)
    weight !== nothing ? (A∇ = weight*∇) : (A∇ = ∇)
    return mixed_norm(D, N1, N2; T=T, proj_opt=proj_opt)∘A∇
end

# Gradient mixed norm (2-D, batch)

function gradient_norm_batch(N1::Number, N2::Number, n::NTuple{D,Int64}, nc::Int64, nb::Int64, h::NTuple{2,S}; T::DataType=Float32, weight::Union{Nothing,AbstractLinearOperator}=nothing, proj_opt::Union{Nothing,ProjOptions}=nothing) where {S<:Real,D}
    ∇ = gradient_operator_batch(n, nc, nb, h; T=T)
    weight !== nothing ? (A∇ = weight*∇) : (A∇ = ∇)
    return mixed_norm_batch(D, N1, N2; T=T, proj_opt=proj_opt)∘A∇
end

# Pareto-search routines

pareto_search_proj21(ptn::AbstractArray{T,N}, ε::T; maxiter::Union{Nothing,Int64}=nothing, xrtol::Union{Nothing,T}=nothing) where {T<:Real,N} = T(solve(ZeroProblem(λ -> obj_pareto_search_proj21(λ, ptn, ε), (T(0), maximum(ptn))), Roots.Brent(); xreltol=isnothing(xrtol) ? eps(T) : xrtol, maxevals=isnothing(maxiter) ? length(ptn) : maxiter))

obj_pareto_search_proj21(λ::T, ptn::AbstractArray{T,N}, ε::T) where {T<:Real,N} = sum(Flux.relu.(ptn.-λ))-ε

# Algebraic utils

ptdot(v1::AbstractArray{T,N}, v2::AbstractArray{T,N}) where {T,N} = sum(v1.*conj.(v2); dims=N)
ptnorm1(p::AbstractArray{CT,N}; η::T=T(0)) where {T<:Real,N,CT<:RealOrComplex{T}} = sum(abs.(p).+η; dims=N)
ptnorm2(p::AbstractArray{CT,N}; η::T=T(0)) where {T<:Real,N,CT<:RealOrComplex{T}} = sqrt.(sum(abs.(p).^2 .+η^2; dims=N))
ptnormInf(p::AbstractArray{CT,N}; η::T=T(0)) where {T<:Real,N,CT<:RealOrComplex{T}} = sqrt.(maximum(abs.(p).+η; dims=N))
function ptnorm2_batch(p::AbstractArray{CT,4}; η::T=T(0)) where {T<:Real,CT<:RealOrComplex{T}}
    nx,ny,nc,nb = size(p)
    p = reshape(p, nx,ny,2,div(nc,2)*nb)
    return reshape(sqrt.(abs.(p[:,:,1:1,:]).^2+abs.(p[:,:,2:2,:]).^2 .+η^2), nx,ny,div(nc,2),nb)
end
function ptnorm2_batch(p::AbstractArray{CT,5}; η::T=T(0)) where {T<:Real,CT<:RealOrComplex{T}}
    nx,ny,nz,nc,nb = size(p)
    p = reshape(p, nx,ny,nz,2,div(nc,2)*nb)
    return reshape(sqrt.(abs.(p[:,:,1:1,:]).^2+abs.(p[:,:,2:2,:]).^2+abs.(p[:,:,3:3,:]).^2 .+η^2), nx,ny,div(nc,2),nb)
end
norm21(v::AbstractArray{CT,N}; η::T=T(0)) where {T<:Real,N,CT<:RealOrComplex{T}} = sum(ptnorm2(v; η=η))
norm22(v::AbstractArray{CT,N}; η::T=T(0)) where {T<:Real,N,CT<:RealOrComplex{T}} = sqrt(sum(ptnorm2(v; η=η).^2))
norm2Inf(v::AbstractArray{CT,N}; η::T=T(0)) where {T<:Real,N,CT<:RealOrComplex{T}} = maximum(ptnorm2(v; η=η))
function norm21_batch(v::AbstractArray{CT,4}; η::T=T(0)) where {T<:Real,CT<:RealOrComplex{T}}
    _,_,nc,nb = size(v)
    return reshape(sum(ptnorm2_batch(v; η=η); dims=(1,2)), div(nc,2), nb)
end
function norm21_batch(v::AbstractArray{CT,5}; η::T=T(0)) where {T<:Real,CT<:RealOrComplex{T}}
    _,_,_,nc,nb = size(v)
    return reshape(sum(ptnorm2_batch(v; η=η); dims=(1,2,3)), div(nc,2), nb)
end

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
"""
Object representing a possible root inside a given `interval`.
`status` is either `:unknown` or `:unique`.
If `status` is `:unique` then we know that there is a unique root of the
function in question inside the given `interval`.
"""
struct Root{T}
    interval::T
    status::Symbol
end

interval(x::Root) = x.interval

# use root_status since just `status` is too generic
root_status(x::Root) = x.status

show(io::IO, root::Root) = print(io, "Root($(root.interval), :$(root.status))")

isunique(root::Root{T}) where {T} = (root.status == :unique)

⊆(a::Interval, b::Root) = a ⊆ b.interval   # the Root object has the interval in the first entry
⊆(a::Root, b::Root) = a.interval ⊆ b.interval

big(a::Root) = Root(big(a.interval), a.status)


text:
module TestWithBase

using Base: splat
using Test
using ThreadsX

args_and_kwargs(args...; kwargs...) = args, (; kwargs...)

inc(x) = x + 1

mapi(args...) = map(args...)

raw_testdata = """
collect(1:10)
collect(Float64, 1:10)
collect(x for x in 1:10 if isodd(x))
collect(Float64, (x for x in 1:10 if isodd(x)))
collect(inc(x) for x in 1:10)
collect(Float64, (inc(x) for x in 1:10))
collect(x * y for x in 1:10, y in 11:20)
collect(Float64, (x * y for x in 1:10, y in 11:20))
collect(x * y for x in 1:0, y in 11:20)
collect(Float64, (x * y for x in 1:0, y in 11:20))
map(inc, 1:10)
map(inc, (x for x in 1:10 if isodd(x)))
map(inc, Float64[])
map(inc, ones(3, 3))
map(inc, ones(3, 0))
map(inc, ones(0, 3))
map(*, 1:10, 11:20)
map(*, ones(3, 3), ones(3, 3))
map(*, ones(3, 0), ones(3, 0))
map(*, ones(0, 3), ones(0, 3))
map(splat(*), Iterators.product(1:10, 11:20))
map(splat(*), Iterators.product(1:0, 11:20))
reduce(+, 1:10)
reduce(+, 1:0)
reduce(+, Bool[])
reduce(+, (x^2 for x in 1:10))
reduce(+, (x for x in 1:10 if isodd(x)))
reduce(+, (y for x in 1:10 for y in 1:x))
mapreduce(inc, +, 1:10)
mapreduce(*, +, 1:10, 11:20)
mapreduce(*, +, 1:0, 11:10)
mapreduce(inc, +, (x^2 for x in 1:10))
mapreduce(inc, +, (x for x in 1:10 if isodd(x)))
mapreduce(inc, +, (y for x in 1:10 for y in 1:x))
sum(1:10)
sum(1:0)
sum(x -> x^2, 1:10)
sum(x^2 for x in 1:10)
sum(x for x in 1:10 if isodd(x))
sum(y for x in 1:10 for y in 1:x)
prod(1:10)
prod(1:0)
prod(x -> x + 1, 1:10)
count(isodd.(1:10))
count(isodd, 1:10)
count(Bool[])
count(isodd, 1:0)
maximum(1:10)
maximum(inc, 1:10)
minimum(1:10)
minimum(inc, 1:10)
any(fill(false, 10))
any([fill(false, 10); true])
any([])
all(fill(true, 10))
all([fill(true, 10); false])
all([])
findfirst([fill(false, 10); true])
findfirst(iseven, 1:10)
findfirst(reshape([fill(false, 10); true; true], 3, 4))
findfirst(==(5), reshape(1:12, 3, 4))
findfirst([])
findfirst(identity, [])
findlast([fill(true, 10); false])
findlast(isodd, 1:10)
findlast(reshape([fill(true, 10); false; false], 3, 4))
findlast(isodd, reshape(1:12, 3, 4))
findlast([])
findlast(identity, [])
findall([fill(false, 10); true])
findall(iseven, 1:10)
findall(reshape([fill(false, 10); true; true], 3, 4))
findall(==(5), reshape(1:12, 3, 4))
findall([])
findall(identity, [])
extrema(1:10)
extrema(sin, 1:10)
issorted(0:9)
issorted(0:9; rev=true)
issorted(0:9; by=_ -> 1)
issorted(reverse(0:9))
issorted(reverse(0:9); rev=true)
issorted(reverse(0:9); by=_ -> 1)
issorted([])
issorted([]; rev=true)
issorted([]; by=_ -> 1)
unique([1, 2, 6, 2])
unique(Real[1, 1.0, 2])
unique(x -> x^2, [1, -1, -3, 4, 3])
Set([1, -1, -3, 4, 3])
Set(x^2 for x in [1, -1, -3, 4, 3])
"""

# An array of `(label, (f, args, kwargs))`
testdata = mapreduce(vcat, split(raw_testdata, "\n", keepempty = false)) do x
    @debug "Parsing: $x"
    f, rest = split(x, "(", limit = 2)
    ex = Meta.parse("DUMMY($rest")
    ex.args[1] = args_and_kwargs
    tests = Any[@eval ($x, ($(Symbol(f)), $ex...))]
    if f == "map"
        y = replace(x, r"^map" => "mapi")
        push!(tests, @eval ($y, (mapi, $ex...)))
    end
    return tests
end

@testset "$label" for (label, (f, args, kwargs)) in testdata
    g = getproperty(ThreadsX, nameof(f))
    @testset "default basesize" begin
        @test g(args...; kwargs...) == f(args...; kwargs...)
    end
    @testset for basesize in 1:3
        @test g(args...; kwargs..., basesize = basesize) == f(args...; kwargs...)
    end
end

function test_all_implementations(test, name)
    @testset "Base" begin
        test(getproperty(Base, name))
    end
    @testset "ThreadsX" begin
        f = getproperty(ThreadsX, name)
        @testset "default basesize" begin
            test(f)
        end
        @testset for basesize in 1:3
            test((args...) -> f(args...; basesize = basesize))
        end
    end
end

@testset "foreach(x -> ys[x] = x^2, 1:5)" begin
    test_all_implementations(:foreach) do foreach
        xs = 1:5
        ys = zero(xs)
        foreach(xs) do x
            ys[x] = x^2
        end
        @test ys == xs .^ 2
    end
end

@testset "foreach((i, x) -> ys[i] = x^2, eachindex(ys, xs), xs)" begin
    test_all_implementations(:foreach) do foreach
        xs = 11:15
        ys = zero(xs)
        foreach(eachindex(ys, xs), xs) do i, x
            ys[i] = x^2
        end
        @test ys == xs .^ 2
    end
end

@testset "foreach(..., product(1:2, 1:3))" begin
    test_all_implementations(:foreach) do foreach
        xs = Iterators.product(1:2, 1:3)
        ys = fill((-1, -1), 2, 3)
        foreach(xs) do I
            ys[I...] = I
        end
        @test ys == Tuple.(xs)
    end
end

@testset "foreach(..., product(1:2, 1:3, 1:4))" begin
    test_all_implementations(:foreach) do foreach
        xs = Iterators.product(1:2, 1:3, 1:4)
        ys = fill((-1, -1, -1), 2, 3, 4)
        foreach(xs) do I
            ys[I...] = I
        end
        @test ys == Tuple.(xs)
    end
end

@testset "map!(x -> x^2, ys, xs)" begin
    test_all_implementations(:map!) do map!
        xs = 11:15
        ys = zero(xs)
        map!(x -> x^2, ys, xs)
        @test ys == xs .^ 2
    end
end

@testset "map!(+, ys, xs1, xs2)" begin
    test_all_implementations(:map!) do map!
        xs1 = 1:5
        xs2 = xs1 .* 10
        ys = zero(xs1)
        map!(+, ys, xs1, xs2)
        @test ys == xs1 .+ xs2
    end
end

@testset "map!(+, dest, matrix, matrix')" begin
    test_all_implementations(:map!) do map!
        matrix = reshape(1:9, 3, 3)
        dest = zero(matrix)
        map!(+, dest, matrix, matrix')
        @test dest == matrix .+ matrix'
    end
end

end  # module


text:
# Optimization of bin counts for histograms, heatmaps, hexbin plots, etc.
#
# I'm using the penalized maximum-likelihood method proposed in
#   Birge, L, and Rozenholc, Y. (2006) How many bins should be put in a regular
#   histogram?
#
# There has been quite a bit written on this problem, but there are a number of
# methods that all seem to give good results with little difference. Birge's
# method is simple (it's just AIC with an extra logarithmic term), has a decent
# theoretical justification, and is general enough to apply to multidimensional
# and non-regular bin selecetion problems. Though, the penalty they use was
# optimized for regular histograms, so may need to be tweaked.
#
# The Birge penalty is
#    penalty(D) = D - 1 + log(D)^2.5
# where D is the number of bins. The 2.5 constant was arrived at emperically by
# optimizing over samples from example density functions.
#

# Penalized log-likelihood function for a histogram with d regular bins.
#
# Args:
#   d: Number of bins in the histogram.
#   n: Number of sample (which should equal sum(bincounts[1:d])).
#   bincounts: An array giving the number occurrences in each bin.
#   binwidth: Width of each bin in the histogram.
#
# Returns:
#   Log-likelihood with Birge's penalty applied.
#
function bincount_pll(d::Int, n::Int, bincounts::Vector{Int}, binwidth::Float64)
    ll = 0
    for i in 1:d
        if bincounts[i] > 0
            ll += bincounts[i] * log(bincounts[i] / (n * binwidth))
        end
    end
    ll - (d - 1 + log(d)^2.5)
end

# Optimize the number of bins for a regular one dimensional histogram.
#
# Args:
#   xs: A sample.
#
# Returns:
#   A tuple of the form (d, bincounts), where d gives the optimal number of
#   bins, and bincounts is an array giving the number of occurances in each bin.
#
function choose_bin_count_1d(xs::AbstractVector, min_bin_count=1, max_bin_count=150)
    n = length(xs)
    if n <= 1
        return 1, Int[0]
    end

    x_min, x_max = Gadfly.concrete_minimum(xs), Gadfly.concrete_maximum(xs)
    span = x_max - x_min

    d_min = min_bin_count
    d_max = max_bin_count
    bincounts = zeros(Int, d_max)

    d_best = d_min
    pll_best = -Inf

    # Brute force optimization: since the number of bins has to be reasonably
    # small to plot, this is pretty quick and very simple.
    for d in d_min:d_max
        binwidth = span / d
        bincounts[1:d] = 0

        for x in xs
            if !isconcrete(x)
                continue
            end
            bincounts[max(1, min(d, int(ceil((x - x_min) / binwidth))))] += 1
        end

        pll = bincount_pll(d, n, bincounts, binwidth)

        if pll > pll_best
            d_best = d
            pll_best = pll
        end
    end

    bincounts[1:d_best] = 0
    binwidth = span / d_best
    for x in xs
        if !isconcrete(x)
            continue
        end
        bincounts[max(1, min(d_best, int(ceil((x - x_min) / binwidth))))] += 1
    end

    (d_best, bincounts)
end

# Optimize the number of bins for regular two dimensional histograms.
#
# Args:
#   xs: Dimension one data.
#   ys: Dimension two data.
#
# Returns:
#   A tuple of the form (dx, dy, bincounts), where dx, dy gives the number of
#   bins in each respective dimension and bincounts is a dx by dy matrix giving
#   the count in each bin.
#
function choose_bin_count_2d(xs::AbstractVector, ys::AbstractVector,
                             xminbincount::Int, xmaxbincount::Int,
                             yminbincount::Int, ymaxbincount::Int)

    # For two demensions, I'm just going to optimize the marginal bin counts.
    # This might not be optimal, but its simple and fast.

    x_min, x_max = Gadfly.concrete_minimum(xs), Gadfly.concrete_maximum(xs)
    y_min, y_max = Gadfly.concrete_minimum(ys), Gadfly.concrete_maximum(ys)

    dx, _ = choose_bin_count_1d(xs, xminbincount, xmaxbincount)
    dy, _ = choose_bin_count_1d(ys, yminbincount, ymaxbincount)

    # bin widths
    wx = (x_max - x_min) / dx
    wy = (y_max - y_min) / dy

    bincounts = zeros(Int, (dy, dx))
    for (x, y) in zip(xs, ys)
        if !Gadfly.isconcrete(x) || !Gadfly.isconcrete(y)
            continue
        end

        i = max(1, min(dx, int(ceil((x - x_min) / wx))))
        j = max(1, min(dy, int(ceil((y - y_min) / wy))))
        bincounts[j, i] += 1
    end

    (dy, dx, bincounts)
end

# Optimize the number of bins for hexagonal 2d histograms.
#
# Args:
#   xs: Dimension one data.
#   ys: Dimension two data.
#
# Returns:
#   A tuple of the form (size, bincounts), where 'size' is the hexagon size and
#   bincounts is an arrays storing the count for each hexagon, using axial
#   coordinates.
#
function choose_hex_bin_count(xs::AbstractVector, ys::AbstractVector)
    # TODO: this should probably be a thing
end


text:
using DesignModels
using Documenter

makedocs(;
    modules=[DesignModels],
    authors="azzaare <jf@baffier.fr> and contributors",
    repo="https://github.com/azzaare/DesignModels.jl/blob/{commit}{path}#L{line}",
    sitename="DesignModels.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://azzaare.github.io/DesignModels.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/azzaare/DesignModels.jl.git",
)


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

@testset "SimData" begin
    @test_throws ArgumentError SimData(1M10, (:a, :b, :c), rand(10, 2))
    let nms = (:a, :b), dta = rand(20, 2), sd = SimData(2000Q1, nms, copy(dta)), dta2 = rand(size(dta)...)
        @test firstdate(sd) == 2000Q1
        @test lastdate(sd) == 2000Q1 + 20 - 1
        @test frequencyof(sd) == Quarterly
        @test sd isa AbstractMatrix
        @test size(sd) == size(dta)
        # integer indexing must be identical to dta
        for i in axes(dta, 2)
            @test sd[:,i] == dta[:,i]
        end
        for j in axes(dta, 1)
            @test sd[j,:] == dta[j,:]
        end
        for i in eachindex(dta)
            @test sd[i] == dta[i]
        end
        # set indexing with integers
        for i in eachindex(dta2)
            sd[i] = dta2[i]
        end
        @test sd[:] == dta2[:]
        for i in axes(dta, 1)
            sd[i,:] = dta[i,:]
        end
        @test sd[:] == dta[:]
        for i in axes(dta2, 2)
            sd[:,i] = dta2[:,i]
        end
        @test sd[:] == dta2[:]
        # access by dot notation (for the columns)
        @test propertynames(sd) == nms
        @test sd.a isa TSeries
        @test sd.b isa TSeries
        @test_throws BoundsError sd.c
        @test sd.a.values == dta2[:,1]
        sd.a[:] = dta[:,1]
        sd.b[:] = dta[:,2]
        @test sd[:] == dta[:]
        @test sd[:a] isa TSeries && sd[:a].values == dta[:,1]
        @test sd["b"] isa TSeries && sd["b"].values == dta[:,2]
        #
        sd.a = dta[:,1]
        sd.b = dta[:,2]
        @test sd[:] == dta[:]
        #
        sd[:] = dta[:]
        sd.a = sd.b
        @test sd[:,1] == sd[:,2]
        sd.a = zeros(size(dta, 1))
        @test sum(abs, sd.a.values) == 0
        @test_throws BoundsError sd.a = ones(length(sd.a) + 5)
        # access to rows by MIT
        sd[:] = dta[:]
        @test sd[2000Q1] isa NamedTuple{nms,NTuple{length(nms),Float64}}
        for (i, idx) in enumerate(mitrange(sd))
            @test [values(sd[idx])...] == dta[i,:]
        end
        sd[2000Q1] = zeros(size(dta, 2))
        @test sum(abs, sd[1,:]) == 0
        @test_throws BoundsError sd[1999Q4]
        @test_throws DimensionMismatch sd[2000Q3] = zeros(size(dta, 2) + 1)
        sd[2000Q2] = (b = 7.0, a = 1.5)
        @test sd[2,:] == [1.5, 7.0]
        sd[2004Q4] = (a = 5.0,)
        @test sd[end,1] == 5.0
        @test_throws BoundsError sd[2004Q4] = (a = 5.0, c = 1.1)
        @test_throws ArgumentError sd[2000Y]
        @test_throws BoundsError sd[1999Q4] = (a = 5.0, b = 1.1)
        @test_throws BoundsError sd[2010Q4] = rand(2)
        #
        @test sd[2001Q1:2002Q4] isa SimData
        @test sd[2001Q1:2002Q4] == sd[5:12,:]
        sd[2001Q1:2002Q4] = 1:16
        @test sd[5:12,:] == reshape(1.0:16.0, :, 2)
        @test_throws BoundsError sd[1111Q4:2002Q1]
        @test_throws BoundsError sd[2002Q1:2200Q2]
        @test_throws DimensionMismatch sd[2001Q1:2002Q4] = 1:17
        # assign new column
        sd1 = hcat(sd, c=sd.a + 3.0)
        @test sd1[nms] == sd
        @test sd1[(:a,:c)] == sd1[:,[1,3]]
        # access with 2 args MIT and Symbol
        @test sd[2001Q2, (:a, :b)] isa NamedTuple
        let foo = sd[2001Q2:2002Q1, (:a, :b)]
            @test foo isa SimData
            @test size(foo) == (4, 2)
            @test firstdate(foo) == 2001Q2
        end
        @test_throws BoundsError sd[1999Q1, (:a,)]
        @test_throws BoundsError sd[2001Q1:2001Q2, (:a, :c)]
        @test_throws Exception sd.c = 5
        @test similar(sd) isa typeof(sd)
        @test_throws BoundsError sd[1999Q1:2000Q4] = zeros(8, 2)
        @test_throws BoundsError sd[2004Q1:2013Q4, [:a, :b]]
        # setindex with two
        sd[2001Q1, (:b,)] = 3.5
        @test sd[5,2] == 3.5
        sd[2000Q1:2001Q4, (:b,)] = 3.7
        @test all(sd[1:8,2] .== 3.7)
        @test_throws BoundsError sd[1999Q1:2000Q4, (:a, :b)] = 5.7
        @test_throws BoundsError sd[2000Q1:2000Q4, (:a, :c)] = 5.7
    end
end

@testset "SimData show" begin
    # test the case when column labels are longer than the numbers
    let io = IOBuffer(), sd =  SimData(1U, (:verylongandsuperboringnameitellya, :anothersuperlongnamethisisridiculous, :a), rand(20, 3) .* 100)
        show(io, sd)
        lines = readlines(seek(io, 0))
        # labels longer than 10 character are abbreviated with '…' at the end
        @test length(split(lines[2], '…')) == 3
    end
    let io = IOBuffer() , sd = SimData(1U, (:alpha, :beta), zeros(24, 2))
        show(io, sd)
        lines = readlines(seek(io, 0))
        lens = length.(lines)
        # when labels are longer than the numbers, the entire column stretches to fit the label
        @test lens[2] == lens[3]
    end
    nrow = 24
    letters = Symbol.(['a':'z'...])
    for (nlines, fd) = zip([3, 4, 5, 6, 7, 8, 22, 23, 24, 25, 26, 30], Iterators.cycle((qq(2010, 1), mm(2010, 1), yy(2010), ii(1))))
        for ncol = [2,5,10,20]
            # display size if nlines × 80
            # data size is nrow × ncol
            # when printing data there are two header lines - summary and column names
            io = IOBuffer()
            sd = SimData(fd, tuple(letters[1:ncol]...), rand(nrow, ncol))
            show(IOContext(io, :displaysize => (nlines, 80)), MIME"text/plain"(), sd)
            lines = readlines(seek(io, 0))
            @test length(lines) == max(3, min(nrow + 2, nlines - 3))
            @test maximum(length, lines) <= 80
            io = IOBuffer()
            show(IOContext(io, :limit => false), sd)
            lines = readlines(seek(io, 0))
            @test length(lines) == nrow + 2
        end
    end
end


text:
using Distributions, Random, LinearAlgebra
# This is a bad design change it
Base.@kwdef struct EKI{A,B,C,D,E}
    Γ::A
    Σ::B
    N::C
    Δt::D
    J::E
end

⊗(a, b) = a * b'

"""
eki_step!(u, forward_map, y̅, J, ξ, Γ, Δt)

# Description
Take a single step of the Ensemble Kalman Inversion (EKI) algorithm
# Langevin formulation
du = Cᵘᵖ Γ⁻¹(y- G(u)) dt + Cᵘᵖ Γ^(-1/2) dW
# Fokker-Planck formulation
∂ₜρ = - ∇ᵤ⋅([y- G(u)]ᵀΓ⁻¹Cov(u, G(u)) + 1/2 Tr(Cov(u, G)Γ⁻¹Cov(u, G) Hᵤ(ρ) )
where the covariances are calculated using the instantaneous density
and Hᵤ is the Hessian of ρ. Note the location of the diffusivity
# Heuristically it is not ∇⋅(κ∇ρ) but rather Trace(κ∇⊗∇ρ)

# Arguments
`u`: Vector{S} where S. OVERWRITTEN. Empirical prior distribution
`forward_map`: Function. Operates on objects of type S
`data`: Vector or Array
`J`: int, number of ensemble members
`ξ`: multivariate distribution,
`Γ`: matrix corresponding to regularization
`Δt`: timestep size

# Return
nothing
"""
function eki_step!(u, forward_map, data, J, ξ, Γ, Δt)
    u̅ = mean(u)
    G = forward_map.(u) # error handling needs to go here
    G̅ = mean(G)

    Cᵘᵖ = (u[1] - u̅) ⊗ (G[1] - G̅)
    Cᵖᵖ = (G[1] - G̅) ⊗ (G[1] - G̅)
    for j in 2:J
        Cᵘᵖ += (u[j] - u̅) ⊗ (G[j] - G̅)
        Cᵖᵖ += (G[j] - G̅) ⊗ (G[j] - G̅)
    end
    Cᵘᵖ *= 1 / J
    Cᵖᵖ *= 1 / J

    y = [data + rand(ξ) for j in 1:J]
    residual = y - G

    factored_Cᵖᵖ = cholesky(Symmetric(Cᵖᵖ + 1 / Δt * Γ))
    for j in 1:J
        u[j] += Cᵘᵖ * (factored_Cᵖᵖ \ residual[j])
    end
    return nothing
end

function eki_loop!(u, forward_map, y̅, J, Γ, Δt, N)
    ξ = MvNormal(1 / Δt * Γ)
    for n in 1:N
        eki_step!(u, forward_map, y̅, J, ξ, Γ, Δt)
    end
end

function calibrate!(x, G, y, method::EKI)
    (; J, Γ, Δt, N) = method
    eki_loop!(G, x, y, J, Γ, Δt, N)
end

function two_player_eki_step!(u, forward_map, data, J, ξ, Γ, Δt)
    u̅ = mean(u)
    G = forward_map.(u) # error handling needs to go here
    G̅ = mean(G)

    Cᵘᵘ = (u[1] - u̅) ⊗ (u[1] - u̅)
    Cᵘᵖ = (u[1] - u̅) ⊗ (G[1] - G̅)
    Cᵖᵖ = (G[1] - G̅) ⊗ (G[1] - G̅)
    for j in 2:J
        Cᵘᵘ += (u[j] - u̅) ⊗ (u[j] - u̅)
        Cᵘᵖ += (u[j] - u̅) ⊗ (G[j] - G̅)
        Cᵖᵖ += (G[j] - G̅) ⊗ (G[j] - G̅)
    end
    Cᵘᵘ *= 1 / J
    Cᵘᵖ *= 1 / J
    Cᵖᵖ *= 1 / J

    y = [data + rand(ξ) for j in 1:J]
    residual = y - G

    factored_Cᵖᵖ = cholesky(Symmetric(Cᵖᵖ + 1 / Δt * Γ))
    for j in 1:J
        u[j] += Cᵘᵖ * (factored_Cᵖᵖ \ residual[j])
    end

    for j in div(J, 2):J
        u[j] += logpdf(MvNormal(u̅, Cᵘᵘ), u[j])
    end

    return nothing
end

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
module LDBDatasets

export textures_label,
       get_textures_dataset,
       get_mnist_dataset,
       get_1d_dataset,
       get_dataset

import HTTP, Tar
import CodecZlib: GzipDecompressorStream
import BufferedStreams: BufferedInputStream
import FileIO: load
import StatsBase: sample
import WaveletsExt: ClassData, generateclassdata
using MLDatasets

include("textures_dict.jl")
const textures_url = "https://sipi.usc.edu/database/textures.tar.gz"

"""
    extract_tgz_from_url(link, [dir])

Decompress and extract files from a .tar.gz file that is obtained from `link` into the
directory `dir`.

# Argummets
- `link::String`: URL link to download .tar.gz file.
- `dir::String`: (Default: `./data/`) Directory to save decompressed and extracted files.
"""
function extract_tgz_from_url(link::String, dir::String = "./data/")
    r = HTTP.request(:GET, link)
    Tar.extract(GzipDecompressorStream(BufferedInputStream(r.body)), dir)
end

"""
    get_textures_dataset(list, train_size, test_size, [dir])

Obtains the textures dataset.

# Arguments
- `list::AbstractVector{T} where T<:AbstractString`: List of texture files to generate
  classes of data. See the full list of texture files in the `textures_label` variable in
  `textures_dict.jl`.
- `train_size::Union{S, AbstractVector{S}} where S<:Integer`: Train size. If `train_size` is
  given as a scalar, then each class in `list` will have `train_size` number of data.
  Similarly, if `train_size` is a vector, then each entry is the train size to the
  corresponding entry in `list`.
- `test_size::Union{S, AbstractVector{S}} where S<:Integer`: Test size. If `test_size` is
  given as a scalar, then each class in `list` will have `test_size` number of data.
  Similarly, if `test_size` is a vector, then each entry is the test size to the
  corresponding entry in `list`.
- `dir::T where T<:AbstractString`: (Default: `"./data/textures/`) Directory to save
  decompressed and extracted textures files.

# Returns
Output in the form of `(train_x, train_y), (test_x, test_y)`
"""
function get_textures_dataset(list::AbstractVector{T},
                              train_size::AbstractVector{S},
                              test_size::AbstractVector{S},
                              dir::T = "./data/textures/") where
                             {T<:AbstractString, S<:Integer}
    @assert length(list) == length(train_size) == length(test_size)
    # Downloads the textures dataset
    if ~isdir(dir) || "textures" ∉ readdir(dir)
        extract_tgz_from_url(textures_url, dir)
    end
    # Get train and test samples
    train_x = []
    train_y = []
    test_x = []
    test_y = []
    for (texture, train_sz, test_sz) in zip(list, train_size, test_size)
        filename = textures_label[texture]
        path = dir * "textures/" * filename
        (tr_x, tr_y), (te_x, te_y) = generate_texture_samples(path, texture, train_sz, test_sz)
        push!(train_x, tr_x)
        push!(train_y, tr_y)
        push!(test_x, te_x)
        push!(test_y, te_y)
    end
    train_x = cat(train_x..., dims = 3)
    train_y = vcat(train_y...)
    test_x = cat(test_x..., dims = 3)
    test_y = vcat(test_y...)
    return (train_x, train_y), (test_x, test_y)
end

function get_textures_dataset(list::AbstractVector{T}, train_size::S, test_size::S,
    dir::T = "./data/textures/") where
    {T<:AbstractString, S<:Integer}
    # Each class in list will have `train_size` number of training data, and `test_size`
    # number of test data
    n = length(list)
    train_size = repeat([train_size], n)
    test_size = repeat([test_size], n)
    return get_textures_dataset(list, train_size, test_size, dir)
end

"""
    generate_texture_samples(path, texture, train_size, test_size)

Generate samples of size 128 x 128 for class `texture`.

# Arguments
- `path::T where T<:AbstractString`: Path of the `texture` file.
- `texture::T where T<:AbstractString`: The corresponding texture name.
- `train_size::S where S<:Integer`: Train size.
- `test_size::S where S<:Integer`: Test size.

# Returns
Output in the form of `(train_x, train_y), (test_x, test_y)`
"""
function generate_texture_samples(path::T, texture::T, train_size::S, test_size::S) where
                                 {T<:AbstractString, S<:Integer}
    # Load image and convert to elements to Float64
    img = load(path) |> x -> convert.(Float64, x)
    # Size of image
    m, n = size(img)
    # Get Cartesian Indices
    indices = CartesianIndices((1:(m-127), 1:(n-127)))
    # Randomly sample the cartesian indices
    samples = sample(indices, train_size+test_size, replace = false)
    # Generate subimage
    x = Array{Float64,3}(undef, (128,128,train_size+test_size))
    for (k,s) in enumerate(samples)
        i = s[1]
        j = s[2]
        x[:,:,k] = img[i:(i+127),j:(j+127)]
    end
    # Split to train and test sets
    train_x = x[:,:,1:train_size]
    train_y = repeat([texture], train_size)
    test_x = x[:,:,(train_size+1):end]
    test_y = repeat([texture], test_size)
    return (train_x, train_y), (test_x, test_y)
end

"""
    get_mnist_dataset([dir], [pad_widths])

Obtains the MNIST dataset.

# Arguments
- `train_size::Union{S, AbstractVector{S}} where S<:Integer`: (Default: 10000) Train size.
  If `train_size` is given as a scalar, then each class in MNIST will have at most
  `train_size` number of data. Similarly, if `train_size` is a vector, then each entry is
  the train size to the corresponding entry in 0,1,...,9.
- `test_size::Union{S, AbstractVector{S}} where S<:Integer`: (Default: 10000) Test size. If
  `test_size` is given as a scalar, then each class in MNIST will have at most `test_size`
  number of data. Similarly, if `test_size` is a vector, then each entry is the test size to
  the corresponding entry in 0,1,...,9.
- `pad_widths::Integer`: (Default: `2`) Width of zero padding. Used to reshape the size of
  images into that of power of 2 for full wavelet decomposition capabilities.
- `dir::T where T<:AbstractString`: (Default: `"./data/textures/`) Directory to save
  decompressed and extracted textures files.

# Returns
Output in the form of `(train_x, train_y), (test_x, test_y)`
"""
function get_mnist_dataset(train_size::AbstractVector{T} = 10000,
                           test_size::AbstractVector{T} = 10000,
                           pad_widths::T = 2,
                           dir::String = "./data/MNIST/") where T<:Integer
    # Downloads the MNIST dataset
    isdir(dir) || MNIST.download(dir, i_accept_the_terms_of_use = true)
    # Unpack the downloaded data into variables
    train_x, train_y = MNIST.traindata(Float64, dir = dir)
    test_x, test_y = MNIST.testdata(Float64, dir = dir)
    # Zero padding to make 28x28 images into 32x32 (using default settings)
    train_x = cat(train_x, dims = 3) |> x -> zero_padding(x, pad_widths)
    test_x = cat(test_x, dims = 3) |> x -> zero_padding(x, pad_widths)
    # Subsample dataset
    train_x, train_y = get_subsample(train_x, train_y, train_size)
    test_x, test_y = get_subsample(test_x, test_y, test_size)
    return (train_x, train_y), (test_x, test_y)
end

function get_mnist_dataset(train_size::T = 10000, test_size::T = 10000, pad_widths::T = 2,
                           dir::String = "./data/MNIST/") where T<:Integer
    return get_mnist_dataset(repeat([train_size], 10), repeat([test_size], 10), pad_widths, dir)
end

"""
    get_subsample(X, y, size)

Subsamples the data from MNIST.

# Arguments
- `X::AbstractArray{T} where T`: Data.
- `y::AbstractVector{S} where S`: Corresponding classes.
- `size::Vector{R} where R<:Integer`: Subsample size for each category.
"""
function get_subsample(X::AbstractArray{T}, y::AbstractVector{S}, size::Vector{R}) where
                      {T, S, R<:Integer}
    # Allocate variable to store list of subsamples
    Xsub = []
    ysub = []
    # Generate subsamples
    classes = unique(y)
    @assert length(size) == length(classes)
    for (i,c) in enumerate(classes)
        idx = findall(x -> (x==c), y)
        nsubsamples = min(length(idx), size[i])
        push!(Xsub, X[:,:,idx[1:nsubsamples]])
        push!(ysub, y[idx[1:nsubsamples]])
    end
    # Concatenate the subsamples
    Xsub = cat(Xsub..., dims = 3)
    ysub = vcat(ysub...)
    return Xsub, ysub
end

"""
    zero_padding(data, [pad_widths])

Zero-padding on `data` of size ``(m, n, k)`` into size ``(m + 2d, n + 2d, k)`` where ``d``
is the padding width.

# Arguments
- `data::AbstractArray{T,3} where T`: Array of size (m, n, k).
- `pad_widths::Integer`: (Default: `2`) Width of zero padding. Used to reshape the size of
  images into that of power of 2 for full wavelet decomposition capabilities.

# Returns
- `padded::Array{T,3}`: Padded `data`.
"""
function zero_padding(data::AbstractArray{T,3}, pad_widths::Integer = 2) where T
    m, n, k = size(data)
    padded = zeros(T, (m+2*pad_widths, n+2*pad_widths, k))
    padded[(pad_widths+1):(end-pad_widths), (pad_widths+1):(end-pad_widths), :] = data
    return padded
end

"""
    get_1d_dataset(type, train_size, test_size)

Wrapper for `WaveletsExt.generateclassdata` for generating 1D class data.

# Arguments
- `type::Symbol`: The type of class data, `:tri` (for triangular waveforms) or `:cbf` (for
  cylinder-bell-funnel).
- `train_size::Union{S, AbstractVector{S}} where S<:Integer`: Train size. If `train_size` is
  given as a scalar, then each class in `list` will have `train_size` number of data.
  Similarly, if `train_size` is a vector, then each entry is the train size to the
  corresponding entry in `list`.
- `test_size::Union{S, AbstractVector{S}} where S<:Integer`: Test size. If `test_size` is
  given as a scalar, then each class in `list` will have `test_size` number of data.
  Similarly, if `test_size` is a vector, then each entry is the test size to the
  corresponding entry in `list`.

# Returns
Output in the form of `(train_x, train_y), (test_x, test_y)`
"""
function get_1d_dataset(type::Symbol, train_size::T, test_size::T) where T<:Integer
    return get_1d_dataset(type, repeat([train_size], 3), repeat([test_size], 3))
end

function get_1d_dataset(type::Symbol,
                         train_size::AbstractVector{T},
                         test_size::AbstractVector{T}) where T<:Integer
    @assert length(train_size) == length(test_size) == 3
    train_x, train_y = generateclassdata(ClassData(type, train_size...))
    test_x, test_y = generateclassdata(ClassData(type, test_size...))
    return (train_x, train_y), (test_x, test_y)
end

get_tri_dataset(args...) = get_1d_dataset(:tri, args...)
get_cbf_dataset(args...) = get_1d_dataset(:cbf, args...)

"""
    get_dataset(dataset, args...)

Collect/generate dataset.

# Arguments
- `dataset::Symbol`: Type of class data. Supported types are:
    - `:mnist`: MNIST dataset (2D).
    - `:textures`: Textures dataset (2D).
    - `:tri`: Triangular waveforms (1D).
    - `:cbf`: Cylinder-bell-funnel dataset (1D).
- `args...`: Additional argument for generating dataset. See:
    - `:mnist`: [`get_mnist_dataset`](@ref).
    - `:textures`: [`get_textures_dataset`](@ref).
    - `:tri`: [`get_1d_dataset`](@ref).
    - `:cbf`: [`get_1d_dataset`](@ref).

# Returns
Output in the form of `(train_x, train_y), (test_x, test_y)`
"""
function get_dataset(dataset::Symbol, args...)
    @assert dataset ∈ [:mnist, :textures, :tri, :cbf]
    return @eval $(Symbol("get_$(dataset)_dataset"))($args...)
end

end # module


text:
using Plots
Plots.GRBackend()

function flood(f::Array;shift=(0.,0.),cfill=:RdBu_11,clims=(),kv...)
    if length(clims)==2
        @assert clims[1]<clims[2]
        @. f=min(clims[2],max(clims[1],f))
    else
        clims = (minimum(f),maximum(f))
    end
    Plots.contour(axes(f,1).+shift[1],axes(f,2).+shift[2],f',
        linewidth=0,linecolor=:black,
        fill = (true,palette(cfill)),
        clims = clims, aspect_ratio=:equal; kv...)
end

addbody(x,y;c=:black) = Plots.plot!(Shape(x,y), c=c, legend=false)
function μbody_plot!(flow::Flow{m};levels=[0.5],lines=:black) where {m}
    @inside flow.σ[I] = sum(WaterLily.ϕ(i,CartesianIndex(I,i),flow.μ₀) for i ∈ 1:m)/m
    BC!(flow.σ)
    contour!(flow.σ';levels,lines)
end

function sim_gif!(sim;duration=1,step=0.1,verbose=true,
                  remeasure=false,μbody=false,size=(700,600),kv...)
    t₀ = round(sim_time(sim))
    t = range(t₀,t₀+duration;step)
    gr(show=false;size)
    @time @gif for tᵢ in t
        sim_step!(sim,tᵢ;remeasure)
        @inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
        flood(sim.flow.σ;shift=(-0.5,-0.5),clims=(-5,5),kv...)
        μbody && μbody_plot!(sim.flow)
        verbose && println("tU/L=",round(tᵢ,digits=4),
            ", Δt=",round(sim.flow.Δt[end],digits=3))
    end
    return
end

function vectors(f::Array)
    x = axes(f,1)' .* ones(size(f,2))
    y = ones(size(f,2))' .* axes(f,2)
    u,v = f[:,:,1]',f[:,:,2]'
    Plots.quiver(x[:],y[:],quiver=(u[:],v[:]),
                aspect_ratio=:equal)
end


text:
module ch11
# Dictionaries
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#chap11

# -------------------------
# A Dictionary is a Mapping
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_a_dictionary_is_a_mapping
# -------------------------

eng2sp = Dict()
eng2sp["one"] = "uno"

# Create new Dict
eng2sp = Dict("one"=>"uno", "two"=>"dos", "three"=>"tres")
eng2sp["two"]

try
    eng2sp["four"]
catch e
    if e isa KeyError
        println("Missing key")
        println(e)
    else
        println(e)
    end
end

length(eng2sp)

# keys
keylist = keys(eng2sp)
for key in keylist
    println(key)
end

# put keys into an array
x = collect(keylist)

# IN operator
"one" ∈ keylist

# values
valuelist = values(eng2sp)
for value in values(eng2sp)
    println(value)
end

# --------------------------------------
# Dictionary as a Collection of Counters
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#dictionary_collection_counters
# --------------------------------------
function histogram(str)
    dic = Dict()
    for char in str
        if char ∉ keys(dic)
            dic[char] = 1
        else
            dic[char] += 1
        end
    end

    return dic
end

dic = histogram("audreyandryan")

dic['a']

# get to access Dict
get(dic, 'a', 0)
get(dic, 'z', 0)

# ------------------------
# Looping and Dictionaries
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_looping_and_dictionaries
# ------------------------
function printhist(dict)
    for key in keys(dict)
        println(key, " ", dict[key])
    end
end
printhist(dic)

# --------------
# Reverse Lookup
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_reverse_lookup
# --------------

# isequal
func = isequal(3)
func(4)
func(3)

findall(isequal(3), dic)

# -----------------------
# Dictionaries and Arrays
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_dictionaries_and_arrays
# -----------------------

# inverts a dictionary
function invertdict(dict)
    idict = Dict()
    for key in keys(dict)
        value = dict[key]
        if value ∈ keys(idict)
            push!(idict[value], key)
        else
            idict[value] = [key] # make it an array
        end
    end

    return idict
end

idict = invertdict(dic)
#idict[2]
@show(idict[2])


text:
isblock(x) = length(x) == 1 && x[1] isa Expr && x[1].head == :block
function flatten_expr!(x)
    isb = isblock(x)
    if isb
        x = MacroTools.striplines(x[1])
        filter!(z->z isa Symbol || z.head != :line, x.args)
        x = (x.args...,)
    end
    x
end
function build_expr(head::Symbol, args)
    ex = Expr(head)
    append!(ex.args, args)
    ex
end

"""
    get_variables(O) -> Vector{Union{Sym, Term}}

Returns the variables in the expression. Note that the returned variables are
not wrapped in the `Num` type.

# Examples
```julia
julia> @parameters t
(t,)

julia> @variables x y z(t)
(x, y, z(t))

julia> ex = x + y + sin(z)
(x + y) + sin(z(t))

julia> Symbolics.get_variables(ex)
3-element Vector{Any}:
 x
 y
 z(t)
```
"""
get_variables(e::Num, varlist=nothing) = get_variables(value(e), varlist)
get_variables!(vars, e, varlist=nothing) = vars

is_singleton(e::Term) = operation(e) isa Sym
is_singleton(e::Sym) = true
is_singleton(e) = false

get_variables!(vars, e::Number, varlist=nothing) = vars

function get_variables!(vars, e::Symbolic, varlist=nothing)
    if is_singleton(e)
        if isnothing(varlist) || any(isequal(e), varlist)
            push!(vars, e)
        end
    else
        foreach(x -> get_variables!(vars, x, varlist), arguments(e))
    end
    return (vars isa AbstractVector) ? unique!(vars) : vars
end

function get_variables!(vars, e::Equation, varlist=nothing)
  get_variables!(vars, e.rhs, varlist)
end

get_variables(e, varlist=nothing) = get_variables!([], e, varlist)

# Sym / Term --> Symbol
Base.Symbol(x::Union{Num,Symbolic}) = tosymbol(x)
tosymbol(x; kwargs...) = x
tosymbol(x::Sym; kwargs...) = nameof(x)
tosymbol(t::Num; kwargs...) = tosymbol(value(t); kwargs...)

"""
    diff2term(x::Term) -> Term
    diff2term(x) -> x

Convert a differential variable to a `Term`. Note that it only takes a `Term`
not a `Num`.

```julia
julia> Symbolics.diff2term(Symbolics.value(D(D(x))))
xˍtt(t)
```
"""
function diff2term(O)
    istree(O) || return O
    if is_derivative(O)
        (x, t, order) = flatten_differential(O)
        return lower_varname(x, t, order)
    end
    return Term{Real}(operation(O), diff2term.(arguments(O)))
end

function flatten_differential(O::Term)
    @assert is_derivative(O) "invalid differential: $O"
    is_derivative(arguments(O)[1]) || return (arguments(O)[1], operation(O).x, 1)
    (x, t, order) = flatten_differential(arguments(O)[1])
    isequal(t, operation(O).x) || throw(ArgumentError("non-matching differentials on lhs: $t, $(operation(O).x)"))
    return (x, t, order + 1)
end

"""
    tosymbol(x::Union{Num,Symbolic}; states=nothing, escape=true) -> Symbol

Convert `x` to a symbol. `states` are the states of a system, and `escape`
means if the target has escapes like `val"y⦗t⦘"`. If `escape` then it will only
output `y` instead of `y⦗t⦘`.

# Examples

```julia
julia> @parameters t; @variables z(t)
(z(t),)

julia> Symbolics.tosymbol(z)
Symbol("z⦗t⦘")
```
"""
function tosymbol(t::Term; states=nothing, escape=true)
    if operation(t) isa Sym
        if states !== nothing && !(t in states)
            return nameof(operation(t))
        end
        op = nameof(operation(t))
        args = arguments(t)
    elseif operation(t) isa Differential
        term = diff2term(t)
        op = Symbol(operation(term))
        args = arguments(term)
    else
        @goto err
    end

    return escape ? Symbol(op, "⦗", join(args, ", "), "⦘") : op
    @label err
    error("Cannot convert $t to a symbol")
end

function lower_varname(var::Symbolic, idv, order)
    order == 0 && return var
    name = string(nameof(operation(var)))
    underscore = 'ˍ'
    idx = findlast(underscore, name)
    append = string(idv)^order
    if idx === nothing
        newname = Symbol(name, underscore, append)
    else
        nidx = nextind(name, idx)
        newname = Symbol(name[1:idx], name[nidx:end], append)
    end
    return Sym{symtype(operation(var))}(newname)(arguments(var)[1])
end

function lower_varname(t::Symbolic, iv)
    var, order = var_from_nested_derivative(t)
    lower_varname(var, iv, order)
end
lower_varname(t::Sym, iv) = t

"""
    makesym(x::Union{Num,Symbolic}, kwargs...) -> Sym

`makesym` takes the same arguments as [`tosymbol`](@ref), but it converts a
`Term` in the form of `x(t)` to a `Sym` in the form of `x⦗t⦘`.

# Examples
```julia
julia> @parameters t; @variables x(t)
(x(t),)

julia> Symbolics.makesym(x)
x⦗t⦘
```
"""
makesym(t::Symbolic; kwargs...) = Sym{symtype(t)}(tosymbol(t; kwargs...))
makesym(t::Num; kwargs...) = makesym(value(t); kwargs...)

var_from_nested_derivative(x, i=0) = (missing, missing)
var_from_nested_derivative(x::Term,i=0) = operation(x) isa Differential ? var_from_nested_derivative(arguments(x)[1],i+1) : (x,i)
var_from_nested_derivative(x::Sym,i=0) = (x,i)


text:
module exportData
include("../../Alessio_Data/loadAlessioData.jl");
using .loadAlessioData
export getExpData

function getExpData(ind; Vind=1, tInds=1:125)

        # get reference data
        basepath = "../../Alessio_Data/2011-05-28_Rec78-103_Pace_Apex/"

        # print some diagnostic notices
        print("Loading data for $(ind)...\t")

        # load in some data
        expStates = loadFile(basepath, ind)

        print("Done.\n")

        # form the space and time masks
        spaceMasks = [spatialMask(state) for state in expStates]
        timeMasks  = [temporalMask(state) for state in expStates]
        mutualSpaceMask = mutualMask(spaceMasks)
        mutualTimeMask = mutualMask(timeMasks)

        # determine mutualSpaceMask centroid for square center
        center = spatialMaskCentroid(mutualSpaceMask)
        samplePoint = Int.(center./0.06)
        t, V = sampleData(expStates, samplePoint, mutualTimeMask, mutualSpaceMask)
        APD, DI, APA = loadAlessioData.analyzeVoltage(t, V[1]; threshold=0.25)
        ii = 1:min(length(APD),length(DI))
        BCL = sum(APD[ii] .+ DI[ii])/length(APD[ii])

        # select EPI for data to fit?
        data = collect(transpose(V[Vind][tInds]));

        return t[tInds], data, BCL
end

function getExpData(ind; Vind=1, xInds=64:64, tInds=1:125)

        # get reference data
        basepath = "../../Alessio_Data/2011-05-28_Rec78-103_Pace_Apex/"

        # print some diagnostic notices
        print("Loading data for $(ind)...\t")

        # load in some data
        expStates = loadFile(basepath, ind)

        print("Done.\n")

        # form the space and time masks
        spaceMasks = [spatialMask(state) for state in expStates]
        timeMasks  = [temporalMask(state) for state in expStates]
        mutualSpaceMask = mutualMask(spaceMasks)
        mutualTimeMask = mutualMask(timeMasks)

        # determine mutualSpaceMask centroid for square center
        samplePoints = [[x,x] for x in xInds]
        BCLs = []
        data = []
        ts   = []
        for samplePoint in samplePoints
        	print("Sampling at $(samplePoint) index...\t")
        	t, V = sampleData(expStates, samplePoint, mutualTimeMask, mutualSpaceMask)
        	print("Done\n")

        	push!(ts, t)

		APD, DI, APA = loadAlessioData.analyzeVoltage(t, V[1]; threshold=0.25)
		ii = 1:min(length(APD),length(DI))
		BCL = sum(APD[ii] .+ DI[ii])/length(APD[ii])
		push!(BCLs, BCL)

		# select EPI for data to fit?
		if isempty(data)
			data = collect(transpose(V[Vind][tInds]))
		else
			data = vcat(data, collect(transpose(V[Vind][tInds])))
		end
	end
	t = ts[1]
	BCL = sum(BCLs)/length(BCLs)

	return t[tInds], data, BCL
end

end


text:
q1 = nominal_configuration(biped)
v1 = zeros(biped.nq)

h = 0.01
T = 100

s = Simulator(biped, T, h=h)
status = simulate!(s, q1, v1)

vis = Visualizer()
open(vis)
visualize!(vis, s)