text:
export computestat, computestat!, computestat_parallel, computestat_parallel!,
       PowerSpectrum, CrossSpectrum, Correlation, Coherency,
       Coherence, MeanPhaseDiff, PLV, PPC, PLI, PLI2Unbiased, WPLI,
       WPLI2Debiased, JammalamadakaR, JuppMardiaR,
       HurtadoModulationIndex, UniformScores, JackknifeSurrogates,
       MultiJackknifeSurrogates, jackknife_bias, jackknife_var,
       Jackknife, JackknifeCorrelation, Bootstrap, Permutation, GroupMean

#
# Utilities
#

typealias AbstractVecOrMat{T} Union{AbstractVector{T}, AbstractMatrix{T}}
typealias StridedVecOrMat{T} Union{StridedVector{T}, StridedMatrix{T}}

ntrials(X::AbstractVecOrMat) = size(X, 1)
nchannels(X::AbstractVecOrMat) = size(X, 2)

# Upper triangular symmetric/hermitian multiply, ignoring lower triangular part
Ac_mul_A!{T<:Complex}(out::StridedMatrix{T}, X::StridedVecOrMat{T}) =
    BLAS.herk!('U', 'C', real(one(T)), X, real(zero(T)), out)
At_mul_A!{T}(out::StridedMatrix{T}, X::StridedVecOrMat{T}) =
    BLAS.syrk!('U', 'C', real(one(T)), X, real(zero(T)), out)
Ac_mul_A!{T}(out::StridedMatrix{T}, X::StridedVecOrMat{T}) = At_mul_A!(out, X)

# Normalize to unit circle
function unitnormalize!{T<:Complex}(out::AbstractArray{T}, X::AbstractArray{T})
    length(out) == length(X) || error(DimensionMismatch("work size mismatch"))
    for i = 1:length(X)
        @inbounds out[i] = X[i]/abs(X[i])
    end
    out
end

# Convert uncentered covariance to coherency or coherence
function cov2coh!{T<:Real}(out::Union{AbstractMatrix{Complex{T}}, AbstractMatrix{T}},
                           XXc::Union{AbstractMatrix{Complex{T}}, AbstractMatrix{T}}=out,
                           f=Base.IdFun())
    size(XXc, 2) == size(XXc, 1) || error(ArgumentError("XXc must be square"))
    size(out, 1) == size(out, 2) == size(XXc, 2) ||
        error(DimensionMismatch("output size mismatch"))

    for j = 1:size(XXc, 2)
        out[j, j] = x = 1/sqrt(real(XXc[j, j]))
        for i = 1:j-1
            @inbounds out[i, j] = f(XXc[i, j])*(real(out[i, i])*x)
        end
    end
    for j = 1:size(XXc, 2)
        @inbounds out[j, j] = one(T)
    end
    out
end
function cov2coh!{T<:Real}(out::Union{AbstractMatrix{Complex{T}}, AbstractMatrix{T}},
                           X::AbstractVecOrMat, Y::AbstractVecOrMat,
                           Xtmp::AbstractVecOrMat{T}, Ytmp::AbstractVecOrMat{T},
                           XYc::Union{AbstractMatrix{Complex{T}}, AbstractMatrix{T}}=out,
                           f=Base.IdFun())
    (size(XYc, 1) == length(Xtmp) && size(XYc, 2) == length(Ytmp)) ||
        error(DimensionMismatch("work size mismatch"))
    (size(out, 1) == size(XYc, 1) && size(out, 2) == size(XYc, 2)) ||
        error(DimensionMismatch("output size mismatch"))

    sumabs2!(Xtmp, X)
    sumabs2!(Ytmp, Y)
    for i = 1:length(Xtmp)
        Xtmp[i] = 1/sqrt(Xtmp[i])
    end
    for i = 1:length(Ytmp)
        Ytmp[i] = 1/sqrt(Ytmp[i])
    end

    for j = 1:size(XYc, 2), i = 1:size(XYc, 1)
        @inbounds out[i, j] = f(XYc[i, j])*(Xtmp[i]*Ytmp[j])
    end
    out
end

# Work array for X or Y in two-argument form
cov2coh_work{T}(X::AbstractVector{Complex{T}}) = Array(T, 1)
cov2coh_work{T}(X::AbstractMatrix{Complex{T}}) = Array(T, 1, nchannels(X))
cov2coh_work{T}(X::AbstractVector{T}) = Array(T, 1)
cov2coh_work{T}(X::AbstractMatrix{T}) = Array(T, 1, nchannels(X))

# Fastest possible unsafe submatrix
immutable ConjFun <: Base.Func{1} end
call(::ConjFun, x) = conj(x)

# Check that X and Y have the same number of trials
chkXY(X, Y) = size(X, 1) == size(Y, 1) ||
    error(DimensionMismatch("X and Y must have same number of trials"))
chkout(out, X) = size(out, 1) == size(out, 2) == size(X, 2) ||
    error(DimensionMismatch("output size mismatch"))
chkout(out, X, Y) = (size(out, 1) == size(X, 2) && size(out, 2) == size(Y, 2)) ||
    error(DimensionMismatch("output size mismatch"))
chkinput(out, X) = chkout(out, X)
chkinput(out, X, Y) = (chkXY(X, Y); chkout(out, X, Y))

#
# Types and basics
#

abstract Statistic

abstract PairwiseStatistic <: Statistic
abstract NormalizedPairwiseStatistic{Normalized} <: PairwiseStatistic

# General definitions of allocoutput
Base.eltype{T,S}(s::PairwiseStatistic, X::AbstractArray{T}, Y::AbstractArray{S}) =
    promote_type(eltype(s, X), eltype(s, Y))
allocoutput{T}(t::PairwiseStatistic, X::AbstractVecOrMat{T}) =
    zeros(eltype(t, X), size(X, 2), size(X, 2))
allocoutput{T}(t::PairwiseStatistic, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    zeros(eltype(t, X, Y), size(X, 2), size(Y, 2))
outputdims{R<:PairwiseStatistic}(t::Type{R}) = 2

#
# Diagonal values
#
diagval{T<:Statistic}(::Type{T}) = nothing

#
# Handling of n-d arrays
#
# X is channels x trials x ...
# output is channels x channels x ...
computestat(t::Statistic, X::AbstractArray) =
    computestat!(t, allocoutput(t, X), allocwork(t, X), X)
computestat(t::Statistic, X::AbstractArray, Y::AbstractArray) =
    computestat!(t, allocoutput(t, X, Y), allocwork(t, X, Y), X, Y)
computestat_parallel(t::Statistic, X::AbstractArray) =
    computestat_parallel!(t, allocoutput(t, X), allocwork(t, X), X)
computestat_parallel(t::Statistic, X::AbstractArray, Y::AbstractArray) =
    computestat_parallel!(t, allocoutput(t, X, Y), allocwork(t, X, Y), X, Y)

allocwork{T}(t::Statistic, X::AbstractVecOrMat{T}) =
    error("allocwork($(typeof(t)), X) not defined")
allocwork{T}(t::Statistic, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    error("allocwork($(typeof(t)), X, Y) not defined")
allocwork{T}(t::Statistic, X::AbstractArray{T}) = allocwork(t, unsafe_view(X, :, :, 1))
function allocwork(t::Statistic, X::AbstractArray, Y::AbstractArray)
    allocwork(t, unsafe_view(X, :, :, 1), unsafe_view(Y, :, :, 1))
end

@generated function allocoutput{T,N}(t::Statistic, X::AbstractArray{T,N}, Y::AbstractArray{T,N}=X)
    :(zeros(eltype(t, X, Y), size(X, 2), size(Y, 2), $([:(size(X, $i)) for i = 3:N]...)))
end

computestat!{T}(t::Statistic, out::AbstractArray, work, X::AbstractVecOrMat{T}) =
    error("computestat! not defined for $(typeof(t)) with single input matrix, or incorrect work $(typeof(work)) or output $(typeof(out)) given")
computestat!{T}(t::Statistic, out::AbstractArray, work, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    error("computestat! not defined for $(typeof(t)) with two input matrices, or incorrect work $(typeof(work)) or output $(typeof(out)) given")

function checknd(t::Statistic, out::AbstractArray, X::AbstractArray)
    isempty(X) && throw(ArgumentError("X is empty"))
    trail = Base.trailingsize(X, 3)
    Base.trailingsize(out, outputdims(typeof(t))+1) == trail || throw(DimensionMismatch("output size mismatch"))
    trail
end

function checknd(t::Statistic, out::AbstractArray, X::AbstractArray, Y::AbstractArray)
    (isempty(X) || isempty(Y)) && throw(ArgumentError("X or Y is empty"))
    trail = Base.trailingsize(X, 3)
    Base.trailingsize(Y, 3) == trail || throw(DimensionMismatch("X and Y must have same trailing dimensions"))
    Base.trailingsize(out, outputdims(typeof(t))+1) == trail || throw(DimensionMismatch("output size mismatch"))
    trail
end

# Generate an unsafe_view of the output for the ith independent input
@generated function output_view(t::Statistic, out::AbstractArray, i::Union{Int, UnitRange{Int}})
    :(unsafe_view(out, $(fill(:(:), outputdims(t))...), i))
end

# nd computestat!
function computestat!{T,V}(t::Statistic, out::AbstractArray, work::V, X::AbstractArray{T})
    for i = 1:checknd(t, out, X)
        computestat!(t, output_view(t, out, i), work, unsafe_view(X, :, :, i))
    end
    out
end

function computestat!{T,V}(t::Statistic, out::AbstractArray, work::V, X::AbstractArray{T}, Y::AbstractArray{T})
    for i = 1:checknd(t, out, X, Y)
        computestat!(t, output_view(t, out, i), work, unsafe_view(X, :, :, i), unsafe_view(Y, :, :, i))
    end
    out
end

# parallel nd computestat!
function chunkranges(x::Int, n::Int)
    parts = Array(UnitRange{Int}, n)
    chunklen = div(x, n)
    curind = 1
    for i = 1:n
        chunklen = div(x - curind + 1, n-i+1)
        parts[i] = curind:curind+chunklen-1
        curind = curind+chunklen
    end
    parts
end

function computestat_parallel_complete!(t, out, refs, ranges)
    for i = 1:length(refs)
        res = fetch(refs[i])
        isa(res, Exception) && rethrow(res)
        copy!(output_view(t, out, ranges[i]), res)
    end
    out
end

function computestat_parallel!(t::Statistic, out, work, X::AbstractArray)
    np = max(nprocs()-1, 1)
    p = np == 1 ? (1:1) : 2:np+1
    ranges = chunkranges(checknd(t, out, X), np)
    n = length(ranges)

    refs = [remotecall(p[i], computestat, t, X[:, :, ranges[i]]) for i = 1:length(ranges)]
    computestat_parallel_complete!(t, out, refs, ranges)
end

function computestat_parallel!{T}(t::Statistic, out, work, X::AbstractArray{T}, Y::AbstractArray{T})
    np = max(nprocs()-1, 1)
    p = np == 1 ? (1:1) : 2:np+1
    ranges = chunkranges(checknd(t, out, X, Y), np)
    n = length(ranges)

    refs = [remotecall(p[i], computestat, t, X[:, :, ranges[i]], Y[:, :, ranges[i]]) for i = 1:length(ranges)]
    computestat_parallel_complete!(t, out, refs, ranges)
end

#
# Normalized statistics
#

@generated function normalized(t::NormalizedPairwiseStatistic{false})
    :($(t.name.name){true}())
end

allocwork{T<:Real}(t::NormalizedPairwiseStatistic{false}, X::AbstractVecOrMat{Complex{T}}) =
    (Array(Complex{T}, size(X, 1), size(X, 2)), allocwork(normalized(t), X))
allocwork{T<:Real}(t::NormalizedPairwiseStatistic{false}, X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) =
    (Array(Complex{T}, size(X, 1), size(X, 2)), Array(Complex{T}, size(X, 1), size(X, 2)), allocwork(normalized(t), X, Y))
computestat!{T<:Real}(t::NormalizedPairwiseStatistic{false}, out::AbstractArray,
                      work::Tuple{Matrix{Complex{T}}, Any}, X::AbstractVecOrMat{Complex{T}}) =
    computestat!(normalized(t), out, work[2], unitnormalize!(work[1], X))
computestat!{T<:Real}(t::NormalizedPairwiseStatistic{false}, out::AbstractArray,
                      work::Tuple{Matrix{Complex{T}}, Matrix{Complex{T}}, Any},
                      X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) =
    computestat!(normalized(t), out, work[3], unitnormalize!(work[1], X), unitnormalize!(work[2], Y))

#
# GroupMean
#

immutable GroupMean{R<:Statistic} <: Statistic
    transform::R
    indices::Vector{Int}
    offsets::Vector{Int}
    nchannels::Int
end
ngroups(x::GroupMean) = length(x.offsets)-1
outputdims{R<:GroupMean}(t::Type{R}) = 1

function GroupMean(t::Statistic, nchannels::Int, pairs::Vector{Vector{Tuple{Int,Int}}})
    indices = Int[]
    offsets = Int[1]
    for ipairgroup = 1:length(pairs)
        pairgroup = pairs[ipairgroup]
        for ipair in 1:size(pairgroup, 1)
            ch1, ch2 = pairgroup[ipair]
            (ch1 > nchannels || ch2 > nchannels || ch1 < 1 || ch2 < 1) && throw(ArgumentError("channel indices must be in range [1, nchannels"))
            if ch1 > ch2
                ch1, ch2 = ch2, ch1
            end
            push!(indices, (ch2-1)*nchannels+ch1)
        end
        push!(offsets, length(indices)+1)
    end
    GroupMean(t, indices, offsets, nchannels)
end

@generated function allocoutput{T,R<:PairwiseStatistic,N}(t::GroupMean{R}, X::AbstractArray{T,N}, Y::AbstractArray{T,N}=X)
    :(zeros(eltype(t.transform, X, Y), ngroups(t), $([:(size(X, $i)) for i = 3:N]...)))
end

allocwork{T}(t::GroupMean, X::AbstractVecOrMat{T}) =
    (allocoutput(t.transform, X), allocwork(t.transform, X))
allocwork{T}(t::GroupMean, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    (allocoutput(t.transform, X, Y), allocwork(t.transform, X, Y))

function grouppairs!{T}(t::GroupMean, out::AbstractVector{T}, chval::AbstractMatrix{T})
    length(out) == ngroups(t) || throw(ArgumentError("invalid output size"))
    @inbounds for ipairgroup = 1:ngroups(t)
        m = zero(T)
        @simd for ipair = t.offsets[ipairgroup]:t.offsets[ipairgroup+1]-1
            pairindex = t.indices[ipair]
            m += chval[pairindex]
        end
        out[ipairgroup] = m/(t.offsets[ipairgroup+1]-t.offsets[ipairgroup])
    end
    out
end

function computestat!{T}(t::GroupMean, out::AbstractArray, work::Tuple{Matrix, Any},
                         X::AbstractVecOrMat{T})
    chval, twork = work
    computestat!(t.transform, chval, twork, X)
    grouppairs!(t, out, chval)
    out
end
function computestat!{T}(t::GroupMean, out::AbstractArray, work::Tuple{Matrix, Any},
                         X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T})
    chval, twork = work
    computestat!(t.transform, chval, twork, X, Y)
    grouppairs!(t, out, chval)
    out
end

#
# (Inefficient) generic jackknife surrogate computation
#

abstract AbstractJackknifeSurrogates{R<:Statistic} <: Statistic

# Jackknifing single samples
immutable JackknifeSurrogates{R<:Statistic} <: AbstractJackknifeSurrogates{R}
    transform::R
end
Base.similar(t::JackknifeSurrogates, stat::Statistic) = JackknifeSurrogates(stat)

# Jackknifing multiple samples (e.g. trials)
# This is NOT a delete-n jackknife; it creates only div(ntrials(X), n)
# samples
immutable MultiJackknifeSurrogates{R<:Statistic} <: AbstractJackknifeSurrogates{R}
    transform::R
    n::Int
end
Base.similar(t::MultiJackknifeSurrogates, stat::Statistic) = MultiJackknifeSurrogates(stat, t.n)

immutable JackknifeSurrogatesOutput{T<:StridedArray,S<:StridedArray}
    trueval::T
    surrogates::S
end

# Number of elements to be deleted at a time
jnn(::JackknifeSurrogates) = 1
jnn(t::MultiJackknifeSurrogates) = t.n

function copyrows!(out, outrow, X, rowinds)
    for j = 1:size(X, 2)
        i = 1
        for x in rowinds
            @inbounds out[outrow+i-1, j] = X[x, j]
            i += 1
        end
    end
    out
end

@generated function allocoutput{T,N}(t::AbstractJackknifeSurrogates, X::AbstractArray{T,N}, Y::AbstractArray{T,N}=X)
    :(JackknifeSurrogatesOutput(zeros(eltype(t.transform, X, Y), size(X, 2), size(Y, 2), $([:(size(X, $i)) for i = 3:N]...)),
                                zeros(eltype(t.transform, X, Y), div(size(X, 1), jnn(t)), size(X, 2), size(Y, 2), $([:(size(X, $i)) for i = 3:N]...))))
end

function allocwork{T}(t::AbstractJackknifeSurrogates, X::AbstractVecOrMat{T})
    Xsurrogate = Array(T, size(X, 1)-jnn(t), size(X, 2))
    (Xsurrogate, allocoutput(t.transform, Xsurrogate), allocwork(t.transform, X), allocwork(t.transform, Xsurrogate))
end
function computestat!{T}(t::AbstractJackknifeSurrogates, out::JackknifeSurrogatesOutput, work,
                         X::AbstractVecOrMat{T})
    isa(work, Tuple{Matrix{T}, Any, Any, Any}) || throw(ArgumentError("invalid work object"))
    stat = t.transform
    ntrials, nch = size(X)
    njn = div(ntrials, jnn(t))
    Xsurrogate, outtmp, worktrue, worksurrogate = work
    trueval = out.trueval
    surrogates = out.surrogates

    (size(Xsurrogate, 1) == ntrials - jnn(t) && size(Xsurrogate, 2) == nch) || throw(ArgumentError("invalid work object"))
    ntrials % jnn(t) == 0 || throw(DimensionMismatch("ntrials not evenly divisible by $(jnn(t))"))
    chkinput(trueval, X)

    fill!(surrogates, NaN)

    # True value
    computestat!(stat, trueval, worktrue, X)

    # First jackknife
    copyrows!(Xsurrogate, 1, X, jnn(t)+1:size(X, 1))
    computestat!(stat, outtmp, worksurrogate, Xsurrogate)
    copy!(unsafe_view(surrogates, 1, :, :), outtmp)

    # Subsequent jackknifes
    for i = 2:njn
        copyrows!(Xsurrogate, (i-2)*jnn(t)+1, X, (i-2)*jnn(t)+1:(i-1)*jnn(t))
        computestat!(stat, outtmp, worksurrogate, Xsurrogate)
        copy!(unsafe_view(surrogates, i, :, :), outtmp)
    end

    out
end

function allocwork{T}(t::AbstractJackknifeSurrogates, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T})
    Xsurrogate = Array(T, size(X, 1)-jnn(t), size(X, 2))
    Ysurrogate = Array(T, size(Y, 1)-jnn(t), size(Y, 2))
    (Xsurrogate, Ysurrogate, allocoutput(t.transform, Xsurrogate, Ysurrogate),
     allocwork(t.transform, X, Y), allocwork(t.transform, Xsurrogate, Ysurrogate))
end
function computestat!{T}(t::AbstractJackknifeSurrogates, out::JackknifeSurrogatesOutput, work,
                         X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T})
    isa(work, Tuple{Matrix{T}, Matrix{T}, Any, Any, Any}) || throw(ArgumentError("invalid work object"))
    stat = t.transform
    ntrials, nchX = size(X)
    njn = div(ntrials, jnn(t))
    nchY = size(Y, 2)
    Xsurrogate, Ysurrogate, outtmp, worktrue, worksurrogate = work
    trueval = out.trueval
    surrogates = out.surrogates

    (size(Xsurrogate, 1) == size(Ysurrogate, 1) == ntrials - jnn(t) &&
     size(Xsurrogate, 2) == nchX &&
     size(Ysurrogate, 2) == nchY) || throw(ArgumentError("invalid work object"))
    ntrials % jnn(t) == 0 || throw(DimensionMismatch("ntrials not evenly divisible by $(jnn(t))"))
    ntrials > 0 || error("X is empty")
    chkinput(trueval, X, Y)

    fill!(surrogates, NaN)

    # True value
    computestat!(stat, trueval, worktrue, X, Y)

    # First jackknife
    copyrows!(Xsurrogate, 1, X, jnn(t)+1:size(X, 1))
    copyrows!(Ysurrogate, 1, Y, jnn(t)+1:size(X, 1))
    computestat!(stat, outtmp, worksurrogate, Xsurrogate, Ysurrogate)
    copy!(unsafe_view(surrogates, 1, :, :), outtmp)

    # Subsequent jackknifes
    for i = 2:njn
        copyrows!(Xsurrogate, (i-2)*jnn(t)+1, X, (i-2)*jnn(t)+1:(i-1)*jnn(t))
        copyrows!(Ysurrogate, (i-2)*jnn(t)+1, Y, (i-2)*jnn(t)+1:(i-1)*jnn(t))
        computestat!(stat, outtmp, worksurrogate, Xsurrogate, Ysurrogate)
        copy!(unsafe_view(surrogates, i, :, :), outtmp)
    end

    out
end

# Only unit normalize once
allocwork{R<:NormalizedPairwiseStatistic{false},T<:Real}(t::AbstractJackknifeSurrogates{R}, X::AbstractVecOrMat{Complex{T}}) =
    (Array(Complex{T}, size(X, 1), size(X, 2)), allocwork(similar(t, normalized(t.transform)), X))
allocwork{R<:NormalizedPairwiseStatistic{false},T<:Real}(t::AbstractJackknifeSurrogates{R}, X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) =
    (Array(Complex{T}, size(X, 1), size(X, 2)), Array(Complex{T}, size(Y, 1), size(Y, 2)), allocwork(similar(t, normalized(t.transform)), X, Y))
computestat!{R<:NormalizedPairwiseStatistic{false},T<:Real}(t::AbstractJackknifeSurrogates{R}, out::JackknifeSurrogatesOutput,
                                                            work, X::AbstractVecOrMat{Complex{T}}) =
    computestat!(similar(t, normalized(t.transform)), out, work[2], unitnormalize!(work[1], X))
computestat!{R<:NormalizedPairwiseStatistic{false},T<:Real}(t::AbstractJackknifeSurrogates{R}, out::JackknifeSurrogatesOutput,
                                                            work, X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) =
    computestat!(similar(t, normalized(t.transform)), out, work[3], unitnormalize!(work[1], X), unitnormalize!(work[2], Y))

# Estimate of variance from jackknife surrogates
function _jackknife_var!{T}(out::AbstractArray{T}, surrogates::AbstractArray{T})
    n = size(surrogates, 1)
    for j = 1:Base.trailingsize(surrogates, 2)
        # Compute mean
        m = zero(eltype(surrogates))
        @simd for i = 1:n
            @inbounds m += surrogates[i, j]
        end
        m /= n

        # Compute variance
        v = zero(eltype(out))
        @simd for i = 1:n
            @inbounds v += abs2(surrogates[i, j] - m)
        end

        @inbounds out[j] = v*(n-1)/n
    end
    out
end
function jackknife_var!{T}(out::AbstractArray{T}, jn::JackknifeSurrogatesOutput)
    size(out) == size(jn.trueval) || throw(ArgumentError("output size mismatch"))
    _jackknife_var!(out, jn.surrogates)
    out
end
jackknife_var(jn::JackknifeSurrogatesOutput) = jackknife_var!(similar(jn.trueval), jn)

# Estimate of (first-order) bias from jackknife surrogates
function jackknife_bias!(out::Matrix, jn::JackknifeSurrogatesOutput)
    trueval = jn.trueval
    surrogates = jn.surrogates
    n = size(surrogates, 1)

    for k = 1:size(surrogates, 3), j = 1:size(surrogates, 2)
        v = zero(eltype(out))
        @simd for i = 1:n
            @inbounds v += surrogates[i, j, k]
        end
        @inbounds out[j, k] = (n - 1)*(v/n - trueval[j, k])
    end
    out
end
jackknife_bias(jn::JackknifeSurrogatesOutput) = jackknife_bias!(similar(jn.trueval), jn)

# Jackknifing of n-d arrays
function computestat!(t::AbstractJackknifeSurrogates, out::JackknifeSurrogatesOutput, work, X::AbstractArray)
    !isempty(X) || error(ArgumentError("X is empty"))
    for i = 1:Base.trailingsize(X, 3)
        computestat!(t, JackknifeSurrogatesOutput(unsafe_view(out.trueval, :, :, i), unsafe_view(out.surrogates, :, :, :, i)),
                     work, unsafe_view(X, :, :, i))
    end
    out
end
function computestat!{T}(t::AbstractJackknifeSurrogates, out::JackknifeSurrogatesOutput, work, X::AbstractArray{T}, Y::AbstractArray{T})
    !isempty(X) || error(ArgumentError("X is empty"))
    for i = 1:Base.trailingsize(X, 3)
        computestat!(t, JackknifeSurrogatesOutput(unsafe_view(out.trueval, :, :, i), unsafe_view(out.surrogates, :, :, :, i)),
                     work, unsafe_view(X, :, :, i), unsafe_view(Y, :, :, i))
    end
    out
end
checknd(t::AbstractJackknifeSurrogates, out::JackknifeSurrogatesOutput, X::AbstractArray) = checknd(t.transform, out.trueval, X)
checknd(t::AbstractJackknifeSurrogates, out::JackknifeSurrogatesOutput, X::AbstractArray, Y::AbstractArray) = checknd(t.transform, out.trueval, X, Y)
function computestat_parallel_complete!(t::AbstractJackknifeSurrogates, out::JackknifeSurrogatesOutput, refs, ranges)
    for i = 1:length(refs)
        res = fetch(refs[i])
        isa(res, Exception) && rethrow(res)
        copy!(unsafe_view(out.trueval, :, :, ranges[i]), res.trueval)
        copy!(unsafe_view(out.surrogates, :, :, :, ranges[i]), res.surrogates)
    end
    out
end

#
# Computation of just the jackknife trueval and variance
#

immutable Jackknife{R<:Statistic,S} <: Statistic
    transform::S
end
Jackknife(t::Statistic) =
    Jackknife{typeof(t),JackknifeSurrogates{typeof(t)}}(JackknifeSurrogates(t))
Jackknife(t::Statistic, n::Int) =
    Jackknife{typeof(t),MultiJackknifeSurrogates{typeof(t)}}(MultiJackknifeSurrogates(t, n))

immutable JackknifeOutput{T<:StridedArray}
    trueval::T
    var::T
end

@generated function allocoutput{T,N}(t::Jackknife, X::AbstractArray{T,N}, Y::AbstractArray{T,N}=X)
    :(JackknifeOutput(zeros(eltype(t.transform.transform, X, Y), size(X, 2), size(Y, 2), $([:(size(X, $i)) for i = 3:N]...)),
                      zeros(eltype(t.transform.transform, X, Y), size(X, 2), size(Y, 2), $([:(size(X, $i)) for i = 3:N]...))))
end

allocwork{T}(t::Jackknife, X::AbstractVecOrMat{T}) =
    (allocwork(t.transform, X), zeros(eltype(t.transform.transform, X), div(size(X, 1), jnn(t.transform)), size(X, 2), size(X, 2)))
allocwork{T}(t::Jackknife, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    (allocwork(t.transform, X, Y), zeros(eltype(t.transform.transform, X, Y), div(size(X, 1), jnn(t.transform)), size(X, 2), size(Y, 2)))

function computestat!{T,S}(t::Jackknife, out::JackknifeOutput, work::Tuple{Any, Array{S,3}},
                                 X::AbstractVecOrMat{T})
    surrogates = JackknifeSurrogatesOutput(out.trueval, work[2])
    computestat!(t.transform, surrogates, work[1], X)
    jackknife_var!(out.var, surrogates)
    out
end
function computestat!{T,S}(t::Jackknife, out::JackknifeOutput, work::Tuple{Any, Array{S,3}},
                                 X::AbstractVecOrMat{T},
                                 Y::AbstractVecOrMat{T})
    surrogates = JackknifeSurrogatesOutput(out.trueval, work[2])
    computestat!(t.transform, surrogates, work[1], X, Y)
    jackknife_var!(out.var, surrogates)
    out
end

# Jackknifing of n-d arrays
function computestat!{S,T}(t::Jackknife, out::JackknifeOutput, work::S, X::AbstractArray{T})
    !isempty(X) || error(ArgumentError("X is empty"))
    for i = 1:Base.trailingsize(X, 3)
        surrogates = JackknifeSurrogatesOutput(unsafe_view(out.trueval, :, :, i), work[2])
        computestat!(t.transform, surrogates, work[1], unsafe_view(X, :, :, i))
        jackknife_var!(unsafe_view(out.var, :, :, i), surrogates)
    end
    out
end
function computestat!{S,T}(t::Jackknife, out::JackknifeOutput, work::S, X::AbstractArray{T}, Y::AbstractArray{T})
    !isempty(X) || error(ArgumentError("X is empty"))
    for i = 1:Base.trailingsize(X, 3)
        surrogates = JackknifeSurrogatesOutput(unsafe_view(out.trueval, :, :, i), work[2])
        computestat!(t.transform, surrogates, work[1], unsafe_view(X, :, :, i), unsafe_view(Y, :, :, i))
        jackknife_var!(unsafe_view(out.var, :, :, i), surrogates)
    end
    out
end
checknd(t::Jackknife, out::JackknifeOutput, X::AbstractArray) = checknd(t.transform.transform, out.trueval, X)
checknd(t::Jackknife, out::JackknifeOutput, X::AbstractArray, Y::AbstractArray) = checknd(t.transform.transform, out.trueval, X, Y)
function computestat_parallel_complete!(t::Jackknife, out::JackknifeOutput, refs, ranges)
    for i = 1:length(refs)
        res = fetch(refs[i])
        isa(res, Exception) && rethrow(res)
        copy!(unsafe_view(out.trueval, :, :, ranges[i]), res.trueval)
        copy!(unsafe_view(out.var, :, :, ranges[i]), res.trueval)
    end
    out
end

#
# Correlation of jackknife surrogates with another statistic
#

immutable JackknifeCorrelation{R<:Statistic,S,T<:AbstractFloat,F} <: PairwiseStatistic
    transform::S
    y::Vector{T}
    yssq::T
    f::F

    function JackknifeCorrelation(transform, y, f)
        μy = mean(y)
        yssq = zero(eltype(y))
        centeredy = zeros(eltype(y), length(y))
        @inbounds @simd for i = 1:length(y)
            centeredy[i] = y[i] - μy
            yssq += abs2(y[i] - μy)
        end
        new(transform, centeredy, yssq, f)
    end
end
JackknifeCorrelation(t::PairwiseStatistic, y::AbstractVector, f=Base.IdFun()) =
    JackknifeCorrelation{typeof(t),JackknifeSurrogates{typeof(t)},eltype(y),typeof(f)}(JackknifeSurrogates(t), float(y), f)
JackknifeCorrelation(t::PairwiseStatistic, n::Int, y::AbstractVector, f=Base.IdFun()) =
    JackknifeCorrelation{typeof(t),MultiJackknifeSurrogates{typeof(t)},eltype(y),typeof(f)}(MultiJackknifeSurrogates(t, n), float(y), f)

Base.eltype{T<:Real}(::JackknifeCorrelation, X::AbstractArray{Complex{T}}) = T
Base.eltype{T<:Real}(::JackknifeCorrelation, X::AbstractArray{T}) = T

function allocwork(t::JackknifeCorrelation, X::AbstractVecOrMat)
    outtype = eltype(t.transform.transform, X)
    outtype <: Real || throw(ArgumentError("JackknifeCorrelation defined only for statistics with real outputs"))
    (zeros(outtype, div(size(X, 1), jnn(t.transform)), size(X, 2), size(X, 2)), allocwork(t.transform, X))
end
function allocwork{T}(t::JackknifeCorrelation, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T})
    outtype = eltype(t.transform.transform, X)
    outtype <: Real || throw(ArgumentError("JackknifeCorrelation defined only for statistics with real outputs"))
    (zeros(outtype, div(size(X, 1), jnn(t.transform)), size(X, 2), size(Y, 2)), allocwork(t.transform, X, Y))
end

function jackknife_cor!{T<:Real}(t::JackknifeCorrelation, out::AbstractVecOrMat{T}, surrogates::AbstractArray{T,3})
    y = t.y
    @inbounds for itrail = 1:Base.trailingsize(surrogates, 2)
        μx = zero(T)
        @simd for isurrogate = 1:size(surrogates, 1)
            surrogates[isurrogate, itrail] = t.f(surrogates[isurrogate, itrail])
            μx += surrogates[isurrogate, itrail]
        end
        μx /= size(surrogates, 1)
        cov = zero(T)
        xssq = zero(T)
        @simd for isurrogate = 1:size(surrogates, 1)
            centeredx = surrogates[isurrogate, itrail] - μx
            cov += centeredx*y[isurrogate]
            xssq += abs2(centeredx)
        end
        out[itrail] = -cov/sqrt(xssq*t.yssq)
    end
    out
end

function computestat!{T<:Real}(t::JackknifeCorrelation, out::AbstractMatrix{T}, work::Tuple{Array{T,3},Any},
                               X::AbstractVecOrMat)
    length(t.y) == div(size(X, 1), jnn(t.transform)) || throw(DimensionMismatch("length of y does not match number of surrogates"))
    surrogates = JackknifeSurrogatesOutput(out, work[1])
    computestat!(t.transform, surrogates, work[2], X)
    jackknife_cor!(t, out, work[1])
end
function computestat!{T<:Real,S}(t::JackknifeCorrelation, out::AbstractMatrix{T}, work::Tuple{Array{T,3},Any},
                                 X::AbstractVecOrMat{S}, Y::AbstractVecOrMat{S})
    length(t.y) == div(size(X, 1), jnn(t.transform)) || throw(DimensionMismatch("length of y does not match number of surrogates"))
    surrogates = JackknifeSurrogatesOutput(out, work[1])
    computestat!(t.transform, surrogates, work[2], X, Y)
    jackknife_cor!(t, out, work[1])
end

#
# Jackknife{GroupMean}
#

@generated function allocoutput{T,R<:Statistic,N}(t::Jackknife{GroupMean{R}}, X::AbstractArray{T,N}, Y::AbstractArray{T,N}=X)
    :(JackknifeOutput(zeros(eltype(t.transform.transform.transform, X, Y), ngroups(t.transform.transform), $([:(size(X, $i)) for i = 3:N]...)),
                      zeros(eltype(t.transform.transform.transform, X, Y), ngroups(t.transform.transform), $([:(size(X, $i)) for i = 3:N]...))))
end

allocwork{R<:Statistic}(t::Jackknife{GroupMean{R}}, X::AbstractVecOrMat) =
    (allocwork(JackknifeSurrogates(t.transform.transform.transform), X),
     allocoutput(JackknifeSurrogates(t.transform.transform.transform), X),
     zeros(eltype(t.transform.transform.transform, X), size(X, 1), ngroups(t.transform.transform)))
allocwork{T,R<:Statistic}(t::Jackknife{GroupMean{R}}, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    (allocwork(JackknifeSurrogates(t.transform.transform.transform), X, Y),
     allocoutput(JackknifeSurrogates(t.transform.transform.transform), X, Y),
     zeros(eltype(t.transform.transform.transform, X), size(X, 1), ngroups(t.transform.transform)))

function grouppairs!{T,R<:Statistic}(t::Jackknife{GroupMean{R}}, group_trueval::AbstractArray{T},
                                     group_var::AbstractArray{T}, group_surrogates::AbstractMatrix{T},
                                     jn::JackknifeSurrogatesOutput)
    fill!(group_surrogates, zero(T))
    jnt = t.transform.transform

    # First compute means of groups and jackknife surrogates
    for ipairgroup = 1:ngroups(jnt)
        m = zero(T)
        for ipair = jnt.offsets[ipairgroup]:jnt.offsets[ipairgroup+1]-1
            pairindex = jnt.indices[ipair]
            m += jn.trueval[pairindex]
            @simd for isurrogate = 1:size(jn.surrogates, 1)
                group_surrogates[isurrogate, ipairgroup] += jn.surrogates[isurrogate, pairindex]
            end
        end
        group_trueval[ipairgroup] = m/(jnt.offsets[ipairgroup+1]-jnt.offsets[ipairgroup])
    end

    # Now compute jackknife variances
    _jackknife_var!(group_var, group_surrogates)

    # Adjust for taking the means
    for ipairgroup = 1:ngroups(jnt)
        group_var[ipairgroup] /= abs2(jnt.offsets[ipairgroup+1]-jnt.offsets[ipairgroup])
    end
    nothing
end

function computestat!{R<:Statistic,S}(t::Jackknife{GroupMean{R}}, out::JackknifeOutput,
                                      work::Tuple{Any, JackknifeSurrogatesOutput{Array{S,2}, Array{S,3}}, Array{S,2}},
                                      X::AbstractVecOrMat)
    computestat!(JackknifeSurrogates(t.transform.transform.transform), work[2], work[1], X)
    grouppairs!(t, out.trueval, out.var, work[3], work[2])
    out
end
function computestat!{R<:Statistic,T,S}(t::Jackknife{GroupMean{R}}, out::JackknifeOutput,
                                              work::Tuple{Any, JackknifeSurrogatesOutput{Array{S,2}, Array{S,3}}, Array{S,2}},
                                              X::AbstractVecOrMat{T},
                                              Y::AbstractVecOrMat{T})
    computestat!(JackknifeSurrogates(t.transform.transform.transform), work[2], work[1], X, Y)
    grouppairs!(t, out.trueval, out.var, work[3], work[2])
    out
end

# Group jackknifing of n-d arrays
function computestat!{R<:Statistic,S}(t::Jackknife{GroupMean{R}}, out::JackknifeOutput,
                                      work::Tuple{Any, JackknifeSurrogatesOutput{Array{S,2}, Array{S,3}}, Array{S,2}}, X::AbstractArray)
    !isempty(X) || throw(ArgumentError("X is empty"))
    for i = 1:Base.trailingsize(X, 3)
        computestat!(JackknifeSurrogates(t.transform.transform.transform), work[2], work[1], unsafe_view(X, :, :, i))
        grouppairs!(t, unsafe_view(out.trueval, :, i), unsafe_view(out.var, :, i), work[3], work[2])
    end
    out
end
function computestat!{R<:Statistic,T,S}(t::Jackknife{GroupMean{R}}, out::JackknifeOutput,
                                        work::Tuple{Any, JackknifeSurrogatesOutput{Array{S,2}, Array{S,3}}, Array{S,2}}, X::AbstractArray{T}, Y::AbstractArray{T})
    !isempty(X) || throw(ArgumentError("X is empty"))
    for i = 1:Base.trailingsize(X, 3)
        computestat!(JackknifeSurrogates(t.transform.transform.transform), work[2], work[1], unsafe_view(X, :, :, i), unsafe_view(Y, :, :, i))
        grouppairs!(t, unsafe_view(out.trueval, :, i), unsafe_view(out.var, :, i), work[3], work[2])
    end
    out
end

#
# Bootstrap and Permutation definitions and common functionality
#

immutable Bootstrap{R<:Statistic} <: Statistic
    transform::R
    indices::Matrix{Int32}    # Sample weights for each bootstrap. ntrials x nbootstraps
end
Bootstrap{T<:Integer}(t::Statistic, indices::Matrix{T}) = Bootstrap(t, convert(Matrix{Int32}, indices))
Bootstrap(transform::Statistic, ntrials::Int, nbootstraps::Int) =
    Bootstrap(transform, rand(Int32(1):Int32(ntrials), ntrials, nbootstraps))

immutable Permutation{S} <: Statistic
    transform::S
    indices::Matrix{Int32}
end
function Permutation(transform::Statistic, ntrials::Int, nperm::Int)
    indices = Array(Int32, ntrials, nperm)
    v = [1:ntrials;]
    for iperm = 1:nperm
        indices[:, iperm] = shuffle!(v)
    end
    Permutation(transform, indices)
end

outputdims{R<:Statistic}(t::Union{Type{Bootstrap{R}}, Type{Permutation{R}}}) = outputdims(R)+1
@generated function allocoutput{T,N}(t::Union{Bootstrap, Permutation}, X::AbstractArray{T,N}, Y::AbstractArray{T,N}=X)
    :(zeros(eltype(t.transform, X, Y), size(X, 2), size(Y, 2), size(t.indices, 2), $([:(size(X, $i)) for i = 3:N]...)))
end
@generated function allocoutput{T,R<:Statistic,N}(t::Union{Bootstrap{GroupMean{R}}, Permutation{GroupMean{R}}}, X::AbstractArray{T,N}, Y::AbstractArray{T,N}=X)
    :(zeros(eltype(t.transform.transform, X, Y), ngroups(t.transform), size(t.indices, 2), $([:(size(X, $i)) for i = 3:N]...)))
end

# Copy `indices` from X to out
function copyindices!(out, X, indices, iindices)
    @inbounds for j = 1:size(X, 2)
        @simd for i = 1:size(indices, 1)
            @inbounds out[i, j] = X[indices[i, iindices], j]
        end
    end
    out
end

#
# Bootstrapping
#

allocwork{T}(t::Bootstrap, X::AbstractVecOrMat{T}) =
    (Array(T, size(X, 1), size(X, 2)), allocwork(t.transform, X))
function bootstrap!{T,U}(t::Statistic, indices::Matrix{Int32}, out::AbstractArray,
                         Xbootstrap::Matrix{T}, work::U, X::AbstractVecOrMat{T})
    size(X, 1) == size(indices, 1) || throw(ArgumentError("number of bootstrap trials does not match number of trials"))
    (size(X, 1) == size(Xbootstrap, 1) &&
     size(X, 2) == size(Xbootstrap, 2)) || throw(ArgumentError("invalid work object"))

    for ibootstrap = 1:size(indices, 2)
        copyindices!(Xbootstrap, X, indices, ibootstrap)
        computestat!(t, output_view(t, out, ibootstrap), work, Xbootstrap)
    end
    out
end
computestat!{R<:Statistic,T}(t::Bootstrap{R}, out::AbstractArray,
                             work::Tuple{Matrix{T}, Any},
                             X::AbstractVecOrMat{T}) =
    bootstrap!(t.transform, t.indices, out, work[1], work[2], X)

allocwork{T}(t::Bootstrap, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    (Array(T, size(X, 1), size(X, 2)), Array(T, size(Y, 1), size(Y, 2)), allocwork(t.transform, X, Y))
function bootstrap!{T,U}(t::Statistic, indices::Matrix{Int32}, out::AbstractArray,
                         Xbootstrap::Matrix{T}, Ybootstrap::Matrix{T}, work::U,
                         X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T})
    size(X, 1) == size(indices, 1) || throw(ArgumentError("number of bootstrap trials does not match number of trials"))
    (size(X, 1) == size(Xbootstrap, 1) &&
     size(Y, 1) == size(Ybootstrap, 1) &&
     size(X, 2) == size(Xbootstrap, 2) &&
     size(Y, 2) == size(Ybootstrap, 2)) || throw(ArgumentError("invalid work object"))

    for ibootstrap = 1:size(indices, 2)
        copyindices!(Xbootstrap, X, indices, ibootstrap)
        copyindices!(Ybootstrap, Y, indices, ibootstrap)
        computestat!(t, output_view(t, out, ibootstrap), work, Xbootstrap, Ybootstrap)
    end
    out
end
computestat!{R<:Statistic,T}(t::Bootstrap{R}, out::AbstractArray,
                             work::Tuple{Matrix{T}, Matrix{T}, Any},
                             X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    bootstrap!(t.transform, t.indices, out, work[1], work[2], work[3], X, Y)

# Only unit normalize once
computestat!{R<:NormalizedPairwiseStatistic{false},T}(t::Bootstrap{R}, out::AbstractArray,
                                                      work::Tuple{Matrix{T}, Any},
                                                      X::AbstractVecOrMat{T}) =
    bootstrap!(normalized(t.transform), t.indices, out, work[1], work[2][2], unitnormalize!(work[2][1], X))
computestat!{R<:NormalizedPairwiseStatistic{false},T}(t::Bootstrap{R}, out::AbstractArray,
                                                      work::Tuple{Matrix{T}, Matrix{T}, Any},
                                                      X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    bootstrap!(normalized(t.transform), t.indices, out, work[1], work[2], work[3][3],
               unitnormalize!(work[3][1], X), unitnormalize!(work[3][2], Y))

#
# Permutation
#

allocwork{T}(t::Permutation, X::AbstractVecOrMat{T}) =
    (Array(T, size(X, 1), size(X, 2)), allocwork(t.transform, X, X))
# Separate allocwork method here to avoid allocating a second array for normalizing X
allocwork{T,R<:NormalizedPairwiseStatistic{false}}(t::Permutation{R}, X::AbstractVecOrMat{T}) =
    (Array(T, size(X, 1), size(X, 2)), Array(T, size(X, 1), size(X, 2)), allocwork(normalized(t.transform), X, X))
function permutation!{T,V}(t::Statistic, indices::Matrix{Int32}, out::AbstractArray, perm::Array{T,2}, work::V,
                           X::AbstractVecOrMat{T})
    size(X, 1) == size(indices, 1) || throw(ArgumentError("number of permutation trials does not match number of trials"))
    (size(X, 1) == size(perm, 1) &&
     size(X, 2) == size(perm, 2)) || throw(ArgumentError("invalid work object"))

    for iperm = 1:size(indices, 2)
        copyindices!(perm, X, indices, iperm)
        computestat!(t, output_view(t, out, iperm), work, perm, X)
    end
    out
end
computestat!{R<:Statistic,T}(t::Permutation{R}, out::AbstractArray, work::Tuple{Matrix{T}, Any},
                             X::AbstractVecOrMat{T}) =
    permutation!(t.transform, t.indices, out, work[1], work[2], X)

allocwork{T}(t::Permutation, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    (Array(T, size(X, 1), size(X, 2)), allocwork(t.transform, X, Y))
function permutation!{T,V}(t::Statistic, indices::Matrix{Int32}, out::AbstractArray, perm::Array{T,2}, work::V,
                           X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T})
    size(X, 1) == size(indices, 1) || throw(ArgumentError("number of permutation trials does not match number of trials"))
    (size(X, 1) == size(perm, 1) &&
     size(X, 2) == size(perm, 2)) || throw(ArgumentError("invalid work object"))

    for iperm = 1:size(indices, 2)
        copyindices!(perm, X, indices, iperm)
        computestat!(t, output_view(t, out, iperm), work, perm, Y)
    end
    out
end
computestat!{R<:Statistic,T}(t::Permutation{R}, out::AbstractArray, work::Tuple{Matrix{T}, Any},
                             X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    permutation!(t.transform, t.indices, out, work[1], work[2], X, Y)

# Only unit normalize once
computestat!{R<:NormalizedPairwiseStatistic{false},T}(t::Permutation{R}, out::AbstractArray,
                                                      work::Tuple{Matrix{T}, Matrix{T}, Any},
                                                      X::AbstractVecOrMat{T}) =
    permutation!(normalized(t.transform), t.indices, out, work[2], work[3], unitnormalize!(work[1], X))
computestat!{R<:NormalizedPairwiseStatistic{false},T}(t::Permutation{R}, out::AbstractArray,
                                                      work::Tuple{Matrix{T}, Any},
                                                      X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    permutation!(normalized(t.transform), t.indices, out, work[1], work[2][3],
               unitnormalize!(work[2][1], X), unitnormalize!(work[2][2], Y))

#
# Efficient permutation for JackknifeCorrelation
#

function perm_jackknife_cor!{T<:Real}(t::JackknifeCorrelation, out::AbstractArray{T}, surrogates::AbstractArray{T,3}, permy::Matrix)
    y = t.y
    trailsize = Base.trailingsize(surrogates, 2)
    @inbounds for itrail = 1:trailsize
        μx = zero(T)
        @simd for isurrogate = 1:size(surrogates, 1)
            surrogates[isurrogate, itrail] = t.f(surrogates[isurrogate, itrail])
            μx += surrogates[isurrogate, itrail]
        end
        μx /= size(surrogates, 1)

        xssq = zero(T)
        @simd for isurrogate = 1:size(surrogates, 1)
            surrogates[isurrogate, itrail] -= μx
            xssq += abs2(surrogates[isurrogate, itrail])
        end

        for iperm = 1:size(permy, 2)
            cov = zero(T)
            @simd for isurrogate = 1:size(surrogates, 1)
                cov += surrogates[isurrogate, itrail]*permy[isurrogate, iperm]
            end
            out[(iperm-1)*trailsize+itrail] = -cov/sqrt(xssq*t.yssq)
        end
    end
    out
end

function invpermy(y, indices)
    size(indices, 1) == length(y) || throw(DimensionMismatch("length of indices does not match length of y"))
    permy = zeros(eltype(y), size(indices))
    for j = 1:size(indices, 2)
        @simd for i = 1:size(indices, 1)
            @inbounds permy[indices[i, j], j] = y[i]
        end
    end
    permy
end

immutable JackknifeCorrelationPermutationWork{T<:Real,S}
    yvals::Matrix{T}
    twork::S
end

allocwork{R<:JackknifeCorrelation}(t::Permutation{R}, X::AbstractVecOrMat) =
    JackknifeCorrelationPermutationWork(invpermy(t.transform.y, t.indices), allocwork(t.transform, X))
function computestat!{R<:JackknifeCorrelation}(t::Permutation{R}, out::AbstractArray, work::JackknifeCorrelationPermutationWork,
                                               X::AbstractVecOrMat)
    indices = t.indices
    yvals = work.yvals
    twork = work.twork
    njn = div(size(X, 1), jnn(t.transform.transform))
    njn == size(indices, 1) || throw(ArgumentError("number of permutation trials does not match number of trials"))
    njn == length(t.transform.y) || throw(DimensionMismatch("length of y does not match number of surrogates"))

    surrogates = JackknifeSurrogatesOutput(output_view(t.transform.transform.transform, out, 1), twork[1])
    computestat!(t.transform.transform, surrogates, twork[2], X)
    perm_jackknife_cor!(t.transform, out, twork[1], yvals)
end

allocwork{R<:JackknifeCorrelation,T}(t::Permutation{R}, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    JackknifeCorrelationPermutationWork(invpermy(t.transform.y, t.indices), allocwork(t.transform, X, Y))
function computestat!{R<:JackknifeCorrelation,S}(t::Permutation{R}, out::AbstractArray, work::JackknifeCorrelationPermutationWork,
                                                 X::AbstractVecOrMat{S}, Y::AbstractVecOrMat{S})
    indices = t.indices
    yvals = work.yvals
    twork = work.twork
    njn = div(size(X, 1), jnn(t.transform.transform))
    njn == size(indices, 1) || throw(ArgumentError("number of permutation trials does not match number of trials"))
    njn == length(t.transform.y) || throw(DimensionMismatch("length of y does not match number of surrogates"))

    surrogates = JackknifeSurrogatesOutput(output_view(t.transform.transform.transform, out, 1), twork[1])
    computestat!(t.transform.transform, surrogates, twork[2], X, Y)
    perm_jackknife_cor!(t.transform, out, twork[1], yvals)
end

#
# Transforms
#

include("transforms/spectrum.jl")
include("transforms/correlation.jl")
include("transforms/coherence.jl")
include("transforms/plv.jl")
include("transforms/pli.jl")
include("transforms/jammalamadaka.jl")
include("transforms/juppmardia.jl")
include("transforms/modulationindex.jl")
include("transforms/uniformscores.jl")


text:
struct Node{R <: AbstractFloat, N}
    ind::CartesianIndex{N}
    val::R
end

<(a::Node, b::Node) = a.val < b.val

ext: jl


text:
# benchmarking

using ArrayViews

# auxiliary functions

mps(a::AbstractArray, rtimes::Int, et::Float64) = length(a) * 1.0e-6 * rtimes / et

myrepr(r::Range) = repr(r)
myrepr(r::Real) = repr(r)
myrepr(r::Colon) = ":"

# timing functions

function time_view1d{T}(a::AbstractArray{T}, rtimes::Int)
    # warm up
    n::Int = length(a)
    for i = 1:n
        @inbounds a[i]
    end

    # benchmark
    s = zero(T)
    et = @elapsed for t = 1 : rtimes
        for i = 1:n
            @inbounds s += a[i]
        end
    end
    return (et, s)
end

function time_view2d{T}(a::AbstractArray{T}, rtimes::Int)
    # warm up
    m::Int = size(a, 1)
    n::Int = size(a, 2)
    for j = 1:n, i = 1:m
        @inbounds a[i,j]
    end

    # benchmark
    s = zero(T)
    et = @elapsed for t = 1 : rtimes
        for j = 1:n, i = 1:m
            @inbounds s += a[i,j]
        end
    end
    return (et, s)
end

function time_view3d{T}(a::AbstractArray{T}, rtimes::Int)
    # warm up
    m::Int = size(a, 1)
    n::Int = size(a, 2)
    k::Int = size(a, 3)
    for l = 1:k, j = 1:n, i = 1:m
        @inbounds a[i,j,l]
    end

    # benchmark
    s = zero(T)
    et = @elapsed for t = 1 : rtimes
        for l = 1:k, j = 1:n, i = 1:m
            @inbounds s += a[i,j,l]
        end
    end
    return (et, s)
end

function perf_view(a::Array, i1; rtimes::Int=100000)
    et_s, _ = time_view1d(sub(a, i1), rtimes)
    et_v, _ = time_view1d(view(a, i1), rtimes)
    et_u, _ = time_view1d(unsafe_view(a, i1), rtimes)

    v = view(a, i1)
    mps_s = mps(v, rtimes, et_s)
    mps_v = mps(v, rtimes, et_v)
    mps_u = mps(v, rtimes, et_u)

    @printf("%-24s:   %10.3f MPS   %10.3f MPS (%6.3fx)   %10.3f MPS (%6.3fx)\n",
        "[$(myrepr(i1))]", mps_s, mps_v, mps_v / mps_s, mps_u, mps_u / mps_s)
end

function perf_view(a::Array, i1, i2; rtimes::Int=100000)
    et_s, _ = time_view2d(sub(a, i1, i2), rtimes)
    et_v, _ = time_view2d(view(a, i1, i2), rtimes)
    et_u, _ = time_view2d(unsafe_view(a, i1, i2), rtimes)

    v = view(a, i1, i2)
    mps_s = mps(v, rtimes, et_s)
    mps_v = mps(v, rtimes, et_v)
    mps_u = mps(v, rtimes, et_u)

    @printf("%-24s:   %10.3f MPS   %10.3f MPS (%6.3fx)   %10.3f MPS (%6.3fx)\n",
        "[$(myrepr(i1)), $(myrepr(i2))]", mps_s, mps_v, mps_v / mps_s, mps_u, mps_u / mps_s)
end

function perf_view(a::Array, i1, i2, i3; rtimes::Int=100000)
    et_s, _ = time_view3d(sub(a, i1, i2, i3), rtimes)
    et_v, _ = time_view3d(view(a, i1, i2, i3), rtimes)
    et_u, _ = time_view3d(unsafe_view(a, i1, i2, i3), rtimes)

    v = view(a, i1, i2, i3)
    mps_s = mps(v, rtimes, et_s)
    mps_v = mps(v, rtimes, et_v)
    mps_u = mps(v, rtimes, et_u)

    @printf("%-24s:   %10.3f MPS   %10.3f MPS (%6.3fx)   %10.3f MPS (%6.3fx)\n",
        "[$(myrepr(i1)), $(myrepr(i2)), $(myrepr(i3))]", mps_s, mps_v, mps_v / mps_s, mps_u, mps_u / mps_s)
end

# benchmarks

println("Indexing                      sub              view                       unsafe_view")
println("--------------------------------------------------------------------------------------------------")

const a1 = rand(1024)

gc_disable()

println("1D views")
perf_view(a1, :)
perf_view(a1, 1:1024)
perf_view(a1, 1:2:1024)
println()

gc_enable()

const a2 = rand(32, 32)

gc_disable()

println("2D views")
perf_view(a2, :, :)
perf_view(a2, :, 1:32)
perf_view(a2, :, 1:2:32)

perf_view(a2, 1:30, :)
perf_view(a2, 1:30, 1:32)
perf_view(a2, 1:30, 1:2:32)

perf_view(a2, 1:2:30, :)
perf_view(a2, 1:2:30, 1:32)
perf_view(a2, 1:2:30, 1:2:32)
println()

gc_enable()

const a3 = rand(16, 8, 8)

gc_disable()

println("3D views")
perf_view(a3, :, :, :)
perf_view(a3, :, :, 1:8)
perf_view(a3, :, :, 1:2:8)

perf_view(a3, :, 1:8, :)
perf_view(a3, :, 1:8, 1:8)
perf_view(a3, :, 1:8, 1:2:8)

perf_view(a3, :, 1:2:8, :)
perf_view(a3, :, 1:2:8, 1:8)
perf_view(a3, :, 1:2:8, 1:2:8)

perf_view(a3, 1:16, :, :)
perf_view(a3, 1:16, :, 1:8)
perf_view(a3, 1:16, :, 1:2:8)

perf_view(a3, 1:16, 1:8, :)
perf_view(a3, 1:16, 1:8, 1:8)
perf_view(a3, 1:16, 1:8, 1:2:8)

perf_view(a3, 1:16, 1:2:8, :)
perf_view(a3, 1:16, 1:2:8, 1:8)
perf_view(a3, 1:16, 1:2:8, 1:2:8)

perf_view(a3, 1:2:16, :, :)
perf_view(a3, 1:2:16, :, 1:8)
perf_view(a3, 1:2:16, :, 1:2:8)

perf_view(a3, 1:2:16, 1:8, :)
perf_view(a3, 1:2:16, 1:8, 1:8)
perf_view(a3, 1:2:16, 1:8, 1:2:8)

perf_view(a3, 1:2:16, 1:2:8, :)
perf_view(a3, 1:2:16, 1:2:8, 1:8)
perf_view(a3, 1:2:16, 1:2:8, 1:2:8)
println()

gc_enable()


text:
"""
Used in gp.jl. Data pre- / post-processing. Takes a ProfileData object and prepares it for use in GP.
"""

abstract type Problem end

abstract type SequentialProblem <: Problem end
abstract type ResidualProblem <: Problem end
abstract type SlackProblem <: Problem end

default_modify_predictor_fn(x, 𝒟, time_index) = x

struct Sequential <: SequentialProblem
    type::String # "T" or "dT" or "wT" or "KPP" or "TKE"
    parameters # Parameters(...) (see OceanTurb documentation)
    modify_predictor_fn
    convective_adjust
end
Sequential(type; parameters=nothing, modify_predictor_fn=default_modify_predictor_fn, convective_adjust=false) = Sequential(type, parameters, modify_predictor_fn, convective_adjust)

struct Residual <: ResidualProblem
    type::String # "KPP" or "TKE"
    parameters # Parameters(...) (see OceanTurb documentation)
    modify_predictor_fn
    convective_adjust
end
Residual(type; parameters=nothing, modify_predictor_fn=default_modify_predictor_fn, convective_adjust=false) = Residual(type, parameters, modify_predictor_fn, convective_adjust)

struct Slack <: SlackProblem
    type::String # "KPP" or "TKE"
    parameters # Parameters(...) (see OceanTurb documentation)
    modify_predictor_fn
    convective_adjust
end
Slack(type; parameters=nothing, modify_predictor_fn=default_modify_predictor_fn, convective_adjust=false) = Slack(type, parameters, modify_predictor_fn, convective_adjust)

mutable struct Sequential_dT <: SequentialProblem
    variable::String #"T" or "wT"
    Δt::Number # assumes constant time interval between all timesteps
    scaling::Scaling
    Sequential_dT(variable, Δt) = new(variable, Δt) # incomplete initialization
end

mutable struct Sequential_T <: SequentialProblem
    variable::String # "T" or "wT"
    scaling::Scaling
    Sequential_T(variable) = new(variable) # incomplete initialization
end

mutable struct Sequential_wT <: SequentialProblem
    variable::String # "T" or "wT"
    scaling::Scaling
    Sequential_wT(variable) = new(variable) # incomplete initialization
end

mutable struct Sequential_KPP <: SequentialProblem
    variable::String
    physics_data
    evolve_physics_model_fn # Function
    Δt # assumes constant time interval between all timesteps
    scaling::Scaling
    Sequential_KPP(variable, physics_data, evolve_physics_model_fn, Δt) = new(variable, physics_data, evolve_physics_model_fn, Δt) # incomplete initialization
end

mutable struct Sequential_TKE <: SequentialProblem
    variable::String
    physics_data
    evolve_physics_model_fn # Function
    Δt # assumes constant time interval between all timesteps
    scaling::Scaling
    Sequential_TKE(variable, physics_data, evolve_physics_model_fn, Δt) = new(variable, physics_data, evolve_physics_model_fn, Δt) # incomplete initialization
end

mutable struct Residual_KPP <: ResidualProblem
    variable::String
    physics_data
    evolve_physics_model_fn # Function
    Δt # assumes constant time interval between all timesteps
    scaling::Scaling
    Residual_KPP(variable, physics_data, evolve_physics_model_fn, Δt) = new(variable, physics_data, evolve_physics_model_fn, Δt) # incomplete initialization
end

mutable struct Residual_TKE <: ResidualProblem
    variable::String
    physics_data
    evolve_physics_model_fn # Function
    Δt # assumes constant time interval between all timesteps
    scaling::Scaling
    Residual_TKE(variable, physics_data, evolve_physics_model_fn, Δt) = new(variable, physics_data, evolve_physics_model_fn, Δt) # incomplete initialization
end

mutable struct Slack_KPP <: SlackProblem
    variable::String
    physics_data::Array
    scaling::Scaling
    Slack_KPP(variable, physics_data) = new(variable, physics_data) # incomplete initialization
end

mutable struct Slack_TKE <: SlackProblem
    variable::String
    physics_data::Array
    scaling::Scaling
    Slack_TKE(variable, physics_data) = new(variable, physics_data) # incomplete initialization
end

"""
get_v(problem)
# Description
    Return data array from the OceananigansData object associated with the variable of interest.
# Arguments
- 'problem': (Problem).     what mapping you wish to evaluate with the model. (Sequential("T"), Sequential("wT"))
"""
function get_v(problem, les)
    if problem.type == "wT"
        return les.wT
    elseif problem.type in ["T", "dT", "KPP", "TKE"]
        return les.T
    else
        throw(error())
    end
end

"""
get_problem(problem::Problem, data::OceananigansData, timeseries)
# Description
    Creates an instance of a Problem struct depending on the type of mapping.
    Returns the problem and the correct data array from the OceananigansData objectssociated with the variable of interest.
# Arguments
- 'problem': (Problem).      what mapping you wish to evaluate with the model. (Sequential("T"), Sequential("wT"))
- 'les': (OceananigansData). LES simulation output
- 'v': (Array).              Array of values
- 'N²': (Number)             initial buoyancy stratification
- 'D': (Number)              number of grid points in the vertical
- 'Nt': (Integer)            number of time steps
"""
function get_problem(problem, les, v, N², D, Nt)

    Δt = les.t[2]-les.t[1]

    if typeof(problem) <: SequentialProblem

        if problem.type == "dT"
            Δt /= N² # Δt' = Δt / N²
            return Sequential_dT("T", Δt) # Sequential_dT(variable, Δt)

        elseif problem.type == "T"
            return Sequential_T("T")

        elseif problem.type == "wT"
            return Sequential_wT("wT")

        elseif problem.type == "KPP"
            f = closure_free_convection_kpp(problem.parameters, D, Δt, les, n_steps=1)
            evolve_physics_model_fn(T⁰) = f(; T⁰=T⁰)[:,2]
            kpp_data = Array{Array{Float64,1},1}(UndefInitializer(), Nt)
            kpp_data[1] = custom_avg(v[:,1], D)
            for i in 1:Nt-1
                kpp_data[i+1] = evolve_physics_model_fn(v[:,i])
            end
            return Sequential_KPP("T", kpp_data, evolve_physics_model_fn, Δt)

        elseif problem.type == "TKE"
            # Use the LES profile at time index i to predict time index i+1 using TKE
            f = closure_free_convection_tke(problem.parameters, D, Δt, les, n_steps=1)
            evolve_physics_model_fn2(T⁰) = f(; T⁰=T⁰)[:,2]
            tke_data = Array{Array{Float64,1},1}(UndefInitializer(), Nt)
            tke_data[1] = custom_avg(v[:,1], D)
            for i in 1:Nt-1
                tke_data[i+1] = evolve_physics_model_fn2(v[:,i])
            end
            return Sequential_TKE("T", tke_data, evolve_physics_model_fn2, Δt)

        else; throw(error())
        end

    elseif typeof(problem) <: ResidualProblem

        if problem.type == "KPP"
            # Use the LES profile at time index i to predict time index i+1 using KPP
            f = closure_free_convection_kpp(problem.parameters, D, Δt, les, n_steps=1)
            evolve_physics_model_fn3(T⁰) = f(; T⁰=T⁰)[:,2]
            kpp_data = Array{Array{Float64,1},1}(UndefInitializer(), Nt)
            kpp_data[1] = custom_avg(v[:,1], D)
            for i in 1:Nt-1
                kpp_data[i+1] = evolve_physics_model_fn3(v[:,i])
            end
            return Residual_KPP("T", kpp_data, evolve_physics_model_fn3, Δt)

        elseif problem.type == "TKE"
            # Use the LES profile at time index i to predict time index i+1 using TKE
            f = closure_free_convection_tke(problem.parameters, D, Δt, les, n_steps=1)
            evolve_physics_model_fn4(T⁰) = f(; T⁰=T⁰)[:,2]
            tke_data = Array{Array{Float64,1},1}(UndefInitializer(), Nt)
            tke_data[1] = custom_avg(v[:,1], D)
            for i in 1:Nt-1
                tke_data[i+1] = evolve_physics_model_fn4(v[:,i])
            end
            return Residual_TKE("T", tke_data, evolve_physics_model_fn4, Δt)

        else; throw(error())
        end

    elseif typeof(problem) <: SlackProblem

        if problem.type == "KPP"
            # Predict the full evolution of the temperature profile from the initial time step using KPP
            f = closure_free_convection_kpp_full_evolution(problem.parameters, D, Δt, les)
            kpp_data = f()
            kpp_data = [kpp_data[:,i] for i in 1:Nt]
            return Slack_KPP("T", kpp_data)

        elseif problem.type == "TKE"
            # Predict the full evolution of the temperature profile from the initial time step using TKE
            f = closure_free_convection_tke_full_evolution(problem.parameters, D, Δt, les)
            tke_data = f()
            tke_data = [tke_data[:,i] for i in 1:Nt]
            return Slack_TKE("T", tke_data)

        else; throw(error())
        end

    else; throw(error())
    end

end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
module MeasureSystems

#   This file is part of MeasureSystems.jl. It is licensed under the MIT license
#   MeasureSystems Copyright (C) 2021 Michael Reed

import Base: @pure
import UnitSystems
import UnitSystems: UnitSystem, Systems, Constants, Physics, Convert
import UnitSystems: Coupling, measure, unit, universe
export UnitSystems, Measure, measure, cache

# measure

using Measurements
struct Measure{N} end
const measure_cache = Measurement{Float64}[]
@pure measure(::Measure{N}) where N = measure_cache[N]
function cache(M::Measurement{Float64})
    N = findfirst(x->x==M,measure_cache)
    if isnothing(N)
        push!(measure_cache,M)
        N = Base.length(measure_cache)
    end
    return Measure{N}()
end
Base.show(io::IO,M::Measure{N}) where N = show(io,measure(M))

# unit systems

@pure mass(U::UnitSystem,S::UnitSystem) = electronmass(U,S)
@pure electronmass(𝘩,R∞) = αinv^2*R∞*2𝘩/𝘤
@pure electronmass(𝘩,R∞,C::Coupling) = inv(finestructure(C))^2*R∞*2𝘩/𝘤
@pure planckmass(U::UnitSystem,C::Coupling=universe(U)) = electronmass(U,C)/√coupling(C)
@pure planck(U::UnitSystem,C::Coupling=universe(U)) = 2π*planckreduced(U,C)
@pure newton(U::UnitSystem,C::Coupling=universe(U)) = lightspeed(U,C)*planckreduced(U,C)/planckmass(U,C)^2
@pure charge(U::UnitSystem,C::Coupling=universe(U)) = sqrt(2planck(U)/(permeability(U)/finestructure(U))/(lightspeed(U)*rationalization(U)*lorentz(U)^2))

for unit ∈ (:coupling,:finestructure,:electronunit,:protonunit,:protonelectron)
    @eval @pure $unit(C::Coupling) = UnitSystems.$unit(C)
    @eval @pure $unit(U::UnitSystem) = UnitSystems.$unit(universe(U))
end
for unit ∈ (:boltzmann,:planckreduced,:lightspeed,:permeability,:electronmass,:rationalization,:lorentz,:molarmass)
    @eval @pure $unit(U::UnitSystem) = UnitSystems.$unit(U)
end
for unit ∈ (:boltzmann,:planckreduced,:lightspeed,:permeability,:electronmass,:molarmass)
    @eval @pure $unit(U::UnitSystem,C::Coupling) = $unit(U)
end
for unit ∈ Constants
    @eval @pure $unit(U::UnitSystem,S::UnitSystem) = unit($unit(S)/$unit(U))
end
for unit ∈ Convert
    @eval begin
        @pure @inline $unit(v::Real,U::UnitSystem) = $unit(v,U,Metric)
        @pure @inline $unit(v::Real,U::UnitSystem,S::UnitSystem) = (u=$unit(U,S);isone(u) ? v : v/u)
        @pure @inline $unit(v::Real,U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ},S::UnitSystem{kB,ħ,𝘤,μ₀,mₑ}) where {kB,ħ,𝘤,μ₀,mₑ} = v
        @pure @inline $unit(U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ},S::UnitSystem{kB,ħ,𝘤,μ₀,mₑ}) where {kB,ħ,𝘤,μ₀,mₑ} = 1
    end
    if unit ∉ (Constants...,:permittivity,:charge,:magneticflux,:impedance,:conductance)
        @eval @pure @inline $unit(U::UnitSystem) = $unit(U,Metric)
    end
end
for unit ∈ (Systems...,Constants...,Physics...,Convert...)
    unit ∉ (:length,:time) && @eval export $unit
end

# common conversion factors

import UnitSystems: atm,g₀,lbm,slug,ft,ftUS,rankine,kelvin
for CAL ∈ (:cal,:calₜₕ,:cal₄,:cal₁₀,:cal₂₀,:calₘ,:calᵢₜ)
    KCAL = Symbol(:k,CAL)
    @eval import UnitSystems: $CAL, $KCAL
end

# fundamental constants, αinv = (34259-1/4366.8123)/250 # 137.036 exactly?

const mP = measurement("0.00000002176434(24)")
const αinv = measurement("137.035999084(21)")
const μₑᵤ = 1/measurement("1822.888486209(53)")
const μₚᵤ = measurement("1.007276466621(53)")
const R∞ = measurement("10973731.5681601(210)")
const KJ2014 = measurement("483597.8525(30)")*1e9
const RK2014 = measurement("25812.8074555(59)")

import UnitSystems: ΔνCs,Kcd,NA,kB,𝘩,𝘤,𝘦,ħ,Rᵤ,αL,RK1990,KJ1990,ħ1990
const mₑ,μ₀ = electronmass(𝘩,R∞),2𝘩/𝘤/αinv/𝘦^2 # ≈ 4π*(1e-7+5.5e-17), exact charge
const δμ₀,μₚₑ,αG,ħ2014 = μ₀-4π*1e-7,μₚᵤ/μₑᵤ,(mₑ/mP)^2,2/RK2014/KJ2014^2/π
const mₑ1990,mₑ2014 = electronmass(2π*ħ1990,R∞),electronmass(2π*ħ2014,R∞)

# engineering units # Thomson: αL = 1/2

const Universe = Coupling{cache(αG),cache(1/αinv),cache(μₑᵤ),cache(μₚᵤ)}()
const Gauss = UnitSystem{cache(1e10*Rᵤ*mₑ/μₑᵤ),1e7*ħ,100𝘤,1,cache(1000mₑ),4π,0.01/𝘤}()
const LorentzHeaviside = UnitSystem{cache(1e10*Rᵤ*mₑ/μₑᵤ),1e7*ħ,100𝘤,1,cache(1000mₑ),1,0.01/𝘤}()
const Thomson = UnitSystem{cache(1e10*Rᵤ*mₑ/μₑᵤ),1e7*ħ,100𝘤,1,cache(1000mₑ),4π,1/2}()
const Kennelly = UnitSystem{cache(Rᵤ*mₑ/μₑᵤ/0.001),ħ,𝘤,1e-7,cache(mₑ),4π}() # ?
const ESU = UnitSystem{cache(1e10*Rᵤ*mₑ/μₑᵤ),1e7*ħ,100𝘤,(100𝘤)^-2,cache(1000mₑ),4π}()
const ESU2019 = UnitSystem{1e7*kB,1e7*ħ,100𝘤,cache(1e3*μ₀/𝘤^2),cache(1000mₑ)}()
const EMU = UnitSystem{cache(1e10*Rᵤ*mₑ/μₑᵤ),1e7*ħ,100𝘤,1,cache(1000mₑ),4π}()
const EMU2019 = UnitSystem{1e7*kB,1e7*ħ,100𝘤,cache(1e7*μ₀),cache(1000mₑ)}()
const MTS = UnitSystem{cache(1e6*Rᵤ*mₑ/μₑᵤ),1000ħ,𝘤,4π/1e4,cache(mₑ/1000)}()
const Mixed = UnitSystem{cache(Rᵤ*mₑ/μₑᵤ/0.001),ħ,𝘤,cache(μ₀),cache(mₑ)}()
const Metric = UnitSystem{cache(Rᵤ*mₑ/μₑᵤ/0.001),ħ,𝘤,4π*1e-7,cache(mₑ)}()
const SI1976 = UnitSystem{cache(8.31432mₑ/μₑᵤ/0.001),ħ,𝘤,4π*1e-7,cache(mₑ)}()
const SI2019 = UnitSystem{kB,ħ,𝘤,cache(μ₀),cache(mₑ)}()
const CODATA = UnitSystem{cache(Rᵤ*mₑ2014/μₑᵤ/0.001),cache(ħ2014),𝘤,cache(2RK2014/𝘤/αinv),cache(mₑ2014)}()
const Conventional = UnitSystem{cache(Rᵤ*mₑ1990/μₑᵤ/0.001),ħ1990,𝘤,cache(2RK1990/𝘤/αinv),cache(mₑ1990)}()
const English = UnitSystem{kB*rankine/slug/ft^2,ħ/slug/ft^2,𝘤/ft,4π,cache(mₑ/slug)}()
const EnglishUS = UnitSystem{cache(1000Rᵤ*mₑ/μₑᵤ*rankine/slug/ftUS^2),ħ/slug/ftUS^2,𝘤/ftUS,4π,cache(mₑ/slug)}()

# astronomical units

const GMsun = measurement("1.32712440018(9)")*1e20
const GMearth = measurement("3.986004418(8)")*1e14
const GMjupiter = measurement("1.26686534(9)")*1e17

import UnitSystems: au,LD,day,pc,ly
const GG = 𝘤*ħ/mP^2; const mₛ = GMsun/GG; const Jₛ = mₛ*au^2/day^2; export mₛ,Jₛ,au,day
const IAU = UnitSystem{cache(Rᵤ*mₑ/μₑᵤ/0.001/Jₛ),cache(ħ/day/Jₛ),day*𝘤/au,cache(4π*1e-7*day^2/Jₛ),cache(mₑ/mₛ)}()

# aliased & humorous units

const mf = mass(90/lbm,Metric,English); const Jf = mf*(201.168/14day)^2
const FFF = UnitSystem{cache(1000Rᵤ*mₑ/μₑᵤ*rankine/Jf),cache(ħ/14day/Jf),14day*𝘤/201.168,0,cache(mₑ/mf)}()
const units, US, SI, MKS, temp = UnitSystem, UnitSystem, SI2019, Metric, temperature
const CGS, CGS2019, CGSm, CGSe, HLU = Gauss, EMU2019, EMU, ESU, LorentzHeaviside

# natural units

const Planck = UnitSystem{1,1,1,1,cache(√(4π*αG))}()
const PlanckGauss = UnitSystem{1,1,1,4π,cache(√αG)}()
const Stoney = UnitSystem{1,cache(αinv),1,4π,cache(√(αG*αinv))}()
const Hartree = UnitSystem{1,1,cache(αinv),cache(4π/αinv^2),1}()
const Rydberg = UnitSystem{1,1,cache(2αinv),cache(π/αinv^2),1/2}()
const Schrodinger = UnitSystem{1,1,cache(αinv),cache(4π/αinv^2),cache(√(αG*αinv))}()
const Electronic = UnitSystem{1,cache(αinv),1,4π,1}()
const Natural = UnitSystem{1,1,1,1,1}()
const NaturalGauss = UnitSystem{1,1,1,4π,1}()
const QCD = UnitSystem{1,1,1,1,cache(1/μₚₑ)}()
const QCDGauss = UnitSystem{1,1,1,4π,cache(1/μₚₑ)}()
const QCDoriginal = UnitSystem{1,1,1,cache(4π/αinv),cache(1/μₚₑ)}()

# physical constants

@pure electronmass(U::typeof(Planck),C::Coupling) = sqrt(4π*coupling(C))
@pure electronmass(U::typeof(PlanckGauss),C::Coupling) = sqrt(coupling(C))
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,cache(√(αG*αinv))},C::Coupling) where {kB,ħ,𝘤,μ₀} = sqrt(coupling(C)/finestructure(C))
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,cache(1/μₚₑ)},C::Coupling) where {kB,ħ,𝘤,μ₀} = 1/protonelectron(C)
@pure permeability(U::UnitSystem{kB,ħ,𝘤,cache(4π/αinv^2)},C::Coupling) where {kB,ħ,𝘤} = 4π*finestructure(C)^2
@pure permeability(U::UnitSystem{kB,ħ,𝘤,cache(π/αinv^2)},C::Coupling) where {kB,ħ,𝘤} = π*finestructure(C)^2
@pure lightspeed(U::UnitSystem{kB,ħ,cache(αinv)},C::Coupling) where {kB,ħ} = 1/finestructure(C)
@pure lightspeed(U::UnitSystem{kB,ħ,cache(2αinv)},C::Coupling) where {kB,ħ} = 2/finestructure(C)
@pure planckreduced(U::UnitSystem{kB,cache(αinv)},C::Coupling) where kB = 1/finestructure(C)

@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,cache(mₑ)},C::Coupling) where {kB,ħ,μ₀} = electronmass(planck(U),R∞,C)
@pure electronmass(U::UnitSystem{kB,ħ,100𝘤,μ₀,cache(1000mₑ)},C::Coupling) where {kB,ħ,μ₀} = 1000electronmass(SI,C)
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,cache(mₑ/1000)},C::Coupling) where {kB,ħ,μ₀} = electronmass(SI,C)/1000
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,cache(mₑ2014)},C::Coupling) where {kB,ħ,μ₀} = electronmass(planck(U),R∞,C)
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,cache(mₑ1990)},C::Coupling) where {kB,ħ,μ₀} = electronmass(planck(U),R∞,C)
@pure electronmass(U::UnitSystem{kB,ħ,𝘤/ftUS,μ₀,cache(mₑ/slug)},C::Coupling) where {kB,ħ,μ₀} = electronmass(SI,C)/slug
@pure permeability(U::UnitSystem{kB,ħ,𝘤,cache(μ₀)},C::Coupling) where {kB,ħ,𝘤} = finestructure(C)*2𝘩/𝘤/𝘦^2
@pure permeability(U::typeof(ESU2019),C::Coupling) = 1e3*permeability(SI,C)/𝘤^2
@pure permeability(U::typeof(EMU2019),C::Coupling) = 1e7*permeability(SI,C)
@pure permeability(U::typeof(CODATA),C::Coupling) = 2RK2014*finestructure(C)/𝘤
@pure permeability(U::typeof(Conventional),C::Coupling) = 2RK1990*finestructure(C)/𝘤

@pure molarmass(U::UnitSystem{kB},C::Coupling=universe(U)) = NA*electronmass(U,C)/electronunit(C)
@pure molarmass(U::UnitSystem{1e7*kB},C::Coupling=universe(U)) = 1000molarmass(SI2019,C)
@pure molarmass(U::UnitSystem{1e3*kB},C::Coupling=universe(U)) = molarmass(SI2019,C)/1000
@pure molarmass(U::UnitSystem{cache(boltzmann(MTS))}) = molarmass(CGS)/1e6
@pure molarmass(U::UnitSystem{cache(boltzmann(CGS))}) = molarmass(Natural)
@pure molarmass(U::UnitSystem{cache(boltzmann(FFF))}) = molarmass(Natural)
@pure molarmass(U::UnitSystem{boltzmann(English)},C::Coupling=universe(U)) = 1000molarmass(SI2019,C)
@pure molarmass(U::UnitSystem{cache(boltzmann(EnglishUS))}) = molarmass(Natural)
@pure molarmass(U::UnitSystem{cache(boltzmann(IAU))}) = 1/1000mₛ

@pure luminousefficacy(U::UnitSystem{1}) = 1
@pure luminousefficacy(U::UnitSystem) = power(Kcd,SI2019,U)

@pure universe(::typeof(FFF)) = Universe
for unit ∈ UnitSystems.Systems
    if unit ∉ (:Natural,:NaturalGauss)
        @eval universe(::typeof($unit)) = Universe
    end
end

const dir = dirname(pathof(UnitSystems))
include("$dir/kinematic.jl")
include("$dir/electromagnetic.jl")
include("$dir/thermodynamic.jl")
include("$dir/physics.jl")
include("$dir/systems.jl")

end # module


text:
function createcluster!(ising::Ising, site::CartesianIndex, spin_val::Int)
    ising.cluster[site] = true
    ising.spins[site] = -spin_val

    for nn_table in ising.neighbours_table
        new_site = nn_table[site]
        if ising.spins[new_site] ≡ -ising.spins[site] && ising.cluster[new_site] ≡ false
            if rand() < 1 - exp(-2 * ising.β)
                createcluster!(ising, new_site, spin_val)
            end
        end
    end
end

function wolff!(ising::Ising)
    fill!(ising.cluster, false)
    rand_site = rand(ising.sites)
    spin_val = ising.spins[rand_site]
    createcluster!(ising, rand_site, spin_val)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

real-valued complex polynomials
=#

import Hypatia.PolyUtils

# get interpolation for a predefined complex poly
function get_interp_data(
    R::Type{Complex{T}},
    poly_name::Symbol,
    halfdeg::Int,
    ) where {T <: Real}
    (n, f, gs, g_halfdegs, true_min) = complex_poly_data[poly_name]
    (points, Ps) = PolyUtils.interpolate(R, halfdeg, n, gs, g_halfdegs)
    interp_vals = f.(points)
    return (interp_vals, Ps, true_min)
end

# merge with real polys when complex polyvars are allowed in DynamicPolynomials:
# real-valued complex polynomials
complex_poly_data = Dict{Symbol, NamedTuple}(
    :abs1d => (n = 1,
        f = (z -> 1 + sum(abs2, z)),
        gs = Function[],
        g_halfdegs = Int[],
        true_min = 1,
        ),
    :absunit1d => (n = 1,
        f = (z -> 1 + sum(abs2, z)),
        gs = [z -> 1 - sum(abs2, z)],
        g_halfdegs = [1],
        true_min = 1,
        ),
    :negabsunit1d => (n = 1,
        f = (z -> -sum(abs2, z)),
        gs = [z -> 1 - sum(abs2, z)],
        g_halfdegs = [1],
        true_min = -1,
        ),
    :absball2d => (n = 2,
        f = (z -> 1 + sum(abs2, z)),
        gs = [z -> 1 - sum(abs2, z)],
        g_halfdegs = [1],
        true_min = 1,
        ),
    :absbox2d => (n = 2,
        f = (z -> 1 + sum(abs2, z)),
        gs = [z -> 1 - abs2(z[1]), z -> 1 - abs2(z[2])],
        g_halfdegs = [1, 1],
        true_min = 1,
        ),
    :negabsbox2d => (n = 2,
        f = (z -> -sum(abs2, z)),
        gs = [z -> 1 - abs2(z[1]), z -> 1 - abs2(z[2])],
        g_halfdegs = [1, 1],
        true_min = -2,
        ),
    :denseunit1d => (n = 1,
        f = (z -> 1 + 2real(z[1]) + abs(z[1])^2 + 2real(z[1]^2) +
            2real(z[1]^2 * conj(z[1])) + abs(z[1])^4),
        gs = [z -> 1 - abs2(z[1])],
        g_halfdegs = [1],
        true_min = 0,
        ),
    )


text:
module TestUtil

include("preamble.jl")
using LyceumCore: genpath

@testset "mkgoodpath" begin
    with_tempdir() do dir
        mkpath("dir1/dir2")
        touch("f1.txt")
        touch(joinpath("dir1", "f2.txt"))
        touch(joinpath("dir1", "dir2", "f3.txt"))
        touch(joinpath("dir1", "dir2", "f3_1.txt"))

        @test genpath("dir1/dir2") == joinpath("dir1", "dir2_1")
        @test genpath("f1.txt", force=true) == "f1.txt"
        @test !isfile("f1.txt")
        @test genpath("dir1/dir2/f3.txt") == joinpath("dir1", "dir2", "f3_2.txt")
        @test genpath("dir1/dir2/f3.txt", sep="___") == joinpath("dir1", "dir2", "f3___1.txt")
    end
end

end # module


text:
macro submodel(expr)
    return quote
        _evaluate($(esc(expr)), $(esc(:__varinfo__)), $(esc(:__context__)))
    end
end

macro submodel(prefix, expr)
    return quote
        _evaluate(
            $(esc(expr)),
            $(esc(:__varinfo__)),
            PrefixContext{$(esc(Meta.quot(prefix)))}($(esc(:__context__))),
        )
    end
end


text:
# ---
# title: 132. Palindrome Partitioning II
# id: problem132
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/palindrome-partitioning-ii/description/>
# hidden: true
# ---
#
# Given a string `s`, partition `s` such that every substring of the partition
# is a palindrome.
#
# Return _the minimum cuts needed_ for a palindrome partitioning of `s`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aab"
#     Output: 1
#     Explanation: The palindrome partitioning ["aa","b"] could be produced using 1 cut.
#
#
# **Example 2:**
#
#
#
#     Input: s = "a"
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: s = "ab"
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 2000`
#   * `s` consists of lower-case English letters only.
#
#
## @lc code=start
using LeetCode

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


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct SecurityGroupNetworkInterface <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    securityRuleAssociations::Any # spec type: Union{ Nothing, SecurityRuleAssociations } # spec name: securityRuleAssociations

    function SecurityGroupNetworkInterface(;id=nothing, securityRuleAssociations=nothing)
        o = new()
        validate_property(SecurityGroupNetworkInterface, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(SecurityGroupNetworkInterface, Symbol("securityRuleAssociations"), securityRuleAssociations)
        setfield!(o, Symbol("securityRuleAssociations"), securityRuleAssociations)
        o
    end
end # type SecurityGroupNetworkInterface

const _property_map_SecurityGroupNetworkInterface = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("securityRuleAssociations")=>Symbol("securityRuleAssociations"))
const _property_types_SecurityGroupNetworkInterface = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("securityRuleAssociations")=>"SecurityRuleAssociations")
Base.propertynames(::Type{ SecurityGroupNetworkInterface }) = collect(keys(_property_map_SecurityGroupNetworkInterface))
Swagger.property_type(::Type{ SecurityGroupNetworkInterface }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_SecurityGroupNetworkInterface[name]))}
Swagger.field_name(::Type{ SecurityGroupNetworkInterface }, property_name::Symbol) =  _property_map_SecurityGroupNetworkInterface[property_name]

function check_required(o::SecurityGroupNetworkInterface)
    true
end

function validate_property(::Type{ SecurityGroupNetworkInterface }, name::Symbol, val)
end


text:
abstract type AbstractCCARule <: AbstractODRule end

struct CCA{T} <: AbstractCCARule
    rule::T
    radius::Int
end

"""
    TCA(code; radius=1)

Returns a ```CCA``` object given a specific code and radius.
"""
function CCA(rule; radius=1)
    CCA(rule, radius)
end

function (cca::CCA)(starting_array)

    nextgen = evolution(starting_array, cca.rule, cca.radius)
end

function c_state_reader(neighborhood, radius)

    sum(neighborhood)/length(neighborhood)
end

function evolution(cell, rule, radius::Number)

    neighborhood_size = radius*2+1
    output = zeros(length(cell))
    cell = vcat(cell[end-neighborhood_size÷2+1:end], cell, cell[1:neighborhood_size÷2])

    for i=1:length(cell)-neighborhood_size+1
        output[i] = modf(c_state_reader(cell[i:i+neighborhood_size-1], radius)+rule)[1]
    end

    output
end

function evolution(cell, rule, radius::Tuple)

    neighborhood_size = sum(radius)+1
    output = zeros(length(cell))
    cell = vcat(cell[end-neighborhood_size÷2+1:end], cell, cell[1:neighborhood_size÷2])

    for i=1:length(cell)-neighborhood_size+1
        output[i] = modf(c_state_reader(cell[i:i+neighborhood_size-1], radius)+rule)[1]
    end

    output
end


text:
include("get_div_grad.jl")

function test_symmlq()
  symmlq_tol = 1.0e-5

  # 1. Symmetric and positive definite systems.
  #
  # Cubic spline matrix.
  n = 10;
  A = spdiagm(-1 => ones(n-1), 0 => 4*ones(n), 1 => ones(n-1))
  b = A * [1:n;]

  (x, xcg, stats) = symmlq(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("SYMMLQ: Relative residual: %8.1e\n", resid)
  @test(resid <= symmlq_tol)
  @test(stats.solved)

  # Symmetric indefinite variant.
  A = A - 3 * I
  b = A * [1:n;]
  (x, xcg, stats) = symmlq(A, b, itmax=n+1)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("SYMMLQ: Relative residual: %8.1e\n", resid)
  @test(resid <= symmlq_tol)
  @test(stats.solved)

  # Code coverage.
  (x, xcg, stats) = symmlq(Matrix(A), b)
  show(stats)

  # Sparse Laplacian (CG point will terminate sooner).
  A = get_div_grad(16, 16, 16)
  b = ones(size(A, 1))
  (x, xcg, stats) = symmlq(A, b, atol=1e-12, rtol=1e-12)
  r = b - A * x
  resid = norm(r) / norm(b)
  @show stats
  @printf("SYMMLQ: Relative residual: %8.1e\n", resid)
  @test(resid <= symmlq_tol)
  @test(stats.solved)

  # Symmetric indefinite variant, almost singular.
  A = A - 5 * I
  (x, xcg, stats) = symmlq(A, b)
  r = b - A * x
  resid = norm(r) / norm(b)
  @printf("SYMMLQ: Relative residual: %8.1e\n", resid)
  @test(resid <= 100 * symmlq_tol)
  @test(stats.solved)

  # Test b == 0
  (x, xcg, stats) = symmlq(A, zeros(size(A,1)))
  @test x == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = spdiagm(-1 => ones(Int, n-1), 0 => 4*ones(Int, n), 1 => ones(Int, n-1))
  b = A * [1:n;]
  (x, xcg, stats) = symmlq(A, b)
  @test stats.solved

  # Test error estimate
  A = Matrix(get_div_grad(8, 8, 8))
  b = ones(size(A, 1))
  λest = (1-1e-10)*eigmin(A)
  x_exact = A\b
  (x, xcg, stats) = symmlq(A, b, λest=λest)
  err = norm(x_exact - x)
  errcg = norm(x_exact - xcg)
  @printf("SYMMLQ    : true error: %8.1e\n", err)
  @printf("SYMMLQ-CG : true error: %8.1e\n", errcg)
  @printf("SYMMLQ    : err up-bnd : %8.1e\n", stats.errors[end])
  @printf("SYMMLQ-CG : err up-bnd : %8.1e\n", stats.errorscg[end])
  @test( err <= stats.errors[end] )
  @test( errcg <= stats.errorscg[end])
  (x, xcg, stats) = symmlq(A, b, λest=λest, window=1)
  @printf("SYMMLQ    : err up-bnd : %8.1e\n", stats.errors[end])
  @printf("SYMMLQ-CG : err up-bnd : %8.1e\n", stats.errorscg[end])
  @test( err <= stats.errors[end] )
  @test( errcg <= stats.errorscg[end])
  (x, xcg, stats) = symmlq(A, b, λest=λest, window=5)
  @printf("SYMMLQ    : err up-bnd : %8.1e\n", stats.errors[end])
  @printf("SYMMLQ-CG : err up-bnd : %8.1e\n", stats.errorscg[end])
  @test( err <= stats.errors[end] )
  @test( errcg <= stats.errorscg[end])
end

test_symmlq()


text:
using LinearAlgebra
using Polynomials
using SemiLagrangian: Hermite, PrecalHermite, L, Lprim, K, H, bplus, bminus, _getpolylagrange

function getbp(i,rp,sp)
    res=big(1//1)
    for j=rp:sp
        if j != i
            res /= (i-j)
            if j != 0
                res *= (-j)
            end
        end
    end
    return res
end

function test_precalhermite(ord)
    ph = PrecalHermite(ord)
    d = div(ord,2)
    @test ph.rplus == -d
    @test ph.splus == d+1
    @test ph.rminus == -d-1
    @test ph.sminus == d

    sbp=0//1
    sbm=0//1

    for i=-d:d+1
        @test L(ph,i) == _getpolylagrange(i+d,ord,-d)
        @test Lprim(ph, i) == derivative(L(ph,i))(i)
        @test K(ph, i) == _getpolylagrange(i+d,ord,-d)^2 *Polynomial([-i, 1//1])
        @test H(ph, i) == _getpolylagrange(i+d,ord,-d)^2 * (1-2*Lprim(ph,i)*Polynomial([-i,1//1]))
        if i != 0
            @test bplus(ph, i) == getbp(i,ph.rplus,ph.splus)
            @test bminus(ph, -i) == -getbp(i, ph.rplus, ph.splus)
        end
        sbp += bplus(ph,i)
        sbm += bminus(ph, -i)
    end

    @test sbp == 0
    @test sbm == 0

    if ord == 3
        @test ph.bplus == [-1//3, -1//2, 1, -1//6]
    elseif ord == 5
        @test ph.bplus == [1//20,-1//2,-1//3,1//1,-1//4,1//30]
    end

end
function test_base_hermite(order)
    herm = Hermite(order, Rational{BigInt})
    ord = div(order,2)+1
    tab = rationalize.(BigInt, rand(ord+1), tol = 1 / 1000000)
    fct = Polynomial(tab)
    dec = div(order, 2)
    for i = 1:3
        value = rationalize.(BigInt, rand(), tol = 1 / 1000000)
        res = 0 // 1
        for j = 0:order
            res += herm.tabfct[j+1](value) * fct(j - dec)
        end
        @test res == fct(value)
    end
end
function test_base_hermite2d(order)
    herm = Hermite(order, Rational{BigInt})
    ord = div(order,2)+1
    tab = rationalize.(BigInt, rand(ord+1, ord+1), tol = 1 / 1000000)
    fct = Pol2(tab)
    dec = div(order, 2)
    for i = 1:1
        val_x = rationalize.(BigInt, rand(), tol = 1 / 1000000)
        val_y = rationalize.(BigInt, rand(), tol = 1 / 1000000)
        res = 0 // 1
        for j = 0:order, k = 0:order
            res += herm.tabfct[j+1](val_x) * herm.tabfct[k+1](val_y) * fct(j - dec, k - dec)
        end
        resf = fct(val_x, val_y)
        #        println("order = $order res-fct=$(res-resf)")
        @test res == resf
    end
end
@time @testset "Hermite precal" begin
    for ord=3:2:21
        test_precalhermite(ord)
    end
end
@time @testset "test base hermite interpolation" begin
    for order = 5:4:41
        test_base_hermite(order)
    end
end
@time @testset "test base hermite interpolation2d" begin
    for order = 5:4:41
        test_base_hermite2d(order)
    end
end


text:
x = 20

println(2x+2)
println(2(x)+2)
println(2(x+1))
println((2)x+2)

println(2+3x)
println(2*3x)
println(2^3x)
println(2(div(x,2)))

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc})   # time: 0.32043454
    Base.precompile(Tuple{typeof(org),IOBuffer,Table,Int64})   # time: 0.27505028
    Base.precompile(Tuple{typeof(org),IOBuffer,OrgDoc})   # time: 0.24133022
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,UnorderedList,Int64})   # time: 0.16556492
    Base.precompile(Tuple{typeof(consume),Type{Planning},SubString{String}})   # time: 0.14544973
    let fbody = try __lookup_kwbody__(which(parseorg, (SubString{String},Dict{Char, Vector{DataType}},Vector{DataType},))) catch missing end
    if !ismissing(fbody)
        precompile(fbody, (Bool,Bool,Int64,typeof(parseorg),SubString{String},Dict{Char, Vector{DataType}},Vector{DataType},))
    end
end   # time: 0.14305766
    Base.precompile(Tuple{typeof(html),IOBuffer,OrgDoc})   # time: 0.11243505
    Base.precompile(Tuple{Type{Table},Vector{Union{Nothing, SubString{String}}}})   # time: 0.09191036
    Base.precompile(Tuple{typeof(consume),Type{Timestamp},SubString{String}})   # time: 0.09087488
    Base.precompile(Tuple{typeof(consume),Type{Citation},SubString{String}})   # time: 0.069843575
    Base.precompile(Tuple{typeof(org),IOBuffer,UnorderedList,Int64})   # time: 0.0667007
    Base.precompile(Tuple{typeof(postprocess!),OrgDoc})   # time: 0.0661829
    Base.precompile(Tuple{Type{Section},Vector{Union{Nothing, SubString{String}}}})   # time: 0.055297952
    Base.precompile(Tuple{typeof(org),IOBuffer,AffiliatedKeywordsWrapper,Int64})   # time: 0.054353792
    Base.precompile(Tuple{typeof(term),IOBuffer,TextMarkup{SubString{String}},Vector{String}})   # time: 0.051143266
    Base.precompile(Tuple{typeof(term),IOContext{IOBuffer},OrgDoc,TextMarkup{Vector{Object}}})   # time: 0.045892756
    Base.precompile(Tuple{typeof(layouttable),IOBuffer,Function,Table,Dict{Char, Char},Int64})   # time: 0.04563425
    Base.precompile(Tuple{typeof(parse),Type{OrgDoc},String})   # time: 0.044654246
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,Table,Int64})   # time: 0.04446829
    Base.precompile(Tuple{Type{TableRow},Vector{Union{Nothing, SubString{String}}}})   # time: 0.041455567
    Base.precompile(Tuple{typeof(org),IOContext{IOBuffer},TextMarkup{Vector{Object}}})   # time: 0.038945496
    Base.precompile(Tuple{typeof(term),IOContext{IOBuffer},OrgDoc,SourceBlock,Int64})   # time: 0.038509842
    Base.precompile(Tuple{typeof(consume),Type{TextPlainForced},SubString{String}})   # time: 0.036461886
    let fbody = try __lookup_kwbody__(which(parseorg, (SubString{String},Dict{Char, Vector{<:Type}},Vector{Type},))) catch missing end
    if !ismissing(fbody)
        precompile(fbody, (Bool,Bool,Int64,typeof(parseorg),SubString{String},Dict{Char, Vector{<:Type}},Vector{Type},))
    end
end   # time: 0.035671674
    Base.precompile(Tuple{typeof(term),IOContext{IOBuffer},OrgDoc,LaTeXFragment})   # time: 0.035494458
    Base.precompile(Tuple{typeof(consume),Type{AffiliatedKeywordsWrapper},SubString{String}})   # time: 0.03425677
    Base.precompile(Tuple{typeof(org),IOContext{IOBuffer},UnorderedList,Int64})   # time: 0.031304054
    Base.precompile(Tuple{typeof(consume),Type{TextPlain},SubString{String}})   # time: 0.030219704
    Base.precompile(Tuple{typeof(org),IOBuffer,Paragraph,Int64})   # time: 0.029114475
    Base.precompile(Tuple{typeof(consume),Type{Comment},SubString{String}})   # time: 0.024979195
    Base.precompile(Tuple{typeof(consume),Type{AffiliatedKeyword},SubString{String}})   # time: 0.02246017
    isdefined(Org, Symbol("#98#104")) && Base.precompile(Tuple{getfield(Org, Symbol("#98#104")),Vector{String}})   # time: 0.021567201
    let fbody = try __lookup_kwbody__(which(parseobjects, (Type,SubString{String},))) catch missing end
    if !ismissing(fbody)
        precompile(fbody, (Bool,Bool,typeof(parseobjects),Type,SubString{String},))
    end
end   # time: 0.020833066
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,Keyword{SubString{String}},Int64})   # time: 0.017364282
    Base.precompile(Tuple{typeof(consume),Type{LaTeXFragment},SubString{String}})   # time: 0.01619035
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,SourceBlock,Int64})   # time: 0.01617535
    Base.precompile(Tuple{typeof(iterate),OrgIterator})   # time: 0.014937996
    Base.precompile(Tuple{Type{PropertyDrawer},Vector{Union{Nothing, SubString{String}}}})   # time: 0.014385041
    isdefined(Org, Symbol("#itemconsume#59")) && Base.precompile(Tuple{getfield(Org, Symbol("#itemconsume#59")),SubString{String},SubString{String}})   # time: 0.013828709
    Base.precompile(Tuple{typeof(consume),Type{FootnoteDefinition},SubString{String}})   # time: 0.013215547
    Base.precompile(Tuple{typeof(consume),Type{TextMarkup},SubString{String}})   # time: 0.012036697
    Base.precompile(Tuple{typeof(html),IOBuffer,ExportBlock})   # time: 0.011866605
    Base.precompile(Tuple{typeof(org),IOBuffer,AffiliatedKeyword{Nothing},Int64})   # time: 0.011750292
    Base.precompile(Tuple{typeof(org),IOBuffer,RegularLink})   # time: 0.010200925
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,SpecialBlock,Int64})   # time: 0.010012363
    isdefined(Org, Symbol("#100#106")) && Base.precompile(Tuple{getfield(Org, Symbol("#100#106")),Vector{Vector{Int64}}})   # time: 0.0099559
    Base.precompile(Tuple{typeof(consume),Type{InlineSourceBlock},SubString{String}})   # time: 0.009951612
    Base.precompile(Tuple{typeof(consume),Type{FootnoteReference},SubString{String}})   # time: 0.009520893
    isdefined(Org, Symbol("#100#106")) && Base.precompile(Tuple{getfield(Org, Symbol("#100#106")),Vector{Union{Nothing, Vector{Int64}}}})   # time: 0.009133301
    Base.precompile(Tuple{typeof(org),IOBuffer,SourceBlock,Int64})   # time: 0.008511166
    Base.precompile(Tuple{typeof(consume),Type{RegularLink},SubString{String}})   # time: 0.007196862
    Base.precompile(Tuple{typeof(consume),Type{List},SubString{String}})   # time: 0.007170842
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,QuoteBlock,Int64})   # time: 0.00685181
    Base.precompile(Tuple{typeof(consume),Type{Item},SubString{String}})   # time: 0.00668202
    Base.precompile(Tuple{typeof(org),IOContext{IOBuffer},RegularLink})   # time: 0.006290038
    Base.precompile(Tuple{Type{Block},Vector{Union{Nothing, SubString{String}}}})   # time: 0.006034944
    Base.precompile(Tuple{typeof(term),IOContext{IOBuffer},OrgDoc,RegularLink})   # time: 0.005997927
    Base.precompile(Tuple{typeof(consume),Type{GreaterBlock},SubString{String}})   # time: 0.005965764
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,AffiliatedKeywordsWrapper,Int64})   # time: 0.005771555
    Base.precompile(Tuple{Type{Keyword},Vector{Union{Nothing, SubString{String}}}})   # time: 0.005761046
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,CenterBlock,Int64})   # time: 0.005582146
    isdefined(Org, Symbol("#124#125")) && Base.precompile(Tuple{getfield(Org, Symbol("#124#125")),IOBuffer,LaTeXFragment})   # time: 0.005432231
    Base.precompile(Tuple{Type{Heading},Vector{Union{Nothing, SubString{String}}}})   # time: 0.005305904
    isdefined(Org, Symbol("#textobjupto#70")) && Base.precompile(Tuple{getfield(Org, Symbol("#textobjupto#70")),Int64})   # time: 0.005131729
    Base.precompile(Tuple{typeof(org),IOContext{IOBuffer},SourceBlock,Int64})   # time: 0.005061035
    Base.precompile(Tuple{typeof(html),IOBuffer,Table})   # time: 0.005050735
    Base.precompile(Tuple{typeof(html),IOBuffer,UnorderedList})   # time: 0.004805607
    Base.precompile(Tuple{typeof(parse),Type{LinkPath},SubString{String},Bool})   # time: 0.00473084
    Base.precompile(Tuple{typeof(consume),Type{Clock},SubString{String}})   # time: 0.004381719
    Base.precompile(Tuple{typeof(consume),Type{Entity},SubString{String}})   # time: 0.003819155
    Base.precompile(Tuple{typeof(org),IOBuffer,QuoteBlock,Int64})   # time: 0.003732675
    Base.precompile(Tuple{typeof(iterate),OrgIterator,Vector{Tuple}})   # time: 0.003555848
    let fbody = try __lookup_kwbody__(which(forwardsbalenced, (SubString{String},Int64,Int64,))) catch missing end
    if !ismissing(fbody)
        precompile(fbody, (Dict{Char, Char},Vector{Char},Vector{Char},Vector{Char},typeof(forwardsbalenced),SubString{String},Int64,Int64,))
    end
end   # time: 0.003452854
    Base.precompile(Tuple{typeof(org),IOBuffer,ExportBlock,Int64})   # time: 0.003173403
    Base.precompile(Tuple{typeof(html),IOBuffer,TextMarkup{SubString{String}}})   # time: 0.002870525
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,Paragraph,Int64})   # time: 0.002809988
    Base.precompile(Tuple{typeof(org),IOBuffer,Drawer,Int64})   # time: 0.00278488
    Base.precompile(Tuple{typeof(term),IOBuffer,OrgDoc,ExportBlock,Int64})   # time: 0.002739602
    Base.precompile(Tuple{typeof(html_tagpair),String,Pair{String, String}})   # time: 0.002626423
    Base.precompile(Tuple{typeof(consume),Type{ExportSnippet},SubString{String}})   # time: 0.002619359
    Base.precompile(Tuple{typeof(term),IOContext{IOBuffer},OrgDoc,Paragraph,Int64})   # time: 0.002616784
    Base.precompile(Tuple{typeof(org),IOBuffer,LaTeXFragment})   # time: 0.002530007
    Base.precompile(Tuple{typeof(org),IOBuffer,Keyword{SubString{String}},Int64})   # time: 0.002526863
    Base.precompile(Tuple{typeof(org),IOBuffer,TextMarkup{SubString{String}}})   # time: 0.002374013
    Base.precompile(Tuple{typeof(consume),Type{Drawer},SubString{String}})   # time: 0.002340451
    Base.precompile(Tuple{typeof(org),IOBuffer,ExportSnippet})   # time: 0.002330339
    Base.precompile(Tuple{typeof(html),IOBuffer,RegularLink})   # time: 0.002330261
    Base.precompile(Tuple{typeof(consume),Type{NodeProperty},SubString{String}})   # time: 0.002107064
    Base.precompile(Tuple{typeof(term),IOBuffer,RegularLink,Vector{String}})   # time: 0.002088849
    Base.precompile(Tuple{typeof(org),IOBuffer,Comment,Int64})   # time: 0.002087208
    isdefined(Org, Symbol("#65#71")) && Base.precompile(Tuple{getfield(Org, Symbol("#65#71")),Char})   # time: 0.001969018
    Base.precompile(Tuple{typeof(html),IOBuffer,Keyword{SubString{String}}})   # time: 0.001919663
    Base.precompile(Tuple{typeof(org),IOBuffer,ExampleBlock,Int64})   # time: 0.001773046
    Base.precompile(Tuple{typeof(html),IOBuffer,SourceBlock})   # time: 0.001771655
    Base.precompile(Tuple{typeof(org),IOContext{IOBuffer},InlineSourceBlock})   # time: 0.001747723
    Base.precompile(Tuple{typeof(html),IOBuffer,PlainLink})   # time: 0.001667706
    Base.precompile(Tuple{typeof(term),IOBuffer,TextPlain{SubString{String}},Vector{String}})   # time: 0.00163107
    Base.precompile(Tuple{typeof(html),IOBuffer,TextMarkup{Vector{Object}}})   # time: 0.001566183
    Base.precompile(Tuple{typeof(term),IOContext{IOBuffer},OrgDoc,TextPlain{SubString{String}}})   # time: 0.001551615
    Base.precompile(Tuple{typeof(term),IOContext{IOBuffer},OrgDoc,TextPlain{String}})   # time: 0.001542971
    Base.precompile(Tuple{typeof(html),IOBuffer,TextPlain{SubString{String}}})   # time: 0.001469463
    Base.precompile(Tuple{typeof(parse),Type{TableRow},SubString{String}})   # time: 0.00145073
    Base.precompile(Tuple{typeof(html),IOBuffer,ExampleBlock})   # time: 0.001345151
    Base.precompile(Tuple{typeof(parse),Type{Section},SubString{String}})   # time: 0.001341194
    Base.precompile(Tuple{typeof(html),IOBuffer,TextPlain{String}})   # time: 0.001281813
    Base.precompile(Tuple{typeof(org),IOBuffer,TextMarkup{Vector{Object}}})   # time: 0.001265663
    isdefined(Org, Symbol("#124#125")) && Base.precompile(Tuple{getfield(Org, Symbol("#124#125")),IOBuffer,Entity})   # time: 0.001249151
    Base.precompile(Tuple{typeof(html),IOBuffer,SpecialBlock})   # time: 0.001218596
    Base.precompile(Tuple{typeof(html),IOBuffer,LaTeXFragment})   # time: 0.00120709
    Base.precompile(Tuple{typeof(html),IOBuffer,DescriptiveList})   # time: 0.001175334
    Base.precompile(Tuple{typeof(html),IOBuffer,OrderedList})   # time: 0.001161889
    Base.precompile(Tuple{typeof(iterate),Item,Int64})   # time: 0.001062694
    Base.precompile(Tuple{typeof(*),TextPlain{String},TextPlain{SubString{String}}})   # time: 0.001045501
end

_precompile_()


text:
"""
    readELP(fname::AbstractString)

Import electrode data in the BESA elp format`

#### Arguments

* `fname` = location and name of file to import
* `verbose` = should the function output user feedback

name, a, b = readELP(verbose=false)
"""
function readELP(fname::AbstractString = Pkg.dir("Private", "test", "data", "Standard-10-10-Cap81.elp"); verbose::Bool=false)

    debug("Electrode file: $(fname)")

    df = readtable(fname, header = false, separator = ' ')

    if verbose
        println("Electrodes:  $(length(df[:x1]))")
    end

    return df[:x1].data, df[:x2].data, df[:x3].data
end


text:
module Day11

using SparseArrays
using ..Utils

neighbors_matrix(n, m) = begin
    bidiag(k) = spdiagm(-1 => ones(Bool, k-1), 1 => ones(Bool, k-1))
    tridiag(k) = spdiagm(-1 => ones(Bool, k-1), 0 => ones(Bool, k), 1 => ones(Bool, k-1))

    kron(bidiag(m), tridiag(n)) .| kron(tridiag(m), bidiag(n))
end

function step!(octopuses::Matrix{Int8})
    n, m = size(octopuses)
    octopuses = vec(octopuses)
    neighbors = neighbors_matrix(n, m)

    # Increase each by 1
    octopuses .+= 1

    # Any which are above 9
    already_flashed = spzeros(Bool, length(octopuses))

    flashing = sparse(octopuses .> 9)
    already_flashed .|= flashing
    while any(flashing)
        octopuses .+= neighbors * flashing
        flashing = sparse(octopuses .> 9) .⊻ already_flashed
        already_flashed .|= flashing
    end
    octopuses[already_flashed] .= 0
    return sum(already_flashed)
end

function part1(octopuses)
    sum = 0
    for i in 1:100
        sum += step!(octopuses)
    end
    return sum
end

function part2(octopuses)
    i = 0
    while true
        i += 1
        count = step!(octopuses)
        if count == length(octopuses)
            return i
        end
    end
end

function run()
    test = """5483143223
2745854711
5264556173
6141336146
6357385478
4167524645
2176841721
6882881134
4846848554
5283751526"""

    test_octopuses = parse_array(test)
    @assert part1(test_octopuses) == 1656
    @assert part2(parse_array(test)) == 195

    input = """4764745784
4643457176
8322628477
7617152546
6137518165
1556723176
2187861886
2553422625
4817584638
3754285662"""

    octopuses = parse_array(input)
    @show part1(octopuses)
    octopuses = parse_array(input)
    @show part2(octopuses)
end

end


text:
@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    const glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    const glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
end

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)