text:
module Grid

# export Log, UniformGrid, tauGrid, fermiKGrid, boseKGrid

using StaticArrays

@enum GridType LOG UNIFORM UNILOG

struct Coeff{T <: AbstractFloat}
    bound::SVector{2,T}
    idx::SVector{2,T}
    λ::T
    a::T
    b::T

    function Coeff{T}(bound, idx, λ, dense2sparse::Bool) where {T}
        λ /= abs(idx[2] - idx[1])
        if dense2sparse == false
            bound = (bound[2], bound[1])
            idx = (idx[2], idx[1])
            λ = -λ
        end
        _l1, _l2 = T(1.0), exp(λ * (idx[2] - idx[1]))
        b = (bound[2] - bound[1]) / (_l2 - _l1)
        a = (bound[1] * _l2 - bound[2] * _l1) / (_l2 - _l1)
        return new{T}(bound, idx, λ, a, b)
    end
end

struct UniLog{T <: AbstractFloat}
    bound::SVector{2,T}
    idx::SVector{2,T}
    isopen::SVector{2,Bool}
    M::Int
    N::Int
    λ::T
    d2s::Bool

    function UniLog{T}(bound, init, minterval::T,M::Int,N::Int, dense2sparse::Bool,isopen = @SVector[false,true]) where {T}
        @assert N*minterval<bound[2]-bound[1]
        Nidx = (M+1)*N
        if isopen[2]==true
            Nidx = Nidx-1
        end
        idx = @SVector[init, init + Nidx ]
        λ = (minterval*N/(bound[2]-bound[1]))^(1.0/M)
        return new{T}(bound, idx,isopen, M,N, λ, dense2sparse)
    end
end

function _grid(l::UniLog{T}, i) where {T}
    head = !l.isopen[1] ? l.idx[1] : (l.idx[1]-1)
    tail = !l.isopen[2] ? l.idx[2] : (l.idx[2]+1)
    if l.d2s
        i_n = (i - head)%l.N;
        i_m = (i - head)÷l.N;
        if i_m==0
            return l.bound[1]+l.λ^l.M/l.N*(l.bound[2]-l.bound[1])*i_n
        else
            return l.bound[1]+l.λ^(l.M+1-i_m)*(l.bound[2]-l.bound[1])+(l.λ^(l.M-i_m)-l.λ^(l.M-i_m+1))/l.N*i_n*(l.bound[2]-l.bound[1])
        end
    end
    if !l.d2s
        i_n = (tail - i)%l.N;
        i_m = (tail - i)÷l.N;
        if i_m==0
            return l.bound[2]-l.λ^(l.M)/l.N*(l.bound[2]-l.bound[1])*i_n
        else
            return l.bound[2]-l.λ^(l.M+1-i_m)*(l.bound[2]-l.bound[1])-(l.λ^(l.M-i_m)-l.λ^(l.M-i_m+1))/l.N*i_n*(l.bound[2]-l.bound[1])
        end
    end
end

function _floor(l::UniLog{T}, x::T) where {T}
    head = !l.isopen[1] ? l.idx[1] : (l.idx[1]-1)
    tail = !l.isopen[2] ? l.idx[2] : (l.idx[2]+1)
    length = l.bound[2] - l.bound[1]

    norm = l.d2s ? ((x - l.bound[1])/length) : ((l.bound[2] - x)/length)
    if norm <= 0
        pos = l.d2s ? l.idx[1] : l.idx[2]
        return Base.floor(Int,pos)
    elseif norm >=1
        pos = l.d2s ? l.idx[2] : l.idx[1]
        return Base.floor(Int,pos)
    end

    i_m= Base.floor(log(norm)/log(l.λ))
    if i_m>=l.M
        if l.d2s
            result = Base.floor(head+norm/l.λ^(l.M)*l.N)
            return (l.isopen[1]&&result==head) ? (result+1) : result
        else
            result = Base.floor(tail-norm/l.λ^(l.M)*l.N)
            return (l.isopen[2]&&result==tail) ? (result-1) : result
        end
    end
    i_n= Base.floor((norm-l.λ^(i_m+1))/(l.λ^(i_m)-l.λ^(i_m+1))*l.N)
    if l.d2s
        result = Base.floor(head+(l.M-i_m)*l.N+i_n)
        return (l.isopen[2]&&result==tail) ? result-1 : result
    else
        result = Base.floor((tail-(l.M-i_m)*l.N-i_n)-1)
        return (l.isopen[1]&&result==head) ? result+1 : result
    end
end

function _grid(l::Coeff{T}, pos) where {T}
    return l.a + l.b * exp(l.λ * (T(pos) - l.idx[1]))
end

function _floor(l::Coeff{T}, x::T) where {T}
    pos = l.idx[1] + log((x - l.a) / l.b) / l.λ
    return Base.floor(Int, pos)
end

function checkOrder(grid)
    for idx = 2:length(grid)
        @assert grid[idx - 1] < grid[idx] "The grid at $idx is not in the increase order: \n$grid"
    end
end

struct Log{T,SIZE,SEG} # create a log grid of the type T with SIZE grids and SEG of segments
    grid::MVector{SIZE,T}
    size::Int
    head::T
    tail::T
    segment::SVector{SEG,Int} # ends of each segments
    coeff::SVector{SEG,Coeff{T}}
    isopen::SVector{2,Bool}

    function Log{T,SIZE,SEG}(coeff, range, isopen) where {T <: AbstractFloat,SIZE,SEG}
        @assert SIZE > 1 "Size must be large than 1"
        for ri = 2:SEG
            @assert range[ri - 1][end] + 1 == range[ri][1] "ranges must be connected to each other"
        end
        @assert range[1][1] == 1 "ranges should start with the idx 1"
        @assert range[end][end] == SIZE "ranges ends with $(range[end][end]), expected $SIZE"

        grid, segment = [], []
        for s = 1:SEG
            push!(segment, range[s][end])
            for idx in range[s]
                push!(grid, _grid(coeff[s], idx))
            end
        end
        head, tail = grid[1], grid[end]
        isopen[1] && (grid[1] += eps(T) * 1e4)
        isopen[2] && (grid[end] -= eps(T) * 1e4)
        checkOrder(grid)
        return new{T,SIZE,SEG}(grid, SIZE, head, tail, segment, coeff, isopen)
    end
end

function Base.floor(grid::Log{T,SIZE,2}, x) where {T,SIZE}
    (grid.head <= x <= grid.tail) || error("$x is out of the uniform grid range!")
    segment = grid.segment
    if x < grid[2]
        return 1
    elseif x < grid[segment[1]]
        return _floor(grid.coeff[1], x)
    elseif x < grid[segment[1] + 1]
        return segment[1]
    elseif x < grid[end - 1]
        return _floor(grid.coeff[2], x)
    else
        return SIZE - 1
    end
end

function Base.floor(grid::Log{T,SIZE,3}, x) where {T,SIZE}
    (grid.head <= x <= grid.tail) || error("$x is out of the uniform grid range!")
    segment = grid.segment
    if x < grid[2]
        return 1
    elseif x < grid[segment[1]]
        return _floor(grid.coeff[1], x)
    elseif x < grid[segment[1] + 1]
        return segment[1]
    elseif x < grid[segment[2]]
        return _floor(grid.coeff[2], x)
    elseif x < grid[segment[2] + 1]
        return segment[2]
    elseif x < grid[end - 1]
        return _floor(grid.coeff[3], x)
    else
        return SIZE - 1
    end
end

Base.getindex(grid::Log, i) = grid.grid[i]
Base.firstindex(grid::Log) = 1
Base.lastindex(grid::Log) = grid.size

"""
    Uniform{Type,SIZE}

Create a uniform Grid with a given type and size

# Member:
- `β`: inverse temperature
- halfLife: the grid is densest in the range (0, halfLife) and (β-halfLife, β)
- size: the Grid size
- grid: vector stores the grid
- size: size of the grid vector
- head: grid head
- tail: grid tail
- δ: distance between two grid elements
- isopen: if isopen[1]==true, then grid[1] will be slightly larger than the grid head. Same for the tail.
"""
struct Uniform{T,SIZE}
    grid::SVector{SIZE,T}
    size::Int
    head::T
    tail::T
    δ::T
    isopen::SVector{2,Bool}

    """
        Uniform{Type,SIZE}(head, tail, isopen)

    Create a uniform Grid with a given type and size

    # Arguments:
     - head: the starting point of the grid
     - tail: the end of the grid
     - isopen: if isopen[1]==true, then grid[1]=head+eps; If isopen[2]==true, then grid[2]=tail-eps. Otherwise, grid[1]==head / grid[2]==tail
    """
    function Uniform{T,SIZE}(head, tail, isopen) where {T <: AbstractFloat,SIZE}
        @assert SIZE > 1 "Size must be large than 1"
        grid = Array(LinRange(T(head), T(tail), SIZE))
        isopen[1] && (grid[1] += eps(T))
        isopen[2] && (grid[end] -= eps(T))
        return new{T,SIZE}(grid, SIZE, head, tail, (tail - head) / (SIZE - 1), isopen)
end
end

function Base.floor(grid::Uniform, x)

    (grid.head <= x <= grid.tail) || error("$x is out of the uniform grid range!")

    if grid[2] <= x < grid[end - 1]
        return floor(Int, (x - grid.head) / grid.δ + 1)
    elseif x < grid[2]
        return 1
    else
        return grid.size - 1
end
end

Base.getindex(grid::Uniform, i) = grid.grid[i]
Base.firstindex(grid::Uniform) = 1
Base.lastindex(grid::Uniform) = grid.size

struct UniLogs{T<:AbstractFloat,SIZE,SEG}
    grid::MVector{SIZE,T}
    size::Int
    head::T
    tail::T
    unilogs::SVector{SEG,UniLog{T}}
    segment::SVector{SEG,T} # ends of each segments
    segindex::SVector{SEG,Int}
    isopen::SVector{2,Bool}

    function UniLogs{T,SIZE,SEG}(bounds, minterval::T,M::Int,N::Int, Isopen = @SVector[true,true], issparse = @SVector[false,false]) where {T<:AbstractFloat,SIZE,SEG}
        @assert SEG > 0
        size = (M+1)*N*SEG + 1
        @assert SIZE == size

        grid, segment,segindex = [], [],[]
        unilogs = []
        if issparse[1]==false
            for s = 1:SEG
                if s%2==1
                    bound = @SVector[ bounds[(s+1)÷2], (bounds[(s+1)÷2+1]+bounds[(s+1)÷2])/2]
                    if s == SEG && issparse[2]
                        bound = @SVector[ bounds[(s+1)÷2],bounds[(s+1)÷2+1]]
                    end
                    init = 1 + (M+1)*N*(s-1)
                    push!(segment,bound[2])
                    push!(segindex, 1 + (M+1)*N*(s))
                    isopen = @SVector[false, true]
                    if s == SEG && issparse[2]
                        isopen = @SVector[false, false]
                    end
                    g = UniLog{T}(bound,init,minterval,M,N,true,isopen)
                    push!(unilogs,g)
                    for idx=g.idx[1]:g.idx[2]
                        push!(grid, _grid(g, idx))
                    end
                else
                    bound = @SVector[ (bounds[s÷2+1]+bounds[s÷2])/2,bounds[s÷2+1]]
                    init = 1 + (M+1)*N*(s-1)
                    push!(segment, bound[2])
                    push!(segindex, 1 + (M+1)*N*(s))
                    isopen = @SVector[false,s==SEG ? false : true]
                    g = UniLog{T}(bound,init,minterval,M,N,false,isopen)
                    push!(unilogs,g)
                    for idx=g.idx[1]:g.idx[2]
                        push!(grid, _grid(g, idx))
                    end
                end
            end
        else
            for s = 1:SEG
                if s%2==1
                    bound = @SVector[(bounds[(s+1)÷2+1]+bounds[(s+1)÷2])/2,bounds[(s+1)÷2+1]]
                    if s == 1
                        bound = @SVector[ bounds[1], bounds[2]]
                    end
                    init = 1 + (M+1)*N*(s-1)
                    push!(segment,bound[2])
                    push!(segindex, 1 + (M+1)*N*(s))
                    isopen = @SVector[false, s==SEG ? false : true]
                    g = UniLog{T}(bound,init,minterval,M,N,false,isopen)
                    push!(unilogs,g)
                    for idx=g.idx[1]:g.idx[2]
                        push!(grid, _grid(g, idx))
                    end
                else
                    bound = @SVector[bounds[s÷2+1], (bounds[s÷2+1]+bounds[s÷2+2])/2]
                    if s == SEG && issparse[2]
                        bound = @SVector[bounds[s÷2+1], bounds[(s+2)÷2+1]]
                    end
                    init = 1 + (M+1)*N*(s-1)
                    push!(segment, bound[2])
                    push!(segindex, 1 + (M+1)*N*(s))
                    isopen = @SVector[false,s==SEG ? false : true]
                    g = UniLog{T}(bound,init,minterval,M,N,true,isopen)
                    push!(unilogs,g)
                    for idx=g.idx[1]:g.idx[2]
                        push!(grid, _grid(g, idx))
                    end
                end
            end
        end

        head, tail = grid[1], grid[end]
        Isopen[1] && (grid[1] += eps(T) * 1e4)
        Isopen[2] && (grid[end] -= eps(T) * 1e4)
        checkOrder(grid)
        return new{T,size,SEG}(grid, size, head, tail,unilogs, segment,segindex, Isopen)
    end
end

function Base.floor(grid::UniLogs{T,SIZE,SEG}, x) where {T,SIZE,SEG}
#    (grid.head <= x <= grid.tail) || error("$x is out of the uniform grid range!")
    segment = grid.segment
    seg = 0
    for i=1:SEG
        if x<segment[i] && seg==0
            seg = i
        end
    end
    if seg == 0
        seg = SEG
    end
    result = _floor(grid.unilogs[seg],x)
    return Base.floor(Int,result==grid.size ? result-1 : result)
end

Base.getindex(grid::UniLogs, i) = grid.grid[i]
Base.firstindex(grid::UniLogs) = 1
Base.lastindex(grid::UniLogs) = grid.size

"""
    tau(β, halfLife, size::Int, type = Float64)

Create a logarithmic Grid for the imaginary time, which is densest near the 0 and β

#Arguments:
- `β`: inverse temperature
- halfLife: the grid is densest in the range (0, halfLife) and (β-halfLife, β)
- size: the Grid size
"""
@inline function tau(β, halfLife, size::Int, type=Float64)
    size = Int(size)
    c1 = Grid.Coeff{type}([0.0, 0.5β], [1.0, 0.5size + 0.5], 1.0 / halfLife, true)
    r1 = 1:Int(0.5size)

    c2 = Grid.Coeff{type}([0.5β, β], [0.5size + 0.5, size], 1.0 / halfLife, false)
    r2 = (Int(0.5size) + 1):size
    return Log{type,size,2}([c1, c2], [r1, r2], [true, true])
end

"""
    tauUL(β, minterval, M::Int, N::Int, type = Float64)

Create a unilog Grid for the imaginary time, which is densest near the 0 and β
the grid has 2 segments, 2*(M+1)*N+1 points.

#Arguments:
- `β`: inverse temperature
- minterval: the minimal interval of the grid
- M: (M+1) is the number of the 'log' grid of each segment of the grids, M>=0
- N: each interval of the 'log' grid is further divided uniformly into N part, N>=1
"""
@inline function tauUL(β, minterval, M::Int,N::Int, type=Float64)
    seg = 2
    size = (M+1)*N*seg+1
    bounds = @SVector[0.0,β]
    return UniLogs{Float64,size,seg}(bounds,minterval,M,N)
end

@inline function tauUL(para, M::Int,N::Int, type=Float64)
    return tauUL(para.β, para.τ_min, M, N, type)
end
"""
    fermiK(Kf, maxK, halfLife, size::Int, kFi = floor(Int, 0.5size), type = Float64)

Create a logarithmic fermionic K Grid, which is densest near the Fermi momentum ``k_F``

#Arguments:
- Kf: Fermi momentum
- maxK: the upper bound of the grid
- halfLife: the grid is densest in the range (Kf-halfLife, Kf+halfLife)
- size: the Grid size
- kFi: index of Kf
"""
@inline function fermiK(Kf, maxK, halfLife, size::Int, kFi=floor(Int, 0.5size), type=Float64)
    size = Int(size)
    c1 = Grid.Coeff{type}([0.0, Kf], [1.0, kFi], 1.0 / halfLife, false)
    r1 = 1:kFi

    c2 = Grid.Coeff{type}([Kf, maxK], [kFi - 1.0, size], 1.0 / halfLife, true)
    r2 = kFi + 1:size
    return Log{type,size,2}([c1, c2], [r1, r2], [true, false])
end

"""
    fermiKUL(Kf, maxK, minterval, M, N,  type = Float64)

Create a unilog fermionic K Grid, which is densest near the Fermi momentum `k_F`
the grid has 2 segments, 2*(M+1)*N+1 points.

#Arguments:
- Kf: Fermi momentum
- maxK: the upper bound of the grid
- minterval: the minimal interval of the grid
- M: (M+1) is the number of the 'log' grid of each segment of the grids, M>=0
- N: each interval of the 'log' grid is further divided uniformly into N part, N>=1
"""

@inline function fermiKUL(Kf, maxK, minterval, M::Int,N::Int, type=Float64)
    seg = 2
    size = (M+1)*N*seg+1
    bounds = @SVector[0.0,Kf,maxK]
    return UniLogs{Float64,size,seg}(bounds,minterval,M,N,[true,false],[true,true])
end

@inline function fermiKUL(para, M::Int,N::Int, type=Float64)
    return fermiKUL(para.kF, para.k_max, para.k_min, M, N, type)
end

"""
    boseK(Kf, maxK, halfLife, size::Int, kFi = floor(Int, 0.5size), twokFi = floor(Int, 2size / 3), type = Float64)

Create a logarithmic bosonic K Grid, which is densest near the momentum `0` and `2k_F`

#Arguments:
- Kf: Fermi momentum
- maxK: the upper bound of the grid
- halfLife: the grid is densest in the range (0, Kf+halfLife) and (2Kf-halfLife, 2Kf+halfLife)
- size: the Grid size
- kFi: index of Kf
- twokFi: index of 2Kf
"""
@inline function boseK(
    Kf,
    maxK,
    halfLife,
    size::Int,
    kFi=floor(Int, size / 3),
    twokFi=floor(Int, 2size / 3),
    type=Float64,
)
    size = Int(size)
    λ = 1.0 / halfLife
    c1 = Grid.Coeff{type}([0.0, Kf], [1.0, kFi], λ, true)
    r1 = 1:kFi

    c2 = Grid.Coeff{type}([Kf, 2.0 * Kf], [kFi - 1.0, twokFi], λ, false)
    r2 = kFi + 1:twokFi

    c3 = Grid.Coeff{type}([2.0 * Kf, maxK], [twokFi - 1.0, size], λ, true)
    r3 = twokFi + 1:size

    K = Grid.Log{type,size,3}([c1, c2, c3], [r1, r2, r3], [true, false])
return K
end

"""
    boseKUL(Kf, maxK, minterval, M, N,  type = Float64)

Create a unilog bosonic K Grid, which is densest near the momentum `0` and `2k_F`
the grid has 3 segments, 3*(M+1)*N+1 points.

#Arguments:
- Kf: Fermi momentum
- maxK: the upper bound of the grid
- minterval: the minimal interval of the grid
- M: (M+1) is the number of the 'log' grid of each segment of the grids, M>=0
- N: each interval of the 'log' grid is further divided uniformly into N part, N>=1
"""
@inline function boseKUL(
    Kf,
    maxK,
    minterval,
    M,
    N,
    type=Float64
)
    seg = 3
    size = (M+1)*N*seg+1
    bounds = @SVector[0.0,2Kf,maxK]
    return UniLogs{Float64,size,seg}(bounds,minterval,M,N,[true,false],[false,true])
end

@inline function boseKUL(para, M::Int,N::Int, type=Float64)
    return boseKUL(para.kF, para.k_max, para.k_min, M, N, type)
end

include("interpolate.jl")

end


text:
# ---
# title: 605. Can Place Flowers
# id: problem605
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Greedy
# link: <https://leetcode.com/problems/can-place-flowers/description/>
# hidden: true
# ---
#
# You have a long flowerbed in which some of the plots are planted, and some are
# not. However, flowers cannot be planted in **adjacent** plots.
#
# Given an integer array `flowerbed` containing `0`'s and `1`'s, where `0` means
# empty and `1` means not empty, and an integer `n`, return _if_ `n` new flowers
# can be planted in the `flowerbed` without violating the no-adjacent-flowers
# rule.
#
#
#
# **Example 1:**
#
#
#
#     Input: flowerbed = [1,0,0,0,1], n = 1
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: flowerbed = [1,0,0,0,1], n = 2
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= flowerbed.length <= 2 * 104`
#   * `flowerbed[i]` is `0` or `1`.
#   * There are no two adjacent flowers in `flowerbed`.
#   * `0 <= n <= flowerbed.length`
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
# Fill in the gradients . No defaults since almost always called internally to custom rule
function generate_perturbation_derivatives!(m, p_d, p_f, cache::AbstractSolverCache{1};
                                            settings = PerturbationSolverSettings())
    @assert cache.p_d_symbols == collect(Symbol.(keys(p_d)))

    p = isnothing(p_f) ? p_d : order_vector_by_symbols(merge(p_d, p_f), m.mod.p_symbols)

    # Fill in derivatives
    ret = evaluate_first_order_functions_p!(m, cache, settings, p)
    (ret == :Success) || return ret
    ret = solve_first_order_p!(m, cache, settings)
    maybe_call_function(settings.solve_first_order_p_callback, ret, m, cache, settings)
    return ret
end

function generate_perturbation_derivatives!(m, p_d, p_f, cache::AbstractSolverCache{2};
                                            settings = PerturbationSolverSettings())
    @assert cache.p_d_symbols == collect(Symbol.(keys(p_d)))

    p = isnothing(p_f) ? p_d : order_vector_by_symbols(merge(p_d, p_f), m.mod.p_symbols)

    # Fill in derivatives, first by calling the first-order
    # Fill in derivatives
    ret = evaluate_first_order_functions_p!(m, cache, settings, p)
    (ret == :Success) || return ret
    ret = solve_first_order_p!(m, cache, settings)
    maybe_call_function(settings.solve_first_order_p_callback, ret, m, cache, settings)

    # 2nd Order calculations
    ret = evaluate_second_order_functions_p!(m, cache, settings, p)
    (ret == :Success) || return ret
    ret = solve_second_order_p!(m, cache, settings)
    maybe_call_function(settings.solve_second_order_p_callback, ret, m, cache, settings)
    return ret
end

function evaluate_first_order_functions_p!(m, c, settings, p)
    (settings.print_level > 2) &&
        println("Evaluating first-order function derivatives into cache")
    try
        @unpack y, x = c  # Precondition: valid (y, x) steady states
        isnothing(c.H_p) ||
            fill_array_by_symbol_dispatch(m.mod.H_p!, c.H_p, c.p_d_symbols, y, x, p) #not required if steady_state_p!
        fill_array_by_symbol_dispatch(m.mod.H_yp_p!, c.H_yp_p, c.p_d_symbols, y, x, p)
        fill_array_by_symbol_dispatch(m.mod.H_y_p!, c.H_y_p, c.p_d_symbols, y, x, p)
        fill_array_by_symbol_dispatch(m.mod.H_xp_p!, c.H_xp_p, c.p_d_symbols, y, x, p)
        fill_array_by_symbol_dispatch(m.mod.H_x_p!, c.H_x_p, c.p_d_symbols, y, x, p)
        fill_array_by_symbol_dispatch(m.mod.Γ_p!, c.Γ_p, c.p_d_symbols, p)
        isnothing(c.Ω_p) ||
            fill_array_by_symbol_dispatch(m.mod.Ω_p!, c.Ω_p, c.p_d_symbols, p)
    catch e
        if !is_linear_algebra_exception(e)
            (settings.print_level > 2) && println("Rethrowing exception")
            rethrow(e)
        else
            settings.print_level == 0 || display(e)
            return :Failure # generic failure
        end
    end
    return :Success  # no failing code-paths yet.
end

function evaluate_second_order_functions_p!(m, c, settings, p)
    (settings.print_level > 2) &&
        println("Evaluating second-order function derivatives into cache")
    try
        @unpack y, x = c  # Precondition: valid (y, x) steady states
        fill_array_by_symbol_dispatch(m.mod.Ψ_p!, c.Ψ_p, c.p_d_symbols, y, x, p)
    catch e
        if !is_linear_algebra_exception(e)
            (settings.print_level > 2) && println("Rethrowing exception")
            rethrow(e)
        else
            settings.print_level == 0 || display(e)
            return :Failure # generic failure
        end
    end
    return :Success  # no failing code-paths yet.
end

function solve_first_order_p!(m, c, settings)
    @unpack n_x, n_y, n_ϵ = m
    n = n_x + n_y
    n_p = length(c.p_d_symbols)
    (settings.print_level > 2) && println("Solving first order derivatives of perturbation")

    buff = c.first_order_solver_p_buffer
    try
        if isnothing(m.mod.ȳ_p!) && isnothing(m.mod.x̄_p!)
            # Zeroth-order derivatives if not provided
            # Calculating c.y_p, c.x_p
            A_zero = [c.H_y + c.H_yp c.H_x + c.H_xp]
            # H_p is a vector of vectors
            x_zeroth = A_zero \ -hcat(c.H_p...) # (47)
            c.y_p .= x_zeroth[1:n_y, :]
            c.x_p .= x_zeroth[(n_y + 1):n, :]
        end

        # The derivatives
        for i in 1:n_p
            c.Σ_p[i] .= Symmetric(c.Γ_p[i] * c.Γ' + c.Γ * c.Γ_p[i]')
        end

        # Write equation (52) as E + AX + CXD = 0, a generalized Sylvester equation
        # first-order derivatives
        buff.R .= vcat(c.g_x * c.h_x, c.g_x, c.h_x, c.I_x)
        buff.A .= [c.H_y c.H_xp + c.H_yp * c.g_x]

        # i.e. C = [c.H_yp zeros(n, n_x)]
        buff.C[:, 1:n_y] .= c.H_yp
        ws = GeneralizedSylvesterWs(n, n, n_x, 1)

        # Initialize
        for i in 1:n_p
            # p-specific Sylvester preparation
            buff.bar[1:n_y] = c.y_p[i]
            buff.bar[(n_y + 1):(2 * n_y)] = c.y_p[i]
            buff.bar[(2 * n_y + 1):(2 * n_y + n_x)] = c.x_p[i]
            buff.bar[(2 * n_y + n_x + 1):end] = c.x_p[i]

            buff.dH[:, 1:n_y] = c.H_yp_p[i]
            buff.dH[:, (n_y + 1):(2 * n_y)] = c.H_y_p[i]
            buff.dH[:, (2 * n_y + 1):(2 * n_y + n_x)] = c.H_xp_p[i]
            buff.dH[:, (2 * n_y + n_x + 1):end] = c.H_x_p[i]

            for j in 1:n
                buff.dH[j, :] += c.Ψ[j] * buff.bar
            end
            mul!(buff.E, buff.dH, buff.R, -1.0, 0.0)

            # solves AXB + CXD = E
            # Sylvester
            generalized_sylvester_solver!(buff.A, buff.C, c.h_x, buff.E, 1, ws)
            c.g_x_p[i] .= buff.E[1:n_y, :]
            c.h_x_p[i] .= buff.E[(n_y + 1):n, :]

            # Q weighted derivatives
            c.C_1_p[i] .= c.Q * vcat(c.g_x_p[i], zeros(n_x, n_x))
            c.A_1_p[i] .= c.h_x_p[i]

            # V derivatives
            tmp = c.h_x_p[i] * Array(c.V) * c.h_x'
            c.V_p[i] .= lyapd(c.h_x, c.η * c.Σ_p[i] * c.η' + tmp + tmp')
            # B derivatives
            c.B_p[i] .= c.η * c.Γ_p[i]
        end

        #@exfiltrate  # flip on to see intermediate calculations.  TURN OFF BEFORE PROFILING

    catch e
        if !is_linear_algebra_exception(e)
            (settings.print_level > 2) && println("Rethrowing exception")
            rethrow(e)
        else
            settings.print_level == 0 || display(e)
            return :Failure # generic failure
        end
    end
    return :Success
end

function solve_second_order_p!(m, c, settings)
    @unpack n_x, n_y, n_ϵ, n_z = m
    n = n_x + n_y
    n_p = length(c.p_d_symbols)
    (settings.print_level > 2) &&
        println("Solving second order derivatives of perturbation")
    buff = c.second_order_solver_p_buffer

    try
        # General Prep
        buff.A .= [c.H_y c.H_xp + c.H_yp * c.g_x]
        buff.C[:, 1:n_y] .= c.H_yp
        ws = GeneralizedSylvesterWs(n, n, n_x, 2)
        buff.E .= zeros(n, n_x * n_x)
        buff.R .= vcat(c.g_x * c.h_x, c.g_x, c.h_x, c.I_x)
        buff.A_σ .= [c.H_yp + c.H_y c.H_xp + c.H_yp * c.g_x]
        buff.R_σ .= vcat(c.g_x, zeros(n_y, n_x), c.I_x, zeros(n_x, n_x))

        buff.gh_stack[1:n_y, :] .= reshape(c.g_xx, n_y, n_x * n_x)
        buff.gh_stack[n_y+1:end, :] .= reshape(c.h_xx, n_x, n_x * n_x)
        buff.g_xx_flat .= reshape(c.g_xx, n_y, n_x * n_x)
        for i in 1:n
            for j in 1:n_x
                buff.Ψ_x_sum[j][i] .= c.Ψ_xp[j][i] .+ c.Ψ_x[j][i]
            end
            for j in 1:n_y
                buff.Ψ_y_sum[j][i] .= c.Ψ_yp[j][i] .+ c.Ψ_y[j][i]
            end
        end

        tmp1 = similar(buff.R')
        tmp2 = similar(c.h_x)
        for i in 1:n_p
            # Prep for _xx_p
            # Compute the total derivatives
            buff.bar[1:n_y] = c.y_p[i]
            buff.bar[(n_y + 1):(2 * n_y)] = c.y_p[i]
            buff.bar[(2 * n_y + 1):(2 * n_y + n_x)] = c.x_p[i]
            buff.bar[(2 * n_y + n_x + 1):end] = c.x_p[i]
            buff.dH[:, 1:n_y] = c.H_yp_p[i]
            buff.dH[:, (n_y + 1):(2 * n_y)] = c.H_y_p[i]
            buff.dH[:, (2 * n_y + 1):(2 * n_y + n_x)] = c.H_xp_p[i]
            buff.dH[:, (2 * n_y + n_x + 1):end] = c.H_x_p[i]
            for j in 1:n
                buff.dH[j, :] += c.Ψ[j] * buff.bar
            end
            for j in 1:n
                buff.dΨ[j] .= c.Ψ_p[i][j]
            end
            for j in 1:n_y
                if (c.y_p[i][j] != 0)
                    for k in 1:n
                        buff.dΨ[k] .+= buff.Ψ_y_sum[j][k] .* c.y_p[i][j]
                    end
                end
            end
            for j in 1:n_x
                if (c.x_p[i][j] != 0)
                    for k in 1:n
                        buff.dΨ[k] .+= buff.Ψ_x_sum[j][k] .* c.x_p[i][j]
                    end
                end
            end

            # Constants: (60)
            buff.R_p[1:n_y, :] = c.g_x_p[i] * c.h_x + c.g_x * c.h_x_p[i]
            buff.R_p[(n_y + 1):(2 * n_y), :] = c.g_x_p[i]
            buff.R_p[(2 * n_y + 1):(2 * n_y + n_x), :] = c.h_x_p[i]
            # Flip the sign of E for Sylvester input
            fill!(buff.E, 0.0)
            for j in 1:n
                mul!(tmp1, buff.R_p', c.Ψ[j])
                mul!(tmp2, tmp1, buff.R)
                buff.E[j, :] .-= vec(tmp2)
                buff.E[j, :] .-= vec(tmp2')
                mul!(tmp1, buff.R', buff.dΨ[j])
                mul!(tmp2, tmp1, buff.R)
                buff.E[j, :] .-= vec(tmp2)
                # buff.E[j, :] .= -vec(buff.R_p' * c.Ψ[j] * buff.R + buff.R' * c.Ψ[j] * buff.R_p + buff.R' * buff.dΨ[j] * buff.R)
            end
            # Constants: (56)
            kron!(buff.kron_h_x, c.h_x, c.h_x)
            buff.E .-= buff.dH[:, 1:n_y] * buff.g_xx_flat * buff.kron_h_x # Plug (57) in (56)
            tmp = c.H_yp * buff.g_xx_flat
            kron!(buff.kron_h_x, c.h_x_p[i], c.h_x)
            mul!(buff.E, tmp, buff.kron_h_x, -1.0, 1.0) # Plug (58) in (56), step 1
            kron!(buff.kron_h_x, c.h_x, c.h_x_p[i])
            mul!(buff.E, tmp, buff.kron_h_x, -1.0, 1.0) # Plug (58) in (56), step 2
            buff.E .-= hcat(buff.dH[:, (n_y + 1):(2 * n_y)],
            buff.dH[:, 1:n_y] * c.g_x +
            c.H_yp * c.g_x_p[i] +
            buff.dH[:, (2 * n_y + 1):(2 * n_y + n_x)]) * buff.gh_stack # Plug (59) in (56)

            # Solve the Sylvester equations (56)
            generalized_sylvester_solver!(buff.A, buff.C, c.h_x, buff.E, 2, ws)
            copyto!(c.g_xx_p[i], buff.E[1:n_y, :]) # Reshaping into n_y * n_x * n_x
            copyto!(c.h_xx_p[i], buff.E[(n_y + 1):end, :]) # # Reshaping into n_x * n_x * n_x

            # Prep for _σσ_p
            # Solve _σσ_p
            R_σ_p = vcat(c.g_x_p[i], zeros(n_y + n_x * 2, n_x))
            η_sq_p = c.η * c.Σ_p[i] * c.η'
            C_σ = -(buff.dH[:, 1:n_y] + buff.dH[:, (n_y + 1):(2 * n_y)]) * c.g_σσ
            C_σ -= (buff.dH[:, 1:n_y] * c.g_x + c.H_yp * c.g_x_p[i] + buff.dH[:, (2 * n_y + 1):(2 * n_y + n_x)]) * c.h_σσ
            # C_σ = -hcat(buff.dH[:, 1:n_y] + buff.dH[:, (n_y + 1):(2 * n_y)],
            #             buff.dH[:, 1:n_y] * c.g_x +
            #             c.H_yp * c.g_x_p[i] +
            #             buff.dH[:, (2 * n_y + 1):(2 * n_y + n_x)]) * vcat(c.g_σσ, c.h_σσ) # Plug (65) in (64), flip the sign to solve (64)
            C_σ -= (buff.dH[:, 1:n_y] * buff.g_xx_flat + c.H_yp * buff.E[1:n_y, :]) * vec(c.η_Σ_sq)# (67), 2nd line
            C_σ -= (c.H_yp * buff.g_xx_flat) * vec(η_sq_p) # (67), 3rd line, second part
            for j in 1:n
                mul!(tmp1, R_σ_p', c.Ψ[j])
                mul!(tmp2, tmp1, buff.R_σ)
                C_σ[j] -= dot(tmp2, c.η_Σ_sq)
                C_σ[j] -= dot(tmp2', c.η_Σ_sq)
                mul!(tmp1, buff.R_σ', buff.dΨ[j])
                mul!(tmp2, tmp1, buff.R_σ)
                C_σ[j] -= dot(tmp2, c.η_Σ_sq)
                # C_σ[j] -= dot((R_σ_p' * c.Ψ[j] * buff.R_σ + buff.R_σ' * c.Ψ[j] * R_σ_p + buff.R_σ' * buff.dΨ[j] * buff.R_σ), c.η_Σ_sq) # (67), 1st line
                mul!(tmp1, buff.R_σ', c.Ψ[j])
                mul!(tmp2, tmp1, buff.R_σ)
                C_σ[j] -= dot(tmp2, η_sq_p)
                # C_σ[j] -= dot((buff.R_σ' * c.Ψ[j] * buff.R_σ), η_sq_p) # (67), 3rd line, first part
            end

            # Solve _σσ_p
            X_σ = buff.A_σ \ C_σ # solve (64)
            c.g_σσ_p[:, i] .= X_σ[1:n_y]
            c.h_σσ_p[:, i] .= X_σ[(n_y + 1):end]
            fill!(c.C_2_p[i], 0.0) # reset as we need to use `+=`
            for j in 1:n_z
                for k in 1:n_y
                    if (c.Q[j, k] != 0)
                        c.C_2_p[i][j, :, :] .+= 0.5 * c.Q[j, k] * c.g_xx_p[i][k, :, :]
                    end
                end
            end
            c.A_2_p[i] .= 0.5 * c.h_xx_p[i]
        end

        c.C_0_p .= 0.5 * c.Q * vcat(c.g_σσ_p, zeros(n_x, n_p))
        c.A_0_p .= 0.5 * c.h_σσ_p

    catch e
        settings.print_level == 0 || display(e)
        return :Failure # generic failure
    end
    return :Success
end

function ChainRulesCore.rrule(::typeof(generate_perturbation), m::PerturbationModel,
                              p_d::NamedTuple{DFieldsType,DTupleType}, p_f, order::Val{1};
                              cache = SolverCache(m, Val(1), p_d),
                              settings = PerturbationSolverSettings()) where {DFieldsType,
                                                                              DTupleType}
    (settings.print_level > 2) && println("Calculating generate_perturbation primal ")
    sol = generate_perturbation(m, p_d, p_f, Val(1); cache, settings)
    if (sol.retcode == :Success)
        grad_ret = generate_perturbation_derivatives!(m, p_d, p_f, cache)
        if (grad_ret != :Success)
            sol = FirstOrderPerturbationSolution(:GradientFailure, m, cache)
        end
    end

    function generate_perturbation_pb(Δsol)
        (settings.print_level > 2) && println("Calculating generate_perturbation pullback")
        Δp = (p_d === nothing) ? nothing : zeros(length(p_d))
        if (sol.retcode == :Success) & (p_d !== nothing)
            c = cache # temp to avoid renaming everything
            n_p_d = length(p_d)
            if (~iszero(Δsol.A))
                for i in 1:n_p_d
                    Δp[i] += dot(c.h_x_p[i], Δsol.A)
                end
            end
            if (~iszero(Δsol.g_x))
                for i in 1:n_p_d
                    Δp[i] += dot(c.g_x_p[i], Δsol.g_x)
                end
            end
            if (~iszero(Δsol.C))
                for i in 1:n_p_d
                    Δp[i] += dot(c.C_1_p[i], Δsol.C)
                end
            end
            if (~isnothing(Δsol.x_ergodic))
                if ((Δsol.x_ergodic != NoTangent()) & (Δsol.x_ergodic != ZeroTangent()))
                    for i in 1:n_p_d
                        tmp = c.V.L \ c.V_p[i] / c.V.U
                        tmp[diagind(tmp)] /= 2.0
                        t1 = c.V.L * LowerTriangular(tmp)
                        # Cholesky by default stores the U part, but that information is lost when passing back
                        # from MvNormal. Therefore, we have to extract the components out manually
                        Δp[i] += dot(t1', UpperTriangular(Δsol.x_ergodic.C.factors)) # dot(c.V_p[i], Δsol.x_ergodic) is the original logic
                    end
                end
            end
            if (~iszero(Δsol.Γ))
                for i in 1:n_p_d
                    Δp[i] += dot(c.Γ_p[i], Δsol.Γ)
                end
            end
            if (~iszero(Δsol.B))
                for i in 1:n_p_d
                    Δp[i] += dot(c.B_p[i], Δsol.B)
                end
            end
            if (~isnothing(Δsol.D)) # D is a Distribution object which complicates stuff here
                if ((Δsol.D != NoTangent()) & (Δsol.D != ZeroTangent()))
                    # Only supports diagonal matrices for now.
                    for i in 1:n_p_d
                        Δp[i] += dot(c.Ω_p[i], Δsol.D.σ)
                    end
                end
            end
            if (~iszero(Δsol.x))
                for i in 1:n_p_d
                    Δp[i] += dot(c.x_p[i], Δsol.x)
                end
            end
            if (~iszero(Δsol.y))
                for i in 1:n_p_d
                    Δp[i] += dot(c.y_p[i], Δsol.y)
                end
            end
        end
        Δp_nt = NamedTuple{DFieldsType,DTupleType}(tuple(Δp...))  # turn tuple into named tuple in the same order
        return NoTangent(), NoTangent(),
               Tangent{NamedTuple{DFieldsType,DTupleType},
                       NamedTuple{DFieldsType,DTupleType}}(Δp_nt), NoTangent(), NoTangent()
    end
    # keep the named tuple the same
    return sol, generate_perturbation_pb
end

function ChainRulesCore.rrule(::typeof(generate_perturbation), m::PerturbationModel,
                              p_d::NamedTuple{DFieldsType,DTupleType}, p_f, order::Val{2};
                              cache = SolverCache(m, Val(2), p_d),
                              settings = PerturbationSolverSettings()) where {DFieldsType,
                                                                              DTupleType}
    (settings.print_level > 2) && println("Calculating generate_perturbation primal ")
    sol = generate_perturbation(m, p_d, p_f, Val(2); cache, settings)
    if (sol.retcode == :Success)
        grad_ret = generate_perturbation_derivatives!(m, p_d, p_f, cache)
        if (grad_ret != :Success)
            sol = SecondOrderPerturbationSolution(:GradientFailure, m, cache)
        end
    end
    function generate_perturbation_pb(Δsol)
        (settings.print_level > 2) && println("Calculating generate_perturbation pullback")
        Δp = (p_d === nothing) ? nothing : zeros(length(p_d))
        if (sol.retcode == :Success) & (p_d !== nothing)
            c = cache # temp to avoid renaming everything
            n_p_d = length(p_d)
            if (~iszero(Δsol.A_1))
                for i in 1:n_p_d
                    Δp[i] += dot(c.A_1_p[i], Δsol.A_1)
                end
            end
            if (~iszero(Δsol.g_x))
                for i in 1:n_p_d
                    Δp[i] += dot(c.g_x_p[i], Δsol.g_x)
                end
            end
            if (~iszero(Δsol.C_1))
                for i in 1:n_p_d
                    Δp[i] += dot(c.C_1_p[i], Δsol.C_1)
                end
            end
            if (~iszero(Δsol.Γ))
                for i in 1:n_p_d
                    Δp[i] += dot(c.Γ_p[i], Δsol.Γ)
                end
            end
            if (~iszero(Δsol.B))
                for i in 1:n_p_d
                    Δp[i] += dot(c.B_p[i], Δsol.B)
                end
            end
            if (~iszero(Δsol.A_2))
                for i in 1:n_p_d
                    Δp[i] += dot(c.A_2_p[i], Δsol.A_2)
                end
            end
            if (~iszero(Δsol.g_xx))
                for i in 1:n_p_d
                    Δp[i] += dot(c.g_xx_p[i], Δsol.g_xx)
                end
            end
            if (~iszero(Δsol.C_2))
                for i in 1:n_p_d
                    Δp[i] += dot(c.C_2_p[i], Δsol.C_2)
                end
            end
            if (~isnothing(Δsol.D)) # D is a Distribution object which complicates stuff here
                if ((Δsol.D != NoTangent()) & (Δsol.D != ZeroTangent()))
                    # Only supports diagonal matrices for now.
                    for i in 1:n_p_d
                        Δp[i] += dot(c.Ω_p[i], Δsol.D.σ)
                    end
                end
            end
            if (~iszero(Δsol.x))
                for i in 1:n_p_d
                    Δp[i] += dot(c.x_p[i], Δsol.x)
                end
            end
            if (~iszero(Δsol.y))
                for i in 1:n_p_d
                    Δp[i] += dot(c.y_p[i], Δsol.y)
                end
            end
            # Currently the results for the second-order correction terms are not vector of vectors, but just arrays
            if (~iszero(Δsol.C_0))
                Δp += c.C_0_p' * Δsol.C_0
            end
            if (~iszero(Δsol.g_σσ))
                Δp += c.g_σσ_p' * Δsol.g_σσ
            end
            if (~iszero(Δsol.A_0))
                Δp += c.A_0_p' * Δsol.A_0
            end
        end

        Δp_nt = NamedTuple{DFieldsType,DTupleType}(tuple(Δp...)) # turn tuple into named tuple in the same order
        return NoTangent(), NoTangent(),
               Tangent{NamedTuple{DFieldsType,DTupleType},
                       NamedTuple{DFieldsType,DTupleType}}(Δp_nt), NoTangent(), NoTangent()
    end

    return sol, generate_perturbation_pb
end


text:
using BenchmarkTools
using SolidStateDetectors
using Plots
using Unitful

using SolidStateDetectors: rb_even, rb_odd

### Inner loop of SOR
begin
    T = Float32
    sim = Simulation{T}(SSD_examples[:CGD]);

    calculate_electric_potential!( sim, depletion_handling = false,
                                        convergence_limit = 1e-4,
                                        max_tick_distance = 1u"mm",
                                        missing,
                                        # refinement_limits = [0.2, 0.1, 0.05, 0.02],
                                        refinement_limits = missing,
                                        )

    # plot(
    #     plot(sim.point_types, y = 0),
    #     plot(sim.electric_potential, y = 0),
    # )

    pssrb = SolidStateDetectors.PotentialSimulationSetupRB( sim.detector,
                                                            sim.electric_potential.grid,
                                                            # Grid(sim, max_tick_distance = 1u"mm"),
                                                            SolidStateDetectors.material_properties[SolidStateDetectors.materials["vacuum"]],
                                                            sim.electric_potential.data);

    nthreads = 1 #Base.Threads.nthreads()
    is_weighting_potential = Val{false}()
    depletion_handling = Val{false}()
    only2d = Val{false}()
    update_even_points = Val{true}()
    update_uneven_points = Val{false}()
end

begin
    @info "Grid size: $(size(sim.electric_potential.data))"
    # SolidStateDetectors.outerloop!(pssrb, nthreads, update_even_points, depletion_handling, is_weighting_potential, only2d)
    for dp in (false, true)
        depletion_handling = Val{dp}()
        @info "Depletion handling: $dp"
        for nt in (1, 2, 4, 8, 16, 32, 64)
            @info "N Threads: $nt"
            @btime SolidStateDetectors.outerloop!($pssrb, $nt, $update_even_points, $depletion_handling, $is_weighting_potential, $only2d)
        end
    end
end

# middleloop!
begin
    idx3 = 2
    even_points = true
    update_even_points = Val(even_points)
    is_weighting_potential = Val(false)
    rb_tar_idx, rb_src_idx = even_points ? (rb_even, rb_odd) : (rb_odd, rb_even)
    depletion_handling = Val{false}()
    idx3iseven = Val(iseven(idx3))

    SolidStateDetectors.middleloop!(idx3, rb_tar_idx, rb_src_idx, pssrb,
                         update_even_points, depletion_handling, is_weighting_potential, only2d, idx3iseven)

    @code_warntype SolidStateDetectors.middleloop!(idx3, rb_tar_idx, rb_src_idx, pssrb,
                         update_even_points, depletion_handling, is_weighting_potential, only2d, idx3iseven)

    @btime SolidStateDetectors.middleloop!($idx3, $rb_tar_idx, $rb_src_idx, $pssrb,
                         $update_even_points, $depletion_handling, $is_weighting_potential, $only2d, $idx3iseven)
end

# load_weights_for_innerloop! & innerloop!
begin
    depletion_handling = Val{false}()

    idx3_is_even = Val(iseven(idx3))

    iz = idx3
    inz = iz - 1

    pwwzr        = pssrb.geom_weights[3].weights[1, inz]
    pwwzl        = pssrb.geom_weights[3].weights[2, inz]
    pwΔmpz       = pssrb.geom_weights[3].weights[3, inz]
    Δz_ext_inv_r = pssrb.geom_weights[3].weights[4, inz + 1]
    Δz_ext_inv_l = pssrb.geom_weights[3].weights[4, inz]

    line_weights = Array{T, 2}(undef, size(pssrb.potential, 1) - 2, 6)

    iy = 2
    iny = iy - 1

    pwwyr  = pssrb.geom_weights[2].weights[1, iny]
    pwwyl  = pssrb.geom_weights[2].weights[2, iny]
    pwΔmpy = pssrb.geom_weights[2].weights[3, iny]
    pwΔmpy_pwΔmpz = pwΔmpy * pwΔmpz
    Δy_ext_inv_r_pwΔmpz  = pssrb.geom_weights[2].weights[4, iny + 1] * pwΔmpz
    Δy_ext_inv_l_pwΔmpz  = pssrb.geom_weights[2].weights[4, iny]     * pwΔmpz
    Δz_ext_inv_r_pwΔmpy = Δz_ext_inv_r * pwΔmpy
    Δz_ext_inv_l_pwΔmpy = Δz_ext_inv_l * pwΔmpy

    pwwyr_pwwzr = pwwyr * pwwzr
    pwwyr_pwwzl = pwwyr * pwwzl
    pwwyl_pwwzr = pwwyl * pwwzr
    pwwyl_pwwzl = pwwyl * pwwzl

    SolidStateDetectors.load_weights_for_innerloop!(line_weights, pssrb, iy, iny, iz, inz,
                update_even_points, idx3_is_even,
                pwwzr, pwwzl, pwwyr, pwwyl,
                pwwyr_pwwzr, pwwyr_pwwzl, pwwyl_pwwzr, pwwyl_pwwzl,
                pwΔmpy_pwΔmpz,
                Δy_ext_inv_r_pwΔmpz, Δy_ext_inv_l_pwΔmpz, Δz_ext_inv_r_pwΔmpy, Δz_ext_inv_l_pwΔmpy)

    @btime SolidStateDetectors.load_weights_for_innerloop!($line_weights, $pssrb, $iy, $iny, $iz, $inz,
                $update_even_points, $idx3_is_even,
                $pwwzr, $pwwzl, $pwwyr, $pwwyl,
                $pwwyr_pwwzr, $pwwyr_pwwzl, $pwwyl_pwwzr, $pwwyl_pwwzl,
                $pwΔmpy_pwΔmpz,
                $Δy_ext_inv_r_pwΔmpz, $Δy_ext_inv_l_pwΔmpz, $Δz_ext_inv_r_pwΔmpy, $Δz_ext_inv_l_pwΔmpy)

    SolidStateDetectors.innerloop!(line_weights, pssrb, iy, iny, iz, inz, rb_tar_idx, rb_src_idx,
                update_even_points, idx3_is_even,
                depletion_handling, is_weighting_potential, only2d)

    @code_warntype SolidStateDetectors.innerloop!(line_weights, pssrb, iy, iny, iz, inz, rb_tar_idx, rb_src_idx,
                update_even_points, idx3_is_even,
                depletion_handling, is_weighting_potential, only2d)

    @btime SolidStateDetectors.innerloop!($line_weights, $pssrb, $iy, $iny, $iz, $inz, $rb_tar_idx, $rb_src_idx,
                $update_even_points, $idx3_is_even,
                $depletion_handling, $is_weighting_potential, $only2d)

    @code_llvm SolidStateDetectors.innerloop!(line_weights, pssrb, iy, iny, iz, inz, rb_tar_idx, rb_src_idx,
                update_even_points, idx3_is_even,
                depletion_handling, is_weighting_potential, only2d)
end


text:
"""
    pcr(X, Y, weights = ones(size(X, 1)); nlv)
Principal component regression (PCR) with a SVD factorization.
* `X` : X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.
* `nlv` : Nb. latent variables (LVs) to compute.

`X` and `Y` are internally centered.
The model is computed with an intercept.

## Examples

```julia
n = 10 ; p = 5 ; q = 2
X = rand(n, p)
Y = rand(n, q)

nlv = 3
fm = pcr(X, Y, w; nlv = nlv) ;
Jchemo.coef(fm).B
Jchemo.coef(fm).int
Jchemo.coef(pcr(X, Y[:, 2], w; nlv = nlv)).B
Jchemo.coef(fm; nlv = 1).B

nlv = 3
fm = pcr(X, Y, w; nlv = nlv) ;
Jchemo.predict(fm, X[1:3, :]).pred
Jchemo.predict(plskern(X, Y, w; nlv = nlv), X[1:3, :]).pred

zfm = fm.fm_pca ;
res = summary(zfm, X) ;
pnames(res)
res.explvar

fm.T
Jchemo.transform(zfm, X[1:3, :])
```
"""
function pcr(X, Y, weights = ones(size(X, 1)); nlv)
    pcr!(copy(X), copy(Y), weights; nlv = nlv)
end

function pcr!(X, Y, weights = ones(size(X, 1)); nlv)
    Y = ensure_mat(Y)
    weights = mweights(weights)
    ymeans = colmeans(Y, weights)
    center!(Y, ymeans)
    fm = pcasvd!(X, weights; nlv = nlv)
    D = Diagonal(fm.weights)
    beta = inv(fm.T' * D * fm.T) * fm.T' * D * Y
    # first term = Diagonal(1 ./ fm.sv[1:nlv].^2) if T is D-orthogonal
    # This is the case for the actual version (pcasvd)
    Pcr(fm, fm.T, fm.P, beta', fm.xmeans, ymeans, weights)
end


text:
using Pkg
Pkg.activate(".")
using LTMSim
using DataFrames
using SimpleHypergraphs
using Statistics
using Plots
using PyPlot
using Random
using Serialization
using LaTeXStrings

h = hg_load("/Users/carminespagnuolo/Dropbox/LTMSim.jl/data/got.hgf")

# he -> degree
degrees = Dict{Int, Int}()

for he=1:nhe(h)
    push!(
        degrees,
        he => length(getvertices(h, he))
    )
end

sorted_degrees = sort(collect(degrees), by=x->x[2], rev=true)

y = map(x->x[2], sorted_degrees)

##############
# distribution fitting
##############
f = fit()

clf()
plt.figure(figsize=(7,3.5))

plt.scatter(range(1, stop=nhe(h)), y, s=.2)

ylabel("Hyperedge size", fontsize="xx-large", fontstyle="italic", labelpad=10) #, fontweight="demibold"
xlabel(L"$n$", fontsize="xx-large", fontweight="bold", labelpad=10, position=(0,100))

plt.xticks(fontsize="x-large") #range(0, length(ticks) * 2, step=2),
plt.yticks(fontsize="x-large")

plt.tight_layout()

gcf()

PyPlot.savefig("res/paper/got/distribution.png")


text:
using JLD2
using DataFrames
using IcanProj

"""
Summaries results from vic into single tables.
"""
function summaries_results(path_sim, path_eval, watersheds, resolutions)

    # Load watershed names

    # @load Pkg.dir("IcanProj", "data", "wsh_info.jld2") wsh_info

    wsh_name = get_wsh_name(wsh_info)

    for watershed in watersheds

        for resolution in resolutions

            @show resolution

            # Read results in simulation directory

            df_vic, res_vic = read_all_results(path_sim, watershed, resolution, wsh_name)

            # Save results to evaluation directory

            path = joinpath(path_eval, watershed)

            mkpath(path)

            writetable(joinpath(path, "res_$(resolution).csv"), df_vic)

            @save joinpath(path, "res_$(resolution).jld2") res_vic

        end

    end

end

# Process data

opt = get_options()

path_sim = opt["target_folder"]

path_eval = opt["eval_folder"]

watersheds = opt["stat_sel"]

resolutions = opt["resolutions"]

summaries_results(path_sim, path_eval, watersheds, resolutions)


text:
__precompile__()
module COMBO
export combo

using PyCall

const combo = PyCall.PyNULL()

function __init__()
    copy!(combo, pyimport("combo"))
end

end # module


text:
"""
    conv_bn(kernelsize, inplanes, outplanes, activation = relu;
            rev = false,
            stride = 1, pad = 0, dilation = 1, groups = 1, [bias, weight, init],
            initβ = Flux.zeros32, initγ = Flux.ones32, ϵ = 1f-5, momentum = 1f-1)

Create a convolution + batch normalization pair with ReLU activation.

# Arguments
- `kernelsize`: size of the convolution kernel (tuple)
- `inplanes`: number of input feature maps
- `outplanes`: number of output feature maps
- `activation`: the activation function for the final layer
- `rev`: set to `true` to place the batch norm before the convolution
- `stride`: stride of the convolution kernel
- `pad`: padding of the convolution kernel
- `dilation`: dilation of the convolution kernel
- `groups`: groups for the convolution kernel
- `bias`, `weight`, `init`: initialization for the convolution kernel (see [`Flux.Conv`](#))
- `initβ`, `initγ`: initialization for the batch norm (see [`Flux.BatchNorm`](#))
- `ϵ`, `momentum`: batch norm parameters (see [`Flux.BatchNorm`](#))
"""
function conv_bn(kernelsize, inplanes, outplanes, activation = relu;
                 rev = false,
                 initβ = Flux.zeros32, initγ = Flux.ones32, ϵ = 1f-5, momentum = 1f-1,
                 kwargs...)
  layers = []

  if rev
    activations = (conv = activation, bn = identity)
    bnplanes = inplanes
  else
    activations = (conv = identity, bn = activation)
    bnplanes = outplanes
  end

  push!(layers, Conv(kernelsize, Int(inplanes) => Int(outplanes), activations.conv; kwargs...))
  push!(layers, BatchNorm(Int(bnplanes), activations.bn;
                          initβ = initβ, initγ = initγ, ϵ = ϵ, momentum = momentum))

  return rev ? reverse(layers) : layers
end

"""
    cat_channels(x, y)

Concatenate `x` and `y` along the channel dimension (third dimension).
Equivalent to `cat(x, y; dims=3)`.
Convenient binary reduction operator for use with `Parallel`.
"""
cat_channels(x, y) = cat(x, y; dims = 3)

"""
    skip_projection(inplanes, outplanes, downsample = false)

Create a skip projection
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: the number of output feature maps
- `downsample`: set to `true` to downsample the input
"""
skip_projection(inplanes, outplanes, downsample = false) = downsample ?
  Chain(conv_bn((1, 1), inplanes, outplanes, identity; stride = 2, bias = false)...) :
  Chain(conv_bn((1, 1), inplanes, outplanes, identity; stride = 1, bias = false)...)

# array -> PaddedView(0, array, outplanes) for zero padding arrays
"""
    skip_identity(inplanes, outplanes[, downsample])

Create a identity projection
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: the number of output feature maps
- `downsample`: this argument is ignored but it is needed for compatibility with [`resnet`](#).
"""
function skip_identity(inplanes, outplanes)
  if outplanes > inplanes
    return Chain(MaxPool((1, 1), stride = 2),
                 y -> cat(y, zeros(eltype(y),
                                   size(y, 1),
                                   size(y, 2),
                                   outplanes - inplanes, size(y, 4)); dims = 3))
  else
    return identity
  end
end
skip_identity(inplanes, outplanes, downsample) = skip_identity(inplanes, outplanes)

"""
    mlpblock(planes, hidden_planes; dropout = 0., dense = Dense, activation = gelu)

Feedforward block used in many vision transformer-like models.

# Arguments
- `planes`: Number of dimensions in the input and output.
- `hidden_planes`: Number of dimensions in the intermediate layer.
- `dropout`: Dropout rate.
- `dense`: Type of dense layer to use in the feedforward block.
- `activation`: Activation function to use.
"""
function mlpblock(planes, hidden_planes; dropout = 0., dense = Dense, activation = gelu)
  Chain(dense(planes, hidden_planes, activation), Dropout(dropout),
        dense(hidden_planes, planes, activation), Dropout(dropout))
end

"""
    Attention(in => out)
    Attention(qkvlayer)

Self attention layer used by transformer models. Specify the `in` and `out` dimensions,
or directly provide a `qkvlayer` that maps an input the queries, keys, and values.
"""
struct Attention{T}
  qkv::T
end

Attention(dims::Pair{Int, Int}) = Attention(Dense(dims.first, dims.second * 3; bias = false))

@functor Attention

function (attn::Attention)(x::AbstractArray{T, 3}) where T
  q, k, v = chunk(attn.qkv(x), 3, dims = 1)
  scale = convert(T, sqrt(size(q, 1)))
  score = softmax(batched_mul(batched_transpose(q), k) / scale)
  attention = batched_mul(v, score)

  return attention
end

struct MHAttention{S, T}
  heads::S
  projection::T
end

"""
    MHAttention(in, hidden, nheads; dropout = 0.0)

Multi-head self-attention layer used in many vision transformer-like models.

# Arguments
- `in`: Number of dimensions in the input.
- `hidden`: Number of dimensions in the intermediate layer.
- `nheads`: Number of attention heads.
- `dropout`: Dropout rate for the projection layer.
"""
function MHAttention(in, hidden, nheads; dropout = 0.)
  if (nheads == 1 && hidden == in)
    return Attention(in => in)
  end
  inheads, innerheads = chunk(1:in, nheads), chunk(1:hidden, nheads)
  heads = Parallel(vcat, [Attention(length(i) => length(o)) for (i, o) in zip(inheads, innerheads)]...)
  projection = Chain(Dense(hidden, in), Dropout(dropout))

  MHAttention(heads, projection)
end

@functor MHAttention

function (mha::MHAttention)(x)
  nheads = length(mha.heads.layers)
  xhead = chunk(x, nheads, dims = 1)
  return mha.projection(mha.heads(xhead...))
end

"""
    PatchEmbedding(patch_size)
    PatchEmbedding(patch_height, patch_width)

Patch embedding layer used by many vision transformer-like models to split the input image into patches.
"""
struct PatchEmbedding
  patch_height::Int
  patch_width::Int
end

PatchEmbedding(patch_size) = PatchEmbedding(patch_size, patch_size)

function (p::PatchEmbedding)(x)
  h, w, c, n = size(x)
  hp, wp = h ÷ p.patch_height, w ÷ p.patch_width
  xpatch = reshape(x, hp, p.patch_height, wp, p.patch_width, c, n)

  return reshape(permutedims(xpatch, (1, 3, 5, 2, 4, 6)), p.patch_height * p.patch_width * c,
                 hp * wp, n)
end

@functor PatchEmbedding

"""
    ViPosEmbedding(embedsize, npatches; init = (dims) -> rand(Float32, dims))

Positional embedding layer used by many vision transformer-like models.
"""
struct ViPosEmbedding{T}
  vectors::T
end

ViPosEmbedding(embedsize, npatches; init = (dims::NTuple{2, Int}) -> rand(Float32, dims)) =
  ViPosEmbedding(init((embedsize, npatches)))

(p::ViPosEmbedding)(x) = x .+ p.vectors

@functor ViPosEmbedding

"""
    ClassTokens(dim; init = Flux.zeros32)

Appends class tokens to an input with embedding dimension `dim` for use in many vision transformer models.
"""
struct ClassTokens{T}
  token::T
end

ClassTokens(dim::Integer; init = Flux.zeros32) = ClassTokens(init(dim, 1, 1))

function (m::ClassTokens)(x)
  tokens = repeat(m.token, 1, 1, size(x, 3))
  return hcat(tokens, x)
end

@functor ClassTokens


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of issues.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Issue #65 - Data type in filters with Tables.jl" begin
    table = (A = Int64.(1:20), B = Float64.(1:20))
    rowfilter1(data,i) = i <= div(length(first(data)),2)
    rowfilter2(data,i) = i <= div(length(data),2)

    result1 = pretty_table(String, table;                 filters_row = (rowfilter1,))
    result2 = pretty_table(String, Tables.columns(table); filters_row = (rowfilter1,))
    result3 = pretty_table(String, Tables.rows(table);    filters_row = (rowfilter2,))

    expected = """
┌───────┬─────────┐
│     A │       B │
│ Int64 │ Float64 │
├───────┼─────────┤
│     1 │     1.0 │
│     2 │     2.0 │
│     3 │     3.0 │
│     4 │     4.0 │
│     5 │     5.0 │
│     6 │     6.0 │
│     7 │     7.0 │
│     8 │     8.0 │
│     9 │     9.0 │
│    10 │    10.0 │
└───────┴─────────┘
"""

    @test result1 == expected
    @test result2 == expected
    @test result3 == expected

    table = (
        A = Int64.(1:10),
        B = Float64.(1:10),
        C = Int64.(1:10),
        D = Float64.(1:10)
    )
    colfilter1(data,j) = j <= div(length(data),2)
    colfilter2(data,j) = j <= div(length(first(data)),2)

    result1 = pretty_table(String, table;                 filters_col = (colfilter1,))
    result2 = pretty_table(String, Tables.columns(table); filters_col = (colfilter1,))
    result3 = pretty_table(String, Tables.rows(table);    filters_col = (colfilter2,))

    @test result1 == expected
    @test result2 == expected
    @test result3 == expected
end

struct SimpleTable{T}
    data::Matrix{T}
end

Tables.istable(::SimpleTable) = true
Tables.columnaccess(::SimpleTable) = true
Tables.columnnames(x::SimpleTable) = [Symbol(i) for i = 1:size(x.data, 2)]
Tables.columns(x::SimpleTable) = x
Tables.getcolumn(x::SimpleTable, i::Symbol) = tuple(x.data[parse(Int,string(i)),:]...)

table = SimpleTable([10.0^(i+j) for i in 1:10, j in 1:5])

@testset "Issue #90 - Tables.jl returning tuples as columns" begin
    expected = """
┌──────────┬──────────┬──────────┬──────────┬────────┐
│        1 │        2 │        3 │        4 │      5 │
├──────────┼──────────┼──────────┼──────────┼────────┤
│    100.0 │   1000.0 │  10000.0 │ 100000.0 │  1.0e6 │
│   1000.0 │  10000.0 │ 100000.0 │    1.0e6 │  1.0e7 │
│  10000.0 │ 100000.0 │    1.0e6 │    1.0e7 │  1.0e8 │
│ 100000.0 │    1.0e6 │    1.0e7 │    1.0e8 │  1.0e9 │
│    1.0e6 │    1.0e7 │    1.0e8 │    1.0e9 │ 1.0e10 │
└──────────┴──────────┴──────────┴──────────┴────────┘
"""

    result = pretty_table(String, table)
    @test result == expected
end


text:
module AnnotatedArrays

export
	AnnotatedArray,
	array,
	annotationnames,
	annotate!,
	annot,
	inset,
	showall

include("annotatedarray.jl")

end


text:
##
using IFGF
using LinearAlgebra
using StaticArrays
using Random
using Plots
plotlyjs()
Random.seed!(1)

includet(joinpath(IFGF.PROJECT_ROOT,"test","testutils.jl"))

k    = 8π
λ    = 2π/k
ppw  = 10
dx   = λ/ppw

Xpts = Ypts = sphere_uniform_surface_mesh(dx)
npts = length(Xpts)

K    = HelmholtzKernel(k)
T    = ComplexF64

I   = rand(1:npts,100)
B   = randn(T,npts)
exa = [sum(K( Xpts[i],Ypts[j])*B[j] for j in 1:npts) for i in I]

# trees
splitter   = Trees.CardinalitySplitter(;nmax=200)

# cone size function
ds_func = IFGF.cone_domain_size_func(k,(1.0,π/2,π/2))

# loop over interpolation order
pvec = 2:8
er   = []
for p in pvec
    ptuple  = (p,p,p)
    A       = IFGFOp(K,Xpts,Ypts;splitter,p=ptuple,ds_func,threads=true)
    C       = zeros(T,npts)
    mul!(C,A,B,1,0;threads=true)
    ee = norm(C[I]-exa,2) / norm(exa,2)
    push!(er,ee)
end

print(
    "Uniform points on surfaces of sphere with parameters:
        r   = 1
        k   = $k
        n   = $npts
        ppw = $ppw"
)

fig = plot(pvec,er,yscale=:log10,m=:circle,
    ylabel="approximate ℓ₂ relative error",
    xlabel="p",
    label="sphere")

dir = @__DIR__
savefig(fig,joinpath(dir,"sphere_p_convergence.png"))

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
@testset "SSAD" begin

    dummy_data = reshape(collect(1:20), 2, 10)
    pools = fill(:U, 10)
    ssad = SVDD.SSAD(dummy_data, pools)

    @testset "create" begin
        @test ssad.state == SVDD.model_created
        @test_throws SVDD.ModelStateException SVDD.predict(ssad, dummy_data)
    end

    @testset "initialize" begin
        C1 = C2 = 1.0
        gamma = 0.5
        init_strategy = SVDD.FixedParameterInitialization(MLKernels.GaussianKernel(gamma), (C1, C2))

        SVDD.initialize!(ssad, init_strategy)
        @test ssad.state == SVDD.model_initialized
        @test ssad.data == dummy_data
        @test ssad.C1 ≈ 1.0
        @test ssad.C2 ≈ 1.0
        @test size(ssad.K) == (10, 10)
        @test all(diag(ssad.K) .≈ 1.0)
    end

    @testset "fit" begin
        SVDD.fit!(ssad, TEST_SOLVER)
        @test ssad.state == SVDD.model_fitted
        @test length(ssad.alpha_values) == size(dummy_data,2)
        @test ssad.ρ > 0.0
    end

    @testset "change pools 1" begin
        pools[1] = :Lin
        SVDD.set_pools!(ssad, pools)
        @test ssad.state == SVDD.model_initialized
        SVDD.fit!(ssad, TEST_SOLVER)
        @test ssad.state == SVDD.model_fitted
    end

    @testset "change pools 2" begin
        pools = fill(:U, 10)
        pools[1] = :Lout
        SVDD.set_pools!(ssad, pools)
        @test ssad.state == SVDD.model_initialized
        SVDD.fit!(ssad, TEST_SOLVER)
        @test ssad.state == SVDD.model_fitted
    end

    @testset "change pools 3" begin
        pools = fill(:U, 10)
        pools[1] = :Lout
        pools[10] = :Lin
        SVDD.set_pools!(ssad, pools)
        @test ssad.state == SVDD.model_initialized
        SVDD.fit!(ssad, TEST_SOLVER)
        @test ssad.state == SVDD.model_fitted
    end

    @testset "predict" begin
        predictions = SVDD.predict(ssad, dummy_data)
        @test length(predictions) == size(dummy_data, 2)
        @test predictions[1] == SVDD.predict(ssad, hcat([1.0, 2.0]))[1]
    end

    @testset "params" begin
        @test all(map(x->SVDD.is_valid_param_value(ssad, Val{:C1}, x), [0.0, 0.1, 1.0]))
        @test all(map(x->SVDD.is_valid_param_value(ssad, Val{:C2}, x), [0.0, 0.1, 1.0]))
        @test !any(map(x->SVDD.is_valid_param_value(ssad, Val{:C1}, x), [-0.1, 1.1]))
        @test !any(map(x->SVDD.is_valid_param_value(ssad, Val{:C2}, x), [-0.1, 1.1]))
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ a31887be-d60b-41a8-a561-3ff0f6491e7b
using CSV, DataFrames

# ╔═╡ 1214f8e0-d687-11eb-296c-977d3ccd3c36
md"""
# Normalizing data tables

The popularity of the [tidyverse](https://github.com/tidyverse/tidyverse) collection of packages for R has led to many public data sets being presented as tidy data frames.

## Redundancies and possible inconsistency

As an organizing principle, having a tidy data frame is desirable but it often results in redundancies in the data, as we will see below.
Redundancies are usually not a problem in terms of the amount of storage taken up; storage, even on a public cloud, is very low cost now.
However, the redundancies present the possibility of accidental inconsistencies, which we do want to avoid.

Those working with relational data bases often emphasize the use of a *normalized* database where redundant information is eliminated and the table to be used, e.g. for fitting models, is generated by *joining* several tables of non-redundant information together.

The [DataFrames.jl package](https://github.com/JuliaData/DataFrames.jl) for Julia provides the ability to join tables flexibly and efficiently.  In fact, much of the recent development work on that package, primarily by Bogumił Kamiński, has been in making joins more efficient.

### General approach

The general approach is to isolate those covariates or columns that are characteristics of the subject in a *subject table*.
In terms of the experimental design these covariates are *between-subject* in that they vary between the subjects but not within a subject.
For example, demographic variables are usually between-subject.

Similarly the covariates that are characteristic of the items are *between-item* and are stored in an *item table*.
Word frequency might be an example of a *between-item* covariate.

The *observation table* is a tidy table containing one row for each observation of the response plus the subject identifier, the item identifier and any covariates that vary within-subject and within-item.

When between-subject or between-item covariates are required for exploratory analysis or model fitting the model table is generated from the observation table, the subject table and the item table by a join.

### Data set

First load the packages to be used, read the .csv file, and convert it to a DataFrame.
"""

# ╔═╡ 1ad1a572-22e4-4c9d-b893-95f270c795d0
orig = CSV.read("../../Dorothea_Pregla/SPL_decl.csv", DataFrame; missingstring="NA")

# ╔═╡ 96410f16-b1a8-4207-802b-abb7eeb1d919
size(orig)

# ╔═╡ 07967240-57aa-42ac-aa4f-c74d83055986
describe(orig)

# ╔═╡ 07cbb2a4-7c37-4a31-8b10-88f478b0b4f8
length(unique(orig.subj))

# ╔═╡ 3fd79669-b8fb-40e7-ba5d-6d9339e01f50
md"""
These data were provided by Dorothea Pregla.
Her description of the columns is:
- subj: subject number
- item: number of the item (number of items: n=20)
- condition: a (canonical) and b (non-canonical), we did not use a latin square design, i.e., all participants saw both conditions of all items
- region: region in the sentence (0, 1, 2, 3, question)
  - 0: sentence onset ("hier")
  - 1: verb ("tröstet")
  - 2: NP1 ("der Tiger")
  - 3: adverb ("gerade")
  - question: NP2 + picture selection ("den Esel" + picture selection)
- subj_status: healthy control (HC) or individual with aphasia (IWA)
- sentence: the respective auditory presented sentence
- word: the respective word of the sentence
- rt: listening time for the respective region
- acc: picture selection accuracy
- age, sex, years of education: demographic data of the participant
- aphasia_type: Aachen aphasia test syndrome, NA for control participants
- digit_symbol_substitution: WAIS subtest, number of correct items in two minutes
- digit span, block span: WMS-R subtests, forward and backward raw scores and longest span
- picture span: working memory test, forward and backward raw scores and longest span
- lexdec: lexical decision task, mean reaction times in the task

Many of these columns are characteristics of the subject, of whom there are 71.
These characteristics are repeated to fill out the 14185 rows of the corresponding column in the table.  This is the nature of the redundancy.

Furthermore, if there was some kind of an error and, say, a subject's age was recorded inconsistently in one or more rows of this table, this would not be automatically detected in this format.

A way to enforce consistency is to create one table for the subjects with 71 rows and all the subject-specific columns recorded only once in there.

Those who collected the data would know what columns are subject-specific and we could procede to `select` those columns, along with the `subj` column, and apply `unique` to get the subject table.
We should, of course, check that the number of unique rows in the subject table is exactly 71, the number of subjects.

For the sake of illustration, we will show a way of determining which columns are properties of the subject for which we will use a `GroupedDataFrame`.
"""

# ╔═╡ fb224e36-4427-43f9-a435-d8a7c7c880a3
gdf = groupby(orig, :subj)

# ╔═╡ a90307e7-764c-41e1-a77a-7f6fc30c275b
md"""
We are trying to find the columns that are constant for each subject.
For example, the `subj_status` should be constant for each subject but the `condition` is not.

One way of checking if a vector is constant is to compare all the values to the first one.
A simple version of an `isconstant` function is
"""

# ╔═╡ 13fc105c-3e6e-492b-859c-08f4801c4bdb
isconstant1(v) = all(==(first(v)), v)

# ╔═╡ d3e8956c-87fc-406c-9099-ed78a46d98fd
md"""
where `==(first(v))` creates a function of the form x -> x == first(v) but with a bit less overhead.

This works for simple cases
"""

# ╔═╡ 2e43c5ac-ba4b-4214-863b-cc5fa081169b
isconstant1(first(gdf).age)

# ╔═╡ 934fe224-91fb-41f1-808a-62ea9fc63146
isconstant1(first(gdf).item)

# ╔═╡ 0dc9be9d-28af-4062-945a-61e529200647
md"""
but it has (at least) two problems on edge cases.  First, what if `v` is empty?  The call to `first(v)` will fail.
"""

# ╔═╡ c55be9cf-83e8-4d36-b2c8-01763cfd2020
isconstant1([])

# ╔═╡ 07689442-b01f-4727-84c9-921ee669d7db
md"""
It is reasonably straightforward to avoid this situation by checking first if `v` is empty, in which case we will return `true`, arguing that a vector with no elements is constant.
"""

# ╔═╡ 2ebdf2c4-3234-4e1a-a122-10603b77e4ef
function isconstant2(v)
	isempty(v) || all(==(first(v)), v)
end

# ╔═╡ 757ba155-6dc5-49c5-beea-826ec8efa560
isconstant2([])

# ╔═╡ 0021dfba-e466-4ff4-8f3d-f521e47bbbed
md"""
The second edge case involves missing data.  Comparisons with `missing` always return `missing` not `true` or `false`.
"""

# ╔═╡ 652cf724-f1fa-46c9-be97-8b52c71b606c
isconstant2(first(gdf).aphasia_type)

# ╔═╡ 24dcaf32-e907-41ab-a7b9-099f6b329d2f
md"""
If all of the elements of `v` are missing we conclude that `v` is constant otherwise if the comparison returns `missing` then it is not constant.
We write our final version of `isconstant` using `coalesce`, which is a function that returns the first non-missing value in its argument list.
"""

# ╔═╡ d85cbb0d-6671-4ee0-a013-d5927a33b1bb
function isconstant(v)
	return isempty(v) || all(ismissing, v) || coalesce(all(==(first(v)), v), false)
end

# ╔═╡ 56222937-e45e-4e7a-b49b-d5fb9bace1a7
md"Next we create a method that takes a `GroupedDataFrame` and a name `nm` and iterates over all the subdataframes in `gdf`."

# ╔═╡ a237208d-1456-4fc6-95e8-c1bd10d20c70
isconstant(gdf::GroupedDataFrame, nm) = all(sdf -> isconstant(sdf[!, nm]), gdf)

# ╔═╡ a1c38547-9d80-472c-af21-0e9db3dd41d7
isconstant(gdf, :subj_status)

# ╔═╡ 6bb74de8-d8ea-4674-8419-35a33d8f7533
isconstant(gdf, :item)

# ╔═╡ dbef0ab4-e103-438d-93c2-98d4bd4920ad
isconstant(gdf, :lexdec)

# ╔═╡ 0045e306-8123-4af0-a6e5-70407ef9e6fb
md"""
Finally, we iterate over all the subdataframes in the `GroupedDataFrame` and all the column names, checking this property.
"""

# ╔═╡ be1ec22f-c3f4-447e-8bf4-a301db064b9f
subj_cols = filter(nm -> isconstant(gdf, nm), names(orig))

# ╔═╡ 11d893ed-d4f3-43a6-b0bf-1746336936dd
md"""
Create `subjtable` and check that it has the correct number of rows (71).
"""

# ╔═╡ 242b0b14-4157-4335-b30c-295c2374d5a4
subjtable = unique(select(orig, subj_cols))

# ╔═╡ ef685a16-c24f-4f8b-958d-3f4eae9ec76e
md"""
Next we check if there are columns that are constant across items.
"""

# ╔═╡ 65f93aff-5453-4f90-93d4-d6ad1686abd8
item_cols = let
	gdfitem = groupby(orig, :item)
	filter(nm -> isconstant(gdfitem, nm), names(orig))
end

# ╔═╡ 77a7b199-adf2-426d-8df5-c9327b61aab2
md"""
So `item` is the only column that is constant within items.

Finally we create a response table by removing the columns from `subjtable`, except for `subj` itself.
"""

# ╔═╡ 8e4cadd0-b5f9-42b5-a8b7-bc28551a35e6
responsetable = select(orig, :subj, setdiff(names(orig), subj_cols))

# ╔═╡ 92aacaba-f831-4a8e-a023-d5920de6f0b4
names(responsetable)

# ╔═╡ 1479ed0c-8d07-4fbf-a897-05e5a937127d
md"""
To regenerate the original table (with, perhaps, a different ordering of the columns), we use a `leftjoin` on the `responsetable` and `subjtable` on `:subj`.
"""

# ╔═╡ 6bc71c0f-1c0b-4995-a245-6044a2d275f4
begin
	newtable = leftjoin(responsetable, subjtable, on = :subj)
	describe(newtable)
end

# ╔═╡ 827ac2ca-70e7-4ff1-b414-b8310669fdf4
md"""
## General approach

In most cases of "subject/item" data like this we expect to have 3 tables: a subject table that contains the subject identifier and any columns that are characteristics of the subjects, an item table, and an observation table that contains the subject, the item, any within-subject/within-item experimental conditions and the response.
In this case because the item table is trivial, containing only the item column itself, it can be omitted.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

[compat]
CSV = "~0.8.5"
DataFrames = "~1.1.1"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.6.2"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.31.0"

[[deps.Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[deps.DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "66ee4fe515a9294a8836ef18eea7239c6ac3db5e"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.1.1"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.9"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InvertedIndices]]
deps = ["Test"]
git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.0.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.0"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.2.1"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "ffae887d0f0222a19c406a11c3831776d1383e3d"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.3"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "8ed4a3ea724dac32670b062be3ef1c1de6773ae8"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.4.4"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─1214f8e0-d687-11eb-296c-977d3ccd3c36
# ╠═a31887be-d60b-41a8-a561-3ff0f6491e7b
# ╠═1ad1a572-22e4-4c9d-b893-95f270c795d0
# ╠═96410f16-b1a8-4207-802b-abb7eeb1d919
# ╠═07967240-57aa-42ac-aa4f-c74d83055986
# ╠═07cbb2a4-7c37-4a31-8b10-88f478b0b4f8
# ╟─3fd79669-b8fb-40e7-ba5d-6d9339e01f50
# ╠═fb224e36-4427-43f9-a435-d8a7c7c880a3
# ╟─a90307e7-764c-41e1-a77a-7f6fc30c275b
# ╠═13fc105c-3e6e-492b-859c-08f4801c4bdb
# ╟─d3e8956c-87fc-406c-9099-ed78a46d98fd
# ╠═2e43c5ac-ba4b-4214-863b-cc5fa081169b
# ╠═934fe224-91fb-41f1-808a-62ea9fc63146
# ╟─0dc9be9d-28af-4062-945a-61e529200647
# ╠═c55be9cf-83e8-4d36-b2c8-01763cfd2020
# ╟─07689442-b01f-4727-84c9-921ee669d7db
# ╠═2ebdf2c4-3234-4e1a-a122-10603b77e4ef
# ╠═757ba155-6dc5-49c5-beea-826ec8efa560
# ╟─0021dfba-e466-4ff4-8f3d-f521e47bbbed
# ╠═652cf724-f1fa-46c9-be97-8b52c71b606c
# ╟─24dcaf32-e907-41ab-a7b9-099f6b329d2f
# ╠═d85cbb0d-6671-4ee0-a013-d5927a33b1bb
# ╟─56222937-e45e-4e7a-b49b-d5fb9bace1a7
# ╠═a237208d-1456-4fc6-95e8-c1bd10d20c70
# ╠═a1c38547-9d80-472c-af21-0e9db3dd41d7
# ╠═6bb74de8-d8ea-4674-8419-35a33d8f7533
# ╠═dbef0ab4-e103-438d-93c2-98d4bd4920ad
# ╟─0045e306-8123-4af0-a6e5-70407ef9e6fb
# ╠═be1ec22f-c3f4-447e-8bf4-a301db064b9f
# ╟─11d893ed-d4f3-43a6-b0bf-1746336936dd
# ╠═242b0b14-4157-4335-b30c-295c2374d5a4
# ╟─ef685a16-c24f-4f8b-958d-3f4eae9ec76e
# ╠═65f93aff-5453-4f90-93d4-d6ad1686abd8
# ╟─77a7b199-adf2-426d-8df5-c9327b61aab2
# ╠═8e4cadd0-b5f9-42b5-a8b7-bc28551a35e6
# ╠═92aacaba-f831-4a8e-a023-d5920de6f0b4
# ╟─1479ed0c-8d07-4fbf-a897-05e5a937127d
# ╠═6bc71c0f-1c0b-4995-a245-6044a2d275f4
# ╟─827ac2ca-70e7-4ff1-b414-b8310669fdf4
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
This file contains all the functions that are required to do basic processing of
climate data, such as averaging of data into hourly, 3-hourly and daily bins, or
in the creation of Hovmöller diagrams and frequency bins.

This function assumes that the datasets (for the averaging/summing of data)
follow the following formats:
    - dimensions in the form (longitude,latitude,time)
    - contains data for a single day, for daily and hourly sums

"""

function convert2hourly(data::Array,isavg=false)
    nlon = size(data,1); nlat = size(data,2); nt = size(data,3); dt = Int(nt/24);
    data = reshape(data,(nlon,nlat,dt,:));
    if isavg; data = mean(data,dims=3);
    else;     data = sum(data,dims=3);
    end
end

function convert23hourly(data::Array,isavg=false)
    nlon = size(data,1); nlat = size(data,2); nt = size(data,3); dt = Int(nt/8);
    data = reshape(data,(nlon,nlat,dt,:));
    if isavg; data = mean(data,dims=3);
    else;     data = sum(data,dims=3);
    end
end

function convert2daily(data::Array,isavg=false)
    if isavg; data = mean(data,dims=3);
    else;     data = sum(data,dims=3);
    end
end

function convert2regionmean(data::Array,isavg=true)
    if isavg; data = mean(data,dims=(1,2));
    else;     data = sum(data,dims=(1,2));
    end
end

function convert2zonalmean(data::Array,isavg=true)
    if isavg; data = mean(data,dims=1);
    else;     data = sum(data,dims=1);
    end
end

function convert2meridionalmean(data::Array,isavg=true)
    if isavg; data = mean(data,dims=2);
    else;     data = sum(data,dims=2);
    end
end


text:
##
const ZEROS_DOC = """
    @zeros(args...)

Call `zeros(numbertype, args...)`, where `numbertype` is the datatype selected with [`@init_parallel_kernel`](@ref) and the function `zeros` is chosen to be compatible with the package for parallelization selected with [`@init_parallel_kernel`](@ref) (zeros for Threads and CUDA.zeros for CUDA).
"""
@doc ZEROS_DOC
macro zeros(args...) check_initialized(); esc(_zeros(args...)); end

##
const ONES_DOC = """
    @ones(args...)

Call `ones(numbertype, args...)`, where `numbertype` is the datatype selected with [`@init_parallel_kernel`](@ref) and the function `ones` is chosen to be compatible with the package for parallelization selected with [`@init_parallel_kernel`](@ref) (ones for Threads CUDA.ones for CUDA).
"""
@doc ONES_DOC
macro ones(args...) check_initialized(); esc(_ones(args...)); end

##
const RAND_DOC = """
    @rand(args...)

Call `rand(numbertype, args...)`, where `numbertype` is the datatype selected with [`@init_parallel_kernel`](@ref) and the function `rand` is chosen/implemented to be compatible with the package for parallelization selected with [`@init_parallel_kernel`](@ref).
"""
@doc RAND_DOC
macro rand(args...) check_initialized(); esc(_rand(args...)); end

## MACROS FORCING PACKAGE, IGNORING INITIALIZATION

macro zeros_cuda(args...)    check_initialized(); esc(_zeros(args...; package=PKG_CUDA)); end
macro ones_cuda(args...)     check_initialized(); esc(_ones(args...; package=PKG_CUDA)); end
macro rand_cuda(args...)     check_initialized(); esc(_rand(args...; package=PKG_CUDA)); end
macro zeros_threads(args...) check_initialized(); esc(_zeros(args...; package=PKG_THREADS)); end
macro ones_threads(args...)  check_initialized(); esc(_ones(args...; package=PKG_THREADS)); end
macro rand_threads(args...)  check_initialized(); esc(_rand(args...; package=PKG_THREADS)); end

## ALLOCATOR FUNCTIONS

function _zeros(args...; package::Symbol=get_package())
    numbertype = get_numbertype()
    if     (package == PKG_CUDA)    return :(CUDA.zeros($numbertype, $(args...)))
    elseif (package == PKG_THREADS) return :(Base.zeros($numbertype, $(args...)))
    else                            @KeywordArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
end

function _ones(args...; package::Symbol=get_package())
    numbertype = get_numbertype()
    if     (package == PKG_CUDA)    return :(CUDA.ones($numbertype, $(args...)))
    elseif (package == PKG_THREADS) return :(Base.ones($numbertype, $(args...)))
    else                            @KeywordArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
end

function _rand(args...; package::Symbol=get_package())
    numbertype = get_numbertype()
    if     (package == PKG_CUDA)    return :(CUDA.CuArray(rand($numbertype, $(args...))))
    elseif (package == PKG_THREADS) return :(Base.rand($numbertype, $(args...)))
    else                            @KeywordArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
end