text:
"""
    LogLogistic(a,b)

The *log-logistic distribution* with shape `a` and scale `b` has probability density function

```math
f(x; a, b)=\\frac{(a / b )(x/b)^{a-1}}{1+(x/b)^a)^2}}
```

```julia
LogLogistic()       # LogLogistic distribution with zero location and unit scale, i.e. LogLogistic(0, 1)
LogLogistic(a)      # LogLogistic distribution with location a and unit scale, i.e. LogLogistic(a, 1)
LogLogistic(a, b)   # LogLogistic distribution with location a and scale b

params(d)       # Get the parameters, i.e. (a, b)
shape(d)     # Get the location parameter, i.e. a
scale(d)        # Get the scale parameter, i.e. b
```

External links

* [Log-Logistic distribution on Wikipedia](https://en.wikipedia.org/wiki/Log-logistic_distribution)

"""
struct LogLogistic <: ContinuousUnivariateDistribution
    a
    b
end

#### Parameters

shape(d::LogLogistic) = d.a
scale(d::LogLogistic) = d.b

params(d::LogLogistic) = (d.a, d.b)
#@inline partype(d::LogLogistic{T}) where {T<:Real} = T

#### Statistics

function mean(d::LogLogistic)
    a, b = params(d)
    if a>1
        return (b*pi/a)/sin(pi/a)
    else
        return missing
    end
end
median(d::Logistic) = d.b
#mode(d::Logistic) = d.b
#	{\displaystyle \alpha \left({\frac {\beta -1}{\beta +1}}\right)^{1/\beta }}\alpha\left(\frac{\beta-1}{\beta+1}\right)^{1/\beta}

# std(d::Logistic) = π * d.θ / sqrt3
# var(d::Logistic) = (π * d.θ)^2 / 3
# skewness(d::Logistic{T}) where {T<:Real} = zero(T)
# kurtosis(d::Logistic{T}) where {T<:Real} = T(6)/5

# entropy(d::Logistic) = log(d.θ) + 2

#### Evaluation

# zval(d::Logistic, x::Real) = (x - d.μ) / d.θ
# xval(d::Logistic, z::Real) = d.μ + z * d.θ

function pdf(d::LogLogistic, x)
    a, b = params(d)
    return ((a/b)*(x/b )^(a-1))/((1+(x/b)^a))^2
end

function logpdf(d::LogLogistic, x)
    a, b = params(d)
    return log(a)-log(b)+(a-1)*log(x/b) - 2*log1p((x/b)^a)
end

function cdf(d::LogLogistic, x)
    a, b = params(d)
    return (x^a) / (b^a +x^a)
end

function ccdf(d::LogLogistic, x)
    a, b = params(d)
    return one(x) - (x^a) / (b^a +x^a)
end

# logcdf(d::Logistic, x::Real) = -log1pexp(-zval(d, x))
# logccdf(d::Logistic, x::Real) = -log1pexp(zval(d, x))
function quantile(d::LogLogistic,p)
    a, b = params(d)
    return b * ((p)/(1-p))^inv(a)
end

function rand(rng::AbstractRNG, d::LogLogistic)
    p = rand(rng)
    return quantile(d, p)
end

# quantile(d::Logistic, p::Real) = xval(d, logit(p))
# cquantile(d::Logistic, p::Real) = xval(d, -logit(p))
# invlogcdf(d::Logistic, lp::Real) = xval(d, -logexpm1(-lp))
# invlogccdf(d::Logistic, lp::Real) = xval(d, logexpm1(-lp))

# function gradlogpdf(d::Logistic, x::Real)
#     e = exp(-zval(d, x))
#     ((2e) / (1 + e) - 1) / d.θ
# end

# mgf(d::Logistic, t::Real) = exp(t * d.μ) / sinc(d.θ * t)

# function cf(d::Logistic, t::Real)
#     a = (π * t) * d.θ
#     a == zero(a) ? complex(one(a)) : cis(t * d.μ) * (a / sinh(a))
# end

max_stars_repo_name: abikoushi/IntervalCensored

max_issues_repo_name: abikoushi/IntervalCensored

max_forks_repo_name: abikoushi/IntervalCensored


text:
module TimeSeriesA5P2

#Use this to print
#=
using Weave
codeName = "TimeSeriesA5P2"
weave(Pkg.dir("$(pwd())\\$(codeName).jl"),
  informat="script",
  out_path = "$(pwd())\\output\\$(codeName)_Appendix.html",
  doctype = "md2html")
=#

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD, ForwardDiff,
  Gadfly, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  Measures, Formatting

importall CTIO, CTReg, CTStat, CTVAR

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const YAHOO_DATE_FORMAT = "m/d/yyyy"
const FRED_DATE_FORMAT = "m/d/yyyy"
const PPI_DATE_FORMAT = "m/d/yyyy"
const NIPA_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"
const FRED_NAME = "CPI"
const PPI_NAME = "PPI"
const NIPA_NAME = "NIPA_GDP_UNRATE"

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 2
const NUM_FORMAT = "{:.$(DECIMALS)f}"
const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]
const COLORS = ["DodgerBlue", "OrangeRed", "Green",
  "Brown", "Black", "BlueViolet"]

n2s(x::T where T<:Real) = num2Str(x, 3, Ints=true)

#preprocess PPI data
function preProcessNIPA(dropTop::Bool = true)::Void

  nipaDF::DataFrame = readtable("$DATA_PATH\\$NIPA_NAME.csv")

  #read and pre-process the data
  dateFormat::String = NIPA_DATE_FORMAT
  if :Date ∈ names(nipaDF)
    rename!(nipaDF, :Date, :DATE)
  elseif :_Date ∈ names(nipaDF)
    rename!(nipaDF, :_Date, :DATE)
  elseif :_DATE ∈ names(nipaDF)
    rename!(nipaDF, :_DATE, :DATE)
  end

  nipaDF[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(nipaDF[:,:DATE])::DataVector{Date}
  sort!(nipaDF, cols = [:DATE])

  #difference
  #println("got here")
  #get GDP growth
  nipaDF[:,:gdp] = Vector{Float64}(nipaDF[:,:gdp]) #convert to float from int
  nipaDF[:,:Δlgdp] = similar(nipaDF[:,:gdp])
  nipaDF[2:end,:Δlgdp] = (log).(nipaDF[2:end,:gdp]) .- (log).(nipaDF[1:(end-1),:gdp])
  nipaDF[:,:unrate] /= 100.

  #get change in employement
  nipaDF[:,:Δunrate] = similar(nipaDF[:,:unrate])
  nipaDF[2:end,:Δunrate] = (nipaDF[2:end,:unrate]) .- (nipaDF[1:(end-1),:unrate])

  stream::IOStream = open("$DATA_PATH\\$NIPA_NAME.jls", "w")
  serialize(stream, nipaDF[2:end,:])
  close(stream)

  return nothing
end

#to get the df
function getDF(pathStr::String, refreshData::Bool, preProcessFunc::Function)::DataFrame

  #re-process the data if needed or desired
  if !isfile(pathStr) || refreshData
    preProcessFunc()
  end

  #read the binary file
  stream::IOStream = open(pathStr)
  DF::DataFrame = deserialize(stream)
  close(stream)

  return DF
end

#Calculates the summary statistics
function summaryA5P2(nipaDF::DataFrame, LHS::Vector{Symbol}, RHS::Vector{Symbol}, est::VAREstimate;
    title="Estimate")
  #label the table
  K::Int = length(LHS)
  colNames::Vector{Vector{String}} =
    [["\$\\phi_0\$"; (String).(RHS); "N"; "\$R^2\$"; "F"]]
  numCols = length(colNames[end])

  rowNames::Vector{String} = ((i::Int)->isodd(i)?(String(LHS[(i+1)÷2])):"").(1:(2K))
  numRows::Int = length(rowNames)
  content::Vector{Vector{String}} = [Vector{String}(numCols) for i::Int ∈ 1:numRows]

  #run the tests and record the results
  for k::Int ∈ 1:K
    r::Int = 2(k-1) + 1 #pointer to the beginning of each double row
    content[r][1] = n2s(est.var.ϕ₀[k])
    content[r+1][1] = est.coeffErrors[k] ≠ nothing ? "($(n2s(est.coeffErrors[k][1])))" : ""
    for c::Int ∈ 2:(numCols-3)
      content[r][c] = n2s(est.var.ϕ₁[k, c-1])
      content[r+1][c] = est.coeffErrors[k] ≠ nothing ? "($(n2s(est.coeffErrors[k][c])))" : ""
    end

    content[r][end-2] = est.N[k] ≠ nothing ? "$(n2s(est.N[k]))" : ""
    content[r+1][end-2] = ""
    content[r][end-1] = est.MSE[k] ≠ nothing ? "$(n2s(1.-est.MSE[k]/est.MST[k]))" : ""
    content[r+1][end-1] = ""
    content[r][end] = est.MSE[k] ≠ nothing ? "$(n2s(est.MSR[k]/est.MSE[k]))" : ""
    content[r+1][end] = ""

  end

  #write the table
  summaryTable::String = texTable( title,
    """See Tex File""", #caption
    colNames, #colNames
    Vector{String}(),#contentRowNames
    Vector{Matrix{String}}(), #content
    rowNames, #descRowNames
    content, #descContent
    Vector{String}(), #notes
  )

  return summaryTable
end
#problem specific graphs
function visualizeA5P2(nipaDF::DataFrame, predictions::Vector{VARPrediction})::Void

  ACFRange::Range = 1:24
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  Gadfly.push_theme(:default)
  outDF = deepcopy(nipaDF)
  completecases!(outDF)

  push!(plotNames, "NIPA Summary Graphs")
  push!(plots, Vector{Plot}())
  for series::Symbol ∈ [:Δlgdp, :unrate]
    push!(plots[end],
      plot(outDF, x=:DATE, y=series,  Geom.line, Geom.point,
        Guide.xlabel("Date"),Guide.ylabel("$series"),
        Guide.title("Summary $series")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  push!(plotNames, "UnrestΔlGDP")
  push!(plots, plot(
    layer(x=predictions[1].dates, y=predictions[1].Y[1,:],  Geom.line,
      Theme(default_color=COLORS[1])),
    layer(x=predictions[1].dates, y=predictions[1].YHat[1,:],  Geom.line,
      Theme(default_color=COLORS[2])),
    layer(x=predictions[2].dates, y=predictions[2].YHat[1,:],  Geom.line,
      Theme(default_color=COLORS[3])),
    Guide.manual_color_key("Legend", ["Actual", "t+1 Prediction", "t+4 Prediction"], COLORS[1:3]),
    Guide.xlabel("Date"),Guide.ylabel("ΔlGDP"),
    Guide.title("Unrestricted VAR Prediction of GDP")))

    push!(plotNames, "UnrestUNRATE")
    push!(plots, plot(
      layer(x=predictions[1].dates, y=predictions[1].Y[2,:],  Geom.line,
        Theme(default_color=COLORS[1])),
      layer(x=predictions[1].dates, y=predictions[1].YHat[2,:],  Geom.line,
        Theme(default_color=COLORS[2])),
      layer(x=predictions[2].dates, y=predictions[2].YHat[2,:],  Geom.line,
        Theme(default_color=COLORS[3])),
      Guide.manual_color_key("Legend", ["Actual", "t+1 Prediction", "t+4 Prediction"], COLORS[1:3]),
      Guide.xlabel("Date"),Guide.ylabel("Unemployment"),
      Guide.title("Unrestricted VAR Prediction of Unemployment")))

    push!(plotNames, "RestrictedΔlGDP")
    push!(plots, plot(
      layer(x=predictions[3].dates, y=predictions[3].Y[1,:],  Geom.line,
        Theme(default_color=COLORS[1])),
      layer(x=predictions[3].dates, y=predictions[3].YHat[1,:],  Geom.line,
        Theme(default_color=COLORS[2])),
      layer(x=predictions[4].dates, y=predictions[4].YHat[1,:],  Geom.line,
        Theme(default_color=COLORS[3])),
      Guide.manual_color_key("Legend", ["Actual", "t+1 Prediction", "t+4 Prediction"], COLORS[1:3]),
      Guide.xlabel("Date"),Guide.ylabel("ΔlGDP"),
      Guide.title("Restricted VAR Prediction of GDP")))

  #write the graphs
  for i::Int ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 9inch, 7inch), plots[i])
    #display(plots[i])
  end

  return nothing
end

function lagDF!(DF::DataFrame, target::Symbol, lags::Int = 1; name= Symbol("$(target)L$lags"))::DataFrame
  DF[:,name] = similar(DF[:,target])
  DF[:,name] .= NA
  DF[(lags+1):end,name] = DF[1:(end-lags),target]

  return DF
end

#script for the problem
function A5P2Script(;refreshData::Bool = true)::Void

  LHS::Vector{Symbol} = [:Δlgdp, :unrate, :ΔlgdpL1, :unrateL1]
  RHS::Vector{Symbol} = [:ΔlgdpL1, :unrateL1, :ΔlgdpL2, :unrateL2]
  nipaDF::DataFrame = getDF("$DATA_PATH\\$NIPA_NAME.jls", refreshData, preProcessNIPA)

  for s ∈ LHS[1:2]
    lagDF!(nipaDF,s, 1)
    lagDF!(nipaDF,s, 2)
  end

  ######first, estimate the unrestricted model
  est::VAREstimate = VAREstimate(nipaDF, LHS, RHS)
  tables::Vector{String} = Vector{String}()
  push!(tables, summaryA5P2(nipaDF, LHS, RHS, est, title="Unrestricted VAR 2"))

  #get the predicted values
  validXY::Vector{Bool} = completecases(nipaDF[RHS]) .& completecases(nipaDF[LHS])
  X::Matrix{Float64} = Matrix{Float64}(nipaDF[validXY, RHS])'
  Y::Matrix{Float64} = Matrix{Float64}(nipaDF[validXY, LHS])'
  dates::Vector{Date} = nipaDF[validXY, :DATE]

  #container for holding all predictions
  predictions::Vector{VARPrediction} = Vector{VARPrediction}()

  push!(predictions, VARPrediction(est.var, dates, X, Y, τ = 1))
  push!(predictions, VARPrediction(est.var, dates, X, Y, τ = 4))

  ##################Now do the restricted model
  LHS = [:Δlgdp, :Δunrate, :ΔlgdpL1]
  RHS = [:ΔlgdpL1, :ΔunrateL1, :ΔlgdpL2]
  lagDF!(nipaDF,:Δunrate, 1)

  #estimate the VAR
  est = VAREstimate(nipaDF, LHS, RHS)
  push!(tables, summaryA5P2(nipaDF, LHS, RHS, est, title="Restricted VAR 2"))

  #get the predicted values
  validXY = completecases(nipaDF[RHS]) .& completecases(nipaDF[LHS])
  X = Matrix{Float64}(nipaDF[validXY, RHS])'
  Y = Matrix{Float64}(nipaDF[validXY, LHS])'
  dates= nipaDF[validXY, :DATE]

  push!(predictions, VARPrediction(est.var, dates, X, Y, τ = 1))
  push!(predictions, VARPrediction(est.var, dates, X, Y, τ = 4))

  #make the grpahs
  visualizeA5P2(nipaDF, predictions)

  #write out the tables
  writeTables2File(tables,
    HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
    outName = "TimeSeries A5P2.tex")

  return nothing
end

@time begin
# Uncomment to run
#  A5P2Script(refreshData = true)
end

end

ext: jl

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function read_output(file_path::String)
	file_stream::IOStream = open(file_path, "r")
	file_string::String = readstring(file_stream)
	close(file_stream)

	number_of_samples::Int64 = read_key(file_string, "number_of_samples", Int64)
	d2::Array{Float64, 1} = read_key(file_string, "d2", Array{Float64, 1})
	S2::Array{Float64, 1} = read_key(file_string, "S2", Array{Float64, 1})
	d3::Array{Float64, 1} = read_key(file_string, "d3", Array{Float64, 1})
	theta3::Array{Float64, 1} = read_key(file_string, "theta3", Array{Float64, 1})
	S3::Array{Float64, 1} = read_key(file_string, "S3", Array{Float64, 1})
	t_exec::Float64 = read_key(file_string, "execution_time", Float64)

	return (
		number_of_samples,
		d2,
		S2,
		d3,
		theta3,
		S3,
		t_exec)
end


text:
#---------------------------------------------------------------------#
#This function computes the Initial & Analytic Solutions.
#Written by F.X. Giraldo on 10/2003
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function exact_solution(coord,Npoin,c,DFloat)

    #Initialize
    qe=zeros(DFloat,Npoin,1)
    fe=zeros(DFloat,Npoin,1)
    cc=c*π

    #Generate Grid Points
    for I=1:Npoin
        x=coord[1,I]
        y=coord[2,I]
        qe[I]=sin(cc*x)*sin(cc*y)
        fe[I]=-2*cc^2*sin(cc*x)*sin(cc*y)
    end #I

    return (qe,fe)

end


text:
module LoadBalancing

include("./../constraints/geometric/constants.jl")
include("./1_constraints.jl")
using JuMP
using .Constants
using .LoadBalancingConstraints

export cons_lb_load_balancing

function cons_lb_load_balancing(m)

    m = lb_constraints(m)

    return m
end

end


text:
include("generic_functions.jl")
#example
output_instance_path = "..\\data_ROPF"
output_decomposition_path = "..\\data_sdp"
list_generation = ["generation"]
generation_files_path = "..\\data_ROPF"
max_time = 600 #10 minutes
network_name = "14nodes_disconnect"

f = open("results_train_$(network_name)_v+-0.02_BB$(max_time).csv", "w")
write(f, "Numero ; UB ; LB ; Gap (%)\n")
for i in 0:2499
    # network_name, numero = split(instance_name, '-')
    numero = "train_$i"
    instance_name = network_name*"-"*numero
    matpower_instance_path = "..\\mp_data\\$(network_name)\\$(numero).mat"
    for generation in list_generation
        construct_dat_file_ROPF(instance_name, matpower_instance_path, output_instance_path)
        # generate_clique_decomposition(instance_name, matpower_instance_path, output_decomposition_path)
        ROPF = ROPF_infos(instance_name,
        matpower_instance_path,
        output_instance_path,
        "cholesky",
        output_decomposition_path,
        generation,
        generation_files_path)
        UB, LB = solve1_refgen(ROPF)
        println("UB=$UB ; LB=$LB")
        UB_BB, LB = solve2_refgen(ROPF, max_time)
        println("After B&B")
        println("UB=$UB_BB ; LB=$LB")
        if UB_BB != Inf
            if UB_BB < UB
                cp("solutions\\BB_refgen_solution_$(instance_name).csv", "solutions_mp_data\\refgen_solution_$(instance_name).csv")
            else
                cp("solutions\\solve1_refgen_solution_$(instance_name).csv", "solutions_mp_data\\refgen_solution_$(instance_name).csv")
            end
        end
        gap = (min(UB, UB_BB)-LB)/(min(UB, UB_BB))
        write(f, "$numero ; $(min(UB, UB_BB)) ; $LB ; $gap \n")
    end
end
close(f)


text:
using RandomExtensions, Random, SparseArrays
using Random: Sampler, gentype

@testset "Distributions" begin
    # Normal/Exponential
    @test rand(Normal()) isa Float64
    @test rand(Normal(0.0, 1.0)) isa Float64
    @test rand(Normal(0, 1)) isa Float64
    @test rand(Normal(0, 1.0)) isa Float64
    @test rand(Exponential()) isa Float64
    @test rand(Exponential(1.0)) isa Float64
    @test rand(Exponential(1)) isa Float64
    @test rand(Normal(Float32)) isa Float32
    @test rand(Exponential(Float32)) isa Float32
    @test rand(Normal(ComplexF64)) isa ComplexF64

    # pairs/complexes
    @test rand(make(Pair, 1:3, Float64)) isa Pair{Int,Float64}
    @test rand(make(Pair{Int8}, 1:3, Float64)) isa Pair{Int8,Float64}
    @test rand(make(Pair{Int8,Float32}, 1:3, Float64)) isa Pair{Int8,Float32}
    @test rand(make(Pair{X,Float32} where X, 1:3, Float64)) isa Pair{Int,Float32}
    @test rand(Pair{Int,Float64}) isa Pair{Int,Float64}

    z = rand(make(Complex, 1:3, 6:9))
    @test z.re ∈ 1:3
    @test z.im ∈ 6:9
    @test z isa Complex{Int}
    z = rand(make(ComplexF64, 1:3, 6:9))
    @test z.re ∈ 1:3
    @test z.im ∈ 6:9
    @test z isa ComplexF64
    for (C, R) in ((Complex, Int), (ComplexF64, Float64), (Complex{Int}, Int))
        z = rand(make(C, 1:3))
        @test z.re ∈ 1:3
        @test z.im ∈ 1:3
        @test z isa Complex{R}
    end
    @test rand(make(Complex, 1:3, Float64)) isa Complex{Float64} # promote_type should be used

    @test rand(ComplexF64) isa ComplexF64

    @test rand(make(Complex,Int), 3) isa Vector{Complex{Int}}
    @test rand(make(Complex,1:3), 3) isa Vector{Complex{Int}}

    # Uniform
    @test rand(Uniform(Float64)) isa Float64
    @test rand(Uniform(1:10)) isa Int
    @test rand(Uniform(1:10)) ∈ 1:10
    @test rand(Uniform(Int)) isa Int

    # Bernoulli
    @test rand(Bernoulli()) ∈ (0, 1)
    @test rand(Bernoulli(1)) == 1
    @test rand(Bernoulli(0)) == 0
    # TODO: do the math to estimate proba of failure:
    @test 620 < count(rand(Bernoulli(Bool, 0.7), 1000)) < 780
    for T = (Bool, Int, Float64, ComplexF64)
        r = rand(Bernoulli(T))
        @test r isa T
        @test r ∈ (0, 1)
        r = rand(Bernoulli(T, 1))
        @test r == 1
    end

    # Categorical
    n = rand(1:9)
    @test rand(Categorical(n)) ∈ 1:9
    @test all(∈(1:9), rand(Categorical(n), 10))
    @test rand(Categorical(n)) isa Int
    c = Categorical(Float64(n))
    @test rand(c) isa Float64
    @test rand(c) ∈ 1:9
    c = Categorical([1, 7, 2])
    # cf. Bernoulli tests
    @test 620 < count(==(2), rand(c, 1000)) < 780
    @test rand(c) isa Int
    @test rand(Categorical{Float64}((1, 2, 3, 4))) isa Float64

    @test_throws ArgumentError Categorical(())
    @test_throws ArgumentError Categorical([])
    @test_throws ArgumentError Categorical(x for x in 1:0)
end

const rInt8 = typemin(Int8):typemax(Int8)
const spString = Sampler(MersenneTwister, String)

const RNGS = ([]                   => "global RNG",
              [MersenneTwister(0)] => "MersenneTwister",
              [RandomDevice()]     => "RandomDevice")

@testset "Containers $name" for (rng, name) in RNGS
    # Array
    for T = (Int, Int8)
        for (A, AT) = ((Array, Int8), (Array{T}, T), (Vector, Int8), (Vector{T}, T))
            @inferred rand(rng..., Int8.(1:9), A, 10)
            a = rand(rng..., Int8.(1:9), A, 10)
            @test a isa Vector{AT}
            @test all(in(1:9), a)
            @inferred rand(rng..., Int8, A, 10)
            a = rand(rng..., Int8, A, 10)
            @test a isa Vector{AT}
            @test all(in(rInt8), a)
        end
    end

    # Set
    for S = (Set{Int}, Set, BitSet)
        s = rand(rng..., 1:99, S, 10)
        @test s isa (S === BitSet ? BitSet : Set{Int})
        @test length(s) == 10
        @test rand(s) ∈ 1:99
    end
    for s = (Set([1, 2]), BitSet([1, 2]))
        @test s === rand!(s)
        @test s != Set([1, 2]) # very unlikely
        @test length(s) == 2
        @test s === rand!(s, 3:9) <= Set(3:9)
        @test length(s) == 2
    end
    @test rand(rng..., Pair{Int,Float64}, Set, 3) isa Set{Pair{Int,Float64}}
    @test rand(rng..., Pair{Int,Float64}, Set{Pair}, 3) isa Set{Pair}

    # BitSet
    s = rand(make(BitSet, 1:10, 3))
    @test s isa BitSet
    @test length(s) == 3
    @test s <= Set(1:10)

    @testset "default_sampling(::BitSet) == Int8" begin
        Random.seed!(0)
        rand!(s)
        @test s <= Set(rInt8)
        Random.seed!(0)
        @test s == rand(BitSet, 3)
    end

    # Dict
    for s = (rand(rng..., make(Pair, 1:99, 1:99), Dict, 10),
             rand(rng..., make(Pair, 1:99, 1:99), Dict{Int,Int}, 10))
        @test s isa Dict{Int,Int}
        @test length(s) == 10
        p = rand(s)
        @test p.first ∈ 1:99
        @test p.second ∈ 1:99
    end
    s = Dict(1=>2, 2=>1)
    @test s === rand!(s)
    @test length(s) == 2
    @test first(s).first ∉ (1, 2) # extremely unlikely
    rand!(s, make(Pair, 3:9, Int))
    @test length(s) == 2
    @test first(s).first ∈ 3:9

    d = rand(rng..., Pair{Int,Float64}, Dict, 3)
    @test d isa Dict{Int,Float64}
    dd = rand!(rng..., d, Pair{Int,Int8})
    @test dd === d
    delt = pop!(d)
    @test delt isa Pair{Int,Float64}
    @test delt[2] ∈ rInt8
    @test rand(rng..., Pair{Int,Float64}, Dict{Any,Any}, 3) isa Dict{Any,Any}

    # sparse
    @test rand(rng..., Float64, .5, 10) isa SparseVector{Float64}
    @test rand(rng..., Float64, .5, (10,)) isa SparseVector{Float64}
    @test rand(rng..., Float64, SparseVector, .5, 10) isa SparseVector{Float64}
    @test rand(rng..., Float64, SparseVector, .5, (10,)) isa SparseVector{Float64}

    @test rand(rng..., .5, 10) isa SparseVector{Float64}
    @test rand(rng..., .5, (10,)) isa SparseVector{Float64}
    @test rand(rng..., SparseVector, .5, 10) isa SparseVector{Float64}
    @test rand(rng..., SparseVector, .5, (10,)) isa SparseVector{Float64}

    @test rand(rng..., Int, .5, 10) isa SparseVector{Int}
    @test rand(rng..., Int, .5, (10,)) isa SparseVector{Int}
    @test rand(rng..., Int, SparseVector, .5, 10) isa SparseVector{Int}
    @test rand(rng..., Int, SparseVector, .5, (10,)) isa SparseVector{Int}

    @test rand(rng..., Float64, .5, 10, 3) isa SparseMatrixCSC{Float64}
    @test rand(rng..., Float64, .5, (10, 3)) isa SparseMatrixCSC{Float64}
    @test rand(rng..., Float64, SparseMatrixCSC, .5, 10, 3) isa SparseMatrixCSC{Float64}
    @test rand(rng..., Float64, SparseMatrixCSC, .5, (10, 3)) isa SparseMatrixCSC{Float64}

    @test rand(rng..., .5, 10, 3) isa SparseMatrixCSC{Float64}
    @test rand(rng..., .5, (10, 3)) isa SparseMatrixCSC{Float64}
    @test rand(rng..., SparseMatrixCSC, .5, 10, 3) isa SparseMatrixCSC{Float64}
    @test rand(rng..., SparseMatrixCSC, .5, (10, 3)) isa SparseMatrixCSC{Float64}

    @test rand(rng..., Int, .5, 10, 3) isa SparseMatrixCSC{Int}
    @test rand(rng..., Int, .5, (10, 3)) isa SparseMatrixCSC{Int}
    @test rand(rng..., Int, SparseMatrixCSC, .5, 10, 3) isa SparseMatrixCSC{Int}
    @test rand(rng..., Int, SparseMatrixCSC, .5, (10, 3)) isa SparseMatrixCSC{Int}

    # BitArray
    for S = ([], [Bool], [Bernoulli()])
        @test rand(rng..., S..., BitArray, 10) isa BitVector
        @test rand(rng..., S..., BitVector, 10) isa BitVector
        @test_throws MethodError rand(rng..., S..., BitVector, 10, 20) isa BitVector
        @test rand(rng..., S..., BitArray, 10, 3) isa BitMatrix
        @test rand(rng..., S..., BitMatrix, 10, 3) isa BitMatrix
        @test_throws MethodError rand(rng..., S..., BitVector, 10, 3) isa BitMatrix
    end

    # String
    s = rand(rng..., String)
    @test s isa String
    @test length(s) == 8
    s = rand(rng..., String, 10)
    @test s isa String
    @test length(s) == 10
    s = rand(rng..., "asd", String)
    @test length(s) == 8
    @test Set(s) <= Set("asd")
    @test_throws ArgumentError rand(rng..., String, 2, 3)
    @test_throws ArgumentError rand(rng..., String, 2, 3, 4)
    @test_throws ArgumentError rand(rng..., String, 2, 3, 4, 5)

    # Tuple
    s = rand(rng..., Int, NTuple{3})
    @test s isa NTuple{3,Int}
    s = rand(rng..., 1:3, NTuple{3})
    @test s isa NTuple{3,Int}
    @test all(in(1:3), s)
    s = rand(rng..., 1:3, NTuple{3,Int8})
    @test s isa NTuple{3,Int8}
    @test all(in(1:3), s)

    s = rand(rng..., NTuple{3, Int8})
    @test s isa NTuple{3,Int8}

    s = rand(rng..., Tuple{Int8, UInt8})
    @test s isa Tuple{Int8, UInt8}
    s = rand(rng..., 1:3, Tuple{Int8, UInt8})
    @test s isa Tuple{Int8, UInt8}
    @test all(in(1:3), s)

    s = rand(rng..., 1:3, Tuple, 4)
    @test s isa NTuple{4,Int}
    @test all(in(1:3), s)
    s = rand(rng..., Tuple, 4)
    @test s isa NTuple{4,Float64}

    s = rand(rng..., NTuple{3})
    @test s isa NTuple{3,Float64}

    s = rand(rng..., NTuple{N,UInt8} where N, 3)
    @test s isa NTuple{3,UInt8}
    s = rand(rng..., 1:3, NTuple{N,UInt8} where N, 3)
    @test s isa NTuple{3,UInt8}
    @test all(in(1:3), s)
end

@testset "Rand $name" for (rng, name) in RNGS
    for XT = zip(([Int], [1:3], []), (Int, Int, Float64))
        X, T = XT
        r = Rand(rng..., X...)
        @test collect(Iterators.take(r, 10)) isa Vector{T}
        @test r() isa T
        @test r(2, 3) isa Matrix{T}
        @test r(.3, 2, 3) isa SparseMatrixCSC{T}
    end
    for d = (Uniform(1:10), Uniform(Int))
        @test collect(Iterators.take(d, 10)) isa Vector{Int}
    end
end

struct PairDistrib <: RandomExtensions.Distribution{Pair}
end

Random.rand(rng::AbstractRNG, ::Random.SamplerTrivial{PairDistrib}) = 1=>2

@testset "allow abstract Pair when generating a Dict" begin
    d = rand(PairDistrib(), Dict, 1)
    @test d == Dict(1=>2)
    @test typeof(d) == Dict{Any,Any}
end

@testset "some tight typing" begin
    UI = Random.UInt52()
    @test eltype(rand(MersenneTwister(), Random.Sampler(MersenneTwister, UI), .6, 1, 0)) == UInt64
    @test eltype(rand(UI, Set, 3)) == UInt64
    @test eltype(rand(Uniform(UI), 3)) == UInt64
    a = rand(make(Pair, Int, UI))
    @test fieldtype(typeof(a), 2) == UInt64
end

#=
@testset "rand(::Pair)" begin
    @test rand(1=>3) ∈ (1, 3)
    @test rand(1=>2, 3) isa Vector{Int}
    @test rand(1=>'2', 3) isa Vector{Union{Char, Int}}
end
=#

@testset "rand(::AbstractFloat)" begin
    # check that overridden methods still work
    m = MersenneTwister()
    for F in (Float16, Float32, Float64, BigFloat)
        @test rand(F) isa F
        sp = Random.Sampler(MersenneTwister, RandomExtensions.CloseOpen01(F))
        @test rand(m, sp) isa F
        @test 0 <= rand(m, sp) < 1
        for (CO, (l, r)) = (CloseOpen  => (<=, <),
                            CloseClose => (<=, <=),
                            OpenOpen   => (<,  <),
                            OpenClose  => (<,  <=))
            f = rand(CO(F))
            @test f isa F
            @test l(0, f) && r(f, 1)
        end
        F ∈ (Float64, BigFloat) || continue # only types implemented in Random
        sp = Random.Sampler(MersenneTwister, RandomExtensions.CloseOpen12(F))
        @test rand(m, sp) isa F
        @test 1 <= rand(m, sp) < 2
    end
    @test CloseOpen(1,   2)          === CloseOpen(1.0, 2.0)
    @test CloseOpen(1.0, 2)          === CloseOpen(1.0, 2.0)
    @test CloseOpen(1,   2.0)        === CloseOpen(1.0, 2.0)
    @test CloseOpen(1.0, Float32(2)) === CloseOpen(1.0, 2.0)
    @test CloseOpen(big(1), 2) isa CloseOpen{BigFloat}

    for CO in (CloseOpen, CloseClose, OpenOpen, OpenClose)
        @test_throws ArgumentError CO(1, 1)
        @test_throws ArgumentError CO(2, 1)

        @test CO(Float16(1), 2) isa CO{Float16}
        @test CO(1, Float32(2)) isa CO{Float32}
    end
end

@testset "rand(::Type{<:Tuple})" begin
    for types in ([Base.BitInteger_types..., Float16, Float32, Float64, BigFloat, Char, Bool],
                  [Int, UInt64, Char]) # more repetitions
        tlist = rand(types, rand(0:10))
        T = Tuple{tlist...}
        @test rand(T) isa Tuple{tlist...}
    end
    @test rand(Tuple{}) === ()
    sp = Sampler(MersenneTwister, Tuple)
    @test gentype(sp) == Tuple{}
    @test rand(sp) == ()
    sp = Sampler(MersenneTwister, NTuple{3})
    @test gentype(sp) == NTuple{3,Float64}
    @test rand(sp) isa NTuple{3,Float64}
    sp = Sampler(MersenneTwister, Tuple{Int8,UInt8})
    @test gentype(sp) == Tuple{Int8,UInt8}
    @test rand(sp) isa Tuple{Int8,UInt8}
end

@testset "rand(make(Tuple, ...))" begin
    s = rand([Char, Int, Float64, Bool, 1:3, "abcd", Set([1, 2, 3])], rand(0:10))
    @test rand(make(Tuple, s...)) isa Tuple{Random.gentype.(s)...}
    # explicit test for corner case:
    @test rand(make(Tuple)) == ()
    @test rand(make(Tuple{})) == ()

    t = rand(make(Tuple, 1:3, Char, Int))
    @test t[1] ∈ 1:3
    @test t[2] isa Char
    @test t[3] isa Int && t[3] ∉ 1:3 # extremely unlikely

    t = rand(make(Tuple{Int8,Char,Int128}, 1:3, Char, Int8))
    @test t[1] isa Int8 && t[1] ∈ 1:3
    @test t[2] isa Char
    @test t[3] isa Int128 && t[3] ∈ rInt8
    @test_throws ArgumentError make(Tuple{Int}, 1:3, 1:3)
    @test_throws ArgumentError make(Tuple{Int,Int,Int}, 1:3, 2:4)

    @test rand(make(Tuple, spString, String)) isa Tuple{String,String}

    @test rand(make(Tuple{Int8,Int8})) isa Tuple{Int8,Int8}
    @test rand(make(Tuple{Int8,UInt})) isa Tuple{Int8,UInt}

    # make(Tuple, s, n)
    s = rand(make(Tuple, 1:3, 4))
    @test s isa NTuple{4,Int}
    @test all(in(1:3), s)
    s = rand(make(Tuple, Int8, 4))
    @test s isa NTuple{4,Int8}
    s = rand(make(Tuple, 4))
    @test s isa NTuple{4,Float64}
end

@testset "rand(make(NTuple{N}/Tuple{...}, x))" begin
    s, N = rand([Char, Int, Float64, Bool, 1:3, "abcd", Set([1, 2, 3])]), rand(0:10)
    T = Random.gentype(s)
    rand(make(NTuple{N}, s)) isa NTuple{N,T}
    @test rand(make(NTuple{3}, spString)) isa NTuple{3,String}
    @test rand(make(NTuple{3,UInt8}, 1:3)) isa NTuple{3,UInt8}
    @test rand(make(Tuple{Integer,Integer}, 1:3)) isa Tuple{Int,Int}
    r = rand(make(Tuple{AbstractFloat,AbstractFloat}, 1:3))
    @test r isa Tuple{Float64,Float64}
    @test all(∈(1.0:3.0), r)

    r = rand(make(Tuple{AbstractFloat,Integer}, 1:3))
    @test r isa Tuple{Float64,Int64}
    @test all(in(1:3), r)

    r = rand(make(NTuple{3}))
    @test r isa NTuple{3,Float64}

    r = rand(make(NTuple{N,UInt8} where N, 3))
    @test r isa NTuple{3,UInt8}
    r = rand(make(NTuple{N,UInt8} where N, 1:3, 3))
    @test r isa NTuple{3,UInt8}
    @test all(in(1:3), r)
    r = rand(make(NTuple{N,UInt8} where N, UInt8, 3))
    @test r isa NTuple{3,UInt8}
end

@testset "NamedTuple" begin
    for t = (rand(make(NamedTuple)), rand(NamedTuple))
        @test t == NamedTuple()
    end

    for t = (rand(make(NamedTuple{(:a,)})), rand(NamedTuple{(:a,)}))
        @test t isa NamedTuple{(:a,), Tuple{Float64}}
    end
    for t = (rand(make(NamedTuple{(:a,),Tuple{Int}})),
             rand(NamedTuple{(:a,),Tuple{Int}}))
        @test t isa NamedTuple{(:a,), Tuple{Int}}
    end

    t = rand(make(NamedTuple{(:a,)}, 1:3))
    @test t isa NamedTuple{(:a,), Tuple{Int}}
    @test t.a ∈ 1:3
    t = rand(make(NamedTuple{(:a,),Tuple{Float64}}, 1:3))
    @test t isa NamedTuple{(:a,), Tuple{Float64}}
    @test t.a ∈ 1:3

    for t = (rand(make(NamedTuple{(:a, :b)})),
             rand(NamedTuple{(:a, :b)}))
        @test t isa NamedTuple{(:a, :b), Tuple{Float64,Float64}}
    end
    for t = (rand(make(NamedTuple{(:a, :b),Tuple{Int,UInt8}})),
             rand(NamedTuple{(:a, :b),Tuple{Int,UInt8}}))
        @test t isa NamedTuple{(:a, :b), Tuple{Int,UInt8}}
    end
    t = rand(make(NamedTuple{(:a, :b)}, 1:3))
    @test t isa NamedTuple{(:a, :b), Tuple{Int,Int}}
    @test t.a ∈ 1:3 && t.b ∈ 1:3
    t = rand(make(NamedTuple{(:a, :b),Tuple{Float64,UInt8}}, 1:3))
    @test t isa NamedTuple{(:a, :b), Tuple{Float64,UInt8}}
    @test t.a ∈ 1:3 && t.b ∈ 1:3

    # as container
    @test rand(1:3, NamedTuple{(:a,)}) isa NamedTuple{(:a,), Tuple{Int}}
    @test rand(1:3, NamedTuple{(:a,), Tuple{Float64}}) isa NamedTuple{(:a,), Tuple{Float64}}
    @test rand(1:3, NamedTuple{(:a, :b)}) isa NamedTuple{(:a, :b), Tuple{Int,Int}}
    @test rand(1:3, NamedTuple{(:a, :b), Tuple{Float64,Float64}}) isa NamedTuple{(:a, :b), Tuple{Float64,Float64}}
end

@testset "rand(make(String, ...))" begin
    b = UInt8['0':'9';'A':'Z';'a':'z']

    for (s, c, n) in [(rand(String), b, 8),
                      (rand(make(String, 3)), b, 3),
                      (rand(make(String, "asd")), "asd", 8),
                      (rand(make(String, 3, "asd")), "asd", 3),
                      (rand(make(String, "qwe", 3)), "qwe", 3)]

        @test s ⊆ map(Char, c)
        @test length(s) == n
    end
    @test rand(make(String, Char)) isa String
    @test rand(make(String, 3, Char)) isa String
    @test rand(make(String, Sampler(MersenneTwister, ['a', 'b', 'c']), 10)) isa String
end

@testset "rand(make(Set/BitSet, ...))" begin
    for (S, SS, (low, high)) = ((Set{Int}, Set{Int}, (typemin(Int), typemax(Int))),
                                (Set,      Set,      (0, 1)),
                                (BitSet,   BitSet,   (typemin(Int8), typemax(Int8))))
        for (k, l) = ([1:9] => 1:9, [Int8] => rInt8, [] => ())
            s = rand(make(S, k..., 3))
            @test s isa (SS === Set ? (l == () ? Set{Float64} : Set{eltype(l)}) : SS)
            @test length(s) == 3
            if l == ()
                @test all(x -> low <= x <= high, s)
            else
                @test all(in(l), s)
            end
        end
        rand(make(S, Sampler(MersenneTwister, 1:99), 9)) isa Union{BitSet, Set{Int}}
    end
end

@testset "rand(make(Dict, ...))" begin
    for BD     = (Dict, Base.ImmutableDict),
        (D, S) = (BD{Int16,Int16}     => [],
                  BD{Int16}           => [Pair{Int8,Int16}],
                  BD{K,Int16} where K => [Pair{Int16,Int8}],
                  BD                  => [Pair{Int16,Int16}])

        d = rand(make(D, S..., 3))
        @test d isa BD{Int16,Int16}
        @test length(d) == 3
    end
end

@testset "rand(make(Array/BitArray, ...))" begin
    for (T, Arr) = (Bool => BitArray, Float64 => Array{Float64}),
        k = ([], [T], [Bernoulli(T, 0.3)]),
        (d, dim) = ([(6,)]              => 1,
                    [(2,3)]             => 2,
                    [6]                 => 1,
                    [2, 3]              => 2,
                    [Int8(2), Int16(3)] => 2),
        A = (T == Bool ?
             (BitArray, BitArray{dim}) :
             (Array, Array{Float64}, Array{Float64,dim}, Array{U,dim} where U))

        s = rand(make(A, k..., d...))
        @test s isa  Arr{dim}
        @test length(s) == 6
    end
    @test_throws MethodError rand(make(Matrix, 2))
    @test_throws MethodError rand(make(Vector, 2, 3))
    @test_throws MethodError rand(make(BitMatrix, 2))
    @test_throws MethodError rand(make(BitVector, 2, 3))

    @test rand(make(Array, spString, 9)) isa Array{String}
    @test rand(make(BitArray, Sampler(MersenneTwister, [0, 0, 0, 1]), 9)) isa BitArray
    # TODO: below was testing without explicit `Array` as 1st argument, so this test
    # may now be obsolete, redundant with tests above
    for dims = ((), (2, 3), (0x2, 3), (2,), (0x2,)),
        s    = ([], [Int], [1:3])

        T = s == [] ? Float64 : Int
        if dims != ()
            a = rand(make(Array, s..., dims...))
            @test a isa Array{T,length(dims)}
            if s == [1:3]
                @test all(in(1:3), a)
            end
        end
        if dims isa Dims
            a = rand(make(Array, s..., dims))
            @test a isa Array{T,length(dims)}
            if s == [1:3]
                @test all(in(1:3), a)
            end
        end
    end
end

@testset "rand(make(Sparse...))" begin
    for k = ([], [Float64], [Bernoulli(Float64, 0.3)]),
        (d, dim) = ([(6,)]              => 1,
                    [(2,3)]             => 2,
                    [6]                 => 1,
                    [2, 3]              => 2,
                    [Int8(2), Int16(3)] => 2)

        typ = dim == 1 ? SparseVector : SparseMatrixCSC

        s = rand(make(typ, k..., 0.3, d...))
        @test s isa (dim == 1 ? SparseVector{Float64,Int} :
                                SparseMatrixCSC{Float64,Int})
        @test length(s) == 6
    end
    @test rand(make(SparseVector, spString, 0.3, 9)) isa SparseVector{String}
    @test rand(make(SparseVector, make(SparseMatrixCSC, 1:9, 0.3, 2, 3), .1, 4)) isa SparseVector{SparseMatrixCSC{Int64,Int64},Int64}
end

@testset "rand(make(default))" begin
    @test rand(make()) isa Float64
    @test rand(make(1:3)) isa Int
    @test rand(make(1:3)) ∈ 1:3
    @test rand(make(Float64)) isa Float64
end

@testset "rand(T => x) & rand(T => (x, y, ...))" begin
    @test rand(Complex => Int) isa Complex{Int}
    @test rand(Pair => (String, Int8)) isa Pair{String,Int8}
    @test_throws MethodError rand(1=>2) # calls rand(make(1, 2))

    @test rand(Complex => Int, 3) isa Vector{Complex{Int}}
    @test rand(Pair => (String, Int8), Set, 3) isa Set{Pair{String,Int8}}

    nt = rand(NTuple{4} => Complex => 1:3)
    @test nt isa NTuple{4,Complex{Int64}}
end

##  make(x1, xs...)

struct MyType
    x
end

RandomExtensions.maketype(x::MyType, y) = eltype(x.x:y)

Base.rand(rng::AbstractRNG,
          x::Random.SamplerTrivial{<:RandomExtensions.Make2{<:Integer, MyType}}) =
              rand(rng, x[][1].x:x[][2])

@testset "rand(make(CustomType(), ...))" begin
    m = MyType(3)
    @test rand(make(m, 3)) == 3
    @test rand(make(m, big(5))) isa BigInt
    @test rand(make(m, big(5))) ∈ 3:5
    a = rand(make(m, big(6)), 5)
    @test a isa Vector{BigInt}
    @test length(a) == 5
    @test all(∈(3:6), a)
end

## Make getindex

@testset "Make getindex" begin
    m = make(Pair, 1:2, Bool, make(String, 3))
    @test m isa RandomExtensions.Make3
    @test m[1] == 1:2
    @test m[0x2] == Bool
    @test m[3] isa RandomExtensions.Make2
    @test m[1:2] == (m[1], m[2])
    @test m[[2, 3]] == m[2:end] == (m[2], m[3])
    @test m[big(3)][1:end] == m[3][1:2]
end

## @rand

run_once = false

struct Die
    n::Int
end

Base.eltype(::Type{Die}) = Int

struct DieT{T}
    n::T
end

Base.eltype(::Type{DieT{T}}) where {T} = T

@testset "@rand" begin
    must_run = !run_once
    global run_once = true

    if must_run
        # rng0 to be sure rng is not accessed in `esc`aped body of @rand
        rng0 = MersenneTwister()

        d = Die(6)
        @rand function rand(d::Die)
            7
        end
        @test rand(d) == 7
        @rand function (d::Die) 7 end
        @test rand(d) == 7
        @rand (d::Die) -> 7
        @test rand(d) == 7

        # redefinition
        @rand rand(d::Die) = rand(1:d.n)
        @test rand(d) ∈ 1:6
        @test rand(rng0, d) ∈ 1:6
        @test all(∈(1:6), rand(rng0, d, 10))
        @test eltype(rand(d, 3)) == Int

        @rand function (d::Die) rand(1:d.n) end
        @test rand(d) ∈ 1:6
        @test rand(rng0, d) ∈ 1:6
        @test all(∈(1:6), rand(rng0, d, 10))
        @test eltype(rand(d, 3)) == Int

        @rand (d::Die) -> rand(1:d.n)
        @test rand(d) ∈ 1:6
        @test rand(rng0, d) ∈ 1:6
        @test all(∈(1:6), rand(rng0, d, 10))
        @test eltype(rand(d, 3)) == Int

        # redefinition (multiple inner samplers)
        @rand rand(d::Die) = rand(10:10) + rand(1:d.n)
        @test rand(d) ∈ 11:16
        @test all(∈(11:16), rand(d, 10))

        @rand function (d::Die)
            rand(10:10) + rand(1:d.n)
        end
        @test rand(d) ∈ 11:16
        @test all(∈(11:16), rand(d, 10))

        @rand (d::Die) -> rand(10:10) + rand(1:d.n)
        @test rand(d) ∈ 11:16
        @test all(∈(11:16), rand(d, 10))

        # redefinition access to argument not within rand call
        @rand rand(d::Die) = rand(1:1) + d.n
        @test all(==(7), rand(d, 10))

        # redefinition back to SamplerTrivial
        @rand rand(d::Die) = d.n
        @test all(==(6), rand(d, 100))

        # redefinition (Val(Inf) iff rand calls with 2+ arguments)
        @rand rand(d::Die) = (rand("asd") + rand("asd", 3)[1]; 1)
        s = Sampler(MersenneTwister, d, Val(1))
        @test s.data[1] isa Random.SamplerSimple{String}
        @test s.data[2] isa Random.SamplerSimple{Vector{Char}} # "proof" of Val(Inf)
        s = Sampler(MersenneTwister, d, Val(Inf))
        @test s.data[1] isa Random.SamplerSimple{Vector{Char}}
        @test s.data[2] isa Random.SamplerSimple{Vector{Char}}

        # test esc-correctness
        VAR = 100
        @rand rand(d::Die) = rand(VAR+1:VAR+d.n) - VAR
        @test all(∈(1:6), rand(d, 100))

        # with type parameters
        d = DieT(6)

        @rand rand(d::DieT{T}) where {T} = 1
        @test rand(d) == 1

        @rand function rand(d::DieT{T}) where {T}
            2
        end
        @test rand(d) == 2

        @rand function (d::DieT{T}) where {T}
            3
        end
        @test rand(d) == 3

        @rand ((d::DieT{T}) where {T}) -> 4
        @test rand(d) == 4

        @rand rand(d::DieT{Int}) = 0
        @test rand(d) == 0

        d = DieT(0x0)
        @rand rand(d::DieT{T}) where {T<:UInt8} = rand(typemin(T):typemax(T))
        @test rand(d) isa UInt8

        d = DieT(true)
        @rand rand(d::DieT{T}) where {T<:Bool} =
            T(typemin(T) + rand(typemin(T):typemax(T)))
        @test rand(d) isa Bool
    end
end

module TestAtRand
# only using RandomExtensions, not Random, to check we don't depend on
# some imported names like e.g. SamplerTrivial
using RandomExtensions
using ..RandomExtensionsTests: @test, @testset

struct Die n end

@testset "@rand in module" begin
    @rand rand(d::Die) = d.n # SamplerTrivial
    @test rand(Die(123)) == 123
    @test rand(RandomDevice(), Die(123)) == 123
    @rand rand(d::Die) = rand(1:d.n) # SamplerSimple
    @test rand(Die(1)) == 1
    @test rand(RandomDevice(), Die(1)) == 1
end

end # module TestAtRand


text:
using Lint: StaticTypeAnalysis

@testset "Static Type" begin

@test StaticTypeAnalysis.canequal(Int, Int)
@test StaticTypeAnalysis.canequal(Int, Float64)
@test StaticTypeAnalysis.canequal(Int, Real)
@test StaticTypeAnalysis.canequal(Symbol, Symbol)
@test !StaticTypeAnalysis.canequal(Int, Symbol)
@test !StaticTypeAnalysis.canequal(Int, String)
@test !StaticTypeAnalysis.canequal(Vector{Char}, String)

@test StaticTypeAnalysis.infertype(String, (Int,)) == Union{}
@test StaticTypeAnalysis.infertype(String, (String,)) == String
@test StaticTypeAnalysis.infertype(String, (Vector{UInt8},)) == String

@test StaticTypeAnalysis.eltype(Tuple{Int,Int}) == Int
@test StaticTypeAnalysis.eltype(Tuple{Int,String}) == Any
@test StaticTypeAnalysis.eltype(Tuple{}) == Union{}
@test StaticTypeAnalysis.eltype(Tuple) == Any
@test StaticTypeAnalysis.eltype(Vector{Int}) == Int
@test StaticTypeAnalysis.eltype(Array{Int}) == Int

@test StaticTypeAnalysis.length(Tuple) === nothing
@test StaticTypeAnalysis.length(Array) === nothing
@test StaticTypeAnalysis.length(Vector{Int}) === nothing
@test StaticTypeAnalysis.length(Array{Int}) === nothing
@test StaticTypeAnalysis.length(Pair) == 2
@test StaticTypeAnalysis.length(Tuple{Int,String}) == 2
@test StaticTypeAnalysis.length(NTuple{10, Integer}) == 10

@test StaticTypeAnalysis.typeof_nth(Tuple{Int,String}, 1) == Int
@test StaticTypeAnalysis.typeof_nth(Tuple{Int,String}, 2) == String
@test StaticTypeAnalysis.typeof_nth(Pair{Int,String}, 2) == String
@test StaticTypeAnalysis.typeof_nth(Pair{Int,String}, 1) == Int
@test StaticTypeAnalysis.typeof_nth(Vector{Int}, 1) == Int
@test StaticTypeAnalysis.typeof_nth(Tuple, 1) == Any
@test StaticTypeAnalysis.typeof_nth(Tuple{Int}, 2) == Union{}

end


text:
perms = [1 2 3 4
         3 1 2 4
         2 3 1 4
         1 4 2 3
         2 1 4 3
         4 2 1 3
         3 4 1 2
         1 3 4 2
         4 1 3 2
         4 3 2 1
         2 4 3 1
         3 2 4 1];
nodes = [-1 -1 -1
          1 -1 -1
         -1  1 -1
         -1 -1  1
          1  1  1];
elems = [1 2 3 4
         3 4 2 5];

Nperms = size(perms)[1]

println("\$MeshFormat")
println("2.2 0 8")
println("\$EndMeshFormat")
println("\$Nodes")
println(Nperms^2 * size(nodes)[1])
for p1 = 1:Nperms
  for p2 = 1:Nperms
    p = (p1-1) * Nperms + (p2-1)
    for r = 1:size(nodes)[1]
      @printf("%3d", r + p * size(nodes)[1])
      for c = 1:size(nodes)[2]
        @printf(" %3d", nodes[r,c] + 3*p)
      end
      @printf("\n")
    end
  end
end
println("\$EndNodes")
println("\$Elements")
println(Nperms^2 * size(elems)[1])
for p1 = 1:Nperms
  for p2 = 1:Nperms
    p = (p1-1) * Nperms + (p2-1)
    for r = 1:size(elems)[1]
      pp = (r % 2 == 1)  ? p1 : p2
      @printf("%3d 4 2 0 1", r + p * size(elems)[1])
      for c = 1:size(elems)[2]
        @printf(" %3d", elems[r,perms[pp,c]] + p * size(nodes)[1])
      end
      @printf("\n")
    end
  end
end
println("\$EndElements")


text:
@doc """
Enumerates all connected neutral networks for NKqLandscapes and all
approximately neutral networks for NKLandscapes, and summarizes their
properties.  The objective is to see if there are "giant" connected neutral
networks that contain almost all genotypes at a given fitness level as
predicted by the "percolation" model of Sergey Gavrilets for uncorrelated
landscapes.  (Reference needed.) An integer representation of genotypes
is used extensively.  The functions `gtoi` and
`itog` convert between the integer representation and the
array representation used elsewhere in NKLandscapes.jl.  If the arity
of the landscape is 2, then the integer representation can be
considered as a bit-string representation, and some functions require
that the arity of the landscape is 2.
"""
using NKLandscapes
using DataStructures

@doc """ function neutral_net_summary(counts, num_intervals)
Returns a list of information about the connected neutral nets at each non-empty fitness level
Each list element is a 5-tuple with the following components
*  The fit level (empty fit levels are not included)
*  The lower bound of the fit level
*  The number of genotypes at the fit level
*  The maximum size of a neutral net at the fit level
*  The ratio of maximum size neutral net over the number of genotypes
Note that if this ratio is close to 1.0, then a high proportion of the genotypes at the fitness level
  are in a "giant component" neutral net.
"""
function neutral_net_summary(counts, num_intervals)
  fit_increment = 1.0/num_intervals
  nn3 = map(x->x[3],counts)
  lb = minimum(nn3)
  ub = maximum(nn3)
  summary = []
  for i = lb:ub
    filtered_nn = [y[2] for y in filter(x->x[3]==i,counts)]
    if length(filtered_nn) > 0
      max_nn = maximum(filtered_nn)
      sum_nn = sum(filtered_nn)
      push!(summary,(i, i*fit_increment, sum_nn, max_nn, max_nn/sum_nn))
    end
  end
  return summary
end

@doc """ function neutral_net_summary(ls::NKLandscape)

Produces the same results as the previous function, but starting with an NK or NKq landscape.
"""

function neutral_net_summary(ls::Landscape, num_intervals::Int64=ls.n)
  lsf = lsfits(ls);
  flevs = fitlevs(ls, num_intervals, lsf);
  nn = neutralnets(ls,flevs)
  ncounts = netcounts(nn,flevs)
  neutral_net_summary(ncounts,num_intervals)
end

@doc """
Prints the the list of fitness levels returned by neutral_net_summary in CSV format.
"""
function print_neutral_nets(summary)
  @printf("fit lev\tfit_lb\t total count\t   max count\tpct max\n")
  for s in summary
    @printf("%2d\t%.3f\t%12d\t%12d\t%.5f\n",s[1],s[2],convert(Int,s[3]),convert(Int,s[4]),s[5])
  end
end