text:
#=
   https://www.allendowney.com/blog/2018/10/21/the-game-of-ur-problem/
   """
   Here’s a probability puzzle to ruin your week.

   In the Royal Game of Ur, players advance tokens along a track with 14 spaces.
   To determine how many spaces to advance, a player rolls 4 dice with 4 sides. Two corners
   on each die are marked; the other two are not. The total number of marked corners —
   which is 0, 1, 2, 3, or 4 — is the number of spaces to advance.

   For example, if the total on your first roll is 2, you could advance a token to space 2.
   If you roll a 3 on the next roll, you could advance the same token to space 5.

   Suppose you have a token on space 13. How many rolls did it take to get there?
   """

   See:
   https://www.allendowney.com/blog/lions-and-tigers-and-bears/

   Allen Downey's solution:
   http://nbviewer.jupyter.org/github/AllenDowney/ThinkBayes2/blob/master/solutions/game_of_ur_soln.ipynb?flush=true

   cf ~/blog/game_of_ur_problem.blog
      ~/webppl/game_of_ur_problem.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function game_of_ur_problem2()

    function roll1(a)
        ss = sum(a)
        if ss == 13
            return  a
        elseif ss > 13
            return []
        else
            t = rand(DiscreteUniform(0,4))
            return roll1(vcat(a,t))
        end
    end

    # Simpler
    function roll(s,len)
        if s == 13
            return  len
        elseif s > 13
            return 0
        else
            t = rand(DiscreteUniform(0,4))
            return roll(s+t,len+1)
        end
    end

    len ~ Dirac(roll(0,0))
    # Skip length 0 events
    true ~ Dirac(len > 0)

end

model = game_of_ur_problem2()

num_chns = 4

# chns = sample(model, Prior(), MCMCThreads(), 10_000, num_chns)

# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(20), 1_000)
chns = sample(model, SMC(1000), 40_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

println("Distribution of len")
show_var_dist_pct(chns, :len)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@testset "Real  Ylm" begin

##
import ACE
using JuLIP.Testing
using LinearAlgebra, StaticArrays, BenchmarkTools, Test, Printf
using ACE.SphericalHarmonics
using ACE.SphericalHarmonics: dspher_to_dcart, PseudoSpherical,
               cart2spher, spher2cart, RSHBasis, index_y
using JuLIP: evaluate, evaluate_d, evaluate_ed

verbose = false

##
@info("Testing consistency of Real and Complex SH; Condon-Shortley convention")

function test_r2c(L, cY, rY)
   cYt = similar(cY)
   for l = 0:L
      m = 0
      i = index_y(l, m)
      cYt[i] = rY[i]
      for m = 1:l
         i_p = index_y(l, m)
         i_m = index_y(l, -m)
         # test the expressions
         #  Y_l^m    =      1/√2 (Y_{lm} - i Y_{l,-m})
         #  Y_l^{-m} = (-1)^m/√2 (Y_{lm} + i Y_{l,-m})
         cYt[i_p] = (1/sqrt(2)) * (rY[i_p] - im * rY[i_m])
         cYt[i_m] = (-1)^m * (1/sqrt(2)) * (rY[i_p] + im * rY[i_m])
      end
   end
   return cY ≈ cYt
end

maxL = 20
cSH = SHBasis(maxL)
rSH = RSHBasis(maxL)

for nsamples = 1:30
   R = ACE.Random.rand_sphere()
   cY = evaluate(cSH, R)
   rY = evaluate(rSH, R)
   print_tf(@test test_r2c(maxL, cY, rY))
end
println()

##

@info("Test: check derivatives of real spherical harmonics")
for nsamples = 1:30
   R = @SVector rand(3)
   SH = RSHBasis(5)
   Y, dY = evaluate_ed(SH, R)
   DY = Matrix(transpose(hcat(dY...)))
   errs = []
   verbose && @printf("     h    | error \n")
   for p = 2:10
      h = 0.1^p
      DYh = similar(DY)
      Rh = Vector(R)
      for i = 1:3
         Rh[i] += h
         DYh[:, i] = (evaluate(SH, SVector(Rh...)) - Y) / h
         Rh[i] -= h
      end
      push!(errs, norm(DY - DYh, Inf))
      verbose && @printf(" %.2e | %.2e \n", h, errs[end])
   end
   success = (minimum(errs[2:end]) < 1e-3 * maximum(errs[1:3])) || (minimum(errs) < 1e-10)
   print_tf(@test success)
end
println()

end

ext: jl


text:
module TestSingleton

using InitialValues: Singleton
using Test

@testset begin
    s = Singleton(nothing)
    @test ndims(s) == 0
    @test ndims(typeof(s)) == 0
    @test size(s) == ()
    @test axes(s) == ()
    @test length(s) == 1
    @test eltype(s) === Nothing
    @test s[] === nothing
    @test s[CartesianIndex()] === nothing
    @test collect(s)::Vector{Nothing} == [nothing]
end

end  # module


text:
module CollegeStratBase

using Dates, DocStringExtensions, Formatting, PrettyTables
using FilesLH, LatexLH, StructLH

include("types.jl");
include("constants.jl");
include("helpers.jl");
include("school_groups.jl");
include("directories.jl");
include("notation.jl");
include("time_units.jl");
include("unit_conversions.jl");
include("display.jl");
include("testing.jl");

# Types
export Double, TimeInt, ncInt, TypeInt, CollInt, SchoolInt, GridInt

# School groups
export AbstractSchoolGroups, Schooling3
export ed_idx, ed_label, ed_labels, ed_symbol, ed_symbols, ed_suffixes, ed_suffix, n_school

# Directories
export relBaseDir, base_dir, julia_dir, notation_dir, notation_copy_dir, project_dir, paper_dir, test_dir
export computer_out_dir, computer_mat_dir, computer_log_dir, computer_json_dir, global_comparison_dir
export copy_file

# Notation
export SymTable, symbol_table, reload_symbol_table, lsymbol, ldescription, symbol_entry, copy_symbol_table_from_dropbox
export write_notation_preamble, notation_preamble_path, write_notation_summary
export RegrIntercept, GpaLabel

# Display
export format_number, format_dollars, format_vector, chain_strings
export current_time, show_text_table, show_matrix, fpath_to_show, calibrated_string
export settings_list, settings_table

# Unit conversions
export hours_per_week_to_mtu, hours_data_to_mtu, hours_mtu_to_data, validate_mtu, mtu_to_hours_per_week, per_year_to_per_week
export model_to_data_courses, data_to_model_courses
export dollars_data_to_model, dollars_model_to_data

# Helpers
export present_value, pv_factor, save_text_file, data_model_labels

# Debugging
export dbgLow, dbgMedium, dbgHigh

# Testing
export test_header, test_divider

end # module


text:
# Matern 1/2 isotropic covariance Function

"""
    Mat12Iso <: MaternISO

Isotropic Matern 1/2 kernel (covariance)
```math
k(x,x') = σ^2 \\exp(-|x-y|/ℓ)
```
with length scale ``ℓ`` and signal standard deviation ``σ``.
"""
mutable struct Mat12Iso{T} <: MaternIso where {T<:Real}
    "Length scale"
    ℓ::T
    "Signal variance"
    σ2::T
    "Priors for kernel parameters"
    priors::Array
end

"""
Matern 1/2 isotropic covariance Function

    Mat12Iso(ll::T, lσ::T)

# Arguments
  - `ll::Real`: length scale (given on log scale)
  - `lσ::Real`: signal standard deviation (given on log scale)
"""
Mat12Iso(ll::T, lσ::T) where T = Mat12Iso{T}(exp(ll), exp(2 * lσ), [])

function set_params!(mat::Mat12Iso, hyp::AbstractVector)
    length(hyp) == 2 || throw(ArgumentError("Matern 1/2 covariance function has two parameters, received $(length(hyp))."))
    mat.ℓ, mat.σ2 = exp(hyp[1]), exp(2 * hyp[2])
end

get_params(mat::Mat12Iso{T}) where T = T[log(mat.ℓ), log(mat.σ2) / 2]
get_param_names(mat::Mat12Iso) = [:ll, :lσ]
num_params(mat::Mat12Iso) = 2

cov(mat::Mat12Iso, r::Number) = mat.σ2 * exp(-r / mat.ℓ)

@inline dk_dll(mat::Mat12Iso, r::Real) = r / mat.ℓ * cov(mat, r)


text:
using Plugins
import Plugins
using Test

abstract type AbstractState1 end
Plugins.TemplateStyle(::Type{AbstractState1}) = Plugins.ImmutableStruct()
abstract type AbstractState2 end
abstract type AbstractState3 end

struct Fielder1 <: Plugin end
Plugins.register(Fielder1)
struct Fielder2 <: Plugin end
Plugins.register(Fielder2)

_numberinit() =  0
_numberinit(p) = p
_any() = nothing
_dict() = Dict()
_dict(p) = Dict(string(p) => p)

Plugins.customfield(plugin::Fielder1, ::Type{AbstractState1}) = Plugins.FieldSpec("field1_1", Int64, _numberinit)
Plugins.customfield(plugin::Fielder1, ::Type{AbstractState2}) = Plugins.FieldSpec("field1_2", Any, _any)
Plugins.customfield(plugin::Fielder2, ::Type{AbstractState1}) = Plugins.FieldSpec(:field2_1, Dict{String, Any}, _dict)
Plugins.customfield(plugin::Fielder2, ::Type{AbstractState2}) = Plugins.FieldSpec(:field2_2, Fielder2)

Plugins.customfield(plugin::Fielder2, ::Type{AbstractState3}) = Plugins.FieldSpec(:field2_2, "this_should_throw")

abstract type AppState end

struct CustomTemplate <: Plugins.TemplateStyle end
Plugins.TemplateStyle(::Type{AppState}) = CustomTemplate()

Plugins.typedef(::CustomTemplate, spec) = quote
    evaltest = true
    mutable struct TYPE_NAME <: $(spec.parent_type)
        $(Plugins.structfields(spec))
    end;
    $(Plugins.default_constructor(spec))
    TYPE_NAME
end

Plugins.customfield(plugin::Fielder1, ::Type{AppState}) = Plugins.FieldSpec("field1", Int64, _numberinit)
Plugins.customfield(plugin::Fielder2, ::Type{AppState}) = Plugins.FieldSpec("field2", Float64, _numberinit)

mutable struct CustomFieldsApp{TCustomState}
    state::TCustomState
    function CustomFieldsApp(plugins, hookfns, stateargs...)
        stack = PluginStack(plugins, hookfns)
        state_type = customtype(stack, :AppStateImpl, AppState)
        return new{state_type}(Base.invokelatest(state_type, stateargs...))
    end
end

function testinjected(app::CustomFieldsApp{TState}, stateargs...) where TState
    app.state = TState(stateargs...) # No need to invokelatest
end

abstract type ErrState end

struct ErrTemplate <: Plugins.TemplateStyle end
Plugins.TemplateStyle(::Type{ErrState}) = ErrTemplate()

Plugins.typedef(::ErrTemplate, spec) = quote
    evaltest = true
    dfgdfg
    TYPE_NAME
end

@testset "Plugins.jl custom fields" begin
    @test isnothing(Plugins.customfield(Fielder1(), AbstractArray)) == true

    stack = PluginStack([Fielder1, Fielder2])
    res1 = Plugins.customfields(stack, AbstractState1)
    @test res1.allok == true
    @test res1.results[1] == Plugins.FieldSpec(:field1_1, Int64, _numberinit)
    @test res1.results[2] == Plugins.FieldSpec(:field2_1, Dict{String, Any}, _dict)
    res2 = Plugins.customfields(stack, AbstractState2)
    @test res2.allok == true
    @test res2.results[1] == Plugins.FieldSpec(:field1_2, Any, _any)
    @test res2.results[2] == Plugins.FieldSpec(:field2_2, Fielder2)

    s1 = customtype(stack, :State1, AbstractState1; unique_name=false)
    @test startswith(string(s1), "State1") == true
    s1i = s1(0, Dict())
    @test s1i.field1_1 == 0
    @test_throws Exception s1i.field1_1 = 43
    @test s1i.field2_1 isa Dict

    s1i2 = s1()
    @test s1i2.field1_1 == 0
    @test_throws Exception s1i2.field1_1 = 43
    @test s1i2.field2_1 isa Dict
    @test length(s1i2.field2_1) == 0

    s1i3 = s1(42)
    @test s1i3.field1_1 == 42
    @test_throws Exception s1i3.field1_1 = 43
    @test s1i3.field2_1 isa Dict
    @test length(s1i3.field2_1) == 1
    @test s1i3.field2_1["42"] == 42

    s2 = customtype(stack, :State2, AbstractState2)
    s2i = s2(0, Fielder2())
    @test s2i.field1_2 == 0
    @test s2i.field2_2 isa Fielder2

    s2i2 = s2()
    @test isnothing(s2i2.field1_2)
    @test s2i.field2_2 isa Fielder2

    # Same params yield to same type
    s2_1 = customtype(stack, :State2, AbstractState2)
    @test s2_1 == s2
    @test length(string(s2_1)) > length(string(:State2))

    # different params lead to different types
    s2_2 = customtype(stack, :State2, AbstractState1)
    @test s2_2 != s2_1
    @test length(string(s2_2)) > length(string(:State2))

    @test_throws Exception customtype(stack, :State2, AbstractState3)

    app = CustomFieldsApp([Fielder1, Fielder2], [], 0, 0.0)
    @test app.state.field1 === 0
    @test app.state.field2 === 0.0
    app.state.field1 = 43
    @test app.state.field1 === 43
    @test_throws Exception app.state.field1 = "43"
    @test evaltest == true
    testinjected(app, 43, 43.0)
    @test app.state.field1 === 43

    @test s1().field1_1 == 0
    @test s1().field2_1 isa Dict{String, Any}

    errstack = PluginStack([])
    @test_throws Exception customtype(errstack, :ErrState2, ErrState)

    # Parametetric type
    st1 = customtype(stack, :StateT1, AbstractState1, [:T1])
    @test typeof(st1) === UnionAll
    st1i = st1{Int}(0, Dict())
    @test st1i.field1_1 == 0
    @test_throws Exception st1i.field1_1 = 43
    @test st1i.field2_1 isa Dict
end


text:
include("restrictedSequencingNode.jl")
include("restrictedDD.jl")
include("roughRelaxedBounding.jl")


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct LogQuantNormal <: LogNormalQuantDist
    mu::NestedFloat
    sigma::NestedFloat
    q::NestedFloat
end

function logquantnormal(mu::Float64, sigma::Float64, q::Float64)::Float64
    return exp(quantnormal(mu, sigma, q))
end

lognormalquant(x::LogQuantNormal, mu, sigma, q) = logquantnormal(mu, sigma, q)


text:
# returns the parameters for the limiting Toeplitz
function givenstail(t₀::Real,t₁::Real)
    @assert t₀^2-4t₁^2≥0
    s = (-t₀ + sqrt(t₀^2-4t₁^2))/(2t₁)
    l⁰ = (t₀ + sqrt(t₀^2-4t₁^2))/2
    # if s∞^2 > 1
    #     s∞ = (t₀ + sqrt(t₀^2-4t₁^2))/(2t₁)
    #     l0 = (t₀ - sqrt(t₀^2-4t₁^2))/2
    # end
    c = -sqrt(1-s^2)
    γ¹ = t₁*c
    γ⁰ = c*t₀ + s*γ¹
    l¹ = 2t₁  # = c*γ¹ - st₁
    l² = -t₁*s
    c,s,l⁰,l¹,l²,γ¹,γ⁰
end

ql(a,b,t₀,t₁) = ql!(copy(a),copy(b),t₀,t₁)

function ql!(a,b,t₀,t₁)
    if t₀^2 < 4t₁^2
        error("A QL decomposition only exists outside the essential spectrum")
    end
    # The Givens rotations coming from infinity (with parameters c∞ and s∞) leave us with the almost triangular
    # a[n-1]  b[n-1]   0    0    0
    # b[n-1]   a[n]   t₁    0    0
    #   0       γ¹      γ⁰    0    0
    #   0      l2     l1   l0    0
    #   0       0     l2   l1   l0

    if t₀ < 0
        # we want positive on the diagonal
        Q,L=ql(-a,-b,-t₀,-t₁)
        return -Q,L
    end

    # match dimensions in paper: a has n perts and b has n-1 perts
    n = max(length(a), length(b)+1)
    if n > length(a)
        m = length(a)
        a = resize!(a,n)
        a[m+1:end] .= t₀
    end

    if n > length(b)+1
        m = length(b)
        b = resize!(b,n-1)
        b[m+1:end] .= t₁
    end
    c∞,s∞,l⁰,l¹,l²,γ¹∞,γ⁰∞ = givenstail(t₀,t₁)
    # use recurrence for c. If we have a_0,…,a_N,t0,t0…, then
    # we only need c_-1,c_0,c_1,…,c_{N-1}.
    c=Array{eltype(c∞)}(undef,n)
    s=Array{eltype(c∞)}(undef,n-1)

    # ranges from 1 to N
    γ¹ = Array{eltype(c∞)}(undef,n-1)

    # ranges from 0 to N
    γ⁰ = Array{eltype(c∞)}(undef,n)
    nrms = Array{eltype(c∞)}(undef,n-1)

    γ⁰[n] = c∞*a[n] + s∞*γ¹∞  # k = N

    if n ≠ 1
        γ¹[n-1] = c∞*b[n-1]  # k = N

        k=n-1
        nrms[k] = 1/sqrt(γ⁰[k+1]^2+b[k]^2)
        c[k+1] = γ⁰[k+1]*nrms[k]  # K = N-1
        s[k] = -b[k]*nrms[k] # K = N-1

        @inbounds for k=n-2:-1:1
            γ¹[k] = c[k+2]*b[k]  # k = N-1
            γ⁰[k+1] = c[k+2]*a[k+1] + s[k+1]*γ¹[k+1]  # k = N
            nrms[k] = 1/sqrt(γ⁰[k+1]^2+b[k]^2)
            c[k+1] = γ⁰[k+1]*nrms[k]  # K = N-1
            s[k] = -b[k]*nrms[k] # K = N-1
        end

        γ⁰[1] = c[2]*a[1] + s[1]*γ¹[1]  # k = 0
    end

    c[1] = γ⁰[1] ≥ 0 ? 1 : -1  # k = -1, can't use sign because of 0 case

    Q = HessenbergUnitary(Val{'L'},true,c,s,c∞,s∞)

    L = BandedMatrix{eltype(c∞)}(undef, (n+1,n), (2,0))

    L[1,1] = abs(γ⁰[1]) - l⁰
    @views L[band(0)][2:end] .=  1 ./nrms .- l⁰
    @views L[band(-1)][1:end-1] .=  c[2:end].*γ¹ .- s.*a[1:end-1] .- l¹
    view(L,band(-1))[end] = c∞*γ¹∞ - s∞*a[end] - l¹
    @views L[band(-2)][1:end-1] .= (-).(s[2:end].*b[1:end-1]) .- l²

    if n ≠ 1
        view(L,band(-2))[end] = -s∞*b[end] - l²
    end

    Q, ToeplitzOperator([l¹,l²],[l⁰]) + FiniteOperator(L,ℓ⁰,ℓ⁰)
end

discreteeigs(J::SymTriPertToeplitz) =
    2*J.b*discreteeigs(0.5*(J.dv .- J.a)/J.b,0.5*J.ev/J.b) .+ J.a

connectioncoeffsoperator(J::SymTriPertToeplitz) =
    connectioncoeffsoperator(0.5*(J.dv .- J.a)/J.b,0.5*J.ev/J.b)

# Spectral map takes SequenceSpace to the space of the diagonal Fun
struct SpectralMap{CC,QQ,RS,T} <: Operator{T}
    n::Int  # number of extra eigenvalues
    C::CC
    Q::QQ
    rangespace::RS
end

SpectralMap(n::Int,C::Operator{T},Q::Operator{T},rs) where {T} =
    SpectralMap{typeof(C),typeof(Q),typeof(rs),T}(n,C,Q,rs)

domainspace(::SpectralMap) = SequenceSpace()
rangespace(S::SpectralMap) = S.rangespace

function ldiv_coefficients(S::SpectralMap,v::AbstractVector;opts...)
    # leave first entries n alone
    r = copy(v)
    r[S.n+1:end] .= ldiv_coefficients(S.C,v[S.n+1:end])
    ldiv_coefficients(S.Q,r;opts...)
end

function mul_coefficients(S::SpectralMap,v::AbstractVector;opts...)
    r = mul_coefficients(S.Q,v;opts...)
    # leave first entries n alone
    r[S.n+1:end] .= mul_coefficients(S.C,r[S.n+1:end])
    r
end

function getindex(S::SpectralMap,k::Integer,j::Integer)
    v = mul_coefficients(S,[zeros(j-1);1])
    k ≤ length(v) && return v[k]
    zero(eltype(S))
end

isbanded(S::SpectralMap) = true
function bandwidths(S::SpectralMap)
    bi = bandwidths(S.Q)
    bi[1],bi[2]+bandwidth(S.C,2)
end

# InvSpectralMap is the inverse of a Spectral map
struct InvSpectralMap{CC,QQ,DS,T} <: Operator{T}
    n::Int  # number of extra eigenvalues
    C::CC
    Q::QQ
    domainspace::DS
end

InvSpectralMap(n::Int,C::Operator{T},Q::Operator{T},ds) where {T} =
    InvSpectralMap{typeof(C),typeof(Q),typeof(ds),T}(n,C,Q,ds)

domainspace(S::InvSpectralMap) = S.domainspace
rangespace(::InvSpectralMap) = SequenceSpace()

function mul_coefficients(S::InvSpectralMap,v::AbstractVector;opts...)
    # leave first entries n alone
    r = copy(v)
    r[S.n+1:end] .= ldiv_coefficients(S.C,v[S.n+1:end])
    ldiv_coefficients(S.Q,r;opts...)
end

function ldiv_coefficients(S::InvSpectralMap,v::AbstractVector;opts...)
    r = mul_coefficients(S.Q,v;opts...)
    # leave first entries n alone
    r[S.n+1:end] .= mul_coefficients(S.C,r[S.n+1:end])
    r
end

function getindex(S::InvSpectralMap,k::Integer,j::Integer)
    v = mul_coefficients(S,[zeros(j-1);1])
    k ≤ length(v) && return v[k]
    zero(eltype(S))
end

isbanded(S::InvSpectralMap) = false
bandwidths(S::InvSpectralMap) = (∞,∞)

inv(S::SpectralMap) = InvSpectralMap(S.n,S.C,S.Q,S.rangespace)
inv(S::InvSpectralMap) = SpectralMap(S.n,S.C,S.Q,S.domainspace)

eigen(Jin::SymTriPertToeplitz) = eigfromguess(Jin,discreteeigs(Jin))
function eigfromguess(Jin::SymTriPertToeplitz,approxeigs)
    Qret=Array{HessenbergUnitary{'U',Float64}}(undef, 0)
    λapprox=sort(approxeigs)

    ctsspec = ApproxFun.Segment(Jin.a-2Jin.b,Jin.a+2Jin.b)

    J=Jin

    if length(λapprox) == 0
        C=connectioncoeffsoperator(J)

        D=Multiplication(Fun(identity,Ultraspherical(1,ctsspec)),Ultraspherical(1,ctsspec))

        U=SpaceOperator(C,ℓ⁰,domainspace(D))
        return D,U
    end

    λ=Vector{Float64}()

    tol=1E-14
    for k=1:length(λapprox)
        μ=λapprox[k]

        Q,L=ql(J-μ*I)
        push!(Qret,deflate(Q',k-1))
        J=L*Q+μ*I

         while abs(J[1,2]) > tol
             # μ=J[1,1] DO NOT DO THIS. IF MU IS NOT ACCURATE, J[1,1] CAN BE AN INVALID SHIFT (MW)
             Q,L=ql(J-μ*I)
             J=L*Q+μ*I
             push!(Qret,deflate(Q',k-1))
         end

        push!(λ,J[1,1])
        J=J[2:end,2:end]
    end

    if length(λ) == 1
         Q=Qret[1]

         D=Multiplication(Fun(identity,PointSpace(λ)⊕Ultraspherical(1,ctsspec)),PointSpace(λ)⊕Ultraspherical(1,ctsspec))

         U=SpectralMap(length(λ),connectioncoeffsoperator(J),Q,domainspace(D))
         return D,U
    else
        Q=BandedUnitary(reverse!(Qret))
        D=Multiplication(Fun(identity,PointSpace(λ)⊕Ultraspherical(1,ctsspec)),PointSpace(λ)⊕Ultraspherical(1,ctsspec))

        U=SpectralMap(length(λ),connectioncoeffsoperator(J),Q,domainspace(D))
        return D,U
    end
end


text:
using JuMP, Gurobi, Distributions, StatsFuns, StatsBase

const gurobi_env = Gurobi.Env()

const numVariables = 100
const numJCC = 100
const conRHS = (numVariables*1.0)^2

xi_mean = zeros(Float64,numVariables)
xi_covariance = 0.5*ones(Float64,numJCC,numJCC)

for j = 1:numVariables
	xi_mean[j] = j*1.0/numVariables
end

for i = 1:numJCC
	xi_covariance[i,i] = 1.0
end

xi_cov_chol_tmp = cholfact(xi_covariance)
xi_cov_chol = xi_cov_chol_tmp[:L]

xi_mean_mat = zeros(Float64,numJCC,numVariables)
for j = 1:numVariables
	xi_mean_mat[:,j] = xi_mean[j]*ones(numJCC)
end

# tailored implementation of multivariate normal distribution
function myMvNormal(numSamples::Int64)

	xi = zeros(Float64,numJCC,numVariables,numSamples)

	xi_tmp = rand(Normal(0.0,1.0),numJCC,numVariables,numSamples)
	for samp = 1:numSamples
		xi[:,:,samp] = xi_cov_chol*xi_tmp[:,:,samp] + xi_mean_mat
	end

	return xi
end

# generate random samples from multivariate normal distribution
function generateRandomSamples(numSamples::Int64)

	xi = myMvNormal(numSamples)

	return xi
end

srand(1234)

# generate samples to estimate actual risk level of solution
const numAnalyticalSamples = 100000
const xi_analytical = generateRandomSamples(numAnalyticalSamples)

srand()

# compute conservative bound on the probability that
# constraints aren't satisfied
function boundRiskLevel(x::Array{Float64},reliabilityLevel::Float64=1E-06)

	const numSamples = size(xi_analytical,3)
	simpleEst::Bool = false

	numViolated::Int64 = checkScenarioConstraints(x,xi_analytical)

	riskLevel::Float64 = Inf

	if(simpleEst)
		riskLevel = numViolated*1.0/numSamples
		return riskLevel
	else
		if(numViolated == 0)
			riskLevel = -log(reliabilityLevel)/numSamples
			return riskLevel
		end

		 # termination criterion for bisection
		gamma_tolerance::Float64 = 0.1/numSamples

		# use tail bounds on binomial distribution to estimate lower and upper bounds
		gamma_low::Float64 = (numViolated - sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_low < 0.0)
			gamma_low = 0.0
		end
		gamma_up::Float64 = (numViolated + sqrt(-numSamples*log(reliabilityLevel)/2.0))/numSamples
		if(gamma_up > 1.0)
			gamma_up = 1.0
		end

		gamma::Float64 = 0.0
		func_value::Float64 = 0.0

		# use bisection to solve the nonlinear equation
		while (gamma_up - gamma_low > gamma_tolerance)

			gamma = (gamma_low + gamma_up)/2.0
			func_value = -reliabilityLevel
			for i = 1:numViolated
				tmp_func = i*log(gamma) + (numSamples-i)*log(1-gamma) + lfact(numSamples) - lfact(numSamples-i) - lfact(i)
				func_value += exp(tmp_func)
				if(func_value > 0)
					break
				end
			end
			if(func_value > 0)
				gamma_low = (gamma_low + gamma_up)/2.0
			else
				gamma_up = (gamma_low + gamma_up)/2.0
			end

		end

		riskLevel = gamma_up
		return riskLevel
	end

end

# determines if all constraints are satisfied for a given single scenario
# returns one if scenario constraints are all satisfied, zero otherwise
function checkScenarioConstraints(x::Array{Float64},xi::Array{Float64})

	numSamples::Int64 = size(xi,3)

	numViolated::Int64 = 0
	for samp = 1:numSamples
		for i = 1:numJCC
			con::Float64 = evaluateConstraint(i,x,xi[i,:,samp])
			if(con > 0)
				numViolated += 1
				break
			end
		end
	end

	return numViolated
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables xi, and constraint scalings
function evaluateConstraint(index::Int64,x::Array{Float64},xi::Array{Float64})

	con::Float64 = norm(xi.*x)^2 - conRHS

	return con
end

NumSamples = 50

# solution file name
const baseDirName = "C:/Users/rkannan/Desktop/SA for CCP/Scenario approximation/experiments/normopt_noniid_scenapprox/"

for samp = 1:NumSamples

	const dirName = baseDirName * string(samp) * "/"

	const fileName = dirName * "solution.txt"

	const outputFile = dirName * "trueRiskLevels.txt"
	open(outputFile, "w") do f
	end

	const timeFile = dirName * "trueRiskTime.txt"
	open(timeFile, "w") do f
	end

	soln = readdlm(fileName)

	const numPoints = size(soln,1)

	for iter = 1:numPoints
		tic()
		trueRiskLevel::Float64 = boundRiskLevel(soln[iter,:])
		riskTime = toq()
		println("Risk level #",iter," : ",trueRiskLevel,"  time: ",riskTime)
		open(outputFile, "a") do f
			write(f,"$trueRiskLevel \n")
		end
		open(timeFile, "a") do f
			write(f,"$riskTime \n")
		end
	end

end

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
import AQFED.Random
import AQFED
import Random123
import Random: MersenneTwister

function simulation(r, n)
    spot = 1.0
    vol = 1.0
    strike = 1.0
    sum = 0.0
    sqrttte = 1.0
    @simd for i = 1:n
        #z = randn(r, Float64)
        u = rand(r,Float64)
    z = AQFED.Math.norminv(u) #faster than with AS241, esp when combined with exp call.
        sum += max(spot * exp(vol * z * sqrttte - 0.5 * (vol * sqrttte)^2) - strike, 0)
    end
    return sum / n
end

function simsum(r, n)
    u = Vector{Float64}(undef, n)
    rand!(r, u)
    return sum(u)
end

@testset "BlackSim" begin
    refValue = AQFED.Black.blackScholesFormula(true, 1.0, 1.0, 1.0, 1.0, 1.0)
    gens = [
        MersenneTwister(20130129),
        AQFED.Random.MersenneTwister64(UInt64(20130129)),
        AQFED.Random.Well1024a(UInt32(20130129)),
        AQFED.Random.Mixmax17(UInt64(20130129)),
        AQFED.Random.Mixmax240(UInt64(20130129)),
        AQFED.Random.MRG63k3a(),
        Random123.Philox4x(UInt64, (20130129, 20100921), 10),
        AQFED.Random.Chacha8SIMD(),
        AQFED.Random.Blabla8(),
    ]
    for gen in gens
        value = simulation(gen, 64 * 1024)
        println(typeof(gen), " ", value, " ", refValue, " ", value - refValue)
        @test isapprox(refValue, value, atol = 1e-2)
    end
end

@testset "MT64FirstNumbers" begin
    init = UInt64.([0x12345, 0x23456, 0x34567, 0x45678])
    gen = AQFED.Random.MersenneTwister64(init)
    println("1000 outputs of UInt64")
    for i = 1:1000
        ui = rand(gen, UInt64)
        # print(" " ,ui)
        # if (i % 5 == 4)
        #     println()
        # end
    end
    #numbers comming from reference C implementation
    referenceNumbers = [
        0.35252031,
        0.51052342,
        0.79771733,
        0.39300273,
        0.27216673,
        0.72151068,
        0.43144703,
        0.38522290,
        0.20270676,
        0.58227313,
        0.80812143,
        0.83767297,
        0.92401619,
        0.84065425,
        0.00852052,
        0.13975395,
        0.35250930,
        0.71196972,
        0.14627395,
        0.17775331,
        0.61046382,
        0.49623272,
        0.23292425,
        0.25038837,
        0.04380664,
        0.43275994,
        0.74540936,
        0.33830700,
        0.68832616,
        0.68744230,
        0.63626548,
        0.85932936,
        0.37089670,
        0.50756304,
        0.69925960,
        0.83481025,
        0.09053196,
        0.09523253,
        0.17783108,
        0.78027239,
        0.70071054,
        0.51879252,
        0.83027285,
        0.92895011,
        0.72144803,
        0.18868644,
        0.83655674,
        0.20358945,
        0.99852143,
        0.88340103,
        0.46729949,
        0.96993433,
        0.00162682,
        0.46829774,
        0.59080423,
        0.54921999,
        0.42516462,
        0.54952196,
        0.99534722,
        0.04473888,
        0.71139235,
        0.91881407,
        0.33781561,
        0.45746234,
        0.78292126,
        0.69206723,
        0.66175448,
        0.07091147,
        0.18179208,
        0.38168454,
        0.38819527,
        0.42452711,
        0.22732724,
        0.16191307,
        0.36842667,
        0.13060083,
        0.68833248,
        0.60498705,
        0.19195304,
        0.26628584,
        0.17030858,
        0.23892426,
        0.38430236,
        0.28034283,
        0.76069020,
        0.21560653,
        0.78101667,
        0.90847812,
        0.06467974,
        0.18487868,
        0.23570471,
        0.29475460,
        0.65563767,
        0.10066446,
        0.57272419,
        0.88731391,
        0.60650995,
        0.96346079,
        0.32940100,
        0.29977746,
        0.03798193,
        0.18026822,
        0.22402746,
        0.45480119,
        0.98114604,
        0.25800668,
        0.94362433,
        0.17901062,
        0.36019313,
        0.45933644,
        0.68309457,
        0.28175454,
        0.00774729,
        0.77054527,
        0.99723413,
        0.59807532,
        0.10294164,
        0.32429228,
        0.54928986,
        0.18410980,
        0.08441555,
        0.14230333,
        0.58892064,
        0.94030475,
        0.35378784,
        0.77584320,
        0.71222448,
        0.83565208,
        0.47309248,
        0.23810761,
        0.74408520,
        0.08891527,
        0.09729786,
        0.38377368,
        0.05092308,
        0.69065638,
        0.10449489,
        0.45050670,
        0.92209534,
        0.80083714,
        0.27902692,
        0.26897142,
        0.50650468,
        0.80111472,
        0.54590012,
        0.96406097,
        0.63779553,
        0.81054357,
        0.75369248,
        0.47473037,
        0.89100315,
        0.89395984,
        0.09985519,
        0.34087631,
        0.22293557,
        0.24375510,
        0.31764191,
        0.04076993,
        0.06160830,
        0.41333434,
        0.11883030,
        0.04548820,
        0.01008040,
        0.25336184,
        0.07325432,
        0.49860151,
        0.07148695,
        0.89483338,
        0.87054457,
        0.15116809,
        0.59650469,
        0.47487776,
        0.43490298,
        0.36684681,
        0.16470796,
        0.76865078,
        0.42920071,
        0.20545481,
        0.87615922,
        0.80332404,
        0.36462506,
        0.49571309,
        0.51904488,
        0.15534589,
        0.43719893,
        0.16562157,
        0.37290862,
        0.91842631,
        0.21310523,
        0.87849154,
        0.18532269,
        0.81713354,
        0.52182344,
        0.51845619,
        0.96261204,
        0.18758718,
        0.68897600,
        0.61484764,
        0.46752993,
        0.05865458,
        0.11614359,
        0.90386866,
        0.45781805,
        0.70649579,
        0.50917048,
        0.21210656,
        0.97818608,
        0.00788342,
        0.61375222,
        0.67366318,
        0.24197878,
        0.66177985,
        0.10463932,
        0.67390799,
        0.50025262,
        0.88332650,
        0.77966851,
        0.13403622,
        0.54357114,
        0.97664854,
        0.06540961,
        0.24013176,
        0.67234032,
        0.91347883,
        0.35486839,
        0.87207865,
        0.43036581,
        0.23652488,
        0.81238450,
        0.72058432,
        0.42239916,
        0.80265764,
        0.03552838,
        0.61939480,
        0.50972420,
        0.21053832,
        0.59952743,
        0.36821802,
        0.45659617,
        0.12529468,
        0.76941623,
        0.99878168,
        0.08602783,
        0.81825937,
        0.39350710,
        0.86090923,
        0.36090230,
        0.75628888,
        0.45036982,
        0.44602266,
        0.20595631,
        0.62241953,
        0.36777732,
        0.47523727,
        0.50248178,
        0.73570362,
        0.48237781,
        0.45590948,
        0.73580783,
        0.96403851,
        0.94586342,
        0.48819868,
        0.48102038,
        0.94618182,
        0.90279924,
        0.78396650,
        0.85182389,
        0.92149394,
        0.32679198,
        0.83554856,
        0.28320609,
        0.34598409,
        0.82090005,
        0.40177958,
        0.38888785,
        0.77873931,
        0.23297931,
        0.75329335,
        0.30770340,
        0.71417540,
        0.68939065,
        0.36577776,
        0.50784857,
        0.50928090,
        0.02552055,
        0.85999075,
        0.26692089,
        0.01402799,
        0.67550392,
        0.48305605,
        0.74608351,
        0.63408891,
        0.58904230,
        0.44337996,
        0.42174728,
        0.74041679,
        0.72719148,
        0.19801992,
        0.66263633,
        0.10381594,
        0.32818760,
        0.68369661,
        0.56076212,
        0.68681921,
        0.91616269,
        0.39836106,
        0.39685027,
        0.97507945,
        0.91010563,
        0.27447360,
        0.95538357,
        0.76758522,
        0.60091060,
        0.37734461,
        0.82948248,
        0.06598078,
        0.50147615,
        0.08417763,
        0.18910044,
        0.51661735,
        0.55011011,
        0.64888175,
        0.82986845,
        0.15126656,
        0.92649390,
        0.25494941,
        0.73275293,
        0.94184393,
        0.84755226,
        0.45921936,
        0.72934054,
        0.43722403,
        0.34305596,
        0.10827860,
        0.29026676,
        0.01935431,
        0.46668573,
        0.83247509,
        0.26349603,
        0.01938542,
        0.43222250,
        0.18109983,
        0.29337450,
        0.16721917,
        0.94751650,
        0.67795254,
        0.56666228,
        0.20699452,
        0.23247262,
        0.19138610,
        0.73495506,
        0.85893600,
        0.83411526,
        0.93689655,
        0.91804752,
        0.99352333,
        0.03207550,
        0.28386071,
        0.48029543,
        0.18736013,
        0.31736452,
        0.72542230,
        0.57530912,
        0.04229918,
        0.84798296,
        0.21886935,
        0.98655615,
        0.52243102,
        0.22611020,
        0.42975741,
        0.21726739,
        0.10912048,
        0.96684473,
        0.01092456,
        0.12461901,
        0.57989070,
        0.39848707,
        0.06330277,
        0.62826828,
        0.01159081,
        0.23157320,
        0.64690912,
        0.44876902,
        0.04463930,
        0.18933780,
        0.21284518,
        0.61363480,
        0.67144845,
        0.38625586,
        0.75719122,
        0.40361050,
        0.26708873,
        0.54534727,
        0.90174015,
        0.58654140,
        0.44885346,
        0.35505544,
        0.65317830,
        0.26074572,
        0.39472912,
        0.54366914,
        0.75020660,
        0.76113614,
        0.24595582,
        0.03941247,
        0.60356153,
        0.23615721,
        0.01603475,
        0.72432457,
        0.39837424,
        0.04195329,
        0.81561058,
        0.34208440,
        0.00513953,
        0.92826234,
        0.11410393,
        0.86692030,
        0.25238726,
        0.98258626,
        0.53353856,
        0.72269001,
        0.71850984,
        0.66829681,
        0.03540769,
        0.01676450,
        0.23557835,
        0.78758497,
        0.85969589,
        0.14673207,
        0.28013860,
        0.17796942,
        0.69924087,
        0.44663597,
        0.62112513,
        0.44079883,
        0.48995231,
        0.18411497,
        0.18440877,
        0.74016388,
        0.28845694,
        0.22969080,
        0.76851164,
        0.15551473,
        0.28980810,
        0.40906710,
        0.47619039,
        0.72611392,
        0.55802939,
        0.69365597,
        0.85736313,
        0.83343150,
        0.21324760,
        0.45327806,
        0.33053855,
        0.98198279,
        0.53279389,
        0.76877035,
        0.20548656,
        0.37065042,
        0.59026910,
        0.67418036,
        0.23585843,
        0.98156397,
        0.27849804,
        0.56198954,
        0.68752287,
        0.30073445,
        0.69348664,
        0.72515585,
        0.40629047,
        0.09320027,
        0.24334978,
        0.91407662,
        0.97226538,
        0.33904970,
        0.01717092,
        0.60155725,
        0.03001652,
        0.50979706,
        0.80531036,
        0.17450719,
        0.84984399,
        0.00498130,
        0.51636405,
        0.14080868,
        0.62289701,
        0.07853030,
        0.70567541,
        0.79844050,
        0.63766566,
        0.03559031,
        0.40994535,
        0.08423996,
        0.00389626,
        0.50608347,
        0.19622681,
        0.90537903,
        0.75458034,
        0.75102094,
        0.81491673,
        0.92925931,
        0.38074332,
        0.54817053,
        0.72593246,
        0.02146791,
        0.57990460,
        0.87921074,
        0.59913886,
        0.66726893,
        0.24269154,
        0.73344575,
        0.71826052,
        0.92313935,
        0.05212996,
        0.93771536,
        0.69489385,
        0.57581887,
        0.48106155,
        0.06808800,
        0.33633940,
        0.69142320,
        0.46566781,
        0.70654143,
        0.16541368,
        0.76257631,
        0.82777900,
        0.62958327,
        0.34757935,
        0.10891487,
        0.79912728,
        0.01156543,
        0.23111261,
        0.58535640,
        0.87461956,
        0.21723454,
        0.80409615,
        0.33169686,
        0.72800785,
        0.31218099,
        0.13729737,
        0.41637635,
        0.01234597,
        0.58313811,
        0.66746028,
        0.05105595,
        0.14930937,
        0.56044864,
        0.76196851,
        0.98800104,
        0.37075949,
        0.88740864,
        0.40697115,
        0.96598278,
        0.86013661,
        0.85386784,
        0.23986516,
        0.39027464,
        0.59593927,
        0.00161530,
        0.31768197,
        0.65702729,
        0.66461914,
        0.62937471,
        0.92120758,
        0.87578958,
        0.37539860,
        0.59182615,
        0.12092214,
        0.55130437,
        0.86365117,
        0.38725162,
        0.28757657,
        0.42803199,
        0.39014405,
        0.50253853,
        0.85306128,
        0.92018995,
        0.71421618,
        0.54236780,
        0.96221157,
        0.22956898,
        0.96519876,
        0.06694102,
        0.11915854,
        0.01354308,
        0.24720070,
        0.71671739,
        0.00604305,
        0.65012352,
        0.71151390,
        0.46616159,
        0.99228224,
        0.20684576,
        0.62941006,
        0.84535326,
        0.30678993,
        0.55264568,
        0.50094784,
        0.39409122,
        0.15479416,
        0.36536318,
        0.51925656,
        0.65567178,
        0.67255519,
        0.55089659,
        0.42194295,
        0.27172413,
        0.79540954,
        0.71594806,
        0.88372598,
        0.29179452,
        0.66411306,
        0.57064687,
        0.42494633,
        0.73389255,
        0.12097313,
        0.53338622,
        0.38493233,
        0.79348021,
        0.01851341,
        0.58594454,
        0.88396240,
        0.04410730,
        0.67419924,
        0.62770011,
        0.64644200,
        0.40335135,
        0.17952644,
        0.55564678,
        0.56643922,
        0.37715015,
        0.87092180,
        0.56726159,
        0.34011210,
        0.13661819,
        0.11474177,
        0.93930097,
        0.48549077,
        0.28484289,
        0.13374371,
        0.40966056,
        0.73662873,
        0.37355323,
        0.65216092,
        0.27372469,
        0.56032082,
        0.14882684,
        0.95462890,
        0.17090266,
        0.92374766,
        0.98368259,
        0.68448367,
        0.02872548,
        0.68598279,
        0.04601084,
        0.17170501,
        0.08906644,
        0.23730372,
        0.02929037,
        0.38566261,
        0.68957569,
        0.53021050,
        0.44200157,
        0.32085701,
        0.72520053,
        0.17454174,
        0.19676599,
        0.88243877,
        0.87030228,
        0.15124486,
        0.78670160,
        0.51731632,
        0.56674531,
        0.20910664,
        0.84962640,
        0.05220467,
        0.91783159,
        0.19138968,
        0.68126378,
        0.79574471,
        0.14910848,
        0.28030331,
        0.98067264,
        0.31263980,
        0.67448964,
        0.69266650,
        0.40033551,
        0.22789781,
        0.78317066,
        0.55815261,
        0.11247054,
        0.47337901,
        0.46310033,
        0.53192452,
        0.56164078,
        0.41750378,
        0.43880622,
        0.69739327,
        0.11092778,
        0.18333765,
        0.67222441,
        0.12789170,
        0.88316806,
        0.37891271,
        0.14935268,
        0.64522185,
        0.93902079,
        0.62481092,
        0.21794927,
        0.71535266,
        0.62169579,
        0.65147153,
        0.01411645,
        0.96413465,
        0.01021578,
        0.50605180,
        0.51595053,
        0.03308040,
        0.01497870,
        0.07809658,
        0.35743383,
        0.58079701,
        0.11785557,
        0.89568677,
        0.38793964,
        0.37117709,
        0.13994133,
        0.11032813,
        0.99998594,
        0.06695042,
        0.79774786,
        0.11093584,
        0.23879095,
        0.85918615,
        0.16109636,
        0.63479696,
        0.75023359,
        0.29061187,
        0.53764772,
        0.30652318,
        0.51387302,
        0.81620973,
        0.82433610,
        0.18302488,
        0.79048957,
        0.07598187,
        0.27887732,
        0.37061042,
        0.36441016,
        0.93736882,
        0.77480946,
        0.02269132,
        0.40309874,
        0.16427650,
        0.13969296,
        0.57605029,
        0.00242426,
        0.56626691,
        0.84390990,
        0.87455806,
        0.12321023,
        0.87561663,
        0.60431578,
        0.35880839,
        0.50426282,
        0.50697689,
        0.06631164,
        0.14976092,
        0.89356018,
        0.91473662,
        0.04235237,
        0.50073724,
        0.75969690,
        0.91743994,
        0.79352335,
        0.58078351,
        0.91819984,
        0.53520520,
        0.18267367,
        0.05608828,
        0.68315721,
        0.27264599,
        0.41245634,
        0.69706222,
        0.69666203,
        0.08967342,
        0.64081905,
        0.22576796,
        0.69315628,
        0.53981640,
        0.76059129,
        0.56712344,
        0.94318621,
        0.44081094,
        0.31699284,
        0.29477911,
        0.80069824,
        0.28366921,
        0.96718081,
        0.85345644,
        0.11681215,
        0.47600710,
        0.33448255,
        0.31217271,
        0.35469241,
        0.59511650,
        0.49583692,
        0.48922303,
        0.20215259,
        0.60159380,
        0.17882055,
        0.77601258,
        0.71020391,
        0.41833503,
        0.71522856,
        0.87534517,
        0.43703394,
        0.43056077,
        0.64828071,
        0.43069441,
        0.39356849,
        0.32063367,
        0.92788963,
        0.16878266,
        0.56762591,
        0.56042446,
        0.84958464,
        0.79408949,
        0.08220340,
        0.13922856,
        0.82529019,
        0.27134959,
        0.00278080,
        0.66192389,
        0.01782933,
        0.95404763,
        0.50787645,
        0.85320521,
        0.83690362,
        0.83771227,
        0.46268665,
        0.31716742,
        0.01716647,
        0.68264674,
        0.01789888,
        0.30446846,
        0.14942271,
        0.26982182,
        0.74933947,
        0.50394161,
        0.78444542,
        0.40009256,
        0.40333422,
        0.16627342,
        0.01898760,
        0.04221829,
        0.77960213,
        0.66230976,
        0.56015996,
        0.49535426,
        0.38536259,
        0.40406773,
        0.99930568,
        0.00857945,
        0.16158390,
        0.64805163,
        0.20237524,
        0.59106326,
        0.76968277,
        0.96887042,
        0.29264851,
        0.97373775,
        0.16767633,
        0.33014482,
        0.27426548,
        0.10947014,
        0.75920652,
        0.37757457,
        0.13125207,
        0.00826451,
        0.96684342,
        0.69362226,
        0.22763554,
        0.20717541,
        0.42112268,
        0.22803038,
        0.33481806,
        0.14968742,
        0.71598558,
        0.55126711,
        0.64518015,
        0.65170197,
        0.89103003,
        0.72728361,
        0.24485454,
        0.09410780,
        0.79818029,
        0.54212409,
        0.17790462,
        0.64442619,
        0.62193511,
        0.51193256,
        0.02848781,
        0.05719604,
        0.45795152,
        0.03219332,
        0.28310254,
        0.85746127,
        0.64890240,
        0.20658356,
        0.50946422,
        0.80432490,
        0.08354468,
        0.09222723,
        0.67455943,
        0.44638771,
        0.76366629,
        0.99677267,
        0.89311242,
        0.11627279,
        0.09181302,
        0.44767077,
        0.16448724,
        0.26005539,
        0.28670391,
        0.52465703,
        0.43598116,
        0.41869096,
        0.98043420,
        0.01497272,
        0.51791571,
        0.61825308,
        0.85503436,
        0.63025655,
        0.02719292,
        0.09865668,
        0.30321729,
        0.56998039,
        0.14946350,
        0.64823918,
        0.19931639,
        0.14623555,
        0.54169913,
        0.68944135,
        0.73551005,
        0.46743658,
        0.04109096,
        0.26625801,
        0.09537298,
        0.98207890,
        0.58109721,
        0.70793680,
        0.84379365,
        0.42774726,
        0.12653597,
        0.08566633,
        0.53366781,
        0.33960092,
        0.11036831,
        0.84464510,
        0.16493476,
        0.92493443,
        0.87640673,
        0.52727644,
        0.57181349,
        0.65071340,
        0.00978637,
        0.31700693,
        0.69148222,
        0.85063311,
        0.06781819,
        0.30794534,
        0.65541667,
        0.16400484,
        0.06886223,
        0.96227205,
        0.09633060,
        0.34513153,
        0.31013900,
        0.78165882,
        0.39583699,
        0.86327936,
        0.69269199,
        0.11016575,
        0.67358419,
        0.81775427,
        0.50052824,
        0.30068582,
        0.16606837,
        0.62243724,
        0.47863741,
        0.68796498,
        0.31526949,
        0.41180883,
        0.23022147,
        0.82342139,
        0.83003381,
        0.53571829,
        0.41081533,
        0.48600142,
    ]
    println("\n1000 outputs of Float64")
    for i = 1:1000
        z = rand(gen, Float64)
        # print(" " ,z)
        # if (i % 5 == 4)
        #     println()
        # end
        @test isapprox(referenceNumbers[i], z; atol = 1e-7)
    end
end


text:
#=
Nibble operations
=================

Copyright (c) 2017 Ben J. Ward & Luis Yanes

This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.
=#

"""
    nibble_capacity{T<:Unsigned}(::Type{T})

Returns the number of nibbles that an unsigned integer of type `T`
holds. This is essentially twice the size of the type (in bytes).
"""
@inline function nibble_capacity(::Type{T}) where {T<:Unsigned}
    return sizeof(T) * 2
end

@inline function bitpair_capacity(::Type{T}) where {T<:Unsigned}
    return sizeof(T) * 4
end

"""
    enumerate_nibbles{T<:Unsigned}(x::T)

Count the number of set bits in each nibble (aligned 4 bit segments) of an
unsigned integer `x`.

E.g. An input of:

0100 0010 0001 0110 1100 1110 1101 1111

Would result in:

0001 0001 0001 0010 0010 0011 0011 0100

This is used to identify different occurances of certain bit patterns.
"""
@inline function enumerate_nibbles(x::T) where {T<:Unsigned}
    x = x - ((x >>> 1) & repeatpattern(T, 0x55))
    return (x & repeatpattern(T, 0x33)) + ((x >>> 2) & repeatpattern(T, 0x33))
end


text:
##
function badvar1(x::Vector{Float64})
ex2 = 0.0
ex = 0.0
n = length(x)
for i=1:n
  ex2 = ex2 + x[i]^2
  ex = ex + x[i]
end
return 1.0/(n-1)*(ex2 - (ex)^2/n)
end
@show badvar1([0.0,1,2])
@show badvar1([0.0,1,2].+1.0e8)
##
using Random
Random.seed!(2)
x = randn(100)
@show badvar1(x)
@show badvar1(x.+1.0e8)


text:
using Documenter, ImageBinarization
makedocs(sitename="Documentation",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "Home" => "index.md",
        "Function Reference" => "reference.md"
    ]
)
deploydocs(repo = "github.com/JuliaImages/ImageBinarization.jl.git")


text:
# Molecular Dynamics Result Abstract Interface

"""
    MolecularDynamicsResult
Abstract type to be extended by all concrete structs representing the result of a molecular dynamics simulation.
"""
abstract type MolecularDynamicsResult end

# -----------------------------------------------------------------------------
# Simulation Configuration
# -----------------------------------------------------------------------------

"""
    get_time_range(result::MolecularDynamicsResult)::AbstractVector{Unitful.Time}

Extract the unit-anotated time range of the simulation from the simulation result.

An implementer of this API should implement a method of this function for their custom result type.
"""
get_time_range(result::MolecularDynamicsResult) = throw(UnimplementedError(:get_time_range, result))
Base.length(result::MolecularDynamicsResult) = length(get_time_range(result::MolecularDynamicsResult))

"""
    get_num_bodies(result::MolecularDynamicsResult)::Integer

Extract the number of bodies in the simulation from the simulation result.

An implementer of this API should implement a method of this function for their custom result type.
"""
get_num_bodies(result::MolecularDynamicsResult) = throw(UnimplementedError(:get_num_bodies, result))

"""
    get_bounding_box(result::MolecularDynamicsResult)::SVector{3, SVector{3, Unitful.Length}}

Extract the unit-anotated bounding box of the simulation from the simulation result in an `AtomBase`-compatible format.

An implementer of this API should implement a method of this function for their custom result type.
"""
get_bounding_box(result::MolecularDynamicsResult) = throw(UnimplementedError(:get_bounding_box, result))

"""
    get_boundary_conditions(result::MolecularDynamicsResult)::SVector{3, BoundaryCondition}

Extract the boundary conditions of the simulation from the simulation result in an `AtomBase`-compatible format.

An implementer of this API should implement a method of this function for their custom result type.
"""
get_boundary_conditions(result::MolecularDynamicsResult) = throw(UnimplementedError(:get_boundary_conditions, result))

"""
    reference_temperature(result::MolecularDynamicsResult)::Union{Unitful.Temperature,Missing}

Extract the unit-anotated reference temperature of the simulation from the simulation result.
If there is no thermostat with a reference temperature in this simulation, return missing.

An implementer of this API should implement a method of this function for their custom result type if it supports thermostats.
If not implmented, the default implemenation just returns missing.
"""
reference_temperature(result::MolecularDynamicsResult) = missing

# -----------------------------------------------------------------------------
# Time-Dependent System Data
# -----------------------------------------------------------------------------

"""
    get_time(result::MolecularDynamicsResult, t::Integer)::Unitful.Time

Extract the unit-anotated time of the simulation at a particular timestep from the simulation result.
Automatically defaults to the end of the simulation when `t` is not passed.

The default implementation extracts the time from the result of `get_time_range`.
"""
get_time(result::MolecularDynamicsResult, t::Integer) = get_time_range(result)[t]
get_time(result::MolecularDynamicsResult) = get_time(result, length(result))

"""
    get_positions(result::MolecularDynamicsResult, t::Integer)::AbstractVector{SVector{3, Unitful.Length}}

Extract the unit-anotated position of each particle in the system at a particular timestep from the simulation result in an `AtomBase`-compatible format.
The returned positions should be normalized such that they lie within the bounding box of the system, even if the system is periodic.
Automatically defaults to the end of the simulation when `t` is not passed.

An implementer of this API should implement a method of this function for their custom result type.
"""
get_positions(result::MolecularDynamicsResult, t::Integer) = throw(UnimplementedError(:get_positions, result))
get_positions(result::MolecularDynamicsResult) = get_positions(result, length(result))

"""
    get_velocities(result::MolecularDynamicsResult, t::Integer)::AbstractVector{SVector{3, Unitful.Velocity}}

Extract the unit-anotated velocity of each particle in the system at a particular timestep from the simulation result in an `AtomBase`-compatible format.
Automatically defaults to the end of the simulation when `t` is not passed.

An implementer of this API should implement a method of this function for their custom result type.
"""
get_velocities(result::MolecularDynamicsResult, t::Integer) = throw(UnimplementedError(:get_velocities, result))
get_velocities(result::MolecularDynamicsResult) = get_velocities(result, length(result))

"""
    get_particles(result::MolecularDynamicsResult, t::Integer)::AbstractVector{Atom}

Extract the position of each particle in the system at a particular timestep from the simulation result in an `AtomBase`-compatible format.
Automatically defaults to the end of the simulation when `t` is not passed.

An implementer of this API should implement a method of this function for their custom result type.
"""
get_particles(result::MolecularDynamicsResult, t::Integer) = throw(UnimplementedError(:get_particles, result))
get_particles(result::MolecularDynamicsResult) = get_particles(result, length(result))

"""
    get_system(result::MolecularDynamicsResult, t::Integer)::AbstractSystem{3}

Extract the underlying system at a particular timestep from the simulation result.
Automatically defaults to the end of the simulation when `t` is not passed.

The default implementation combines the results of `get_particles`, `get_bounding_box`, `get_boundary_conditions`, and `get_time` at the timestep to create a `FlexibleSystem` wrapped in a `DyanmicSystem`.
An implementer of this API could implement a method of this function for their custom result type if it supports a more efficient way to extract the system.
"""
get_system(result::MolecularDynamicsResult, t::Integer) = FlexibleSystem(get_particles(result, t), get_bounding_box(result), get_boundary_conditions(result))
get_system(result::MolecularDynamicsResult) = get_system(result, length(result))

# -----------------------------------------------------------------------------
# Simulation Analysis
# -----------------------------------------------------------------------------

"""
    temperature(result::MolecularDynamicsResult, t::Integer)::Unitful.Temperature

Extract the unit-anotated temperature of the simulation at a particular timestep from the simulation result.
Automatically defaults to the end of the simulation when `t` is not passed.

An implementer of this API should implement a method of this function for their custom result type.
"""
temperature(result::MolecularDynamicsResult, t::Integer) = throw(UnimplementedError(:temperature, result))
temperature(result::MolecularDynamicsResult) = temperature(result, length(result))

"""
    kinetic_energy(result::MolecularDynamicsResult, t::Integer)::Unitful.Energy

Extract the unit-anotated kinetic energy of the simulation at a particular timestep from the simulation result.
Automatically defaults to the end of the simulation when `t` is not passed.

An implementer of this API should implement a method of this function for their custom result type.
"""
kinetic_energy(result::MolecularDynamicsResult, t::Integer) = throw(UnimplementedError(:kinetic_energy, result))
kinetic_energy(result::MolecularDynamicsResult) = kinetic_energy(result, length(result))

"""
    potential_energy(result::MolecularDynamicsResult, t::Integer)::Unitful.Energy

Extract the unit-anotated potential energy of the simulation at a particular timestep from the simulation result.
Automatically defaults to the end of the simulation when `t` is not passed.

An implementer of this API should implement a method of this function for their custom result type.
"""
potential_energy(result::MolecularDynamicsResult, t::Integer) = throw(UnimplementedError(:potential_energy, result))
potential_energy(result::MolecularDynamicsResult) = potential_energy(result, length(result))

"""
    total_energy(result::MolecularDynamicsResult, t::Integer)::Unitful.Energy

Extract the unit-anotated total energy of the simulation at a particular timestep from the simulation result.
Automatically defaults to the end of the simulation when `t` is not passed.

The default implementation simply sums the results of `kinetic_energy` and `potential_energy` at the timestep.
An implementer of this API could implement a method of this function for their custom result type if it supports a more efficient way to calculate this quantity.
"""
total_energy(result::MolecularDynamicsResult, t::Integer) = kinetic_energy(result, t) + potential_energy(result, t)
total_energy(result::MolecularDynamicsResult) = total_energy(result, length(result))

"""
    rdf(result::MolecularDynamicsResult, start::Integer, stop::Integer)::Tuple{AbstractVector{Unitful.Length},AbstractVector{Real}}

Calculate the radial distribution function from the simulation result.

To include only a portion of the timesteps for faster computation, increase the start parameter or decrease the stop parameter; by default, the full time range will be iterated over.
The result is a named tuple of vectors:
    - r: unit-annotated interparticle radial distance bins
    - g: distribution value of each bin

The default implementation provided uses the provided implementations of `get_num_bodies`, `get_bounding_box`, `get_boundary_conditions`, and `get_positions`.
An implementor of the API could use a built in implementation if one that fits this spec is available.
"""
function rdf(result::MolecularDynamicsResult, start::Integer, stop::Integer)
    n = get_num_bodies(result)
    box = get_bounding_box(result)
    boundary_conditions = get_boundary_conditions(result)

    # TODO: support more boundary conditions
    @assert hcat(box...) == box[1][1] * I(3)
    @assert all(b isa Periodic for b ∈ boundary_conditions)

    d = Distances.PeriodicEuclidean([box[i][i] for i ∈ 1:3])
    L = box[1][1]

    maxbin = 1000
    z = zero(typeof(L))
    radius = 0.5 * L
    dr = radius / maxbin
    hist_threads = zeros(maxbin, nthreads())
    @threads for t ∈ start:stop
        pos = get_positions(result, t)
        for i ∈ 1:n
            for j ∈ i+1:n
                @inbounds r = evaluate(d, pos[i], pos[j])
                if z < r < radius
                    bin = ceil(Int, r / dr)
                    @inbounds hist_threads[bin, threadid()] += 2
                end
            end
        end
    end
    hist = sum(eachcol(hist_threads))

    c = 4 / 3 * π * n / L^3
    rs = zeros(typeof(dr), maxbin)
    gr = zeros(maxbin)
    tlen = length(start:stop)
    for bin ∈ 1:maxbin
        rlower = (bin - 1) * dr
        rupper = rlower + dr
        nideal = c * (rupper^3 - rlower^3)
        @inbounds rs[bin] = rlower + dr / 2
        @inbounds gr[bin] = (hist[bin] / (tlen * n)) / nideal
    end

    (; r=rs, g=gr)
end
rdf(result::MolecularDynamicsResult, start::Integer) = rdf(result, start, length(result))
rdf(result::MolecularDynamicsResult) = rdf(result, 1, length(result))


text:
@testset "Accuracy" begin

m = 5
σ = 2.0

@testset "High-level NFFT" begin
  @info "High-level NFFT"
  N = (32,32)
  eps =  1e-3
  D = length(N)

  M = prod(N)
  x = rand(Float64,D,M) .- 0.5

  fHat = rand(Float64,M) + rand(Float64,M)*im
  f = ndft_adjoint(x, N, fHat,)
  fApprox = nfft_adjoint(x, N, fHat, reltol=1e-9)
  e = norm(f[:] - fApprox[:]) / norm(f[:])
  @debug "error adjoint nfft "  e
  @test e < eps

  gHat = ndft(x, f)
  gHatApprox = nfft(x, f, reltol=1e-9)
  e = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
  @debug "error nfft "  e
  @test e < eps

  f = ndft_adjoint(x, N, ComplexF32.(fHat))
  fApprox = nfft_adjoint(x, N, ComplexF32.(fHat))
  e = norm(f[:] - fApprox[:]) / norm(f[:])
  @debug "error adjoint nfft "  e
  @test e < eps

  gHat = ndft(x, ComplexF32.(f))
  gHatApprox = nfft(x, ComplexF32.(f))
  e = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
  @debug "error nfft "  e
  @test e < eps
end

@testset "NFFT in multiple dimensions" begin
    for (u,N) in enumerate([(256,), (30,32), (10,12,14), (6,6,6,6)])
      for (pre, storeDeconv, blocking) in zip([ NFFT.LINEAR, NFFT.LINEAR, NFFT.LINEAR, NFFT.FULL, NFFT.TENSOR, NFFT.POLYNOMIAL, NFFT.POLYNOMIAL],
                                           [false, true, false, false, false, false, false],
                                           [true, true, false, false, true, true, false])
        eps = [1e-7, 1e-7, 1e-3, 1e-6, 1e-4]
        for (l,window) in enumerate([:kaiser_bessel, :cosh_type, :gauss, :kaiser_bessel_rev, :spline])
            D = length(N)
            @info "Testing $D, window=$window, pre=$pre, storeDeconv=$storeDeconv, block=$blocking"

            M = prod(N)
            x = rand(Float64,D,M) .- 0.5
            p = plan_nfft(x, N; m, σ, window, precompute = pre, storeDeconvolutionIdx = storeDeconv,
                         fftflags = FFTW.ESTIMATE, blocking)

            pNDFT = NDFTPlan(x, N)

            fHat = rand(Float64,M) + rand(Float64,M)*im
            f = adjoint(pNDFT) * fHat
            fApprox = adjoint(p) * fHat
            e = norm(f[:] - fApprox[:]) / norm(f[:])
            @debug "error adjoint nfft "  e
            @test e < eps[l]

            gHat = pNDFT * f
            gHatApprox = p * f
            e = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
            @debug "error nfft "  e
            @test e < eps[l]
        end
      end
    end
end

@testset "Abstract sampling points" begin
  @info "Abstract sampling points"
    M, N = rand(100:2:200, 2)
    x = range(-0.4, stop=0.4, length=M)
    p = plan_nfft(x, N, fftflags = FFTW.ESTIMATE)
end

@testset "Directional NFFT $D dim" for D in 2:3 begin
    # NFFT along a specified dimension should give the same result as
    # running a 1D NFFT on every slice along that dimension
        eps = 1e-4
        N = tuple( 2*rand(4:8,D)... )
        M = prod(N)
        for d in 1:D
            @info "Testing in $D dimensions directional NFFT along dim=$d"
            x = rand(M) .- 0.5

            f = rand(ComplexF64,N)
            p_dir = plan_nfft(x, N, dims=d)
            fHat_dir = p_dir * f
            g_dir = adjoint(p_dir) * fHat_dir

            p = plan_nfft(x, N[d])
            fHat = similar(fHat_dir)
            g = similar(g_dir)

            sz = size(fHat)
            Rpre = CartesianIndices( sz[1:d-1] )
            Rpost = CartesianIndices( sz[d+1:end] )
            for Ipost in Rpost, Ipre in Rpre
                idx = [Ipre, :, Ipost]
                fview = f[idx...]
                fHat[idx...] = p * vec(fview)

                fHat_view = fHat_dir[idx...]
                g[idx...] = adjoint(p) * vec(fHat_view)
            end

            e = norm( fHat_dir[:] - fHat[:] )
            @test e < eps

            e = norm( g_dir[:] - g[:] ) / norm(g[:])
            @test e < eps
        end
    end
end

@testset "Directional NFFT $D dim" for D in 3:3 begin
  # NFFT along  specified dimensions (1:2 or 2:3) should give the same result as
  # running a 2D NFFT on every slice along that dimensions
      eps = 1e-4
      N = tuple( 2*rand(4:8,D)... )
      M = prod(N)
      for d in 1:(D-1)
          dims = d:(d+1)
          @info "Testing in $D dimensions directional NFFT along dim=$(dims)"
          x = rand(2,M) .- 0.5

          f = rand(ComplexF64, N)
          p_dir = plan_nfft(x, N, dims=dims)
          fHat_dir = p_dir * f
          g_dir = adjoint(p_dir) * fHat_dir

          p = plan_nfft(x, N[dims])
          fHat = similar(fHat_dir)
          g = similar(g_dir)

          Rpre = CartesianIndices( N[1:dims[1]-1] )
          Rpost = CartesianIndices( N[(dims[end]+1):end] )
          for Ipost in Rpost, Ipre in Rpre
              idxf = [Ipre, :, :, Ipost]
              idxfhat = [Ipre, :, Ipost]

              fview = f[idxf...]
              fHat[idxfhat...] = p * fview

              fHat_view = fHat_dir[idxfhat...]
              g[idxf...] = adjoint(p) * fHat_view
          end

          e = norm( fHat_dir[:] - fHat[:] )
          @test e < eps

          e = norm( g_dir[:] - g[:] ) / norm(g[:])
          @test e < eps
      end
  end
end

end


text:
#=
Part of this work is derived from astropy/photutils and kbarbary/sep. The relevant derivations
are considered under a BSD 3-clause license. =#
using LazySets: AbstractHyperrectangle,
                EmptySet,
                HPolygon,
                Hyperrectangle,
                LinearMap,
                area,
                constraints_list,
                intersection,
                linear_map,
                matrix,
                set

using Rotations
using StaticArrays

function circular_overlap_single_subpixel(xmin, ymin, xmax, ymax, r, subpixels)
    frac = 0
    dx = (xmax - xmin) / subpixels
    dy = (ymax - ymin) / subpixels
    r2 = r^2

    x = xmin - 0.5dx
    for i in 1:subpixels
        x += dx
        y = ymin - 0.5dy
        for j in 1:subpixels
            y += dy
            if x^2 + y^2 < r2
                frac += 1
            end
        end
    end
    return frac / subpixels^2
end

"""Area of overlap between a rectangle and a circle"""
function circular_overlap_single_exact(xmin, ymin, xmax, ymax, r)
    R = float(typeof(xmin))
    r ≤ 0 && return zero(R)
    if 0 ≤ xmin
        0 ≤ ymin && return circular_overlap_core(xmin, ymin, xmax, ymax, r)
        0 ≥ ymax && return circular_overlap_core(-ymax, xmin, -ymin, xmax, r)
        return (circular_overlap_single_exact(xmin, ymin, xmax, 0, r) +
                circular_overlap_single_exact(xmin, 0, xmax, ymax, r))
    elseif 0 ≥ xmax
        0 ≤ ymin && return circular_overlap_core(-xmax, ymin, -xmin, ymax, r)
        0 ≥ ymax && return circular_overlap_core(-xmax, -ymax, -xmin, -ymin, r)
        return (circular_overlap_single_exact(xmin, ymin, xmax, 0, r) +
                circular_overlap_single_exact(xmin, 0, xmax, ymax, r))
    else
        0 ≤ ymin && return (circular_overlap_single_exact(xmin, ymin, 0, ymax, r) +
                            circular_overlap_single_exact(0, ymin, xmax, ymax, r))
        0 ≥ ymax && return (circular_overlap_single_exact(xmin, ymin, 0, ymax, r) +
                            circular_overlap_single_exact(0, ymin, xmax, ymax, r))
        return (circular_overlap_single_exact(xmin, ymin, 0, 0, r) +
                circular_overlap_single_exact(0, ymin, xmax, 0, r) +
                circular_overlap_single_exact(xmin, 0, 0, ymax, r) +
                circular_overlap_single_exact(0, 0, xmax, ymax, r))
    end
end

"""Core of circular overlap routine"""
function circular_overlap_core(xmin, ymin, xmax, ymax, r)
    R = float(typeof(xmin))
    xmin^2 + ymin^2 > r^2 && return zero(R)
    xmax^2 + ymax^2 < r^2 && return R((xmax - xmin) * (ymax - ymin))

    d1 = sqrt(xmax^2 + ymin^2)
    d2 = sqrt(xmin^2 + ymax^2)
    if d1 < r && d2 < r
        x1, y1 = sqrt(r^2 - ymax^2), ymax
        x2, y2 = xmax, sqrt(r^2 - xmax^2)
        area = ((xmax - xmin) * (ymax - ymin) -
                area_triangle(x1, y1, x2, y2, xmax, ymax) +
                area_arc(x1, y1, x2, y2, r))
    elseif d1 < r
        x1, y1 = xmin, sqrt(r^2 - xmin^2)
        x2, y2 = xmax, sqrt(r^2 - xmax^2)
        area = (area_arc(x1, y1, x2, y2, r) +
                area_triangle(x1, y1, x1, ymin, xmax, ymin) +
                area_triangle(x1, y1, x2, ymin, x2, y2))
    elseif d2 < r
        x1, y1 = sqrt(r^2 - ymin^2), ymin
        x2, y2 = sqrt(r^2 - ymax^2), ymax
        area = (area_arc(x1, y1, x2, y2, r) +
                area_triangle(x1, y1, xmin, y1, xmin, ymax) +
                area_triangle(x1, y1, xmin, y2, x2, y2))
    else
        x1, y1 = sqrt(r^2 - ymin^2), ymin
        x2, y2 = xmin, sqrt(r^2 - xmin^2)
        area = (area_arc(x1, y1, x2, y2, r) +
                area_triangle(x1, y1, x2, y2, xmin, ymin))
    end

    return R(area)
end

####################################
# Some geometric helpers

"""Area of a triangle defined by three vertices"""
area_triangle(x0, y0, x1, y1, x2, y2) = 0.5abs(x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1))

"""
Area of a circular segment above a chord between two points with circle radius `r`

[Reference](http://mathworld.wolfram.com/CircularSegment.html)
"""
function area_arc(x0, y0, x1, y1, r)
    a = sqrt((x1 - x0)^2 + (y1 - y0)^2)
    θ = 2 * asin(0.5 * a / r)
    return r^2 * (θ - sin(θ)) / 2
end

# is (x, y) contained within triangle specified by three points.
function inside_triangle(x, y, x1, y1, x2, y2, x3, y3)
    c = (y1 > y) ⊻ (y2 > y) && x < (x2 - x1) * (y - y1) / (y2 - y1) + x1
    c += (y2 > y) ⊻ (y3 > y) &&  x < (x3 - x2) * (y - y2) / (y3 - y2) + x2
    c += (y3 > y) ⊻ (y1 > y) && x < (x1 - x3) * (y - y3) / (y1 - y3) + x3

    # true if c is odd
    return c % 2 == 1
end

####################################
# Elliptical routines

"""
    inside_ellipse(x, y, h, k, cxx, cyy, cxy)

x: x coordinate of the test point
y: y coordinate of the test point
h: x coordinate of the center of ellipse
k: y coordinate of the center of ellipse
cxx, cyy, cxy: coefficients of equation of ellipse

Utility function to find whether a point is inside ellipse or not.

If point inside ellipse: Returns true else returns false

General equation of ellipse:
    cxx * (x - h)^2 + cxy * (x - h) * (y - k) + cyy * (y - k)^2 = 1
"""
@inline inside_ellipse(x, y, h, k, cxx, cyy, cxy) = cxx * (x - h)^2 + cxy * (x - h) * (y - k) + cyy * (y - k)^2  - 1 < 0

function elliptical_overlap_single_subpixel(xmin, ymin, xmax, ymax, cxx, cyy, cxy, subpixels)
    frac = 0
    dx = (xmax - xmin) / subpixels
    dy = (ymax - ymin) / subpixels

    x = xmin - 0.5dx
    for i in 1:subpixels
        x += dx
        y = ymin - 0.5dy
        for j in 1:subpixels
            y += dy
            if inside_ellipse(x, y, 0, 0, cxx, cyy, cxy)
                frac += 1
            end
        end
    end
    return frac / subpixels^2
end

function sort_order(d1, d2, d3)
    if d1 < d2
        if d2 < 3
            return SA[1, 2, 3]
        elseif d1 < d3
            return SA[1, 3, 2]
        else
            return SA[3, 1, 2]
        end
    else
        if d1 < d3
            return SA[2, 1, 3]
        elseif d2 < d3
            return SA[2, 3, 1]
        else
            return SA[3, 2, 1]
        end
    end
end

# overlap between triangle and unit circle
function triangle_unitcircle_overlap(x1, y1, x2, y2, x3, y3)
    # distances
    d1 = x1^2 + y1^2
    d2 = x2^2 + y2^2
    d3 = x3^2 + y3^2

    # order by distances
    ds = SA[d1, d2, d3]
    # order = sortperm(ds)
    order = sortperm(ds)
    ds = ds[order]
    x1, x2, x3 = SA[x1, x2, x3][order]
    y1, y2, y3 = SA[y1, y2, y3][order]

    # which are inside circle
    inside = map(d -> d < 1, ds)

    # on circle
    on = map(d -> d ≈ 1, ds)

    # triangle is completely inside circle
    if inside[3] || on[3]
        return area_triangle(x1, y1, x2, y2, x3, y3)

    # if vertex 1 or 2 are on the edge, then dot product with 3 to determine intersection
    elseif inside[2] || on[2]
        intersect13 = !on[1] || x1 * (x3 - x1) + y1 * (y3 - y1) < 0
        intersect23 = !on[2] || x2 * (x3 - x2) + y2 * (y3 - y2) < 0
        if intersect13 && intersect23 && !on[2]
            point1 = circle_segment_single2(x1, y1, x3, y3)
            point2 = circle_segment_single2(x2, y2, x3, y3)

            return (area_triangle(x1, y1, x2, y2, point1...) +
                    area_triangle(x2, y2, point1..., point2...) +
                    area_arc(point1..., point2..., 1))
        elseif intersect13
            point1 = circle_segment_single2(x1, y1, x3, y3)

            return (area_triangle(x1, y1, x2, y2, point1...) +
                    area_arc(x2, y2, point1..., 1))
        elseif intersect23
            point2 = circle_segment_single2(x2, y2, x3, y3)

            return (area_triangle(x1, y1, x2, y2, point2...) +
                    area_arc(x1, y1, point2..., 1))
        else
            return area_arc(x1, y1, x2, y2, 1)
        end
    elseif on[1]
        return 0.0
    elseif inside[1]
        point1, point2 = circle_segment(x2, y2, x3, y3)
        point3 = circle_segment_single2(x1, y1, x2, y2)
        point4 = circle_segment_single2(x1, y1, x3, y3)

        if point1[1] > 1 # no intersection
            # check if (x1, y2) and origin are on opposite sides of segment
            if (((0 - point3[2]) * (point4[1] - point3[1]) > (point4[2] - point3[2]) * (0 - point3[1])) !=
                ((y1 - point3[2]) * (point4[1] - point3[1]) > (point4[2] - point3[2]) * (x1 - point3[1])))
                return (area_triangle(x1, y1, point3..., point4...) + π - area_arc(point3..., point4..., 1))
            else
                return area_triangle(x1, y1, point3..., point4...) + area_arc(point3..., point4..., 1)
            end
        else
            # ensure that point1 is the point closest to (x2, y2)
            if (point2[1] - x2)^2 + (point2[2] - y2)^2 < (point1[1] - x2)^2 + (point1[2] - y2)^2
                point1, point2 = point2, point1
            end

            return (area_triangle(x1, y1, point3..., point1...) +
                    area_triangle(x1, y1, point1..., point2...) +
                    area_triangle(x1, y1, point2..., point4...) +
                    area_arc(point1..., point3..., 1) +
                    area_arc(point2..., point4..., 1))
        end
    else
        point1, point2 = circle_segment(x1, y1, x2, y2)
        point3, point4 = circle_segment(x2, y2, x3, y3)
        point5, point6 = circle_segment(x3, y3, x1, y1)

        if point1[1] ≤ 1
            xp = (point1[1] + point2[1]) / 2
            yp = (point1[2] + point2[2]) / 2
            return (triangle_unitcircle_overlap(x1, y1, x3, y3, xp, yp) +
                    triangle_unitcircle_overlap(x2, y2, x3, y3, xp, yp))
        elseif point3[1] ≤ 1
            xp = (point3[1] + point4[1]) / 2
            yp = (point3[2] + point4[2]) / 2
            return (triangle_unitcircle_overlap(x3, y3, x1, y1, xp, yp) +
                    triangle_unitcircle_overlap(x2, y2, x1, y1, xp, yp))
        elseif point5[1] ≤ 1
            xp = (point5[1] + point6[1]) / 2
            yp = (point5[2] + point6[2]) / 2
            return (triangle_unitcircle_overlap(x1, y1, x2, y2, xp, yp) +
                    triangle_unitcircle_overlap(x3, y3, x2, y2, xp, yp))
        else
            return inside_triangle(0, 0, x1, y1, x2, y2, x3, y3) ? π : 0.0
        end
    end
end

# intersection of a segment with the unit cirlce
function circle_segment(x1, y1, x2, y2)
    point1, point2 = circle_line(x1, y1, x2, y2)

    if ((point1[1] > x1 && point1[1] > x2) || (point1[1] < x1 && point1[1] < x2) ||
        (point1[2] > y1 && point1[2] > y2) || (point1[2] < y1 && point1[2] < y2))
        point1 = (2.0, 2.0)
    end

    if ((point2[1] > x1 && point2[1] > x2) || (point2[1] < x1 && point2[1] < x2) ||
        (point2[2] > y1 && point2[2] > y2) || (point2[2] < y1 && point2[2] < y2))
        point2 = (2.0, 2.0)
    end

    return point1[1] > 1 && point2[1] < 2 ?  (point1, point2) : (point2, point1)
end

# closest intersection of a line with the unit cirlce
function circle_segment_single2(x1, y1, x2, y2)
    point1, point2 = circle_line(x1, y1, x2, y2)
    dx1 = abs(point1[1] - x2)
    dy1 = abs(point1[2] - y2)
    dx2 = abs(point2[1] - x2)
    dy2 = abs(point2[2] - y2)

    if dx1 > dy1
        return dx1 > dx2 ? point2 : point1
    else
        return dy1 > dy2 ? point2 : point1
    end
end

# intersection of a line defined by two points with a unit cirlce
function circle_line(x1, y1, x2, y2)
    dx = x2 - x1
    dy = y2 - y1

    if dx ≈ 0 && dy ≈ 0
        return (2.0, 2.0), (2.0, 2.0)
    elseif abs(dx) > abs(dy)
        # find slope and intercept
        a = dy / dx
        b = y1 - a * x1

        # find quadratic determinant
        δ2 = 1 + a^2 - b^2
        if δ2 > 0 # real solutions
            δ = sqrt(δ2)
            _x1 = (-a * b - δ) / (1 + a^2)
            _y1 = a * _x1 + b
            _x2 = (-a * b + δ) / (1 + a^2)
            _y2 = a * _x2 + b
            return (_x1, _y1), (_x2, _y2)
        else
            return (2.0, 2.0), (2.0, 2.0)
        end
    else
        # find slope and intercept
        a = dx / dy
        b = x1 - a * y1

        # find quadratic determinant
        δ2 = 1 + a^2 - b^2
        if δ2 > 0 # real solutions
            δ = sqrt(δ2)
            _y1 = (-a * b - δ) / (1 + a^2)
            _x1 = a * _y1 + b
            _y2 = (-a * b + δ) / (1 + a^2)
            _x2 = a * _y2 + b
            return (_x1, _y1), (_x2, _y2)
        else
            return (2.0, 2.0), (2.0, 2.0)
        end
    end
end

function elliptical_overlap_exact(xmin, ymin, xmax, ymax, a, b, θ)
    sint, cost = sincosd(-θ)

    scale = a * b

    # reproject ellipse
    x1 = (xmin * cost - ymin * sint) / a
    y1 = (xmin * sint + ymin * cost) / b
    x2 = (xmax * cost - ymin * sint) / a
    y2 = (xmax * sint + ymin * cost) / b
    x3 = (xmax * cost - ymax * sint) / a
    y3 = (xmax * sint + ymax * cost) / b
    x4 = (xmin * cost - ymax * sint) / a
    y4 = (xmin * sint + ymax * cost) / b

    return scale * (triangle_unitcircle_overlap(x1, y1, x2, y2, x3, y3) +
                    triangle_unitcircle_overlap(x1, y1, x4, y4, x3, y3))
end

####################################
# Rectangular routines

function rectangular_overlap_single_subpixel(x0, y0, x1, y1, w, h, θ, subpixels)
    dx = (x1 - x0) / subpixels
    dy = (y1 - y0) / subpixels

    frac = 0
    x = x0 - 0.5dx
    for i in 1:subpixels
        x += dx
        y = y0 - 0.5dy
        for j in 1:subpixels
            y += dy
            if inside_rectangle(x, y, w, h, θ)
                frac += 1
            end
        end
    end

    return frac / subpixels^2

end

# see https://math.stackexchange.com/questions/69099/equation-of-a-rectangle
"""intersection with rectangular using implicit Lamé curve"""
function inside_rectangle(x, y, w, h, θ)
    # transform into frame of rectangle
    u, v = RotMatrix{2}(θ) \ SA[x, y]
    return abs(u) < w / 2 && abs(v) < h / 2
end

function rectangular_overlap_exact(xmin, ymin, xmax, ymax, w, h, θ)
    R = RotMatrix{2}(deg2rad(θ))
    aper = R * Hyperrectangle(zeros(SVector{2}), SA[w / 2, h / 2])
    dy = ymax - ymin
    dx = xmax - xmin
    pix = Hyperrectangle(SA[dx / 2 + xmin, dy / 2 + ymin], SA[dx / 2, dy / 2])
    return intersection_area(pix, aper)
end

# contrib: https://nbviewer.jupyter.org/github/mforets/escritoire/blob/master/2020/Week11/2D_intersection.ipynb
# the idea is to use the intersection for HPolygon, which is faster than the fallback because the
# normal vectors to the half-spaces are sorted in counter--clockwise fashion
function intersection_area(X::AbstractHyperrectangle{N},
    Y::LinearMap{N,<:AbstractHyperrectangle{N}}) where {N}
    X_clist = X |> constraints_list
    X_poly = HPolygon(X_clist, sort_constraints = true, prune = false, check_boundedness = false)

    Y_clist = linear_map(matrix(Y), set(Y)) |> constraints_list
    Y_poly = HPolygon(Y_clist, sort_constraints = true, prune = false, check_boundedness = false)
    inter = intersection(X_poly, Y_poly)

    return inter isa EmptySet ? 0.0 : area(inter)
end


text:
# ROCm-specific runtime libraries

## GPU runtime library

# reset the runtime cache from global scope, so that any change triggers recompilation
GPUCompiler.reset_runtime()

function signal_exception()
    flag_ptr = get_global_pointer(Val(:__global_exception_flag), Int64)
    unsafe_store!(flag_ptr, 1)
    # TODO: Can we sync here?
    endpgm() # stop this wavefront
    return
end

function report_exception(ex)
    # Add kernel ID and exception string to exception ring buffer
    ring_ptr = get_global_pointer(Val(:__global_exception_ring), LLVMPtr{ExceptionEntry,AS.Global})
    ring_ptr = unsafe_load(ring_ptr)
    our_signal = _completion_signal()
    prev = UInt64(1)
    while prev != UInt64(0)
        # Try to write to this slot, and skip if we fail (because another wavefront wrote first)
        prev = atomic_cas!(reinterpret(LLVMPtr{UInt64,AS.Global}, ring_ptr), UInt64(0), our_signal)
        if prev == UInt64(0)
            # We get a ReadOnlyMemoryError without this copy because the data is pinned to the device
            # TODO: Don't use an expensive host malloc
            ex_len = string_length(ex)::Csize_t
            copy_buf = reinterpret(LLVMPtr{UInt8,AS.Global}, malloc(ex_len))
            memcpy!(copy_buf, reinterpret(LLVMPtr{UInt8,AS.Global}, ex), ex_len)
            Base.unsafe_store!(reinterpret(LLVMPtr{UInt64,AS.Global}, ring_ptr+sizeof(UInt64)), UInt64(copy_buf))
            break
        elseif prev == UInt64(1)
            # Tail slot, give up
            break
        end
        ring_ptr += sizeof(ExceptionEntry)
    end
    return
end

report_oom(sz) = @rocprintf("ERROR: Out of dynamic GPU memory (trying to allocate %i bytes)\n", sz)

function device_string_to_host(ex)
    ex_ptr = reinterpret(LLVMPtr{UInt8,1}, ex)
    ex_len = string_length(ex_ptr)
    ex_str = reinterpret(LLVMPtr{UInt8,1}, malloc(Csize_t(ex_len+1)))
    if reinterpret(UInt64, ex_str) == 0
        @rocprintf("Device-to-host string conversion failed\n")
        return reinterpret(Cstring, 0)
    end
    memcpy!(ex_str, ex_ptr, ex_len)
    unsafe_store!(ex_str+ex_len, UInt8(0))
    reinterpret(Cstring, ex_str)
end

function report_exception_name(ex)
    # Pass argument in host buffer
    ex_str = device_string_to_host(ex)
    @rocprintf("""
        ERROR: a %s was thrown during kernel execution.
        Stacktrace:
        """, ex_str)
    free(reinterpret(Ptr{Cvoid}, ex_str))
    return
end

function report_exception_frame(idx, func, file, line)
    # Pass arguments in host buffers
    func_str = device_string_to_host(func)
    file_str = device_string_to_host(file)
    @rocprintf(" [%i] %s at %s:%i\n", idx, func_str, file_str, line)
    free(reinterpret(Ptr{Cvoid}, func_str))
    free(reinterpret(Ptr{Cvoid}, file_str))
    return
end

## ROCm device library

function load_device_libs(dev_isa, ctx)
    device_libs_path === nothing && return

    isa_short = replace(dev_isa, "gfx"=>"")
    device_libs = LLVM.Module[]
    bitcode_files = (
        "hc.amdgcn.bc",
        "hip.amdgcn.bc",
        "irif.amdgcn.bc",
        "ockl.amdgcn.bc",
        "oclc_isa_version_$isa_short.amdgcn.bc",
        "opencl.amdgcn.bc",
        "ocml.amdgcn.bc",
    )

    for file in bitcode_files
        ispath(joinpath(device_libs_path, file)) || continue
        name, ext = splitext(file)
        file_path = joinpath(device_libs_path, file)
        lib = parse(LLVM.Module, read(file_path), ctx)
        for f in LLVM.functions(lib)
            attrs = function_attributes(f)
            delete!(attrs, StringAttribute("target-features"))
        end
        push!(device_libs, lib)
    end

    @assert !isempty(device_libs) "No device libs detected!"
    return device_libs
end

function link_device_libs!(mod::LLVM.Module, dev_isa::String, undefined_fns)
    ufns = undefined_fns
    # TODO: only link if used
    # TODO: make these globally/locally configurable
    ctx = LLVM.context(mod)
    libs = load_device_libs(dev_isa, ctx)::Vector{LLVM.Module}
    link_oclc_defaults!(mod, dev_isa, ctx)
    for lib in libs
        # override libdevice's triple and datalayout to avoid warnings
        triple!(lib, triple(mod))
        datalayout!(lib, datalayout(mod))
    end

    GPUCompiler.link_library!(mod, libs)
    dispose.(libs)
end

function link_oclc_defaults!(mod::LLVM.Module, dev_isa::String, ctx; finite_only=false,
                             unsafe_math=false, correctly_rounded_sqrt=true, daz=false)
    # link in some defaults for OCLC knobs, to prevent undefined variable errors
    lib = LLVM.Module("OCLC", ctx)
    triple!(lib, triple(mod))
    datalayout!(lib, datalayout(mod))

    isa_short = replace(dev_isa, "gfx"=>"")
    for (name,value) in (
            "__oclc_ISA_version"=>parse(Int32, isa_short),
            "__oclc_finite_only_opt"=>Int32(finite_only),
            "__oclc_unsafe_math_opt"=>Int32(unsafe_math),
            "__oclc_correctly_rounded_sqrt32"=>Int32(correctly_rounded_sqrt),
            "__oclc_daz_opt"=>Int32(daz))
        gvtype = convert(LLVMType, typeof(value), ctx)
        gv = GlobalVariable(lib, gvtype, name, 4)
        init = ConstantInt(Int32(0), ctx)
        initializer!(gv, init)
        unnamed_addr!(gv, true)
        constant!(gv, true)
        linkage!(gv, LLVM.API.LLVMWeakAnyLinkage)
    end

    link!(mod, lib)
end


text:
function find_files_by_ext(startdir, ext = [""])
    foundfiles = String[]
    for (root, dirs, files) in walkdir(startdir)
        for file in files
            if last(splitext(file)) ∈ ext
                push!(foundfiles, joinpath(root, file))
            end
        end
    end
    return normpath.(foundfiles)
end

find_files_by_ext(startdir, ext::AbstractString) = find_files_by_ext(startdir, [ext])

function apply_citeproc(startdir)
    bib = first(find_files_by_ext(startdir, ".bib"))
    csl = first(find_files_by_ext(startdir, ".csl"))
    html = find_files_by_ext(startdir, ".html")
    for file in html
        run(`pandoc $file -o $file --bibliography $bib --csl $csl --metadata link-citations=true --filter pandoc-citeproc -f markdown --to html `)
    end
    return
end


text:
export backtracking_linesearch

function backtracking_linesearch(model_orig, q, dobs, f_prev, g, proj; alpha=1f0, tau=.1f0, c1=1f-4, maxiter=10, verbose=false)

    # evaluate FWI objective as function of step size alpha
    function objective(alpha,p)
        model.m = proj(model_orig.m + alpha*reshape(p,model.n))

        # Set up linear operator and calculate data residual
        info = JUDI.TimeModeling.Info(prod(model.n), dobs.nsrc, JUDI.TimeModeling.get_computational_nt(q.geometry,dobs.geometry,model))
        F = JUDI.TimeModeling.judiModeling(info,model,q.geometry,dobs.geometry)
        dpred = F*q
        return .5f0*norm(dpred - dobs)^2
    end

    model = deepcopy(model_orig)    # don't modify original model
    p = -g/norm(g,Inf)  # normalized descent direction
    f_new = objective(alpha,p)
    iter = 1
    verbose == true && println("    Iter LS: ", iter, "; ", f_new, " <= ", f_prev + c1*alpha*dot(g,p), "; alpha: ", alpha)

    # sufficient decrease (Armijo) condition
    while f_new > f_prev + c1*alpha*dot(g,p) && iter < maxiter
        alpha *= tau
        f_new = objective(alpha,p)
        iter += 1
        verbose == true && println("    Iter LS: ", iter, "; ", f_new, " <= ", f_prev + c1*alpha*dot(g,p), "; alpha: ", alpha)
    end
    return alpha*p
end

function backtracking_linesearch(model_orig, q, dobs, f_prev, g, proj, objective; alpha=1f0, tau=.1f0, c1=1f-4, maxiter=10, verbose=false)

    model = deepcopy(model_orig)    # don't modify original model
    p = -g/norm(g,Inf)  # normalized descent direction
    model.m = proj(model_orig.m + alpha*reshape(p,model_orig.n))
    f_new = objective(model, q, dobs; compute_gradient=false)
    println("    LS: ",f_new)
    iter = 1
    verbose == true && println("    Iter LS: ", iter, "; ", f_new, " <= ", f_prev + c1*alpha*dot(g,p), "; alpha: ", alpha)

    # sufficient decrease (Armijo) condition
    while f_new > f_prev + c1*alpha*dot(g,p) && iter < maxiter
        alpha *= tau
        model.m = proj(model_orig.m + alpha*reshape(p,model_orig.n))
        f_new = objective(model, q, dobs; compute_gradient=false)
        println("    LS: ",f_new)
        iter += 1
        verbose == true && println("    Iter LS: ", iter, "; ", f_new, " <= ", f_prev + c1*alpha*dot(g,p), "; alpha: ", alpha)
    end
    return alpha*p
end


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file include the SudokuTester module's source code, runs the tests, and writes a report.

include("sudoku_tester.jl")

using Test

@testset verbose=true "Sudoku puzzle" begin

@testset "Common definitions" begin
    SudokuTester.puzzle_2()
    SudokuTester.puzzle_rank_array()
    SudokuTester.puzzle_3_entry()
    SudokuTester.solvable_puzzle_2_construction()
    SudokuTester.solvable_puzzle_2_modification()
    SudokuTester.solvable_puzzle_2_assignment()
end

@testset "Permutations helpers" begin
    SudokuTester.block_permutations()
    SudokuTester.intrablock_permutations()
    SudokuTester.mirroring()
end

@testset "Composition helpers" begin
    SudokuTester.symbol_swap()
    SudokuTester.random_puzzle()
end

@testset "Compostition statistics" begin
    SudokuTester.random_permutations(100000)
    SudokuTester.random_permutations(200000)
    SudokuTester.random_permutations(400000)
end

@testset "Validity checks" begin
    SudokuTester.valid_subarray()
    SudokuTester.valid_puzzle_2()
    SudokuTester.satisfies_puzzle_2()
end

@testset "Solution helpers" begin
    SudokuTester.solve_example_puzzles()
    SudokuTester.solve_random_puzzles()
end

end # end top-level testset

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
module RecordStores

export recordstore

import BSON
import ZipFile
using BangBang: push!!
using Base: open_flags

addfile(writer::ZipFile.Writer, name) = ZipFile.addfile(writer, name)
files(x) = x.files
filename(x) = x.name

struct DirWriter
    path::String

    function DirWriter(path)
        mkdir(path)
        return new(path)
    end
end

addfile(writer::DirWriter, name) = joinpath(writer.path, name)

Base.close(::DirWriter) = nothing
Base.flush(::DirWriter) = nothing

struct DirReader
    path::String
end

files(reader::DirReader) = joinpath.(reader.path, readdir(reader.path))
filename(fullpath::AbstractString) = basename(fullpath)
Base.close(::DirReader) = nothing

const _archivers = (
    zip = (writer = ZipFile.Writer, reader = ZipFile.Reader),
    dir = (writer = DirWriter, reader = DirReader),
)

function guess_archiver(path)
    if endswith(path, ".zip")
        return :zip
    elseif endswith(path, Base.Filesystem.path_separator) || isdir(path)
        return :dir
    end
    error(
        "Cannot guess `archiver` for path `$path`.",
        " Please specify `archiver` keyword argument.",
    )
end

"""
    recordstore(path; archiver) -> store

# Keyword Arguments
- `archiver ∈ (:zip, :dir)`
"""
recordstore(path::AbstractString; archiver::Symbol = guess_archiver(path)) =
    RecordStore(path; pairs(_archivers[archiver])...)

struct RecordStore
    path::String
    writer
    reader
end

RecordStore(path; writer, reader) = RecordStore(path, writer, reader)

function Base.open(store::RecordStore; read=nothing, write=nothing)
    flags = open_flags(read=read, write=write)
    flags.write && return RecordWriter(store)
    flags.read && return RecordReader(store)
    throw(ArgumentError("Unsupported flags: $flags"))
end

function Base.open(store::RecordStore, mode::AbstractString)
    mode == "w" && return open(store, write=true)
    mode == "r" && return open(store, read=true)
    throw(ArgumentError("Unsupported mode: $mode"))
end

Base.read(store::RecordStore) = open(read, store)

struct RecordWriter
    path::String
    writer
    counter::typeof(Ref(0))
end

RecordWriter(path::AbstractString, writer) = RecordWriter(path, writer, Ref(0))
RecordWriter(store::RecordStore) = RecordWriter(store.path, store.writer(store.path))

Base.close(w::RecordWriter) = close(w.writer)
Base.flush(w::RecordWriter) = flush(w.writer)

maybeclose(_) = nothing
maybeclose(io::IO) = close(io)

function Base.write(w::RecordWriter, obj)
    i = w.counter[] += 1
    f = addfile(w.writer, "$i.bson")
    obj = deepcopy(obj)  # https://github.com/MikeInnes/BSON.jl/issues/26
    try
        BSON.bson(f, obj)
    finally
        maybeclose(f)  # flush
    end
    return w
end

struct RecordReader
    path::String
    reader
end

RecordReader(path::AbstractString, reader) = RecordReader(path, reader)
RecordReader(store::RecordStore) = RecordReader(store.path, store.reader(store.path))

Base.close(r::RecordReader) = close(r.reader)

fileid(x) = parse(Int, chop(filename(x); tail = length(".bson")))
sortedfiles(r::RecordReader) = sort!(collect(files(r.reader)); by = fileid)

function Base.read(r::RecordReader)
    return mapfoldl(BSON.load, push!!, sortedfiles(r); init = Union{}[])
end

Base.iterate(r::RecordReader) = iterate(r, (sortedfiles(r), 1))
function Base.iterate(r::RecordReader, (files, i))
    i > length(files) && return nothing
    return (BSON.load(files[i]), (files, i + 1))
end

Base.IteratorSize(::Type{<:RecordReader}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:RecordReader}) = Base.EltypeUnknown()

end # module


text:
using Test
using NCDatasets

filename = tempname()
# The mode "c" stands for creating a new file (clobber)
ds = NCDataset(filename,"c")

# define the dimension "lon" and "lat" with the size 10 and 11 resp.
ds.dim["lon"] = 26

v = defVar(ds,"var_string",String,("lon",))
v1 = defVar(ds,"var_string1",String,("lon",))

data = ["$(Char(Int('A')+i))" for i = 0:25]

v[:] = data

# element-wise
for i in 1:length(data)
    v1[i] = data[i]
end

close(ds)

ds = NCDataset(filename)
vs = ds["var_string"]
vs1 = ds["var_string"]

@test vs[:] == data
@test vs1[:] == data

@test vs[2:3] == data[2:3]

@test vs[2] == data[2]

# ignore additional indices
@test vs[2,1,1] == data[2,1,1]

# ignore additional indices
@test vs[2:3,1,1] == data[2:3,1,1]

close(ds)


text:
## IndexDicts
abstract type AbstractIndexDict{K, V} <: AbstractDict{K, V} end

makekeys(::Type{UnitRange{K}}, start::K, stop::K) where {K} = start : stop
function makekeys(::Type{Base.OneTo{K}}, start::K, stop::K) where {K}
    @boundscheck start === K(1) || error()
    Base.OneTo(stop)
end

"""
$(TYPEDEF)

An associative type whose keys are an `AbstractUnitRange`, and whose values are stored in
a `Vector`. `IndexDict` is an ordered associative collection, with the order determined by key range.
The nature of the keys enables very fast lookups and stores.

# Examples

```julia-repl
julia> IndexDict(2 : 4, [4, 5, 6])
RigidBodyDynamics.CustomCollections.IndexDict{Int64,UnitRange{Int64},Int64} with 3 entries:
  2 => 4
  3 => 5
  4 => 6

julia> IndexDict{Int32, UnitRange{Int32}}(i => 3 * i for i in Int32[4, 2, 3])
RigidBodyDynamics.CustomCollections.IndexDict{Int32,UnitRange{Int32},Int64} with 3 entries:
  2 => 6
  3 => 9
  4 => 12
```
"""
struct IndexDict{K, KeyRange<:AbstractUnitRange{K}, V} <: AbstractIndexDict{K, V}
    keys::KeyRange
    values::Vector{V}
end

Base.broadcastable(x::IndexDict) = Ref(x)

"""
$(TYPEDEF)

Like [`IndexDict`](@ref), but contains an additional `Bool` dirty bit to be used in algorithms involving cached data.
"""
mutable struct CacheIndexDict{K, KeyRange<:AbstractUnitRange{K}, V} <: AbstractIndexDict{K, V}
    keys::KeyRange
    values::Vector{V}
    dirty::Bool
    function CacheIndexDict{K, KeyRange, V}(keys::KeyRange, values::Vector{V}) where {K, V, KeyRange<:AbstractUnitRange{K}}
        @boundscheck length(keys) == length(values) || error("Mismatch between keys and values.")
        new{K, KeyRange, V}(keys, values, true)
    end
end

setdirty!(d::CacheIndexDict) = (d.dirty = true; nothing)
isdirty(d::CacheIndexDict) = d.dirty

# Constructors
for IDict in (:IndexDict, :CacheIndexDict)
    @eval begin
        function $IDict{K, KeyRange, V}(keys::KeyRange) where {K, KeyRange<:AbstractUnitRange{K}, V}
            $IDict{K, KeyRange, V}(keys, Vector{V}(undef, length(keys)))
        end

        function $IDict{K, KeyRange, V}(kv::Vector{Pair{K, V}}) where {K, KeyRange<:AbstractUnitRange{K}, V}
            if !issorted(kv, by = first)
                sort!(kv; by = first)
            end
            start, stop = if isempty(kv)
                K(1), K(0)
            else
                first(first(kv)), first(last(kv))
            end
            keys = makekeys(KeyRange, start, stop)
            for i in eachindex(kv)
                keys[i] === first(kv[i]) || error()
            end
            values = map(last, kv)
            $IDict{K, KeyRange, V}(keys, values)
        end

        function $IDict{K, KeyRange}(kv::Vector{Pair{K, V}}) where {K, KeyRange<:AbstractUnitRange{K}, V}
            $IDict{K, KeyRange, V}(kv)
        end

        function $IDict{K, KeyRange, V}(itr) where {K, KeyRange<:AbstractUnitRange{K}, V}
            kv = Pair{K, V}[]
            for x in itr
                push!(kv, convert(K, first(x)) => convert(V, last(x)))
            end
            $IDict{K, KeyRange, V}(kv)
        end

        function $IDict{K, KeyRange}(itr) where {K, KeyRange<:AbstractUnitRange{K}}
            kv = [convert(K, first(x)) => last(x) for x in itr]
            $IDict{K, KeyRange}(kv)
        end
    end
end

@inline Base.isempty(d::AbstractIndexDict) = isempty(d.values)
@inline Base.length(d::AbstractIndexDict) = length(d.values)
@inline Base.iterate(d::AbstractIndexDict, i = 1) = i > length(d) ? nothing : (d.keys[i] => d.values[i], i + 1)
@inline Base.keys(d::AbstractIndexDict{K}) where {K} = d.keys
@inline Base.values(d::AbstractIndexDict) = d.values
@inline Base.haskey(d::AbstractIndexDict, key) = key ∈ d.keys
@inline keyindex(key::K, keyrange::Base.OneTo{K}) where {K} = Int(key)
@inline keyindex(key::K, keyrange::UnitRange{K}) where {K} = Int(key - first(keyrange) + 1)
@propagate_inbounds Base.getindex(d::AbstractIndexDict{K}, key::K) where {K} = d.values[keyindex(key, d.keys)]
@propagate_inbounds Base.setindex!(d::AbstractIndexDict{K}, value, key::K) where {K} = d.values[keyindex(key, d.keys)] = value
@propagate_inbounds Base.get(d::AbstractIndexDict{K}, key::K, default) where {K} = d[key]