text:
# Pillai et al. style Meta-foodchain modified from original code by Sebastian Kusch

using Fastnet

na = ["Producer","Consumer","Predator","Empty"]

present1 = (1,2,3)
present2 = (2,3)
present3 = 3

L10=LinkType(1,4)
L20=LinkType(2,4)
L30=LinkType(3,4)
L21=LinkType(2,1)
L31=LinkType(3,1)
L32=LinkType(3,2)

linklist=[L10,L20,L30,L21,L31,L32]
la=["1-0","2-0","3-0","2-1","3-1","3-2"]

routes1=(1,2,3)
routes2=(4,5)
routes3=6

# parameters
dispersal=[1.0, 0.5, 0.25]
mortality=[1.0, 0.1, 0.01]

#network
z=5
n=10000
k=(n*z)÷2
net=FastNet(n,k,4,linklist; nodealias=na,linkalias=la)
randomgraph!(net)
rand_init = rand(1:4,n)                          # generate random patch state
for i=1:n                                        # Assign patch state at random
    nodestate_f!(net,node_f(net,i),rand_init[i])
end

function rates!(rates,t)    # This functins computes the rates of processes
    rates[1]=mortality[1]*countnodes(net,present1)             # species 1 death
    rates[2]=mortality[2]*countnodes(net,present2)                # species 2 death
    rates[3]=mortality[3]*countnodes(net,present3)                     # species 3 death
    rates[4]=dispersal[1]*countlinks(net,routes1)             # species 1 migration
    rates[5]=dispersal[2]*countlinks(net,routes2)                  # species 2 migration
    rates[6]=dispersal[3]*countlinks(net,routes3)                # species 3 migration
    nothing
end

function loss1()        # This is what we do when a primary producer dies
    nde=randomnode(net,present1)          # Find a random node of state 1,2 or 3
    nodestate_f!(net,nde,4)                  # Set the state of the node to empty patch
end

function loss2()        # This is what we do when a species of state 2 dies
    nde=randomnode(net,present2)             # Find a random node of state 2 or 3
    nodestate_f!(net,nde,1)                  # Set the state of the node to S1
end

function loss3()                     # This is what we do when a species of state 3 dies
    nde=randomnode(net,present3)                  # Find a random node of state 3
    nodestate!(net,nde,2)                  # Set the state of the node to S2
end

function disperse1()       # This is what we do when a species 1 migration process is triggered
    alink=randomlink(net,routes1)    # Find a random species 1 migration link
    src = linksrc(net,alink)
    dst = linkdst(net,alink)
    if nodestate(net,src) == 4    # determine migration destination in undirected link and update to S1
        nodestate!(net,src,1)
    else
        nodestate!(net,dst,1)
    end
end

function disperse2()       # This is what we do when a species 2 migration process is triggered
    alink=randomlink(net,routes2)         # Find a random species 2 migration link
    src = linksrc(net,alink)
    dst = linkdst(net,alink)
    if nodestate(net,src) == 1  # determine migration destination in undirected link and update to S2
        nodestate!(net,src,2)
    else
        nodestate!(net,dst,2)
    end
end

function disperse3()       # This is what we do when a species 3 migration process is triggered
    alink=randomlink(net,routes3)                # Find a random species 3 migration link
    src = linksrc(net,alink)
    dst = linkdst(net,alink)
    if nodestate(net,src) == 2 # determine migration destination in undirected link and update to S3
        nodestate!(net,src,3)
    else
        nodestate!(net,dst,3)
    end
end

# Simulation
sim=FastSim(net,rates!,[loss1,loss2,loss3,disperse1,disperse2,disperse3])   # initialize the simulation
@time runsim!(sim,10,2)


text:
using Test
import BenchmarkTools

@testset "solve: zero allocations" begin
    fs = (sin, cos, x -> -sin(x))
    x0 = (3, 4)
    Ms = (
        Order1(),
        Order2(),
        Order5(),
        Order8(),
        Order16(),
        Roots.Order1B(),
        Roots.Order2B(),
        Roots.Bisection(),
        Roots.A42(),
        Roots.AlefeldPotraShi(),
        Roots.Brent(),
        Roots.Newton(),
        Roots.Halley(),
        Roots.Schroder(),
    ) # not Order0(), FalsePosition()

    for M in Ms
        @test BenchmarkTools.@ballocated(solve(ZeroProblem($fs, $x0), $M)) == 0
    end
end

@testset "simple: zero allocations" begin
    @test BenchmarkTools.@ballocated(Roots.bisection(sin, 3, 4)) == 0
    @test BenchmarkTools.@ballocated(Roots.secant_method(sin, 3)) == 0
    @test BenchmarkTools.@ballocated(Roots.muller(sin, 2.9, 3.0, 3.1)) == 0
    @test BenchmarkTools.@ballocated(Roots.newton((sin, cos), 3)) == 0
    @test BenchmarkTools.@ballocated(Roots.dfree(sin, 3)) == 0
end

ext: jl


text:
using ADVUMPS
using ADVUMPS:magofβ,magofdβ,num_grad
using Plots
using Zygote
using JLD2
using FileIO

# magnetisation
begin
    magplot = plot()
    for D = [2,4,8,16,32]
        mag = []
        for β = 0.41:0.0025:0.48
            printstyled("D = $(D), β = $(β) \n"; bold=true, color=:red)

            env = vumps_env(Ising(),β,D)
            mag = [mag;magnetisation(env, Ising(), β)]

            chkp_file = "./data/Ising()_β$(round(β+0.0025,digits=4))_D$(D).jld2"
            if isfile(chkp_file) == false
                save(chkp_file, "env", env)
            end
        end
        β = 0.41:0.0025:0.48
        magplot = plot!(β,mag,seriestype = :scatter,title = "magnetisation", label = "VUMPS D = $(D)", lw = 3)
    end

    tmag = []
    for β = 0.41:0.001:0.48
        tmag = [tmag;magofβ(Ising(), β)]
    end
    β = 0.41:0.001:0.48
    magplot = plot!(β,tmag,label = "exact", lw = 2)
    xlabel!("β")
    ylabel!("mag")
    savefig(magplot,"./plot/2Disingmag.svg")
end

# dmag/dβ
begin
    dmagplot = plot()
    for D = [2,4,8]
        dmag = []
        for β = 0.1:0.05:1.0
            printstyled("D = $(D), β = $(β) \n"; bold=true, color=:red)
            function foo(β)
                env = vumps_env(Ising(),β,D)
                Zygote.@ignore begin
                    chkp_file = "./data/Ising()_β$(round(β+0.05,digits=2))_D$(D).jld2"
                    if isfile(chkp_file) == false
                        save(chkp_file, "env", env)
                    end
                end
                magnetisation(env,Ising(),β)
            end
            dmag = [dmag;Zygote.gradient(foo,β)[1]]
        end
        β = 0.1:0.05:1.0
        dmagplot = plot!(β,dmag,seriestype = :scatter,title = "dmag/dβ", label = "ADVUMPS D = $(D)", lw = 2)
    end
    tdmag = []
    for β = 0.1:0.01:1.0
        tdmag = [tdmag;magofdβ(Ising(), β)]
    end
    β = 0.1:0.01:1.0
    dmagplot = plot!(β,tdmag, label = "exact", lw = 2)
    xlabel!("β")
    ylabel!("dmag/dβ")
    savefig(dmagplot,"./plot/2Disingdmag.svg")
end

# energy
begin
    engplot = plot()
    for D = [2,4,8]
        eng = []
        for β = 0.1:0.05:1.0
            printstyled("D = $(D), β = $(β) \n"; bold=true, color=:red)
            function foo(β)
                env = vumps_env(Ising(),β,D)
                Zygote.@ignore begin
                    chkp_file = "./data/Ising()_β$(round(β+0.05,digits=2))_D$(D).jld2"
                    if isfile(chkp_file) == false
                        save(chkp_file, "env", env)
                    end
                end
                -log(Z(env))
            end
            eng = [eng;Zygote.gradient(foo,β)[1]]
        end
        β = 0.1:0.05:1.0
        engplot = plot!(β,eng,seriestype = :scatter,title = "-dlog(Z)/dβ", label = "ADVUMPS D = $(D)", lw = 2)
    end
    β = 0.1:0.01:1.0
    tengmag = []
    for β = 0.1:0.01:1.0
        env = vumps_env(Ising(),β,8)
        chkp_file = "./data/Ising()_β$(round(β+0.01,digits=2))_D$(8).jld2"
        if isfile(chkp_file) == false
            save(chkp_file, "env", env)
        end
        tengmag = [tengmag;energy(env,Ising(),β)]
    end
    engplot = plot!(β,tengmag, label = "VUMPS D = 8", lw = 2)
    xlabel!("β")
    ylabel!("energy")
    savefig(engplot,"./plot/2Disingene.svg")
end


text:
function show(io::IO, p::MPoly)
    print(io, typeof(p))
    print(io, '(')
    print(io, p)
    print(io, ')')
end

function printmonomial(io::IO, m, vars)
    first = true
    for i = 1:length(m)
        if m[i] >= 1 || m[i] <= -1
            if first
                first = false
            else
                print(io, '*')
            end
            print(io, vars[i])
            if m[i] >= 2 || m[i] <= -1
                print(io, '^')
                print(io, m[i] > 0 ? m[i] : "($(m[i]))")
            end
        end
    end
end

function print{T}(io::IO, p::MPoly{T})
    first = true
    for (m, c) in p
        if first
            if T <: Complex
                if m == zeros(Int, length(m))
                    print(io, "$c")
                else
                    print(io, "($c)*")
                end
            else
                if m == zeros(Int, length(m)) || c != one(T)
                    print(io, c)
                end
            end
            printmonomial(io, m, vars(p))
            first = false
        else
            if T <: Complex
                if m == zeros(Int, length(m))
                    print(io, " + $c")
                else
                    print(io, " + ($c)*")
                end
            else
                if c >= zero(T)
                    if m != zeros(Int, length(m)) && c == one(T)
                        print(io, " + ")
                    else
                        print(io, " + $c")
                    end
                else
                    print(io, " - $(-c)")
                end
            end
            printmonomial(io, m, vars(p))
        end
    end
    if first
        print(io, zero(T))
    end
end

showcompact(io::IO, p::MPoly) = print(io, p)


text:
#
# fft.jl -
#
# Implementation of FFT and circulant convolution operators.
#
#------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (C) 2017-2021, Éric Thiébaut.
# Copyright (C) 2015-2016, Éric Thiébaut, Jonathan Léger & Matthew Ozon.
#

module FFTs

# Be nice with the caller: re-export `fftshift` and `ifftshift` but not `fft`,
# `ifft`, etc. as the `FFTOperator` is meant to replace them.
export
    CirculantConvolution,
    FFTOperator,
    fftfreq,
    fftshift,
    goodfftdim,
    goodfftdims,
    ifftshift,
    rfftdims

using ..Foundations
using ..LazyAlgebra
using ..LazyAlgebra:
    @certify, bad_argument, bad_size, compose
import ..LazyAlgebra:
    adjoint, apply!, vcreate, MorphismType, mul!,
    input_size, input_ndims, input_eltype,
    output_size, output_ndims, output_eltype,
    identical

import Base: *, /, \, inv, show

using ArrayTools

import AbstractFFTs: Plan, fftshift, ifftshift

using FFTW
import FFTW: fftwNumber, fftwReal, fftwComplex, FFTWPlan, cFFTWPlan, rFFTWPlan

# All planning flags.
const PLANNING = (FFTW.ESTIMATE | FFTW.MEASURE | FFTW.PATIENT |
                  FFTW.EXHAUSTIVE | FFTW.WISDOM_ONLY)

# The time needed to allocate temporary arrays is negligible compared to the
# time taken to compute a FFT (e.g., 5µs to allocate a 256×256 array of double
# precision complexes versus 1.5ms to compute its FFT).  We therefore do not
# store any temporary arrays in the FFT operator.  Only the FFT plans are
# cached in the operator.

#------------------------------------------------------------------------------
# Extend LazyAlgebra framework for FFTW plans.
#
# This simplify a lot the implementation of FFT and circulant convolution
# operators without loss of performances.

macro checksize(name, arg, dims)
    return quote
        size($(esc(arg))) == $(esc(dims)) ||
            bad_size($(esc(name)), " must have dimensions ", $(esc(dims)))
    end
end

input_size(P::FFTWPlan) = P.sz
output_size(P::FFTWPlan) = P.osz
#input_strides(P::FFTWPlan) = P.istride
#output_strides(P::FFTWPlan) = P.ostride
flags(P::FFTWPlan) = P.flags

destroys_input(A::FFTWPlan) =
    (flags(A) & (FFTW.PRESERVE_INPUT|FFTW.DESTROY_INPUT)) == FFTW.DESTROY_INPUT

preserves_input(A::FFTWPlan) =
    (flags(A) & (FFTW.PRESERVE_INPUT|FFTW.DESTROY_INPUT)) == FFTW.PRESERVE_INPUT

# Extend `vcreate` for FFTW plans.  Rationale: result must be of predictible type
# and checking input argument is skipped (this will be done by `apply!`).
#
# Create result for an in-place complex-complex forward/backward FFT
# transform.
function vcreate(::Type{Direct},
                 A::cFFTWPlan{Complex{T},K,true,N},
                 x::StridedArray{Complex{T},N},
                 scratch::Bool) where {T<:fftwReal,K,N}
    return (scratch && isa(x, Array) ? x : Array{Complex{T}}(undef, output_size(A)))
end

# Create result for an out-of-place complex-complex forward/backward FFT
# transform.
function vcreate(::Type{Direct},
                 A::cFFTWPlan{Complex{T},K,false,N},
                 x::StridedArray{Complex{T},N},
                 scratch::Bool) where {T<:fftwReal,K,N}
    return Array{Complex{T}}(undef, output_size(A))
end

# Create result for a real-complex or a complex-real forward/backward FFT
# transform.  The result is necessarily a new array whatever the `scratch`
# flag.
function vcreate(::Type{Direct},
                 A::rFFTWPlan{T,K,false,N},
                 x::StridedArray{T,N},
                 scratch::Bool) where {T<:fftwReal,K,N}
    return Array{Complex{T}}(undef, output_size(A))
end

function vcreate(::Type{Direct},
                 A::rFFTWPlan{Complex{T},K,false,N},
                 x::StridedArray{Complex{T},N},
                 scratch::Bool) where {T<:fftwReal,K,N}
    return Array{T}(undef, output_size(A))
end

# Extend `apply!` for FFTW plans.  We want to compute:
#
#    y = α⋅F⋅x + β⋅y
#
# with as few temporaries as possible.  If β = 0, then there are no needs
# to save the contents of y which can be used directly for the output of
# the transform.  Extra checks are required to make sure the contents x is
# not damaged unless scratch is true.  It tuns out that the implementation
# depends on the type of transform so several versions are coded below.

# Apply in-place complex-complex forward/backward FFT transform.
function apply!(α::Number,
                ::Type{Direct},
                A::cFFTWPlan{Complex{T},K,true,N},
                x::StridedArray{Complex{T},N},
                scratch::Bool,
                β::Number,
                y::StridedArray{Complex{T},N}) where {T<:fftwReal,N,K}
    @checksize "argument" x  input_size(A)
    @checksize "result"   y output_size(A)
    if α == 0
        vscale!(y, β)
    elseif β == 0
        mul!(y, A, vscale!(y, α, x))
    elseif scratch
        vcombine!(y, α, mul!(x, A, x), β, y)
    else
        z = copy(x)
        vcombine!(y, α, mul!(z, A, z), β, y)
    end
    return y
end

# Apply out-of-place complex-complex forward/backward FFT transform.
function apply!(α::Number,
                ::Type{Direct},
                A::cFFTWPlan{Complex{T},K,false,N},
                x::StridedArray{Complex{T},N},
                scratch::Bool,
                β::Number,
                y::StridedArray{Complex{T},N}) where {T<:fftwReal,N,K}
    @checksize "argument" x  input_size(A)
    @checksize "result"   y output_size(A)
    if α == 0
        vscale!(y, β)
    elseif β == 0
        safe_mul!(y, A, x, scratch && x !== y)
        α == 1 || vscale!(y, α)
    else
        vcombine!(y, α, safe_mul(A, x, scratch), β, y)
    end
    return y
end

# Apply real-to-complex forward transform.  The transform is necessarily
# out-of-place.
function apply!(α::Number,
                ::Type{Direct},
                A::rFFTWPlan{T,K,false,N},
                x::StridedArray{T,N},
                scratch::Bool,
                β::Number,
                y::StridedArray{Complex{T},N}) where {T<:fftwReal,K,N}
    @checksize "argument" x  input_size(A)
    @checksize "result"   y output_size(A)
    if α == 0
        vscale!(y, β)
    elseif β == 0
        safe_mul!(y, A, x, scratch)
        α == 1 || vscale!(y, α)
    else
        vcombine!(y, α, safe_mul(A, x, scratch), β, y)
    end
    return y
end

# Apply complex-to-real (c2r) backward transform. Preserving input is not
# possible for multi-dimensional c2r transforms so we must copy the input
# argument x.
function apply!(α::Number,
                ::Type{Direct},
                A::rFFTWPlan{Complex{T},K,false,N},
                x::StridedArray{Complex{T},N},
                scratch::Bool,
                β::Number,
                y::StridedArray{T,N}) where {T<:fftwReal,K,N}
    @checksize "argument" x  input_size(A)
    @checksize "result"   y output_size(A)
    if α == 0
        vscale!(y, β)
    elseif β == 0
        safe_mul!(y, A, x, scratch)
        α == 1 || vscale!(y, α)
    else
        vcombine!(y, α, safe_mul(A, x, scratch), β, y)
    end
    return y
end

"""
```julia
safe_mul!(dest, A, src, scratch=false) -> dest
```

overwrite `dest` with the result of applying operator `A` to `src` and
returns `dest`.  Unless `scratch` is true, it is guaranteed that `src` is
preserved which may involve making a temporary copy of it.

See also [`safe_mul`](@ref).

"""
function safe_mul!(dest::StridedArray{Complex{T},N},
                   A::cFFTWPlan{Complex{T},K,inplace,N},
                   src::StridedArray{Complex{T},N},
                   scratch::Bool = false) where {T<:fftwReal,K,inplace,N}
    _safe_mul!(dest, A, src, scratch)
end

function safe_mul!(dest::StridedArray{Complex{T},N},
                   A::rFFTWPlan{T,K,inplace,N},
                   src::StridedArray{T,N},
                   scratch::Bool = false) where {T<:fftwReal,K,inplace,N}
    _safe_mul!(dest, A, src, scratch)
end

function safe_mul!(dest::StridedArray{T,N},
                   A::rFFTWPlan{Complex{T},K,inplace,N},
                   src::StridedArray{Complex{T},N},
                   scratch::Bool = false) where {T<:fftwReal,K,inplace,N}
    _safe_mul!(dest, A, src, scratch)
end

function _safe_mul!(dest::StridedArray, A::FFTWPlan,
                    src::StridedArray{T,N}, scratch::Bool) where {T,N}
    if scratch || preserves_input(A)
        mul!(dest, A, src)
    else
        mul!(dest, A, copy(src))
    end
    return dest
end

"""
```julia
safe_mul(A, x, scratch=false)
```

yields the result of applying operator `A` to `x`.  Unless `scratch` is
true, it is guaranteed that input `x` is preserved which may involve making
a temporary copy of it.

See also [`safe_mul!`](@ref).

"""
function safe_mul(A::cFFTWPlan{Complex{T},K,inplace,N},
                  x::StridedArray{Complex{T},N},
                  scratch::Bool = false) where {T<:fftwReal,K,inplace,N}
    y = Array{Complex{T},N}(undef, output_size(A))
    safe_mul!(y, A, x, scratch)
end

function safe_mul(A::rFFTWPlan{T,K,inplace,N},
                  x::StridedArray{T,N},
                  scratch::Bool = false) where {T<:fftwReal,K,inplace,N}
    y = Array{Complex{T},N}(undef, output_size(A))
    safe_mul!(y, A, x, scratch)
end

function safe_mul(A::rFFTWPlan{Complex{T},K,inplace,N},
                  x::StridedArray{Complex{T},N},
                  scratch::Bool = false) where {T<:fftwReal,K,inplace,N}
    y = Array{T,N}(undef, output_size(A))
    safe_mul!(y, A, x, scratch)
end

#------------------------------------------------------------------------------
# FFT operator.

"""
```julia
FFTOperator(A) -> F
```

yields an FFT operator suitable for computing the fast Fourier transform of
arrays similar to `A`.  The operator can also be specified by the
real/complex floating-point type of the elements of the arrays to transform
and their dimensions:

```julia
FFTOperator(T, dims) -> F
```

where `T` is one of `Float64`, `Float32` (for a real-complex FFT),
`Complex{Float64}`, `Complex{Float32}` (for a complex-complex FFT) and
`dims` gives the dimensions of the arrays to transform (by the `Direct` or
`InverseAdjoint` operation).

The interest of creating such an operator is that it caches the ressources
necessary for fast computation of the FFT and can be therefore *much*
faster than calling `fft`, `rfft`, `ifft`, etc.  This is especially true on
small arrays.  Keywords `flags` and `timelimit` may be used to specify
planning options and time limit to create the FFT plans (see
http://www.fftw.org/doc/Planner-Flags.html).  The defaults are
`flags=FFTW.MEASURE` and no time limit.

An instance of `FFTOperator` is a linear mapping which can be used as any
other mapping:

```julia
F*x     # yields the FFT of x
F'*x    # yields the adjoint FFT applied to x, that is the backward FFT of x
F\\x     # yields the inverse FFT of x
```

See also: [`fft`](@ref), [`plan_fft`](@ref), [`bfft`](@ref),
          [`plan_bfft`](@ref), [`rfft`](@ref), [`plan_rfft`](@ref),
          [`brfft`](@ref), [`plan_brfft`](@ref).

"""
struct FFTOperator{T<:fftwNumber,  # element type of input
                   N,              # number of dimensions
                   C<:fftwComplex, # element type of output
                   F<:Plan{T},     # type of forward plan
                   B<:Plan{C}      # type of backward plan
                   } <: LinearMapping
    ncols::Int             # number of input elements
    inpdims::NTuple{N,Int} # input dimensions
    outdims::NTuple{N,Int} # output dimensions
    forward::F             # plan for forward transform
    backward::B            # plan for backward transform
end

# Real-to-complex FFT.
function FFTOperator(::Type{T},
                     dims::NTuple{N,Int};
                     timelimit::Real = FFTW.NO_TIMELIMIT,
                     flags::Integer = FFTW.MEASURE) where {T<:fftwReal,N}
    # Check arguments and build dimension list of the result of the forward
    # real-to-complex (r2c) transform.
    planning = check_flags(flags)
    ncols = check_size(dims)
    zdims = rfftdims(dims)

    # Compute the plans with suitable FFTW flags.  The forward transform (r2c)
    # must preserve its input, while the backward transform (c2r) may destroy
    # it (in fact there are no input-preserving algorithms for
    # multi-dimensional c2r transforms implemented in FFTW, see
    # http://www.fftw.org/doc/Planner-Flags.html).
    forward = plan_rfft(Array{T}(undef, dims);
                        flags = (planning | FFTW.PRESERVE_INPUT),
                        timelimit = timelimit)
    backward = plan_brfft(Array{Complex{T}}(undef, zdims), dims[1];
                          flags = (planning | FFTW.DESTROY_INPUT),
                          timelimit = timelimit)

    # Build operator.
    F = typeof(forward)
    B = typeof(backward)
    return FFTOperator{T,N,Complex{T},F,B}(ncols, dims, zdims,
                                           forward, backward)
end

# Complex-to-complex FFT.
function FFTOperator(::Type{T},
                     dims::NTuple{N,Int};
                     timelimit::Real = FFTW.NO_TIMELIMIT,
                     flags::Integer = FFTW.MEASURE) where {T<:fftwComplex,N}
    # Check arguments.  The input and output of the complex-to-complex
    # transform have the same dimensions.
    planning = check_flags(flags)
    ncols = check_size(dims)
    temp = Array{T}(undef, dims)

    # Compute the plans with suitable FFTW flags.  For maximum efficiency, the
    # transforms are always applied in-place and thus cannot preserve their
    # inputs.
    forward = plan_fft!(temp; flags = (planning | FFTW.DESTROY_INPUT),
                        timelimit = timelimit)
    backward = plan_bfft!(temp; flags = (planning | FFTW.DESTROY_INPUT),
                          timelimit = timelimit)

    # Build operator.
    F = typeof(forward)
    B = typeof(backward)
    return FFTOperator{T,N,T,F,B}(ncols, dims, dims, forward, backward)
end

@callable FFTOperator

# Constructor for dimensions not specified as a tuple.
FFTOperator(T::Type{<:fftwNumber}, dims::Integer...; kwds...) =
    FFTOperator(T, dims; kwds...)

# The following 2 definitions are needed to avoid ambiguities.
FFTOperator(T::Type{<:fftwReal}, dims::Tuple{Vararg{Integer}}; kwds...) =
    FFTOperator(T, to_size(dims); kwds...)
FFTOperator(T::Type{<:fftwComplex}, dims::Tuple{Vararg{Integer}}; kwds...) =
    FFTOperator(T, to_size(dims); kwds...)

# Constructor for transforms applicable to a given array.
FFTOperator(A::DenseArray{T,N}; kwds...) where {T<:fftwNumber,N} =
    FFTOperator(T, size(A); kwds...)

# Traits:
MorphismType(::FFTOperator{<:Complex}) = Endomorphism()

ncols(A::FFTOperator) = A.ncols
ncols(A::Adjoint{<:FFTOperator}) = ncols(unveil(A))
ncols(A::Inverse{<:FFTOperator}) = ncols(unveil(A))
ncols(A::InverseAdjoint{<:FFTOperator}) = ncols(unveil(A))

input_size(A::FFTOperator) = A.inpdims # FIXME: input_size(A.forward)
input_size(A::FFTOperator, i::Integer) = get_dimension(input_size(A), i)
output_size(A::FFTOperator) = A.outdims
output_size(A::FFTOperator, i::Integer) = get_dimension(output_size(A), i)
input_ndims(A::FFTOperator{T,N,C}) where {T,N,C} = N
output_ndims(A::FFTOperator{T,N,C}) where {T,N,C} = N
input_eltype(A::FFTOperator{T,N,C}) where {T,N,C} = T
output_eltype(A::FFTOperator{T,N,C}) where {T,N,C} = C

# 2 FFT operators can be considered the same if they operate on arguments with
# the same element type and the same dimensions.  If the types do not match,
# the matching method is the one which return false, so it is only needed to
# implement the method for two arguments with the same types (omitting the type
# of the plans as it is irrelevant here).
identical(A::FFTOperator{T,N,C}, B::FFTOperator{T,N,C}) where {T,N,C} =
    (input_size(A) == input_size(B))

show(io::IO, A::FFTOperator) = print(io, "FFT")

# Impose the following simplifying rules:
#     inv(F) = n\F'
#     ==> F⋅F' = F'⋅F = n⋅Id
#     ==> inv(F⋅F') = inv(F'⋅F) = inv(F)⋅inv(F') = inv(F')⋅inv(F) = n\Id
*(A::Adjoint{F}, B::F) where {F<:FFTOperator} =
    (identical(unveil(A), B) ? ncols(A)*Id : compose(A, B))
*(A::F, B::Adjoint{F}) where {F<:FFTOperator} =
    (identical(A, unveil(B)) ? ncols(A)*Id : compose(A, B))
*(A::InverseAdjoint{F}, B::Inverse{F}) where {F<:FFTOperator} =
    (identical(unveil(A), unveil(B)) ? (1//ncols(A))*Id : compose(A, B))
*(A::Inverse{F}, B::InverseAdjoint{F}) where {F<:FFTOperator} =
    (identical(unveil(A), unveil(B)) ? (1//ncols(A))*Id : compose(A, B))

function vcreate(P::Type{<:Union{Direct,InverseAdjoint}},
                 A::FFTOperator{T,N,C},
                 x::DenseArray{T,N},
                 scratch::Bool) where {T,N,C}
    vcreate(Direct, A.forward, x, scratch)
end

function vcreate(P::Type{<:Union{Adjoint,Inverse}},
                 A::FFTOperator{T,N,C},
                 x::DenseArray{C,N},
                 scratch::Bool) where {T,N,C}
    vcreate(Direct, A.backward, x, scratch)
end

#
# In principle, FFTW plans can be applied to strided arrays (StridedArray) but
# this imposes that the arguments have the same strides.  So for now, we choose
# to restrict arguments to arrays with contiguous elements (DenseArray).
#

function apply!(α::Number,
                ::Type{Direct},
                A::FFTOperator{T,N,C},
                x::DenseArray{T,N},
                scratch::Bool,
                β::Number,
                y::DenseArray{C,N}) where {T,N,C}
    return apply!(α, Direct, A.forward, x, scratch, β, y)
end

function apply!(α::Number,
                ::Type{Adjoint},
                A::FFTOperator{T,N,C},
                x::DenseArray{C,N},
                scratch::Bool,
                β::Number,
                y::DenseArray{T,N}) where {T,N,C}
    return apply!(α, Direct, A.backward, x, scratch, β, y)
end

function apply!(α::Number,
                ::Type{Inverse},
                A::FFTOperator{T,N,C},
                x::DenseArray{C,N},
                scratch::Bool,
                β::Number,
                y::DenseArray{T,N}) where {T,N,C}
    return apply!(α/ncols(A), Direct, A.backward, x, scratch, β, y)
end

function apply!(α::Number,
                ::Type{InverseAdjoint},
                A::FFTOperator{T,N,C},
                x::DenseArray{T,N},
                scratch::Bool,
                β::Number,
                y::DenseArray{C,N}) where {T,N,C}
    return apply!(α/ncols(A), Direct, A.forward, x, scratch, β, y)
end

#------------------------------------------------------------------------------
# Circulant convolution.

struct CirculantConvolution{T<:fftwNumber,N,
                            C<:fftwComplex,
                            F<:Plan{T},
                            B<:Plan{C}} <: LinearMapping
    dims::NTuple{N,Int}  # input/output dimensions
    zdims::NTuple{N,Int} # complex dimensions
    mtf::Array{C,N}      # modulation transfer function
    forward::F           # plan for forward transform
    backward::B          # plan for backward transform
end

@callable CirculantConvolution

# Traits:
MorphismType(::CirculantConvolution) = Endomorphism()

# Basic methods for a linear operator on Julia's arrays.
input_size(H::CirculantConvolution) = H.dims
output_size(H::CirculantConvolution) = H.dims
input_size(H::CirculantConvolution, i::Integer) = get_dimension(H.dims, i)
output_size(H::CirculantConvolution, i::Integer) = get_dimension(H.dims, i)
input_ndims(H::CirculantConvolution{T,N}) where {T,N} = N
output_ndims(H::CirculantConvolution{T,N}) where {T,N} = N
input_eltype(H::CirculantConvolution{T,N}) where {T,N} = T
output_eltype(H::CirculantConvolution{T,N}) where {T,N} = T

# Basic methods for an array.
Base.eltype(H::CirculantConvolution{T,N}) where {T,N} = T
Base.size(H::CirculantConvolution{T,N}) where {T,N} =
    ntuple(i -> H.dims[(i ≤ N ? i : i - N)], 2*N)
Base.size(H::CirculantConvolution{T,N}, i::Integer) where {T,N} =
    (i < 1 ? bad_dimension_index() : i ≤ N ? H.dims[i] : i ≤ 2N ? H.dims[i-N] : 1)
Base.ndims(H::CirculantConvolution{T,N}) where {T,N} = 2*N

"""
# Circulant convolution operator

The circulant convolution operator `H` is defined by:

```julia
H  = (1/n)*F'*Diag(mtf)*F
```

with `n` the number of elements, `F` the discrete Fourier transform operator
and `mtf` the modulation transfer function.

The operator `H` can be created by:

```julia
H = CirculantConvolution(psf; flags=FFTW.MEASURE, timelimit=Inf, shift=false)
```

where `psf` is the point spread function (PSF).  Note that the PSF is assumed
to be centered according to the convention of the discrete Fourier transform.
You may use `ifftshift` or the keyword `shift` if the PSF is geometrically
centered:

```julia
H = CirculantConvolution(ifftshift(psf))
H = CirculantConvolution(psf, shift=true)
```

The following keywords can be specified:

* `shift` (`false` by default) indicates whether to apply `ifftshift` to `psf`.

* `normalize` (`false` by default) indicates whether to divide `psf` by the sum
  of its values.  This keyword is only available for real-valued PSF.

* `flags` is a bitwise-or of FFTW planner flags, defaulting to `FFTW.MEASURE`.
  If the operator is to be used many times (as in iterative methods), it is
  recommended to use at least `flags=FFTW.MEASURE` (the default) which
  generally yields faster transforms compared to `flags=FFTW.ESTIMATE`.

* `timelimit` specifies a rough upper bound on the allowed planning time, in
  seconds.

The operator can be used as a regular linear operator: `H(x)` or `H*x` to
compute the convolution of `x` and `H'(x)` or `H'*x` to apply the adjoint of
`H` to `x`.

For a slight improvement of performances, an array `y` to store the result of
the operation can be provided:

```julia
apply!(y, [P=Direct,] H, x) -> y
apply!(y, H, x)
apply!(y, H', x)
```

If provided, `y` must be at a different memory location than `x`.

""" CirculantConvolution

function CirculantConvolution(psf::AbstractArray{T,N};
                              kwds...) where {T<:fftwNumber,N}
    CirculantConvolution(copy(psf); kwds...)
end

# Create a circular convolution operator for real arrays.
function CirculantConvolution(psf::DenseArray{T,N};
                              flags::Integer = FFTW.MEASURE,
                              normalize::Bool = false,
                              shift::Bool = false,
                              kwds...) where {T<:fftwReal,N}
    # Check arguments and compute dimensions.
    planning = check_flags(flags)
    n = length(psf)
    dims = size(psf)
    zdims = rfftdims(dims)

    # Allocate array for the scaled MTF, this array also serves as a workspace
    # for planning operations which may destroy their input.
    mtf = Array{Complex{T}}(undef, zdims)

    # Compute the plans with suitable FFTW flags.  The forward transform (r2c)
    # must preserve its input, while the backward transform (c2r) may destroy
    # it (in fact there are no input-preserving algorithms for
    # multi-dimensional c2r transforms).
    forward = safe_plan_rfft(psf; flags = (planning | FFTW.PRESERVE_INPUT), kwds...)
    backward = plan_brfft(mtf, dims[1]; flags = (planning | FFTW.DESTROY_INPUT),
                          kwds...)

    # Compute the scaled MTF *after* computing the plans.
    mul!(mtf, forward, (shift ? ifftshift(psf) : psf))
    if normalize
        sum = mtf[1]
        sum <= 0 && bad_argument("cannot normalize: sum(PSF) ≤ 0")
        sum != 1 && vscale!(mtf, 1/sum)
    end

    # Build operator.
    F = typeof(forward)
    B = typeof(backward)
    CirculantConvolution{T,N,Complex{T},F,B}(dims, zdims, mtf, forward, backward)
end

# Create a circular convolution operator for complex arrays (see
# docs/convolution.md for explanations).
function CirculantConvolution(psf::DenseArray{T,N};
                              flags::Integer = FFTW.MEASURE,
                              normalize::Bool = false,
                              shift::Bool = false,
                              kwds...) where {T<:fftwComplex,N}
    # Check arguments and get dimensions.
    @certify normalize == false "normalizing a complex PSF has no sense"
    planning = check_flags(flags)
    n = length(psf)
    dims = size(psf)

    # Allocate array for the scaled MTF, this array also serves as a workspace
    # for planning operations which may destroy their input.
    mtf = Array{T}(undef, dims)

    # Compute the plans with FFTW flags suitable for out-of-place forward
    # transform and in-place backward transform.
    forward = plan_fft(mtf; flags = (planning | FFTW.PRESERVE_INPUT), kwds...)
    backward = plan_bfft!(mtf; flags = (planning | FFTW.DESTROY_INPUT), kwds...)

    # Compute the MTF *after* computing the plans.
    mul!(mtf, forward, (shift ? ifftshift(psf) : psf))

    # Build the operator.
    F = typeof(forward)
    B = typeof(backward)
    CirculantConvolution{T,N,T,F,B}(dims, dims, mtf, forward, backward)
end

"""

`safe_plan_rfft(x; kwds...)` yields a FFTW plan for computing the real to
complex fast Fourier transform of `x`.  This method is the same as `plan_rfft`
except that it makes sure that `x` is preserved.

"""
function safe_plan_rfft(x::AbstractArray{T,N}; flags::Integer = FFTW.MEASURE,
                        kwds...) where {T<:fftwReal,N}
    planning = (flags & PLANNING)
    if isa(x, StridedArray) && (planning == FFTW.ESTIMATE ||
                                planning == FFTW.WISDOM_ONLY)
        return plan_rfft(x; flags=flags, kwds...)
    else
       return plan_rfft(Array{T}(undef, size(x)); flags=flags, kwds...)
    end
end

function vcreate(::Type{<:Operations},
                 H::CirculantConvolution{T,N},
                 x::AbstractArray{T,N},
                 scratch::Bool) where {T<:fftwNumber,N}
    return Array{T,N}(undef, H.dims)
end

function apply!(α::Number,
                P::Type{<:Union{Direct,Adjoint}},
                H::CirculantConvolution{Complex{T},N,Complex{T}},
                x::AbstractArray{Complex{T},N},
                scratch::Bool,
                β::Number,
                y::AbstractArray{Complex{T},N}) where {T<:fftwReal,N}
    @certify !Base.has_offset_axes(x, y)
    if α == 0
        @certify size(y) == H.dims
        vscale!(y, β)
    else
        n = length(x)
        if β == 0
            # Use y as a workspace.
            mul!(y, H.forward, x) # out-of-place forward FFT of x in y
            _apply!(y, α/n, P, H.mtf) # in-place multiply y by mtf/n
            mul!(y, H.backward, y) # in-place backward FFT of y
        else
            # Must allocate a workspace.
            z = Array{Complex{T}}(undef, H.zdims) # allocate temporary
            mul!(z, H.forward, x) # out-of-place forward FFT of x in z
            _apply!(z, α/n, P, H.mtf) # in-place multiply z by mtf/n
            mul!(z, H.backward, z) # in-place backward FFT of z
            vcombine!(y, 1, z, β, y)
        end
    end
    return y
end

function apply!(α::Number,
                P::Type{<:Union{Direct,Adjoint}},
                H::CirculantConvolution{T,N,Complex{T}},
                x::AbstractArray{T,N},
                scratch::Bool,
                β::Number,
                y::AbstractArray{T,N}) where {T<:fftwReal,N}
    @certify !Base.has_offset_axes(x, y)
    if α == 0
        @certify size(y) == H.dims
        vscale!(y, β)
    else
        n = length(x)
        z = Array{Complex{T}}(undef, H.zdims) # allocate temporary
        mul!(z, H.forward, x) # out-of-place forward FFT of x in z
        _apply!(z, α/n, P, H.mtf) # in-place multiply z by mtf/n
        if β == 0
            mul!(y, H.backward, z) # out-of-place backward FFT of z in y
        else
            w = Array{T}(undef, H.dims) # allocate another temporary
            mul!(w, H.backward, z) # out-of-place backward FFT of z in y
            vcombine!(y, 1, w, β, y)
        end
    end
    return y
end

"""
```julia
_apply!(arr, α, P, mtf)
```

stores in `arr` the elementwise multiplication of `arr` by `α*mtf` if `P` is
`Direct` or by `α*conj(mtf)` if `P` is `Adjoint`.  An error is thrown if the
arrays do not have the same dimensions.  It is assumed that `α ≠ 0`.

"""
function _apply!(arr::AbstractArray{Complex{T},N},
                 α::Number, ::Type{Direct},
                 mtf::AbstractArray{Complex{T},N}) where {T,N}
    @certify axes(arr) == axes(mtf)
    if α == 1
        @inbounds @simd for i in eachindex(arr, mtf)
            arr[i] *= mtf[i]
        end
    else
        alpha = promote_multiplier(α, T)
        @inbounds @simd for i in eachindex(arr, mtf)
            arr[i] *= alpha*mtf[i]
        end
    end
end

function _apply!(arr::AbstractArray{Complex{T},N},
                 α::Number, ::Type{Adjoint},
                 mtf::AbstractArray{Complex{T},N}) where {T,N}
    @certify axes(arr) == axes(mtf)
    if α == 1
        @inbounds @simd for i in eachindex(arr, mtf)
            arr[i] *= conj(mtf[i])
        end
    else
        alpha = promote_multiplier(α, T)
        @inbounds @simd for i in eachindex(arr, mtf)
            arr[i] *= alpha*conj(mtf[i])
        end
    end
end

#------------------------------------------------------------------------------
# Utilities.

"""

`check_flags(flags)` checks whether `flags` is an allowed bitwise-or
combination of FFTW planner flags (see
http://www.fftw.org/doc/Planner-Flags.html) and returns the filtered flags.

"""
function check_flags(flags::Integer)
    planning = flags & PLANNING
    flags == planning ||
        bad_argument("only FFTW planning flags can be specified")
    return UInt32(planning)
end

"""

`get_dimension(dims, i)` yields the `i`-th dimension in tuple of integers
`dims`.  Like for broadcasting rules, it is assumed that the length of
all dimensions after the last one are equal to 1.

"""
get_dimension(dims::NTuple{N,Int}, i::Integer) where {N} =
    (i < 1 ? bad_dimension_index() : i ≤ N ? dims[i] : 1)
# FIXME: should be in ArrayTools
bad_dimension_index() = error("invalid dimension index")

"""
```julia
goodfftdim(len)
```

yields the smallest integer which is greater or equal `len` and which is a
multiple of powers of 2, 3 and/or 5.  If argument is an array dimesion list
(i.e. a tuple of integers), a tuple of good FFT dimensions is returned.

Also see: [`goodfftdims`](@ref), [`rfftdims`](@ref), [`FFTOperator`](@ref).

"""
goodfftdim(len::Integer) = goodfftdim(Int(len))
goodfftdim(len::Int) = nextprod([2,3,5], len)

"""
```julia
goodfftdims(dims)
```

yields a list of dimensions suitable for computing the FFT of arrays whose
dimensions are `dims` (a tuple or a vector of integers).

Also see: [`goodfftdim`](@ref), [`rfftdims`](@ref), [`FFTOperator`](@ref).

"""
goodfftdims(dims::Integer...) = map(goodfftdim, dims)
goodfftdims(dims::Union{AbstractVector{<:Integer},Tuple{Vararg{Integer}}}) =
    map(goodfftdim, dims)

"""
```julia
rfftdims(dims)
```

yields the dimensions of the complex array produced by a real-complex FFT of a
real array of size `dims`.

Also see: [`goodfftdim`](@ref), [`FFTOperator`](@ref).

"""
rfftdims(dims::Integer...) = rfftdims(dims)
rfftdims(dims::NTuple{N,Integer}) where {N} =
    ntuple(d -> (d == 1 ? (Int(dims[d]) >>> 1) + 1 : Int(dims[d])), Val(N))
# Note: The above version is equivalent but much faster than
#     ((dims[1] >>> 1) + 1, dims[2:end]...)
# which is not optimized out by the compiler.

"""
### Generate Discrete Fourier Transform frequency indexes or frequencies

Syntax:

```julia
k = fftfreq(dim)
f = fftfreq(dim, step)
```

With a single argument, the function returns a vector of `dim` values set with
the frequency indexes:

```
k = [0, 1, 2, ..., n-1, -n, ..., -2, -1]   if dim = 2*n
k = [0, 1, 2, ..., n,   -n, ..., -2, -1]   if dim = 2*n + 1
```

depending whether `dim` is even or odd.  These rules are compatible to what is
assumed by `fftshift` (which to see) in the sense that:

```
fftshift(fftfreq(dim)) = [-n, ..., -2, -1, 0, 1, 2, ...]
```

With two arguments, `step` is the sample spacing in the direct space and the
result is a floating point vector with `dim` elements set with the frequency
bin centers in cycles per unit of the sample spacing (with zero at the start).
For instance, if the sample spacing is in seconds, then the frequency unit is
cycles/second.  This is equivalent to:

```
fftfreq(dim)/(dim*step)
```

See also: [`FFTOperator`](@ref), [`fftshift`](@ref).

"""
function fftfreq(_dim::Integer)
    dim = Int(_dim)
    n = div(dim, 2)
    f = Array{Int}(undef, dim)
    @inbounds begin
        for k in 1:dim-n
            f[k] = k - 1
        end
        for k in dim-n+1:dim
            f[k] = k - (1 + dim)
        end
    end
    return f
end

function fftfreq(_dim::Integer, step::Real)
    dim = Int(_dim)
    scl = Cdouble(1/(dim*step))
    n = div(dim, 2)
    f = Array{Cdouble}(undef, dim)
    @inbounds begin
        for k in 1:dim-n
            f[k] = (k - 1)*scl
        end
        for k in dim-n+1:dim
            f[k] = (k - (1 + dim))*scl
        end
    end
    return f
end

end # module


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 3: FKBP12-mTOR
"""

# ╔═╡ 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 = 5
	nsamples = 2
	ncells = 10

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

	experimentresults = load(datapath)["experimentresults"]

	medianmeasurements = Array{Float64,3}(undef, ncells, nreplicates, nsamples)
	montecarlomeasurements = Array{Float64,3}(undef, ncells, nreplicates, nsamples)
	positivecontrolmontecarlomeasurements = Array{Float64,3}(undef, ncells, nreplicates, nsamples)
	negativecontrolmontecarlomeasurements = Array{Float64,3}(undef, ncells, nreplicates, nsamples)

	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[1][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] = samplemedianresults
		montecarlomeasurements[:, :, i] = samplelessthan10results
		positivecontrolmontecarlomeasurements[:, :, i] = samplepositivecontrollessthan10results
		negativecontrolmontecarlomeasurements[:, :, i] = samplenegativecontrollessthan10results
	end

	normalizedmontecarlomeasurements = (montecarlomeasurements .- negativecontrolmontecarlomeasurements) ./ (positivecontrolmontecarlomeasurements .- negativecontrolmontecarlomeasurements)
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 = -Rap, 2 = +Rap)
"""

# ╔═╡ 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
"""

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

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

	plot(p1, p2, layout = grid(1, 2), legend = :none, plot_title = "FKBP12-mTOR")
end

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

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

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

# ╔═╡ f55402f2-6a25-11eb-18df-59b6da713a4d
mediansflat = [medianmeasurements[:,:,1] medianmeasurements[:,:,2]]

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

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

# ╔═╡ 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 one quite severe. But since an ANOVA will be used, must be careful not to delete points from a balanced design. Therefore, no outliers will be removed.

However, those two outliers in [1] Replicate 1, 3, and 5, and [2] Replicate 2 are so extreme that I will windsorize them.
"""

# ╔═╡ 6da539c0-5907-11eb-0f68-c7a5f1ee90ca
begin
	medianmeasurementsw = copy(medianmeasurements)
	mediansflatw = copy(mediansflat)
	medianmeasurementsw[4,1,1] = medianmeasurementsw[5,1,1]
	mediansflatw[4,1] = mediansflatw[5,1]
	medianmeasurementsw[2,3,1] = medianmeasurementsw[5,3,1]
	mediansflatw[2,3] = mediansflatw[5,3]
	medianmeasurementsw[9,5,1] = medianmeasurementsw[4,5,1]
	mediansflatw[9,5] = mediansflatw[4,5]
	medianmeasurementsw[4,2,2] = medianmeasurementsw[10,2,2]
	mediansflatw[4,7] = mediansflatw[10,7]
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 fairly 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(mediansflatw[:,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 values were within 1.96 (5% level).
"""

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

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

# ╔═╡ 56cb4410-6f0a-11eb-2680-b598165ea2f1
std(medianresult.crossedcellmeans, dims = 2)

# ╔═╡ fd2d4a40-57ab-11eb-2688-53b009a37680
md"""
#### Interpret main effects
Rapamycin caused a detectable but small decrease in median distance.
"""

# ╔═╡ 895181b0-59f2-11eb-3d4e-1d92833b0e72
md"""
#### Figure
"""

# ╔═╡ 95fd3850-59f2-11eb-0914-a999ed3b9e5f
let
	import StatsPlots.mm
	fkbp12_mtor_median = [medianmeasurementsw[:,:,1] |> vec; medianmeasurementsw[:,:,2] |> vec]
	groups = repeat([1,2], inner = 50)
	boxplot(groups, fkbp12_mtor_median, outliers=false,
			legend=:none,
			seriescolor=:white,
		    left_margin=2mm,
			size=(256, 512),
			line=(2, 0.75),
		    xgrid=:none,
			xaxis=("Rapamycin", (1:2, ["-Rap", "+Rap"])),
			yaxis=("Median distance (nm)", [0, 600]))
	fkbp12_mtor_median_means = dropdims(mean(medianmeasurementsw, dims=1), dims=1)
	dotplot!([1,1,1,1,1,2,2,2,2,2], fkbp12_mtor_median_means |> vec, mode = :none, label="", marker=(4, 0.75, :rect, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray]), stroke(0)))
	dotplot!(groups, fkbp12_mtor_median, mode = :density, label="", marker=(2, 0.5, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray], inner=10), stroke(0)))
end

# ╔═╡ 1715afae-6aee-11eb-3c61-a3902bf7e827
md"""Create paper-quality version and save it:"""

# ╔═╡ 15ee46b0-6aee-11eb-0dce-e39b244d6d3c
let
	import StatsPlots.mm
	fkbp12_mtor_median = [medianmeasurementsw[:,:,1] |> vec; medianmeasurementsw[:,:,2] |> vec]
	groups = repeat([1,2], inner = 50)
	boxplot(groups, fkbp12_mtor_median, outliers=false,
			legend=:none,
		    left_margin=25mm,
			seriescolor=:white,
		    guidefontsize=52,
            tickfontsize=48,
			size=(1024, 2048),
			line=(8, 0.75),
		    xgrid=:none,
			xaxis=("Rapamycin", (1:2, ["-Rap", "+Rap"])),
			yaxis=("Median distance (nm)", [0, 600]))
	fkbp12_mtor_median_means = dropdims(mean(medianmeasurementsw, dims=1), dims=1)
	dotplot!([1,1,1,1,1,2,2,2,2,2], fkbp12_mtor_median_means |> vec, mode = :none, label="", marker=(12, 0.75, :rect, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray]), stroke(0)))
	dotplot!(groups, fkbp12_mtor_median, mode = :density, label="", marker=(8, 0.5, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray], inner=10), stroke(0)))
	savefig(joinpath(outputdir, "fkbp12-mTOR-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], xaxis = ("Replicates", [1, 2, 3, 4, 5]), yaxis = ("Fraction associated"))
	p2 = boxplot(montecarlomeasurements[:, :, 2], xaxis = ("Replicates", [1, 2, 3, 4, 5]), yaxis = ("Fraction associated"))

	plot(p1, p2, layout = grid(1, 2), legend = :none, plot_title = "FKBP12-mTOR")
end

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

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

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

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

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

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

# ╔═╡ 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 [1] Replicate 3 and 4 outliers.
"""

# ╔═╡ 0498eae0-59e7-11eb-028c-47a1e64c4e5e
begin
	montecarlomeasurementsw = copy(montecarlomeasurements)
	montecarloflatw = copy(montecarloflat)
	montecarlomeasurementsw[8,3,1] = montecarlomeasurementsw[7,3,1]
	montecarlomeasurementsw[8,4,1] = montecarlomeasurementsw[6,4,1]
	montecarloflatw[8,3] = montecarloflatw[7,3]
	montecarloflatw[8,4] = montecarloflatw[6,4]
end;

# ╔═╡ 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 somewhat 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, (1,3,2)),
	[fixed, subject],
    factornames = ["Rapamycin", "Replicate"],
)

# ╔═╡ 2553a3a0-57a9-11eb-247d-a161d85d4410
md"""
#### Interpret main effects
Rapamycin has a large effect and is a bit above the significance threshold.
"""

# ╔═╡ 38571c12-59f3-11eb-337f-13164dbd552b
md"""
#### Figure
"""

# ╔═╡ 475acc70-59f3-11eb-08a1-99e272fbc933
let
	import StatsPlots.mm
	fkbp12_mtor_montecarlo = [montecarlomeasurementsw[:,:,1] |> vec; montecarlomeasurementsw[:,:,2] |> vec]
	groups = repeat([1,2], inner = 50)
	boxplot(groups, fkbp12_mtor_montecarlo, outliers=false,
			legend=:none,
			seriescolor=:white,
		    left_margin=3mm,
			line=(2, 0.75),
		    size=(256, 512),
		    xgrid=:none,
			xaxis=("Rapamycin", (1:2, ["-Rap", "+Rap"])),
			yaxis=("Fraction associated", [0, 0.25]))
	fkbp12_mtor_montecarlo_means = dropdims(mean(montecarlomeasurementsw, dims=1), dims=1)
	dotplot!([1,1,1,1,1,2,2,2,2,2], fkbp12_mtor_montecarlo_means |> vec, mode = :none, label="", marker=(4, 0.75, :rect, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray]), stroke(0)))
	dotplot!(groups, fkbp12_mtor_montecarlo, mode = :density, label="", marker=(2, 0.5, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray], inner=10), stroke(0)))
end

# ╔═╡ 9f7009e0-6aef-11eb-0a49-a15c24bee94d
md"""Create paper-quality version and save it:"""

# ╔═╡ c3d837e0-6aee-11eb-1840-053a9dca4891
let
	import StatsPlots.mm
	fkbp12_mtor_montecarlo = [montecarlomeasurementsw[:,:,1] |> vec; montecarlomeasurementsw[:,:,2] |> vec]
	groups = repeat([1,2], inner = 50)
	boxplot(groups, fkbp12_mtor_montecarlo, outliers=false,
			legend=:none,
		    left_margin=25mm,
			seriescolor=:white,
		    guidefontsize=52,
            tickfontsize=48,
			size=(1024, 2048),
			line=(8, 0.75),
		    xgrid=:none,
			xaxis=("Rapamycin", (1:2, ["-Rap", "+Rap"])),
			yaxis=("Fraction associated", [0, 0.25]))
	fkbp12_mtor_montecarlo_means = dropdims(mean(montecarlomeasurementsw, dims=1), dims=1)
	dotplot!([1,1,1,1,1,2,2,2,2,2], fkbp12_mtor_montecarlo_means |> vec, mode = :none, label="", marker=(12, 0.75, :rect, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray]), stroke(0)))
	dotplot!(groups, fkbp12_mtor_montecarlo, mode = :density, label="", marker=(8, 0.5, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray], inner=10), stroke(0)))
	savefig(joinpath(outputdir, "fkbp12-mTOR-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], xaxis = ("Replicates", [1, 2, 3, 4, 5]), yaxis = ("Fraction bound"))
	p2 = boxplot(normalizedmontecarlomeasurements[:, :, 2], xaxis = ("Replicates", [1, 2, 3, 4, 5]), yaxis = ("Fraction bound"))

	plot(p1, p2, layout = grid(1, 2), legend = :none, plot_title = "FKBP12-mTOR")
end

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

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

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

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

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

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

However, I will Windsorize the worst outliers of [1] Replicate 4 and 5.
"""

# ╔═╡ 49094280-66a6-11eb-1261-7fe254f02b6d
begin
	normalizedmontecarlomeasurementsw = copy(normalizedmontecarlomeasurements)
	normmontecarloflatw = copy(normmontecarloflat)
	normalizedmontecarlomeasurementsw[8,4,1] = normalizedmontecarlomeasurementsw[5,4,1]
	normalizedmontecarlomeasurementsw[6,5,1] = normalizedmontecarlomeasurementsw[9,5,1]
	normmontecarloflatw[8,4] = normmontecarloflatw[5,4]
	normmontecarloflatw[6,5] = normmontecarloflatw[9,5]
end;

# ╔═╡ 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 somewhat 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(normmontecarloflatw)

# ╔═╡ 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(normmontecarloflatw[:,i]) for i ∈ axes(normmontecarloflatw, 2)]

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

# ╔═╡ 80f22522-575d-11eb-211b-59fe55838280
md"""
From skewness and kurtosis statistics, one was greater in value than 1.96 (5% level), [1] Replicate 4, with skewness 1.78 and kurtosis 2.36.
"""

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

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

# ╔═╡ 408e5d30-5764-11eb-124f-4f2630d1d908
md"""
#### Interpret main effects
Rapamycin had a moderate effect but was not significant (2.13% to 5.11%).
"""

# ╔═╡ 8c39ec6e-5a72-11eb-3b9b-1d45fbdeaa8c
md"""
#### Figure
"""

# ╔═╡ 9cd79dc0-5a72-11eb-12a5-cd1e42c92297
let
	import StatsPlots.mm
	fkbp12_mtor_montecarlo = [normalizedmontecarlomeasurements[:,:,1] |> vec; normalizedmontecarlomeasurements[:,:,2] |> vec]
	groups = repeat([1,2], inner = 50)
	boxplot(groups, fkbp12_mtor_montecarlo, outliers=false,
			legend=:none,
			seriescolor=:white,
		    left_margin=3mm,
			line=(3, 0.75),
		    size=(256, 512),
		    xgrid=:none,
			xaxis=("Rapamycin", (1:2, ["-Rap", "+Rap"])),
			yaxis=("Median distance (nm)", [-0.1,0.2]))
	fkbp12_mtor_montecarlo_means = dropdims(mean(normalizedmontecarlomeasurements, dims=1), dims=1)
	dotplot!([1,1,1,1,1,2,2,2,2,2], fkbp12_mtor_montecarlo_means |> vec, mode = :none, label="", marker=(6, 0.75, :rect, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray]), stroke(0)))
	dotplot!(groups, fkbp12_mtor_montecarlo, mode = :density, label="", marker=(4, 0.5, repeat([:orange, :darkblue, :darkred, :darkgreen, :darkgray], inner=10), stroke(0)))
end

# ╔═╡ f55ff670-6bbd-11eb-3650-8585a350ef2c
md"""
## Example localization maps
"""

# ╔═╡ 8e3320f0-6c24-11eb-1e67-ef196dc98679
let
	# Exp 1 (FKBP12-mTOR)
	# A3
	insetx, insety = [14600, 16000], [16300, 17700]
	r = experimentresults[1][5][1][9]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "1 - A9 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "1 - A9 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

# ╔═╡ b0ecbc82-6c8f-11eb-14ee-41bf2e845adc
let
	# Exp 1 (FKBP12-mTOR)
	# B4
	insetx, insety = [19000, 20400], [14300, 15700]
	r = experimentresults[1][5][2][8]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "1 - B8 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "1 - B8 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

# ╔═╡ eb152440-6caf-11eb-34bd-a9fc47f16a4d
let
	# Plot all
	r = localizationsplot.(experimentresults[1][j][i][k] for k ∈ 1:10 for i ∈ 1:2 for j ∈ 1:5)
	p = plot(r..., size=(2048, 2048), fmt=:png)
	savefig(joinpath(outputdir, "fkbp12-mTOR-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(normmontecarloflatw)

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

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

# ╔═╡ 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
# ╟─f55402f2-6a25-11eb-18df-59b6da713a4d
# ╟─25253910-5516-11eb-0a7a-1d73ed51d4a1
# ╟─d8c96c5e-5517-11eb-0861-056d5e425f04
# ╟─87eb8f70-5518-11eb-1639-f91b4bcb5c9b
# ╠═6da539c0-5907-11eb-0f68-c7a5f1ee90ca
# ╟─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
# ╟─265006c0-55c1-11eb-2ba8-3d628d096704
# ╠═c101d900-57ab-11eb-04ab-77cd4e0ff4b3
# ╠═56cb4410-6f0a-11eb-2680-b598165ea2f1
# ╟─fd2d4a40-57ab-11eb-2688-53b009a37680
# ╟─895181b0-59f2-11eb-3d4e-1d92833b0e72
# ╟─95fd3850-59f2-11eb-0914-a999ed3b9e5f
# ╟─1715afae-6aee-11eb-3c61-a3902bf7e827
# ╠═15ee46b0-6aee-11eb-0dce-e39b244d6d3c
# ╟─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
# ╟─c18cc050-56ad-11eb-261c-4177aa31ca7c
# ╠═0498eae0-59e7-11eb-028c-47a1e64c4e5e
# ╟─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
# ╟─38571c12-59f3-11eb-337f-13164dbd552b
# ╟─475acc70-59f3-11eb-08a1-99e272fbc933
# ╟─9f7009e0-6aef-11eb-0a49-a15c24bee94d
# ╠═c3d837e0-6aee-11eb-1840-053a9dca4891
# ╟─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
# ╠═49094280-66a6-11eb-1261-7fe254f02b6d
# ╟─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
# ╟─408e5d30-5764-11eb-124f-4f2630d1d908
# ╟─8c39ec6e-5a72-11eb-3b9b-1d45fbdeaa8c
# ╠═9cd79dc0-5a72-11eb-12a5-cd1e42c92297
# ╟─f55ff670-6bbd-11eb-3650-8585a350ef2c
# ╠═8e3320f0-6c24-11eb-1e67-ef196dc98679
# ╠═b0ecbc82-6c8f-11eb-14ee-41bf2e845adc
# ╠═eb152440-6caf-11eb-34bd-a9fc47f16a4d
# ╟─9c92b290-56ae-11eb-2595-85845fe95f0e
# ╟─2623d480-56af-11eb-2ab8-0d0afd8c45dc
# ╟─a655a120-56ae-11eb-30be-459d7c3067dc
# ╟─7922d1d0-56b0-11eb-261c-5f87ca587394


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -11.08168115564476 <= q <= 12.629895565298366)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.8013086439641905 + -0.577380120181807*1/(1 + exp(-(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4])))) + -0.9073468716494122*1/(1 + exp(-(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4])))) + -0.34890000179172853*1/(1 + exp(-(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))))))) + 1/(1 + exp(-(-0.6400914085711236 + 0.5999574270201746*1/(1 + exp(-(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4])))) + -0.3302091884595031*1/(1 + exp(-(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4])))) + 0.9262651374928534*1/(1 + exp(-(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))))))) + 1/(1 + exp(-(0.7171327203675486 + 0.9873737246227612*1/(1 + exp(-(-0.15161094843442724 + 0.3901524114721284*$(x[1]) + -0.7620628362671806*$(x[2]) + 0.5531251925343885*$(x[3]) + 0.8877751701397409*$(x[4])))) + -0.03055409136192777*1/(1 + exp(-(-0.9517329762307463 + 0.7324043290359472*$(x[1]) + -0.6011940599882695*$(x[2]) + 0.8132214951552479*$(x[3]) + 0.2016959253369457*$(x[4])))) + 0.20254351925716962*1/(1 + exp(-(0.5700564024677015 + -0.6508434167160906*$(x[1]) + 0.3531203702734036*$(x[2]) + 0.9592602519851474*$(x[3]) + -0.2457804299770685*$(x[4]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module ChainRules

using Base.Broadcast: materialize, materialize!, broadcasted, Broadcasted, broadcastable
using ChainRulesCore
using Compat
using LinearAlgebra
using LinearAlgebra.BLAS
using Random
using RealDot: realdot
using Statistics

# Basically everything this package does is overloading these, so we make an exception
# to the normal rule of only overload via `ChainRulesCore.rrule`.
import ChainRulesCore: rrule, frule

if VERSION < v"1.3.0-DEV.142"
    # In prior versions, the BLAS submodule also exported `dot`, which caused a conflict
    # with its parent module. To get around this, we can simply create a hard binding for
    # the one we want to use without qualification.
    import LinearAlgebra: dot
end

# numbers that we know commute under multiplication
const CommutativeMulNumber = Union{Real,Complex}

include("rulesets/Core/core.jl")

include("rulesets/Base/utils.jl")
include("rulesets/Base/nondiff.jl")
include("rulesets/Base/base.jl")
include("rulesets/Base/fastmath_able.jl")
include("rulesets/Base/evalpoly.jl")
include("rulesets/Base/array.jl")
include("rulesets/Base/arraymath.jl")
include("rulesets/Base/indexing.jl")
include("rulesets/Base/sort.jl")
include("rulesets/Base/mapreduce.jl")

include("rulesets/Statistics/statistics.jl")

include("rulesets/LinearAlgebra/utils.jl")
include("rulesets/LinearAlgebra/blas.jl")
include("rulesets/LinearAlgebra/lapack.jl")
include("rulesets/LinearAlgebra/dense.jl")
include("rulesets/LinearAlgebra/norm.jl")
include("rulesets/LinearAlgebra/matfun.jl")
include("rulesets/LinearAlgebra/structured.jl")
include("rulesets/LinearAlgebra/symmetric.jl")
include("rulesets/LinearAlgebra/factorization.jl")

include("rulesets/Random/random.jl")

end # module