text:
# ---
# title: 1594. Maximum Non Negative Product in a Matrix
# id: problem1594
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Greedy
# link: <https://leetcode.com/problems/maximum-non-negative-product-in-a-matrix/description/>
# hidden: true
# ---
#
# You are given a `rows x cols` matrix `grid`. Initially, you are located at the
# top-left corner `(0, 0)`, and in each step, you can only **move right  or
# down** in the matrix.
#
# Among all possible paths starting from the top-left corner `(0, 0)` and ending
# in the bottom-right corner `(rows - 1, cols - 1)`, find the path with the
# **maximum non-negative product**. The product of a path is the product of all
# integers in the grid cells visited along the path.
#
# Return the  _maximum non-negative product   **modulo**  _`109 + 7`.  _If the
# maximum product is **negative** return _`-1`.
#
# **Notice that the modulo is performed after getting the maximum product.**
#
#
#
# **Example 1:**
#
#
#
#     Input: grid = [[-1,-2,-3],
#                    [-2,-3,-3],
#                    [-3,-3,-2]]
#     Output: -1
#     Explanation: It's not possible to get non-negative product in the path from (0, 0) to (2, 2), so return -1.
#
#
# **Example 2:**
#
#
#
#     Input: grid = [[ **1** ,-2,1],
#                    [ **1** , **-2** ,1],
#                    [3, **-4** , **1** ]]
#     Output: 8
#     Explanation: Maximum non-negative product is in bold (1 * 1 * -2 * -4 * 1 = 8).
#
#
# **Example 3:**
#
#
#
#     Input: grid = [[ **1** , 3],
#                    [ **0** , **-4** ]]
#     Output: 0
#     Explanation: Maximum non-negative product is in bold (1 * 0 * -4 = 0).
#
#
# **Example 4:**
#
#
#
#     Input: grid = [[ **1** , 4,4,0],
#                    [ **-2** , 0,0,1],
#                    [ **1** , **-1** , **1** , **1** ]]
#     Output: 2
#     Explanation: Maximum non-negative product is in bold (1 * -2 * 1 * -1 * 1 * 1 = 2).
#
#
#
#
# **Constraints:**
#
#   * `1 <= rows, cols <= 15`
#   * `-4 <= grid[i][j] <= 4`
#
#
## @lc code=start
using LeetCode

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


text:
using Documenter, Hashpipe

makedocs(sitename="Hashpipe Documentation",
         modules = [Hashpipe],
)

deploydocs(
    repo = "github.com/max-Hawkins/Hashpipe.jl.git",
    devbranch = "main"
)

ext: jl


text:
abstract type Coordinates end
struct Polar <: Coordinates end
struct Cartesian <: Coordinates end

circ(r) = if r < 1 1 else 0 end
circ(::Polar, r,θ) = circ(r)
circ(::Cartesian, x,y) = circ(x^2+y^2)
"""
    circ(x,y)

Circular mask of radius 1.
"""
circ(x,y) = circ(Cartesian(), x, y)

Base.convert(::Type{Coordinates}, s::Symbol) = begin
    if s in [:cartesian :cart :xy]
        Cartesian()
    elseif s in [:polar :pol :rθ]
        Polar()
    else
        error("$s not convertible in coordinates.")
    end
end


text:
# Gondzio's multiple centrality correctors method

function update_rxs!(rxs_l, rxs_u, Hmin, Hmax, x_m_l_αΔp, u_m_x_αΔp, s_l_αΔp, s_u_αΔp, nlow, nupp)

    @inbounds @simd for i=1:nlow
        rxs_l[i] = s_l_αΔp[i] * x_m_l_αΔp[i]
        if Hmin <= rxs_l[i] <= Hmax
            rxs_l[i] = 0
        elseif rxs_l[i] < Hmin
            rxs_l[i] -= Hmin
        else
            rxs_l[i] -= Hmax
        end
        if rxs_l[i] > Hmax
            rxs_l[i] = Hmax
        end
    end
    @inbounds @simd for i=1:nupp
        rxs_u[i] = -s_u_αΔp[i]*u_m_x_αΔp[i]
        if Hmin <= -rxs_u[i] <= Hmax
            rxs_u[i] = 0
        elseif -rxs_u[i] < Hmin
            rxs_u[i] += Hmin
        else
            rxs_u[i] += Hmax
        end
        if rxs_u[i] < -Hmax
            rxs_u[i] = -Hmax
        end
    end
end

function multi_centrality_corr!(dda :: DescentDirectionAllocsPC{T}, pad :: PreallocatedData{T}, pt :: Point{T}, α_pri :: T, α_dual :: T,
                                itd :: IterData{T}, fd :: Abstract_QM_FloatData{T}, id :: QM_IntData, cnts :: Counters, res :: Residuals{T},
                                T0 :: DataType) where {T<:Real}

    iter_c = 0 # current number of correction iterations
    corr_flag = true #stop correction if false
    # for storage issues Δ_aff = Δp  and Δ_cc = Δm
    dda.Δxy_aff .= itd.Δxy
    dda.Δs_l_aff .= itd.Δs_l
    dda.Δs_u_aff .= itd.Δs_u
    @inbounds while iter_c < cnts.kc && corr_flag
        # Δp = Δ_aff + Δ_cc
        δα, γ, βmin, βmax = T(0.1), T(0.1), T(0.1), T(10)
        α_p2, α_d2 = min(α_pri + δα, one(T)), min(α_dual + δα, one(T))
        update_pt_aff!(dda.x_m_l_αΔ_aff, dda.u_m_x_αΔ_aff, dda.s_l_αΔ_aff, dda.s_u_αΔ_aff, dda.Δxy_aff, dda.Δs_l_aff, dda.Δs_u_aff,
                        itd.x_m_lvar, itd.uvar_m_x, pt.s_l, pt.s_u, α_p2, α_d2, id.ilow, id.iupp)
        μ_p = compute_μ(dda.x_m_l_αΔ_aff, dda.u_m_x_αΔ_aff, dda.s_l_αΔ_aff, dda.s_u_αΔ_aff, id.nlow, id.nupp)

        σ = (μ_p / itd.μ)^3
        Hmin, Hmax = βmin * σ * itd.μ, βmax * σ * itd.μ

        # corrector-centering step
        update_rxs!(dda.rxs_l, dda.rxs_u, Hmin, Hmax, dda.x_m_l_αΔ_aff, dda.u_m_x_αΔ_aff, dda.s_l_αΔ_aff, dda.s_u_αΔ_aff, id.nlow, id.nupp)
        itd.Δxy .= 0
        itd.Δxy[id.ilow] .+= dda.rxs_l ./ itd.x_m_lvar
        itd.Δxy[id.iupp] .+= dda.rxs_u ./ itd.uvar_m_x
        out = solver!(pad, dda, pt, itd, fd, id, res, cnts, T0, :cc)
        itd.Δs_l .= @views .-(dda.rxs_l .+ pt.s_l .* itd.Δxy[id.ilow]) ./ itd.x_m_lvar
        itd.Δs_u .= @views (dda.rxs_u .+ pt.s_u .* itd.Δxy[id.iupp]) ./ itd.uvar_m_x

        itd.Δxy .+= dda.Δxy_aff
        itd.Δs_l .+= dda.Δs_l_aff
        itd.Δs_u .+= dda.Δs_u_aff
        α_p2, α_d2 = compute_αs(pt.x, pt.s_l, pt.s_u, fd.lvar, fd.uvar, itd.Δxy, itd.Δs_l, itd.Δs_u, id.nvar)

        if α_p2 >= α_pri + γ*δα && α_d2 >= α_dual + γ*δα
            iter_c += 1
            dda.Δxy_aff .= itd.Δxy
            dda.Δs_l_aff .= itd.Δs_l
            dda.Δs_u_aff .= itd.Δs_u
            α_pri, α_dual = α_p2, α_d2
        else
            itd.Δxy .= dda.Δxy_aff
            itd.Δs_l .= dda.Δs_l_aff
            itd.Δs_u .= dda.Δs_u_aff
            corr_flag = false
        end
    end
    return α_pri, α_dual
end

# function to determine the number of centrality corrections (Gondzio's procedure)
function nb_corrector_steps(J_colptr, ncon, nvar, T)
    Ef, Es, rfs = 0, 16 * nvar, zero(T) # 14n = ratio tests and vector initializations
    @inbounds @simd for j=1:ncon+nvar
        lj = (J_colptr[j+1]-J_colptr[j])
        Ef += lj^2
        Es += lj
    end
    rfs = T(Ef / Es)
    if rfs <= 10
        kc = 0
    elseif 10 < rfs <= 30
        kc = 1
    elseif 30 < rfs <= 50
        kc = 2
    elseif rfs > 50
        kc = 3
    else
        p = Int(rfs / 50)
        kc = p + 2
        if kc > 10
            kc = 10
        end
    end
    return kc
end


text:
#Local function used to transform data via appropriate influence functions for the case where bootstrap method is tapered block
function apply_influence_function(x::Vector{T}, bi::BootInput)::Vector{Float64} where {T<:Number}
	if bi.flevel1 == mean
		x_if = x - mean(x, 1)
	elseif bi.flevel1 == sum
		x_if = length(x) * (x - mean(x))
	else
		error("Tapered block bootstrap only implemented for a limited number of cases with known influence functions")
	end
	return(x_if)
end
(apply_influence_function(x::Vector{Vector{T}}, bi::BootInput)::Vector{Vector{Float64}}) where {T<:Number} = Vector{Float64}[ apply_influence_function(x[k], bi) for k = 1:length(x) ]

#Local function used to weight data for the case where bootstrap method is tapered block
function dboot_kernel_weights(bi::BootInput)::Vector{Float64}
    bL = Int(ceil(bi.blocklength))
	kernelInput = Float64[ (1 / bL) * (n - 0.5) for n = 1:bL ]
    kernelWeight = kernel_func_pp2002(kernelInput, bi.bootmethod)
	normTerm = sqrt(bL) / norm(kernelWeight, 2)
    kernelWeight .*= normTerm
	return(kernelWeight)
end
function dboot_weight!(x::Vector{Vector{T}}, bi::BootInput)::Vector{Vector{Float64}} where {T<:Number}
    bL = Int(ceil(bi.blocklength))
    bL <= 1 && return(x)
    w = dboot_kernel_weights(bm)
	length(w) != bL && error("Logic fail. Incorrect length output from dboot_kernel_weights function. Please file an issue.")
    (num_repeat, num_remain) = divrem(size(x, 1), bL)
    wLong = vcat(repeat(w, outer=num_repeat), w[1:num_remain])
    for n = 1:length(x)
        x[n] = wLong .* x[n]
    end
	return(x)
end
#Local function for the two kernel functions proposed in Paparoditis and Politis (2002) "The tapered block bootstrap for general statistics from stationary sequences"
function kernel_func_pp2002_trap(x::Float64)::Float64
    p = 0.43 #Optimal value from PP (2002)
    x < 0 && return(0.0)
    x < p && return(x / p)
    x < 1 - p && return(1.0)
    x < 1 && return((1 - x) / p)
    return(0.0)
end
kernel_func_pp2002_trap(x::Vector{Float64})::Vector{Float64} = Float64[ kernel_func_pp2002_trap(x[k]) for k = 1:length(x) ]
function kernel_func_pp2002_smooth(x::Float64)::Float64
    p = 1.3 #Optimal value from PP (2002)
    (0.0 <= x <= 1.0) && return(1 - abs(2*x - 1)^p)
    return(0.0)
end
kernel_func_pp2002_smooth(x::Vector{Float64})::Vector{Float64} = Float64[ kernel_func_pp2002_smooth(x[k]) for k = 1:length(x) ]
function kernel_func_pp2002(x::Vector{Float64}, bm::BootTapered)::Vector{Float64}
    bm.kernelfunction == :trap && return(kernel_func_pp2002_trap(x))
    bm.kernelfunction == :smooth && return(kernel_func_pp2002_smooth(x))
    error("Invalid kernel function of $(bm.kernelfunction)")
end


text:
# Reference : Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application
using IntervalArithmetic, ForwardDiff
import Base: +, -, *, /, ^, sqrt, exp, log, sin, cos, tan, asin, acos, atan
import IntervalArithmetic: mid, interval

function slope(f::Function, x::Interval, c::Number)
    f(slope_var(x, c)).s
end

struct Slope{T}
    x::Interval{T}  # Interval on which slope is evaluated
    c::Interval{T}  # Point about which slope is evaluated (Interval to get bounded rounding errors)
    s::Interval{T}  # Variable which propogates the slope information
end

Slope(c) = Slope(c, c, 0)
Slope(a, b, c) = Slope(promote(convert(Interval, a), b, c)...)

function slope_var(v::Number)
    Slope(v, v, 1)
end

function slope_var(v::Interval, c::Number)
    Slope(v, c, 1)
end

function interval(u::Slope)
    u.x
end

function mid(u::Slope)
    u.c
end

function slope(u::Slope)
    u.s
end

function +(u::Slope, v::Slope)
    Slope(u.x + v.x, u.c + v.c, u.s + v.s)
end

function -(u::Slope, v::Slope)
    Slope(u.x - v.x, u.c - v.c, u.s - v.s)
end

function *(u::Slope, v::Slope)
    Slope(u.x * v.x, u.c * v.c, u.s * v.c + u.x * v.s)
end

function /(u::Slope, v::Slope)
    Slope(u.x / v.x, u.c / v.c, (u.s - (u.c / v.c) * v.s) / v.x)
end

function +(u, v::Slope)
    Slope(u + v.x, u + v.c, v.s)
end

function -(u, v::Slope)
    Slope(u - v.x, u - v.c, -v.s)
end

function *(u, v::Slope)
    Slope(u * v.x, u * v.c, u * v.s)
end

function /(u, v::Slope)
    Slope(u / v.x, u / v.c, -(u / v.c) * (v.s / v.x))
end

+(v::Slope, u) = u + v

*(v::Slope, u) = u * v

function -(u::Slope, v)
    Slope(u.x - v, u.c - v, u.s)
end

function -(u::Slope)
    Slope(-u.x, -u.c, -u.s)
end

function /(u::Slope, v)
    Slope(u.x / v, u.c / v, u.s / v)
end

function sqr(u::Slope)
    Slope(u.x ^ 2, u.c ^ 2, (u.x + u.c) * u.s)
end

function ^(u::Slope, k::Integer)
    if k == 0
        return Slope(1)
    elseif k == 1
        return u
    elseif k == 2
        return sqr(u)
    else
        hxi = interval(u.x.lo) ^ k
        hxs = interval(u.x.hi) ^ k
        hx = hull(hxi, hxs)

        if (k % 2 == 0) && (0 ∈ u.x)
            hx = interval(0, hx.hi)
        end

        hc = u.c ^ k

        i = u.x.lo - u.c.lo
        s = u.x.hi - u.c.hi

        if ((i == 0) || (s == 0) || (k % 2 == 1 && zero(u.x) ⪽ u.x))
            h1 = k * (u.x ^ (k - 1))
        else
            if k % 2 == 0 || u.x.lo >= 0
                h1 = interval((hxi.hi - hc.lo) / i, (hxs.hi - hc.lo) / s)
            else
                h1 = interval((hxs.lo - hc.hi) / s, (hxi.lo - hc.hi) / i)
            end
        end
        return Slope(hx, hc, h1 * u.s)
    end
end

function sqrt(u::Slope)
    Slope(sqrt(u.x), sqrt(u.c), u.s / (sqrt(u.x) + sqrt(u.c)))
end

function exp(u::Slope)
    hx = exp(u.x)
    hc = exp(u.c)

    i = u.x.lo - u.c.lo
    s = u.x.hi - u.c.hi

    if (i == 0 || s == 0)
        h1 = hx
    else
        h1 = interval((hx.lo - hc.lo) / i, (hx.hi - hc.hi) / s)
    end

    Slope(hx, hc, h1 * u.s)
end

function log(u::Slope)
    hx = log(u.x)
    hc = log(u.c)

    i = u.x.lo - u.c.lo
    s = u.x.hi - u.c.hi

    if (i == 0 || s == 0)
        h1 = 1 / u.x
    else
        h1 = interval((hx.hi - hc.hi) / s, (hx.lo - hc.lo) / i)
    end
    Slope(hx, hc, h1 * u.s)
end

function sin(u::Slope) # Using derivative to upper bound the slope expansion for now
    hx = sin(u.x)
    hc = sin(u.c)
    hs = cos(u.x)
    Slope(hx, hc, hs)
end

function cos(u::Slope) # Using derivative to upper bound the slope expansion for now
    hx = cos(u.x)
    hc = cos(u.c)
    hs = -sin(u.x)
    Slope(hx, hc, hs)
end

function tan(u::Slope) # Using derivative to upper bound the slope expansion for now
    hx = tan(u.x)
    hc = tan(u.c)
    hs = (sec(u.x)) ^ 2
    Slope(hx, hc, hs)
end

function asin(u::Slope)
    hx = asin(u.x)
    hc = asin(u.c)
    hs = 1 / sqrt(1 - (u.x ^ 2))
    Slope(hx, hc, hs)
end

function acos(u::Slope)
    hx = acos(u.x)
    hc = acos(u.c)
    hs = -1 / sqrt(1 - (u.x ^ 2))
    Slope(hx, hc, hs)
end

function atan(u::Slope)
    hx = atan(u.x)
    hc = atan(u.c)
    hs = 1 / 1 + (u.x ^ 2)
    Slope(hx, hc, hs)
end


text:
using BlockArrays, FillArrays, LazyArrays, Test
import BlockArrays: SubBlockIterator, BlockIndexRange, Diagonal

@testset "broadcast" begin
    @testset "BlockArray" begin
        A = BlockArray(randn(6), 1:3)

        @test BlockArrays.BroadcastStyle(typeof(A)) == BlockArrays.BlockStyle{1}()
        bc = Base.broadcasted(exp,A)
        @test axes(A) === axes(bc) === axes(similar(bc, Float64))
        @test exp.(A) == exp.(Vector(A))
        @test axes(A) === axes(exp.(A))

        @test A+A isa BlockArray
        @test axes(A + A,1).lasts == axes(A .+ A,1).lasts == axes(A,1).lasts
        @test axes(A .+ 1,1).lasts == axes(A,1).lasts

        A = BlockArray(randn(6,6), 1:3,1:3)

        @test BlockArrays.BroadcastStyle(typeof(A)) == BlockArrays.BlockStyle{2}()

        @test exp.(A) == exp.(Matrix(A))
        @test axes(A) == axes(exp.(A))

        @test axes(A + A) == axes(A .+ A) == axes(A)
        @test axes(A .+ 1) == axes(A)
    end

    @testset "PseudoBlockArray" begin
        A = PseudoBlockArray(randn(6), 1:3)

        @test BlockArrays.BroadcastStyle(typeof(A)) == BlockArrays.PseudoBlockStyle{1}()

        @test exp.(A) == exp.(Vector(A))
        @test axes(A,1).lasts == axes(exp.(A),1).lasts

        @test A+A isa PseudoBlockArray
        @test axes(A + A,1).lasts == axes(A .+ A,1).lasts == axes(A,1).lasts
        @test axes(A .+ 1,1).lasts == axes(A,1).lasts

        B = PseudoBlockArray(randn(6,6), 1:3,1:3)

        @test BlockArrays.BroadcastStyle(typeof(B)) == BlockArrays.PseudoBlockStyle{2}()

        @test exp.(B) == exp.(Matrix(B))
        @test axes(B) == axes(exp.(B))

        @test axes(B + B) == axes(B .+ B) == axes(B)
        @test axes(B .+ 1) == axes(B)
        @test axes(A .+ 1 .+ B) == axes(B)
        @test A .+ 1 .+ B == Vector(A) .+ 1 .+ B == Vector(A) .+ 1 .+ Matrix(B)
    end

    @testset "Mixed" begin
        A = BlockArray(randn(6), 1:3)
        B = PseudoBlockArray(randn(6), 1:3)

        @test A + B isa BlockArray
        @test B + A isa BlockArray

        @test axes(A + B,1).lasts == axes(A,1).lasts

        C = randn(6)

        @test A + C isa BlockVector{Float64}
        @test C + A isa BlockVector{Float64}
        @test B + C isa PseudoBlockVector{Float64}
        @test C + B isa PseudoBlockVector{Float64}

        @test blocksize(A+C) == blocksize(C+A) == blocksize(A)
        @test blocksize(B+C) == blocksize(C+B) == blocksize(B)

        A = BlockArray(randn(6,6), 1:3, 1:3)
        D = Diagonal(ones(6))
        @testset "avoid stack overflow in _generic_block_broadcast" begin
            @test Base.BroadcastStyle(typeof(view(D, Block(1)[1:2], Block(1)[1:2]))) isa Base.Broadcast.DefaultArrayStyle{2}
        end
        @test blocksize(A + D) == blocksize(A)
        @test blocksize(B .+ D) == (3,1)
    end

    @testset "Mixed block sizes" begin
        A = BlockArray(randn(6), 1:3)
        B = BlockArray(randn(6), fill(2,3))
        @test blocksize(A+B) == (5,)

        A = BlockArray(randn(6,6), 1:3, 1:3)
        B = BlockArray(randn(6,6), fill(2,3), fill(3,2))

        @test blocksize(A+B) == (5,3)
    end

    @testset "UnitRange" begin
        n = 3
        x = mortar([1:4n, 1:n])
        @test eltype(x.blocks) <: UnitRange
        y = 1:length(x)
        z = randn(size(x))
        x2 = vcat(x.blocks...)
        y2 = copy(y)
        z2 = copy(z)
        @test (@. z = x + y + z; z) == (@. z2 = x2 + y2 + z2; z2)
    end

    @testset "Special broadcast" begin
        v = mortar([1:3,4:7])
        @test broadcast(+, v) isa BlockVector{Int,Vector{UnitRange{Int}}}
        @test broadcast(+, v) == v
        @test broadcast(-, v) isa BlockVector{Int,Vector{StepRange{Int,Int}}}
        @test broadcast(-, v) == -v == -Vector(v)
        @test broadcast(+, v, 1) isa BlockVector{Int,Vector{UnitRange{Int}}}
        @test broadcast(+, v, 1) == Vector(v).+1
        @test broadcast(*, 2, v) isa BlockVector{Int,Vector{StepRange{Int,Int}}}
        @test broadcast(*, 2, v) == 2Vector(v)
    end

    @testset "BlockVector broadcast" begin
        a = BlockArray(randn(3),[1,2])
        b = broadcast(+, Ref(1), a)
        @test b isa BlockArray{eltype(a),1}
        @test Base.axes1(a) == Base.axes1(b)
        @test Vector(b) == broadcast(+, Ref(1), Vector(a))
        a = BlockArray(randn(16),[1,3,5,7])
        b = broadcast(+, Ref(1), a)
        @test b isa BlockArray{eltype(a),1}
        @test Base.axes1(a) == Base.axes1(b)
        @test Vector(b) == broadcast(+, Ref(1), Vector(a))
    end

    @testset "special axes" begin
        A = BlockArray(randn(6), Ones{Int}(6))
        B = BlockArray(randn(6), Ones{Int}(6))
        @test axes(A+B,1) === axes(A,1)

        C = BlockArray(randn(6), (BlockArrays._BlockedUnitRange(1,2:6),))
        @test axes(A+C,1) === BlockArrays._BlockedUnitRange(1,1:6)
    end

    @testset "Views" begin
        A = BlockArray(randn(6), 1:3)
        @test Base.BroadcastStyle(typeof(view(A, Block(2)))) isa Base.Broadcast.DefaultArrayStyle{1}
        V = view(A, Block.(2:3))
        @test Base.BroadcastStyle(typeof(V)) isa BlockArrays.BlockStyle{1}
        @test V .+ 1 isa BlockArray
        @test 1 .+ V isa BlockArray
        @test V .+ 1 == 1 .+ V == A[Block.(2:3)] .+ 1
        @test -V isa BlockArray
        @test -V == -A[Block.(2:3)]
    end

    @testset "Fill broadcast" begin
        a = BlockArray(randn(6), Ones{Int}(6))
        @test a .+ a == Vector(a) .+ a == Vector(a) .+ Vector(a)
        @test axes(a) ≡ axes(a .+ a)

        @test a .+ (1:6) == (1:6) .+ a

        b = BlockArray(randn(6), (BlockArrays._BlockedUnitRange(1, 2:6),))
        @test b .+ b == Vector(b) .+ b == Vector(b) .+ Vector(b)
        @test axes(b) ≡ axes(b .+ b)
        @test axes(a .+ b,1) ≡ BlockArrays._BlockedUnitRange(1, 1:6)

        @test a .+ PseudoBlockArray(b) == PseudoBlockArray(a) .+ b

        A = BlockArray(randn(6), 1:3)
        @test blockisequal(axes(A .* Zeros(6)), axes(A .* zeros(6)))
        @test blockisequal(axes(A .* Ones(6)), axes(A .* ones(6)))
        @test blockisequal(axes(A .* Zeros(axes(A))), axes(Zeros(axes(A)) .* A), axes(A .* zeros(6)))
        @test blockisequal(axes(A .* Ones(axes(A))), axes(Ones(axes(A)) .* A), axes(A .* ones(6)))
    end

    @testset "recurrences" begin
        N = 1000
        n = mortar(BroadcastArray(Fill,Base.OneTo(N),Base.OneTo(N)))
        k = mortar(BroadcastArray(Base.OneTo,Base.OneTo(N)))

        @test view(n, Block(5)) ≡ Fill(5,5)
        @test view(k,Block(5)) ≡ Base.OneTo(5)
        a = b = c = 0.0
        # for some reason the following was causing major slowdown. I think it
        # went pass a limit to Base.Broadcast.flatten which caused `bc.f` to have a strange type.
        # bc = Base.Broadcast.instantiate(Base.broadcasted(/, Base.broadcasted(*, k, Base.broadcasted(-, Base.broadcasted(-, k, n), a)), Base.broadcasted(+, 2k, b+c-1)))

        bc = Base.Broadcast.instantiate(Base.broadcasted((k,n,a,b,c) -> k * (k-n-a) / (2k+(b+c-1)), k, n, a, b, c))
        @test axes(n,1) ≡ axes(k,1) ≡ axes(bc)[1] ≡ blockedrange(Base.OneTo(N))
        u = (k .* (k .- n .- a) ./ (2k .+ (b+c-1)))
        @test u == (Vector(k) .* (Vector(k) .- Vector(n) .- a) ./ (2Vector(k) .+ (b+c-1)))
        @test copyto!(u, bc) == (k .* (k .- n .- a) ./ (2k .+ (b+c-1)))
        @test @allocated(copyto!(u, bc)) ≤ 1000
        # not clear why allocatinos so high: all allocations are coming from checking
        # axes

        u = PseudoBlockArray{Float64}(undef, collect(1:N))
        @test copyto!(u, bc) == (k .* (k .- n .- a) ./ (2k .+ (b+c-1)))
        @test @allocated(copyto!(u, bc)) ≤ 1000
    end

    @testset "type inferrence" begin
        u = BlockArray(randn(5), [2,3]);
        @inferred(copyto!(similar(u), Base.broadcasted(exp, u)))
        @test exp.(u) == exp.(Vector(u))
    end

    @testset "adjtrans" begin
        a = PseudoBlockArray(randn(6), [2,3])
        b = BlockArray(a)

        @test Base.BroadcastStyle(typeof(a')) isa BlockArrays.PseudoBlockStyle{2}
        @test Base.BroadcastStyle(typeof(b')) isa BlockArrays.BlockStyle{2}

        @test exp.(a') == exp.(b') == exp.(Vector(a)')
        @test exp.(a') isa PseudoBlockArray
        @test exp.(b') isa BlockArray
    end

    @testset "subarray" begin
        @testset "vector" begin
            a = PseudoBlockArray(randn(6), [2,3])
            b = BlockArray(a)

            v = view(a,Block.(1:2))
            w = view(b,Block.(1:2))

            @test Base.BroadcastStyle(typeof(v)) isa BlockArrays.PseudoBlockStyle{1}
            @test Base.BroadcastStyle(typeof(w)) isa BlockArrays.BlockStyle{1}

            @test exp.(v) == exp.(w) == exp.(Vector(v))
            @test exp.(v) isa PseudoBlockArray
            @test exp.(w) isa BlockArray
        end

        @testset "matrix" begin
            A = PseudoBlockArray(randn(6,3), [2,3],[1,2])
            B = BlockArray(A)

            v = view(A,1:3,Block.(1:2))
            w = view(B,1:3,Block.(1:2))

            @test Base.BroadcastStyle(typeof(v)) isa BlockArrays.PseudoBlockStyle{2}
            @test Base.BroadcastStyle(typeof(w)) isa BlockArrays.BlockStyle{2}

            @test exp.(v) == exp.(w) == exp.(Matrix(v))
            @test exp.(v) isa PseudoBlockArray
            @test exp.(w) isa BlockArray

            @test Base.BroadcastStyle(typeof(view(A,Block.(1:2),Block.(1:2)))) isa BlockArrays.PseudoBlockStyle{2}
            @test Base.BroadcastStyle(typeof(view(B,Block.(1:2),Block.(1:2)))) isa BlockArrays.BlockStyle{2}
        end
    end

    @testset "SubBlockIterator" begin
        A = BlockArray(1:6, 1:3);
        subblock_lasts = axes(A, 1).lasts
        @test subblock_lasts == [1, 3, 6]
        block_lasts = [1, 3, 4, 6]

        for idx in SubBlockIterator(subblock_lasts, block_lasts)
           B = view(A, idx)
           @test !(parent(B) isa BlockArray)
           @test idx isa BlockIndexRange
           @test idx.block isa Block{1}
           @test idx.indices isa Tuple{UnitRange}
        end

        @test [idx.block.n[1] for idx in SubBlockIterator(subblock_lasts, block_lasts)] == [1,2,3,3]

        @test [idx.indices[1] for idx in SubBlockIterator(subblock_lasts, block_lasts)] == [1:1,1:2,1:1,2:3]
    end
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 IoK8sApiCoreV1StorageOSPersistentVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectReference } # spec name: secretRef
    volumeName::Any # spec type: Union{ Nothing, String } # spec name: volumeName
    volumeNamespace::Any # spec type: Union{ Nothing, String } # spec name: volumeNamespace

    function IoK8sApiCoreV1StorageOSPersistentVolumeSource(;fsType=nothing, readOnly=nothing, secretRef=nothing, volumeName=nothing, volumeNamespace=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1StorageOSPersistentVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1StorageOSPersistentVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1StorageOSPersistentVolumeSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        validate_property(IoK8sApiCoreV1StorageOSPersistentVolumeSource, Symbol("volumeName"), volumeName)
        setfield!(o, Symbol("volumeName"), volumeName)
        validate_property(IoK8sApiCoreV1StorageOSPersistentVolumeSource, Symbol("volumeNamespace"), volumeNamespace)
        setfield!(o, Symbol("volumeNamespace"), volumeNamespace)
        o
    end
end # type IoK8sApiCoreV1StorageOSPersistentVolumeSource

const _property_map_IoK8sApiCoreV1StorageOSPersistentVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretRef")=>Symbol("secretRef"), Symbol("volumeName")=>Symbol("volumeName"), Symbol("volumeNamespace")=>Symbol("volumeNamespace"))
const _property_types_IoK8sApiCoreV1StorageOSPersistentVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("readOnly")=>"Bool", Symbol("secretRef")=>"IoK8sApiCoreV1ObjectReference", Symbol("volumeName")=>"String", Symbol("volumeNamespace")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1StorageOSPersistentVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1StorageOSPersistentVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1StorageOSPersistentVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1StorageOSPersistentVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1StorageOSPersistentVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1StorageOSPersistentVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1StorageOSPersistentVolumeSource)
    true
end

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


text:
floyd(n) =
  pprint([join([lpad(j+binomial(i,2), (j==1?0:1)+ndigits(j+binomial(n,2)), " ")
               for j=1:i])
         for i=1:n])

pprint(matrix) = for i = 1:size(matrix,1) println(join(matrix[i,:])) end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct TendencyIntegralModel{M} <: AbstractOceanModel
    ocean::M
    function TendencyIntegralModel(ocean::M) where {M}
        return new{M}(ocean)
    end
end
vars_state(tm::TendencyIntegralModel, ::Prognostic, FT) =
    vars_state(tm.ocean, Prognostic(), FT)
vars_state(tm::TendencyIntegralModel, ::GradientFlux, FT) = @vars()

function vars_state(m::TendencyIntegralModel, ::Auxiliary, T)
    @vars begin
        ∫du::SVector{2, T}
    end
end

function vars_state(m::TendencyIntegralModel, ::UpwardIntegrals, T)
    @vars begin
        ∫du::SVector{2, T}
    end
end

@inline function integral_load_auxiliary_state!(
    m::TendencyIntegralModel,
    I::Vars,
    Q::Vars,
    A::Vars,
)
    I.∫du = A.∫du

    return nothing
end

@inline function integral_set_auxiliary_state!(
    m::TendencyIntegralModel,
    A::Vars,
    I::Vars,
)
    A.∫du = I.∫du

    return nothing
end

init_state_auxiliary!(tm::TendencyIntegralModel, A::Vars, _...) = nothing

function update_auxiliary_state!(
    dg::DGModel,
    tm::TendencyIntegralModel,
    dQ::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    A = dg.state_auxiliary

    # copy tendency vector to aux state for integration
    function f!(::TendencyIntegralModel, dQ, A, t)
        @inbounds begin
            A.∫du = @SVector [dQ.u[1], dQ.u[2]]
        end

        return nothing
    end
    update_auxiliary_state!(f!, dg, tm, dQ, t)

    # compute integral for Gᵁ
    indefinite_stack_integral!(dg, tm, dQ, A, t, elems) # bottom -> top

    return true
end

#-------------------------------------------------------------------------------
struct FlowIntegralModel{M} <: AbstractOceanModel
    ocean::M
    function FlowIntegralModel(ocean::M) where {M}
        return new{M}(ocean)
    end
end
vars_state(fm::FlowIntegralModel, ::Prognostic, FT) =
    vars_state(fm.ocean, Prognostic(), FT)
vars_state(fm::FlowIntegralModel, ::GradientFlux, FT) = @vars()

function vars_state(m::FlowIntegralModel, ::Auxiliary, T)
    @vars begin
        ∫u::SVector{2, T}
    end
end

function vars_state(m::FlowIntegralModel, ::UpwardIntegrals, T)
    @vars begin
        ∫u::SVector{2, T}
    end
end

@inline function integral_load_auxiliary_state!(
    m::FlowIntegralModel,
    I::Vars,
    Q::Vars,
    A::Vars,
)
    I.∫u = Q.u

    return nothing
end

@inline function integral_set_auxiliary_state!(
    m::FlowIntegralModel,
    A::Vars,
    I::Vars,
)
    A.∫u = I.∫u

    return nothing
end

init_state_auxiliary!(fm::FlowIntegralModel, A::Vars, _...) = nothing

function update_auxiliary_state!(
    dg::DGModel,
    fm::FlowIntegralModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    A = dg.state_auxiliary

    # compute vertical integral of u
    indefinite_stack_integral!(dg, fm, Q, A, t, elems) # bottom -> top

    return true
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 09bee220-1249-11ec-3e2a-4fdde9976a59
begin
	import Pkg
    Pkg.activate(mktempdir())
	Pkg.add(url="https://github.com/brendanjohnharris/Catch22.jl#main")
	using Catch22
	Pkg.add(url="https://github.com/tk3369/YeoJohnsonTrans.jl")
	Pkg.develop(path="C:\\Users\\Brendan\\OneDrive - The University of Sydney (Students)\\Honours\\Code\\ParameterInference.jl")
	Pkg.add(path="C:\\Users\\Brendan\\OneDrive - The University of Sydney (Students)\\Honours\\Code\\ParameterInference.jl")
	using ParameterInference, Plots, StatsBase, Statistics, LinearAlgebra, ParameterInference.NonstationaryProcesses, MultivariateStats
	gr(); fourseas!(); # Plots initialisation
	nothing
end

# ╔═╡ 1cbf17f7-8322-4e3c-a895-163a48332205
md"# Baseline PCA"

# ╔═╡ 2f0cb4e5-4424-42d9-8eb7-e8bf128575ad
md"""
We will construct a dataset that demonstrates the impact of dependencies between features on PCA estimates.
"""

# ╔═╡ 0f10dc15-9153-4d8c-a5e5-c0318c459fa1
md"## Parameters"

# ╔═╡ 230c66e5-7898-48af-9544-439677bee137
md"For a simple example, choose to have two 'parameters' that drive our manually constructed 'features': 𝛼 and 𝛽. Both will vary as random Gaussians.:"

# ╔═╡ fee5a67f-1678-45fc-83b3-3ffc187c19b3
nₛ = 1000;

# ╔═╡ c3e707a4-1e35-4fe0-8e6d-88d0eca27fde
𝜂 = 0.25;

# ╔═╡ 5146ae95-d551-49ec-888a-71fc7bf65d14
𝛼, 𝛽 = randn(nₛ, 2) |> eachcol;

# ╔═╡ b646056a-7c34-4cb3-aac9-d255646f6b75
plotcov(x) = covarianceimage(1:5, x, colormode=:raw, docluster=false, colorbar_title="", title="\$|\\Sigma^2|\$");

# ╔═╡ 2fbff3c4-9377-4624-90d6-d4ea7c8da64e
md"# Feature scenario #1"

# ╔═╡ 404bb5c9-74ab-43a9-8791-67ac23b0c0f8
md"""Next, we construct 5 features according to the parameters.
In our first scenario, these are:

| Feature      	| Description  	|
| ----------- 	| ----------- 	|
| 𝑓₁ -- 𝑓₄      | Sensitive to 𝛼, but intercorrelated  |
| 𝑓₅   		 | Sensitive to 𝛽  |
"""

# ╔═╡ 4f01861c-dd52-462f-b340-ceecd3385c8e
md"""### High-dimensional baseline
In which both 𝛼 and 𝛽 are varying:
"""

# ╔═╡ 92670bd1-ec33-4590-9a53-ceca6f25ac0f
𝑓ₕ₁ = hcat(𝛼, 𝛼, 𝛼, 𝛼, 𝛽)';

# ╔═╡ 8629c505-a1e4-4c73-b2d4-44ecc858373d
plotcov(𝑓ₕ₁)

# ╔═╡ c3c3d6c3-49a0-40a6-b41b-dbdfa20a7d6b
md"""### Test data
In which 𝛼 is constrained:
"""

# ╔═╡ 9247ba48-f242-4479-a92d-e994b1610279
𝑓₁ = hcat(𝜂.*[𝛼, 𝛼, 𝛼, 𝛼]..., 𝛽)';

# ╔═╡ c973721b-c20b-4573-a8fa-a9a9282fc87d
plotcov(𝑓₁)

# ╔═╡ 4032a036-842b-462b-b75d-ebaaaf4aded9
md"### True parameters"

# ╔═╡ 6a957b7d-8fff-46c5-b441-01f8b9d45eb7
scatter(𝛽, 𝜂.*𝛼, aspect_ratio=:equal, ylabel="𝛼", xlabel="𝛽", title="True parameter variation",  markersize=2, markercolor=:gray)

# ╔═╡ 5686638b-ba3f-4691-9291-d91dfa510aa7
md"### PCA"

# ╔═╡ 76e702f7-521c-4445-aa6c-de28215d8d44
𝛴² = cov(𝑓₁');

# ╔═╡ b17a0c26-ce25-4e11-84dd-0d4b86371ae5
𝑃 = eigvecs(𝛴²);

# ╔═╡ 69e0799e-18f4-4f19-9bae-44202c2bba5d
𝑓′₁ = 𝑃'*𝑓₁;

# ╔═╡ cf454ccf-e505-4df9-a3cf-4afba2211d4a
scatter(𝑓′₁[5, :], 𝑓′₁[4, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′", title="PCA estimate"); # Scatter significant PCs

# ╔═╡ 8013d3a2-a011-456b-9bfe-3ec781c31aba
scatter!(𝛽, 𝜂.*𝛼, aspect_ratio=:equal, markersize=2, markercolor=:gray)

# ╔═╡ ce024bfd-b892-4a31-a0d1-1d4befc7b80e
md"### Standardised PCA"

# ╔═╡ 542aec21-e623-4e53-ae4b-a7f240bfcb08
𝑓̂₁ = 𝑓₁./std(𝑓₁, dims=2);

# ╔═╡ 1daa139a-23c0-4995-bf73-abfbb7c518dd
𝛴̂² = cov(𝑓̂₁');

# ╔═╡ 705068d7-033f-49d2-bb59-4b40aa514642
𝑃̂ =  fit(MultivariateStats.PCA, 𝑓̂₁).proj; #eigvecs(𝛴̂²);

# ╔═╡ f1a685e5-ab99-4a61-a2d9-41ebf33372b0
𝑓̂′₁ = 𝑃̂'*𝑓₁;

# ╔═╡ 3a7af0db-8e3d-4100-a7c2-dcd37253a9ce
scatter(𝑓̂′₁[2, :], 𝑓̂′₁[1, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′", title="Standardised PCA estimate");

# ╔═╡ 7a28c225-bab8-40ee-906e-0c5349accc84
scatter!(𝛽, 𝜂.*𝛼, aspect_ratio=:equal, markersize=2, markercolor=:gray)

# ╔═╡ 4d619965-6eb7-4528-804c-00bbf1572e60
md"""### Baseline PCA
In this case, whitening:"""

# ╔═╡ 7ddd5530-7590-42a7-8e66-54432f717eec
𝛴²ₕ = cov(𝑓ₕ₁');

# ╔═╡ a95a3d4c-cb8b-4840-bc67-57fb2be00065
m1 = fit(MultivariateStats.PCA, 𝑓ₕ₁; maxoutdim=5, pratio=1.0)

# ╔═╡ 688dc093-c173-4eb4-ad12-4ff2e899b9f9
𝑃ₕ = m1.proj; # eigvecs(𝛴²ₕ)

# ╔═╡ f4081452-9314-47a4-948d-2813bda53507
𝛬ₕ = sqrt(diagm(m1.prinvars)); # Diagonal(sqrt.(abs.(eigvals(𝛴²ₕ))));

# ╔═╡ 5b853331-88c7-4b7f-b26e-8e0f963a5ed5
𝑓′′₁ = (𝑃ₕ*inv(𝛬ₕ))'*𝑓₁;

# ╔═╡ 1ee56732-671e-447e-adb8-cce45902b935
scatter(𝑓′′₁[2, :], .-𝑓′′₁[1, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′");

# ╔═╡ 70d2fe5f-7f7d-4891-8490-5c508c21c6f6
scatter!(𝛽, 𝜂.*𝛼, aspect_ratio=:equal, markersize=2, markercolor=:gray)

# ╔═╡ bb52a2b3-e856-4785-abb9-78ba15382d11
md"# Feature scenario #2"

# ╔═╡ 70d2f196-1604-4cf4-8e38-05dd8c099c83
md"""What if features one to four are sensitive to both 𝛼 and 𝛽? Then the effect of constraining 𝛼 manifests as an increase in correlation between features. For our second scenario, the features are then:

| Feature      	| Description  	|
| ----------- 	| ----------- 	|
| 𝑓₁ -- 𝑓₄      | Sensitive to 𝛼 and 𝛽, but intercorrelated  |
| 𝑓₅   		 | Sensitive to 𝛽  |
"""

# ╔═╡ 009c662a-9940-43ef-b0c9-d55560e6f49b
md"### High-dimensional baseline"

# ╔═╡ 46e7d2a0-8b62-4112-92df-74538f5f16d8
𝑓ₕ₂ = hcat(𝛼+𝛽/4, 𝛼+𝛽/4, 𝛼+𝛽/4, 𝛼+𝛽/4, 𝛽)';

# ╔═╡ 595e5022-9099-45ab-b5b0-b92fedb96ba8
plotcov(𝑓ₕ₂./std(𝑓ₕ₂, dims=2))

# ╔═╡ 0d062a65-adca-42c4-8b1c-84095d5335c5
md"### Test data"

# ╔═╡ 0f6b9065-55ad-4b1e-bcc3-2e177a0e331d
𝑓₂ = hcat((𝜂.*[𝛼, 𝛼, 𝛼, 𝛼] .+ [𝛽, 𝛽, 𝛽, 𝛽]./4)..., 𝛽)';

# ╔═╡ fd506ba8-9a5b-4235-bc43-b730522dce4f
plotcov(𝑓₂./std(𝑓₂, dims=2))

# ╔═╡ 58083319-d649-4c2d-b4b2-c2df7c49ee46
md"### True parameters"

# ╔═╡ 2a9a4929-8755-4237-a454-fda125d548ca
scatter(𝛽, 𝜂.*𝛼 .+ 𝛽./4, aspect_ratio=:equal, ylabel="𝛼", xlabel="𝛽", title="True parameter variation",  markersize=2, markercolor=:gray)

# ╔═╡ 1762e05f-7eec-4411-865a-eba47531a49c
md"### Standardised PCA"

# ╔═╡ feeaa835-28e0-4a24-905b-53347560ec9a
𝑓̂₂ = 𝑓₂./std(𝑓₂, dims=2);

# ╔═╡ 2358f79f-f44d-4e93-b4da-75e7be715882
𝛴̂²₂ = cov(𝑓̂₂');

# ╔═╡ d9fca85a-63b6-4ab5-ab0e-6cd5e68c4e06
𝑃̂₂ = fit(MultivariateStats.PCA, 𝑓̂₂).proj; # eigvecs(𝛴̂²₂);

# ╔═╡ d700b4ae-8bd5-4437-b515-c21d3da70543
𝑓̂′₂ = 𝑃̂₂'*𝑓₂;

# ╔═╡ 06e226a6-0f0f-4ef6-b916-04d51a3789c8
scatter(𝑓̂′₂[1, :], .-𝑓̂′₂[2, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′", title="Standardised PCA estimate");

# ╔═╡ ea513127-e9bb-40a7-9ec0-2beee8fa6179
scatter!(𝛽, 𝜂.*𝛼 .+ 𝛽./4, markersize=2, markercolor=:gray)

# ╔═╡ f0a7f529-19ac-44f7-8ab1-8d72d596a0ef
md"""### Baseline PCA"""

# ╔═╡ f946e3d7-15b9-45a9-951b-b3d77fc589f8
𝛴²ₕ₂ = cov(𝑓ₕ₂');

# ╔═╡ 652b998d-a6b0-4bb7-908e-c3dc7a82764c
m2 = fit(MultivariateStats.PCA, 𝑓ₕ₂);

# ╔═╡ d192149c-1ef4-4761-9e78-879a97abee57
𝑃ₕ₂ = m2.proj; #eigvecs(𝛴²ₕ₂);

# ╔═╡ 14d5e08c-0808-4352-835e-c173c53fd8d7
𝛬ₕ₂ = sqrt(diagm(m2.prinvars)); # Diagonal(sqrt.(abs.(eigvals(𝛴²ₕ₂))))

# ╔═╡ 772d260a-72b3-4c9e-8e65-6368508cfb52
𝑓′′₂ = (𝑃ₕ₂*inv(𝛬ₕ₂))'*𝑓₂;

# ╔═╡ 0f83a03c-890d-4771-b83f-16d8d5db81e9
cov(𝑓′′₂,dims=2)

# ╔═╡ b9a8e7e4-5f9b-462d-845e-b873f50bfc90
scatter(𝑓′′₂[2, :], .-𝑓′′₂[1, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′");

# ╔═╡ f1df27ea-0537-4f4b-afe7-c0692ca7b26d
scatter!(𝛽, 𝜂.*𝛼 .+ 𝛽./4, markersize=2, markercolor=:gray)

# ╔═╡ 2f987bde-a87a-43a8-8653-0c5812c56a70
md"# Feature scenario #3"

# ╔═╡ 1d18e8c7-7eff-4b75-ad67-e61bc22e12ea
md"""What if each feature has some intrinsic noise? This requires incorporating low-dimensional variance into the estimate.
To show this:

| Feature      	| Description  	|
| ----------- 	| ----------- 	|
| 𝑓₁ -- 𝑓₃      | Sensitive to 𝛼, but intercorrelated and with small noise |
| 𝑓₄   		 | Sensitive to 𝛽, but with small noise  |
| 𝑓₅   		 | Small noise only  |
"""

# ╔═╡ 3821aa68-e3b3-4764-a7ff-0b6b0e6bdc3e
𝜂₃ = 0.25;

# ╔═╡ 9d0c3f2b-8cf4-41f9-8fe8-5cd9a19dbd88
𝑓ₕ₃ = hcat(repeat([𝛼.+𝜂₃.*randn(nₛ)],3)..., 𝛽.+𝜂₃.*randn(nₛ), (𝜂₃.-0.01).*randn(nₛ))';

# ╔═╡ b0d86eb8-6f94-424a-adb2-eae8e40be1fd
plotcov(𝑓ₕ₃)

# ╔═╡ e507beed-984e-4f8b-a59b-bb3dcff1416a
md"""### Zero-dimensional baseline
In which no parameters are varying:
"""

# ╔═╡ ea23a142-1c36-4430-a598-2fcf0f2e3b66
𝑓₀ = hcat(repeat([𝜂₃.*randn(nₛ)], 3)..., 𝜂₃.*randn(nₛ), 𝜂₃.*randn(nₛ))';

# ╔═╡ fe0902aa-b55e-4b72-a35e-77db673747c9
plotcov(𝑓₀)

# ╔═╡ 9b92548e-ae54-48de-a5c9-4be88ae11228
md"### Test data"

# ╔═╡ 33167762-e2c6-4aba-98b7-58e28ee35ccb
𝑓₃ = hcat(repeat([0.25.*𝛼 .+ 𝜂₃.*randn(nₛ)], 3)..., 𝛽.+𝜂₃.*randn(nₛ), 𝜂₃.*randn(nₛ))';

# ╔═╡ 35dcd411-a342-4ee2-9146-f45ad39d6341
plotcov(𝑓₃)

# ╔═╡ d2404839-53bf-49d3-a9bf-a0ea81471c43
md"### True parameters"

# ╔═╡ e64a04de-9257-4bbb-84c4-007c532d9c72
scatter(𝛽, 0.25.*𝛼, aspect_ratio=:equal, ylabel="𝛼", xlabel="𝛽", title="True parameter variation",  markersize=2, markercolor=:gray)

# ╔═╡ 710d51ea-8196-4a7f-91fe-1a9bd3f06d01
md"""
### Standardised PCA
"""

# ╔═╡ a4c8677f-6bdc-4bfe-adcc-3130f45b7aaa
𝑓̂₃ = 𝑓₃./std(𝑓₃, dims=2);

# ╔═╡ 5ad01546-7daa-441f-bf2c-60e2a6e6f0c7
𝛴̂²₃ = cov(𝑓̂₃');

# ╔═╡ 36039419-6df3-45f2-ae2c-77cf6e3b405d
𝑃̂₃ = eigvecs(𝛴̂²₃);

# ╔═╡ d046eb57-d4cb-4eda-9a81-00a234af1d57
𝑓̂′₃ = 𝑃̂₃'*𝑓₃;

# ╔═╡ 4e9b0004-da7a-4cb3-afe9-c3365d37471a
scatter(𝑓̂′₃[4, :], 𝑓̂′₃[5, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′", title="Standardised PCA estimate");

# ╔═╡ 5e2281e9-360b-41dd-9f4f-38e14d8aea86
scatter!(𝛽, 0.25.*𝛼, markersize=2, markercolor=:gray)

# ╔═╡ fa43e761-0b22-4903-aa6c-47050d8c27c0
scatter(𝛼, 𝑓̂′₃[5, :])

# ╔═╡ bcd85e11-8bf0-4968-a528-8b56ce14db5b
md"""### Baseline PCA"""

# ╔═╡ 3fe79e11-947e-46e6-a5ef-1bf1b537c857
md"""
The full baseline-corrected covariance is:

$$\Sigma^{\prime 2} = (P \Lambda_h^{-1})^T \left(\Sigma^2 - \Sigma_0^2\right)(P \Lambda_h^{-1}).$$

But this is not a linear transformation of feature space, so some approximations must (?) be made.

An easy one is to treat the zero-dimensional covariance as diagonal. In that case:

$$\Sigma^2 - \Sigma_0^2 \approx V^T \Sigma^2 V,$$
where:

$$V = \sqrt{I - D(\Sigma^2_0) D(\Sigma^{-2})},$$
And $D(A)$ is the diagonal of $A$. The resulting approximation is simply a linear rescaling. Importantly, it satisfies the condition that $V \to 0$ as $\Sigma^2 \to \Sigma^2_0$.

There are other approximations that also allow covariances to vanish in the same limit...
"""

# ╔═╡ 15970d81-af01-443e-b54a-bda505b1b3be
md"""
For now, we will just look at a heuristic approximation:

$$V = I - D(\Sigma') D(\Sigma_0')$$

which gives a total transformation of:

$$F^\prime = FP\frac{I - D(\Sigma^\prime)^{-1}D(\Sigma^\prime_0)}{\Lambda_h - D(\Sigma^\prime_0)}$$

where $P$ is the projection matrix of the high-dimensional baseline and the primes indicate application of $P$.

"""

# ╔═╡ f771887a-77f0-4a4d-9e4c-7e988201197f
𝛴²ₕ₃ = cov(𝑓ₕ₃');

# ╔═╡ b636df98-f935-46b4-bb08-be96cee18884
𝛴²₀ = cov(𝑓₀');

# ╔═╡ c43d6053-4eeb-450a-bca4-bd5d58b179c1
𝛴²₃ = cov(𝑓₃');

# ╔═╡ beafaf38-3753-4e9f-bf6d-69b98439fe18
m3 = fit(MultivariateStats.PCA, 𝑓ₕ₃)

# ╔═╡ fd00b9e4-f2fc-4571-8596-b1e53f8d7286
𝑃ₕ₃ = m3.proj;#eigvecs(𝛴²ₕ₃);

# ╔═╡ 9d13988d-e625-42b5-8e89-86d519027912
𝛬ₕ₃ = sqrt(diagm(m3.prinvars)); #Diagonal(sqrt.(abs.(eigvals(𝛴²ₕ₃))))

# ╔═╡ 4f5c167d-dbde-4fed-bd2d-63b0100e2d1a
𝛴²ₕ′ = diagm(m3.prinvars)

# ╔═╡ 82841686-b078-4ba2-b5ce-d748070437d2
𝛴²₀′ = cov(𝑃ₕ₃'*𝑓₀, dims=2)

# ╔═╡ 8dded5c7-f4ea-4ee5-b0b9-59b21c7f3097
𝛴²′ = cov(𝑃ₕ₃'*𝑓₃, dims=2)

# ╔═╡ cb77a71b-7cba-489d-a51b-690cdd281a4c
𝑉 = I - inv(sqrt(Diagonal(𝛴²′)))*sqrt(Diagonal(𝛴²₀′));

# ╔═╡ b00f965a-18a1-4ccf-85a9-ae79d588cab6
𝑆 = 𝑉/(𝛬ₕ₃ - sqrt(Diagonal(𝛴²₀′)))

# ╔═╡ 0ae23619-0ce3-4e32-9a11-f9a444235c70
𝑓′′₃ = 𝑆'*𝑃ₕ₃'*𝑓₃;

# ╔═╡ ddb9528f-b57c-4333-9347-8a56aa6890c4
#f′′₃ = (𝑓′′₃./std(𝑓′′₃, dims=2))*((sqrt(Diagonal(𝛴²′)) - sqrt(Diagonal(𝛴²₀′)))/())

# ╔═╡ c3067e6d-dfaf-4938-8bbf-4b51514fd286
𝑓′′₀ = (𝑃ₕ₃*inv(𝛬ₕ₃)*𝑉)'*𝛴²₀*(𝑃ₕ₃*inv(𝛬ₕ₃)*𝑉); # Not 0 because the transform is a strong approximation

# ╔═╡ a2f4c245-031b-4f1c-9f44-b3947cbc9795
cov(𝑓′′₃, dims=2)

# ╔═╡ ad7aac3b-303a-4ca6-a44f-f11cba8d115f
𝑓′′₀

# ╔═╡ dbdcd6c4-102d-41cc-b218-100f74355b85
scatter(𝑓′′₃[2, :], .-𝑓′′₃[1, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′", title="Baseline PCA estimate", left_margin=5Plots.mm);

# ╔═╡ 9d6334c1-7107-4fbb-b383-9a37fc56de74
scatter!(.-𝛽, .-0.25.*𝛼, markersize=1, markercolor=:gray)

# ╔═╡ 07e9a417-5adb-4182-a98d-fd618feaf675
scatter(𝛽, 𝑓′′₃[2, :]) # Compairwise to true parameters

# ╔═╡ 008512c1-3a44-42b4-8695-2a9f4d88ef01
scatter(𝛼, 𝑓′′₃[1, :]) # Compairwise to true parameters

# ╔═╡ 2826e49a-5566-4aa7-be9a-73799dfbebc7
scatter(𝑓̂′₃[4, :], 𝑓′′₃[4, :]); # Compare to standardised PCA estimate. Similar, which is good, but properly scaled relative to other PCs (as above).

# ╔═╡ c79d0914-d0db-4f5c-b1f2-a6154d063e78

# ╔═╡ eaae76e1-9326-432a-8fc5-b03c3d93e75a

# ╔═╡ b7a9e831-6ab4-4afd-a5da-9d3b5414245b
fₕ₃ = 𝑓ₕ₃ .+eps().*randn(size(𝑓ₕ₃))

# ╔═╡ 70109d54-cfcd-46c5-8c67-a3b96286ee09
f₃ = 𝑓₃ .+eps().*randn(size(𝑓₃))

# ╔═╡ 9970c6e2-a2c9-442e-8fca-17b8a5c360b9
𝑃ₕ₄ = eigvecs(𝛴²ₕ₃)

# ╔═╡ 7408597f-8bc6-493a-80f6-79a5201d625b
#Λ = Diagonal(sqrt.(abs.(eigvals(𝛴²ₕ₃))))

# ╔═╡ a031cdd3-7e92-45d8-8037-e6695cd588c2
MultivariateStats.PCA

# ╔═╡ 816e841b-3f3b-4b8c-8c4b-4cd93c685ab7
m = fit(MultivariateStats.PCA, 𝑓ₕ₃; maxoutdim=5, pratio=1.0)

# ╔═╡ 923b9503-5a8e-4dab-984d-6ca07d625bb4
Σₕ = cov(fₕ₃, dims=2)

# ╔═╡ 9938d1b8-8f80-45b0-bdff-48cd74e78fa4
P = m.proj

# ╔═╡ 4771733f-83a5-4c13-b0ea-cb54241f7e19
Λ² = diagm(eigvals(Σₕ))

# ╔═╡ d8847644-4d5b-4113-8870-e8c3dba7e1a2
Λ²[Λ².<0] .= eps()

# ╔═╡ 5463ac9d-9bf7-4363-abb9-d1de01a88681
Λ = sqrt(diagm(m.prinvars))

# ╔═╡ 39871429-722a-4745-9a67-62c2961fb268
ffh = P'*𝑓ₕ₃

# ╔═╡ 31a016af-82fa-41a1-b2e8-8adea326c11c
ff = P'*𝑓₃

# ╔═╡ 31cfab8a-e219-484c-9428-1fef8ea16ca5
ff0 = P'*𝑓₀

# ╔═╡ c64a1c15-824a-49d2-ac0f-10573a1472c5
l0 = cov(ff0, dims=2)

# ╔═╡ 8c9e6bd3-c548-4126-a6fb-a100e91f8747
l = cov(ff, dims=2)

# ╔═╡ 24425756-1e29-4e78-9398-6ef1820940b8
lh = Λ^2

# ╔═╡ 4a74c5b1-d19b-4add-a0ca-3a239d8a7c75
l*(I - Diagonal(l0)/Diagonal(l))/(Λ - l0)

# ╔═╡ 3914f4ec-4b70-4f05-9ec1-a9f70223bfd8
inv(Λ)[3, 3]*l[3, 3]*(1-l0[3, 3]/l[3, 3])

# ╔═╡ 94cf8ec3-5c81-4874-a331-aa0f15563808
a = (P*inv(Λ))'*(𝛴²ₕ₃)*(P*inv(Λ))

# ╔═╡ ffb3c3bb-827e-4173-90af-231bc45fef11
𝛴²₃

# ╔═╡ 2a1d988a-21fd-43b1-9701-dcfcf8f3d09b
A = (P*sqrt((I - Diagonal(l0)/Diagonal(l))/(Λ - l0)))'*(𝛴²₃)*P*sqrt((I - Diagonal(l0)/Diagonal(l))/(Λ - l0))

# ╔═╡ f51dc362-e8ab-494b-94f9-7ec63262afe4
B = (P*inv(Λ))'*𝛴²₀*P*inv(Λ)

# ╔═╡ 9ae7ea6f-6f70-4d92-9009-4ed3ccde3945
C = Real.(((I - sqrt(Diagonal(l0))/sqrt(Diagonal(l)))/(sqrt(Λ) - sqrt(l0))))'*ff

# ╔═╡ 5a0bd5e0-e3fb-47dc-8be2-fec23ea865a9
cov(C, dims=2)

# ╔═╡ 2df2e495-5170-4db3-904b-0c272f69b182
D = Real.(sqrt((I - Diagonal(l0)/Diagonal(l))/(Λ - l0)))'*ff0

# ╔═╡ 53142c10-19e9-42ba-b06d-85e326c8c575
cov(D, dims=2)

# ╔═╡ 008b10fc-8e80-4783-ae37-1f18147d6f55
scatter(C[2, :], .-C[3, :], aspect_ratio=:equal, ylabel="𝛼′", xlabel="𝛽′", title="Standardised PCA estimate", left_margin=5Plots.mm)

# ╔═╡ c43d7521-0e3d-403c-83b6-d9cce7675931
scatter!(.-𝛽, .-0.1.*𝛼, markersize=1, markercolor=:gray)

# ╔═╡ Cell order:
# ╠═09bee220-1249-11ec-3e2a-4fdde9976a59
# ╟─1cbf17f7-8322-4e3c-a895-163a48332205
# ╟─2f0cb4e5-4424-42d9-8eb7-e8bf128575ad
# ╟─0f10dc15-9153-4d8c-a5e5-c0318c459fa1
# ╟─230c66e5-7898-48af-9544-439677bee137
# ╠═fee5a67f-1678-45fc-83b3-3ffc187c19b3
# ╠═c3e707a4-1e35-4fe0-8e6d-88d0eca27fde
# ╠═5146ae95-d551-49ec-888a-71fc7bf65d14
# ╠═b646056a-7c34-4cb3-aac9-d255646f6b75
# ╟─2fbff3c4-9377-4624-90d6-d4ea7c8da64e
# ╟─404bb5c9-74ab-43a9-8791-67ac23b0c0f8
# ╟─4f01861c-dd52-462f-b340-ceecd3385c8e
# ╠═92670bd1-ec33-4590-9a53-ceca6f25ac0f
# ╠═8629c505-a1e4-4c73-b2d4-44ecc858373d
# ╟─c3c3d6c3-49a0-40a6-b41b-dbdfa20a7d6b
# ╠═9247ba48-f242-4479-a92d-e994b1610279
# ╠═c973721b-c20b-4573-a8fa-a9a9282fc87d
# ╟─4032a036-842b-462b-b75d-ebaaaf4aded9
# ╠═6a957b7d-8fff-46c5-b441-01f8b9d45eb7
# ╟─5686638b-ba3f-4691-9291-d91dfa510aa7
# ╠═76e702f7-521c-4445-aa6c-de28215d8d44
# ╠═b17a0c26-ce25-4e11-84dd-0d4b86371ae5
# ╠═69e0799e-18f4-4f19-9bae-44202c2bba5d
# ╠═cf454ccf-e505-4df9-a3cf-4afba2211d4a
# ╠═8013d3a2-a011-456b-9bfe-3ec781c31aba
# ╟─ce024bfd-b892-4a31-a0d1-1d4befc7b80e
# ╠═542aec21-e623-4e53-ae4b-a7f240bfcb08
# ╠═1daa139a-23c0-4995-bf73-abfbb7c518dd
# ╠═705068d7-033f-49d2-bb59-4b40aa514642
# ╠═f1a685e5-ab99-4a61-a2d9-41ebf33372b0
# ╠═3a7af0db-8e3d-4100-a7c2-dcd37253a9ce
# ╠═7a28c225-bab8-40ee-906e-0c5349accc84
# ╟─4d619965-6eb7-4528-804c-00bbf1572e60
# ╠═7ddd5530-7590-42a7-8e66-54432f717eec
# ╠═a95a3d4c-cb8b-4840-bc67-57fb2be00065
# ╠═688dc093-c173-4eb4-ad12-4ff2e899b9f9
# ╠═f4081452-9314-47a4-948d-2813bda53507
# ╠═5b853331-88c7-4b7f-b26e-8e0f963a5ed5
# ╠═1ee56732-671e-447e-adb8-cce45902b935
# ╠═70d2fe5f-7f7d-4891-8490-5c508c21c6f6
# ╟─bb52a2b3-e856-4785-abb9-78ba15382d11
# ╟─70d2f196-1604-4cf4-8e38-05dd8c099c83
# ╟─009c662a-9940-43ef-b0c9-d55560e6f49b
# ╠═46e7d2a0-8b62-4112-92df-74538f5f16d8
# ╠═595e5022-9099-45ab-b5b0-b92fedb96ba8
# ╟─0d062a65-adca-42c4-8b1c-84095d5335c5
# ╠═0f6b9065-55ad-4b1e-bcc3-2e177a0e331d
# ╠═fd506ba8-9a5b-4235-bc43-b730522dce4f
# ╟─58083319-d649-4c2d-b4b2-c2df7c49ee46
# ╠═2a9a4929-8755-4237-a454-fda125d548ca
# ╟─1762e05f-7eec-4411-865a-eba47531a49c
# ╠═feeaa835-28e0-4a24-905b-53347560ec9a
# ╠═2358f79f-f44d-4e93-b4da-75e7be715882
# ╠═d9fca85a-63b6-4ab5-ab0e-6cd5e68c4e06
# ╠═d700b4ae-8bd5-4437-b515-c21d3da70543
# ╠═06e226a6-0f0f-4ef6-b916-04d51a3789c8
# ╠═ea513127-e9bb-40a7-9ec0-2beee8fa6179
# ╟─f0a7f529-19ac-44f7-8ab1-8d72d596a0ef
# ╠═f946e3d7-15b9-45a9-951b-b3d77fc589f8
# ╠═652b998d-a6b0-4bb7-908e-c3dc7a82764c
# ╠═d192149c-1ef4-4761-9e78-879a97abee57
# ╠═14d5e08c-0808-4352-835e-c173c53fd8d7
# ╠═772d260a-72b3-4c9e-8e65-6368508cfb52
# ╠═0f83a03c-890d-4771-b83f-16d8d5db81e9
# ╠═b9a8e7e4-5f9b-462d-845e-b873f50bfc90
# ╠═f1df27ea-0537-4f4b-afe7-c0692ca7b26d
# ╟─2f987bde-a87a-43a8-8653-0c5812c56a70
# ╟─1d18e8c7-7eff-4b75-ad67-e61bc22e12ea
# ╠═3821aa68-e3b3-4764-a7ff-0b6b0e6bdc3e
# ╠═9d0c3f2b-8cf4-41f9-8fe8-5cd9a19dbd88
# ╠═b0d86eb8-6f94-424a-adb2-eae8e40be1fd
# ╟─e507beed-984e-4f8b-a59b-bb3dcff1416a
# ╠═ea23a142-1c36-4430-a598-2fcf0f2e3b66
# ╠═fe0902aa-b55e-4b72-a35e-77db673747c9
# ╟─9b92548e-ae54-48de-a5c9-4be88ae11228
# ╠═33167762-e2c6-4aba-98b7-58e28ee35ccb
# ╠═35dcd411-a342-4ee2-9146-f45ad39d6341
# ╟─d2404839-53bf-49d3-a9bf-a0ea81471c43
# ╠═e64a04de-9257-4bbb-84c4-007c532d9c72
# ╟─710d51ea-8196-4a7f-91fe-1a9bd3f06d01
# ╠═a4c8677f-6bdc-4bfe-adcc-3130f45b7aaa
# ╠═5ad01546-7daa-441f-bf2c-60e2a6e6f0c7
# ╠═36039419-6df3-45f2-ae2c-77cf6e3b405d
# ╠═d046eb57-d4cb-4eda-9a81-00a234af1d57
# ╠═4e9b0004-da7a-4cb3-afe9-c3365d37471a
# ╠═5e2281e9-360b-41dd-9f4f-38e14d8aea86
# ╠═fa43e761-0b22-4903-aa6c-47050d8c27c0
# ╟─bcd85e11-8bf0-4968-a528-8b56ce14db5b
# ╟─3fe79e11-947e-46e6-a5ef-1bf1b537c857
# ╠═15970d81-af01-443e-b54a-bda505b1b3be
# ╠═f771887a-77f0-4a4d-9e4c-7e988201197f
# ╠═b636df98-f935-46b4-bb08-be96cee18884
# ╠═c43d6053-4eeb-450a-bca4-bd5d58b179c1
# ╠═beafaf38-3753-4e9f-bf6d-69b98439fe18
# ╠═fd00b9e4-f2fc-4571-8596-b1e53f8d7286
# ╠═9d13988d-e625-42b5-8e89-86d519027912
# ╠═4f5c167d-dbde-4fed-bd2d-63b0100e2d1a
# ╠═82841686-b078-4ba2-b5ce-d748070437d2
# ╠═8dded5c7-f4ea-4ee5-b0b9-59b21c7f3097
# ╠═cb77a71b-7cba-489d-a51b-690cdd281a4c
# ╠═b00f965a-18a1-4ccf-85a9-ae79d588cab6
# ╠═0ae23619-0ce3-4e32-9a11-f9a444235c70
# ╠═ddb9528f-b57c-4333-9347-8a56aa6890c4
# ╠═c3067e6d-dfaf-4938-8bbf-4b51514fd286
# ╠═a2f4c245-031b-4f1c-9f44-b3947cbc9795
# ╠═ad7aac3b-303a-4ca6-a44f-f11cba8d115f
# ╠═dbdcd6c4-102d-41cc-b218-100f74355b85
# ╠═9d6334c1-7107-4fbb-b383-9a37fc56de74
# ╠═07e9a417-5adb-4182-a98d-fd618feaf675
# ╠═008512c1-3a44-42b4-8695-2a9f4d88ef01
# ╠═2826e49a-5566-4aa7-be9a-73799dfbebc7
# ╠═c79d0914-d0db-4f5c-b1f2-a6154d063e78
# ╠═eaae76e1-9326-432a-8fc5-b03c3d93e75a
# ╠═b7a9e831-6ab4-4afd-a5da-9d3b5414245b
# ╠═70109d54-cfcd-46c5-8c67-a3b96286ee09
# ╠═9970c6e2-a2c9-442e-8fca-17b8a5c360b9
# ╠═7408597f-8bc6-493a-80f6-79a5201d625b
# ╠═a031cdd3-7e92-45d8-8037-e6695cd588c2
# ╠═816e841b-3f3b-4b8c-8c4b-4cd93c685ab7
# ╠═923b9503-5a8e-4dab-984d-6ca07d625bb4
# ╠═9938d1b8-8f80-45b0-bdff-48cd74e78fa4
# ╠═4771733f-83a5-4c13-b0ea-cb54241f7e19
# ╠═d8847644-4d5b-4113-8870-e8c3dba7e1a2
# ╠═5463ac9d-9bf7-4363-abb9-d1de01a88681
# ╠═39871429-722a-4745-9a67-62c2961fb268
# ╠═31a016af-82fa-41a1-b2e8-8adea326c11c
# ╠═31cfab8a-e219-484c-9428-1fef8ea16ca5
# ╠═c64a1c15-824a-49d2-ac0f-10573a1472c5
# ╠═8c9e6bd3-c548-4126-a6fb-a100e91f8747
# ╠═24425756-1e29-4e78-9398-6ef1820940b8
# ╠═4a74c5b1-d19b-4add-a0ca-3a239d8a7c75
# ╠═3914f4ec-4b70-4f05-9ec1-a9f70223bfd8
# ╠═94cf8ec3-5c81-4874-a331-aa0f15563808
# ╠═ffb3c3bb-827e-4173-90af-231bc45fef11
# ╠═2a1d988a-21fd-43b1-9701-dcfcf8f3d09b
# ╠═f51dc362-e8ab-494b-94f9-7ec63262afe4
# ╠═9ae7ea6f-6f70-4d92-9009-4ed3ccde3945
# ╠═5a0bd5e0-e3fb-47dc-8be2-fec23ea865a9
# ╠═2df2e495-5170-4db3-904b-0c272f69b182
# ╠═53142c10-19e9-42ba-b06d-85e326c8c575
# ╠═008b10fc-8e80-4783-ae37-1f18147d6f55
# ╠═c43d7521-0e3d-403c-83b6-d9cce7675931


text:
# Accelerated inner product

C = zeros(m,p)

for j=1:m
    for l=1:p
        C[j,l] = dot(conj(view(A,j,:)),view(B,:,l))
    end
end

max_stars_repo_path: Code/Numerical Analysis & Machine Learning/Numerical Linear Algebra/InnerProduct.jl

max_stars_repo_name: BambooFlower/Math-Scripts

max_issues_repo_name: BambooFlower/Math-Scripts

max_forks_repo_name: BambooFlower/Math-Scripts


text:
struct CompoundSample{T, EBS, S} <: EBayesSample{T}
    vec::AbstractVector{EBS}
    probs::AbstractVector{S}
    Z::T
end

response(Z::CompoundSample) = Z.Z
nuisance_parameter(Z::CompoundSample) = (Z.vec, Z.probs)

function Base.isless(a::CompoundSample, b::CompoundSample)
    nuisance_parameter(a) == nuisance_parameter(b) ||
        throw(ArgumentError("Comparison only implemented for homoskedastic compound samples"))
    Base.isless(response(a), response(b))
end

function compound(Zs::AbstractVector{<:EBayesSample})
    n = nobs(Zs)
    Zs_no_data = set_response.(Zs)
    Zs_unique_dict = countmap(Zs_no_data)
    Zs_unique = collect(keys(Zs_unique_dict))
    Zs_unique_prob = values(Zs_unique_dict)./n
    CompoundSample.(Ref(Zs_unique), Ref(Zs_unique_prob), response.(Zs))
end

function compound(Zs::MultinomialSummary{<:EBayesSample})
    n = nobs(Zs)
    Zs_no_data = set_response.(collect(keys(Zs)))
    _mult = fweights(multiplicity(Zs))
    Zs_unique_dict = countmap(Zs_no_data, _mult)
    Zs_unique = collect(keys(Zs_unique_dict))
    Zs_unique_prob = collect(values(Zs_unique_dict))./n
    Zs_compound = CompoundSample.(Ref(Zs_unique), Ref(Zs_unique_prob), response.(collect(keys(Zs))))
    summarize(Zs_compound, _mult)
end

function likelihood_distribution(Z::CompoundSample, μ)
    n = length(Z.vec)
    MixtureModel(likelihood_distribution.(Z.vec, μ), Z.probs)
end

function marginalize(Z::CompoundSample, prior::Distribution)
    n = length(Z.vec)
    MixtureModel(marginalize.(Z.vec, Ref(prior)),  Z.probs)
end

# TODO: # beware a bit of this, hopefully compound won't be used too much
function skedasticity(Zs::AbstractArray{<:CompoundSample,1})
    Empirikos.Homoskedastic()
end

# hard code as missing.
struct HeteroskedasticSamples{EBS <: EBayesSample{Missing}, S}
    vec::AbstractVector{EBS}
    probs::AbstractVector{S}
end

function heteroskedastic(Zs::AbstractVector{<:EBayesSample})
    _comp = compound(Zs)[1]
    HeteroskedasticSamples( _comp.vec, _comp.probs)
end

function heteroskedastic(Zs::MultinomialSummary{<:EBayesSample})
    _comp = collect(keys(compound(Zs)))[1]
    HeteroskedasticSamples( _comp.vec, _comp.probs)
end

function CompoundSample(Zs::HeteroskedasticSamples)
    CompoundSample(Zs.vec, Zs.probs, missing)
end


text:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-------------------------
"""
    mes = meshmeasures(p, t)

compute areas of every element
"""
function meshareas(p, t)
    d12 = p[t[:, 2], :] - p[t[:, 1], :]
    d13 = p[t[:, 3], :] - p[t[:, 1], :]
    mes = (d12[:, 1] .* d13[:, 2] .- d12[:, 2] .* d13[:, 1]) / 2.
    return reshape(abs.(mes), length(mes))
end
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-------------------------
"""
   I, J = indKM_sparse(t)

non zero indices in stiffness and mass P1 sparse matrices
"""
function indKM_sparse(t)
    it1, it2, it3 = t[:, 1], t[:, 2], t[:, 3]
    I = vcat(it1, it2, it3, it2, it3, it1, it1, it2, it3)
    J = vcat(it2, it3, it1, it1, it2, it3, it1, it2, it3)
    return I, J
end
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-------------------------
"""
     SK, SM = assembKM_P12D(p, t)

assemble the P1 stifness and mass matrices coefficients.
"""
function assembKM_P12D(p, t)
    np, dim, nt = size(p, 1), size(p, 2), size(t, 1)
    # compute ar add define arK and arM (not used ?)
    ar = meshareas(p, t)
    arK, arM = ar, ar
    it1, it2, it3 = t[:, 1], t[:, 2], t[:, 3]
    g1x, g1y, g2x, g2y, g3x, g3y = pdetrg(p, t, ar)
    c3 = (g1x .* g2x .+ g1y .* g2y) .* arK
    c1 = (g2x .* g3x .+ g2y .* g3y) .* arK
    c2 = (g3x .* g1x .+ g3y .* g1y) .* arK
    SK = vcat(c3, c1, c2, c3, c1, c2, -c2 - c3, -c3 - c1, -c1 - c2)
    # mass matrix
    aod = arM / 12. # off diagonal element
    ad = 2 * aod    # diagonal element
    SM = vcat(aod, aod, aod, aod, aod, aod, ad, ad, ad)
    return SK, SM
end
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-------------------------
"""
     nothing = assembK_P12D_inplace!(p, t, S)

assemble the P1 stifness matrice coefficients in place. Useful for SparseDiffTools.
"""
function assembK_P12D_inplace!(pin, t, S)
    p = if typeof(pin) <: Vector
        reshape(pin, length(pin) ÷ 2, 2)
    else
        pin
    end
    np, dim, nt = size(p, 1), size(p, 2), size(t, 1)
    # compute ar add define arK and arM (not used ?)
    ar = meshareas(p, t)
    arK, arM = ar, ar
    it1, it2, it3 = t[:, 1], t[:, 2], t[:, 3]
    g1x, g1y, g2x, g2y, g3x, g3y = pdetrg(p, t, ar)
    c3 = (g1x .* g2x .+ g1y .* g2y) .* arK
    c1 = (g2x .* g3x .+ g2y .* g3y) .* arK
    c2 = (g3x .* g1x .+ g3y .* g1y) .* arK
    SK = vcat(c3, c1, c2, c3, c1, c2, -c2 - c3, -c3 - c1, -c1 - c2)
    copy!(S, SK)
    return nothing
end
"""
     nothing = assembM_P12D_inplace!(p, t, S)

assemble the P1 mass matrice coefficients in place. Useful for SparseDiffTools.
"""
function assembM_P12D_inplace!(pin, t, S)
    p = if typeof(pin) <: Vector
        reshape(pin, length(pin) ÷ 2, 2)
    else
        pin
    end
    np, dim, nt = size(p, 1), size(p, 2), size(t, 1)
    # compute ar add define arK and arM (not used ?)
    arM = meshareas(p, t)
    aod = arM / 12. # off diagonal element
    ad = 2 * aod    # diagonal element
    SM = vcat(aod, aod, aod, aod, aod, aod, ad, ad, ad)
    copy!(S, SM)
    return nothing
end
#-------------------------------------------------------------------------------
"""
   g1x, g1y, g2x, g2y, g3x, g3y = pdetrg(p, t, ar)

returns the  the gradient components of the triangle base functions
"""
function pdetrg(p, t, ar)
    # corner point indices
    a1, a2, a3 =  t[:, 1], t[:, 2], t[:, 3]
    # triangle sides
    r23x =  p[a3, 1] .- p[a2, 1]
    r23y =  p[a3, 2] .- p[a2, 2]
    r31x =  p[a1, 1] .- p[a3, 1]
    r31y =  p[a1, 2] .- p[a3, 2]
    r12x =  p[a2, 1] .- p[a1, 1]
    r12y =  p[a2, 2] .- p[a1, 2]
    g1x =  - 0.5 .* r23y ./ ar
    g1y =    0.5 .* r23x ./ ar
    g2x =  - 0.5 .* r31y ./ ar
    g2y =    0.5 .* r31x ./ ar
    g3x =  - 0.5 .* r12y ./ ar
    g3y =    0.5 .* r12x ./ ar
    return g1x, g1y, g2x, g2y, g3x, g3y
end
#-------------------------------------------------------------------------------
"""
    be = btri(T::Array{Int64, 2})

finds boundary edges
"""
function btri(T::Array{Int64, 2})
    # form all faces, non - duplicates are surface triangles
    faces = vcat(T[:, [1, 2]], T[:, [1, 3]], T[:, [2, 3]])
    faces = sort(faces, dims = 2)
    _, ix, jx = uniquerows_int(faces)
    ee = collect(1:(maximum(jx) + 1)) .- 1e-8
    result = fit(Histogram, jx, (1:(maximum(jx) + 1)) .- 1e-8,
                 closed = :right)
    vec = result.weights
    qx = findall(x -> x != 0, vec .- 2)
    return faces[ix[qx], :]
end
#-------------------------------------------------------------------------------
function uniquerows_int(x::Matrix{Int64})
    Vx = [x[k, :] for k = 1:size(x, 1)]
    jj = groupslices(Vx)
    ix, jx = firstinds(jj), similar(jj)
    for k = 1:length(ix)
        jx[jj .== ix[k]] .= k
    end
    Vy = Vx[ix]
    y = [Vy[k][l] for k = 1:length(Vy), l = 1:length(Vy[1])]
    return y, ix, jx
end


text:
"""
    isearlytermination(thread::Thread)

Check if we have met a termination condition for the thread.
"""
function isearlytermination(thread::Thread)
    # Ensure, minimum depth 1 searched.
    # Ensure, minumum 1023+ nodes searched.
    # And that time management is up.
    (thread.ss.depth > 1) && ((thread.ss.nodes & 1023) == 1023) && (elapsedtime(thread.timeman) > thread.timeman.max_time) &&
    (!thread.timeman.isinfinite) && (thread.timeman.depth == 0)
end

"""
    find_best_move()

Probe the tablebase if appropriate, or perform the absearch routine.
"""
function find_best_move(thread::Thread, ttable::TT_Table, depth::Int)::Int
    board = thread.board

    # Probe the syzygy tablebase.
    # At present, it is assumed 5 men tablebases are available.
    # Also assumes we want to always probe.
    if (count(occupied(board)) <= 5)
        res = tb_probe_root(board)
        if res !== TB_RESULT_FAILED
            # Extract the move to the PV and return an evaluation.
            return interpret_syzygy(thread, res)
        end
    end

    ABORT_SIGNAL[] = false
    initTimeManagement!(thread.timeman)

    # If no tablebase result found, we proceed with iterative deepening and absearch.
    eval = iterative_deepening(thread, ttable, depth)
    ABORT_SIGNAL[] = true
    return eval
end

"""
    iterative_deepening(thread::Thread, ttable::TT_Table, max_depth::Int)

Performs the iterative deepening approach to absearch, and forms the main core of the engine routine.
"""
function iterative_deepening(thread::Thread, ttable::TT_Table, max_depth::Int)::Int
    eval = -MATE
    # Below is an idea to reduce history heuristic "noise" between searches.
    # Testing shows some benefit.
    for i in 1:2
        for j in 1:64
            @inbounds @simd for k in 1:64
                thread.history[i][j][k] = fld(thread.history[i][j][k], 2)
            end
        end
    end
    for i in 1:6
        for j in 1:64
            for k in 1:6
                @inbounds @simd for l in 1:64
                    thread.counterhistory[i][j][k][l] = fld(thread.counterhistory[i][j][k][l], 2)
                    thread.followhistory[i][j][k][l] = fld(thread.followhistory[i][j][k][l], 2)
                end
            end
        end
    end
    # We iterate through depths until we either:
    # a) reach the maximum depth, or
    # b) hit a time management constraint.
    best_move = MOVE_NONE
    for depth in 1:max_depth
        eval = aspiration_window(thread, ttable, depth, eval)
        if thread.stop
            thread.pv[1].list[1] = best_move
            break
        else
            best_move = thread.pv[1][1]
        end
        if (!thread.timeman.isinfinite && elapsedtime(thread.timeman) > thread.timeman.max_time) ||
            (thread.timeman.depth !== 0 && depth >= thread.timeman.depth) || istermination(thread.timeman)
            break
        end
    end
    return eval
end

"""
    aspiration_window(thread::Thread, ttable::TT_Table, depth::Int, eval::Int)

When operating within the iterative deepening framework, aspiration windows are used via calls to this function.
"""
function aspiration_window(thread::Thread, ttable::TT_Table, depth::Int, eval::Int)::Int
    δ = 22 + fld(abs(eval), 256)
    if depth >= WINDOW_DEPTH
        α = max(-MATE, eval - δ)
        β = min(MATE, eval + δ)
    end
    return aspiration_window_internal(thread, ttable, depth, eval, -MATE, MATE, δ)
end

# The internals of aspiration_window
function aspiration_window_internal(thread::Thread, ttable::TT_Table, depth::Int, eval::Int, α::Int, β::Int, δ::Int)::Int
    eval = eval
    thread.ss.seldepth = 0
    thread.ss.depth = depth

    while !thread.stop
        eval = absearch(thread, ttable, α, β, depth, 0, false)

        # window cond met
        if ((α < eval < β) || (elapsedtime(thread.timeman) > 2.5)) && !thread.stop
            uci_report(thread, ttable, α, β, eval)
            return eval
        end

        # fail low
        if eval <= α
            β = fld(α + β, 2)
            α = max(-MATE, α - δ)

        # fail high
        elseif eval >= β
            β = min(MATE, β + δ)
        end

        # expand window
        δ += fld(δ, 2)
    end
    return eval
end

# https://www.chessprogramming.org/Quiescence_Search
"""
    qsearch()

Quiescence search function.
"""
function qsearch(thread::Thread, ttable::TT_Table, α::Int, β::Int, ply::Int)::Int
    board = thread.board
    pv = thread.pv

    # ensure pv is clear
    clear!(thread.pv[ply + 1])

    # default val
    tt_eval = -MATE
    tt_move = MOVE_NONE

    thread.ss.seldepth = max(thread.ss.seldepth, ply)
    thread.ss.nodes += 1

    # Check for abort signal, or TimeManagement criteria.
    if isearlytermination(thread) || (ABORT_SIGNAL[] == true)
        thread.stop = true
        return 0
    end

    # Check if the current position is drawn. (Ignoring stalemate)
    if isdrawbymaterial(board) || is50moverule(board)
        return draw_noise(thread)
    end

    # If we have reached the maximum search ply depth, we stop here.
    if ply >= MAX_PLY
        return evaluate(board, thread.pktable)
    end

    # probe the transposition table
    tt_entry = getTTentry(ttable, board.hash)
    if tt_entry !== NO_ENTRY
        tt_eval = Int(tt_entry.eval)
        tt_move = tt_entry.move
        tt_value = Int(ttvalue(tt_entry, ply))
        if (tt_entry.bound == BOUND_EXACT) ||
            ((tt_entry.bound == BOUND_LOWER) && (tt_value >= β)) ||
            ((tt_entry.bound == BOUND_UPPER) && (tt_value <= α))
            return tt_value
        end
    end

    # We may treat checks in qsearch specially.
    if ischeck(board)
        best = -MATE
        margin = 0
    else
        if tt_eval !== -MATE
            eval = thread.evalstack[ply + 1] = tt_eval
        elseif (ply > 0) && (thread.movestack[ply] == NULL_MOVE)
            eval = thread.evalstack[ply + 1] = -thread.evalstack[ply] + 2*TEMPO_BONUS
        else
            eval = thread.evalstack[ply + 1] = evaluate(board, thread.pktable)
        end

        best = eval

        # Evaluation pruning step.
        if eval > α
            α = eval
        end
        if α >= β
            return eval
        end

        # Delta pruning step. / Futility pruning.
        # If a "best case" (but perhaps non-existent) move, plus a small margin, is not enough to raise alpha, we can stop here.
        # Rather not perform this when in check.
        margin = α - eval - Q_FUTILITY_MARGIN
        if !ischeck(board) && (optimistic_move_estimator(board) < margin)
            return eval
        end
    end

    moveorder = thread.moveorders[ply + 1]

    if ischeck(board)
        # We need evasions, so generate quiet moves.
        init_normal_moveorder!(thread, tt_move, ply)
        @inbounds thread.killer1s[ply + 2] = MOVE_NONE
        @inbounds thread.killer2s[ply + 2] = MOVE_NONE
    else
        # We just look at noisy moves.
        # We make sure that the moves satisfy a SEE margin.
        init_noisy_moveorder!(thread, ply, max(1, margin))
    end

    best = qsearch_internal(thread, ttable, α, β, ply, tt_move, best)
    clear!(moveorder)
    return best
end

# The internals of qsearch.
function qsearch_internal(thread::Thread, ttable::TT_Table, α::Int, β::Int, ply::Int, tt_move::Move, best::Int)::Int
    # We can select skipquiets = false during qsearch, as if it's a NOISY type moveorder, quiets are skipped anyway.
    # Otherwise, quiets are needed to generate king check evasions.
    played = 0
    skipquiets = false
    @inbounds pv_current = thread.pv[ply + 1]
    @inbounds pv_future = thread.pv[ply + 2]
    while ((move = selectmove!(thread, ply, skipquiets)) !== MOVE_NONE) && !thread.stop

        if ischeck(thread.board) && (played >= 4) && (best > -MATE + MAX_PLY) && !istactical(thread.board, move)
            skipquiets = true
            continue
        end

        u = apply_move!(thread, move)
        played += 1
        eval = -qsearch(thread, ttable, -β, -α, ply + 1)
        undo_move!(thread, move, u)

        # Check for improvements, and update PV.
        # Ensure we met no stop signals.
        if (eval > best) && (!thread.stop)
            best = eval
            if eval > α
                α = best
                clear!(pv_current)
                push!(pv_current, move)
                updatepv!(pv_current, pv_future)
            end
        end

        # fail high?
        if (α >= β) && !thread.stop
            return best
        end
    end
    return best
end

"""
    absearch()

The main alpha-beta search function, containing all the pruning heuristics.
"""
function absearch(thread::Thread, ttable::TT_Table, α::Int, β::Int, depth::Int, ply::Int, cutnode::Bool)::Int
    board = thread.board
    @inbounds pv_current = thread.pv[ply + 1]
    @inbounds pv_future = thread.pv[ply + 2]

    # Set initial history stats to zero.
    hist = cmhist = fmhist = 0

    # Initial α value for later computation.
    init_α = α
    clear!(pv_current)

    # Is this the root node?
    isroot = iszero(ply)

    # Is this a pvnode?
    pvnode = β !== α + 1

    # Best evaluation defaults to -MATE.
    best = -MATE

    # default tt_eval, tt_move
    tt_eval = -MATE
    tt_move = MOVE_NONE

    # ensure +ve depth
    if depth < 0
        depth = 0
    end

    # If depth reaches 0, we enter quiescence search.
    # If we are in check, we prefer to hold off one more ply.
    if iszero(depth) && !ischeck(board)
        q_eval = qsearch(thread, ttable, α, β, ply)
        return q_eval
    end

    # update thread details
    thread.ss.seldepth = max(thread.ss.seldepth, ply)
    thread.ss.nodes += 1

    # Check for abort signal, or TimeManagement criteria.
    if isearlytermination(thread) || (ABORT_SIGNAL[] == true)
        thread.stop = true
        return evaluate(board, thread.pktable)
    end

    # If this is not the root node, we can check for early exit conditions.
    if isroot == false

        # If we detect a draw, we can return an evaluation of zero.
        if isdrawbymaterial(board) || is50moverule(board) || isrepetition(board)
            return draw_noise(thread)
        end

        # If we reach the max ply, we evaluate the board here and return.
        if ply >= MAX_PLY
            eval = evaluate(board, thread.pktable)
            return eval
        end

        # Check for mating evaluations, and prune for the best mate.
        # This avoids the need to look for sub-optimal continuations.
        if α > -MATE + ply
            mate_α = α
        else
            mate_α = -MATE + ply
        end
        if β < MATE - ply - 1
            mate_β = β
        else
            mate_β = MATE - ply - 1
        end
        if mate_α >= mate_β
            return mate_α
        end
    end

    # Probe the transposition table.
    tt_entry = getTTentry(ttable, board.hash)
    if tt_entry !== NO_ENTRY
        tt_eval = Int(tt_entry.eval)
        tt_value = Int(ttvalue(tt_entry, ply))
        tt_move = tt_entry.move
        if (tt_entry.depth >= depth) && (depth == 0 || (pvnode == false))
            if (tt_entry.bound == BOUND_EXACT) ||
                ((tt_entry.bound == BOUND_LOWER) && (tt_value >= β)) ||
                ((tt_entry.bound == BOUND_UPPER) && (tt_value <= α))
                return tt_value
            end
        end
    end

    # If we have less than 5 pieces, we are in the realms of the Syzygy tablebase.
    # Attempt to probe.
    if (count(occupied(board)) <= 5) && !isroot
        _eval = tb_probe_wdl(board)
        if _eval !== TB_RESULT_FAILED
            thread.ss.tbhits += 1

            # Is the tablebase score losing?
            if iszero(_eval)
                eval = -MATE + MAX_PLY + ply + 1
                tt_bound = BOUND_UPPER

            # Is the tablebase a draw? Blessed / cursed loss and wins are draws
            elseif 1 <= _eval <= 3
                eval = 0
                tt_bound = BOUND_EXACT

            # Else, the tablebase is a win.
            else
                eval = MATE - MAX_PLY - ply - 1
                tt_bound = BOUND_LOWER
            end

            # add to transposition table
            if (tt_bound == BOUND_EXACT) ||
                ((tt_bound == BOUND_LOWER) && (eval >= β)) ||
                ((tt_bound == BOUND_UPPER) && (eval <= α))
                setTTentry!(ttable, board.hash, eval, MOVE_NONE, MAX_PLY - 1, tt_bound)
                return eval
            end
        end
    end

    # Set a static evaluation.
    # If the last move was a NULL move, the evaluation is simply a change in 2*TEMPO_BONUS (flip a sign).
    if tt_eval !== -MATE
        eval = thread.evalstack[ply + 1] = tt_eval
        if eval == 0
            eval = 2 * (thread.ss.nodes & 1) - 1
        end
    elseif (ply > 0) && (thread.movestack[ply] == NULL_MOVE)
        eval = thread.evalstack[ply + 1] = -thread.evalstack[ply] + 2*TEMPO_BONUS
    else
        eval = thread.evalstack[ply + 1] = evaluate(board, thread.pktable)
    end

    # Reset killer moves for the upcoming ply.
    @inbounds thread.killer1s[ply + 2] = MOVE_NONE
    @inbounds thread.killer2s[ply + 2] = MOVE_NONE

    # Razoring.
    # If the evaluation plus a small margin is still below alpha, we drop into the quiescence search.
    if (pvnode === false) && (ischeck(board) === false) && (depth <= RAZOR_DEPTH) && (eval + RAZOR_MARGIN <= α)
        q_eval = qsearch(thread, ttable, α, β, ply)
        return q_eval
    end

    # Check to see if the evaluation history has improved over the last few plies.
    # We can use this to modify some pruning heuristics later on.
    improving = ((ply >= 2) && (thread.evalstack[ply + 1] > thread.evalstack[ply - 1])) ? 2 : 1

    # Beta pruning.
    # If the evaluation minus a margin is still better than beta, we can prune here.
    if (pvnode === false) && (ischeck(board) === false) && (depth <= BETA_PRUNE_DEPTH) && (eval - BETA_PRUNE_MARGIN * (depth - improving + 1) >= β)
        return eval
    end

    # Null move pruning.
    # Check we are greater than depth 2, eval > β.
    # Not done in a pvnode.
    # Not done when in check.
    # Check we have non-pawn material left on the board.
    # Do not do more than 2 b2b null moves.
    if (depth >= 2) && (eval >= β) && !pvnode && !ischeck(board) &&
        (ply > 0 ? (thread.movestack[ply] !== NULL_MOVE) : true)  &&
        (ply > 1 ? (thread.movestack[ply - 1] !== NULL_MOVE) : true) &&
        ((tt_entry === NO_ENTRY) || (tt_eval >= β)) &&
        has_non_pawn_material(board)

        reduction = fld(depth, 5) + 3 + min(fld(eval - β, 150), 3)

        u = apply_null!(thread)

        cand_eval = -absearch(thread, ttable, -β, -β + 1, depth - reduction, ply + 1, !cutnode)
        undo_null!(thread, u)
        if (cand_eval >= β)
            return β
        end
    end

    best_move = MOVE_NONE

    # Set pruning variables used during the main loop.
    futility_margin = FUTILITY_MARGIN * depth
    see_quiet_margin = SEE_QUIET_MARGIN * depth
    see_noisy_margin = SEE_NOISY_MARGIN * depth^2
    skipquiets = false
    played = 0
    num_quiets = 0

    # Probcut pruning.
    # We are not in a PV-node.
    # A potential exists on the board such that it can beat β plus a margin.
    if !pvnode &&  (depth >= PROBCUT_DEPTH) && (abs(β) < (MATE - MAX_PLY)) &&
        (eval >= β || eval + optimistic_move_estimator(board) >= β + PROBCUT_MARGIN)
        probcut_count = 0
        raised_β = min(β + PROBCUT_MARGIN, MATE - MAX_PLY - 1)
        init_noisy_moveorder!(thread, ply, raised_β - eval)
        while ((move = selectmove!(thread, ply, true)) !== MOVE_NONE) && !thread.stop && (probcut_count < 2 + (cutnode ? 2 : 0))
            probcut_count += 1
            u = apply_move!(thread, move)
            # We check that the move holds in a qsearch.
            eval = -qsearch(thread, ttable, -β, -β + 1, ply + 1)
            # If it holds, and is above raised β, we do a fuller search.
            if eval >= raised_β
                eval = -absearch(thread, ttable, -β, -β + 1, depth - 4, ply + 1, !cutnode)
            end
            undo_move!(thread, move, u)

            if (eval >= raised_β)
                clear!(thread.moveorders[ply + 1])
                return eval
            end
        end
        clear!(thread.moveorders[ply + 1])
    end

    # Internal iterative deepening step.
    # Aim is to perform a quick search to fill the transposition table with potential moves.
    if (depth > 6) && (tt_move === MOVE_NONE)
        absearch(thread, ttable, α, β, depth - 7, ply + 1, cutnode)
        tt_entry = getTTentry(ttable, board.hash)
        if tt_entry !== NO_ENTRY
            tt_move = tt_entry.move
        end
    end

    # Init the move ordering
    init_normal_moveorder!(thread, tt_move, ply)
    @inbounds quiets_tried = thread.quietstack[ply + 1]
    @inbounds moveorder = thread.moveorders[ply + 1]
    while ((move = selectmove!(thread, ply, skipquiets)) !== MOVE_NONE) && !thread.stop

        isquiet = !istactical(board, move)

        # If we pick a quiet move, extract the history statistics.
        if isquiet
            hist, cmhist, fmhist = gethistory(thread, move, ply)
            num_quiets += 1
        end

        # Quiet move pruning steps.
        # We require that a non-mating line exists before pruning.
        # Opt not to prune if there is only king + pawn material remaining.
        if isquiet && (best > -MATE + MAX_PLY) && has_non_pawn_material(board)
            # Quiet move futility pruning step, using history.
            if (depth <= FUTILITY_PRUNE_DEPTH) && (eval + futility_margin <= α) &&
                (hist + cmhist + fmhist < FUTILITY_LIMIT[improving])
                skipquiets = true
            end

            # Quiet move futility pruning, using an additional margin.
            if (depth <= FUTILITY_PRUNE_DEPTH) && (eval + futility_margin + FUTILITY_MARGIN_NOHIST <= α)
                skipquiets = true
            end

            # Late move pruning.
            if (depth <= LATE_MOVE_PRUNE_DEPTH) && (num_quiets >= LATE_MOVE_COUNT[improving][depth + 1])
                skipquiets = true
            end

            # Counter move pruning
            if (moveorder.stage > STAGE_COUNTER) && (depth <= COUNTER_PRUNE_DEPTH[improving]) && (cmhist < COUNTER_PRUNE_LIMIT[improving])
                continue
            end

            # Follow up move pruning
            if (moveorder.stage > STAGE_COUNTER) && (depth <= FOLLOW_PRUNE_DEPTH[improving]) && (fmhist < FOLLOW_PRUNE_LIMIT[improving])
                continue
            end
        end

        # Prune moves which fail the static exchange evaluator.
        # Only ran if our best evaluation is not a mating line.
        if (moveorder.stage > STAGE_GOOD_NOISY) && (depth <= SEE_PRUNE_DEPTH) && (best > -MATE + MAX_PLY) &&
            (static_exchange_evaluator(board, move, isquiet ? see_quiet_margin : see_noisy_margin) == false)
            continue
        end

        u = apply_move!(thread, move)
        played += 1
        if isquiet && (quiets_tried.idx <= MAX_QUIET_TRACK)
            push!(quiets_tried, move)
        end

        # Late move reduction calculations.
        # This allows later searches to probe to a lower depth under given scenarios.
        if (depth > 2) && (played > 1) && (isquiet || cutnode || skipquiets)
            reduction = @inbounds LMRTABLE[min(depth, 64)][min(played, 64)]
            # Below are quiet-move specific reduction conditions.
            if isquiet
                if !pvnode
                    reduction += 1
                end
                # If we are not improving, increase the reduction depth.
                if isone(improving)
                    reduction += 1
                end
                # Add reductions for cutnodes
                if cutnode
                    reduction += 2
                end
                # Killer moves, and counter moves, are worth looking at more.
                if moveorder.stage < STAGE_INIT_QUIET
                    reduction -= 1
                end
                # Alter if we are moving our king when in check.
                if ischeck(board) && (type(board[from(move)]) === KING)
                    reduction += 1
                end
                # Adjust on the history
                reduction -= max(-2, min(2, fld(hist + cmhist + fmhist, 5000)))
            elseif (depth <= 7) && (played >= 3)
                reduction += 1
            end
            reduction = min(depth - 1, max(reduction, 1))
        else
            reduction = 1
        end

        # do we need an extension?
        hist_thresh = 14000 + 10*ply^2
        if (ischeck(board) || (isquiet && num_quiets <= 4 && cmhist >= hist_thresh && fmhist >= hist_thresh)) && !isroot
            newdepth = depth + 1
        else
            newdepth = depth
        end

        # perform search, taking into account LMR
        if reduction !== 1
            cand_eval = -absearch(thread, ttable, -α - 1, -α, newdepth - reduction, ply + 1, true)
        end
        if ((reduction !== 1) && (cand_eval > α)) || (reduction == 1 && !(pvnode && played == 1))
            cand_eval = -absearch(thread, ttable, -α - 1, -α, newdepth - 1, ply + 1, !cutnode)
        end
        if (pvnode && (played == 1 || ((cand_eval > α) && (isroot || cand_eval < β))))
            cand_eval = -absearch(thread, ttable, -β, -α, newdepth - 1, ply + 1, false)
        end

        # Revert move.
        undo_move!(thread, move, u)

        # Have we found a better move?
        # Ensure we encountered no stops, and update the PV.
        if (cand_eval > best) && (!thread.stop)
            best = cand_eval
            best_move = move
            if cand_eval > α
                α = cand_eval
                # Update the PV
                clear!(pv_current)
                push!(pv_current, best_move)
                updatepv!(pv_current, pv_future)

                # Did we fail high? Then we can break.
                if α >= β
                    break
                end
            end
        end
    end

    if iszero(played)
        clear!(moveorder)
        clear!(quiets_tried)
        if ischeck(board)
            # add ply to give an indication of the "fastest" mate
            return -MATE + ply
        else
            # Stalemate is a draw.
            return 0
        end
    end

    # Update the history heuristics.
    # This is done when a quiet move fails high.
    if (best >= β) && (best_move !== MOVE_NONE) && !istactical(board, best_move)
        updatehistory!(thread, quiets_tried, ply, depth^2)
    end

    clear!(moveorder)
    clear!(quiets_tried)

    # Update the transposition table.
    # We don't do this at the root node.
    if isroot == false
        tt_bound = best >= β ? BOUND_LOWER : (best > init_α ? BOUND_EXACT : BOUND_UPPER)
        setTTentry!(ttable, board.hash, eval, best_move, depth, tt_bound)
    end
    return best
end

"""
    static_exchange_evaluator(board::Board, move::Move)

Returns true if a move passes a static exchange criteria, false otherwise.
"""
const SEE_VALUES = @SVector [200, 775, 775, 1250, 2400, 15000]

function static_exchange_evaluator(board::Board, move::Move, threshold::Int)
    from_sqr = Int(from(move))
    to_sqr = Int(to(move))

    move_flag = flag(move)

    from_piece = piece(board, from_sqr)
    to_piece = piece(board, to_sqr)
    victim = (move_flag < __KNIGHT_PROMO) ? from_piece : makepiece(PieceType((move_flag >> 12) - 3), board.turn)

    occ = occupied(board)
    occ ⊻= (Bitboard(from_sqr) | Bitboard(to_sqr))

    if move_flag === __ENPASS
        occ ⊻= Bitboard(board.enpass)
    end

    attackers = (pawns(board) & pawnAttacks(!board.turn, to_sqr) & friendly(board)) |
    (pawns(board) & pawnAttacks(board.turn, to_sqr) & enemy(board)) |
    (knightMoves(to_sqr) & knights(board)) |
    (kingMoves(to_sqr) & kings(board))
    if !isempty(bishoplike(board))
        attackers |= (bishopMoves(to_sqr, occ) & bishoplike(board))
    end
    if !isempty(rooklike(board))
        attackers |= (rookMoves(to_sqr, occ) & rooklike(board))
    end

    attackers &= occ

    if isempty(attackers & enemy(board))
        return true
    end

    color = !board.turn

    balance = -threshold

    if to_piece !== BLANK
        @inbounds balance += SEE_VALUES[type(to_piece).val]
    end

    if move_flag >= __KNIGHT_PROMO
        @inbounds balance += SEE_VALUES[type(victim).val] - SEE_VALUES[1]
    end

    if move_flag === __ENPASS
        @inbounds balance += SEE_VALUES[1]
    end

    if balance < 0
        return false
    end

    @inbounds balance -= SEE_VALUES[type(victim).val]

    if balance >= 0
        return true
    end

    while true
        our_attackers = attackers & board[color]

        # if we can't attack, we lose
        if isempty(our_attackers)
            break
        end

        # find weakest piece to recapture
        for i in 1:6
            piecetype = PieceType(i)
            if isempty(our_attackers & board[piecetype]) == false
                victim = piecetype
                break
            end
        end

        # remove our attacker
        occ ⊻= Bitboard(poplsb(our_attackers & board[victim])[1])

        # check for diag moves
        if (victim === PAWN || victim === BISHOP || victim === QUEEN)
            attackers |= bishopMoves(to_sqr, occ) & bishoplike(board)
        end

        # check for rank/file moves
        if (victim === ROOK || victim === QUEEN)
            attackers |= rookMoves(to_sqr, occ) & rooklike(board)
        end

        attackers &= occ

        balance = -balance - 1 - SEE_VALUES[victim.val]
        color = !color

        if balance >= 0
            break
        end

    end
    # if it's your turn, you lost the SEE loop
    if board.turn === color
        return false
    else
        return true
    end
end

# delta pruning
function optimistic_move_estimator(board::Board)
    # assume pawn at minimum
    value = SEE_VALUES[1]

    # find highest val targets
    for i in 5:-1:2
        piecetype = PieceType(i)
        if isempty(board[board.turn] & board[piecetype]) == false
            @inbounds value = SEE_VALUES[i]
            break
        end
    end

    # promo checks
    if isempty(pawns(board) & board[board.turn] & (board.turn == WHITE ? RANK_7 : RANK_2)) == false
        @inbounds value += SEE_VALUES[5] - SEE_VALUES[1]
    end

    return value
end

function updatepv!(pv_current::MoveStack, pv_future::MoveStack)
    for tmp_pv_move in pv_future
        push!(pv_current, tmp_pv_move)
    end
end

function draw_noise(thread::Thread)
    return 2 * (thread.ss.nodes & 1) - 1
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
import Distributions.pdf, Distributions.logpdf, Random.rand, Base.length

struct MixedSupport <: ValueSupport end

"""
    Factored{N} <: Distribution{Multivariate, MixedSupport}
a `Distribution` type that can be used to combine multiple `UnivariateDistribution`'s and sample from them.
Example: it can be used as `prior = Factored(Normal(0,1), Uniform(-1,1))`
"""
struct Factored{N} <: Distribution{Multivariate,MixedSupport}
    p::NTuple{N,UnivariateDistribution}
    Factored(args::UnivariateDistribution...) = new{length(args)}(args)
end
"""
    pdf(d::Factored, x) = begin
Function to evaluate the pdf of a `Factored` distribution object
"""
function pdf(d::Factored{N}, x) where {N}
    s = pdf(d.p[1], x[1])
    for i = 2:N
        s *= pdf(d.p[i], x[i])
    end
    s
end

"""
    logpdf(d::Factored, x) = begin
Function to evaluate the logpdf of a `Factored` distribution object
"""
function logpdf(d::Factored{N}, x) where {N}
    s = logpdf(d.p[1], x[1])
    for i = 2:N
        s += logpdf(d.p[i], x[i])
    end
    s
end

"""
    rand(rng::AbstractRNG, factoreddist::Factored)
function to sample one element from a `Factored` object
"""
rand(rng::AbstractRNG, factoreddist::Factored{N}) where {N} =
    ntuple(i -> rand(rng, factoreddist.p[i]), Val(N))

"""
    length(p::Factored) = begin
returns the number of distributions contained in `p`.
"""
length(p::Factored{N}) where {N} = N

export Factored


text:
module  SummaryInfo

using PopGen
using DataFrames
using Test

cats = @nancycats;

@testset "alleleaverage" begin
    @test length(alleleaverage(cats,  rounding = true)) == 2
    @test alleleaverage(cats,  rounding = true) isa NamedTuple
    @test length(alleleaverage(cats,  rounding = false)) == 2
    @test alleleaverage(cats,  rounding = false) isa NamedTuple
end

@testset "richness" begin
    @test richness(cats, by = "locus") isa DataFrame
    @test size(richness(cats, by = "locus")) == (9, 2)
    @test richness(cats, by = "population") isa DataFrame
    @test size(richness(cats, by = "population")) == (153, 3)
end

@testset "summary F/D/etc. stats" begin
    smry_glob = summary(cats)
    smry_loc = summary(cats, by = "locus")
    @test smry_glob isa DataFrame
    @test smry_loc isa DataFrame
    @test size(smry_glob) == (1, 10)
    @test size(smry_loc) ==(9, 11)
end

end # module


text:
function resample(x::AbstractArray{T,2}, fact::Integer=10, fun::Function=first) where {T <: Real}
    nrow, ncol = size(x)
    nrow2 = fld(nrow, fact)
    ncol2 = fld(ncol, fact)

    # type = ifelse(fun in [first, last], T, Float32)
    out = zeros(T, nrow2, ncol2)

    @inbounds for i in 1:nrow2
        indx = (i-1)*fact + 1 : i*fact
        @inbounds for j in 1:ncol2
            indy = (j-1)*fact + 1 : j*fact
            # out[i, j] = fun(x[indx, indy])
            out[i, j] = fun(view(x, indx, indy));
        end
    end
    out
end

function resample(x::AbstractArray{T,3}, fact::Integer=10, fun::Function=first) where {T <: Real}
    nrow, ncol, ntime = size(x)
    nrow2 = fld(nrow, fact)
    ncol2 = fld(ncol, fact)

    out = zeros(T, nrow2, ncol2, ntime)
    @inbounds for i in 1:ntime
        # out[:, :, i] = resample(view(x, :, :, i), fact, fun)
        out[:, :, i] = resample(@view(x[:, :, i]), fact, fun)
    end
    out
end

function resample(x::GeoArray, fact::Integer=10, fun::Function=first)
    # reconstruct the a.f
    GeoArray(resample(x.A, fact, fun), x.f, x.crs)
end


text:
using .DistributionsAD: TuringDirichlet, TuringWishart, TuringInverseWishart,
                        FillVectorOfUnivariate, FillMatrixOfUnivariate, VectorOfUnivariate,
                        MatrixOfUnivariate, FillVectorOfMultivariate, VectorOfMultivariate

# TuringDirichlet

_clamp(x, ::TuringDirichlet) = _clamp(x, SimplexBijector())
bijector(d::TuringDirichlet) = SimplexBijector()

function link(
    d::TuringDirichlet,
    x::AbstractVecOrMat{<:Real},
    proj::Bool = true,
)
    return SimplexBijector{proj}()(x)
end

function invlink(
    d::TuringDirichlet,
    y::AbstractVecOrMat{<:Real},
    proj::Bool = true,
)
    return inv(SimplexBijector{proj}())(y)
end

function logpdf_with_trans(
    d::TuringDirichlet,
    x::AbstractVecOrMat{<:Real},
    transform::Bool,
)
    ϵ = _eps(eltype(x))
    lp = logpdf(d, x .+ ϵ)
    if transform
        lp -= logabsdetjac(bijector(d), x)
    end
    return lp
end

# TuringWishart

function logpdf_with_trans(
    d::TuringWishart,
    X::AbstractMatrix{<:Real},
    transform::Bool
)
    return _logpdf_with_trans_pd(d, X, transform)
end
function logpdf_with_trans(
    d::TuringWishart,
    X::AbstractArray{<:AbstractMatrix{<:Real}},
    transform::Bool
)
    return mapvcat(X) do x
        _logpdf_with_trans_pd(d, x, transform)
    end
end
link(d::TuringWishart, X::AbstractMatrix{<:Real}) = PDBijector()(X)
invlink(d::TuringWishart, Y::AbstractMatrix{<:Real}) = inv(PDBijector())(Y)
function getlogp(d::TuringWishart, Xcf, X)
    return 0.5 * ((d.df - (dim(d) + 1)) * logdet(Xcf) - tr(d.chol \ X)) - d.c0
end

# TuringInverseWishart

function logpdf_with_trans(
    d::TuringInverseWishart,
    X::AbstractMatrix{<:Real},
    transform::Bool
)
    _logpdf_with_trans_pd(d, X, transform)
end
function logpdf_with_trans(
    d::TuringInverseWishart,
    X::AbstractArray{<:AbstractMatrix{<:Real}},
    transform::Bool
)
    return mapvcat(X) do x
        _logpdf_with_trans_pd(d, x, transform)
    end
end
link(d::TuringInverseWishart, X::AbstractMatrix{<:Real}) = PDBijector()(X)
invlink(d::TuringInverseWishart, Y::AbstractMatrix{<:Real}) = inv(PDBijector())(Y)
function getlogp(d::TuringInverseWishart, Xcf, X)
    Ψ = d.S
    return -0.5 * ((d.df + dim(d) + 1) * logdet(Xcf) + tr(Xcf \ Ψ)) - d.c0
end

# filldist and arraydist

function logpdf_with_trans(
    dist::FillVectorOfUnivariate,
    x::AbstractVector{<:Real},
    istrans::Bool,
)
    return _sum(x) do x
        logpdf_with_trans(dist.v.value, x, istrans)
    end
end
function logpdf_with_trans(
    dist::FillVectorOfUnivariate,
    x::AbstractMatrix{<:Real},
    istrans::Bool,
)
    return vec(sum(logpdf_with_trans(dist.v.value, x, istrans), dims = 1))
end
function link(
    dist::FillVectorOfUnivariate,
    x::AbstractVector{<:Real},
)
    return link(dist.v.value, x)
end
function link(
    dist::FillVectorOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return link(dist.v.value, x)
end
function invlink(
    dist::FillVectorOfUnivariate,
    x::AbstractVector{<:Real},
)
    return invlink(dist.v.value, x)
end
function invlink(
    dist::FillVectorOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return invlink(dist.v.value, x)
end

function logpdf_with_trans(
    dist::FillMatrixOfUnivariate,
    x::AbstractMatrix{<:Real},
    istrans::Bool,
)
    return _sum(x) do x
        logpdf_with_trans(dist.dists.value, x, istrans)
    end
end
function logpdf_with_trans(
    dist::FillMatrixOfUnivariate,
    x::AbstractArray{<:AbstractMatrix{<:Real}},
    istrans::Bool,
)
    return mapvcat(x) do x
        logpdf_with_trans(dist, x, istrans)
    end
end
function link(
    dist::FillMatrixOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return link(dist.dists.value, x)
end
function invlink(
    dist::FillMatrixOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return invlink(dist.dists.value, x)
end

function logpdf_with_trans(
    dist::FillVectorOfMultivariate,
    x::AbstractMatrix{<:Real},
    istrans::Bool,
)
    return sum(logpdf_with_trans(dist.dists.value, x, istrans))
end
function logpdf_with_trans(
    dist::FillVectorOfMultivariate,
    x::AbstractArray{<:AbstractMatrix{<:Real}},
    istrans::Bool,
)
    return mapvcat(x) do x
        logpdf_with_trans(dist, x, istrans)
    end
end
function link(
    dist::FillVectorOfMultivariate,
    x::AbstractMatrix{<:Real},
)
    return link(dist.dists.value, x)
end
function invlink(
    dist::FillVectorOfMultivariate,
    x::AbstractMatrix{<:Real},
)
    return invlink(dist.dists.value, x)
end

function logpdf_with_trans(
    dist::VectorOfMultivariate,
    x::AbstractMatrix{<:Real},
    istrans::Bool,
)
    return _sumeachcol(x, dist.dists) do c, dist
        logpdf_with_trans(dist, c, istrans)
    end
end
function logpdf_with_trans(
    dist::VectorOfMultivariate,
    x::AbstractArray{<:AbstractMatrix{<:Real}},
    istrans::Bool,
)
    return mapvcat(x) do x
        logpdf_with_trans(dist, x, istrans)
    end
end
function link(
    dist::VectorOfMultivariate,
    x::AbstractMatrix{<:Real},
)
    return eachcolmaphcat(x, dist.dists) do c, dist
        link(dist, c)
    end
end
function invlink(
    dist::VectorOfMultivariate,
    x::AbstractMatrix{<:Real},
)
    return eachcolmaphcat(x, dist.dists) do c, dist
        invlink(dist, c)
    end
end

function logpdf_with_trans(
    dist::VectorOfUnivariate,
    x::AbstractVector{<:Real},
    istrans::Bool,
)
    return _sum(dist.v, x) do dist, x
        logpdf_with_trans(dist, x, istrans)
    end
end
function logpdf_with_trans(
    dist::VectorOfUnivariate,
    x::AbstractMatrix{<:Real},
    istrans::Bool,
)
    return mapvcat(1:size(x,2)) do i
        c = view(x, :, i)
        sum(logpdf_with_trans.(dist.v, c, istrans))
    end
end
function link(
    dist::VectorOfUnivariate,
    x::AbstractVector{<:Real},
)
    return mapvcat(dist.v, x) do dist, x
        link(dist, x)
    end
end
function link(
    dist::VectorOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return mapvcat(eachcol(x)) do x
        link(dist, x)
    end
end

function logpdf_with_trans(
    dist::MatrixOfUnivariate,
    x::AbstractMatrix{<:Real},
    istrans::Bool,
)
    return _sum(dist.dists, x) do dist, x
        logpdf_with_trans(dist, x, istrans)
    end
end
function logpdf_with_trans(
    dist::MatrixOfUnivariate,
    x::AbstractArray{<:AbstractMatrix{<:Real}},
    istrans::Bool,
)
    return mapvcat(x) do x
        logpdf_with_trans(dist, x, istrans)
    end
end
function link(
    dist::MatrixOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return mapvcat(dist.dists, x) do dist, x
        link(dist, x)
    end
end
function invlink(
    dist::MatrixOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return mapvcat(dist.dists, x) do dist, x
        invlink(dist, x)
    end
end


text:
using VsRecorder
using VsRecorder.Scenes: TeamView

@testset "Team View" begin
    config = vs_setup(PokemonBattle, language = VsI18n.ZHS, init_descriptors = false)
    ctx = VsRecorderBase.initialize(config)

    expected = """Dragapult (F) @ Life Orb
Ability: Clear Body
Level: 60
- Protect
- Draco Meteor
- Thunderbolt
- Shadow Ball

Tyranitar (M) @ Weakness Policy
Ability: Sand Stream
Level: 55
- Rock Slide
- Crunch
- High Horsepower
- Protect

Grimmsnarl (M) @ Iapapa Berry
Ability: Prankster
Level: 100
- Play Rough
- Sucker Punch
- Fake Out
- Power Swap

Arcanine (M) @ Mago Berry
Ability: Intimidate
Level: 50
- Snarl
- Flamethrower
- Safeguard
- Protect

Gastrodon (F) @ Sitrus Berry
Ability: Storm Drain
Level: 50
- Muddy Water
- Earth Power
- Yawn
- Recover

Corviknight (F) @ Leftovers
Ability: Mirror Armor
Level: 50
- Roost
- Bulk Up
- Brave Bird
- Body Press
"""
    frame = VsFrame(image = VsRecorder.load_data("scenes", "zhs", "team-view.png"))
    parsed = vs_tryparse_scene(TeamView, frame, ctx)
    team = parsed.team
    # Tesseract has slightly different behaviours on different platforms
    @test team.title ∈ ("达拉斯R赛", "達拉斯R賽")
    @test team.author == "露露"
    @test team.rental_code == "0000 0008 53NO V8"
    @test export_team(team) == expected
end


text:
#!/usr/bin/env julia

# Copyright (c) 2016, Jost von Hardenberg - ISAC-CNR, Italy

using RainFARM
using ArgParse
using Compat, Compat.Printf

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--nf", "-n"
            help = "Subdivisions for downscaling [default: 2]"
            arg_type = Int
            default = 2
        "--outfile", "-o", "--out"
            help = "Weights file"
            arg_type = AbstractString
            default = "grid.nc"
        "infile"
            help = "The input file to downscale"
            arg_type = AbstractString
            required = true
        "varname"
            help = "Input variable name"
            arg_type = AbstractString
            required = true
    end
	s.description="Generates a sample file with a grid equal to the input file downscaled by a factor NF"
    return parse_args(s)
end

args = parse_commandline()
nf=args["nf"]
filenc=args["infile"]
outfile=args["outfile"]
varnc=args["varname"]

println("Reading file ",filenc)

(pr,lon_mat,lat_mat)=read_netcdf2d(filenc, varnc);

# Creo la griglia fine
nss=size(pr)
if (length(nss)>=3)
    pr=pr[:,:,1]
end
ns=nss[1];

(lon_f, lat_f)=lon_lat_fine(lon_mat, lat_mat,nf);

println("Output size: ",size(lon_f))

write_netcdf2d(outfile,reshape(pr,ns,ns,1),lon_f,lat_f,varnc,filenc)


text:
@testset "MPolyGcd" begin
  QA, A = PolynomialRing(Hecke.QQ, "A")
  K, a = number_field(A^3 - 2, "a")
  Kx, (x,) = PolynomialRing(K, ["x"])
  g = gcd(-715//2*a^2*x^3 - 715//2*a^2*x^2 - 365525875//2*a^2*x - 365525875//2*a^2, (-715//2*a^2 + 511225)*x^2 - 365525875//2*a^2 + 261351000625)
  @test g == x^2 + 511225
end


text:
@generated function det{T}(A::StaticMatrix{T})
    if size(A) == (1,1)
        return quote
            $(Expr(:meta, :inline))
            @inbounds return A[1]
        end
    elseif size(A) == (2,2)
        return quote
            $(Expr(:meta, :inline))
            @inbounds return A[1]*A[4] - A[3]*A[2]
        end
    elseif size(A) == (3,3)
        return quote
            $(Expr(:meta, :inline))
            #@inbounds a = A[5]*A[9] - A[8]*A[6]
            #@inbounds b = A[8]*A[3] - A[2]*A[9]
            #@inbounds c = A[2]*A[6] - A[5]*A[3]
            #@inbounds return A[1]*a + A[4]*b + A[7]*c

            @inbounds x0 = SVector(A[1], A[2], A[3])
            @inbounds x1 = SVector(A[4], A[5], A[6])
            @inbounds x2 = SVector(A[7], A[8], A[9])
            return vecdot(x0, cross(x1, x2))
        end
    else
        S = typeof((one(T)*zero(T) + zero(T))/one(T))
        return quote # Implementation from Base
            if istriu(A) || istril(A)
                return convert($S, det(UpperTriangular(A)))::$S # Is this a Julia bug that a convert is not type stable??
            end
            AA = convert(Array{$S}, A)
            return det(lufact(AA))
        end
    end
end


text:
ts = TS(df_timetype_index)

@test index(ts) == df_timetype_index[!, :Index]


text:
# This file was generated, do not modify it. # hide
ensemble = machine(ensemble_model, X, y)
estimates = evaluate!(ensemble, resampling=CV())
estimates |> pprint