text:
withenv("LINES"=>24, "COLUMNS"=>80) do
    @testset "sizing" begin
        seed!(1337)
        p = @inferred spy(sprand(10,10,.15))
        @test p isa Plot
        @test_reference(
            "references/spy/default_10x10.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(sprand(10,15,.15))
        @test_reference(
            "references/spy/default_10x15.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(sprand(15,10,.15))
        @test_reference(
            "references/spy/default_15x10.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(sprand(2000,200,.0001))
        @test_reference(
            "references/spy/default_2000x200.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(sprand(200,2000,.0001))
        @test_reference(
            "references/spy/default_200x2000.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(sprand(2000,2000,.1))
        @test_reference(
            "references/spy/default_overdrawn.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(sprandn(200,200,.001))
        @test_reference(
            "references/spy/default_200x200_normal.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        @test_reference(
            "references/spy/default_200x200_normal_nocolor.txt",
            @io2str(show(IOContext(::IO, :color=>false), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(Matrix(sprandn(200,200,.001)))
        @test_reference(
            "references/spy/default_200x200_normal.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(sprandn(200,200,.001), width=10)
        @test_reference(
            "references/spy/default_200x200_normal_small.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(sprandn(200,200,.001), height=5)
        @test_reference(
            "references/spy/default_200x200_normal_small.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = @inferred spy(sprandn(200,200,.001), height=5, width=20)
        @test_reference(
            "references/spy/default_200x200_normal_misshaped.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
    end

    @testset "parameters" begin
        seed!(1337)
        p = @inferred spy(sprandn(200,200,.001), color=:green)
        @test_reference(
            "references/spy/parameters_200x200_green.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
        @test_reference(
            "references/spy/parameters_200x200_green_nocolor.txt",
            @io2str(show(IOContext(::IO, :color=>false), p)),
            render = BeforeAfterFull()
        )
        seed!(1337)
        p = spy(sprandn(200,200,.001), title="Custom Title", canvas=DotCanvas, border=:ascii)
        @test_reference(
            "references/spy/parameters_200x200_dotcanvas.txt",
            @io2str(show(IOContext(::IO, :color=>true), p)),
            render = BeforeAfterFull()
        )
    end
end


text:
module Pairs

export defaultbyfirst, setbyfirst, getbyfirst, setkv, getkv, rmkv

import ..compat_findfirst
const findfirst = compat_findfirst

"""
    setbyfirst(collection, item) -> item

Set `item` in a `collection`.
If `first() of an exisiting item matches `first(item)` it is replaced.
Otherwise the new `item` is inserted at the end of the `collection`.
"""
function setbyfirst(c, item, eq = ==)
    k = first(item)
    i = findfirst(x->eq(first(x), k), c)
    if i === nothing
        push!(c, item)
    else
        c[i] = item
    end
    return item
end

"""
    getbyfirst(collection, key [, default]) -> item

Get `item` from collection where `first(item)` matches `key`.
"""
function getbyfirst(c, k, default=nothing, eq = ==)
    i = findfirst(x->eq(first(x), k), c)
    return i === nothing ? default : c[i]
end

"""
    defaultbyfirst(collection, item)

If `first(item)` does not match match `first()` of any existing items,
insert the new `item` at the end of the `collection`.
"""
function defaultbyfirst(c, item, eq = ==)
    k = first(item)
    if findfirst(x->eq(first(x), k), c) === nothing
        push!(c, item)
    end
    return
end

"""
    setkv(collection, key, value)

Set `value` for `key` in collection of key/value `Pairs`.
"""
setkv(c, k, v) = setbyfirst(c, k => v)

"""
    getkv(collection, key [, default]) -> value

Get `value` for `key` in collection of key/value `Pairs`,
where `first(item) == key` and `value = item[2]`
"""
function getkv(c, k, default=nothing)
    i = findfirst(x->first(x) == k, c)
    return i === nothing ? default : c[i][2]
end

"""
    rmkv(collection, key)

Remove `key` from `collection` of key/value `Pairs`.
"""
function rmkv(c, k, default=nothing)
    i = findfirst(x->first(x) == k, c)
    if !(i === nothing)
        deleteat!(c, i)
    end
    return
end

end # module Pairs

ext: jl


text:
"""
    twosite_variance(AL, AC, VL, VR, O)

```
  ┌────────┐    ┌────────┐
┌─1   AL   2────1   AC   2─┐
│ └────3───┘    └───3────┘ │
│      │            │      │
│    .─3────────────4──.   │
│   (         O         )  │
│    `─1────────────2──'   │
│      │            │      │
│ ┌────3───┐    ┌───3────┐ │
└─1conj(VL)2┐  ┌1conj(VR)2─┘
  └────────┘│  │└────────┘
(1)─────────┘  └─────────(2)
```
"""
function twosite_variance(
    AL::AbstractTensor3,
    AC::AbstractTensor3,
    VL::AbstractTensor3,
    VR::AbstractTensor3,
    O::AbstractTensor4,
)
    @tensoropt ((blp, brp, klp, krp) = 1, (x, y, z, l, r) = χ) twosite_variance[l, r] :=
        AL[x, y, klp] *
        AC[y, z, krp] *
        conj(VL)[x, l, blp] *
        conj(VR)[r, z, brp] *
        O[blp, brp, klp, krp]
end
function twosite_variance(
    AL::AbstractTensor3,
    AC::AbstractTensor3,
    VL::AbstractTensor3,
    O::AbstractTensor4,
)
    @tensoropt ((blp, brp, klp, krp) = 1, (x, y, z, l, r) = χ) twosite_variance[l, r] :=
        AL[x, y, klp] *
        AC[y, z, krp] *
        conj(VL)[x, l, blp] *
        conj(VL)[z, r, brp] * # VR = permutedims(VL, (2, 1, 3))
        O[blp, brp, klp, krp]
end

function vl_from_al(AL::AbstractTensor3)
    P, D = size(AL, 3), size(AL, 2)
    VL = permutedims(reshape(nullspace(_tensor_to_matrix(AL, true)'), D, P, :), (1, 3, 2))
    size(VL, 2) == D * (P - 1) || @warn "failed to solve nullspace in VL"
    VL
end
function vr_from_ar(AR::AbstractTensor3)
    P, D = size(AR, 3), size(AR, 2)
    VR = reshape(nullspace(_tensor_to_matrix(AR, false))', :, D, P)
    size(VR, 1) == D * (P - 1) || @warn "failed to solve nullspace in VR"
    AR
end

function enlargestep(
    TV::AbstractMatrix{U},
    AL::Tensor3{T},
    AR::Tensor3{T},
    AC::Tensor3{T},
    C::Matrix{T},
    VL::Tensor3{T},
    VR::Tensor3{T};
    maxenlarge = size(C, 1) * (size(AC, 3) - 1),
    ϵ = 0.0,
) where {T,U}

    W = promote_type(T, U)
    d, p = size(C, 1), size(AC, 3)

    x, y, z = svd(TV)
    enlarge_d = min(sum(y .> ϵ), maxenlarge)
    enlarge_L = mul_matrix_from_right(VL, x[:, 1:enlarge_d])
    enlarge_R = mul_matrix_from_left(VR, z[:, 1:enlarge_d]')

    AL_new, AR_new, AC_new, C_new = zeros(W, d + enlarge_d, d + enlarge_d, p),
    zeros(W, d + enlarge_d, d + enlarge_d, p),
    zeros(W, d + enlarge_d, d + enlarge_d, p),
    zeros(W, d + enlarge_d, d + enlarge_d)

    AL_new[1:d, 1:d, :] .= AL
    AR_new[1:d, 1:d, :] .= AR
    AC_new[1:d, 1:d, :] .= AC
    C_new[1:d, 1:d] .= C
    AL_new[1:d, d+1:end, :] .= enlarge_L
    AR_new[d+1:end, 1:d, :] .= enlarge_R

    (al = AL_new, ar = AR_new, ac = AC_new, c = C_new), (norm(y), norm(y[enlarge_d+1:end]))
end
function enlargestep(
    TV::AbstractMatrix{U},
    AL::Tensor3{T},
    AC::Tensor3{T},
    C::Matrix{T},
    VL::Tensor3{T};
    maxenlarge = size(C, 1) * (size(AC, 3) - 1),
    ϵ = 0.0,
) where {T,U}

    W = promote_type(T, U)
    d, p = size(C, 1), size(AC, 3)
    VR = permutedims(VL, (2, 1, 3))

    x, y, z = svd(TV)
    enlarge_d = min(sum(y .> ϵ), maxenlarge)
    enlarge_L = mul_matrix_from_right(VL, x[:, 1:enlarge_d])

    AL_new, AC_new, C_new = zeros(W, d + enlarge_d, d + enlarge_d, p),
    zeros(W, d + enlarge_d, d + enlarge_d, p),
    zeros(W, d + enlarge_d, d + enlarge_d)

    AL_new[1:d, 1:d, :] .= AL
    AC_new[1:d, 1:d, :] .= AC
    C_new[1:d, 1:d] .= C
    AL_new[1:d, d+1:end, :] .= enlarge_L

    (al = AL_new, ac = AC_new, c = C_new), (norm(y), norm(y[enlarge_d+1:end]))
end


text:
using DistStat, Random, LinearAlgebra, LoopVectorization

mutable struct COXUpdate
    maxiter::Int
    step::Int
    verbose::Bool
    tol::Real
    function COXUpdate(; maxiter::Int=100, step::Int=10, verbose::Bool=false, tol::Real=1e-10)
        maxiter > 0 || throw(ArgumentError("maxiter must be greater than 0."))
        tol > 0 || throw(ArgumentError("tol must be positive."))
        new(maxiter, step, verbose, tol)
    end
end

"""
    breslow_ind(x)

returns indexes of result of cumulutive sum corresponding to "W". `x` is assumed to be nonincreasing.
"""
function breslow_ind(x::AbstractVector)
    uniq = unique(x)
    lastinds = findlast.(isequal.(uniq), [x])
    invinds = findfirst.(isequal.(x), [uniq])
    lastinds[invinds]
end

mutable struct COXVariables{T, A}
    m::Int # rows, number of subjects
    n::Int # cols, number of predictors
    β::MPIVector{T,A} # [n]
    β_prev::MPIVector{T,A} # r x [n]
    δ::A # indicator for right censoring (0 if censored)
    λ::T # regularization parameter
    t::AbstractVector{T} # a vector containing timestamps, should be in descending order.
    breslow::A
    σ::T # step size. 1/(2 * opnorm(X)^2) for guaranteed convergence
    grad::MPIVector{T,A}
    w::A
    W::A
    W_dist::MPIMatrix{T,A}
    q::A # (1-π)δ, distributed
    eval_obj::Bool
    function COXVariables(X::MPIMatrix{T,A}, δ::A, λ::AbstractFloat,
                        t::A;
                        σ::T=convert(T, 1/(2*opnorm(X; verbose=true)^2)),
                        eval_obj=false) where {T,A}
        m, n = size(X)
        β = MPIVector{T,A}(undef, n)
        β_prev = MPIVector{T,A}(undef, n)
        fill!(β, zero(T))
        fill!(β_prev, zero(T))

        δ = convert(A{T}, δ)
        breslow = convert(A{Int}, breslow_ind(convert(Array, t)))

        grad = MPIVector{T,A}(undef, n)
        w = A{T}(undef, m)
        W = A{T}(undef, m)
        W_dist = MPIMatrix{T,A}(undef, 1, m)
        q = A{T}(undef, m)

        new{T,A}(m, n, β, β_prev, δ, λ, t, breslow, σ, grad, w, W, W_dist, q, eval_obj)
    end
end

function reset!(v::COXVariables{T,A}; seed=nothing) where {T,A}
    fill!(v.β, zero(T))
    fill!(v.β_prev, zero(T))
end

function soft_threshold(x::T, λ::T) ::T where T <: AbstractFloat
    @assert λ >= 0 "Argument λ must be greater than or equal to zero."
    x > λ && return (x - λ)
    x < -λ && return (x + λ)
    return zero(T)
end

function π_δ!(out, w, W_dist, δ, breslow, W_range)
    # fill `out` with zeros beforehand.
    m = length(δ)
    W_base = minimum(W_range) - 1
    W_local = W_dist.localarray
    @avx for i in 1:m
        outi = zero(eltype(w))
        for j in 1:length(W_range)
            outi += ifelse(breslow[i] <= breslow[j + W_base], δ[j + W_base] * w[i] / W_local[j], zero(eltype(w)))
        end
        out[i] = outi
    end
    DistStat.Barrier()
    DistStat.Allreduce!(out)
    return out
end

function get_breslow!(out, cumsum_w, bind)
    out .= cumsum_w[bind]
    out
end

function cox_grad!(out, w, W, W_dist, t, q, X, β, δ, bind)
    T = eltype(β)
    m, n = size(X)
    mul!(w, X, β)
    w .= exp.(w)
    cumsum!(q, w) # q is used as a dummy variable
    get_breslow!(W, q, bind)
    W_dist .= distribute(reshape(W, 1, :))
    fill!(q, zero(eltype(q)))
    π_δ!(q, w, W_dist, δ, bind, W_dist.partitioning[DistStat.Rank()+1][2])
    q .= δ .- q
    mul!(out, transpose(X), q) # ([n]) = (n x [m]) x (m)
    out
end

cox_grad!(v::COXVariables{T,A}, X) where {T,A} = cox_grad!(v.grad, v.w, v.W, v.W_dist, v.t, v.q, X, v.β, v.δ, v.breslow)

function update!(X::MPIArray, u::COXUpdate, v::COXVariables{T,A}) where {T,A}
    #mul!(v.tmp_m_local1, X, v.β) # {m} = {m x [n]} * {[n]}.
    #v.tmp_m_local1 .= exp.(v.tmp_m_local1) # w
    #cumsum!(v.tmp_m_local2, v.tmp_m_local1) # W. TODO: deal with ties.
    #v.tmp_1m .= distribute(reshape(v.tmp_m_local2, 1, :)) # W_dist: distribute W.
    #v.tmp_mm .= v.π_ind .* v.tmp_m_local1 ./ v.tmp_1m # (π_ind .* w) ./ W_dist. computation order is determined for CuArray safety.
    #pd = mul!(v.tmp_m_local1, v.tmp_mm, v.δ) # {m} = {m x [m]} * {m}.
    #v.tmp_m_local2 .= v.δ .- pd # {m}.
    #grad = mul!(v.tmp_n, transpose(X), v.tmp_m_local2) # {[n]} = {[n] x m} * {m}.
    cox_grad!(v, X)
    v.β .= soft_threshold.(v.β .+ v.σ .* v.grad , v.λ) # {[n]}.
end

function get_objective!(X::MPIArray, u::COXUpdate, v::COXVariables{T,A}) where {T,A}
    v.grad .= (v.β .!= 0) # grad is dummy
    nnz = sum(v.grad)

    if v.eval_obj
        v.w .= exp.(mul!(v.w, X, v.β))
        cumsum!(v.q, v.w) # q is dummy
        get_breslow!(v.W, v.q, v.breslow)
        obj = dot(v.δ, mul!(v.q, X, v.β) .- log.(v.W)) .- v.λ .* sum(abs.(v.β))
        return false, (obj, nnz)
    else
        v.grad .= abs.(v.β_prev .- v.β)
        return false, (maximum(v.grad), nnz)
    end
end

function cox_one_iter!(X::MPIArray, u::COXUpdate, v::COXVariables)
    copyto!(v.β_prev, v.β)
    update!(X, u, v)
end

function loop!(X::MPIArray, u, iterfun, evalfun, args...)
    converged = false
    t = 0
    while !converged && t < u.maxiter
        t += 1
        iterfun(X, u, args...)
        if t % u.step == 0
            converged, monitor = evalfun(X, u, args...)
            if DistStat.Rank() == 0
                println(t, ' ', monitor)
            end
        end
    end
end

function cox!(X::MPIArray, u::COXUpdate, v::COXVariables)
    loop!(X, u, cox_one_iter!, get_objective!, v)
end

include("cmdline.jl")
opts = parse_commandline_cox()
if DistStat.Rank() == 0
    println("world size: ", DistStat.Size())
    println(opts)
end

m = opts["rows"]
n = opts["cols"]
iter = opts["iter"]
interval = opts["step"]
T = Float64
A = Array
using CUDA
if opts["gpu"]
    A = CuArray

    # All the GPU-related functions here
    function π_δ_kernel!(out, w, W_dist, δ, breslow, W_range)
        # fill `out` with zeros beforehand.
        idx_x = (blockIdx().x-1) * blockDim().x + threadIdx().x
        stride_x = blockDim().x * gridDim().x
        W_base = minimum(W_range) - 1
        for i in idx_x:stride_x:length(out)
            for j in W_range
                @inbounds if breslow[i] <= breslow[j]
                    out[i] += δ[j] * w[i] / W_dist[j - W_base]
                end
            end
        end
    end

    function π_δ!(out::CuArray, w::CuArray, W_dist, δ, breslow, W_range)
        numblocks = ceil(Int, length(w)/256)
        CUDA.@sync begin
            @cuda threads=256 blocks=numblocks π_δ_kernel!(out, w, W_dist.localarray, δ, breslow, W_range)
        end
        DistStat.Allreduce!(out)
        out
    end

    function breslow_kernel!(out, cumsum_w, bind)
        idx_x = (blockIdx().x-1) * blockDim().x + threadIdx().x
        stride_x = blockDim().x * gridDim().x
        for i = idx_x: stride_x:length(out)
            out[i]=cumsum_w[bind[i]]
        end
    end

    function get_breslow!(out::CuArray, cumsum_w::CuArray, bind)
        numblocks = ceil(Int, length(out)/256)
        CUDA.@sync begin
            @cuda threads=256 blocks=numblocks breslow_kernel!(out, cumsum_w, bind)
        end
        out
    end
end
if opts["Float32"]
    T = Float32
end
init_opt = opts["init_from_master"]
seed = opts["seed"]
eval_obj = opts["eval_obj"]
censor_rate = opts["censor_rate"]
lambda = opts["lambda"]

X = MPIMatrix{T, A}(undef, m, n)
rand!(X; common_init=init_opt, seed=seed)
δ = convert(A{T}, rand(m) .> censor_rate)
DistStat.Bcast!(δ) # synchronize the choice for δ.

uquick = COXUpdate(;maxiter=2, step=1, verbose=true)
u = COXUpdate(;maxiter=iter, step=interval, verbose=true)
# for simulation run, we just assume that the data are in reversed order of survivial time
t = convert(A{T}, collect(reverse(1:size(X,1))))
v = COXVariables(X, δ, lambda, t; eval_obj=eval_obj)
cox!(X, uquick, v)
reset!(v; seed=seed)
if DistStat.Rank() == 0
    @time cox!(X, u, v)
else
    cox!(X, u, v)
end


text:
"""
`trivial_selection!`: Trivial selection, copying over particles

### Arguments
* `E` - particle ensemble
"""
function trivial_selection!(E::TE) where {TE<:Ensemble}
    @. E.ω̂ = E.ω
    @. E.b̂ = E.b
    @. E.ξ̂ = deepcopy(E.ξ)
    @. E.d̂ = deepcopy(E.d)

    E
end

"""
`repopulate!`: After allocating the number of offspring of each particle, copy
the particles over.

### Arguments
* `E` - particle ensemble
* `B` - bin data structure
### Optional Arguments
"""
function repopulate!(E::TE, B::TB) where {TE<:Ensemble,TB<:Bins}
    n_particles = length(E)
    n_spawned = 0
    # copy over the particles
    for i = 1:n_particles
        # identify the bin of the current particle
        @inbounds bin = E.b[i]
        for k = 1:E.o[i]
            @inbounds E.ξ̂[k+n_spawned] = deepcopy(E.ξ[i])
            @inbounds E.ω̂[k+n_spawned] = B.ν[bin] / B.target[bin]
            @inbounds E.b̂[k+n_spawned] = bin
            @inbounds E.d̂[k+n_spawned] = deepcopy(E.d[i])
        end
        @inbounds n_spawned += E.o[i]
    end
    E, B
end

"""
`uniform_selection!`: Uniformly select particles, ensuring each bin with
positive bin weight has at least one offspring.

### Arguments
* `E` - particle ensemble
* `B` - bin data structure
* `t` - t-th seletion step
### Optional Arguments
* `allocation_resampler=systematic` - resampling scheme amongst bins
* `within_bin_resampler=multinomial` - resampling scheme within bins
"""
function uniform_selection!(E::TE, B::TB, t::Int; allocation_resampler = systematic, within_bin_resampler = multinomial) where {TE<:Ensemble,TB<:Bins}

    # zero out offspring counts
    @. E.o = 0
    @. B.target = 0
    # ensure each nonempty bin has at least one particle
    minimal_bin_allocation!(B)
    n_particles = length(E)
    # number of remaining particles to allocate
    n_allocate = n_particles - sum(B.target)
    try
        # allocate remaining particles
        uniform_bin_allocation!(B, E, n_allocate, allocation_resampler = allocation_resampler)
        # set number of offspring of each particle
        within_bin_allocation!(E, B, within_bin_resampler = within_bin_resampler)
    catch e
        # fall back to trivial allocation if uniform fails
        if e isa DomainError
            @printf("[%d]: TRIVIAL ALLOCATION\n", t)
            trivial_allocation!(E, B)
        else
            rethrow()
        end
    end
    # populate the particles
    repopulate!(E, B)

    E, B
end

"""
`optimal_selection!`: Perform optimal selection of the particles, ensuring each
non empty bin has at least one particle.

### Arguments
* `E` - particle ensemble
* `B` - bin data structure
* `v²` - v² variance function estimator
* `t` - t-th seletion step
### Optional Arguments
* `allocation_resampler=systematic` - resampling scheme amongst bins
* `within_bin_resampler=multinomial` - resampling scheme within bins
"""
function optimal_selection!(E::TE, B::TB, v²::F, t::Int; allocation_resampler = systematic, within_bin_resampler = multinomial) where {TE<:Ensemble,TB<:Bins,F<:Function}

    # zero out offspring counts
    @. E.o = 0
    @. B.target = 0
    # ensure each nonempty bin has at least one particle
    minimal_bin_allocation!(B)
    n_particles = length(E)
    # number of remaining particles to allocate
    n_allocate = n_particles - sum(B.target)
    try
        # allocate remaining particles
        optimal_bin_allocation!(B, E, v², t, n_allocate, allocation_resampler = allocation_resampler);
        # set number of offspring of each particle
        within_bin_allocation!(E, B, within_bin_resampler = within_bin_resampler);
    catch e
        # fall back to trivial allocation if optimal fails
        if e isa DomainError
            @printf("[%d]: TRIVIAL ALLOCATION\n", t)
            trivial_allocation!(E, B)
        else
            rethrow()
        end
    end
    # populate the particles
    repopulate!(E, B);

    E, B
end

"""
`targeted_selection!`: Perform targeted selection of the particles, ensuring each
non empty bin has at least one particle.

### Arguments
* `E` - particle ensemble
* `B` - bin data structure
* `G` - target function
* `t` - t-th seletion step
### Optional Arguments
* `allocation_resampler=systematic` - resampling scheme amongst bins
* `within_bin_resampler=multinomial` - resampling scheme within bins
"""
function targeted_selection!(E::TE, B::TB, G::F, t::Int; allocation_resampler = systematic, within_bin_resampler = multinomial) where {TE<:Ensemble,TB<:Bins,F<:Function}

    # zero out offspring counts
    @. E.o = 0
    @. B.target = 0
    # ensure each nonempty bin has at least one particle
    minimal_bin_allocation!(B)
    n_particles = length(E)
    # number of remaining particles to allocate
    n_allocate = n_particles - sum(B.target)
    try
        # allocate remaining particles
        targeted_bin_allocation!(B, E, G, t, n_allocate, allocation_resampler = allocation_resampler)
        # set number of offspring of each particle
        within_bin_allocation!(E, B, within_bin_resampler=within_bin_resampler)
    catch e
        # fall back to trivial allocation if targeted fails
        if e isa DomainError
            @printf("[%d]: TRIVIAL ALLOCATION\n",t);
            trivial_allocation!(E, B)
        else
            rethrow()
        end
    end
    # populate the particles
    repopulate!(E, B)

    E, B
end

max_stars_repo_name: gideonsimpson/JuWeightedEnsemble

max_issues_repo_name: gideonsimpson/JuWeightedEnsemble

max_forks_repo_name: gideonsimpson/JuWeightedEnsemble


text:
using Citations
using Test

module Demo
    using Citations
    export mul, add

    @declare_citation "Demo Citation Module"

    mul(x,y)=x*y
    @declare_citation mul "-2000BC, Babylonians, Multiplication"

    add(x,y) = x+y
    Citations.cite(::typeof(add)) = "Lost To Time, Addition"

end
using .Demo

@testset "Citations.jl" begin
    @testset "Functions" begin
        @test cite(mul) == "-2000BC, Babylonians, Multiplication"
        @test cite(add) == "Lost To Time, Addition"
    end

    @testset "Modules" begin
        @test cite(Demo) == "Demo Citation Module"
    end

    @testset "Citing julia" begin
        @test match(r"Julia.*Jeff.*Alan.*Stefan.*Viral", cite(Core)) != nothing
    end
end


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..")) # add OceanTurbulenceParameterEstimation to environment stack

using
  Documenter,
  Literate,
  Plots,  # so that Literate.jl does not capture precompilation output
  Distributions,
  OceanTurbulenceParameterEstimation

# Gotta set this environment variable when using the GR run-time on CI machines.
# This happens as examples will use Plots.jl to make plots and movies.
# See: https://github.com/jheinen/GR.jl/issues/278
ENV["GKSwstype"] = "100"

#####
##### Generate examples
#####

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/literated")

examples = [
  # "convective_adjustment_perfect_model_calibration.jl",
  # "convective_adjustment_perfect_model_calibration_uki.jl",
]

for example in examples
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR; flavor = Literate.DocumenterFlavor())
end

#####
##### Build and deploy docs
#####

# Set up a timer to print a space ' ' every 240 seconds. This is to avoid CI
# timing out when building demanding Literate.jl examples.
Timer(t -> println(" "), 0, interval=240)

format = Documenter.HTML(
  collapselevel = 2,
     prettyurls = get(ENV, "CI", nothing) == "true",
      canonical = "https://adelinehillier.github.io/OceanTurbulenceParameterEstimation/dev/",
)

pages = [
    "Home" => "index.md",
    "Installation Instructions" => "installation_instructions.md",
    #=
    "Examples" => [
        "literated/convective_adjustment_perfect_model_calibration.md",
        "literated/convective_adjustment_perfect_model_calibration_uki.md",
        ],
    =#
    "Library" => [
        "Contents"       => "library/outline.md",
        "Public"         => "library/public.md",
        "Private"        => "library/internals.md",
        "Function index" => "library/function_index.md",
        ],
]

makedocs(
   sitename = "OceanTurbulenceParameterEstimation.jl",
    modules = [OceanTurbulenceParameterEstimation],
     format = format,
      pages = pages,
    doctest = true,
     strict = true,
      clean = true,
  checkdocs = :exports
)

deploydocs(        repo = "github.com/adelinehillier/OceanTurbulenceParameterEstimation.jl.git",
               versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
              devbranch = "main",
           push_preview = true
)


text:
module NonLinearFEOperatorsTests

using Test
using Gridap

import Gridap: ∇

# Define manufactured functions
ufun(x) = x[1] + x[2]
ufun_grad(x) = VectorValue(1.0,1.0)
∇(::typeof(ufun)) = ufun_grad
bfun(x) = -(3.0*x[1]+x[2]+1.0)

@law ν(x,u) = (u+1.0)*x[1]
@law dν(x,du) = x[1]*du

# Construct the discrete model
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))

# Construct the FEspace
order = 2
diritag = "boundary"
fespace = H1ConformingFESpace(Float64,model,order,diritag)

# Define test and trial
V = TestFESpace(fespace)
U = TrialFESpace(fespace,ufun)

# Define integration mesh and quadrature
trian = Triangulation(model)
quad = CellQuadrature(trian,degree=3*order-1)

# Define cell field describing the source term
bfield = CellField(trian,bfun)

# Define residual and jacobian
a(u,v,du) = inner( ∇(v), ν(u)*∇(du))
res(u,v) = a(u,v,u) - inner(v,bfield)
jac(u,v,du) = a(u,v,du) + inner(∇(v),dν(du)*∇(u))

# Define Assembler
assem = SparseMatrixAssembler(V,U)

# Define the FEOperator
op = NonLinearFEOperator(res,jac,V,U,assem,trian,quad)

# Define the FESolver
ls = LUSolver()
nls = NLSolver(
  ls,method=:newton,show_trace=false,ftol=1.e-10,iterations=20)
solver = NonLinearFESolver(nls)

# Solve!
uh = solve(solver,op)

# Define exact solution and error
u = CellField(trian,ufun)
e = u - uh

# Define norms to measure the error
l2(u) = inner(u,u)
sh1(u) = inner(∇(u),∇(u))
h1(u) = sh1(u) + l2(u)

# Compute errors
el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))
ul2 = sqrt(sum( integrate(l2(u),trian,quad) ))
uh1 = sqrt(sum( integrate(h1(u),trian,quad) ))

@test el2/ul2 < 1.e-8
@test eh1/uh1 < 1.e-7

#writevtk(trian,"trian",nref=4,cellfields=["uh"=>uh,"u"=>u,"e"=>e])

# Further tests
@test TrialFESpace(op) === U
@test TestFESpace(op) === V

zh = zero(U)
zh, cache = solve!(zh,solver,op)
@test free_dofs(zh) ≈ free_dofs(uh)

zh = zero(U)
zh, solve!(zh,solver,op,cache)
@test free_dofs(zh) ≈ free_dofs(uh)

t_Ω = NonLinearFETerm(res,jac,trian,quad)

b0fun(x) = 0.0
b0(v) = inner(v,b0fun)
t_source = FESource(b0,trian,quad)

op = NonLinearFEOperator(V,U,assem,t_Ω,t_source)
op = NonLinearFEOperator(V,U,t_Ω,t_source)

end # module NonLinearFEOperatorsTests


text:
module SemifieldAlgebra

export Semifield

"""
    struct Semifield{T,A,M,MI,Z,O} <: Number
        val::T
    end

`T` is the value type, `A` is the addition function, `M` is the
multiplication function, `MI` is the muliplication inverse function,
`Z` is the zero element function and `O` is the one element function.
"""
struct Semifield{T,A,M,MI,Z,O} <: Number
    val::T
end

Base.show(io::IO, x::Semifield) = print(io, x.val)
Base.:+(x::Semifield{T,A,M,MI,Z,O}, y::Semifield{T,A,M,MI,Z,O}) where {T,A,M,MI,Z,O} =
    Semifield{T,A,M,MI,Z,O}( A(x.val, y.val) )
Base.:*(x::Semifield{T,A,M,MI,Z,O}, y::Semifield{T,A,M,MI,Z,O}) where {T,A,M,MI,Z,O} =
    Semifield{T,A,M,MI,Z,O}( M(x.val, y.val) )
Base.:/(x::Semifield{T,A,M,MI,Z,O}, y::Semifield{T,A,M,MI,Z,O}) where {T,A,M,MI,Z,O} =
    Semifield{T,A,M,MI,Z,O}( MI(x.val, y.val) )
Base.zero(::Type{Semifield{T,A,M,MI,Z,O}}) where {T,A,M,MI,Z,O} = Semifield{T,A,M,MI,Z,O}(Z)
Base.zero(::Semifield{T,A,M,MI,Z,O}) where {T,A,M,MI,Z,O} = Semifield{T,A,M,MI,Z,O}(Z)
Base.one(::Type{Semifield{T,A,M,MI,Z,O}})  where {T,A,M,MI,Z,O} = Semifield{T,A,M,MI,Z,O}(O)
Base.one(::Semifield{T,A,M,MI,Z,O}) where {T,A,M,MI,Z,O} = Semifield{T,A,M,MI,Z,O}(O)

end


text:
################################################################
# Array index functions
# Not part of P4est, but really usefull to access sc_array from julia
################################################################

function sc_array_index(sc_array_object::sc_array_t, iz)
    @assert iz < sc_array_object.elem_count
    return sc_array_object.array + sc_array_object.elem_size * iz
end

function p4est_tree_array_index(sc_array_object::sc_array_t, it)
    @assert sc_array_object.elem_size == sizeof(p4est_tree_t)
    @assert it in 0:sc_array_object.elem_count
    return Ptr{p4est_tree_t}(sc_array_object.array + sc_array_object.elem_size*it)
end

function p4est_quadrant_array_index(sc_array_object::sc_array_t, it)
    @assert sc_array_object.elem_size == sizeof(p4est_quadrant_t)
    @assert it < sc_array_object.elem_count
    return Ptr{p4est_quadrant_t}(sc_array_object.array + sc_array_object.elem_size*it)
end

function p8est_tree_array_index(sc_array_object::sc_array_t, it)
  @assert sc_array_object.elem_size == sizeof(p8est_tree_t)
  @assert it in 0:sc_array_object.elem_count
  return Ptr{p8est_tree_t}(sc_array_object.array + sc_array_object.elem_size*it)
end

function p8est_quadrant_array_index(sc_array_object::sc_array_t, it)
  @assert sc_array_object.elem_size == sizeof(p8est_quadrant_t)
  @assert it < sc_array_object.elem_count
  return Ptr{p8est_quadrant_t}(sc_array_object.array + sc_array_object.elem_size*it)
end

################################################################
# Array index functions
# Static functions are not part of the API/ABI
################################################################

function sc_array_index(array::Ptr{sc_array_t}, iz)
    sc_array_object = unsafe_wrap(Array, array, 1)[1]
    return sc_array_index(sc_array_object, iz)
end

function p4est_tree_array_index(array::Ptr{sc_array_t}, it)
    sc_array_object = unsafe_wrap(Array, array, 1)[1]
    return p4est_tree_array_index(sc_array_object, it)
end

function p4est_quadrant_array_index(array::Ptr{sc_array_t}, it)
    sc_array_object = unsafe_wrap(Array, array, 1)[1]
    return p4est_quadrant_array_index(sc_array_object, it)
end

function p8est_tree_array_index(array::Ptr{sc_array_t}, it)
  sc_array_object = unsafe_wrap(Array, array, 1)[1]
  return p8est_tree_array_index(sc_array_object, it)
end

function p8est_quadrant_array_index(array::Ptr{sc_array_t}, it)
  sc_array_object = unsafe_wrap(Array, array, 1)[1]
  return p8est_quadrant_array_index(sc_array_object, it)
end

################################################################
# p4est_init function
# This function is part of p4est_base.h header.
# This header is strongly related with p4est_config.h
# and sc_config.h, so not added in automatic generation
# We wrap this binding manually.
################################################################
function p4est_init(log_handler, log_threshold)
    ccall((:p4est_init, p4est_lib), Cvoid, (sc_log_handler_t, Cint), log_handler, log_threshold)
end

################################################################
# Dereference functions
# Not part of P4est, but really usefull to access data from pointers
################################################################

PXEST_NAMED_STRUCTS = Union{p4est_quadrant, p4est_tree, p4est_inspect, p4est_connectivity, p4est, p4est_geometry,
                            p4est_ghost_t, p4est_ghost_exchange, p4est_iter_volume_info, p4est_iter_face_side,
                            p4est_iter_face_info, p4est_iter_corner_side, p4est_iter_corner_info, p4est_lnodes,
                            p4est_lnodes_rank, p4est_lnodes_buffer, p4est_search_local_t, p4est_vtk_context,
                            p6est_connectivity, p2est_quadrant, p6est_ghost, p8est_quadrant, p8est_tree, sc_array,
                            p8est_inspect, p8est_connectivity, p8est, p8est_geometry, p8est_ghost_t,
                            p8est_ghost_exchange, p8est_iter_volume_info, p8est_iter_face_side, p8est_iter_face_info,
                            p8est_iter_edge_side, p8est_iter_edge_info, p8est_iter_corner_side, p8est_iter_corner_info,
                            p8est_lnodes, p8est_lnodes_rank, p8est_lnodes_buffer, p8est_search_local_t, p8est_vtk_context}

function Base.getindex(ptr::Ptr{T}) where {T<:PXEST_NAMED_STRUCTS}
    unsafe_wrap(Array, ptr, 1)[]
end

function Base.getindex(ptr::Ptr{T}, i) where {T<:PXEST_NAMED_STRUCTS}
    unsafe_wrap(Array, ptr, i)[i]
end


text:
"""
	vcongruence(
        V::Lar.Points; epsilon=1e-6
    )::Tuple{Lar.Points, Array{Array{Int,1},1}}

Evaluates the Vertex Congruence for 3D-points ``V ∈ ℳ(3,n)``.

The function determines the points of ``V`` closer than ``epsilon`` and
builds a new Vertex Set made of the representative of each point cluster.

The method returns:
 - the new Vertex Set
 - a map that, for every new vertex, identifies the old vertices it is made of
"""
function vcongruence(V::Lar.Points; epsilon=1e-6)
    vclasses, visited = [], []
    V = convert(Array{Float64,2}, V)
    kdtree = NearestNeighbors.KDTree(V)
    for vidx = 1 : size(V, 2)  if !(vidx in visited)
        nearvs = NearestNeighbors.inrange(kdtree, V[:,vidx], epsilon)
        push!(vclasses, nearvs)
        append!(visited, nearvs)
    end  end
    W = hcat([sum(V[:,class], dims=2)/length(class) for class in vclasses]...)
    return W, vclasses
end

"""
    cellCongruence(
        Delta::Lar.ChainOp,
        inclasses::Array{Array{Int,1},1};
        dim::Int = 0
    )::Tuple{Lar.Cells, Array{Array{Int,1},1}}

Evaluates the Cell Congruence for Cochain ``Delta`` with classes ``inclasses``.

The function determines the new Cochain (as ``Lar.Cells``) built from ``cop``
where the lower order cells are merged according to ``inclasses`` map.

The parameter ``dim``, if specified, represent the order of the cells
(that is the least number of lower order cells a current order cell is made of).

The method returns:
- the new Cochain Operator as ``Lar.Cells``
- a map that, for every new cell, identifies the old cells it is made of
"""

function cellcongruence(Delta, inclasses; dim = 0)
    cellarray = Lar.cop2lar(Delta)
    newcell = Vector(undef, size(Delta,2))
    for (k, class) in enumerate(inclasses)  for e in class
            newcell[e] = k
    end end
    cells = [sort(map(x -> newcell[x], face)) for face in cellarray]
    okcells = [cell for cell in cells if length(Set(cell)) > dim] # non-empty cells
    classes = DefaultOrderedDict{Vector, Vector}([])
    for (k,face) in enumerate(okcells)
        classes[face] == [] ? classes[face] = [k] : append!(classes[face], [k])
    end
    cells = collect(keys(classes))
    outclasses = collect(values(classes))
    return cells, outclasses
end

"""
    chaincongruence(
        W::Lar.Points,
        Delta_0::Lar.ChainOp,
        Delta_1::Lar.ChainOp
        epsilon = 1e-6
    )::Tuple{Lar.Points, Lar.Cells, Lar.Cells}

Cell Congruence Enabling Algorithm for 3D Complexes.

The algorithm performs the identify vertices within ``epsilon`` and coherently
update the Chain Operators.
"""
function chaincongruence(W, Delta_0, Delta_1; epsilon = 1e-6)
    V, vclasses = vcongruence(W, epsilon = epsilon)
    EV, eclasses = cellcongruence(Delta_0, vclasses, dim=1)
    FE, _ = cellcongruence(Delta_1, eclasses, dim=2)
    return V, EV, FE
end

"""
	signedCellCongruence(
		cop::Lar.ChainOp,
		lo_cls::Array{Array{Int,1},1},
		lo_sign::Array{Array{Int8,1},1};
		imp = false,
		d = 0
	)
Evaluates the Cell Congruence for a Cochain ``cop`` with classes ``lo_cls``.
The function determines the new Cochain Operator built from ``cop`` where
the lower order cells are merged according to ``lo_cls`` map.
``lo_sign`` specifies whether a cell must be considered in reverse order.
If optional paramenter ``imp`` is set to ``true`` then FP imprecisions
are taken into account in the sense that lower order cells may have collided.
The parameter ``d`` represent then the order of the cell (that also is the
least number of lower order cells a current order cell is made of).
The method returns:
 - the new Cochain Operator
 - a map that, for every new cell, identifies the old cells it is made of
 - a map that, for every new cell, specify if old cells have changed ordering.
"""
function signedCellCongruence(cop, lo_cls, lo_sign; imp = false, dim = 0)::Tuple{
        Lar.ChainOp,
		Array{Array{Int,1},1},
		Array{Array{Int8,1},1}
    }
#function cellCongruence(
#		cop::Lar.ChainOp,
#		lo_cls::Array{Array{Int,1},1},
#		lo_sign::Array{Array{Int8,1},1};
#		imp = false,
#		dim = 0
#	)::Tuple{
#		Lar.ChainOp,
#		Array{Array{Int,1},1},
#		Array{Array{Int8,1},1}
#	}

	# Collide columns
	copCols = []
	for i = 1 : length(lo_cls)
		col = sum([
			cop[:, lo_cls[i][j]] .* lo_sign[i][j]
			for j = 1 : length(lo_cls[i])
		])
#		if imp
#			#TODO remove zeros
#			if length(col.nzind) > dim  push!(copCols, col)  end
#		end
		push!(copCols, col)
	end

	# Retrieve Matrix Form and extract rows
	cop = hcat(copCols...)
	rows = [cop[row, :] for row = 1 : cop.m]

	# Adjustinfg signs such that the first value always is `-1`
	sign = ones(Int8, cop.m)
	for row = 1 : cop.m  if rows[row].nzval[1] > 0
		rows[row] = -rows[row]
		sign[row] = -1
	end  end

	# Sort rows in order to collide them quickly
	rows_ord = sortperm(rows)
	nrows = unique(rows[rows_ord])
	ho_cls = [Array{Int,1}() for i = 1 : length(nrows)]
	nidx = 1

	# Collide cells with same structure
	for cidx = 1 : cop.m
		if rows[rows_ord[cidx]] != nrows[nidx]  nidx = nidx + 1;  end
		push!(ho_cls[nidx], rows_ord[cidx])
	end

	# Reshaping Sign and Cochain
	ho_sign = [[sign[el] for el in cell] for cell in ho_cls]
#	cop = convert(Lar.ChainOp, hcat(nrows...)')

	return hcat(nrows...)', ho_cls, ho_sign
end


text:
m = page_model()
include("../src/components/CH4emissions.jl")

add_comp!(m, ch4emissions)

set_param!(m, :ch4emissions, :e0_baselineCH4emissions, readpagedata(m, "data/e0_baselineCH4emissions.csv")) #PAGE 2009 documentation pp38
set_param!(m, :ch4emissions, :er_CH4emissionsgrowth, readpagedata(m, "data/er_CH4emissionsgrowth.csv"))

##running Model
run(m)

# Generated data
emissions= m[:ch4emissions,  :e_regionalCH4emissions]

# Recorded data
emissions_compare=readpagedata(m, "test/validationdata/e_regionalCH4emissions.csv")

@test emissions ≈ emissions_compare rtol=1e-3


text:
# # Bukdu v0.4.5
using Bukdu # ApplicationController Conn HTML Router CLI render routes get post
using Bukdu.HTML5.Form # change
                       # form_for label_for
                       # text_area text_input radio_button checkbox file_input
                       # submit
using Documenter.Utilities.DOM: @tags
using HTTP: Multipart
using Base64

struct FormController <: ApplicationController
    conn::Conn
end

function layout(body)
    """

  <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />

Bukdu.HTML5.Form

$body

"""
end

function preview(multipart::Multipart)
    if eof(multipart)
        string("empty")
    else
        @tags div img
        mark(multipart)
        data = read(multipart)
        reset(multipart)
        contents = [
            div(multipart.filename),
            div(repr(length(data)), " bytes"),
        ]
        if startswith(multipart.contenttype, "image")
            base64 = base64encode(data)
            push!(contents, img[:src => string("data:", multipart.contenttype, ";base64, ", base64),
                                :alt => multipart.filename])
        end
        div(contents)
    end
end

function preview(::Nothing)
    ""
end

struct User
end

global changeset = Changeset(User, (name="", famous=false, season="summer", intro="", photo1=nothing, photo2=nothing))

function index(c::FormController)
    global changeset
    @tags div
    form1 = form_for(changeset, (FormController, post_result), method=post, multipart=true) do f
        div(
            div.([
                text_area(f, :name, placeholder="Name", rows="3", cols="50"),
                label_for(checkbox(f, :famous), "Famous"),
                label_for(file_input(f, :photo1), "Photo1"),
                preview(changeset.changes.photo1),
                label_for(file_input(f, :photo2), "Photo2"),
                preview(changeset.changes.photo2),
                submit("Submit"),
            ])
        )
    end
    body = div(form1)
    render(HTML, layout(body))
end

function post_result(c::FormController)
    global changeset
    @tags div a h3 li p
    result = change(changeset, c.params)
    if !isempty(result.changes)
        changeset.changes = merge(changeset.changes, result.changes)
    end
    body = div(
        h3( a[:href => "/"]("back") ),
        p(isempty(result.changes) ?
            "nothing's changed" :
            "✅ changed"),
        li( "result: ", string(result) ),
        li( "global changeset: ", string(changeset) ),
        preview(changeset.changes.photo1),
        preview(changeset.changes.photo2),
    )
    render(HTML, layout(body))
end

if PROGRAM_FILE == basename(@__FILE__)

routes() do
    get("/", FormController, index)
    post("/post_result", FormController, post_result)
end

Bukdu.start(8080)

Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if


text:
## Calculus

Base.sum(f::Fun{JacobiWeight{C,DD,RR,TT}}) where {C<:Chebyshev,DD,RR,TT} =
    sum(setdomain(f,canonicaldomain(f))*fromcanonicalD(f))
linesum(f::Fun{JacobiWeight{C,DD,RR,TT}}) where {C<:Chebyshev,DD,RR,TT} =
    linesum(setdomain(f,canonicaldomain(f))*abs(fromcanonicalD(f)))

for (Func,Len) in ((:(Base.sum),:complexlength),(:linesum,:arclength))
    @eval begin
        function $Func(f::Fun{JacobiWeight{C,DD,RR,TT}}) where {C<:Chebyshev,DD<:IntervalOrSegment,RR,TT}
            tol=1e-10
            d,β,α,n=domain(f),f.space.β,f.space.α,ncoefficients(f)
            n == 0 && return zero(cfstype(f))*$Len(d)
            g=Fun(space(f).space,f.coefficients)
            if β ≤ -1.0 && abs(first(g))≤tol
                $Func(increase_jacobi_parameter(-1,f))
            elseif α ≤ -1.0 && abs(last(g))≤tol
                $Func(increase_jacobi_parameter(+1,f))
            elseif β ≤ -1.0 || α ≤ -1.0
                fs = Fun(f.space.space,f.coefficients)
                return Inf*0.5*$Len(d)*(sign(fs(leftendpoint(d)))+sign(fs(rightendpoint(d))))/2
            elseif β == α == -0.5
                return 0.5*$Len(d)*f.coefficients[1]*π
            elseif β == α == 0.5
                return 0.5*$Len(d)*(n ≤ 2 ? f.coefficients[1]/2 : f.coefficients[1]/2 - f.coefficients[3]/4)*π
            elseif β == 0.5 && α == -0.5
                return 0.5*$Len(d)*(n == 1 ? f.coefficients[1] : f.coefficients[1] + f.coefficients[2]/2)*π
            elseif β == -0.5 && α == 0.5
                return 0.5*$Len(d)*(n == 1 ? f.coefficients[1] : f.coefficients[1] - f.coefficients[2]/2)*π
            else
                c = zeros(cfstype(f),n)
                c[1] = 2^(β+α+1)*gamma(β+1)*gamma(α+1)/gamma(β+α+2)
                if n > 1
                    c[2] = c[1]*(β-α)/(β+α+2)
                    for i=1:n-2
                        c[i+2] = (2(β-α)*c[i+1]-(β+α-i+2)*c[i])/(β+α+i+2)
                    end
                end
                return 0.5*$Len(d)*dotu(f.coefficients,c)
            end
        end
        $Func(f::Fun{JacobiWeight{PS,DD,RR,TT}}) where {PS<:PolynomialSpace,DD,RR,TT} =
            $Func(Fun(f,JacobiWeight(space(f).β, space(f).α, Chebyshev(domain(f)))))
    end
end

function differentiate(f::Fun{JacobiWeight{SS,DD,RR,TT}}) where {SS,DD<:IntervalOrSegment,RR,TT}
    S=f.space
    d=domain(f)
    ff=Fun(S.space,f.coefficients)
    if S.β==S.α==0
        u=differentiate(ff)
        Fun(JacobiWeight(0.,0.,space(u)),u.coefficients)
    elseif S.β==0
        x=Fun(identity,d)
        M=tocanonical(d,x)
        Mp=tocanonicalD(d,leftendpoint(d))
        u=-Mp*S.α*ff +(1-M).*differentiate(ff)
        Fun(JacobiWeight(0.,S.α-1,space(u)),u.coefficients)
    elseif S.α==0
        x=Fun(identity,d)
        M=tocanonical(d,x)
        Mp=tocanonicalD(d,leftendpoint(d))
        u=Mp*S.β*ff +(1+M).*differentiate(ff)
        Fun(JacobiWeight(S.β-1,0.,space(u)),u.coefficients)
    else
        x=Fun(identity,d)
        M=tocanonical(d,x)
        Mp=tocanonicalD(d,leftendpoint(d))
        u=(Mp*S.β)*(1-M).*ff- (Mp*S.α)*(1+M).*ff +(1-M.^2).*differentiate(ff)
        Fun(JacobiWeight(S.β-1,S.α-1,space(u)),u.coefficients)
    end
end

function integrate(f::Fun{JacobiWeight{SS,DD,RR,TT}}) where {SS,DD<:IntervalOrSegment,RR,TT}
    S=space(f)
    # we integrate by solving u'=f
    tol=1e-10
    g=Fun(S.space,f.coefficients)
    if S.β ≈ 0 && S.α ≈ 0
        integrate(g)
    elseif S.β ≤ -1 && abs(first(g)) ≤ tol
        integrate(increase_jacobi_parameter(-1,f))
    elseif S.α ≤ -1 && abs(last(g)) ≤ tol
        integrate(increase_jacobi_parameter(+1,f))
    elseif S.β ≈ -1 && S.α ≈ -1
        error("Implement")
    elseif S.β ≈ -1 && S.α ≈ 0
        p=first(g)  # first value without weight
        fp = Fun(f-Fun(S,[p]),S.space)  # Subtract out right value and divide singularity via conversion
        d=domain(f)
        Mp=tocanonicalD(d,leftendpoint(d))
        integrate(fp) ⊕ Fun(LogWeight(1.,0.,S.space),[p/Mp])
    elseif S.β ≈ -1 && S.α > 0 && isapproxinteger(S.α)
        # convert to zero case and integrate
        integrate(Fun(f,JacobiWeight(S.β,0.,S.space)))
    elseif S.α ≈ -1 && S.β ≈ 0
        p=last(g)  # last value without weight
        fp = Fun(f-Fun(S,[p]),S.space)  # Subtract out right value and divide singularity via conversion
        d=domain(f)
        Mp=tocanonicalD(d,leftendpoint(d))
        integrate(fp) ⊕ Fun(LogWeight(zero(TT),one(TT),S.space),[-p/Mp])
    elseif isapprox(S.α,-1) && S.β > 0 && isapproxinteger(S.β)
        # convert to zero case and integrate
        integrate(Fun(f,JacobiWeight(zero(TT),S.α,S.space)))
    elseif S.β ≈ 0
        D = Derivative(JacobiWeight(S.β, S.α+1, S.space))
        D\f   # this happens to pick out a smooth solution
    elseif S.α ≈ 0
        D = Derivative(JacobiWeight(S.β+1, S.α, S.space))
        D\f   # this happens to pick out a smooth solution
    elseif isapproxinteger(S.β) || isapproxinteger(S.α)
        D = Derivative(JacobiWeight(S.β+1, S.α+1, S.space))
        D\f   # this happens to pick out a smooth solution
    else
        s=sum(f)
        if abs(s)<1E-14
            D=Derivative(JacobiWeight(S.β+1, S.α+1, S.space))
            \(D,f; tolerance=1E-14)  # if the sum is 0 we don't get step-like behaviour
        else
            # we normalized so it sums to zero, and so backslash works
            w = Fun(x->exp(-40x^2),81)
            w1 = Fun(S,coefficients(w))
            w2 = Fun(x->w1(x),domain(w1))
            c  = s/sum(w1)
            v  = f-w1*c
            (c*integrate(w2)) ⊕ integrate(v)
        end
    end
end

function Base.cumsum(f::Fun{JacobiWeight{SS,DD,RR,TT}}) where {SS,DD<:IntervalOrSegmentDomain,RR,TT}
    g=integrate(f)
    S=space(f)

    if (S.β==0 && S.α==0) || S.β>-1
        g-first(g)
    else
        @warn "Function is not integrable at left endpoint.  Returning a non-normalized indefinite integral."
        g
    end
end

## Operators

function jacobiweightDerivative(S::JacobiWeight{<:Any,<:IntervalOrSegment})
    d = domain(S)
    # map to canonical
    Mp=fromcanonicalD(d,leftendpoint(d))
    DD=jacobiweightDerivative(setdomain(S,ChebyshevInterval()))

    return DerivativeWrapper(SpaceOperator(DD.op.op,S,setdomain(rangespace(DD),d))/Mp,1)
end

function jacobiweightDerivative(S::JacobiWeight{<:Any,<:ChebyshevInterval})
    d=domain(S)

    if S.β==S.α==0
        DerivativeWrapper(SpaceOperator(Derivative(S.space),S,JacobiWeight(0.,0.,rangespace(Derivative(S.space)))),1)
    elseif S.β==0
        w=Fun(JacobiWeight(0,1,ConstantSpace(d)),[1.])

        DD=-S.α + w*Derivative(S.space)
        rs=S.α==1 ? rangespace(DD) : JacobiWeight(0.,S.α-1,rangespace(DD))
        DerivativeWrapper(SpaceOperator(DD,S,rs),1)
    elseif S.α==0
        w=Fun(JacobiWeight(1,0,ConstantSpace(d)),[1.])

        DD=S.β + w*Derivative(S.space)
        rs=S.β==1 ? rangespace(DD) : JacobiWeight(S.β-1,0.,rangespace(DD))
        DerivativeWrapper(SpaceOperator(DD,S,rs),1)
    else
        w=Fun(JacobiWeight(1,1,ConstantSpace(d)),[1.])
        x=Fun()

        DD=S.β*(1-x) - S.α*(1+x) + w*Derivative(S.space)
        rs=S.β==1&&S.α==1 ? rangespace(DD) : JacobiWeight(S.β-1,S.α-1,rangespace(DD))
        DerivativeWrapper(SpaceOperator(DD,S,rs),1)
    end
end

Derivative(S::JacobiWeight{SS,DDD}) where {SS,DDD<:IntervalOrSegment} = jacobiweightDerivative(S)

function Derivative(S::JacobiWeight{SS,DD}, k::Integer) where {SS,DD<:IntervalOrSegment}
    if k==1
        Derivative(S)
    else
        D=Derivative(S)
        DerivativeWrapper(TimesOperator(Derivative(rangespace(D),k-1),D),k)
    end
end

## Multiplication

#Left multiplication. Here, S is considered the domainspace and we determine rangespace accordingly.

function Multiplication(f::Fun{<:JacobiWeight},S::JacobiWeight)
    M=Multiplication(Fun(space(f).space,f.coefficients),S.space)
    if space(f).β+S.β==space(f).α+S.α==0
        rsp=rangespace(M)
    else
        rsp=JacobiWeight(space(f).β+S.β,space(f).α+S.α,rangespace(M))
    end
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

function Multiplication(f::Fun, S::JacobiWeight)
    M=Multiplication(f,S.space)
    rsp=JacobiWeight(S.β,S.α,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

function Multiplication(f::Fun{<:JacobiWeight},S::PolynomialSpace)
    M=Multiplication(Fun(space(f).space,f.coefficients),S)
    rsp=JacobiWeight(space(f).β,space(f).α,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

#Right multiplication. Here, S is considered the rangespace and we determine domainspace accordingly.

function Multiplication(S::JacobiWeight, f::Fun{<:JacobiWeight})
    M=Multiplication(Fun(space(f).space,f.coefficients),S.space)
    dsp=canonicalspace(JacobiWeight(S.β-space(f).β,S.α-space(f).α,rangespace(M)))
    MultiplicationWrapper(f,SpaceOperator(M,dsp,S))
end

function Multiplication(S::JacobiWeight, f::Fun)
    M=Multiplication(f,S.space)
    dsp=JacobiWeight(S.β,S.α,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,dsp,S))
end

# function Multiplication{D<:JacobiWeight,T,V,ID<:IntervalOrSegment}(S::Space{V,D},f::Fun{ID,T})
#     M=Multiplication(Fun(f.coefficients,space(f).space),S)
#     dsp=JacobiWeight(-space(f).β,-space(f).α,rangespace(M))
#     MultiplicationWrapper(f,SpaceOperator(M,dsp,S))
# end

## Conversion
for (OPrule,OP) in ((:maxspace_rule,:maxspace),(:union_rule,:union))
    @eval begin
        function $OPrule(A::JacobiWeight,B::JacobiWeight)
            if domainscompatible(A,B) && isapproxinteger(A.β-B.β) && isapproxinteger(A.α-B.α)
                ms=$OP(A.space,B.space)
                if min(A.β,B.β)==0.0 && min(A.α,B.α) == 0.0
                    return ms
                else
                    return JacobiWeight(min(A.β,B.β),min(A.α,B.α),ms)
                end
            end
            NoSpace()
        end
        $OPrule(A::JacobiWeight,B::Space{D}) where {D<:IntervalOrSegmentDomain} = $OPrule(A,JacobiWeight(0.,0.,B))
    end
end

for FUNC in (:hasconversion,:isconvertible)
    @eval begin
        $FUNC(A::JacobiWeight{S1,D},B::JacobiWeight{S2,D}) where {S1,S2,D<:IntervalOrSegmentDomain} =
            isapproxinteger(A.β-B.β) &&
            isapproxinteger(A.α-B.α) && A.β ≥ B.β && A.α ≥ B.α && $FUNC(A.space,B.space)

        $FUNC(A::JacobiWeight{S,D},B::Space{D}) where {S,D<:IntervalOrSegmentDomain} =
            $FUNC(A,JacobiWeight(0.,0.,B))
        $FUNC(B::Space{D},A::JacobiWeight{S,D}) where {S,D<:IntervalOrSegmentDomain} =
            $FUNC(JacobiWeight(0.,0.,B),A)
    end
end

# return the space that has banded Conversion to the other, or NoSpace

function conversion_rule(A::JacobiWeight,B::JacobiWeight)
    if isapproxinteger(A.β-B.β) && isapproxinteger(A.α-B.α)
        ct=conversion_type(A.space,B.space)
        ct==NoSpace() ? NoSpace() : JacobiWeight(max(A.β,B.β),max(A.α,B.α),ct)
    else
        NoSpace()
    end
end

conversion_rule(A::JacobiWeight,B::Space{D}) where {D<:IntervalOrSegmentDomain} = conversion_type(A,JacobiWeight(0,0,B))

# override defaultConversion instead of Conversion to avoid ambiguity errors

defaultConversion(A::Space{<:IntervalOrSegmentDomain,<:Real},B::JacobiWeight{<:Any,<:IntervalOrSegmentDomain}) =
    ConversionWrapper(SpaceOperator(
        Conversion(JacobiWeight(0,0,A),B),
        A,B))
defaultConversion(A::JacobiWeight{<:Any,<:IntervalOrSegmentDomain},B::Space{<:IntervalOrSegmentDomain,<:Real}) =
    ConversionWrapper(SpaceOperator(
        Conversion(A,JacobiWeight(0,0,B)),
        A,B))

function defaultConversion(A::JacobiWeight{<:Any,<:IntervalOrSegmentDomain},B::JacobiWeight{<:Any,<:IntervalOrSegmentDomain})
    if isapprox(A.β,B.β) && isapprox(A.α,B.α)
        ConversionWrapper(SpaceOperator(Conversion(A.space,B.space),A,B))
    elseif isapprox(A.β-B.β, A.space.b-B.space.b) && isapprox(A.α-B.α, A.space.a-B.space.a)
        ConversionWrapper(SpaceOperator(Multiplication(jacobiweight(A.β-B.β,A.α-B.α,domain(A)),A.space),A,B))
    else
        C=JacobiWeight(A.β,A.α,Jacobi(B.space.b+A.β-B.β,B.space.a+A.α-B.α,domain(A)))
        ConversionWrapper(Conversion(C,B)*Conversion(A,C))
    end
end

## Evaluation

function  Base.getindex(op::ConcreteEvaluation{<:JacobiWeight,typeof(leftendpoint)},kr::AbstractRange)
    S=op.space
    @assert op.order ≤ 1
    d=domain(op)

    @assert S.β ≥ 0
    if S.β==0
        if op.order==0
            2^S.α*getindex(Evaluation(S.space,op.x),kr)
        else #op.order ===1
            @assert isa(d,IntervalOrSegment)
            2^S.α*getindex(Evaluation(S.space,op.x,1),kr)-(tocanonicalD(d,leftendpoint(d))*S.α*2^(S.α-1))*getindex(Evaluation(S.space,op.x),kr)
        end
    else
        @assert op.order==0
        zeros(eltype(op), length(kr))
    end
end

function  Base.getindex(op::ConcreteEvaluation{<:JacobiWeight,typeof(rightendpoint)},kr::AbstractRange)
    S=op.space
    @assert op.order<=1
    d=domain(op)

    @assert S.α>=0
    if S.α==0
        if op.order==0
            2^S.β*getindex(Evaluation(S.space,op.x),kr)
        else #op.order ===1
            @assert isa(d,IntervalOrSegment)
            2^S.β*getindex(Evaluation(S.space,op.x,1),kr)+
                (tocanonicalD(d,leftendpoint(d))*S.β*2^(S.β-1))*getindex(Evaluation(S.space,op.x),kr)
        end
    else
        @assert op.order==0
        zeros(eltype(op), length(kr))
    end
end

## Definite Integral

for (Func,Len,Sum) in ((:DefiniteIntegral,:complexlength,:sum),(:DefiniteLineIntegral,:arclength,:linesum))
    ConcFunc = Meta.parse("Concrete"*string(Func))

    @eval begin
        $Func(S::JacobiWeight{SS,D}) where {SS,D<:IntervalOrSegment} = $ConcFunc(S)

        getindex(Σ::$ConcFunc,k::Integer) = eltype(Σ)($Sum(Fun(domainspace(Σ),[zeros(eltype(Σ),k-1);1])))

        function getindex(Σ::$ConcFunc{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT},T},k::Integer) where {LT,D<:IntervalOrSegment,R,T,TT}
            λ = order(domainspace(Σ).space)
            dsp = domainspace(Σ)
            d = domain(Σ)
            C = $Len(d)/2

            if dsp.β==dsp.α==λ-0.5
                k == 1 ? convert(T,C*gamma(λ+one(T)/2)*gamma(one(T)/2)/gamma(λ+one(T))) : zero(T)
            else
                convert(T,$Sum(Fun(dsp,[zeros(T,k-1);1])))
            end
        end

        function getindex(Σ::$ConcFunc{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT},T},kr::AbstractRange) where {LT,D<:IntervalOrSegment,R,T,TT}
            λ = order(domainspace(Σ).space)
            dsp = domainspace(Σ)
            d = domain(Σ)
            C = $Len(d)/2

            if dsp.β==dsp.α==λ-0.5
                T[k == 1 ? C*gamma(λ+one(T)/2)*gamma(one(T)/2)/gamma(λ+one(T)) : zero(T) for k=kr]
            else
                T[$Sum(Fun(dsp,[zeros(T,k-1);1])) for k=kr]
            end
        end

        function bandwidths(Σ::$ConcFunc{JacobiWeight{Ultraspherical{LT,D,R},D,R,TT}}) where {LT,D<:IntervalOrSegment,R,TT}
            λ = order(domainspace(Σ).space)
            β,α = domainspace(Σ).β,domainspace(Σ).α
            if β==α && isapproxinteger(β-0.5-λ) && λ ≤ ceil(Int,β)
                0,2*(ceil(Int,β)-λ)
            else
                0,∞
            end
        end

        function getindex(Σ::$ConcFunc{JacobiWeight{Chebyshev{D,R},D,R,TT},T},k::Integer) where {D<:IntervalOrSegment,R,T,TT}
            dsp = domainspace(Σ)
            d = domain(Σ)
            C = $Len(d)/2

            if dsp.β==dsp.α==-0.5
                k == 1 ? convert(T,C*π) : zero(T)
            else
                convert(T,$Sum(Fun(dsp,[zeros(T,k-1);1])))
            end
        end

        function getindex(Σ::$ConcFunc{JacobiWeight{Chebyshev{D,R},D,R,TT},T},kr::AbstractRange) where {D<:IntervalOrSegment,R,T,TT}
            dsp = domainspace(Σ)
            d = domain(Σ)
            C = $Len(d)/2

            if dsp.β==dsp.α==-0.5
                T[k == 1 ? C*π : zero(T) for k=kr]
            else
                T[$Sum(Fun(dsp,[zeros(T,k-1);1])) for k=kr]
            end
        end

        function bandwidths(Σ::$ConcFunc{JacobiWeight{Chebyshev{D,R},D,R,TT}}) where {D<:IntervalOrSegment,R,TT}
            β,α = domainspace(Σ).β,domainspace(Σ).α
            if β==α && isapproxinteger(β-0.5) && 0 ≤ ceil(Int,β)
                0,2ceil(Int,β)
            else
                0,∞
            end
        end
    end
end


text:
"""
```
chull = convexhull(img)
```
Computes the convex hull of a binary image and returns the vertices of convex hull as a CartesianIndex array.
"""
function convexhull(img::AbstractArray{T,2}) where {T<:Union{Bool,Gray{Bool}}}

    function getboundarypoints(img)
        points = CartesianIndex{2}[]
        for j in axes(img, 2)
            v = Base.view(img, :, j)
            i1 = findfirst(v)
            if i1 != nothing
                i2 = findlast(v)
                push!(points, CartesianIndex(i1, j))
                if i1 != i2
                    push!(points, CartesianIndex(i2, j))
                end
            end
        end
        return points
    end

    function right_oriented(ref, a, b)
        return (a[2] - ref[2]) * (b[1] - ref[1]) - (b[2] - ref[2]) * (a[1] - ref[1]) < 0
    end

    function collinear(a, b, c)
        return (a[2] - c[2]) * (b[1] - c[1]) - (b[2] - c[2]) * (a[1] - c[1]) == 0
    end

    dist2(a, b) = sum(abs2, (a - b).I)

    function angularsort!(points, ref, img)
        last_point = ref

        for i in eachindex(points)
            if points[i] == last_point
                deleteat!(points, i)
                break
            end
        end

        sort!(points; lt=(a, b) -> right_oriented(last_point, a, b))

        i = 0
        while i <= length(points)
            i = i + 1
            if i + 1 > length(points)
                break
            end

            if collinear(last_point, points[i], points[i + 1])
                if dist2(last_point, points[i]) < dist2(last_point, points[i + 1])
                    deleteat!(points, i)
                else
                    deleteat!(points, i + 1)
                end
                i = i - 1
            end
        end
    end

    # Used Graham scan algorithm

    points = getboundarypoints(img)
    last_point = CartesianIndex(map(r -> first(r) - 1, axes(img)))
    for point in points
        if point[2] > last_point[2] || (point[2] == last_point[2] && point[1] > last_point[1])
            last_point = point
        end
    end
    angularsort!(points, last_point, img)

    n_points = length(points)

    if n_points < 3
        error("Not enough points to compute convex hull.")
    end

    convex_hull = CartesianIndex{2}[]
    push!(convex_hull, last_point)
    push!(convex_hull, points[1])
    push!(convex_hull, points[2])

    for i in 3:n_points
        while (
            right_oriented(convex_hull[end], convex_hull[end - 1], points[i]) ||
            collinear(convex_hull[end], convex_hull[end - 1], points[i])
        )
            pop!(convex_hull)
        end
        push!(convex_hull, points[i])
    end

    return convex_hull
end


text:
# cross-correlation module
export clean_up!, clean_up, correlate, phasecorrelate
export coherence!, coherence, deconvolution!, deconvolution, whiten, whiten!

"""
    clean_up!(A,freqmin,freqmax,fs)

Demean, detrend, taper and filter time series.

# Arguments
- `A::AbstractArray`: Time series.
- `fs::Real`: Sampling rate of time series `A` in Hz.
- `freqmin::Real`: Pass band low corner frequency in Hz.
- `freqmax::Real`: Pass band high corner frequency in Hz.
"""
function clean_up!(A::AbstractArray, freqmin::Real, freqmax::Real, fs::Real;
                   corners::Int=4,zerophase::Bool=true,max_length::Real=20.)
    detrend!(A)
    taper!(A,fs,max_length=max_length)
    bandpass!(A,freqmin,freqmax,fs,corners=corners,zerophase=zerophase)
    return nothing
end
clean_up(A::AbstractArray, freqmin::Real, freqmax::Real, fs::Real;
         corners::Int=4, zerophase::Bool=true,max_length::Real=20.) =
         (U = deepcopy(A); clean_up!(U,freqmin,freqmax, fs,
         corners=corners, zerophase=zerophase, max_length=max_length); return U)

clean_up!(C::CorrData,freqmin::Real,freqmax::Real; corners::Int=4,
          zerophase::Bool=true,max_length::Real=20.) = (clean_up!(C.corr,
          freqmin,freqmax,C.fs,corners=corners,zerophase=zerophase,
          max_length=max_length);C.freqmin=Float64(freqmin);
          C.freqmax=Float64(freqmax);return nothing)

clean_up(C::CorrData,freqmin::Real,freqmax::Real; corners::Int=4,
         zerophase::Bool=true,max_length::Real=20.) = (U = deepcopy(C);
         clean_up!(U,freqmin,freqmax,corners=corners,zerophase=zerophase,
         max_length=max_length);return U)

clean_up!(R::RawData,freqmin::Real,freqmax::Real; corners::Int=4,
          zerophase::Bool=true,max_length::Real=20.) = (clean_up!(R.x,freqmin,
          freqmax,R.fs,corners=corners,zerophase=true, max_length=max_length);
          R.freqmin=Float64(freqmin);R.freqmax=Float64(freqmax);return nothing)

clean_up(R::RawData,freqmin::Real,freqmax::Real; corners::Int=4,
       zerophase::Bool=true,max_length::Real=20.) = (U = deepcopy(R);
       clean_up!(U,freqmin,freqmax,corners=corners,zerophase=zerophase,
       max_length=max_length);return U)

"""
    correlate(FFT1, FFT2, N, maxlag, corr_type='cross-correlation')

Cross-correlate ambient noise data in the frequency domain.

Cross-correlation can be done using one of three options:

- Cross-correlation: ``C_{AB}(ω) = u_A(ω) u^∗_B(ω)``
- Coherence: ``C_{AB}(ω) = \frac{u_A(ω) u^∗_B(ω)}{∣ u_A(ω) ∣  ∣ u_B(ω) ∣}``
- Deconvolution: ``C_{AB}(ω) = \frac{u_A(ω) u^∗_B(ω)}{∣ u_B(ω) ∣^2}``

Smoothing of FFTs for coherence and deconvolution should be done before
cross-correlating.

# Arguments
- `FFT1::AbstractArray`: Complex Array of fourier transform of ambient noise data.
- `FFT2::AbstractArray`: Complex Array of fourier transform of ambient noise data.
- `N::Int`: Number of input data points in time domain, equal to `cc_len` * `fs`.
- `maxlag::Int`: Number of data points in cross-correlation to save,
                 e.g. `maxlag = 2000` will save lag times = -2000/fs:2000/fs s.
- `corr_type::String`: Type of correlation: `cross-correlation`, `coherence` or
                       `deconv`.
"""
function correlate(FFT1::AbstractArray{Complex{T}}, FFT2::AbstractArray{Complex{T}},
                      N::Int, maxlag::Int) where T <: AbstractFloat
    # take inverse fft
    corrT = irfft(conj.(FFT1) .* FFT2,N,1)

    # return corr[-maxlag:maxlag]
    t = vcat(0:Int(N  / 2)-1, -Int(N  / 2):-1)
    ind = findall(abs.(t) .<= maxlag)
    newind = fftshift(ind,1)
    return corrT[newind,:]
end

"""
    phasecorrelate(FFT1, FFT2, maxlag)

Phase Cross-correlate (PCC) ambient noise data in the frequency domain.

# Arguments
- `FFT1::AbstractArray`: Complex Array of fourier transform of ambient noise data.
- `FFT2::AbstractArray`: Complex Array of fourier transform of ambient noise data.
- `N::Int`: Number of input data points in time domain, equal to `cc_len` * `fs`.
- `maxlag::Int`: Number of data points in cross-correlation to save,
                 e.g. `maxlag = 2000` will save lag times = -2000/fs:2000/fs s.
- `corr_type::String`: Type of correlation: `cross-correlation`, `coherence` or
                       `deconv`.
"""
function phasecorrelate(FFT1::AbstractArray{Complex{T}}, FFT2::AbstractArray{Complex{T}},
                      N::Int, maxlag::Int) where T <: AbstractFloat
    # take inverse fft
    corrT = real.(ifft(conj.(FFT1) .* FFT2,1))

    # return corr[-maxlag:maxlag]
    t = vcat(0:Int(N  / 2)-1, -Int(N  / 2):-1)
    ind = findall(abs.(t) .<= maxlag)
    newind = fftshift(ind,1)
    return corrT[newind,:]
end

"""
    correlate(FFT1, FFT2, maxlag,corr_type="CC")

Cross-correlate ambient noise data in the frequency domain.

Cross-correlation can be done using one of two options:

- CC: Cross-correlation, i.e. ``C_{AB}(ω) = u_A(ω) u^∗_B(ω)``
- PCC: Phase cross-correlation, see [Ventosa et al., 2019]

# Arguments
- `FFT1::FFTData`: FFTData object of fft'd ambient noise data.
- `FFT2::FFTData`: FFTData object of fft'd ambient noise data.
- `maxlag::Real`: Maximum lag time (in seconds) in cross-correlation to save,
                     e.g. `maxlag = 20.` will save lag times = -20.:20. s.
- `corr_type::String`: Type of correlation: `CC` or `PCC`.
"""
function correlate(FFT1::FFTData, FFT2::FFTData, maxlag::Real;corr_type::String="CC")

    comp = FFT1.name[end] * FFT2.name[end]
    # get intersect of dates; return nothing if no intersect
    inter = intersect(FFT1.t,FFT2.t)
    if length(inter) == 0
        throw(ArgumentError("No common windows for $(FFT1.name)-$(FFT2.name) $(FFT1.id)"))
    end

    ind1 = findall(x -> x ∈ inter, FFT1.t)
    ind2 = findall(x -> x ∈ inter, FFT2.t)
    N = convert(Int,round(FFT1.cc_len * FFT1.fs)) # number of data points
    if uppercase(corr_type) == "CC"
        corr = correlate(@views(FFT1.fft[:,ind1]), @views(FFT2.fft[:,ind2]),
                     N,convert(Int,round(maxlag * FFT1.fs)))
    elseif uppercase(corr_type) == "PCC"
        corr = phasecorrelate(@views(FFT1.fft[:,ind1]), @views(FFT2.fft[:,ind2]),
                     N,convert(Int,round(maxlag * FFT1.fs)))
    else
        throw(ArgumentError("Unrecognized cross-correlation type $corr_type. Options are CC and PCC."))
    end

    rotated = false

    return CorrData(FFT1, FFT2, comp, rotated, corr_type, Float64(maxlag), inter, corr)
end

"""
   whiten!(A, freqmin, freqmax, fs, N; pad=50)

Whiten spectrum of rfft `A` between frequencies `freqmin` and `freqmax`.
Returns the whitened rfft of the time series.

# Arguments
- `A::AbstractArray`: Time series.
- `freqmin::Real`: Pass band low corner frequency.
- `freqmax::Real`: Pass band high corner frequency.
- `fs::Real`: Sampling rate of time series `A`.
- `N::Int`: Number of input time domain samples for each rfft.
- `pad::Int`: Number of tapering points outside whitening band.
"""
function whiten!(A::AbstractArray{Complex{T}}, freqmin::Real,
                 freqmax::Real, fs::Real,N::Int;pad::Int=50) where T <: AbstractFloat
   Nrows,Ncols = size(A)

   # get whitening frequencies
   freqvec = FFTW.rfftfreq(N,fs)
   left = findfirst(x -> x >= freqmin, freqvec)
   right = findfirst(freqmax .<= freqvec)
   low, high = left - pad, right + pad

   if low <= 1
       low = 1
       left = low + pad
   end

   if high > length(freqvec)
       high = length(freqvec)- 1
       right = high - pad
   end

   compzero = complex(T(0))
   padarr = similar(A,T,pad)
   padarr .= T(0.):T(pad-1)
   # left zero cut-off
   A[1:low-1,:] .= compzero
   # left tapering
   A[low:left-1,:] .= cos.(T(pi) ./ T(2) .+ T(pi) ./ T(2) .* padarr ./ pad).^2 .* exp.(im .* angle.(A[low:left-1,:]))
   # pass band
   A[left:right-1,:] .= exp.(im .* angle.(A[left:right-1,:]))
   # right tapering
   A[right:high-1,:] .= cos.(T(pi) ./ T(2) .* padarr ./ pad).^2 .* exp.(im .* angle.(A[right:high-1,:]))
   # right zero cut-off
   A[high:end,:] .= compzero
   return nothing
end
whiten(A::AbstractArray, freqmin::Real, freqmax::Real, fs::Real, N::Int;
    pad::Int=50) = (U = deepcopy(A);
    whiten!(U,freqmin,freqmax,fs,N,pad=pad);
    return U)
"""
   whiten(F, freqmin, freqmax)

Whiten spectrum of FFTData `F` between frequencies `freqmin` and `freqmax`.
Uses real fft to speed up computation.
Returns the whitened (single-sided) fft of the time series.

# Arguments
- `F::FFTData`: FFTData object of fft'd ambient noise data.
- `freqmin::Real`: Pass band low corner frequency.
- `freqmax::Real`: Pass band high corner frequency.
- `pad::Int`: Number of tapering points outside whitening band.
"""
function whiten!(F::FFTData, freqmin::Real, freqmax::Real;pad::Int=50)
    if freqmin < F.freqmin && freqmax > F.freqmax
        @warn "Whitening frequencies ($freqmin, $freqmax Hz) are wider than frequencies
        in FFTData ($(F.freqmin),$(F.freqmax) Hz). Whitening in ($(F.freqmin),$(F.freqmax) Hz) band."
    elseif freqmin < F.freqmin
        @warn "Low whitening frequency $freqmin Hz is lower than minumum frequency
        in FFTData ($(F.freqmin) Hz). Whitening in ($(F.freqmin),$freqmax Hz) band."
    elseif freqmax > F.freqmax
        @warn "High whitening frequency $freqmax Hz is higher than maximum frequency
        in FFTData ($(F.freqmax) Hz). Whitening in ($freqmin,$(F.freqmax) Hz) band."
    end

    N = convert(Int, F.fs * F.cc_len) # number of data points
    freqmin = max(freqmin,F.freqmin) # check for freqmin = 0
    freqmax = min(freqmax,max(F.freqmax,1 / F.cc_len)) # check for freqmax = 0
    whiten!(F.fft, freqmin, freqmax, F.fs, N, pad=pad)
    F.whitened = true
    F.freqmin = Float64(freqmin)
    F.freqmax = Float64(freqmax)
    return nothing
end
whiten(F::FFTData, freqmin::Real, freqmax::Real;pad::Int=50) =
      (U = deepcopy(F); whiten!(U,freqmin,freqmax,pad=pad);return U)

function whiten!(R::RawData,freqmin::Real, freqmax::Real; pad::Int=50)
    if freqmin < R.freqmin && freqmax > R.freqmax
        @warn "Whitening frequencies ($freqmin, $freqmax Hz) are wider than frequencies
        in RawData ($(R.freqmin),$(R.freqmax) Hz). Whitening in ($(R.freqmin),$(R.freqmax) Hz) band."
    elseif freqmin < R.freqmin
        @warn "Low whitening frequency $freqmin Hz is lower than minumum frequency
        in RawData ($(R.freqmin) Hz). Whitening in ($(R.freqmin),$freqmax Hz) band."
    elseif freqmax > R.freqmax
        @warn "High whitening frequency $freqmax Hz is higher than maximum frequency
        in FFTData ($(R.freqmax) Hz). Whitening in ($freqmin,$(R.freqmax) Hz) band."
    end

    N = convert(Int, R.fs * R.cc_len) # number of data points
    freqmin = max(freqmin,R.freqmin) # check for freqmin = 0
    freqmax = min(freqmax,max(R.freqmax,1 / R.cc_len)) # check for freqmax = 0
    FFT = rfft(R.x,1)
    whiten!(FFT,freqmin,freqmax,R.fs, N, pad=pad)
    R.x .= irfft(FFT,N,1)
    R.freqmin = Float64(freqmin)
    R.freqmax = Float64(freqmax)
    R.whitened = true
    return nothing
end
whiten(R::RawData,freqmin::Real, freqmax::Real; pad::Int=50) =
      (U = deepcopy(R); whiten!(U,freqmin,freqmax,pad=pad);return U)

function whiten!(N::NodalData,freqmin::Real, freqmax::Real; pad::Int=50)
    @assert freqmin > 0 "Whitening frequency must be greater than zero."
    @assert freqmax <= N.fs[1] / 2 "Whitening frequency must be less than or equal to Nyquist frequency."
    Npts = size(N.data,1) # number of data points
    FFT = rfft(N.data,1)
    whiten!(FFT,freqmin,freqmax,N.fs[1], Npts, pad=pad)
    N.data .= irfft(FFT,Npts,1)
    return nothing
end
whiten(N::NodalData,freqmin::Real, freqmax::Real; pad::Int=50) =
    (U = deepcopy(N); whiten!(U,freqmin,freqmax,pad=pad);return U)

"""

  coherence!(F,half_win, water_level)

Apply coherence method to FFTData `F`. Where,
``C_{AB}(ω) = \frac{u_A(ω) u^∗_B(ω)}{∣ u_A(ω) ∣  ∣ u_B(ω) ∣}``

# Arguments
- `F::FFTData`: FFTData object of fft'd ambient noise data.
- `half_win::Int`: Number of points in half-window to smooth spectrum.
- `water_level::AbstractFloat`: Regularization parameter for spectral smoothing.
                                0.01 is a common value [Mehta, 2007].
"""
function coherence!(F::FFTData, half_win::Int,
                    water_level::Union{Nothing,AbstractFloat}=nothing)
    smoothF = smooth(abs.(F.fft),half_win)
    if !isnothing(water_level)
        reg = water_level .* mean(abs.(F.fft),dims=1)
        smoothF .+= reg
    end
    F.fft ./= smoothF
    return nothing
end
coherence(F::FFTData,half_win::Int,
          water_level::Union{Nothing,AbstractFloat}=nothing) =
          (U = deepcopy(F);coherence!(U,half_win,water_level);
          return U)

"""

  deconvolution!(F,half_win, water_level)

Apply deconvolution method to FFTData `F`. Where,
``C_{AB}(ω) = \frac{u_A(ω) u^∗_B(ω)}{∣ u_B(ω) ∣^2}``

# Arguments
- `F::FFTData`: FFTData object of fft'd ambient noise data.
- `half_win::Int`: Number of points in half-window to smooth spectrum.
- `water_level::AbstractFloat`: Regularization parameter for spectral smoothing.
                              0.01 is a common value [Mehta, 2007].
"""
function deconvolution!(F::FFTData, half_win::Int,
                        water_level::Union{Nothing,AbstractFloat}=nothing)
    smoothF = smooth(abs.(F.fft).^2,half_win)
    if !isnothing(water_level)
        reg = water_level .* mean(abs.(F.fft).^2,dims=1)
        smoothF .+= reg
    end
    F.fft ./= smoothF
    return nothing
end
deconvolution(F::FFTData,half_win::Int,
              water_level::Union{Nothing,AbstractFloat}=nothing) =
              (U = deepcopy(F);deconvolution!(U,half_win,water_level);
              return U)


text:
using Test

function naive!(C, A, B)
    row, col = size(C)
    colA = size(A)[2]
    for i in 1:row
        for j in 1:col
            for k in 1:colA
                C[i, j] = C[i, j] + A[i, k] * B[k, j]
            end
        end
    end
end

function unroll_k!(C, A, B)
    row, col = size(C)
    colA = size(A)[2]
    for i in 1:row
        for j in 1:col
            for k in 1:2:colA
                C[i, j] = C[i, j] + A[i, k] * B[k, j] + A[i, k+1] * B[k+1, j]
            end
        end
    end
end

function unroll_j!(C, A, B)
    row, col = size(C)
    colA = size(A)[2]
    for i in 1:row
        for j in 1:2:col
            for k in 1:colA
                C[i,   j] = C[i, j] + A[i, k] * B[k, j]
                C[i, j+1] = C[i, j+1] + A[i, k] * B[k, j+1]
            end
        end
    end
end

function unroll_i!(C, A, B)
    row, col = size(C)
    colA = size(A)[2]
    for i in 1:2:row
        for j in 1:col
            for k in 1:colA
                C[i,   j] = C[i,   j] + A[i,   k] * B[k, j]
                C[i+1, j] = C[i+1, j] + A[i+1, k] * B[k, j]
            end
        end
    end
end

function unroll_ij!(C, A, B)
    row, col = size(C)
    colA = size(A)[2]
    for i in 1:2:row
        for j in 1:2:col
            for k in 1:colA
                C[i ,  j] = C[i,   j] + A[i, k] * B[k, j]
                C[i+1, j] = C[i+1, j] + A[i+1, k] * B[k ,j]
                C[i, j+1] = C[i, j+1] + A[i, k] * B[k, j+1]
                C[i+1, j+1] = C[i+1, j+1] + A[i+1, k] * B[k, j+1]
            end
        end
    end
end

function main(N)
    A = reshape(convert.(Float64, collect(1:N^2)), N, N)
    B = reshape(convert.(Float64, collect(1:N^2)), N, N)
    C = zeros(N, N)
    @time refC = A * B
    C = zeros(N, N)
    @time naive!(C, A, B)
    C = zeros(N ,N)
    @time unroll_k!(C, A, B)
    @test C ≈ refC
    C = zeros(N, N)
    @time unroll_i!(C, A, B)
    @test C ≈ refC
    C = zeros(N, N)
    @time unroll_j!(C, A, B)
    @test C ≈ refC
    C = zeros(N, N)
    @time unroll_ij!(C, A, B)
    @test C ≈ refC
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# AUTO GENERATED FILE - DO NOT EDIT

export datetimepicker

"""
    datetimepicker(;kwargs...)

A DateTimePicker component.
DateTimePicker is a datetime input component.
The inputs can be initialized with the `defaultValue` property and the
input date, on ISO format, is specified with the `value` property.
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `value` (String; optional): The input's current value.
- `defaultValue` (String; optional): The input's default value. Accepts values on the ISO format.
- `placeholder` (String; optional): The date input placeholder.
- `timezonePlaceholder` (String; optional): The timezone input placeholder.
- `renderTimezone` (Bool; optional): Defines if the timezone input should be rendered. Defaults to true.
- `style` (Dict | String; optional): The container's style.
- `dateInputStyle` (Dict | String; optional): The date input's style.
- `timezoneInputStyle` (Dict; optional): The timezone input's container style.
"""
function datetimepicker(; kwargs...)
        available_props = Symbol[:id, :value, :defaultValue, :placeholder, :timezonePlaceholder, :renderTimezone, :style, :dateInputStyle, :timezoneInputStyle]
        wild_props = Symbol[]
        return Component("datetimepicker", "DateTimePicker", "dash_cool_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: raulrozza/dash-cool-components

max_issues_repo_name: raulrozza/dash-cool-components

max_forks_repo_name: raulrozza/dash-cool-components


text:
@ghdef mutable struct Tag
    tag::Union{String, Nothing}
    sha::Union{String, Nothing}
    url::Union{HTTP.URI, Nothing}
    message::Union{String, Nothing}
    tagger::Union{Dict, Nothing}
    object::Union{Dict, Nothing}
    verification::Union{Dict, Nothing}
end

namefield(tag::Tag) = tag.sha

@api_default function tag(api::GitHubAPI, repo, tag_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/git/refs/tags/$(name(tag_obj))"; options...)
    return Tag(result)
end

@api_default function tags(api::GitHubAPI, repo; options...)
    result, paged_data = gh_get_paged_json(api, "/repos/$(name(repo))/git/refs/tags"; options...)
    return map(Tag, result), paged_data
end

@api_default function create_tag(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/git/tags"; options...)
    return Tag(result)
end


text:
module FPTControl

##################################
# Fixed Point Transformation Control toolbox #
##################################

# In PID case n=the order of the system +1
# In PD case n=the order of the system
# lamda is a control parameter for tuning the PID or PD
# errors is a vector where you put togather the related errors.
# example for PID n=3 (2nd order system) errors=[error_int,error,error_p], where error=qNominal - q (by defininition nominal Trajectory - Realized Trajectory)
# errror_p is the time derivative of the error
# for PD ,when n=2 (2nd order system) errors=[error,error_p]
function KB(n,lambda,errors,nominal)
  out=0
  for l=0:n-1
    out=out+(factorial(n)/(factorial(l)*factorial(n-l))*errors[l+1]*lambda^(n-l))
  end
  out=out+nominal
  return out
end
# Simple Euler integral (rectangle method)
function Integ(Integral,step,input)
    out=Integral+step*input
    return out
end

# Simple Sigmoid function, mostly you can use any Sigmoid type function, like tanh()
function sigmoid(x)
	s=x/(1+abs(x))
  return s
end
# The Original Deformation function for SISO case (note here I used tanh() for sigmoid function)
function G_SISO(past_input,past_response,desired,Kc,Bc,Ac)
  out=(Kc+past_input)*(1 + Bc * tanh(Ac * (past_response - desired))) - Kc
  return out
end

# The Original Deformation function for MIMO case
function G_MIMO(past_input,past_response,desired,error_limit,Kc,Bc,Ac) # Inputs and Outputs are vectors
#  need control parameters  K B A
	error_norm=vecnorm(past_response - desired,2)
	# If the norm of the error is greater then the limnit compute the deformation
	# (it is not near the Fixed Point)
  if error_norm>error_limit
    e_direction=(past_response-desired)/error_norm
    B_factor= Bc* sigmoid(Ac *error_norm)
    out=(1+B_factor)*past_input +B_factor*Kc*e_direction
  else
    out=past_input # Almost in the Fixed Point
  end
  return out
end

export KB,Integ,G_SISO,G_MIMO
end # module


text:
module CallbackTests

using Photon, Test
using Photon.Layers
using Photon.Callbacks
using Photon.Metrics
using Photon.Losses

function simple_conv_model()
    model = Sequential(
        Conv2D(4, 3, :relu),
        Conv2D(8, 3, :relu),
        MaxPool2D(),
        Dense(16, :relu),
        Dense(10)
    )
    return model
end

function getdata(s=28)
    [(
        KorA(randn(Float32,s,s,1,16)),
        KorA(randn(Float32,10,16))
    ) for i=1:10]
end

function test_callbacks()
    model = simple_conv_model()
    workout = Workout(model, MSELoss())

    data = getdata()
    filename = "myworkout.dat"
    rm(filename, force=true)

    train!(workout, data, epochs=2, cb=AutoSave(:loss, filename))
    @assert isfile(filename)
    rm(filename, force=true)

    train!(workout, data, epochs=2, cb=EpochSave(filename))
    @assert isfile(filename)
    rm(filename, force=true)

    val_data = getdata()
    workout = Workout(model, MSELoss())
    train!(workout, data, val_data, epochs=10, cb=EarlyStop(:val_loss))
    @assert workout.epochs < 10
end

@testset "Callbacks" begin
    test_callbacks()
end

end


text:
@testset "Tracking results" begin
    results = Tracking.TrackingResults(
        TrackingState(GPSL1, 100Hz, 100),
        EarlyPromptLateCorrelator(NumAnts(2)),
        true,
        Tracking.BitBuffer(),
        45dBHz
    )

    @test @inferred(get_carrier_doppler(results)) == 100Hz
    @test @inferred(get_carrier_phase(results)) == 0.0
    @test @inferred(get_code_doppler(results)) == 100Hz / 1540
    @test @inferred(get_code_phase(results)) == 100
    correlator = @inferred get_correlator(results)
    @test correlator == EarlyPromptLateCorrelator(NumAnts(2))
    @test @inferred(get_early(results)) == [0.0, 0.0]
    @test @inferred(get_prompt(results)) == [0.0, 0.0]
    @test @inferred(get_late(results)) == [0.0, 0.0]
    @test @inferred(get_bits(results)) == 0
    @test @inferred(get_num_bits(results)) == 0
    @test @inferred(get_secondary_code_or_bit_found(results)) == false
    @test @inferred(get_cn0(results)) == 45dBHz

    results = Tracking.TrackingResults(
        TrackingState(GPSL1, 100Hz, 100),
        EarlyPromptLateCorrelator(NumAnts(2)),
        false,
        Tracking.BitBuffer(),
        45dBHz
    )
    correlator = @inferred get_correlator(results)
    @test get_early(correlator) == [0, 0]
    @test get_prompt(correlator) == [0, 0]
    @test get_late(correlator) == [0, 0]
end


text:
module StandardLibrary

using ModelingToolkit
using OrderedCollections

export standard_library

@parameters t
D = Differential(t)

@variables E[1:2](t) F[1:2](t)

# Linear resistance (:R)
@parameters R
R_dict = Dict(
    :description => "Generalised Linear Resistor",
    :numports => 1,
    :parameters => OrderedDict(
        R => "Resistance"
    ),
    :equations => [0 ~ E[1] - R*F[1]],
    :defaults => Dict(R => 1.)
)

# Linear capacitor (:C)
@parameters C
@variables q(t)
C_dict = Dict(
    :description => "Generalised Linear Resistor",
    :numports => 1,
    :parameters => OrderedDict(
        C => "Capacitance"
    ),
    :states => OrderedDict(
        q => "Generalised Position"
    ),
    :equations => [
        0 ~ q/C - E[1],
        D(q) ~ F[1]
    ],
    :defaults => Dict(C => 1., q => 0.)
)

# Linear inductance (:I)
@parameters L
@variables p(t)
I_dict = Dict(
    :description => "Generalised Linear Inductor",
    :numports => 1,
    :parameters => OrderedDict(
        L => "Inductance"
    ),
    :states => OrderedDict(
        p => "Generalised Momentum"
    ),
    :equations => [
        0 ~ p/L - F[1],
        D(p) ~ E[1]
    ],
    :defaults => Dict(L => 1., p => 0.)
)

# Source of effort (:Se)
@parameters e
Se_dict = Dict(
    :description => "Effort source",
    :numports => 1,
    :parameters => OrderedDict(
        e => "Effort"
    ),
    :equations => [0 ~ e - E[1]],
    :defaults => Dict(e => 1.)
)

# Source of flow (:Sf)
@parameters f
Sf_dict = Dict(
    :description => "Flow source",
    :numports => 1,
    :parameters => OrderedDict(
        f => "Flow"
    ),
    :equations => [0 ~ f - F[1]],
    :defaults => Dict(f => 1.)
)

# Transformer (:TF)
@parameters n
TF_dict = Dict(
    :description => "Linear Transformer",
    :numports => 2,
    :parameters => OrderedDict(
        n => "Ratio"
    ),
    :equations => [
        0 ~ E[2] - n * E[1],
        0 ~ F[1] - n * F[2]
    ],
    :defaults => Dict(n => 1.)
)

const standard_library = Dict(
    :R => R_dict,
    :C => C_dict,
    :I => I_dict,
    :Se => Se_dict,
    :Sf => Sf_dict,
    :TF => TF_dict
)

end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
@testset "JuMP: constraints" begin
    m = Model(CBLS.Optimizer)

    err = _ -> 1.0
    concept = _ -> true

    @variable(m, X[1:10], DiscreteSet(1:4))

    @constraint(m, X in Error(err))
    @constraint(m, X in Predicate(concept))

    @constraint(m, X in AllDifferent())
    @constraint(m, X in AllEqual())
    @constraint(m, X in AllEqualParam(2))
    @constraint(m, X in AlwaysTrue())
    @constraint(m, X[1:4] in DistDifferent())
    @constraint(m, X[1:2] in Eq())
    @constraint(m, X in Ordered())
end

@testset "JuMP: sudoku 9x9" begin
    m, X = sudoku(3)
    optimize!(m)
    solution_ = value.(X)
    @info solution_summary(m)
    display(solution_, Val(:sudoku))
end

@testset "JuMP: golomb(5)" begin
    m, X = golomb(5)
    optimize!(m)
    @info solution_summary(m)
    @info "JuMP: golomb(5)" value.(X)
end

@testset "JuMP: magic_square(3)" begin
    m, X = magic_square(3)
    optimize!(m)
    @info solution_summary(m)
    @info "JuMP: magic_square(3)" value.(X)
end

@testset "JuMP: n_queens(5)" begin
    m, X = n_queens(5)
    optimize!(m)
    @info solution_summary(m)
    @info "JuMP: n_queens(5)" value.(X)
end

@testset "JuMP: qap(12)" begin
    m, X = qap(12, qap_weights, qap_distances)
    optimize!(m)
    @info solution_summary(m)
    @info "JuMP: qap(12)" value.(X)
end

@testset "JuMP: basic opt" begin
    model = Model(CBLS.Optimizer)
    set_optimizer_attribute(model, "iteration", 100)
    @test get_optimizer_attribute(model, "iteration") == 100
    set_time_limit_sec(model, 5.0)
    @test time_limit_sec(model) == 5.0

    @variable(model, x in DiscreteSet(0:20))
    @variable(model, y in DiscreteSet(0:20))

    @constraint(model, [x,y] in Predicate(v -> 6v[1] + 8v[2] >= 100 ))
    @constraint(model, [x,y] in Predicate(v -> 7v[1] + 12v[2] >= 120 ))

    objFunc = v -> 12v[1] + 20v[2]
    @objective(model, Min, ScalarFunction(objFunc))

    optimize!(model)
    @info solution_summary(model)
    @info "JuMP: basic opt" value(x) value(y) (12*value(x)+20*value(y))
end

@testset "JuMP: Chemical equilibrium" begin
    m, X = chemical_equilibrium(atoms_compounds, elements_weights, standard_free_energy)
    # set_optimizer_attribute(m, "iteration", 10000)
    # set_time_limit_sec(m, 120.0)
    optimize!(m)
    @info solution_summary(m)
    @info "JuMP: $compounds_names ⟺ $mixture_name" value.(X)
end

@testset "JuMP: Scheduling" begin
    model, completion_times, start_times = scheduling(processing_times, due_times)
    optimize!(model)
    @info solution_summary(model)
    @info "JuMP: scheduling" value.(start_times) value.(completion_times) processing_times due_times
    @info (value.(start_times)+processing_times == value.(completion_times))
end

@testset "JuMP: Minimum Cut" begin
    graph = zeros(5, 5)
    graph[1,2] = 1.0
    graph[1,3] = 2.0
    graph[1,4] = 3.0
    graph[2,5] = 1.0
    graph[3,5] = 2.0
    graph[4,5] = 3.0

    model, links = mincut(graph; source=1, sink=5, interdiction=1)
    optimize!(model)
    @info solution_summary(model)
    @info "JuMP: min-1-cut" value.(links)
end


text:
#=
#cd("/home/david/Documents/gitclones/AlphaZero.jl")
cd("/home/davidyang/Documents/gitclones/AlphaZero.jl")
using Pkg
Pkg.activate(".")
using AlphaZero
#include("/home/david/Documents/gitclones/AlphaZero.jl/games/amazon/main.jl")
include("/home/davidyang/Documents/gitclones/AlphaZero.jl/games/amazon/main.jl")

gs = Amazon.GameSpec()
g = GI.init(gs)
GI.render(g)
while !g.finished
  action = g|>GI.actions_mask|>findall|>rand
  println("action taken ",GI.action_string(gs,action))
  GI.play!(g,action)
  GI.render(g)
end

gs = Amazon.GameSpec()
g = GI.init(gs)
GI.render(g)
Int.(GI.vectorize_state(gs,GI.current_state(g))[:,:,4])
GI.play!(g,CartesianIndex(1,2))
GI.render(g)
Int.(GI.vectorize_state(gs,GI.current_state(g))[:,:,5])
GI.play!(g,CartesianIndex(4,5))
GI.render(g)
Int.(GI.vectorize_state(gs,GI.current_state(g))[:,:,6])
=#

cd("/home/davidyang/Documents/gitclones/AlphaZero.jl")
#cd("/home/david/Documents/gitclones/AlphaZero.jl")
#@everywhere using Pkg
#@everywhere Pkg.activate(".")
#@everywhere using AlphaZero
using Pkg
Pkg.activate(".")
using AlphaZero
experiment = Examples.experiments["amazon"]
session = Session(experiment, dir="sessions/amazon2")
resume!(session)

duel1 = session.benchmark[1]
#duel1 equivalent to eval
progress = UI.Log.Progress(session.logger, duel1.sim.num_games)
env = session.env
net() = Network.copy(env.bestnn, on_gpu=duel1.sim.use_gpu, test_mode=true)
simulator = Simulator(net, record_trace) do net
  player = Benchmark.instantiate(duel1.player, env.gspec, net)
  baseline = Benchmark.instantiate(duel1.baseline, env.gspec, net)
  return TwoPlayers(player, baseline)
end
#=
samples, elapsed = @timed simulate(
    simulator, env.gspec, duel1.sim,
    game_simulated=(() -> next!(progress)))
=#
oracles = simulator.make_oracles()
spawn_oracles, done =
    AlphaZero.batchify_oracles(oracles; duel1.sim.num_workers, duel1.sim.batch_size, duel1.sim.fill_batches)

println("debugging play game")
trace = @time play_game(env.gspec, simulator.make_player(oracles), flip_probability=0.0)

flip_probability = 0.0
player = simulator.make_player(oracles)
game = GI.init(env.gspec)
trace = Trace(GI.current_state(game))
#actions, π_target = think(player, game)
env = player.white.mcts
#MCTS.explore!(env,game,player.white.niters)
#η = MCTS.dirichlet_noise(game, env.noise_α)
#MCTS.run_simulation!(env,GI.clone(game),η=η)
#MCTS.state_info(env,GI.current_state(game))
#env.oracle(state)
#Network.evaluate(env.oracle,state)
nn = env.oracle
gspec = Network.game_spec(nn)
actions_mask = GI.actions_mask(GI.init(gspec, state))
x = GI.vectorize_state(gspec, state)
a = Float32.(actions_mask)
xnet, anet = Network.to_singletons.(Network.convert_input_tuple(nn, (x, a)))
#net_output = Network.forward_normalized(nn, xnet, anet)
p, v = Network.forward(nn, xnet)

while true
  if GI.game_terminated(game)
    return trace
  end
  if !iszero(flip_probability) && rand() < flip_probability
    GI.apply_random_symmetry!(game)
  end
  actions, π_target = think(player, game)
  τ = player_temperature(player, game, length(trace))
  π_sample = apply_temperature(π_target, τ)
  a = actions[Util.rand_categorical(π_sample)]
  GI.play!(game, a)
  push!(trace, π_target, GI.white_reward(game), GI.current_state(game))
end


text:
abstract type AbstractMINERVA end

"""
    MINERVA(;L=.8, n_features=10, n_traces=20, c=[0.0])

Model object for the MINERVA model model.

Parameters
- `L`: learning rate or encoding accuracy probability
* `memory`: a matrix containing memory traces where each row corresponds to a trace
* `n_features`: the number of features per trace
* `n_traces`: the number of traces in model
* `c`: a vector of decision criteria

**References**

Hintzman, D. L. (1988). Judgments of frequency and recognition model in a multiple-trace model model.
    Psychological Review, 95(4), 528.
"""
mutable struct MINERVA <: AbstractMINERVA
    L::Float64
    memory::Array{Int,2}
    n_features::Int
    n_traces::Int
    c::Vector{Float64}
end

Broadcast.broadcastable(x::MINERVA) = Ref(x)

function MINERVA(;L=.8, n_features=10, n_traces=20, c=[0.0])
    memory = fill(0, n_traces, n_features)
    return MINERVA(L, memory, n_features, n_traces, c)
end

function stimulus_similarity(model, p)
     rand() <= model.h ? (return p) : (return rand(setdiff(-1:1, p)))
end

"""
    random_stimulus(model::AbstractMINERVA)

Generates random stimulus with a specified number of features in the `model` object.
"""
random_stimulus(model::AbstractMINERVA) = random_stimulus(model.n_features)
random_stimulus(n_features) = rand(-1:1, n_features)

"""
    random_stimulus(model::AbstractMINERVA, n)

Generates `n` random stimuli with a specified number of features in model object. Returns a `n` by
`n_features` matrix.
"""
random_stimulus(model::AbstractMINERVA, n) = random_stimulus(n, model.n_features)
random_stimulus(n_traces, n_features) = rand(-1:1, n_traces, n_features)

random_distractor(model::AbstractMINERVA) = random_stimulus(model::AbstractMINERVA)

function encode(model::AbstractMINERVA, x)
    x == 0 ? (return 0) : nothing
    rand() <= model.L ? (return x) : (return 0)
end

"""
    encode!(model, stimuli)

Encodes a matrix of stimuli with accuracy governed by parameter `L`

* `model`: `MINERVA` model object
* `stimuli`: a matrix in which rows represent stimuli and columns represent features

"""
function encode!(model::AbstractMINERVA, stimuli)
    model.memory .= encode.(model, stimuli)
    return nothing
end

"""
    compute_intensity(model::AbstractMINERVA, probe)

Computes echo intensity using a given memory `probe`.

* `model`: `MINERVA` model object
* `probe`: a vector of feature values representing a memory probe

"""
compute_intensity(model::AbstractMINERVA, probe) = compute_intensity(model.memory, probe)

function compute_intensity(memory, probe)
    return sum(compute_activation(memory, probe))
end

"""
    compute_activation(M, p)

Computes activation value for each trace in matrix `M` for a given probe `p`.

* `M`: a matrix in which rows represent memory traces and columns represent feature values
* `p`: a vector of feature values representing a memory probe

"""
compute_activation(M, p) = compute_similarity(M, p) .^ 3

"""
    compute_activation(M, p)

Computes simularity between memory traces in matrix `M` and memory probe `p` using the dot product.

* `M`: a matrix in which rows represent memory traces and columns represent feature values
* `p`: a vector of feature values representing a memory probe
"""
function compute_similarity(M, p)
    s = M * p
    n = length(p) .- sum((M .== 0) * (p .== 0), dims=2)
    return s ./ n
end


text:
module jLRS

# standard library
using Printf

# third party libraries
using Distributions
using HDF5
using Proj4
using Glob
using AstroLib
using DataFrames
using Dates
using LinearAlgebra
using Random
using SatelliteToolbox
using SQLite
using SparseArrays
using StaticArrays

jLRS_seed = 123456
# Seed the generator at the beginning. So the *same* code with the *same* inputs
# will always produce the same numbers. If even one parameter changes, we will
# get different results!
Random.seed!(jLRS_seed)

# Type definitions required by the rest of the module
include("types.jl")

# Common math-y functions
include("misc.jl")

# Solar data required for many calculations
include("solar.jl")

# Surface samplings (BRDF, land cover, etc.)
include("surface.jl")

# SIF-related functions
include("SIF.jl")

# Functions for OCO-type
include("OCO.jl")

# Functions for Geostationary-type
include("Geostationary.jl")

# Functions for SwathNadir-type
include("SwathNadirOrbiter.jl")

# Main calculations
include("light_response.jl")

# Plotting tools
include("plotting.jl")

@info "You have loaded jLRS - thanks!"

end # module


text:
using Test

using Downhill
using Downhill: argumentvec, gradientvec, step_origin, fnval, fnval_origin
using Downhill: mcholesky!

using LinearAlgebra

using Random

const test_rng = MersenneTwister(123581321)

"""
    rosenbrock!(x, g; b=2)

Return the value of Rosenbrock function in `length(x)` dimensions.
"""
function rosenbrock!(x::AbstractVector, g::AbstractVector; b=2)
    f = zero(eltype(g))
    fill!(g, 0)
    inds = eachindex(x, g)
    for i in 2:last(inds)
        f += (1 - x[i-1])^2 + b * (x[i] - x[i-1]^2)^2
        g[i-1] += 2 * (x[i-1] - 1) + 4 * b * x[i-1] * (x[i-1]^2 - x[i])
        g[i] += 2 * b * (x[i] - x[i-1]^2)
    end
    return f, g
end

OPT_TYPES = (
    SteepestDescent,
    HyperGradDescent,
    FixedRateDescent,
    MomentumDescent,
    NesterovMomentum,
    CGDescent,
    BFGS,
    CholBFGS,
)

@testset "Testing package DescentMethods" begin
# Test that all methods have a constructor
include("constructors.jl")

# Test that all methods have interface functions
include("interface.jl")

# Test linesearch on quadratic functions
include("linesearch.jl")

# Test modified Cholesky decomposition
include("cholesky.jl")

# Test convergence on the Rosenbrock function
include("convergence.jl")
end