text:
@testset "reinterpret" begin
    # Gray
    for sz in ((4,), (4,5))
        a = rand(Gray{N0f8}, sz)
        for T in (Gray{N0f8}, Gray{Float32}, Gray{Float64})
            b = @inferred(convert(Array{T}, a))
            rb = @inferred(reinterpretc(eltype(T), b))
            @test eltype(rb) == eltype(T) && ndims(rb) == length(sz)
            @test size(rb) == sz
            c = copy(rb)
            rc = @inferred(reinterpretc(T, c))
            @test eltype(rc) == T && ndims(rc) == length(sz)
            @test size(rc) == sz
        end
    end
    for sz in ((4,), (4,5))
        # Bool/Gray{Bool}
        b = rand(Bool, sz)
        rb = @inferred(reinterpretc(Gray{Bool}, b))
        @test eltype(rb) == Gray{Bool} && ndims(rb) == length(sz)
        @test size(rb) == sz
        c = copy(rb)
        rc = @inferred(reinterpretc(Bool, c))
        @test eltype(rc) == Bool && ndims(rc) == length(sz)
        @test size(rc) == sz
    end
    for sz in ((4,), (4,5))
        b = Gray24.(reinterpret(N0f8, rand(UInt8, sz)))
        for T in (UInt32, RGB24)
            rb = @inferred(reinterpretc(T, b))
            @test eltype(rb) == T && ndims(rb) == length(sz)
            @test size(rb) == sz
            c = copy(rb)
            rc = @inferred(reinterpretc(Gray24, c))
            @test eltype(rc) == Gray24 && ndims(rc) == length(sz)
            @test size(rc) == sz
        end
    end
    # TransparentGray
    a = rand(AGray{N0f8}, (4,5))
    for T in (AGray{N0f8}, GrayA{Float32}, AGray{Float64})
        b = @inferred(convert(Array{T}, a))
        rb = @inferred(reinterpretc(eltype(T), b))
        @test eltype(rb) == eltype(T) && ndims(rb) == 3
        @test size(rb) == (2,4,5)
        c = copy(rb)
        rc = @inferred(reinterpretc(T, c))
        @test eltype(rc) == T && ndims(rc) == 2
        @test size(rc) == (4,5)
    end
    # Color3
    a = rand(RGB{N0f8}, (4,5))
    for T in (RGB{N0f8}, HSV{Float32}, XYZ{Float64})
        b = @inferred(convert(Array{T}, a))
        rb = @inferred(reinterpretc(eltype(T), b))
        @test eltype(rb) == eltype(T) && ndims(rb) == 3
        @test size(rb) == (3,4,5)
        c = copy(rb)
        rc = @inferred(reinterpretc(T, c))
        @test eltype(rc) == T && ndims(rc) == 2
        @test size(rc) == (4,5)
    end
    for a in (rand(RGB{N0f8}, 4), rand(RGB{N0f8}, (4,5)))
        b = @inferred(reinterpretc(HSV{Float32}, float32.(a)))
        @test eltype(b) == HSV{Float32}
        @test ndims(b) == ndims(a)
    end
    # Transparent color
    a = rand(ARGB{N0f8}, (4,5))
    for T in (ARGB{N0f8}, AHSV{Float32}, AXYZ{Float64})
        b = @inferred(convert(Array{T}, a))
        rb = @inferred(reinterpretc(eltype(T), b))
        @test eltype(rb) == eltype(T) && ndims(rb) == 3
        @test size(rb) == (4,4,5)
        c = copy(rb)
        rc = @inferred(reinterpretc(T, c))
        @test eltype(rc) == T && ndims(rc) == 2
        @test size(rc) == (4,5)
    end
    # XRGB/RGBX
    a = rand(RGB{N0f8}, (4,5))
    for T in (XRGB{N0f8},RGBX{Float32})
        b = @inferred(convert(Array{T}, a))
        rb = @inferred(reinterpretc(eltype(T), b))
        @test eltype(rb) == eltype(T) && ndims(rb) == 3
        @test size(rb) == (4,4,5)
        c = copy(rb)
        rc = @inferred(reinterpretc(T, c))
        @test eltype(rc) == T && ndims(rc) == 2
        @test size(rc) == (4,5)
    end
    a = [RGB(1,0,0) RGB(0,0,1);
         RGB(0,1,0) RGB(1,1,1)]
    @test @inferred(reinterpretc(N0f8, a)) == cat([1 0; 0 1; 0 0], [0 1; 0 1; 1 1]; dims=3)
    b = convert(Array{BGR{N0f8}}, a)
    @test @inferred(reinterpretc(N0f8, b)) == cat([0 0; 0 1; 1 0], [1 1; 0 1; 0 1]; dims=3)
    # RGB24, ARGB32
    for sz in ((4,), (4,5))
        a = rand(UInt32, sz)
        for T in (RGB24, ARGB32)
            b = @inferred(reinterpretc(T, a))
            @test eltype(b) == T && ndims(b) == length(sz)
            @test size(b) == sz
            @test eltype(b) == T
            @test @inferred(reinterpret(UInt32, b)) == a
        end
    end

    # 1d
    a = RGB{Float64}[RGB(1,1,0)]
    af = @inferred(reinterpretc(Float64, a))
    anew = @inferred(reinterpretc(RGB, vec(af)))
    @test anew[1] == a[1]
    @test ndims(anew) == 0

    # #33 and its converse
    a = reinterpretc(BGRA{N0f8}, [0xf0884422])
    @test isa(a, AbstractVector) && a == [BGRA{N0f8}(0.533,0.267,0.133,0.941)]
    @test reinterpretc(UInt32, a) == [0xf0884422]
    @test size(reinterpretc(BGRA{N0f8}, rand(UInt32, 5, 5))) == (5,5)
    @test size(color_view(ARGB32, rand(BGRA{N0f8}, 5, 5))) == (5,5)
    a = reinterpretc(BGRA{N0f8}, reshape([0x22, 0x44, 0x88, 0xf0, 0x01, 0x02, 0x03, 0x04], 4, 2))
    @test a == [BGRA{N0f8}(0.533,0.267,0.133,0.941), BGRA{N0f8}(0.012, 0.008, 0.004, 0.016)]
    @test reinterpret(UInt8, a) == [0x22, 0x44, 0x88, 0xf0, 0x01, 0x02, 0x03, 0x04]
    @test color_view(ARGB32, a) == reinterpretc(ARGB32, [0xf0884422,0x04030201])

    # indeterminate type tests
    a = Array{RGB{AbstractFloat}}(undef, 3)
    @test_throws ErrorException reinterpretc(Float64, a)
    a = Vector{RGB}(undef, 3)
    @test_throws ErrorException reinterpretc(Float64, a)

    # Invalid conversions
    a = rand(UInt8, 4,5)
    ret = @test_throws TypeError reinterpretc(Gray, a)
    a = rand(Int8, 4,5)
    ret = @test_throws TypeError reinterpretc(Gray, a)
end

@testset "convert" begin
    a = [RGB(1,0,0) RGB(0,0,1);
         RGB(0,1,0) RGB(1,1,1)]
    c = @inferred(convert(Array{BGR}, a))
    @test eltype(c) == BGR{N0f8}
    c = @inferred(convert(Array{BGR{Float32}}, a))
    @test eltype(c) == BGR{Float32}
    c = @inferred(convert(Array{Lab}, a))
    @test eltype(c) == Lab{Float32}
    for a in (rand(Float32, (4,5)),
              bitrand(4,5))
        b = @inferred(convert(Array{Gray}, a))
        @test eltype(b) == Gray{eltype(a)}
        b = @inferred(convert(Array{Gray{N0f8}}, a))
        @test eltype(b) == Gray{N0f8}
    end

    # Gray images wrapped by an OffsetArray.
    A = rand(8,8)
    for img in ( Gray.(A),
                 Gray.(N0f8.(A)),
                 Gray.(N0f16.(A)) )
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{Gray{Float32},2,Array{Gray{Float32},2}}, imgo))
        @test eltype(s) == Gray{Float32}
        @test s isa OffsetArray{Gray{Float32},2,Array{Gray{Float32},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end

    for img in ( Gray.(A),
                 Gray.(N0f8.(A)),
                 Gray.(N0f16.(A)) )
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{Gray{N0f8},2,Array{Gray{N0f8},2}}, imgo))
        @test eltype(s) == Gray{N0f8}
        @test s isa OffsetArray{Gray{N0f8},2,Array{Gray{N0f8},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end

    for img in ( Gray.(A),
                 Gray.(N0f8.(A)),
                 Gray.(N0f16.(A)) )
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{Gray{N0f16},2,Array{Gray{N0f16},2}},imgo))
        @test eltype(s) == Gray{N0f16}
        @test s isa OffsetArray{Gray{N0f16},2,Array{Gray{N0f16},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end

    # Color images wrapped by an OffsetArray.
    A = rand(RGB{Float32},8,8)
    for img in ( A,
                 n0f8.(A),
                 n6f10.(A),
                 n4f12.(A),
                 n2f14.(A),
                 n0f16.(A))
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{RGB{N0f8},2,Array{RGB{N0f8},2}},imgo))
        @test eltype(s) == RGB{N0f8}
        @test s isa OffsetArray{RGB{N0f8},2,Array{RGB{N0f8},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end

    A = rand(RGB{Float32},8,8)
    for img in ( A,
                 n0f8.(A),
                 n6f10.(A),
                 n4f12.(A),
                 n2f14.(A),
                 n0f16.(A))
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{RGB{Float32},2,Array{RGB{Float32},2}},imgo))
        @test eltype(s) == RGB{Float32}
        @test s isa OffsetArray{RGB{Float32},2,Array{RGB{Float32},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end
end

@testset "eltype conversion" begin
    @test float32(Float64) == Float32
    @test float32(N0f8)      == Float32
    @test float64(RGB{N0f8}) == RGB{Float64}

    a = [RGB(1,0,0) RGB(0,0,1);
         RGB(0,1,0) RGB(1,1,1)]
    @test eltype(a) == RGB{N0f8}
    @test eltype(n0f8.(a))       == RGB{N0f8}
    @test eltype(n6f10.(a)) == RGB{N6f10}
    @test eltype(n4f12.(a)) == RGB{N4f12}
    @test eltype(n2f14.(a)) == RGB{N2f14}
    @test eltype(n0f16.(a)) == RGB{N0f16}
#    @test eltype(float16.(a)) == RGB{Float16}
    @test eltype(float32.(a)) == RGB{Float32}
    @test eltype(float64.(a)) == RGB{Float64}

    a = N0f8[0.1,0.2,0.3]
    @test eltype(a) == N0f8
    @test eltype(n0f8.(a))       == N0f8
    @test eltype(n6f10.(a)) == N6f10
    @test eltype(n4f12.(a)) == N4f12
    @test eltype(n2f14.(a)) == N2f14
    @test eltype(n0f16.(a)) == N0f16
#    @test eltype(float16.(a)) == Float16
    @test eltype(float32.(a)) == Float32
    @test eltype(float64.(a)) == Float64

    a = OffsetArray(N0f8[0.1,0.2,0.3], -1:1)
    @test eltype(a) == N0f8
    @test eltype(n0f8.(a))       == N0f8
    @test eltype(n6f10.(a)) == N6f10
    @test eltype(n4f12.(a)) == N4f12
    @test eltype(n2f14.(a)) == N2f14
    @test eltype(n0f16.(a)) == N0f16
#    @test eltype(float16.(a)) == Float16
    @test eltype(float32.(a)) == Float32
    @test eltype(float64.(a)) == Float64
    @test axes(float32.(a)) == (-1:1,)

    @testset "float" begin
        # float(::Type) is equivalent to floattype
        @test @inferred(float(RGBA{Float32})) == RGBA{Float32}
        @test @inferred(float(BGR{N0f8})    ) == BGR{Float32}
        @test @inferred(float(Gray{N0f8})   ) == Gray{Float32}

        # `floattype(::Number)` is inconsistent to `float(::NUmber)`, see issue:
        # https://github.com/JuliaMath/FixedPointNumbers.jl/issues/127
        # @test @inferred(float(N0f8)         ) == Float32
        # @test @inferred(float(Bool)         ) == Float32
        # @test @inferred(float(Float32)      ) == Float32
        # @test @inferred(float(Float64)      ) == Float64

        @test @inferred(float(ARGB32)) == ARGB{Float32}
        @test @inferred(float(AGray32)) == AGray{Float32}
        @test @inferred(float(RGB24)) == RGB{Float32}
        @test @inferred(float(Gray24)) == Gray{Float32}

        # float(x::Colorant)
        @test float(oneunit(Gray{N0f8})) == oneunit(Gray{Float32})
        @test float(RGB(oneunit(Gray{N0f8}))) == RGB(oneunit(Gray{Float32}))
    end
end

@testset "ambiguities" begin
    # issue #40
    d = Dict(:a=>1, :b=>2.0)
    @test isa(Dict(k=>v for (k,v) in d), Dict{Symbol,Real})
end

nothing


text:
using StrideArrays, VectorizationBase, ProgressMeter
using StrideArrays: StaticFloat
function jmultpackab!(
  C,
  A,
  B,
  ::Val{W₁},
  ::Val{W₂},
  ::Val{R₁},
  ::Val{R₂},
) where {W₁,W₂,R₁,R₂}
  M, N = size(C)
  K = size(B, 1)
  zc, za, zb = StrideArrays.zstridedpointer.((C, A, B))
  @elapsed(
    StrideArrays.jmultpackAB!(
      zc,
      za,
      zb,
      StaticInt{1}(),
      StaticInt{0}(),
      M,
      K,
      N,
      VectorizationBase.NUM_CORES,
      StaticFloat{W₁}(),
      StaticFloat{W₂}(),
      StaticFloat{R₁}(),
      StaticFloat{R₂}(),
    )
  )
end

function bench_size(
  Cs,
  As,
  Bs,
  ::Val{W₁},
  ::Val{W₂},
  ::Val{R₁},
  ::Val{R₂},
) where {W₁,W₂,R₁,R₂}
  if length(first(Cs)) < length(last(Cs))
    jmultpackab!(
      first(Cs),
      first(As),
      first(Bs),
      Val{W₁}(),
      Val{W₂}(),
      Val{R₁}(),
      Val{R₂}(),
    )
  else
    jmultpackab!(last(Cs), last(As), last(Bs), Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
  end
  gflop = 0.0
  for (C, A, B) ∈ zip(Cs, As, Bs)
    M, K, N = StrideArrays.matmul_sizes(C, A, B)
    # sleep(0.5)
    t = jmultpackab!(C, A, B, Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
    gf = 2e-9M * K * N / t
    gflop += gf
  end
  gflop / length(As)
end
matrix_sizes(s::Int) = (s, s, s)
matrix_sizes(MKN::NTuple{3,Int}) = MKN
size_range(l, u, len) = round.(Int, exp.(range(log(l), stop = log(u), length = len)))
function matrix_range(l, u, len, ::Type{T} = Float64) where {T}
  matrix_range(size_range(l, u, len), T)
end
function matrix_range(S, ::Type{T} = Float64) where {T}
  Alen = 0
  Blen = 0
  Clen = 0
  for s ∈ S
    M, K, N = matrix_sizes(s)
    Alen = max(Alen, M * K)
    Blen = max(Blen, K * N)
    Clen = max(Clen, M * N)
  end
  Abuf = rand(T, Alen)
  Bbuf = rand(T, Blen)
  Cbuf = rand(T, Clen)
  As = Vector{
    Base.ReshapedArray{T,2,SubArray{T,1,Vector{T},Tuple{Base.OneTo{Int}},true},Tuple{}},
  }(
    undef,
    length(S),
  )
  Bs = similar(As)
  Cs = similar(As)
  for (i, s) ∈ enumerate(S)
    M, K, N = matrix_sizes(s)
    As[i] = reshape(view(Abuf, Base.OneTo(M * K)), (M, K))
    Bs[i] = reshape(view(Bbuf, Base.OneTo(K * N)), (K, N))
    Cs[i] = reshape(view(Cbuf, Base.OneTo(M * N)), (M, N))
  end
  Cs, As, Bs
end
function gflop_map(
  Cs,
  As,
  Bs,
  ::Val{W₁},
  ::Val{W₂},
  ::Val{R₁},
  ::Val{R₂},
) where {W₁,W₂,R₁,R₂}
  if length(first(Cs)) < length(last(Cs))
    jmultpackab!(
      first(Cs),
      first(As),
      first(Bs),
      Val{W₁}(),
      Val{W₂}(),
      Val{R₁}(),
      Val{R₂}(),
    )
  else
    jmultpackab!(last(Cs), last(As), last(Bs), Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
  end
  gflops = Vector{Float64}(undef, length(Cs))
  for (i, C, A, B) ∈ zip(eachindex(gflops), Cs, As, Bs)
    M, K, N = StrideArrays.matmul_sizes(C, A, B)
    t = jmultpackab!(C, A, B, Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
    gflops[i] = 2e-9M * K * N / t
  end
  gflops
end

# function gridsearch(
#     CsAsBs = matrix_range(1_500, 10_000, 100), w₁range = 0.012:0.0005:0.013, w₂range = 0.024:0.001:0.026, r₁range = 0.44:0.01:0.5, r₂range = 0.75:0.025:0.8
# )
#     search_space = Iterators.product(w₁range, w₂range, r₁range, r₂range)
#     best = Ref(((0.0,0.0),(0.0,0.0),-Inf))
#     gflop_array = let (Cs,As,Bs) = CsAsBs, iter_prod = search_space, p = Progress(length(iter_prod)), best = best
#         map(iter_prod) do (W₁, W₂, R₁, R₂)
#             gflops = bench_size(Cs, As, Bs, Val(W₁), Val(W₂), Val(R₁), Val(R₂))
#             b = best[]
#             recent = ((M_c, K_c, N_c), (R₂, R₃), gflops)
#             bb = if last(b) > gflops
#                 b
#             else
#                 best[] = recent
#             end
#             ProgressMeter.next!(p, showvalues = [(:Last, recent), (:Best, bb)])
#             gflops
#         end
#     end
#     gflop_array, best
# end

# search_range = (120:24:120, 700:200:1100, 4000:2000:6000, 0.44:0.01:0.47, 0.80:0.025:0.875)
# gflop_array, best = search(search_range...); getindex.(search_range, Tuple(argmax(gflop_array)))

# S, Cs, As, Bs = matrix_range(1500, 10_000, 100);

# search_range = (120:24:120, 1000:100:1100, 5000:1000:5000, 0.44:0.01:0.47, 0.80:0.025:0.825)
# gflop_array, best = search((S,Cs,As,Bs), search_range...); getindex.(search_range, Tuple(argmax(gflop_array)))

# gflops_range = gflop_map(CsConst, AsConst, BsConst, Val{StrideArrays.W₁Default}(), Val{StrideArrays.W₂Default}(), Val{StrideArrays.R₁Default}(), Val{StrideArrays.R₂Default}());

# using StatsBase, UnicodePlots

# summarystats(gflops_range)
# lineplot(SR, gflops_range, title = "Square Matrix GFLOPS", xlabel = "Size", ylabel = "GFLOPS")

# S[10]
# StatsBase.summarystats(@view(gflops_range[30:end]))
# findmin(gflops_range)

# julia> search_range = (120:24:120, 700:200:1100, 4000:2000:6000, 0.40:0.01:0.46, 0.85:0.025:1)
# (120:24:120, 700:200:1100, 4000:2000:6000, 0.4:0.01:0.46, 0.85:0.025:1.0)

# julia> getindex.(search_range, Tuple(argmax(gflop_array)))
# (120, 900, 4000, 0.45, 0.875)

# julia> size(gflop_array)
# (1, 3, 2, 7, 7)

# julia> argmax(gflop_array)
# CartesianIndex(1, 2, 1, 6, 2)

# function genmats(N)
#     A = rand(N,N)
#     B = rand(N,N)
#     C = similar(A); p = StrideArrays.zstridedpointer
#     C, A, B, p(C), p(A), p(B)
# end
# C10_000,A10_000,B10_000,zc10_000,za10_000,zb10_000 = genmats(10_000);
# @time(StrideArrays.jmultpackAB!(zc10_000, za10_000, zb10_000, StaticInt{1}(), StaticInt{0}(), Val(72), Val(1.875), Val(1.05), 10_000,10_000,10_000, StrideArrays.CloseOpen(0, 18), Val(1)))

# search_range = (72:24:120, 1.80:0.025:1.9, 1.00:0.25:1.1);
# search_range = (96:24:144, 1.6:0.05:1.9, 1.0:0.05:1.2);

# julia> search_range = (72:24:120, 1.85:0.025:1.95, 1.1:0.25:1.2);

#   Last:  (120, 583, 5562, (1.875, 1.0), 1553.379336642954)
#   Best:  (120, 607, 5346, (1.8, 1.0), 1559.8622396497135)
# (120, 1.8, 1.0)

# julia> gflop_array, best = search(search_range...); getindex.(search_range, Tuple(argmax(gflop_array)))
# Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:05:24
#   Last:  (120, 560, 5265, 1514.4606921365082)
#   Best:  (96, 738, 3996, 1560.6626281274844)
# (96, 1.85, 1.1)

# julia> best_gflop_array, best_search_range = gflop_array, search_range;

# julia> search_range = (72:24:120, 1.80:0.025:1.9, 1.05:0.25:1.1);

# julia> gflop_array, best = search(search_range...); getindex.(search_range, Tuple(argmax(gflop_array)))
# Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:23:06
#   Last:  (120, 557, 5301, (0.51, 0.91), 1468.106008702092)
#   Best:  (72, 819, 3600, (0.45, 0.91), 1563.9777269382319)
# (72, 0.45, 0.91)

# 72 with (0.4-0.44) x 0.9-9.92 looks good
# julia> permutedims(gflop_array, (2,3,1))
# 7×5×2 Array{Float64, 3}:
# [:, :, 1] =   # (96 x 0.40:0.01:0.46 x 0.9:0.005:0.92)
#  1566.74  1512.84  1521.73  1558.3   1559.93
#  1569.36  1563.85  1563.67  1559.56  1563.38
#  1566.34  1551.83  1562.32  1558.95  1561.91
#  1567.51  1561.19  1560.29  1557.78  1564.33
#  1567.02  1560.08  1547.32  1559.42  1565.47
#  1501.42  1550.91  1484.06  1554.05  1553.47
#  1503.18  1545.86  1549.18  1550.41  1548.85

# [:, :, 2] =   # (72 x 0.40:0.01:0.46 x 0.9:0.005:0.92)
#  1531.42  1521.71  1520.14  1520.35  1520.72
#  1539.86  1536.97  1532.17  1532.33  1536.38
#  1543.98  1541.24  1535.58  1540.67  1536.43
#  1550.4   1543.33  1542.23  1540.41  1548.76
#  1552.08  1544.08  1548.16  1547.86  1546.74
#  1548.65  1487.99  1544.13  1542.94  1541.12
#  1476.7   1512.51  1542.46  1537.65  1535.09

# using BlackBoxOptim

# # Noisy test problem
# using Random
# function randrosenbrock(x)
#   return sum( 100*( x[2:end] .- x[1:end-1].^2 ).^2 .+ ( x[1:end-1] .- 1 ).^2 ) + randexp()
# end
# # # res = compare_optimizers(randrosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 4, MaxFuncEvals = 10_000);
# # # :adaptive_de_rand_1_bin_radiuslimited performed well here and is recomended in `BlackBoxOptim`'s README
# bboptimize(
#     randrosenbrock;
#     SearchRange = fill((-5.0,5.0), 4),
#     Method = :adaptive_de_rand_1_bin_radiuslimited,
#     MaxFuncEvals = 10_000
# )

# using Optim
# Optim.optimize(randrosenbrock, fill(-5.0, 4), fill(5.0, 4), rand(4) .* 10 .- 5, SAMIN(), Optim.Options(iterations = 10_000))
# Optim.optimize(randrosenbrock, rand(4) .* 10 .- 5, ParticleSwarm(lower=fill(-5.0, 4),upper=fill(5.0, 4)), Optim.Options(iterations = 10_000))

T = Float64
min_size = round(Int, sqrt(0.65 * StrideArrays.VectorizationBase.CACHE_SIZE[3] / sizeof(T)))
max_size = round(Int, sqrt(32 * StrideArrays.VectorizationBase.CACHE_SIZE[3] / sizeof(T)))

SR = size_range(max_size, min_size, 100);
const CsConst, AsConst, BsConst = matrix_range(SR, T);

function matmul_objective(params)
  print("Params: ", params, "; ")
  W₁, W₂, R₁, R₂ = params
  # print("(W₁ = $(round(W₁, sigdigits=4)); W₂ = $(round(W₂, sigdigits=4)); R₁ = $(round(R₁, sigdigits=4)); R₂ = $(round(R₂, sigdigits=4))); ")
  gflop = bench_size(CsConst, AsConst, BsConst, Val{W₁}(), Val{W₂}(), Val{R₁}(), Val{R₂}())
  println(gflop)
  -gflop
end

using Optim
hours = 60.0 * 60.0;
days = 24hours;
init = [
  StrideArrays.W₁Default,
  StrideArrays.W₂Default,
  StrideArrays.R₁Default,
  StrideArrays.R₂Default,
]
# init = [0.0125, 0.025, 0.6275, 0.9579]
# init = [0.0060790786747738235, 0.4531988431700635, 0.47560416900859487, 0.6776801310495106]
# opt = Optim.optimize(
#     matmul_objective, [0.001, 0.01, 0.3, 0.4], [0.1, 2.0, 0.9, 0.99], init, SAMIN(),
#     Optim.Options(iterations = 10^6, time_limit = days/8)
# );
opt = Optim.optimize(
  matmul_objective,
  init,
  ParticleSwarm(lower = [0.001, 0.01, 0.3, 0.4], upper = [0.1, 2.0, 0.9, 0.99]),
  Optim.Options(iterations = 10^6, time_limit = days / 8),
);
# (W₁ = 0.0125; W₂ = 0.025; R₁ = 0.6275; R₂ = 0.9579); 1551.5119875598948
# (W₁ = 0.0125; W₂ = 0.025; R₁ = 0.6275; R₂ = 0.8122); 1467.716330273459
# (W₁ = 0.0394; W₂ = 0.025; R₁ = 0.6275; R₂ = 0.9579); 1499.605639963072
# (W₁ = 0.0125; W₂ = 1.187; R₁ = 0.6275; R₂ = 0.9579); 1545.9928365889534

# using BlackBoxOptim

# bbres = bboptimize(
#     matmul_objective;
#     SearchRange = [(0.001, 0.1), (0.01,1.0), (0.25, 0.8), (0.25, 0.99)],
#     Method = :adaptive_de_rand_1_bin_radiuslimited,
#     MaxTime = 60.0 * 60.0 * 4.0 # 4 hours
# )

ext: jl


text:
classicalquantum_insts(complex::Bool, use_EF::Bool) = [
    [(d, complex, use_EF)
    for d in vcat(3, 10:10:40, 50:25:100, 150:50:750)] # includes compile run
    ]

insts = OrderedDict()
insts["nat"] = (nothing, classicalquantum_insts(true, false))
insts["ext"] = (nothing, classicalquantum_insts(true, true))
return (ClassicalQuantum, insts)


text:
@recipe function f(domain::AbstractPlanarDomain, nb::Integer = 100, ni::Integer = 100)
    seriestype := :scatter
    markersize --> 2
    label --> ["Boundary" "Interior"]

    boundary = boundary_points(domain, nb, distribution = :linear)[1]
    interior = interior_points(domain, ni)

    x =
        Float64.(
            hcat(
                [getindex.(boundary, 1); fill(NaN, max(0, ni - nb))],
                [getindex.(interior, 1); fill(NaN, max(0, nb - ni))],
            ),
        )
    y =
        Float64.(
            hcat(
                [getindex.(boundary, 2); fill(NaN, max(0, ni - nb))],
                [getindex.(interior, 2); fill(NaN, max(0, nb - ni))],
            ),
        )

    return x, y
end


text:
# QP Lasso testproblem

using Test, LinearAlgebra, SparseArrays, Random
using COSMO

# generate problem data
rng = Random.MersenneTwister(1313)
n = 15
m = 50 * n

A = sprandn(rng, m, n, 0.5)
vtrue = 1/n * sprand(rng, n, 0.5)
noise = 1/4 * randn(rng, m)
b = Vector(A * vtrue + noise)
λ = 0.2 * norm(A' * b, Inf)

# define lasso problem as QP
A1 = [A zeros(m, n) -Matrix(1.0I, m, m)]
A2 = [-Matrix(1.0I, n, n) Matrix(1.0I, n, n) zeros(n, m);
       Matrix(1.0I, n, n) Matrix(1.0I, n, n) zeros(n, m)]

b1 = -b;
b2 = zeros(2 * n)

P = 2 * Matrix(Diagonal([zeros(2 * n); ones(m)]))# times two to cancel the 1/2 in the obj_val function
q = [zeros(n); λ * ones(n); zeros(m)]

constraint1 = COSMO.Constraint(A1, b1, COSMO.ZeroSet)
constraint2 = COSMO.Constraint(A2, b2, COSMO.Nonnegatives)
constraints = [constraint1; constraint2]

model = COSMO.Model()
assemble!(model, P, q, constraints)

res = COSMO.optimize!(model)

@testset "QP - Lasso" begin
  @test res.status == :Solved
  @test isapprox(res.obj_val, 46.40521553063313, atol=1e-3)
end
nothing


text:
using MatrixHelper
using Base.Test

# write your own tests here
a = rand(21,21)
asp = sortpermrows(a)
as = sortrows(a)

@test a[asp,:] == as
for i in 1:21
	@test selectpermrows(a, i) == asp[i]
end
@test selectpermrows(a, 6:12) == asp[6:12]

asp = sortpermcols(a)
as = sortcols(a)
@test a[:,asp] == as
for i in 1:21
	@test selectpermcols(a, i) == asp[i]
end
@test selectpermcols(a, 6:12) == asp[6:12]


text:
# Autogenerated wrapper script for ripgrep_jll for x86_64-w64-mingw32
export rg

JLLWrappers.@generate_wrapper_header("ripgrep")
JLLWrappers.@declare_executable_product(rg)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        rg,
        "bin\\rg.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


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

using Markdown
using InteractiveUtils

# ╔═╡ b11fdebe-5507-11eb-1001-537a20e51d0d
begin
	let
		env = mktempdir()
		import Pkg
		Pkg.activate(env)
		Pkg.Registry.update()
		Pkg.add("FileIO")
		Pkg.add("JLD2")
		Pkg.add("StatsBase")
		Pkg.add("StatsPlots")
		Pkg.add("SimpleANOVA")
		Pkg.add(path="../..")
	end
	using SMLMAssociationAnalysis_NCB
	using FileIO
	using StatsBase
	using Statistics
	using StatsPlots
	using SimpleANOVA
end

# ╔═╡ 56dcb910-5507-11eb-0b0a-a1c2aef89207
md"""
# Analysis of Experiment 2: p53-MEG3
"""

# ╔═╡ 90d2acb0-5507-11eb-39c4-df93c258dd98
md"""Load packages"""

# ╔═╡ 86502ff0-550d-11eb-1c13-553251c061ee
md"""Load results file and compute normalized values"""

# ╔═╡ 5e3db410-550d-11eb-0331-49a447464b07
begin
	### Load saved data

	nreplicates = 3
	nsamples = 4
	ncells = 10

	outputdir = "../../output"
	datapath = joinpath(outputdir, "results.jld2")

	experimentresults = load(datapath)["experimentresults"]

	medianmeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)
	montecarlomeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)
	positivecontrolmontecarlomeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)
	negativecontrolmontecarlomeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)

	for k ∈ 1:2
		for i ∈ 1:nsamples
			samplemedianresults = Array{Float64,2}(undef, ncells, nreplicates)
			samplelessthan10results = Array{Float64,2}(undef, ncells, nreplicates)
			samplepositivecontrollessthan10results = Array{Float64,2}(undef, ncells, nreplicates)
			samplenegativecontrollessthan10results = Array{Float64,2}(undef, ncells, nreplicates)

			for j ∈ 1:nreplicates
				replicateresults = experimentresults[k][j][i]
				lessthanlimitreplicate = [(x.distances .< 200) .& (x.percentileranks .< 0.1) for x ∈ replicateresults]
				positivecontrollessthanlimitreplicate = [(x.positivecontrol_distances .< 200) .& (x.positivecontrol_percentileranks .< 0.1) for x ∈ replicateresults]
				negativecontrollessthanlimitreplicate = median.([count.([(x.negativecontrol_distances[m] .< 200) .& (x.negativecontrol_percentileranks[m] .< 0.1) for m ∈ eachindex(x.negativecontrol_distances)]) for x ∈ replicateresults])
				lessthan10 = count.(lessthanlimitreplicate) ./ length.(lessthanlimitreplicate)
				positivecontrollessthan10 = count.(positivecontrollessthanlimitreplicate) ./ length.(lessthanlimitreplicate)
				negativecontrollessthan10 = negativecontrollessthanlimitreplicate ./ length.(lessthanlimitreplicate)
				mediandistances = map(x -> x.mediandistance, replicateresults)
				samplemedianresults[:, j] = mediandistances
				samplelessthan10results[:, j] = lessthan10
				samplepositivecontrollessthan10results[:, j] = positivecontrollessthan10
				samplenegativecontrollessthan10results[:, j] = negativecontrollessthan10
			end
			medianmeasurements[:, :, i, k] = samplemedianresults
			montecarlomeasurements[:, :, i, k] = samplelessthan10results
			positivecontrolmontecarlomeasurements[:, :, i, k] = samplepositivecontrollessthan10results
			negativecontrolmontecarlomeasurements[:, :, i, k] = samplenegativecontrollessthan10results
		end
	end

	medianmeasurements = cat(
		cat(medianmeasurements[:, :, 1:2, 1], medianmeasurements[:, :, 3:4, 1], dims = 4),
		cat(medianmeasurements[:, :, 1:2, 2], medianmeasurements[:, :, 3:4, 2], dims = 4),
		dims = 5,
	)
	montecarlomeasurements = cat(
		cat(montecarlomeasurements[:, :, 1:2, 1], montecarlomeasurements[:, :, 3:4, 1], dims = 4),
		cat(montecarlomeasurements[:, :, 1:2, 2], montecarlomeasurements[:, :, 3:4, 2], dims = 4),
		dims = 5,
	)
	positivecontrolmontecarlomeasurements = cat(
		cat(positivecontrolmontecarlomeasurements[:, :, 1:2, 1], positivecontrolmontecarlomeasurements[:, :, 3:4, 1], dims = 4),
		cat(positivecontrolmontecarlomeasurements[:, :, 1:2, 2], positivecontrolmontecarlomeasurements[:, :, 3:4, 2], dims = 4),
		dims = 5,
	)
	negativecontrolmontecarlomeasurements = cat(
		cat(negativecontrolmontecarlomeasurements[:, :, 1:2, 1], negativecontrolmontecarlomeasurements[:, :, 3:4, 1], dims = 4),
		cat(negativecontrolmontecarlomeasurements[:, :, 1:2, 2], negativecontrolmontecarlomeasurements[:, :, 3:4, 2], dims = 4),
		dims = 5,
	)

	normalizedmontecarlomeasurements = (montecarlomeasurements .- negativecontrolmontecarlomeasurements) ./ (positivecontrolmontecarlomeasurements .- negativecontrolmontecarlomeasurements)

	nothing
end

# ╔═╡ 294e6a02-550e-11eb-3450-d9d5989ec08c
md"""
Relevant data are stored in these variables:
- `medianmeasurements` - median distance between the pairs
- `montecarlomeasurements` - fraction of pairs considered associated (p < 0.1, dist < 200 nm)
- `normalizedmontecarlomeasurements` - fraction of pairs considered associated, normalized bassed on simulations of 0% and 100% binding

The data is stored in matrices with these dimensions:
- 1 - cell
- 2 - experimental replicate
- 3 - condition 1 (1 = -Dox, 2 = +Dox)
- 4 - condition 2 (Exp 1: 1 = -Nut, 2 = +Nut; Exp 2: 1 = MEG3, 2 = GAPDH)
- 5 - experiment (1 = p53-Mdm2, 2 = p53-MEG3)
"""

# ╔═╡ 9a94f370-5511-11eb-1f60-7fe83ed11def
md"""
### Planning
Central Limit Theorem - Only 10 measurements per level, but 30 per condition. Might be close enough for a normal distribution.

Group sizes equal

Doxycycline and RNA type are independent.
"""

# ╔═╡ 572ba5f2-5513-11eb-1937-3f03cec8252e
md"""
### Medians
#### Check unusual cases
"""

# ╔═╡ 864b181e-5513-11eb-0304-8983d37a8170
let
	p1 = boxplot(medianmeasurements[:, :, 1, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
	p2 = boxplot(medianmeasurements[:, :, 2, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
	p3 = boxplot(medianmeasurements[:, :, 1, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
	p4 = boxplot(medianmeasurements[:, :, 2, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))

	plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "MDM2-p53")
end

# ╔═╡ 27e24820-5514-11eb-1415-8711a194892d
md"""
Potential outliers in box plots:
"""

# ╔═╡ f0393540-5514-11eb-2f76-3514cc3bc85c
[medianmeasurements[5,2,1,1,2],
 medianmeasurements[[4,7],2,2,1,2],
 medianmeasurements[2,1,1,2,2],
medianmeasurements[6,1,2,2,2]]

# ╔═╡ f02df530-5515-11eb-0d9e-8dbb63ed2023
md"""Evaluating z-scores within each replicate:"""

# ╔═╡ ffa9f6c0-576e-11eb-0f11-3139a34ff859
mediansflat = [medianmeasurements[:, :, 1, 1, 2] medianmeasurements[:, :, 2, 1, 2] medianmeasurements[:, :, 1, 2, 2] medianmeasurements[:, :, 2, 2, 2]]

# ╔═╡ 25253910-5516-11eb-0a7a-1d73ed51d4a1
medianzscores = [zscore(mediansflat[:,i]) for i ∈ axes(mediansflat, 2)]

# ╔═╡ d8c96c5e-5517-11eb-0861-056d5e425f04
md"""
- In [1,1] Replicate 1: All within 1.96
- In [1,1] Replicate 2: All within 1.96
- In [1,1] Replicate 3: All within 1.96
- In [1,2] Replicate 1: 1 value, Cell 5 (10%) just under 2.58. None outside 3.29.
- In [1,2] Replicate 2: 1 value, Cell 4 (10%) outside 1.96. None outside 2.58.
- In [1,2] Replicate 3: All within 1.96
- In [2,1] Replicate 1: 1 value, Cell 2 (10%) outside 1.96 significantly. None outside 2.58.
- In [2,1] Replicate 2: All within 1.96
- In [2,1] Replicate 3: 1 value, Cell 3 (10%) outside 1.96 significantly. None outside 2.58.
- In [2,2] Replicate 1: 1 value, Cell 6 (10%) outside 1.96 significantly. None outside 2.58.
- In [2,2] Replicate 2: 1 value, Cell 7 (10%) outside 1.96. None outside 2.58.
- In [2,2] Replicate 3: 1 value, Cell 8 (10%) outside 1.96. None outside 2.58.
"""

# ╔═╡ 87eb8f70-5518-11eb-1639-f91b4bcb5c9b
md"""
None of these potential outliers are due to data entry error or measurement error. Sampling a different population is possible but unknown. These potential outliers are relatively modest, though two quite severe. Since an ANOVA will be used, must be careful not to delete points from a balanced design. Therefore, no outliers will be removed. The two severe outliers will be Windsorized.
"""

# ╔═╡ 74a25f10-59ab-11eb-3a45-bf3b0554a07a
begin
	medianmeasurementsw = copy(medianmeasurements)
	mediansflatw = copy(mediansflat)
	medianmeasurementsw[5,1,2,1,2] = medianmeasurements[4,1,2,1,2]
	mediansflatw[5,4] = mediansflat[4,4]
	medianmeasurementsw[6,1,2,2,2] = medianmeasurements[8,1,2,2,2]
	mediansflatw[6,10] = mediansflat[8,10]
end;

# ╔═╡ a87aac80-5518-11eb-38d4-5fe90140b820
md"""
#### Check assumption of homoscedasticity
"""

# ╔═╡ f6e23cd0-5518-11eb-226a-c5751eb65b34
md"""
Based on the ZResid/ZPred plot, the data appears mildly heteroscedastic; outliers really stand out.

Conducted Levene's test, a 1-way ANOVA of the absolute deviation between each value and the group mean, by replicate.
"""

# ╔═╡ 0ffdd990-5519-11eb-376b-33be0cc966d7
levene(mediansflatw)

# ╔═╡ 28a8da30-5519-11eb-2fb0-efe693c6e44c
md"""
By Levene's test, the data is very likely heteroscedastic.

Group sizes are equal, but samples are on the smaller side. CIs and significance tests may be affected.
"""

# ╔═╡ 5f8de452-5519-11eb-0a2c-051d3e701dd5
md"""
#### Check assumption of normality
"""

# ╔═╡ c4dbef60-55bd-11eb-1822-49c3473000e8
md"""
There are a couple notable deviations but for the *most* part they appear close to normal. The worst ones had the potential outliers identified previously.
"""

# ╔═╡ cf59a26e-55bd-11eb-1c4b-316a603b9007
medianskewness = [skewness(mediansflat[:,i]) for i ∈ axes(mediansflatw, 2)]

# ╔═╡ 135dbe6e-55be-11eb-3a89-dd3d02c3ad4a
mediankurtosis = [kurtosis(mediansflatw[:,i]) for i ∈ axes(mediansflatw, 2)]

# ╔═╡ 6778a880-55be-11eb-1016-c59335beef0a
md"""
From skewness and kurtosis statistics, all were within 1.96 (5% level).
"""

# ╔═╡ 9bf9aaf0-55be-11eb-17a2-75b3eec550f8
md"""
The primary source of non-normality appear to be outliers, relatively minimal. As a result, I don't intend to do anything special.
"""

# ╔═╡ 265006c0-55c1-11eb-2ba8-3d628d096704
md"""
#### Conduct ANOVA tests
"""

# ╔═╡ c101d900-57ab-11eb-04ab-77cd4e0ff4b3
medianresult = anova(
    permutedims(medianmeasurementsw[:, :, :, :, 2], (1, 3, 2, 4)),
    [fixed, subject, fixed],
    factornames = ["Doxycycline", "Replicate", "RNA"],
)

# ╔═╡ 19934eb0-57ba-11eb-3f45-11c68057c1f7
md"""
#### Check presence of interactions
"""

# ╔═╡ 25bb5660-57ba-11eb-3999-4bbf47d29d8a
plot(medianresult)

# ╔═╡ 342cfff0-57ba-11eb-07d6-2b90717f7c7b
md"""
There seem to be large interactions, but they aren't significant. Since we expect RNAs to be different anyway, will continue to examine ANOVA within each RNA.
"""

# ╔═╡ fd2d4a40-57ab-11eb-2688-53b009a37680
md"""
#### Interpret main effects
RNA identity is not significant and has an insignificant effect size.

Doxycycline is not significant and has an insignificant effect size. As the RNA should be independent, we'll do another ANOVA within each RNA.
"""

# ╔═╡ 5f1ce3b0-55c1-11eb-3d03-b53516cfc976
medianresultMEG3 = anova(
    permutedims(medianmeasurementsw[:, :, :, 1, 2], (1, 3, 2)),
    [fixed, subject],
    factornames = ["Doxycycline", "Replicate"],
)

# ╔═╡ 0fddbb50-5772-11eb-01cc-e57dbec14a6b
medianresultGAPDH = anova(
    permutedims(medianmeasurementsw[:, :, :, 2, 2], (1, 3, 2)),
    [fixed, subject],
    factornames = ["Doxycycline", "Replicate"],
)

# ╔═╡ 2ed8ff90-568e-11eb-10a3-0b97b2425d2f
md"""
#### Interpret main effects of within-RNA ANOVAs

Doxycycline produced an insignificant effect with GAPDH (both in size and probability).

Doxycycline produced a large but not significant effect with MEG3, indicating MEG3 and p53 were brought closer together but not unambiguously. This is probably mostly due to the induced expression of MEG3.
"""

# ╔═╡ 98c74810-6af9-11eb-3988-1f4c77155983
md"""
### Figure
"""

# ╔═╡ 68a010a0-6a3a-11eb-07df-9d60e51f3aab
let
	import StatsPlots.mm
	rnagroups = repeat([1,2], inner = 60)
	doxgroups = repeat([2, 1], inner = 30, outer = 2)
	p53_meg3_median = [medianmeasurementsw[:,:,1,1,2] |> vec; medianmeasurementsw[:,:,2,1,2] |> vec; medianmeasurementsw[:,:,1,2,2] |> vec; medianmeasurementsw[:,:,2,2,2] |> vec]
	groupedboxplot(rnagroups, p53_meg3_median, group = doxgroups, outliers=false,
			label=["- Dox", "+ Dox"],
			guidefontsize=13,
			tickfontsize=12,
			legend=:none,
			top_margin=10mm,
		    left_margin=10mm,
			seriescolor=[:white :lightgray],
			line=(2, 1.0),
			gridopacity=0.3,
			xgrid=:none,
			xaxis=("RNA", (1:2, ["MEG3", "GAPDH"])),
			size=(256,512),
			yaxis=("Median exclusive pairwise distance (nm)", (0,1200)))
	p53_meg3_median_means = dropdims(mean(medianmeasurementsw[:,:,:,:,2], dims=1), dims=1)
	groupeddotplot!([1,1,1,1,1,1,2,2,2,2,2,2], group = [2,2,2,1,1,1,2,2,2,1,1,1], p53_meg3_median_means |> vec, mode = :none, label="", marker=(4, 0.75, :rect, repeat([:orange, :darkblue, :darkred]), stroke(0)))
	groupeddotplot!(rnagroups, p53_meg3_median, group = doxgroups, mode = :density, label="", marker=(2, 0.5, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
end

# ╔═╡ 26985da0-6afa-11eb-2675-69a745a338c4
md"""Create paper-quality version and save it:"""

# ╔═╡ 252c7190-6afa-11eb-00fd-5f527f9fa189
let
	import StatsPlots.mm
	rnagroups = repeat([1,2], inner = 60)
	doxgroups = repeat([2, 1], inner = 30, outer = 2)
	p53_meg3_median = [medianmeasurementsw[:,:,1,1,2] |> vec; medianmeasurementsw[:,:,2,1,2] |> vec; medianmeasurementsw[:,:,1,2,2] |> vec; medianmeasurementsw[:,:,2,2,2] |> vec]
	groupedboxplot(rnagroups, p53_meg3_median, group = doxgroups, outliers=false,
			label=["- Dox", "+ Dox"],
			guidefontsize=52,
			tickfontsize=48,
			legend=:none,
			top_margin=10mm,
		    left_margin=25mm,
			seriescolor=[:white :lightgray],
			line=(8, 1.0),
			gridopacity=0.3,
			xgrid=:none,
			xaxis=("RNA", (1:2, ["MEG3", "GAPDH"])),
			size=(1024,2048),
			yaxis=("Median exclusive pairwise distance (nm)", (0,1200)))
	p53_meg3_median_means = dropdims(mean(medianmeasurementsw[:,:,:,:,2], dims=1), dims=1)
	groupeddotplot!([1,1,1,1,1,1,2,2,2,2,2,2], group = [2,2,2,1,1,1,2,2,2,1,1,1], p53_meg3_median_means |> vec, mode = :none, label="", marker=(12, 0.75, :rect, repeat([:orange, :darkblue, :darkred]), stroke(0)))
	groupeddotplot!(rnagroups, p53_meg3_median, group = doxgroups, mode = :density, label="", marker=(8, 0.5, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
	savefig(joinpath(outputdir, "p53-meg3-medians.png"))
end

# ╔═╡ 25d0ce1e-5691-11eb-2cc8-49980b6b2161
md"""
### Monte Carlo fraction associated
#### Check unusual cases
"""

# ╔═╡ 8c0ebd10-56a9-11eb-3e8c-95b53bf12061
let
	p1 = boxplot(montecarlomeasurements[:, :, 1, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p2 = boxplot(montecarlomeasurements[:, :, 2, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p3 = boxplot(montecarlomeasurements[:, :, 1, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p4 = boxplot(montecarlomeasurements[:, :, 2, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))

	plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "Mdm2-p53")
end

# ╔═╡ faa812c0-56aa-11eb-3b6b-99fe5f200b67
md"""
Potential outliers in box plots:
"""

# ╔═╡ 012a7c50-56ab-11eb-3d50-c9820925f9ac
[montecarlomeasurements[2,3,1,1,2],
 montecarlomeasurements[[7,9],3,1,2,2]]

# ╔═╡ 1637a040-56ac-11eb-2433-49c37383e005
md"""Evaluating z-scores within each replicate:"""

# ╔═╡ 87ec3520-56ac-11eb-2a1a-2981155d4168
montecarloflat = [montecarlomeasurements[:, :, 1, 1, 2] montecarlomeasurements[:, :, 2, 1, 2] montecarlomeasurements[:, :, 1, 2, 2] montecarlomeasurements[:, :, 2, 2, 2]]

# ╔═╡ 9ce1f5f0-56ac-11eb-22a5-d9b5d6a18978
montecarlozscores = [zscore(montecarloflat[:,i]) for i ∈ axes(montecarloflat, 2)]

# ╔═╡ ba7d2760-56ac-11eb-11ca-d543a628c68b
md"""
- In [1,1] Replicate 1: 1 value, Cell 4 (10%) outside 1.96 barely. None outside 2.58.
- In [1,1] Replicate 2: All within 1.96
- In [1,1] Replicate 3: 1 value, Cell 2 (10%) outside 1.96 significantly. None outside 2.58.
- In [1,2] Replicate 1: 1 value, Cell 4 (10%) outside 1.96. None outside 2.58.
- In [1,2] Replicate 2: All within 1.96
- In [1,2] Replicate 3: All within 1.96
- In [2,1] Replicate 1: All within 1.96
- In [2,1] Replicate 2: All within 1.96
- In [2,1] Replicate 3: All within 1.96
- In [2,2] Replicate 1: All within 1.96
- In [2,2] Replicate 2: All within 1.96
- In [2,2] Replicate 3: 1 value, Cell 8 (10%) outside 1.96. None outside 2.58.
"""

# ╔═╡ 03d7ceae-6d78-11eb-3f4b-13cf782576c0
begin
	montecarlomeasurementsw = copy(montecarlomeasurements)
	montecarloflatw = copy(montecarloflat)
	montecarlomeasurementsw[2,3,1,1,2] = montecarlomeasurements[8,3,1,1,2]
	montecarloflatw[2,3] = montecarloflat[8,3]
end;

# ╔═╡ c18cc050-56ad-11eb-261c-4177aa31ca7c
md"""
None of these potential outliers are due to data entry error or measurement error. Sampling a different population is possible but unknown. These potential outliers are relatively modest. In addition, since an ANOVA will be used, must be careful not to delete points from a balanced design. Therefore, no outliers will be removed.

However, I will Windsorize the one severe outlier in [1,1] Replicate 3.
"""

# ╔═╡ 2bba9b02-56ae-11eb-2e30-439dc20273f8
md"""
#### Check assumption of homoscedasticity
"""

# ╔═╡ 0f8a00a0-56af-11eb-2dc6-abd2a45298a4
md"""
Based on the ZResid/ZPred plot, the data appears very heteroscedastic.

Conducted Levene's test, a 1-way ANOVA of the absolute deviation between each value and the group mean, by replicate.
"""

# ╔═╡ b1ae0750-56af-11eb-0714-35b78e197b5e
levene(montecarloflatw)

# ╔═╡ f9676320-56af-11eb-1481-ef9502edded8
md"""
By Levene's test, the data is very likely heteroscedastic.

Group sizes are equal, but samples are on the smaller side. CIs and significance tests may be affected.
"""

# ╔═╡ 4e36f730-56b0-11eb-21e8-0d2e7d0a2562
md"""
#### Check assumption of normality
"""

# ╔═╡ 4423b3d0-56b2-11eb-22be-fd86deabef7f
md"""
There are a couple notable deviations but for the *most* part they appear close to normal.
"""

# ╔═╡ db6922c0-56b2-11eb-2601-5525d3f6bf70
montecarloskewness = [skewness(montecarloflatw[:,i]) for i ∈ axes(montecarloflatw, 2)]

# ╔═╡ 047df370-56b3-11eb-3036-2393031bd967
montecarlokurtosis = [kurtosis(montecarloflatw[:,i]) for i ∈ axes(montecarloflatw, 2)]

# ╔═╡ 3316f790-56b3-11eb-3d9c-e9ebf797c923
md"""
From skewness and kurtosis statistics, none was greater in value than 1.96 (5% level).
"""

# ╔═╡ 6f1da2c0-56b3-11eb-2506-5f9a944735b6
md"""
#### Conduct ANOVA tests
"""

# ╔═╡ c29f8040-57a7-11eb-2cb4-7974bcf74e27
montecarloresult = anova(
    permutedims(montecarlomeasurementsw[:, :, :, :, 2], (1, 3, 2, 4)),
    [fixed, subject, fixed],
    factornames = ["Doxycycline", "Replicate", "RNA"],
)

# ╔═╡ 2553a3a0-57a9-11eb-247d-a161d85d4410
md"""
#### Interpret main effects
RNA identity is significant and has a very large effect size.

Doxycycline is not significant and has a small effect size. As the RNA should be independent, we'll do another ANOVA within each RNA.
"""

# ╔═╡ 755d8510-56b3-11eb-3a0c-53b643ba2053
montecarloresultMEG3 = anova(
    permutedims(montecarlomeasurementsw[:, :, :, 1, 2], (1, 3, 2)),
    [fixed, subject],
    factornames = ["Doxycycline", "Replicate"],
)

# ╔═╡ 28212d70-57a7-11eb-1c94-3fa38021ba8e
montecarloresultGAPDH = anova(
    permutedims(montecarlomeasurementsw[:, :, :, 2, 2], (1, 3, 2)),
    [fixed, subject],
    factornames = ["Doxycycline", "Replicate"],
)

# ╔═╡ 248d1560-5753-11eb-2156-89cfb6742205
md"""
#### Interpret main effects of within-RNA ANOVAs
Doxycycline had a medium effect on MEG3 and p53, although not significant. Doxycyline increased fraction bound.

Doxycyline had no effect on GAPDH and p53 and was highly nonsignificant
"""

# ╔═╡ 174f2590-6a46-11eb-0621-73610de2d1c0
let
	import StatsPlots.mm
	rnagroups = repeat([1,2], inner = 60)
	doxgroups = repeat([2, 1], inner = 30, outer = 2)
	p53_meg3_montecarlo = [montecarlomeasurementsw[:,:,1,1,2] |> vec; montecarlomeasurementsw[:,:,2,1,2] |> vec; montecarlomeasurementsw[:,:,1,2,2] |> vec; montecarlomeasurementsw[:,:,2,2,2] |> vec]
	groupedboxplot(rnagroups, p53_meg3_montecarlo, group = doxgroups, outliers=false,
			label=["- Dox", "+ Dox"],
			guidefontsize=13,
			tickfontsize=12,
			legend=:none,
			top_margin=10mm,
		    left_margin=10mm,
			seriescolor=[:white :lightgray],
			line=(2, 1.0),
			gridopacity=0.3,
			xgrid=:none,
			xaxis=("RNA", (1:2, ["MEG3", "GAPDH"])),
			size=(256,512),
			yaxis=("Fraction associated", (0,0.25)))
	p53_meg3_montecarlo_means = dropdims(mean(montecarlomeasurementsw[:,:,:,:,2], dims=1), dims=1)
	groupeddotplot!([1,1,1,1,1,1,2,2,2,2,2,2], group = [2,2,2,1,1,1,2,2,2,1,1,1], p53_meg3_montecarlo_means |> vec, mode = :none, label="", marker=(4, 0.75, :rect, repeat([:orange, :darkblue, :darkred]), stroke(0)))
	groupeddotplot!(rnagroups, p53_meg3_montecarlo, group = doxgroups, mode = :density, label="", marker=(2, 0.75, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
end

# ╔═╡ 0185bed0-6afb-11eb-279d-b9e008c964c3
md"""Create paper-quality version and save it:"""

# ╔═╡ fecb59c0-6afa-11eb-15fe-cbab4e3f44e1
let
	import StatsPlots.mm
	rnagroups = repeat([1,2], inner = 60)
	doxgroups = repeat([2, 1], inner = 30, outer = 2)
	p53_meg3_montecarlo = [montecarlomeasurementsw[:,:,1,1,2] |> vec; montecarlomeasurementsw[:,:,2,1,2] |> vec; montecarlomeasurementsw[:,:,1,2,2] |> vec; montecarlomeasurementsw[:,:,2,2,2] |> vec]
	groupedboxplot(rnagroups, p53_meg3_montecarlo, group = doxgroups, outliers=false,
			label=["- Dox", "+ Dox"],
			guidefontsize=52,
			tickfontsize=48,
			legend=:none,
		    left_margin=25mm,
		    top_margin=10mm,
			seriescolor=[:white :lightgray],
			line=(8, 1.0),
			gridopacity=0.3,
			xgrid=:none,
			xaxis=("RNA", (1:2, ["MEG3", "GAPDH"])),
			size=(1024,2048),
			yaxis=("Fraction associated", (0,0.25)))
	p53_meg3_montecarlo_means = dropdims(mean(montecarlomeasurementsw[:,:,:,:,2], dims=1), dims=1)
	groupeddotplot!([1,1,1,1,1,1,2,2,2,2,2,2], group = [2,2,2,1,1,1,2,2,2,1,1,1], p53_meg3_montecarlo_means |> vec, mode = :none, label="", marker=(12, 0.75, :rect, repeat([:orange, :darkblue, :darkred]), stroke(0)))
	groupeddotplot!(rnagroups, p53_meg3_montecarlo, group = doxgroups, mode = :density, label="", marker=(8, 0.75, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
	savefig(joinpath(outputdir, "p53-meg3-montecarlo.png"))
end

# ╔═╡ 78883b10-5756-11eb-32cf-b95fdeebda79
md"""
### Normalized Monte Carlo fraction associated
This approach was abandoned because we determined it was over-normalizing the data, but still included here for completeness.

#### Check unusual cases
"""

# ╔═╡ df7dce20-5756-11eb-3796-037e8e3c49f9
let
	p1 = boxplot(normalizedmontecarlomeasurements[:, :, 1, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p2 = boxplot(normalizedmontecarlomeasurements[:, :, 2, 1, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p3 = boxplot(normalizedmontecarlomeasurements[:, :, 1, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p4 = boxplot(normalizedmontecarlomeasurements[:, :, 2, 2, 2], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))

	plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "Mdm2-p53")
end

# ╔═╡ faf5aa60-5756-11eb-3b5b-f1ad5e004af2
md"""
Potential outliers in box plots:
"""

# ╔═╡ ff845a40-5756-11eb-2ce8-d5ba8458e2ae
[normalizedmontecarlomeasurements[4,1,1,1,2],
 normalizedmontecarlomeasurements[9,2,1,1,2],
 normalizedmontecarlomeasurements[5,1,1,2,2],
 normalizedmontecarlomeasurements[7,3,1,2,2]]

# ╔═╡ d2282670-5757-11eb-0bd9-b1aaf0b37d70
md"""Evaluating z-scores within each replicate:"""

# ╔═╡ d99fec80-5757-11eb-39cd-4b0febbf61e3
normmontecarloflat = [normalizedmontecarlomeasurements[:, :, 1, 1, 2] normalizedmontecarlomeasurements[:, :, 2, 1, 2] normalizedmontecarlomeasurements[:, :, 1, 2, 2] normalizedmontecarlomeasurements[:, :, 2, 2, 2]]

# ╔═╡ edee22b2-5757-11eb-317a-d3e4d31486cd
normmontecarlozscores = [zscore(normmontecarloflat[:,i]) for i ∈ axes(normmontecarloflat, 2)]

# ╔═╡ a05efb90-5758-11eb-10ca-fb9fb8e05821
md"""
- In [1,1] Replicate 1: 1 value, Cell 4 (10%) outside 1.96. None outside 2.58.
- In [1,1] Replicate 2: 1 value, Cell 9 (10%) outside 1.96. None outside 2.58.
- In [1,1] Replicate 3: 1 value, Cell 2 (10%) outside 1.96. None outside 2.58.
- In [1,2] Replicate 1: 1 value, Cell 5 (10%) outside 1.96. None outside 2.58.
- In [1,2] Replicate 2: All within 1.96
- In [1,2] Replicate 3: All within 1.96
- In [2,1] Replicate 1: 1 value, Cell 6 (10%) outside 1.96 barely. None outside 2.58.
- In [2,1] Replicate 2: All within 1.96
- In [2,1] Replicate 3: 1 value, Cell 7 (10%) outside 1.96. None outside 2.58.
- In [2,2] Replicate 1: All within 1.96
- In [2,2] Replicate 2: All within 1.96
- In [2,2] Replicate 3: 1 value, Cell 8 (10%) outside 1.96. None outside 2.58.
"""

# ╔═╡ 7323e450-5759-11eb-1c62-7ba711bf418b
md"""
None of these potential outliers are due to data entry error or measurement error. Sampling a different population is possible but unknown. These potential outliers are relatively modest. In addition, since an ANOVA will be used, must be careful not to delete points from a balanced design. Therefore, no outliers will be removed.
"""

# ╔═╡ df4eee90-5759-11eb-02ef-53a0c6f31f07
md"""
#### Check assumption of homoscedasticity
"""

# ╔═╡ 8c2a8790-575b-11eb-34be-99f584826d43
md"""
Based on the ZResid/ZPred plot, the data appears heteroscedastic, although not much relationship between residual and predicted.

Conducted Levene's test, a 1-way ANOVA of the absolute deviation between each value and the group mean, by replicate.
"""

# ╔═╡ 9aade460-575b-11eb-29da-4ff6881ceb60
levene(normmontecarloflat)

# ╔═╡ b8f00430-575b-11eb-0ada-7f0001326feb
md"""
By Levene's test, the data is likely heteroscedastic.

Group sizes are equal, but samples are on the smaller side. CIs and significance tests may be affected.
"""

# ╔═╡ 1326cdce-575c-11eb-0f72-cf6eebe27bdf
md"""
#### Check assumption of normality
"""

# ╔═╡ 3c6921b0-575d-11eb-160f-6f06db5965f3
md"""
There are a couple notable deviations but for the *most* part they appear close to normal.
"""

# ╔═╡ 59314bb0-575d-11eb-2b54-e1bf609e6d36
normmontecarloskewness = [skewness(normmontecarloflat[:,i]) for i ∈ axes(normmontecarloflat, 2)]

# ╔═╡ 682c7090-575d-11eb-3b12-d50d450d3fce
normmontecarlokurtosis = [kurtosis(normmontecarloflat[:,i]) for i ∈ axes(normmontecarloflat, 2)]

# ╔═╡ 80f22522-575d-11eb-211b-59fe55838280
md"""
From skewness and kurtosis statistics, none was greater in value than 1.96 (5% level).
"""

# ╔═╡ 91c676d0-575d-11eb-2fd5-97ca14bfdec0
md"""
#### Conduct ANOVA tests
Regular ANOVA
"""

# ╔═╡ b7624d5e-575d-11eb-1221-974537c3b48e
normmontecarloresult = anova(
    permutedims(normalizedmontecarlomeasurements[:, :, :, :, 2], (1, 3, 2, 4)),
    [fixed, subject, fixed],
    factornames = ["Doxycycline", "Replicate", "RNA"],
)

# ╔═╡ 8f0dad80-575f-11eb-2312-adae92d50e83
md"""
### Check presence of interactions
"""

# ╔═╡ a5489880-575f-11eb-1afa-c3bfaf5dccc8
plot(normmontecarloresult)

# ╔═╡ 4b211122-575f-11eb-0aa9-59fc4e9a2df7
md"""
The lines cross, but the F-statistic is not significant. Interpreting simple effects valid. Continuing to look within each RNA.
"""

# ╔═╡ f325b3c0-57ba-11eb-1cae-17d58dd4ec29
normmontecarloresultMEG3 = anova(
    permutedims(normalizedmontecarlomeasurements[:, :, :, 1, 2], (1, 3, 2)),
    [fixed, subject],
    factornames = ["Doxycycline", "Replicate"],
)

# ╔═╡ 0f2bc640-57bb-11eb-2632-e1058fbc9108
normmontecarloresultGAPDH = anova(
    permutedims(normalizedmontecarlomeasurements[:, :, :, 2, 2], (1, 3, 2)),
    [fixed, subject],
    factornames = ["Doxycycline", "Replicate"],
)

# ╔═╡ 408e5d30-5764-11eb-124f-4f2630d1d908
md"""
#### Interpret main effects
Doxycycline had a small but nonsignificant effect with MEG3 (1.55% to 3.72%). Increasing MEG3 concentration may slightly increase the binding of MEG3 and p53, but evidence is weak.

Doxycyline had no significant effect with GAPDH (0.678% to 1.19%). Doxycycline shouldn't do anything to GAPDH, so this is expected.
"""

# ╔═╡ 795d7b9e-6bbe-11eb-1871-777e334518d5
md"""
## Localization/molecule maps
"""

# ╔═╡ abdc8f60-6bc5-11eb-1203-c36ca41662e8
md"""
#### Examples
"""

# ╔═╡ 892650c0-6bbe-11eb-3154-eb5c97e8c955
let
	# Exp 3 (p53-MEG3)
	# A2
	insetx, insety = [13800, 15200], [20300, 21700]
	r = experimentresults[2][1][1][2]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "3 - A2 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "3 - A2 dSTORM 1400nm.png"))
	p1 = localizationsplot(r, insetbox = [insetx, insety], include_scalebar = true)
	p2 = insetplot(r, insetx, insety, include_scalebar = true)
	plot(p1, p2, layout = grid(1,2), size=(1024,512), fmt = :png)
end

# ╔═╡ b87bb41e-6bcb-11eb-2a50-81d8ce260449
let
	# Exp 3 (p53-MEG3)
	# B2
	insetx, insety = [19800, 21200], [20900, 22300]
	r = experimentresults[2][1][2][2]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "3 - B2 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "3 - B2 dSTORM 1400nm.png"))
	p1 = localizationsplot(r, insetbox = [insetx, insety], include_scalebar = true)
	p2 = insetplot(r, insetx, insety, include_scalebar = true)
	plot(p1, p2, layout = grid(1,2), size=(1024,512), fmt = :png)
end

# ╔═╡ 516bf770-6bcd-11eb-37ab-f9ebc03c3cc4
let
	# Exp 3 (p53-MEG3)
	# C10
	insetx, insety = [8700, 10100], [20900, 22300]
	r = experimentresults[2][1][3][10]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "3 - C10 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "3 - C10 dSTORM 1400nm.png"))
	p1 = localizationsplot(r, insetbox = [insetx, insety], include_scalebar = true)
	p2 = insetplot(r, insetx, insety, include_scalebar = true)
	plot(p1, p2, layout = grid(1,2), size=(1024,512), fmt = :png)
end

# ╔═╡ 3c58ffc0-6bcf-11eb-13be-ddb38df5499d
let
	# Exp 3 (p53-MEG3)
	# D10
	insetx, insety = [16650, 18050], [29600, 31000]
	r = experimentresults[2][1][4][10]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "3 - D10 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "3 - D10 dSTORM 1400nm.png"))
	p1 = localizationsplot(r, insetbox = [insetx, insety], include_scalebar = true)
	p2 = insetplot(r, insetx, insety, include_scalebar = true)
	plot(p1, p2, layout = grid(1,2), size=(1024,512), fmt=:png)
end

# ╔═╡ c50f0010-6ce4-11eb-17df-39ae7df11a90
let
	# Plot all
	r = localizationsplot.(experimentresults[2][j][i][k] for k ∈ 1:10 for i ∈ 1:4 for j ∈ 1:3)
	p = plot(r..., size=(2048, 2048), layout=grid(10,12), fmt=:png)
	savefig(joinpath(outputdir, "p53-meg3-all-localizations.png"))
	p
end

# ╔═╡ 9c92b290-56ae-11eb-2595-85845fe95f0e
md"""
## Functions
"""

# ╔═╡ 2623d480-56af-11eb-2ab8-0d0afd8c45dc
function zresid_zpred(data)
	z = zscore(data)
	zpred = repeat(mean(z, dims = 1), 10)
	zresid = z .- zpred
	zresid, zpred
end

# ╔═╡ a655a120-56ae-11eb-30be-459d7c3067dc
function zresid_zpred_plot(data)
	zresid, zpred = zresid_zpred(data)
	scatter(
		zresid,
		zpred,
		xaxis = ("Standardized Residual (ZResid)"),
		yaxis = ("Standardized Predicted Value (ZPred)"),
		legend = :none,
	)
end

# ╔═╡ c8f83d60-5518-11eb-2643-a30df2dd9fb9
zresid_zpred_plot(mediansflatw)

# ╔═╡ 3a3fa580-56ae-11eb-1f19-f3c06ef2d9d7
zresid_zpred_plot(montecarloflatw)

# ╔═╡ 59522180-575a-11eb-39c4-f9eaa27f7901
zresid_zpred_plot(normmontecarloflat)

# ╔═╡ 7922d1d0-56b0-11eb-261c-5f87ca587394
function qqnormplot(data)
	zresid, _ = zresid_zpred(data)
	p = [qqnorm(zresid[:, i]) for i ∈ 1:12]
	plot(p..., layout = grid(4, 3), legend = :none, plot_title = "MDM2-p53 qqnorm")
end

# ╔═╡ 6f727ed0-5519-11eb-1462-01c83aa975b6
qqnormplot(mediansflat)

# ╔═╡ 5ef3f050-56b0-11eb-31f9-419a4161b67e
qqnormplot(montecarloflatw)

# ╔═╡ 1823a8d0-575c-11eb-0fa1-893a914f69f7
qqnormplot(normmontecarloflat)

# ╔═╡ Cell order:
# ╟─56dcb910-5507-11eb-0b0a-a1c2aef89207
# ╟─90d2acb0-5507-11eb-39c4-df93c258dd98
# ╠═b11fdebe-5507-11eb-1001-537a20e51d0d
# ╟─86502ff0-550d-11eb-1c13-553251c061ee
# ╠═5e3db410-550d-11eb-0331-49a447464b07
# ╟─294e6a02-550e-11eb-3450-d9d5989ec08c
# ╟─9a94f370-5511-11eb-1f60-7fe83ed11def
# ╟─572ba5f2-5513-11eb-1937-3f03cec8252e
# ╟─864b181e-5513-11eb-0304-8983d37a8170
# ╟─27e24820-5514-11eb-1415-8711a194892d
# ╟─f0393540-5514-11eb-2f76-3514cc3bc85c
# ╟─f02df530-5515-11eb-0d9e-8dbb63ed2023
# ╟─ffa9f6c0-576e-11eb-0f11-3139a34ff859
# ╟─25253910-5516-11eb-0a7a-1d73ed51d4a1
# ╟─d8c96c5e-5517-11eb-0861-056d5e425f04
# ╟─87eb8f70-5518-11eb-1639-f91b4bcb5c9b
# ╠═74a25f10-59ab-11eb-3a45-bf3b0554a07a
# ╟─a87aac80-5518-11eb-38d4-5fe90140b820
# ╟─c8f83d60-5518-11eb-2643-a30df2dd9fb9
# ╟─f6e23cd0-5518-11eb-226a-c5751eb65b34
# ╟─0ffdd990-5519-11eb-376b-33be0cc966d7
# ╟─28a8da30-5519-11eb-2fb0-efe693c6e44c
# ╟─5f8de452-5519-11eb-0a2c-051d3e701dd5
# ╟─6f727ed0-5519-11eb-1462-01c83aa975b6
# ╟─c4dbef60-55bd-11eb-1822-49c3473000e8
# ╟─cf59a26e-55bd-11eb-1c4b-316a603b9007
# ╟─135dbe6e-55be-11eb-3a89-dd3d02c3ad4a
# ╟─6778a880-55be-11eb-1016-c59335beef0a
# ╟─9bf9aaf0-55be-11eb-17a2-75b3eec550f8
# ╟─265006c0-55c1-11eb-2ba8-3d628d096704
# ╠═c101d900-57ab-11eb-04ab-77cd4e0ff4b3
# ╟─19934eb0-57ba-11eb-3f45-11c68057c1f7
# ╟─25bb5660-57ba-11eb-3999-4bbf47d29d8a
# ╟─342cfff0-57ba-11eb-07d6-2b90717f7c7b
# ╟─fd2d4a40-57ab-11eb-2688-53b009a37680
# ╠═5f1ce3b0-55c1-11eb-3d03-b53516cfc976
# ╠═0fddbb50-5772-11eb-01cc-e57dbec14a6b
# ╟─2ed8ff90-568e-11eb-10a3-0b97b2425d2f
# ╟─98c74810-6af9-11eb-3988-1f4c77155983
# ╟─68a010a0-6a3a-11eb-07df-9d60e51f3aab
# ╟─26985da0-6afa-11eb-2675-69a745a338c4
# ╠═252c7190-6afa-11eb-00fd-5f527f9fa189
# ╟─25d0ce1e-5691-11eb-2cc8-49980b6b2161
# ╟─8c0ebd10-56a9-11eb-3e8c-95b53bf12061
# ╟─faa812c0-56aa-11eb-3b6b-99fe5f200b67
# ╟─012a7c50-56ab-11eb-3d50-c9820925f9ac
# ╟─1637a040-56ac-11eb-2433-49c37383e005
# ╟─87ec3520-56ac-11eb-2a1a-2981155d4168
# ╟─9ce1f5f0-56ac-11eb-22a5-d9b5d6a18978
# ╟─ba7d2760-56ac-11eb-11ca-d543a628c68b
# ╠═03d7ceae-6d78-11eb-3f4b-13cf782576c0
# ╟─c18cc050-56ad-11eb-261c-4177aa31ca7c
# ╟─2bba9b02-56ae-11eb-2e30-439dc20273f8
# ╟─3a3fa580-56ae-11eb-1f19-f3c06ef2d9d7
# ╟─0f8a00a0-56af-11eb-2dc6-abd2a45298a4
# ╟─b1ae0750-56af-11eb-0714-35b78e197b5e
# ╟─f9676320-56af-11eb-1481-ef9502edded8
# ╟─4e36f730-56b0-11eb-21e8-0d2e7d0a2562
# ╟─5ef3f050-56b0-11eb-31f9-419a4161b67e
# ╟─4423b3d0-56b2-11eb-22be-fd86deabef7f
# ╟─db6922c0-56b2-11eb-2601-5525d3f6bf70
# ╟─047df370-56b3-11eb-3036-2393031bd967
# ╟─3316f790-56b3-11eb-3d9c-e9ebf797c923
# ╟─6f1da2c0-56b3-11eb-2506-5f9a944735b6
# ╠═c29f8040-57a7-11eb-2cb4-7974bcf74e27
# ╟─2553a3a0-57a9-11eb-247d-a161d85d4410
# ╠═755d8510-56b3-11eb-3a0c-53b643ba2053
# ╠═28212d70-57a7-11eb-1c94-3fa38021ba8e
# ╟─248d1560-5753-11eb-2156-89cfb6742205
# ╟─174f2590-6a46-11eb-0621-73610de2d1c0
# ╟─0185bed0-6afb-11eb-279d-b9e008c964c3
# ╠═fecb59c0-6afa-11eb-15fe-cbab4e3f44e1
# ╟─78883b10-5756-11eb-32cf-b95fdeebda79
# ╟─df7dce20-5756-11eb-3796-037e8e3c49f9
# ╟─faf5aa60-5756-11eb-3b5b-f1ad5e004af2
# ╠═ff845a40-5756-11eb-2ce8-d5ba8458e2ae
# ╟─d2282670-5757-11eb-0bd9-b1aaf0b37d70
# ╟─d99fec80-5757-11eb-39cd-4b0febbf61e3
# ╟─edee22b2-5757-11eb-317a-d3e4d31486cd
# ╟─a05efb90-5758-11eb-10ca-fb9fb8e05821
# ╟─7323e450-5759-11eb-1c62-7ba711bf418b
# ╟─df4eee90-5759-11eb-02ef-53a0c6f31f07
# ╟─59522180-575a-11eb-39c4-f9eaa27f7901
# ╟─8c2a8790-575b-11eb-34be-99f584826d43
# ╟─9aade460-575b-11eb-29da-4ff6881ceb60
# ╟─b8f00430-575b-11eb-0ada-7f0001326feb
# ╟─1326cdce-575c-11eb-0f72-cf6eebe27bdf
# ╟─1823a8d0-575c-11eb-0fa1-893a914f69f7
# ╟─3c6921b0-575d-11eb-160f-6f06db5965f3
# ╟─59314bb0-575d-11eb-2b54-e1bf609e6d36
# ╟─682c7090-575d-11eb-3b12-d50d450d3fce
# ╟─80f22522-575d-11eb-211b-59fe55838280
# ╟─91c676d0-575d-11eb-2fd5-97ca14bfdec0
# ╠═b7624d5e-575d-11eb-1221-974537c3b48e
# ╟─8f0dad80-575f-11eb-2312-adae92d50e83
# ╟─a5489880-575f-11eb-1afa-c3bfaf5dccc8
# ╟─4b211122-575f-11eb-0aa9-59fc4e9a2df7
# ╠═f325b3c0-57ba-11eb-1cae-17d58dd4ec29
# ╠═0f2bc640-57bb-11eb-2632-e1058fbc9108
# ╟─408e5d30-5764-11eb-124f-4f2630d1d908
# ╟─795d7b9e-6bbe-11eb-1871-777e334518d5
# ╟─abdc8f60-6bc5-11eb-1203-c36ca41662e8
# ╠═892650c0-6bbe-11eb-3154-eb5c97e8c955
# ╠═b87bb41e-6bcb-11eb-2a50-81d8ce260449
# ╠═516bf770-6bcd-11eb-37ab-f9ebc03c3cc4
# ╠═3c58ffc0-6bcf-11eb-13be-ddb38df5499d
# ╠═c50f0010-6ce4-11eb-17df-39ae7df11a90
# ╟─9c92b290-56ae-11eb-2595-85845fe95f0e
# ╟─2623d480-56af-11eb-2ab8-0d0afd8c45dc
# ╟─a655a120-56ae-11eb-30be-459d7c3067dc
# ╟─7922d1d0-56b0-11eb-261c-5f87ca587394


text:
"""
    dualBinaryTree(ptr, N)
Draw the dual domain's binary partition tree.
"""
function dualBinaryTree(ptr, N)
    n = 0 # number of nodes in the tree
    jmax = length(ptr)+1
    for j = 1:jmax-1
        n += length(ptr[j])
    end
    n += N

    node_ind = (Array{Int64,1})[]
    ind = 1
    for j = 1:jmax-1
        lvl_nodes = (Int64)[]
        for i = 1:length(ptr[j])
            push!(lvl_nodes, ind)
            ind += 1
        end
        push!(node_ind, lvl_nodes)
    end
    push!(node_ind, [ind-1+i for i=1:N])

    BT = Graph(n)
    X = zeros(n,2)
    hv = ["(0,0)" for i = 1:n]
    for j = 1:jmax-1
        for i in 1:length(ptr[j])
            c1 = ptr[j][i][1]
            add_edge!(BT, Edge(node_ind[j][i], node_ind[j+1][c1]))
            X[node_ind[j+1][c1],:] = X[node_ind[j][i],:] + [0,1]
            hv[node_ind[j+1][c1]] = "($(j),$(c1-1))"
            if length(ptr[j][i]) == 2
                c2 = ptr[j][i][2]
                add_edge!(BT, Edge(node_ind[j][i], node_ind[j+1][c2]))
                X[node_ind[j+1][c1],:] = X[node_ind[j][i],:] + [-2^(jmax-j),1]
                X[node_ind[j+1][c2],:] = X[node_ind[j][i],:] + [2^(jmax-j),1]
                hv[node_ind[j+1][c2]] = "($(j),$(c2-1))"
            end
        end
    end
    return X, BT, node_ind, hv
end

function find_deepest_leaf(node, ptr)
    jmax = length(ptr)
    j, k = node
    while j <= jmax
        k = ptr[j][k][1]
        j += 1
    end
    return [j, k]
end


text:
function WilcoxonSignedRankTest()
    return WilcoxonSignedRankTest(())
end

function WilcoxonSignedRankTest(arg0::NaNStrategy, arg1::TiesStrategy)
    return WilcoxonSignedRankTest((NaNStrategy, TiesStrategy), arg0, arg1)
end

function wilcoxon_signed_rank(obj::WilcoxonSignedRankTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "wilcoxonSignedRank", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function wilcoxon_signed_rank_test(obj::WilcoxonSignedRankTest, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jboolean)
    return jcall(obj, "wilcoxonSignedRankTest", jdouble, (Vector{jdouble}, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end


text:
using AbstractAlgebra
using Random

@testset "Input-output AbstractAlgebra" begin

    rng = Random.MersenneTwister(42)

    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:lex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :lex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :degrevlex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y, x*y - y^2, x*y - x^2]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :degrevlex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y, x*y - y^2, x*y - x^2]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :degrevlex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:deglex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y, x*y - y^2, x*y - x^2]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :deglex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:deglex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y, x*y - y^2, x*y - x^2]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :deglex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3//4, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :degrevlex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:lex)
    fs = [x^2*y + 3//4, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :lex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y, x*y - y^2, x*y - x^2]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :degrevlex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y, x*y - y^2, x*y - x^2]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :degrevlex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:deglex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y, x*y - y^2, x*y - x^2]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :deglex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:deglex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y, x*y - y^2, x*y - x^2]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    meta = Groebner.set_metaparameters(ring, :deglex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs
end

@testset "Input-output generic :hasparent" begin

    rng = Random.MersenneTwister(42)

    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:lex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    ring.origring = :hasparent
    meta = Groebner.set_metaparameters(ring, :lex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    ring.origring = :hasparent
    meta = Groebner.set_metaparameters(ring, :degrevlex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:degrevlex)
    fs = [x^2*y + 3//4, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    ring.origring = :hasparent
    meta = Groebner.set_metaparameters(ring, :degrevlex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs

    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:lex)
    fs = [x^2*y + 3//4, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    ring.origring = :hasparent
    meta = Groebner.set_metaparameters(ring, :lex, false, false, :exact, rng)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs, meta)
    @test fsfs == fs
end


text:
include("abstractfeatures.jl")
using DataFrames
using ..ChemistryFeaturization: elements
using ..ChemistryFeaturization.AbstractType: AbstractCodec, AbstractAtoms
using ..ChemistryFeaturization.Codec: SimpleCodec
using ..ChemistryFeaturization.Utils.OrbitalFeatureUtils:
    _name_to_econf, _indexorbital, _econf_to_name
using ..ChemistryFeaturization.Data: valenceshell_conf_df
using SparseArrays

"""
    OrbitalFeatureDescriptor

Describes the orbital configuration of an element.
"""
struct OrbitalFeatureDescriptor <: AbstractEnvironmentFeatureDescriptor
    encoder_decoder::AbstractCodec
    function OrbitalFeatureDescriptor()
        # trim down the DataFrame to have only the Symbol, and valence shell configuration (extracted from "Electronic Structure")
        new(SimpleCodec(default_ofd_encode, default_ofd_decode))
    end
end

function (ofd::OrbitalFeatureDescriptor)(a::AbstractAtoms)
    @assert all([el in valenceshell_conf_df[:, :Symbol] for el in elements(a)]) "All elements must be valid and accounted for in the periodic table!"
    ofd.encoder_decoder(ofd, a)
end

# encode
(sc::SimpleCodec)(ofd::OrbitalFeatureDescriptor, a::AbstractAtoms) = sc.encode_f(ofd, a)

# decode
(sc::SimpleCodec)(
    ofd::OrbitalFeatureDescriptor,
    encoded_features::SparseMatrixCSC{Tv,Ti},
) where {Tv,Ti} = sc.decode_f(ofd, encoded_features)

function output_shape(ofd::OrbitalFeatureDescriptor, sc::SimpleCodec)
    # return the number
end

Base.show(io::IO, efd::OrbitalFeatureDescriptor) =
    print(io, "OrbitalFeatureDescriptor created.")

function Base.show(io::IO, ::MIME"text/plain", efd::OrbitalFeatureDescriptor)
    st = "OrbitalFeatureDescriptor created."
    print(io, st)
end

"""
    default_ofd_encode(ofd::OrbitalFeatureDescriptor, a::AbstractAtoms)

Default encoding scheme for OrbitalFeatureDescriptor, which returns a
sparse matrix, where column_i is a sparse representation of the electronic
configuration of the ith element in the [Atoms](@ref atoms) object.
"""
function default_ofd_encode(ofd::OrbitalFeatureDescriptor, a::AbstractAtoms)
    I = Vector{Int16}()
    J = Vector{Int16}()
    V = Vector{Int16}()
    col = 0
    for i in elements(a)
        # X, Y are equivalent to I, V for a SparseVector
        X, Y = _name_to_econf(i)
        col += 1    # column number

        for i = 1:length(X)
            append!(I, X[i])
            append!(J, col)
            append!(V, Y[i])
        end
    end
    sparse(I, J, V) # create a SparseMatrix with the valence shell configuration of each element as a column vector
end

"""
    default_ofd_decode(
        ofd::OrbitalFeatureDescriptor,
        encoded_features::SparseArrays.AbstractSparseMatrixCSC{Tv,Ti}
    ) where {Tv,Ti}

Default decoding scheme for OrbitalFeatureDescriptor, which returns the
vector of elements encoded, given the encoded sparse matrix.
"""
function default_ofd_decode(
    ofd::OrbitalFeatureDescriptor,
    encoded_features::SparseArrays.AbstractSparseMatrixCSC{Tv,Ti},
) where {Tv,Ti}
    elements = String[]
    for i = 1:size(encoded_features)[2]
        push!(elements, _econf_to_name(encoded_features[:, i]))
    end
    return elements
end


text:
#####
##### `mean`
#####

_denom(x, dims::Colon) = length(x)
_denom(x, dims::Integer) = size(x, dims)
_denom(x, dims) = mapreduce(i->size(x, i), Base.mul_prod, unique(dims), init=1)

# TODO: We have `mean(f, x; dims)` as of 1.3.0-DEV.36
# https://github.com/JuliaDiff/ChainRules.jl/issues/85
function rrule(::typeof(mean), x::AbstractArray{<:Real}; dims=:)
    y_sum, sum_pullback = rrule(sum, x; dims=dims)
    n = _denom(x, dims)
    function mean_pullback(ȳ)
        _, ∂sum_x = sum_pullback(ȳ)
        ∂x = unthunk(∂sum_x) / n
        return (NoTangent(), ∂x)
    end
    return y_sum / n, mean_pullback
end

#####
##### variance
#####

function rrule(
    ::typeof(Statistics.var),
    x::AbstractArray{<:Number};
    corrected::Bool=true,
    dims=:,
    mean=mean(x, dims=dims)
)
    y = Statistics.var(x; corrected=corrected, mean=mean, dims=dims)
    function variance_pullback(dy)
        pre = 2 // (_denom(x, dims) - corrected)
        dx = pre .* unthunk(dy) .* (x .- mean)
        return (NoTangent(), ProjectTo(x)(dx))
    end
    y, variance_pullback
end

function rrule(
    ::typeof(Statistics.std),
    x::AbstractArray{<:Number};
    corrected::Bool=true,
    dims=:,
    mean=mean(x, dims=dims)
)
    y = Statistics.std(x; corrected=corrected, mean=mean, dims=dims)
    function std_pullback(dy)
        pre = 1 // (_denom(x, dims) - corrected)
        dx = pre .* unthunk(dy) .* (x .- mean) ./ y
        return (NoTangent(), ProjectTo(x)(dx))
    end
    y, std_pullback
end


text:
#The m-by-n Hilbert matrix has matrix elements
# H_{ij} = 1/(i+j-1)
export Hilbert, InverseHilbert
"""
[`Hilbert` matrix](http://en.wikipedia.org/wiki/Hilbert_matrix)

```julia
julia> A=Hilbert(5)
Hilbert{Rational{Int64}}(5,5)

julia> Matrix(A)
5x5 Array{Rational{Int64},2}:
 1//1  1//2  1//3  1//4  1//5
 1//2  1//3  1//4  1//5  1//6
 1//3  1//4  1//5  1//6  1//7
 1//4  1//5  1//6  1//7  1//8
 1//5  1//6  1//7  1//8  1//9

julia> Matrix(Hilbert(5))
5x5 Array{Rational{Int64},2}:
 1//1  1//2  1//3  1//4  1//5
 1//2  1//3  1//4  1//5  1//6
 1//3  1//4  1//5  1//6  1//7
 1//4  1//5  1//6  1//7  1//8
 1//5  1//6  1//7  1//8  1//9
```
Inverses are also integer matrices:

```julia
julia> inv(A)
5x5 Array{Rational{Int64},2}:
    25//1    -300//1     1050//1    -1400//1     630//1
  -300//1    4800//1   -18900//1    26880//1  -12600//1
  1050//1  -18900//1    79380//1  -117600//1   56700//1
 -1400//1   26880//1  -117600//1   179200//1  -88200//1
   630//1  -12600//1    56700//1   -88200//1   44100//1
```
"""
struct Hilbert{T} <: AbstractMatrix{T}
    m :: Int
    n :: Int
end
Hilbert(T::Type, m::Integer, n::Integer) = Hilbert{T}(m, n)
Hilbert(T::Type{<:Integer}, m::Integer, n::Integer) = Hilbert(Rational{T}, m, n)
Hilbert(m::Integer, n::Integer) = Hilbert(Int, m, n)
Hilbert(n::Integer) = Hilbert(n, n)

# Define its size
size(H::Hilbert, dim::Integer) = dim==1 ? H.m : dim==2 ? H.n : 1
size(H::Hilbert)= size(H,1), size(H,2)

# Index into a Hilbert matrix
function getindex(H::Hilbert{T}, i::Integer, j::Integer) where {T}
    return one(T)/(i+j-1)
end

# Dense version
Matrix(H::Hilbert) = [H[i,j] for i=1:size(H,1), j=1:size(H,2)]

# Some properties
ishermitian(H::Hilbert) = H.m==H.n
isposdef(H::Hilbert) = H.m==H.n
det(H::Hilbert) = inv(det(inv(H)))

# Inverse of a Hilbert matrix
struct InverseHilbert{T} <: AbstractMatrix{T}
    n :: Int
end
InverseHilbert(T::Type, n::Integer) = InverseHilbert{T}(n)
InverseHilbert(n::Integer) = InverseHilbert(Int, n)

# Define its size
size(A::InverseHilbert, dim::Integer) = dim==1 || dim==2 ? A.n : 1
size(A::InverseHilbert) = size(A,1), size(A,2)

# Index into a inverse Hilbert matrix
function getindex(A::InverseHilbert{T}, i::Integer, j::Integer) where {T}
    return T((-1)^(i+j)*(i+j-1)*binomial(A.n+i-1,A.n-j)*
              binomial(A.n+j-1,A.n-i)*binomial(i+j-2,i-1)^2)
end
# Use bigger numbers in case of bigger types
function getindex(A::InverseHilbert{T}, i::Integer, j::Integer) where {T<:Union{BigInt,BigFloat}}
    N = big(A.n)
    return T((-1)^(i+j)*(i+j-1)*binomial(N+i-1,N-j)*
              binomial(N+j-1,N-i)*binomial(big(i+j-2),i-1)^2)
end

# Explicit formula for the determinant
det(A::InverseHilbert{T}) where {T} = prod(T,(2k+1)*binomial(2k,k)^2 for k=1:A.n-1)

# Dense version
Matrix(A::InverseHilbert) = [A[i,j] for i=1:size(A,1), j=1:size(A,2)]

# Define Base.inv
function inv(H::Hilbert{T}) where {T}
    H.m == H.n || throw(ArgumentError("Works only for square Hilbert matrices."))
    return InverseHilbert(T,H.n)
end
function inv(A::InverseHilbert{T}) where {T}
    HT = promote_type(T,typeof(Rational(one(T))))
    return Hilbert(HT,A.n)
end