text:
# This is all AD-related stuff. If you're looking to understand TemporalGPs, this can be
# safely ignored.

using Zygote: @adjoint, accum, AContext

# This context doesn't allow any globals.
struct NoContext <: Zygote.AContext end

# Stupid implementation to obtain type-stability.
Zygote.cache(cx::NoContext) = (cache_fields=nothing)

# Stupid implementation.
Base.haskey(cx::NoContext, x) = false

Zygote.accum_param(::NoContext, x, Δ) = Δ

nograd_pullback(Δ) = nothing

Zygote._pullback(::AContext, ::typeof(eltype), x) = eltype(x), nograd_pullback

# Hacks to help the compiler out in very specific situations.
Zygote.accum(a::Array{T}, b::Array{T}) where {T<:Real} = a + b

Zygote.accum(a::SArray{size, T}, b::SArray{size, T}) where {size, T<:Real} = a + b

Zygote.accum(a::Tuple, b::Tuple, c::Tuple) = map(Zygote.accum, a, b, c)

function Zygote._pullback(
    ::AContext, ::Type{SArray{S, T, N, L}}, x::NTuple{L, T},
) where {S, T, N, L}
    SArray_pullback(Δ::Nothing) = nothing
    SArray_pullback(Δ::NamedTuple{(:data,)}) = nothing, Δ.data
    SArray_pullback(Δ::SArray{S}) = nothing, Δ.data
    return SArray{S, T, N, L}(x), SArray_pullback
end

function Zygote._pullback(
    ctx::AContext, ::Type{SArray{S, T, N, L}}, x::NTuple{L, Any},
) where {S, T, N, L}
    new_x, convert_pb = Zygote._pullback(ctx, StaticArrays.convert_ntuple, T, x)
    out, pb = Zygote._pullback(ctx, SArray{S, T, N, L}, new_x)
    SArray_pullback(Δ::Nothing) = nothing
    SArray_pullback(Δ::SArray{S}) = SArray_pullback((data=Δ.data,))
    function SArray_pullback(Δ::NamedTuple{(:data,)})
        _, Δnew_x = pb(Δ)
        _, ΔT, Δx = convert_pb(Δnew_x)
        return nothing, ΔT, Δx
    end
    return SArray{S, T, N, L}(x), SArray_pullback
end

Zygote.@adjoint function collect(x::SArray{S, T, N, L}) where {S, T, N, L}
    collect_pullback(Δ::Array) = ((data = ntuple(i -> Δ[i], Val(L)), ), )
    return collect(x), collect_pullback
end

Zygote.@adjoint function vcat(A::SVector{DA}, B::SVector{DB}) where {DA, DB}
    function vcat_pullback(Δ::SVector)
        ΔA = Δ[SVector{DA}(1:DA)]
        ΔB = Δ[SVector{DB}((DA+1):(DA+DB))]
        return ΔA, ΔB
    end
    return vcat(A, B), vcat_pullback
end

# Implementation of the matrix exponential that assumes one doesn't require access to the
# gradient w.r.t. `A`, only `t`. The former is a bit compute-intensive to get at, while the
# latter is very cheap.

time_exp(A, t) = exp(A * t)
ZygoteRules.@adjoint function time_exp(A, t)
    B = exp(A * t)
    return B, Δ->(nothing, sum(Δ .*  (A * B)))
end

# THIS IS A TEMPORARY FIX WHILE I WAIT FOR #445 IN ZYGOTE TO BE MERGED.
# FOR SOME REASON THIS REALLY HELPS...
@adjoint function (::Type{T})(x, sz) where {T <: Fill}
    back(Δ::AbstractArray) = (sum(Δ), nothing)
    back(Δ::NamedTuple) = (Δ.value, nothing)
    return Fill(x, sz), back
end

function Zygote._pullback(::Zygote.AContext, ::typeof(vcat), x::Zeros, y::Zeros)
    vcat_pullback(Δ) = (nothing, nothing, nothing)
    return vcat(x, y), vcat_pullback
end

@adjoint function collect(x::Fill)
    function collect_Fill_back(Δ)
        return ((value=reduce(accum, Δ), axes=nothing),)
    end
    return collect(x), collect_Fill_back
end

@adjoint function step(x::StepRangeLen)
    return step(x), Δ -> ((ref=nothing, step=Δ, len=nothing, offset=nothing),)
end

@adjoint function BlockDiagonal(blocks::Vector)
    function BlockDiagonal_pullback(Δ::NamedTuple{(:blocks,)})
        return (Δ.blocks,)
    end
    return BlockDiagonal(blocks), BlockDiagonal_pullback
end

@adjoint function Base.map(f::Tf, x::Fill) where {Tf}
    y_el, back = Zygote._pullback(__context__, f, x.value)
    function map_Fill_pullback(Δ::Union{NamedTuple, Tangent})
        if Δ isa Tangent
            Δ_ = (value=Δ.value, axes=Δ.axes)
        else
            Δ_ = Δ
        end
        Δf, Δx_el = back(Δ_.value)
        return Δf, (value = Δx_el, axes=nothing)
    end
    return Fill(y_el, size(x)), map_Fill_pullback
end

function Base.map(f::Tf, x1::Fill, x2::Fill) where {Tf}
    @assert size(x1) == size(x2)
    y_el = f(x1.value, x2.value)
    return Fill(y_el, size(x1))
end

function Base.map(f::Tf, x1::Fill, x2::Fill) where {Tf<:Function}
    @assert size(x1) == size(x2)
    y_el = f(x1.value, x2.value)
    return Fill(y_el, size(x1))
end

Zygote.@adjoint function Base.map(f::Tf, x1::Fill, x2::Fill) where {Tf}
    @assert size(x1) == size(x2)
    y_el, back = Zygote._pullback(__context__, f, x1.value, x2.value)
    function map_Fill_pullback(Δ::NamedTuple)
        Δf, Δx1_el, Δx2_el = back(Δ.value)
        return (Δf, (value = Δx1_el, axes=nothing), (value = Δx2_el, axes=nothing))
    end
    return Fill(y_el, size(x1)), map_Fill_pullback
end

@adjoint function Base.getindex(x::Fill, n::Int)
    function getindex_FillArray_pullback(Δ)
        return ((value = Δ, axes = nothing), nothing)
    end
    return x[n], getindex_FillArray_pullback
end

@adjoint function Base.getindex(x::SVector{1}, n::Int)
    getindex_SArray_pullback(Δ) = (SVector{1}(Δ), nothing)
    return x[n], getindex_SArray_pullback
end

@adjoint function Base.getindex(x::SVector{1, 1}, n::Int)
    getindex_pullback(Δ) = (SMatrix{1, 1}(Δ), nothing)
    return x[n], getindex_SArray_pullback
end

#
# AD-free pullbacks for a few things. These are primitives that will be used to write the
# gradients.
#

function cholesky_rrule(Σ::Symmetric{<:Real, <:StridedMatrix})
    C = cholesky(Σ)
    function cholesky_pullback(Δ::NamedTuple)
        U, Ū = C.U, Δ.factors
        Σ̄ = Ū * U'
        Σ̄ = LinearAlgebra.copytri!(Σ̄, 'U')
        Σ̄ = ldiv!(U, Σ̄)
        BLAS.trsm!('R', 'U', 'T', 'N', one(eltype(Σ)), U.data, Σ̄)

        for n in diagind(Σ̄)
            Σ̄[n] /= 2
        end
        return (UpperTriangular(Σ̄),)
    end
    return C, cholesky_pullback
end

function cholesky_rrule(S::Symmetric{<:Real, <:StaticMatrix{N, N}}) where {N}
    C = cholesky(S)
    function cholesky_pullback(Δ::NamedTuple)
        U, ΔU = C.U, Δ.factors
        ΔS = U \ (U \ SMatrix{N, N}(Symmetric(ΔU * U')))'
        ΔS = ΔS - Diagonal(ΔS ./ 2)
        return ((data=SMatrix{N, N}(UpperTriangular(ΔS)), ),)
    end
    return C, cholesky_pullback
end

@adjoint function cholesky(S::Symmetric{<:Real, <:StaticMatrix{N, N}}) where {N}
    return cholesky_rrule(S)
end

function logdet_pullback(C::Cholesky)
    return logdet(C), function(Δ)
        return ((uplo=nothing, info=nothing, factors=Diagonal(2 .* Δ ./ diag(C.factors))),)
    end
end

function Zygote.accum(a::UpperTriangular, b::UpperTriangular)
    return UpperTriangular(Zygote.accum(a.data, b.data))
end

function Zygote.accum(D::Diagonal{<:Real}, U::UpperTriangular{<:Real, <:SMatrix})
    return UpperTriangular(D + U.data)
end

function Zygote.accum(a::Diagonal, b::UpperTriangular)
    return UpperTriangular(a + b.data)
end

Zygote.accum(a::UpperTriangular, b::Diagonal) = Zygote.accum(b, a)

Zygote._symmetric_back(Δ::UpperTriangular{<:Any, <:SArray}, uplo) = Δ
function Zygote._symmetric_back(Δ::SMatrix{N, N}, uplo) where {N}
    if uplo === 'U'
        return SMatrix{N, N}(UpperTriangular(Δ) + UpperTriangular(Δ') - Diagonal(Δ))
    else
        return SMatrix{N, N}(LowerTriangular(Δ) + LowerTriangular(Δ') - Diagonal(Δ))
    end
end

# Temporary hacks.

using Zygote: literal_getproperty, literal_indexed_iterate, literal_getindex

function Zygote._pullback(::NoContext, ::typeof(*), A::Adjoint, B::AbstractMatrix)
    times_pullback(::Nothing) = nothing
    times_pullback(Δ) = nothing, Adjoint(B * Δ'), A' * Δ
    return A * B, times_pullback
end

function Zygote._pullback(::NoContext, ::typeof(literal_getproperty), C::Cholesky, ::Val{:U})
    function literal_getproperty_pullback(Δ)
        return (nothing, (uplo=nothing, info=nothing, factors=UpperTriangular(Δ)))
    end
    literal_getproperty_pullback(Δ::Nothing) = nothing
    return literal_getproperty(C, Val(:U)), literal_getproperty_pullback
end

Zygote.accum(x::Adjoint...) = Adjoint(Zygote.accum(map(parent, x)...))

Zygote.accum(x::NamedTuple{(:parent,)}, y::Adjoint) = (parent=accum(x.parent, y.parent),)

function Zygote.accum(A::UpperTriangular{<:Any, <:SMatrix{P}}, B::SMatrix{P, P}) where {P}
    return Zygote.accum(SMatrix{P, P}(A), B)
end

function Zygote.accum(B::SMatrix{P, P}, A::UpperTriangular{<:Any, <:SMatrix{P}}) where {P}
    return Zygote.accum(B, SMatrix{P, P}(A))
end

function Zygote.accum(a::Tangent{T}, b::NamedTuple) where {T}
    return Zygote.accum(a, Tangent{T}(; b...))
end

Base.:(+)(::Tangent, ::Nothing) = ZeroTangent()

function Base.:(-)(
    A::UpperTriangular{<:Real, <:SMatrix{N, N}}, B::Diagonal{<:Real, <:SVector{N}},
) where {N}
    return UpperTriangular(A.data - B)
end

function _symmetric_back(Δ, uplo)
    L, U, D = LowerTriangular(Δ), UpperTriangular(Δ), Diagonal(Δ)
    return collect(uplo == Symbol(:U) ? U .+ transpose(L) - D : L .+ transpose(U) - D)
end
_symmetric_back(Δ::Diagonal, uplo) = Δ
_symmetric_back(Δ::UpperTriangular, uplo) = collect(uplo == Symbol('U') ? Δ : transpose(Δ))
_symmetric_back(Δ::LowerTriangular, uplo) = collect(uplo == Symbol('U') ? transpose(Δ) : Δ)

function Zygote._pullback(
    ctx::AContext, ::Type{Symmetric}, X::StridedMatrix{<:Real}, uplo=:U,
)
    function Symmetric_pullback(Δ)
        ΔX = Δ === nothing ? nothing : _symmetric_back(Δ, uplo)
        return nothing, ΔX, nothing
    end
    return Symmetric(X, uplo), Symmetric_pullback
end

# function Zygote._pullback(cx::AContext, ::typeof(literal_indexed_iterate), xs::Tuple, ::Val{i}) where i
#   y, b = Zygote._pullback(cx, literal_getindex, xs, Val(i))
#   back(::Nothing) = nothing
#   back(ȳ) = b(ȳ[1])
#   (y, i+1), back
# end

# function Zygote._pullback(cx::AContext, ::typeof(literal_indexed_iterate), xs::Tuple, ::Val{i}, st) where i
#   y, b = Zygote._pullback(cx, literal_getindex, xs, Val(i))
#   back(::Nothing) = nothing
#   back(ȳ) = (b(ȳ[1])..., nothing)
#   (y, i+1), back
# end

# Zygote._pullback(cx::AContext, ::typeof(getproperty), x, f::Symbol) =
#   Zygote._pullback(cx, Zygote.literal_getproperty, x, Val(f))

# Zygote._pullback(cx::AContext, ::typeof(getfield), x, f::Symbol) =
#   Zygote._pullback(cx, Zygote.literal_getproperty, x, Val(f))

# Zygote._pullback(cx::AContext, ::typeof(literal_getindex), x::NamedTuple, ::Val{f}) where f =
#   Zygote._pullback(cx, Zygote.literal_getproperty, x, Val(f))

# Zygote._pullback(cx::AContext, ::typeof(literal_getproperty), x::Tuple, ::Val{f}) where f =
#   Zygote._pullback(cx, Zygote.literal_getindex, x, Val(f))

# function Zygote._pullback(
#     ::AContext,
#     T::Type{<:StructArray{T, N, C} where {T, N, C<:NamedTuple}},
#     x::Union{Tuple, NamedTuple},
# )
#     function StructArray_pullback(Δ::NamedTuple{(:components, )})
#         @show typeof(x), typeof(Δ.components)
#         return (nothing, Δ.components)
#     end
#     return T(x), StructArray_pullback
# end

function Zygote._pullback(::AContext, T::Type{<:StructArray}, x::Tuple)
    function StructArray_pullback(Δ::NamedTuple{(:components, )})
        return (nothing, values(Δ.components))
    end
    return T(x), StructArray_pullback
end

# `getproperty` accesses the `components` field of a `StructArray`. This rule makes that
# explicit.
function Zygote._pullback(
    ctx::AContext, ::typeof(Zygote.literal_getproperty), x::StructArray, ::Val{p},
) where {p}
    value, pb = Zygote._pullback(
        ctx, Zygote.literal_getproperty, getfield(x, :components), Val(p),
    )
    function literal_getproperty_pullback(Δ)
        return nothing, (components=pb(Δ)[2], ), nothing
    end
    return value, literal_getproperty_pullback
end

function time_ad(label::String, f, x...)
    println("primal: ", label)
    return @time f(x...)
end

time_ad(::Val{:disabled}, label::String, f, x...) = f(x...)

function Zygote._pullback(ctx::AContext, ::typeof(time_ad), label::String, f, x...)
    println("Forward: ", label)
    out, pb = @time Zygote._pullback(ctx, f, x...)
    function time_ad_pullback(Δ)
        println("Pullback: ", label)
        Δinputs = @time pb(Δ)
        return (nothing, nothing, Δinputs...)
    end
    return out, time_ad_pullback
end

function Zygote._pullback(ctx::AContext, ::typeof(\), A::Diagonal{<:Real}, x::Vector{<:Real})
    out, pb = Zygote._pullback(ctx, (a, x) -> a .\ x, diag(A), x)
    function ldiv_pullback(Δ)
        if Δ === nothing
            return nothing
        else
            _, Δa, Δx = pb(Δ)
            return nothing, Diagonal(Δa), Δx
        end
    end
    return out, ldiv_pullback
end

function Zygote._pullback(ctx::AContext, ::typeof(\), A::Diagonal{<:Real}, x::Matrix{<:Real})
    out, pb = Zygote._pullback(ctx, (a, x) -> a .\ x, diag(A), x)
    function ldiv_pullback(Δ)
        if Δ === nothing
            return nothing
        else
            _, Δa, Δx = pb(Δ)
            return nothing, Diagonal(Δa), Δx
        end
    end
    return out, ldiv_pullback
end

using Base.Broadcast: broadcasted

function Zygote._pullback(
    ::AContext, ::typeof(broadcasted), ::typeof(\), a::Vector{<:Real}, x::Vector{<:Real},
)
    y = a .\ x
    function broadcast_ldiv_pullback(Δ::Union{Nothing, Vector{<:Real}})
        if Δ === nothing
            return nothing
        else
            return nothing, nothing, -(Δ .* y ./ a), a .\ Δ
        end
    end
    return y, broadcast_ldiv_pullback
end

function Zygote._pullback(
    ::AContext, ::typeof(broadcasted), ::typeof(\), a::Vector{<:Real}, x::Matrix{<:Real},
)
    y = a .\ x
    function broadcast_ldiv_pullback(
        Δ::Union{Nothing, Matrix{<:Real}, Adjoint{<:Real, <:Matrix{<:Real}}},
    )
        if Δ === nothing
            return nothing
        else
            return nothing, nothing, -vec(sum(Δ .* y ./ a; dims=2)), a .\ Δ
        end
    end
    return y, broadcast_ldiv_pullback
end


text:
@doc raw"""
```
Z,p,q = retrieve_surface(algorithm::DiscreteShapeBound, img::AbstractArray, iterations::Int=2000)
```

Same as [`DiscreteShape`](@ref retrieve_surface(::DiscreteShape, ::AbstractArray, ::iterations))
except it has its initial conditions bound by the image as per bellow where
``E`` is the brightness of the image:
```math
\begin{gathered}
Z_{i,j}=\begin{cases}-100E_{i,j} &\text{if } E_{i,j}>0.75\\0 &\text{otherwise}\\
\end{cases}\\\\p,q=\nabla E
\end{gathered}
```
# Example
Compute the heightmap for a synthetic image generated by `generate_surface`.
```julia
using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(), 1, [0.2,0,0.9], radius = 5)

#calculate the heightmap (using 500 iterations)
Z,p,q = retrieve_surface(DiscreteShapeBound(), img, 500)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)
```
# Reference
1. S. Elhabian, "Hands on Shape from Shading", Computer Vision and Image Processing, 2008.
"""
function retrieve_surface(algorithm::DiscreteShapeBound, img::AbstractArray, iterations::Int=2000; smoothness::Int=1000)
    ρ,I,σ,τ = estimate_img_properties(img)
    λ = smoothness
    return retrieve_surface(DiscreteShapeBound(), img, ρ, I, iterations, smoothness=λ)
end

function retrieve_surface(algorithm::DiscreteShapeBound, img::AbstractArray, albedo::Real, illumination_direction::Vector{T} where T <: Real, iterations::Int=2000; smoothness::Int=1000)
    ρ = albedo
    I = illumination_direction
    E = Array{Float64}(img)

    #initialize variables
    p = zeros(Complex{Float64},axes(E))
    q = zeros(Complex{Float64},axes(E))
    R = zeros(Complex{Float64},axes(E))
    Z = zeros(axes(E))

    #apply bounding conditions
    for i in CartesianIndices(E)
        if E[i] > 0.75
            Z[i] = -100*E[i]
        else
            Z[i] = 0.0
        end
    end
    q, p = Array{Complex{Float64}}.(imgradients(Z, KernelFactors.sobel, "replicate"))
    return solve_EulerLagrange(ρ, I, iterations, p, q, R, smoothness, E, Z)
end

ext: jl


text:
type CVSystem
    Pa::Vector{Float64}
    Pv::Vector{Float64}
    V::Vector{Float64}
    E::Vector{Float64}
    Q::Vector{Float64}
    t::Vector{Float64}
    sparams::SolverParams
    mparams::ModelParams

    function CVSystem(nvar=0,old=Dict("a"=>0),restart="no")
        this = new()
        if restart == "no"
            this.mparams = ModelParams();
        elseif restart == "yes"
            mparams = old["system"]["mparams"];
            this.mparams = ModelParams(mparams,restart);
        end
        this.Pa = [];
        this.V = [];
        this.Q = [];
        this.Pv = [];
        this.E = [];
        this.t = [];
        this.sparams = SolverParams(nvar);
        return this
    end
end


text:
"""
$SIGNATURES

Internal function to check that two vectors have matching length.
"""
function check_dims(vec1, vec2)
    if length(vec1) != length(vec2)
        throw(DimensionMismatch("The vectors don't have the same length."))
    end
end

"""
$SIGNATURES

Internal function to unpack kwargs.
"""
dict(;kw...) = isempty(kw) ? nothing : LittleDict{Symbol,Any}(kw)

"""
$SIGNATURES

Check if a function has already been defined on the board and is available
by checking the current functions in reverse order. Note that you could fool
this system but it should be ok in most standard cases.
"""
function isdef(f, b)
    isnothing(b) && return false
    fl = findlast(fc -> fc.name == f.name, b.functions)
    isnothing(fl) && return false
    fc = b.functions[fl]
    return fc.s == f.s
end


text:
module DSUtils

using LinearAlgebra
using FreqTables
using NamedArrays

using Reexport
@reexport using StatsBase
@reexport using DataFrames
@reexport using Plots

export BCDiag, bcdiag
export biasplot, ksplot, rocplot, accuracyplot
export liftcurve, cumliftcurve
export liftable, cumliftable

export onehot, onehot!

export ranks                # equal count binning
export strstd               # string standardization
export sumxm                # element-wise sum of vectors with missing as 0

export kstest               # two sample KS and location
export auroc                # concordance and AUROC
export concordance          # Concordance with bounds function

export infovalue

include("bcdiag.jl")
include("onehotencoder.jl")
include("utils.jl")

end # module


text:
import Base: in, view, size, eltype, names

using Dates
using HDF5

# TODO: remove duplication with Sario constants
const MASK_FILENAME = "masks.h5"
const UNW_FILENAME = "unw_stack.h5"
const CC_FILENAME = "cc_stack.h5"

# dataset names for general 3D stacks
const STACK_DSET = "stack"
const STACK_MEAN_DSET = "stack_mean"
const STACK_FLAT_DSET = "stack_flat"
const STACK_FLAT_SHIFTED_DSET = "stack_flat_shifted"

# Mask file datasets
const GEO_MASK_DSET = "geo"
const GEO_MASK_SUM_DSET = "geo_sum"
const IGRAM_MASK_DSET = "igram"

const DEM_RSC_DSET = "dem_rsc"
const GEOLIST_DSET = "geo_dates"
const INTLIST_DSET = "int_dates"

const REFERENCE_ATTR = "reference"
const REFERENCE_STATION_ATTR = "reference_station"

# Type alias for commonly used compositite type
const Igram = Tuple{Date,Date}

# To check if either date of an igram is contained with a geo date list
Base.in(igram::Igram, geo_date_list::AbstractArray{Date}) =
    ((igram[1] in geo_date_list) || (igram[2] in geo_date_list))

Base.in(date::Date, igram_list::AbstractArray{Igram}) =
    [date in igram for igram in igram_list]

temporal_baseline(igram::Igram) = (igram[2] - igram[1]).value
temporal_baseline(igram_list::Array{Igram}) =
    [temporal_baseline(igram) for igram in igram_list]
span(dates::AbstractArray{Date,1}) = (dates[end] - dates[1]).value
span(igrams::AbstractArray{Igram,1}) = (sort(igrams)[end][2] - sort(igrams)[1][1]).value

_get_day_nums(dts::AbstractArray{Date,1}) = [(d - dts[1]).value for d in dts]

# Extend views to work for datasets
function Base.view(dset::HDF5Dataset, I::Vararg{Any})
    idxs = to_indices(dset, I)
    return dset[idxs...]
end

function Base.view(dset::HDF5Dataset, i::Colon, j::Colon, k)
    # Trick to get the slice view of a 3D dset into 2D instead of (m, n, 1)
    return view(dset[i, j, k], :, :)
end

function Base.view(dset::HDF5Dataset, i::Int, j::Int, k::Colon)
    # Hack for 1D depth slice to have 1D instead of (1, 1, N)
    return view(dset[i, j, k], 1, 1, :)
end

function eltype(h5file::String, dset::String)
    h5open(h5file) do f
        return eltype(f[dset])
    end
end

# Allow the HDF5 `names` to work on just filenam
function names(fname::AbstractString)
    h5open(fname) do f
        return names(f)
    end
end
function names(fname::AbstractString, group::AbstractString)
    h5open(fname) do f
        return names(f[group])
    end
end

function get_chunk(fname::AbstractString, dset::AbstractString)
    h5open(fname) do f
        return HDF5.get_chunk(f[dset])
    end
end

# From https://github.com/JuliaLang/julia/blob/master/base/broadcast.jl#L662
# broadcastable(x::Union{AbstractArray,Number,Ref,Tuple,Broadcasted}) = x
# # Default to collecting iterables — which will error for non-iterables
# broadcastable(x) = collect(x)
Base.broadcastable(x::HDF5Dataset) = read(x)
Base.collect(x::HDF5Dataset) = collect(read(x))
Base.iterate(x::HDF5Dataset) = iterate(read(x))
Base.iterate(x::HDF5Dataset, state) = iterate(read(x), state)

"""Get the number of megabytes of RAM available
Note this is higher than "Sys.free_memory()"
"""
function getmemavail()
    out = Pipe()
    lines = read(`free -m`, String)
    memline = split(lines, '\n')[2]
    return parse(Float64, split(memline)[end])
end


text:
module Sesam

using ModelingToolkit, DifferentialEquations, IfElse
#using Unitful
using DistributionFits: fit_mode_flat
using Distributions
using MTKHelpers

export plant_const, plant_face, plant_face_fluct, plant_face_fluct_fake,
    sesam3, seam3, sesam3CN, plant_sesam_system,
    sesam3_protect, sesam3_revMM
export calculate_β_NR_sesam3, calculate_β_PR_sesam3
#export smoothstep
#export SystemParUpdater

#include("units.jl")
include("util.jl")
include("sesam3_revMM.jl")
include("sesam3.jl")
include("sesam3P.jl")
include("seam3.jl")
include("sesam3_protect.jl")
include("plant_const.jl")
include("plant_face.jl")
include("plant_face_fluct.jl")
include("compose.jl")

end


text:
"""
# Moon

This module defines seismic models for the Earth's moon.
"""
module Moon

using ..SeisModels

"""
	ALL_MODELS

List of Moon model instances as `Symbol`s.

If adding new inbuilt Moon models, append the name of the new model to this list
"""
const ALL_MODELS = (:MOON_WEBER_2011,)

include("weber_2011.jl")

end # module


text:
using DataDeps

register(DataDep(
    "synthetic",
    """
    Dataset: synthetic

    A simulated multivariate normal.
    Classes are determined by the first two dimensions, which correspond to signs of X*beta with beta[1] = 10, beta[2] = -10.
    Covariance structure is as follows:

    Σ[1,1] = 1, Σ[2,2] = 3, Σ[i,i] = 2
    Σ[1,2] = 0.9
    Σ[i,j] = 1e-3*randn()

    Observations: 1000
    Features:     500
    Classes:      2
    """,
    "script", # nothing to download
    "514ef72019b85634859baf459b38f37295d5fcc84d06be7d8e4a71a99a74cbdf";
    fetch_method=function(unused, localdir)
        #
        rng = StableRNG(2000)
        (m, n) = (10^3, 500)

        # covariance matrix
        Σ = Matrix{Float64}(2*I, n, n)
        for j in 1:n, i in j+1:n
            Σ[j,i] = 1e-3*randn(rng)
        end
        Σ[1,1] = 1.0
        Σ[2,2] = 3.0
        Σ[1,2] = 0.9

        L, _ = cholesky(Symmetric(Σ))

        # simulate data
        X = Matrix{Float64}(undef, m, n)
        for i in axes(X, 1)
            @views X[i, :] .= L*randn(rng, n)
        end

        # coefficients
        beta = zeros(n)
        beta[1] = 10.0
        beta[2] = -10.0

        # targets
        y = sign.(X*beta)
        target = map(yi -> yi > 0 ? 'A' : 'B', y)

        local_file = joinpath(localdir, "data.csv")
        df = hcat(DataFrame(target=target), DataFrame(X, :auto))
        CSV.write(local_file, df)
        return local_file
    end
))


text:
module Interpolation
using Images
using Interpolations
import Interpolations.Throw # see https://github.com/JuliaMath/Interpolations.jl/issues/127

using DeformableRegistration.Transformation
using DeformableRegistration.ImageProcessing

export interpolateImage,InterpLinearFast
export interpolateDeformationField

function interpolateImage(
    I::regImage,transformedGrid::scaledArray;
    doDerivative=false,interpolationScheme=InterpLinearFast)

    return interpolateImage(I, transformedGrid.data, doDerivative=doDerivative,
                            interpolationScheme=interpolationScheme)
end

function interpolateImage(
    I::regImage,transformedGrid::Array{Float64,1};
    doDerivative=false,interpolationScheme=InterpLinearFast)

    return interpolateArray(Array(I.data), I.voxelsize, I.shift, transformedGrid, size(I.data);
                            doDerivative=doDerivative,interpolationScheme=interpolationScheme)
end

function interpolateArray(
    data::Array{Float64,2}, voxelsize::Array{Float64,1},
    shift::Array{Float64,1}, transformedGrid::Array{Float64,1}, newDimensions::Tuple{Vararg{Int64}};
    doDerivative=false, interpolationScheme=InterpLinearFast)

  # use faster linear interpolation as default
  if(interpolationScheme == InterpLinearFast)
    return InterpLinearFast(data, voxelsize, shift, transformedGrid, newDimensions, doDerivative=doDerivative)
  end

  # determine number of new points, pixel spacing and spatial domain of the image
  numberOfPoints::Int = ceil(Int64,size(transformedGrid,1)/2)

  imageIntNoScaling = interpolate(data, interpolationScheme, OnCell())

  xRange = shift[1]+voxelsize[1]/2:voxelsize[1]:shift[1]+voxelsize[1]*(size(data,1))
  yRange = shift[2]+voxelsize[1]/2:voxelsize[2]:shift[2]+voxelsize[2]*(size(data,2))

  imageIntTmp = extrapolate(imageIntNoScaling, 0.0)
  imageInt = scale(imageIntTmp, xRange, yRange)

  transformedImage = zeros(numberOfPoints);
  for i=1:numberOfPoints
      transformedImage[i] = imageInt[transformedGrid[i],transformedGrid[numberOfPoints+i]]
  end

  if(doDerivative)
    imageIntTmp = extrapolate(imageIntNoScaling, Throw())
    imageInt = scale(imageIntTmp, xRange, yRange)

    dY_transformedImage = zeros(numberOfPoints)
    dX_transformedImage = zeros(numberOfPoints)
    for i=1:numberOfPoints
        try # zero if outside boundary (extrapolate(..., 0.0) does not work with gradient)
            dX_transformedImage[i],dY_transformedImage[i] = gradient(imageInt,transformedGrid[i],transformedGrid[i+numberOfPoints])
        catch x
            if isa(x, BoundsError)
                #print(".")
            else
                throw(x)
            end
        end
    end
    return reshape(transformedImage, newDimensions), dX_transformedImage, dY_transformedImage
  end
  return reshape(transformedImage, newDimensions)
end

function InterpLinearFast(data::Array{Float64,2}, voxelsize::Array{Float64,1}, shift::Array{Float64,1},
                            transformedGrid::Array{Float64,1}, newDimensions::Tuple{Vararg{Int64}};doDerivative=false)

  # determine number of new points, pixel spacing and spatial domain of the image
  numberOfPoints::Int = convert(Int64,(size(transformedGrid,1)/2))

  # setup variables
  transformedImage = zeros(numberOfPoints)
  dY_transformedImage = zeros(numberOfPoints)
  dX_transformedImage = zeros(numberOfPoints)
  p = zeros(4)
  sizeY= size(data)[2]
  sizeX = size(data)[1]
  x::Float64 = 0.0; y::Float64 = 0.0;
  xf::Int = 0; yf::Int = 0;

  # interpolate image at new points
  for i=1:numberOfPoints

    x = (transformedGrid[i]   - shift[1]) / voxelsize[1] + .5
    y = (transformedGrid[numberOfPoints+i] - shift[2]) / voxelsize[2] + .5

    if ((x <= 0) | (y <= 0) | (x >= sizeX+1) | (y >= sizeY+1))
      continue
    end

    xf = floor(Int64,x)
    yf = floor(Int64,y)

    p[1] = ((xf<1)          | (yf<1))        ?   0.0: data[xf,yf]
    p[2] = ((xf+1>sizeX)    | (yf<1))        ?   0.0: data[xf+1,yf]
    p[3] = ((xf<1)          | (yf+1>sizeY))  ?   0.0: data[xf,yf+1]
    p[4] = ((xf+1>sizeX)    | (yf+1>sizeY))  ?   0.0: data[xf+1,yf+1]

    x = x - xf
    y = y - yf

    transformedImage[i] = (p[1] * (1-x) + p[2] * x) * (1-y) + (p[3] * (1-x) + p[4] * x) * y
    dX_transformedImage[i] = ((p[3] - p[1]) * (1-x) + (p[4] - p[2]) * x) / voxelsize[1]
    dY_transformedImage[i] = ((p[2] - p[1]) * (1-y) + (p[4] - p[3]) * y) / voxelsize[2]

  end

  if(doDerivative)
    return reshape(transformedImage,newDimensions),dY_transformedImage,dX_transformedImage
  else
    return reshape(transformedImage,newDimensions)
  end
end

function interpolateDeformationField(deformationField::scaledArray, newPoints::scaledArray; interpolationScheme = InterpLinearFast)# BSpline(Linear())
    dims = deformationField.dimensions
    deformationFieldX = reshape(deformationField.data[1:prod(dims)], dims[1], dims[2])
    deformationFieldY = reshape(deformationField.data[prod(dims)+1:end], dims[1], dims[2])

    return scaledArray([
        interpolateArray(deformationFieldX, deformationField.voxelsize, deformationField.shift,
            newPoints.data, (prod(newPoints.dimensions),), interpolationScheme = interpolationScheme, doDerivative=false);
        interpolateArray(deformationFieldY, deformationField.voxelsize, deformationField.shift,
            newPoints.data, (prod(newPoints.dimensions),), interpolationScheme = interpolationScheme, doDerivative=false)
        ],newPoints.dimensions, newPoints.voxelsize, newPoints.shift)

end

end


text:
using Test, TSVD, LinearAlgebra, SparseArrays

@testset "test tsvd with m = $m, n = $n, and p = $p" for
    (m, n, p) in ((10, 6, 0.8),
                  (6, 10, 0.8),
                  (10, 10, 0.8),
                  (100, 60, 0.1),
                  (60, 100, 0.1),
                  (100, 100, 0.1))

    mnp = round(Integer, m*n*p)

    for A in (randn(m, n),
              complex.(randn(m, n), randn(m, n)),
              sprandn(m, n, p),
              sparse(rand(1:m, mnp), rand(1:n, mnp), complex.(randn(mnp), randn(mnp)), m, n))

        Uf, sf, Vf = svd(Array(A))

        for k = 1:5
            U, s, V = TSVD.tsvd(A, k)
            @test norm(s - sf[1:k]) < sqrt(eps(real(eltype(A))))*mnp
            @test norm(abs.(U'Uf[:,1:k]) - I) < sqrt(eps(real(eltype(A))))*mnp
            @test norm(abs.(V'Vf[:,1:k]) - I) < sqrt(eps(real(eltype(A))))*mnp

            s, V = TSVD.tsvd2(A, k)
            # tmp = TSVD._teig(TSVD.AtA(A, randn(n)), k, debug = true)
            # @test norm(sqrt(reverse(tmp[1]))[1:k] - sf[1:k]) < sqrt(eps(real(eltype(A))))*mnp
            @test norm(s - sf[1:k]) < sqrt(eps(real(eltype(A))))*mnp
            @test norm(abs.(V'Vf[:,1:k]) - I) < sqrt(eps(real(eltype(A))))*mnp
        end
    end
end

@testset "Issue 9" begin
    data = rand(1:100, 50, 50)
    @test @inferred(TSVD.tsvd(data, 2))[2] ≈ svdvals(data)[1:2]
end

@testset "Issue 25. Degenerate basis" begin
  @test tsvd([1 0 ; 0 1], 1, initvec=[1.0; 0.0])[2] ≈ [1.0]
end


text:
immutable InterpolationData{F,uType,tType,kType,cacheType} <: Function
  f::F
  timeseries::uType
  ts::tType
  ks::kType
  notsaveat_idxs::Vector{Int}
  cache::cacheType
end

immutable CompositeInterpolationData{F,uType,tType,kType,cacheType} <: Function
  f::F
  timeseries::uType
  ts::tType
  ks::kType
  alg_choice::Vector{Int}
  notsaveat_idxs::Vector{Int}
  cache::cacheType
end

(interp::InterpolationData)(tvals,idxs,deriv) = ode_interpolation(tvals,interp,idxs,deriv)
(interp::CompositeInterpolationData)(tvals,idxs,deriv) = ode_interpolation(tvals,interp,idxs,deriv)
(interp::InterpolationData)(val,tvals,idxs,deriv) = ode_interpolation!(val,tvals,interp,idxs,deriv)
(interp::CompositeInterpolationData)(val,tvals,idxs,deriv) = ode_interpolation!(val,tvals,interp,idxs,deriv)

function InterpolationData(id::InterpolationData,f)
  InterpolationData(f,id.timeseries,
                      id.ts,
                      id.ks,
                      id.notsaveat_idxs,
                      id.cache)
end

function CompositeInterpolationData(id::CompositeInterpolationData,f)
  CompositeInterpolationData(f,id.timeseries,
                               id.ts,
                               id.ks,
                               id.alg_choice,
                               id.notsaveat_idxs,
                               id.cache)
end


text:
module PaperAnalysis
    using OceanWaveSpectralFitting, Optim, DSP
    include("nonparametrics.jl")

    const WindSea = JS_BWG_HNE_DL{0,40} # no aliasing and 40m water depth

    "
        movingaverage(vector, half_window_size)

    Compute the moving average of a vector"
    function movingaverage(vector::Vector{T}, half_window_size) where {T}
        out = similar(vector)
        for i = 1:length(vector)
            boundarydistance = min(length(vector) - i, i - 1)
            window = min(half_window_size, boundarydistance)
            @views out[i] = sum(vector[i-window:i+window])/(2window+1)
        end
        return out
    end

    "
        initialisepeak(ts, Δ, ωₗ, ωᵤ)

    Compute the peak frequency and direction for a certain frequency band."
    function initialisepeak(ts::Matrix{Float64}, Δ, ωₗ, ωᵤ)
        P = OceanWaveSpectralFitting.Periodogram(ts,Δ)
        S = P.ordinates[end÷2+1:end]
        Ω = P.Ω[end÷2+1:end]
        peakind = argmax([(ωₗ < ω < ωᵤ) ? real(s[1,1]) : -Inf for (s,ω) ∈ zip(S, Ω)])
        meandir = angle(imag(S[peakind][2,1])+1im*imag(S[peakind][3,1]))
        return Ω[peakind], mod2pi.(meandir.+pi) .- pi
    end

    "Function to generate an initial guess when we have a wind-sea only model."
    function initialguess(::Type{<:WindSea},ts,Δ,cutoffs)
        peak,mdir = initialisepeak(ts,Δ,maximum(cutoffs),4)
        return [0.7, peak, 2, 5, mdir, 4, 2.7, 0.55, 0.26]
    end

    "Type to store the timeseries, initial parameters, cutoff and proposed model."
    struct InitialiseData
        ts::Matrix{Float64}
        Δ::Float64
        cutoff::Float64
        model::Type{<:TimeSeriesModel}
        θ₀::Vector{Float64}
        function InitialiseData(model::Type{<:TimeSeriesModel},ts::Matrix{Float64},Δ::Float64,cutoff)
            θ₀ = initialguess(model,ts,Δ,cutoff)
            new(ts,Δ,cutoff,model,θ₀)
        end
    end

    "Make constraints for a wind-sea only model."
    function makeconstraints(::Type{<:WindSea},cut)
        lxW = [0, cut, 1,  3,-Inf, 0,  0,   0,  0  ]
        uxW = [4, 4,   Inf,8, Inf, 2π, Inf ,Inf,Inf]
        return TwiceDifferentiableConstraints(lxW,uxW)
    end

    "Type to store the results of fitting one choice of cutoff."
    struct FittingResults
        θ_est::Vector{Float64}
        varθ::Matrix{Float64}
        converged::Bool
        function FittingResults(ini,iterations)
            obj = DebiasedWhittleLikelihood(ini.model, ini.ts, ini.Δ; lowerΩcutoff = ini.cutoff, upperΩcutoff = 3.8)
            constraints = makeconstraints(ini.model, ini.cutoff)
            res = optimize(TwiceDifferentiable(Optim.only_fgh!(obj), ini.θ₀), constraints, ini.θ₀, IPNewton(), Optim.Options(iterations = iterations))
            obj = DebiasedWhittleLikelihood(ini.model, ini.ts, ini.Δ; lowerΩcutoff = ini.cutoff, upperΩcutoff = 3.8)
            EH = zeros(length(ini.θ₀),length(ini.θ₀))
            obj(nothing,nothing,EH,res.minimizer)
            varθ = nothing
            try
                varθ = inv(EH)
            catch
                varθ = fill(NaN, size(EH))
            end
            new(res.minimizer, varθ, Optim.converged(res))
        end
    end

    "Type to store fitted models and the initial values."
    struct FittedModel
        ini::InitialiseData
        fits::FittingResults
        function FittedModel(ini::InitialiseData;iterations=3000)
            fits = FittingResults(ini,iterations)
            new(ini,fits)
        end
    end
end


text:
using Documenter

# using Pkg
# Pkg.activate("..")

using FinancialSymbology

makedocs(
    sitename = "FinancialSymbology",
    format = Documenter.HTML(),
    modules = [FinancialSymbology],
    doctest = true
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#

max_stars_repo_name: robfs/figi

max_issues_repo_name: robfs/figi

max_forks_repo_name: robfs/figi


text:
### Day 13 - Problem 2
## A Maze of Twisty Little Cubicles
## Author: Thanasis Georgiou

# Cell type
immutable Cell
  x::Int
  y::Int
  weight::Int
end

# Equality check for cells
function Base.isequal(a::Cell, b::Cell)::Bool
  return a.x == b.x && a.y == b.y
end

# Manhattan distance from cell to cell
function man_distance(a::Cell, b::Cell)::Int
  return abs(a.x - b.x) + abs(a.y - b.y)
end

# Manhattan distance from coords to cell
function man_distance(x::Int, y::Int, b::Cell)::Int
  return abs(x - b.x) + abs(y - b.y)
end

createHeuristicCell(x::Int, y::Int, target::Cell) = Cell(x, y, man_distance(x, y, target))

# -1 For wall, anything else is step count
floor = Array(Int, (50, 50))

# Create floor plan
for x = 0:49, y = 0:49
  magic = x * x + 3 * x + 2 * x * y + y + y * y
  magic += 1358
  bitCount = count(bit -> bit == '1', bits(magic))

 floor[x + 1, y + 1] = isodd(bitCount)? -1 : 0
end

# Find optimal path
target = Cell(31 + 1, 39 + 1, 0)
visited = Array(Cell, 0)
queue = Array(Cell, 0)
from = Dict{Cell, Cell}()
path = Array(Cell, 0)

# Add starting position to queue
push!(queue, createHeuristicCell(2, 2, target))

while length(queue) > 0
  # Sort queue
  sort!(queue, lt = (a, b) -> a.weight < b.weight)

  # Grab most promising cell
  cell = queue[1]
  push!(visited, cell)
  deleteat!(queue, 1)
  steps = floor[cell.x, cell.y]

  if (steps > 50)
    continue
  end

  # Find neighbors
  # For each neighbor check it's a wall. If not, create a cell, push to queue
  # add to dictionary
  neighbors = [
    createHeuristicCell(cell.x + 1, cell.y, target),
    createHeuristicCell(cell.x - 1, cell.y, target),
    createHeuristicCell(cell.x, cell.y + 1, target),
    createHeuristicCell(cell.x, cell.y - 1, target),
  ]
  for neighbor in neighbors
    if neighbor.x > 0 && neighbor.y > 0 && floor[neighbor.x, neighbor.y] == false
      # Check if we have visited this before
      visitedIndex = findfirst(x -> isequal(x, neighbor), visited)
      if visitedIndex == 0
        push!(queue, neighbor)
        from[neighbor] = cell
        floor[neighbor.x, neighbor.y] = steps + 1
      else
        oldCell = visited[visitedIndex]
        if floor[oldCell.x, oldCell.y] > steps + 1
          deleteat!(visited, visitedIndex)
          from[neighbor] = cell
          floor[neighbor.x, neighbor.y] = steps + 1
        end
      end
    end
  end
end

tileCount = count(x -> x > 0 && x <= 50, floor) + 1 # Plus one for starting position
println("Finished execution, reachable tiles: $tileCount")

for x = 1:40
  for y = 1:45
    if floor[x, y] > 0 && floor[x, y] <= 50
      print("x")
    elseif floor[x, y] == -1
      print("█")
    else
      print(" ")
    end
  end
  println()
end


text:
include("../mesh/mesh_clear.jl")

@testset "Meshgrid" begin

end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
println("Starting -- Julia seems to be working.")
println("")

using Random
Random.seed!(1)

###############################################################
println("Let's first create an initial cluster, with 16 star formation events")

include("./Cluster.jl")

lmax = 10

n_disks = 16
n_stars_per_disk = 64
σ = 0.006

m_min = 1.0
m_max = 100.0
γₘ = -2.0

a_min = 1.0
a_max = 100.0
γₐ = 0.0

e_min = 0.0
e_max = 0.3
γₑ = 1.0

clust = Cluster_initial(n_disks, n_stars_per_disk, lmax,
                        σ = σ,
                        m_min = m_min, m_max = m_max, γₘ = γₘ,
                        a_min = a_min, a_max = a_max, γₐ = γₐ,
                        e_min = e_min, e_max = e_max, γₑ = γₑ)

println("We can compute this cluster's binding energy E and spin s")
Etot_norm, s = get_ES_from_clust(clust) # Here, Etot_norm = Etot / (J N^2) = -E
                                        # J = G mₘᵢₙ² / aₘᵢₙ

###############################################################
println("The cluster has been created and its (E, s) computed")
println("Now let's plot it")
# It should reproduce the figure "Initial_cluster.png"

using Plots

my_plot_1 = plot_that_cluster(clust)

display(my_plot_1)

###############################################################
println("")
println("Now let's create the initial DF related to this initial cluster")
# This means we create all the Kₖ = (m, a, e), Nₖ, Lₖ and Jₗ[K, K'], and we translate (E, s) from the initial cluster into (Eₜₒₜ, Lₜₒₜ)

# First we create the m, a, e bins, and fill them with stars i.e. create the N(K) vector
include("./Miscellaneous.jl")

n_pop_m = 10
n_pop_a = 3
n_pop_e = 1

Bins_m, Bins_a, Bins_e, K, L, N = create_KLN(m_min, m_max, γₘ, n_pop_m,
                                            a_min, a_max, γₐ, n_pop_a,
                                            e_min, e_max, γₑ, n_pop_e)

# Then we create the matrix of coupling coefficients J
include("./Coupling_coefficients_interpolated.jl")

J = Coupling_coefficients_interpolated(K, lmax)

# Then we translate the normalized binding energy and spin into the pure energy and angular momentum
include("./Distribution.jl")

Etot, Ltot = get_EL_for_dist(Etot_norm, s,
                                L, N, m_min, a_min)

###############################################################
println("The initial DF realted to the initial cluster has been created")
println("")
println("Let's modify this DF so that it maximizes the entropy S")
# This means we relax the cluster

include("./Solve_by_energy_sequence.jl")
using LinearAlgebra

dist = solve_by_energy_sequence(J, L, N, Etot, Ltot)
println( norm(dist.cost) < 10^(-6) ) # This is to check that the optimization method converged

###############################################################
println("The entropy has been maximized")
println("Finally, let's plot the final DF")
# It should reproduce the figure "Relaxed_cluster.png"

my_plot_2 = plot_that_dist("m", # "m" means that we want to plot the segregation in mass.
                                # "a" or "e" would give the segregations in semi-major axis and eccentricity
                        Etot_norm, s, lmax,
                        dist,
                        Bins_m, n_pop_m,
                        Bins_a, n_pop_a,
                        Bins_e, n_pop_e)

display(my_plot_2)

readline() # This is just so that the figure panels remain open

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe


text:
function calc_width_ngon(a,b)
if length(a) == length(b) && length(a) >= 2
    # SIDES
    v = zeros(length(a)+1);
    v[1] = abs(a[1]+im*b[1]);
    for i = 2:length(a)
        v[i] = abs((a[i]-a[i-1])+im*(b[i]-b[i-1]));
    end
    v[length(a)+1] = abs(a[length(a)]+im*b[length(a)]);
    # HEIGHTS
    h = zeros(length(a)+1);
    # h1
    j = 1;
    while j <= length(a)
        h1 = abs(a[j]*(b[1]-b[j])-b[j]*(a[1]-a[j]))/v[1];
        if h1 > h[1]
            h[1] = h1;
        end
        j = j+1;
    end
    # hi
    for i = 2:length(a)
        h[i] = abs(a[i-1]*b[i]-b[i-1]*a[i])/v[i];
        j = 1;
        while j <= length(a)
            hi = abs((a[i-1]-a[j])*(b[i]-b[j])-(b[i-1]-b[j])*(a[i]-a[j]))/v[i];
            if hi > h[i]
                h[i] = hi;
            end
            j = j+1;
        end
    end
    # hn
    j = 1;
    while j <= length(a)
        hn = abs((a[length(a)]-a[j])*b[j]-(b[length(a)]-b[j])*a[j])/v[length(a)+1];
        if hn > h[length(a)+1]
            h[length(a)+1] = hn;
        end
        j = j+1;
    end
    # WIDTH
    W = minimum(h);
    return W
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
module AutoGrid
export AutoGridConfig, AbstractFuncConfig, FuncConfig, SavingFuncConfig, WrappingFuncConfig, auto_grid

import LinearAlgebra: dot

abstract type AbstractAutoGridConfig end
abstract type AbstractFuncConfig{saving} end

nothingfun(x) = nothing

struct AutoGridConfig{N, FloatT<:Number, FuncT<:Function, XscaleT<:Function, XscaleInvT<:Function, FuncConfigs<:Tuple{Vararg{<:AbstractFuncConfig,N}}, LT<:Number} <: AbstractAutoGridConfig
    func::FuncT
    lims::NTuple{2,LT}
    xscale::XscaleT
    xscale_inv::XscaleInvT
    funcConfigs::FuncConfigs
    start_gridpoints::Int
    maxgridpoints::Int
    min_eps::FloatT
    reltol_dd::FloatT
    abstol_dd::FloatT
end

""" Standard config """
struct FuncConfig{saving, YscaleT<:Function, Diff2T, Diff3T} <: AbstractFuncConfig{saving}
    yscale::YscaleT
    find_peaks::Bool
    diff2::Diff2T
    diff3::Diff3T
    FuncConfig{saving}(yscale; find_peaks=true, diff2=normdiff2, diff3=normdiffdiff2) where saving =
        new{saving, typeof(yscale), typeof(diff2), typeof(diff3)}(yscale, find_peaks, diff2, diff3)
end
FuncConfig(args...; kwargs...) = FuncConfig{true}(args...; kwargs...)

""" Only saves the corresponding data, not used to computing grid """
struct SavingFuncConfig{saving} <: AbstractFuncConfig{saving}
    function SavingFuncConfig{saving}() where saving
        saving === true || throw(ArgumentError("SavingFuncConfig must have \"saving===true\" "))
        return new{true}()
    end
    SavingFuncConfig() = SavingFuncConfig{true}()
end
function Base.getproperty(sfc::SavingFuncConfig, sym::Symbol)
    if sym === :yscale
        return nothingfun
    elseif sym === :find_peaks
        return false
    else
        return getfield(sfc, sym)
    end
end

struct WrappingFuncConfig{saving, Diff2T, Diff3T, T<:Number} <: AbstractFuncConfig{saving}
    # yscale::YscaleT
    find_peaks::Bool
    wraplims::Tuple{T,T}
    diff2::Diff2T
    diff3::Diff3T
    WrappingFuncConfig{saving}(wraplims; find_peaks=true, diff2=normdiff2, diff3=normdiffdiff2) where saving=
        new{saving, typeof(diff2), typeof(diff3), typeof(wraplims[1])}(find_peaks, wraplims, diff2, diff3)
end
WrappingFuncConfig(args...; kwargs...) = WrappingFuncConfig{true}(args...; kwargs...)

function Base.getproperty(sfc::WrappingFuncConfig, sym::Symbol)
    if sym === :yscale
        return identity
    else
        return getfield(sfc, sym)
    end
end

AutoGridConfig(func::FuncT, lims::NTuple{2,LT}, xscale::XscaleT, xscale_inv::XscaleInvT, funcConfigs::FuncConfigs;
                start_gridpoints=30, maxgridpoints=2000,
                min_eps::FloatT=(xscale(lims[2])-xscale_inv(lims[1]))/100000,
                reltol_dd::FloatT=0.1, abstol_dd::FloatT=1e-6) where {N, FloatT<:Number, FuncT<:Function, XscaleT<:Function, XscaleInvT<:Function, FuncConfigs<:Tuple{Vararg{<:AbstractFuncConfig,N}}, LT<:Number} =
    AutoGridConfig{N, FloatT, FuncT, XscaleT, XscaleInvT, FuncConfigs, LT}(func, lims, xscale, xscale_inv, funcConfigs, start_gridpoints, maxgridpoints, min_eps, reltol_dd, abstol_dd)
"""
values, grid = auto_grid(func::Function, N, lims::Tuple, xscale=(identity, identity), FINCCONFIGSDEF; start_gridpoints, maxgridpoints, min_eps, reltol_dd, abstol_dd)
values, grid = auto_grid(agc::AutoGridConfig)

    Arguments:
    Compute a `grid` that tries to capture all features of the function `func` in the range `lim=(xmin,xmax)`.
    `func`: `Function` so that that returns a `Tuple` of values, e.g. `func(x) = (v1,v2,...,vN)` where each `vi` is a scalar or AbstractArray.
        note that the behavior for non scalar vi might be unintuitive.
    The algorithm will try to produce a grid so that features are captured in all of the values, but widely varying scaled might cause problems.
    `xscale`: `Tuple` containing a monotone function and its inverse corresponding to the default axis.
        Example: If you want to produce a plot with logarithmic x-axis, then you should set `xscale=(log10, exp10)`.
    `yscales`: `Tuple` containing one monotone function for each of the values `(v1,v2,...,vN)` corresponding to the default axis.

    Kwargs:
    `start_gridpoints`: The number of initial grid points.
    `maxgridpoints`: A warning will be thrown if this number of gridpoints are reached.
    `min_eps`: Lower limit on how small steps will be taken (in `xscale`).
        Example: with `xscale=(log10, exp10)` and `eps=1e-2`then `log10(grid[k+1])-log10(grid[k]) > (1e-2)/2`, i.e `grid[k+1] > 1.012grid[k]`.
    `reltol_dd = 0.05`, `abstol_dd = 0.01`: The criteria for further refinement is
        `norm(2*y2-y1-y3) < max(reltol_dd*max(norm(y2-y1), norm(y3-y2)), abstol_dd)`,
        where `y1,y2,y3` (in yscale) are 3 equally spaced points (in xscale).

    Output:
    `values`: `Tuple` with one vector per output value of `func`.
    `grid`: `Vector` with the grid corresponding to the values, i.e `func(grid[i])[j] = values[j][i]`

    Example: The following code computes a grid and plots the functions: abs(sinc(x)) and cos(x)+1
        in lin-log and log-log plots respectively.

    func = x -> (abs(sinc(x)),cos(x)+1)
    lims = (0.1,7.0)
    xscale = (log10, exp10)
    yscales = (identity, log10)
    y, x = auto_grid(func, lims, xscale, yscales, min_eps=1e-3)

    plot(x, y[1], m=:o; xscale=:log10, layout=(2,1), yscale=:identity, subplot=1, lab="sinc(x)", size=(800,600))
    plot!(x, y[2], m=:o, xscale=:log10, subplot=2, yscale=:log10, lab="cos(x)+1", ylims=(1e-5,10))

"""
# function auto_grid(func, lims::Tuple, xscale::Tuple, yscales::Tuple;
#     start_gridpoints=30, maxgridpoints=2000,
#     min_eps=(xscale[1](lims[2])-xscale[1](lims[1]))/10000,
#     reltol_dd=0.05, abstol_dd=1e-2)
function auto_grid(agc::AbstractAutoGridConfig)
    num_gridpoints = agc.start_gridpoints
    func = agc.func
    funcConfigs = agc.funcConfigs

    # Linearly spaced grid in xscale
    init_grid = collect(LinRange(agc.xscale(agc.lims[1]), agc.xscale(agc.lims[2]), agc.start_gridpoints))

    x1, x1_id, y1, y1_id = eval_at_point(agc, init_grid[1])

    values_id = tuple_of_vectors(y1_id, funcConfigs)
    grid_id = [x1_id,]

    xend, xend_id, yend, yend_id = eval_at_point(agc, init_grid[end])

    buffer_x = [xend,]
    buffer_x_id = [xend_id,]
    buffer_values = tuple_of_vectors(yend, funcConfigs)
    buffer_values_id = tuple_of_vectors(yend_id, funcConfigs)
    # Fill right buffer with 2,.....,(end-1) in reverse order
    # xend is already added
    for i in (length(init_grid)-1):-1:2
        xi, xi_id, yi, yi_id = eval_at_point(agc, init_grid[i])
        push_buffers!(xi, xi_id, yi, yi_id, buffer_x, buffer_x_id, buffer_values, buffer_values_id, funcConfigs)
    end

    x2, x2_id, y2, y2_id = pop_buffers!(buffer_x, buffer_x_id, buffer_values, buffer_values_id, funcConfigs)
    x3, x3_id, y3, y3_id = pop_buffers!(buffer_x, buffer_x_id, buffer_values, buffer_values_id, funcConfigs)

    # y1, y2, y3, unwrap_refine = potential_unwrap(y1, y2, y3, x1, x2, x3, agc)

    # Inside loop we assume x1,x2,x3 are set to current interval
    # x1, y1 should already be in `grid`, `values_id` etc., x2, x3 not
    unwrap_fail = false
    while true
        stop_refine = num_gridpoints >= agc.maxgridpoints
        stop_refine && @warn "Maximum number of gridpoints reached in refine_grid! at $x2_id, no further refinement will be made. Increase maxgridpoints to get better accuracy." maxlog=1

        y1, y2, y3, unwrap_refine = potential_unwrap(y1, y2, y3, x1, x2, x3, agc)
        deriv_refine = should_refine(y1, y2, y3, x1, x2, x3, agc)
        eps_too_small = (max(abs(x2 - x1),abs(x3 - x2)) < agc.min_eps)
        # If unwrapping unsuccessful and epsilon already too small
        if !unwrap_fail && unwrap_refine && eps_too_small
            # TODO change this to something better, based on derivatives
            @warn "Unwrapping unnuccessful, unwrapping ignored from now"
            unwrap_fail = true
        end
        if !stop_refine && !eps_too_small && (unwrap_refine || (deriv_refine && !unwrap_fail))
            refine_left = refine_right = true
            if 3*abs(x2 - x1) < abs(x3 - x2) # Don't refine a side if one interval is much smaller
                refine_left = false
            elseif abs(x2 - x1) > 3*abs(x3 - x2)
                refine_right = true
            end
            # Rightmost point
            push_buffers!(x3, x3_id, y3, y3_id, buffer_x, buffer_x_id, buffer_values, buffer_values_id, funcConfigs)
            if refine_right
                xm_right = (x2+x3)/2
                # Second to last point
                xi, xi_id, yi, yi_id = eval_at_point(agc, xm_right)
                if refine_left
                    # The new point is not used next iteration
                    push_buffers!(xi, xi_id, yi, yi_id, buffer_x, buffer_x_id, buffer_values, buffer_values_id, funcConfigs)
                else
                    # The new right point (possibly)
                    x3, x3_id, y3, y3_id = xi, xi_id, yi, yi_id
                end
                num_gridpoints += 1
            end
            if refine_left
                xm_left = (x1+x2)/2
                # Old middle point is now last point
                x3, x3_id, y3, y3_id = x2, x2_id, y2, y2_id
                # New middle point
                x2, x2_id, y2, y2_id = eval_at_point(agc, xm_left)
                num_gridpoints += 1
            end
        else
            # Save completed values to output
            push!(grid_id, x2_id)
            tuple_map(save, values_id, y2_id, y2, agc.funcConfigs)

            if isempty(buffer_x) # We are done
                push!(grid_id, x3_id)
                tuple_map(save, values_id, y3_id, y3, agc.funcConfigs)
                break
            end
            # Prepare for next loop
            x1, x1_id, y1, y1_id = x2, x2_id, y2, y2_id
            x2, x2_id, y2, y2_id = x3, x3_id, y3, y3_id
            x3, x3_id, y3, y3_id = pop_buffers!(buffer_x, buffer_x_id, buffer_values, buffer_values_id, funcConfigs)
            # Go to the next point
        end
    end
    return values_id, grid_id
end

function eval_at_point(agc, xi)
    xi_id = agc.xscale_inv(xi)
    yi_id = agc.func(xi_id)
    yi = apply_config_scales(agc.funcConfigs, yi_id)
    return xi, xi_id, yi, yi_id
end

function push_buffers!(xi, xi_id, yi, yi_id, buffer_x, buffer_x_id, buffer_values, buffer_values_id, funcConfigs)
    push!(buffer_x, xi) # Slightly unnessesary alloc
    push!(buffer_x_id, xi_id)
    tuple_map(save, buffer_values, yi, yi, funcConfigs)
    tuple_map(save, buffer_values_id, yi_id, yi_id, funcConfigs)
    return
end

function pop_buffers!(buffer_x, buffer_x_id, buffer_values, buffer_values_id, funcConfigs)
    xi = pop!(buffer_x)
    xi_id = pop!(buffer_x_id)
    yi = pop_tuple_of_vectors!(buffer_values, funcConfigs)
    yi_id = pop_tuple_of_vectors!(buffer_values_id, funcConfigs)
    return xi, xi_id, yi, yi_id
end

apply_scale(fc::AbstractFuncConfig, y_id) = fc.yscale(y_id)
#apply_scale(fc::SavingFuncConfig, y_id) = nothing

function apply_config_scales(funcConfigs, y_id)
    tuple_map(apply_scale, funcConfigs, y_id)
end

# TODO We can scale when saving instead of recomputing scaling
# Svectors should also be faster in general
function should_refine(v1, v2, v3, x1, x2, x3, agc::AbstractAutoGridConfig)
    not_linear = !all(tuple_map(is_almost_linear, v1, v2, v3, x1, x2, x3, agc.funcConfigs, agc))
    if not_linear
        return not_linear
    end
    is_peak = any(tuple_map(potential_peak, v1, v2, v3, agc.funcConfigs, agc))
    return is_peak
end

function potential_peak(y1, y2, y3, fc::AbstractFuncConfig, agc::AbstractAutoGridConfig)
    return fc.find_peaks && dotdiff(y2,y1,y3,y2) < 0
end

@inline dotdiff(v1::T,v2::T,v3::T,v4::T) where T<:Real = (v1-v2)*(v3-v4)
# TODO UNSAFE
@inline function dotdiff(v1::T,v2::T,v3::T,v4::T) where T
    (length(v1) == length(v2) == length(v3) == length(v4)) || throw(DimensionMismatch())
    s = zero(dot(first(v1)-first(v2), first(v3)-first(v4)))
    for i in eachindex(v1)
        @inbounds s += dot(v1[i]-v2[i], v3[i]-v4[i])
    end
    return s
end

function is_almost_linear(y1, y2, y3, x1, x2, x3, fc::AbstractFuncConfig, agc::AbstractAutoGridConfig)
    # We assume that x2-x1 \approx  x3-x2, so we need to check that y2-y1 approx y3-y2
    # Essentially low second derivative compared to derivative, so that linear approximation is good.
    # Second argument to avoid too small steps when derivatives are small
    dd = derivderiv2(y1,y2,y3,x1,x2,x3)
    #rhs = max((agc.reltol_dd)^2*max(deriv2(y2,y1,x2,x1), deriv2(y3,y2,x3,x2)), (agc.abstol_dd)^2)
    d = max(deriv2(y2,y1,x2,x1), deriv2(y3,y2,x3,x2))
    return dd <  agc.reltol_dd*max(d, agc.abstol_dd)/(x3-x1)^2
end
is_almost_linear(y1, ym, y2, x1, x2, x3, fc::SavingFuncConfig, agc::AbstractAutoGridConfig) = true

function potential_unwrap(v1, v2, v3, x1, x2, x3, agc::AbstractAutoGridConfig)
    vnew = tuple_map(potential_unwrap, v1, v2, v3, x1, x2, x3, agc.funcConfigs, agc)
    v1_new, v2_new, v3_new, do_refine =  zip(vnew...)
    return v1_new, v2_new, v3_new, any(do_refine)
end

potential_unwrap(y1, y2, y3, x1, x2, x3, fc::AbstractFuncConfig, agc::AbstractAutoGridConfig) = y1, y2, y3, false

function potential_unwrap(y1, y2, y3, x1, x2, x3, fc::WrappingFuncConfig, agc::AbstractAutoGridConfig)
    mid = +(fc.wraplims...)/2
    width = -(-(fc.wraplims...))
    y2new = y2 - fld(y2-y1 + mid, width) * width
    y3new = y3 - fld(y3-y2new + mid, width) * width
    do_refine = (sqrt(normdiff2(y2new, y1)) > width/2) || (sqrt(normdiff2(y3new, y2new)) > width/2)
    return y1, y2new, y3new, do_refine
end

deriv2(y1,y2,x1,x2) = normdiff2(y1,y2)/abs2(x1-x2)
# This will probably never (almost) be called
@inline normdiff2(y1::Number,y2::Number) = abs2(y1-y2)

function normdiff2(y1::A,y2::A) where {T<:Number, N, A <: Union{NTuple{N,T},AbstractArray{T}}}
    length(y1) == length(y2) || throw(DimensionMismatch())
    sum = abs2(y1[1]-y2[1])
    @inbounds for i in 2:length(y1)
        sum += abs2(y1[i]-y2[i])
    end
    return sum
end
function normdiff2(x,y)
    length(x) == length(y) || throw(DimensionMismatch())
    sum = normdiff2(x[1],y[1])
    @inbounds for i in 2:length(x)
        sum += normdiff2(x[i],y[i])
    end
    return sum
end
# @inline normdiff(x,y) = sqrt(normdiff2(x,y))

function derivderiv2(y1,y2,y3,x1,x2,x3)
    d1 = abs(x1-x2)
    d2 = abs(x2-x3)
    diff = normdiffdiff2(y1,y2,y3,d1,d2)
    dsquare = (d1*d2*(d1+d2))^2/4
    diff/dsquare
end

# This will probably never(almost) be called
@inline normdiffdiff2(x::Number,y::Number,z::Number,d1,d2) = abs2(d2*x+d1*z-(d1+d2)*y)

function normdiffdiff2(x::A,y::A,z::A,d1,d2) where {T<:Number, N, A <: Union{NTuple{N,T},AbstractArray{T}}}
    length(x) == length(y) || throw(DimensionMismatch())
    sum = abs2(d2*x[1]+d1*z[1]-(d1+d2)*y[1])
    @inbounds for i in 2:length(x)
        sum += abs2(d2*x[i]+d1*z[i]-(d1+d2)*y[i])
    end
    return sum
end
function normdiffdiff2(x,y,z,d1,d2)
    length(x) == length(y) || throw(DimensionMismatch())
    sum = normdiffdiff2(x[1],y[1],z[1],d1,d2)
    @inbounds for i in 2:length(x)
        sum += normdiffdiff2(x[i],y[i],z[i],d1,d2)
    end
    return sum
end
# @inline normdiffdiff(x,y,z) = sqrt(normdiffdiff2(x,y,z))

# @generated function apply_tuple3(fs::Tuple{Vararg{<:Any,N}}, x::Tuple{Vararg{<:Any,N}}) where N
#     vec = []
#     for i in 1:N
#         push!(vec, :(fs[$i].(x[$i])))
#     end
#     :(Core.tuple($(vec...)))
# end

@generated function tuple_of_vectors(yi_id::Tuple{Vararg{<:Any,N}}, funcConfigs::Tuple{Vararg{<:AbstractFuncConfig,N}}) where {N}
    vec = []
    for i in 1:N
        # if funcConfigs.types[i] <: AbstractFuncConfig{true}
            push!(vec, :([yi_id[$i],]))
        # else
            # push!(vec, :(nothing))
        # end
    end
    :(Core.tuple($(vec...),))
end

@generated function pop_tuple_of_vectors!(vals::Tuple{Vararg{<:Any,N}}, funcConfigs::Tuple{Vararg{<:AbstractFuncConfig,N}}) where {N}
    vec = []
    for i in 1:N
        # if funcConfigs.types[i] <: AbstractFuncConfig{true}
            push!(vec, :(pop!(vals[$i])))
        # else
            # push!(vec, :(nothing))
        # end
    end
    :(Core.tuple($(vec...),))
end

#save(values, ym_id, funcConfig::AbstractFuncConfig{false}) = nothing
function save(values,  ym_id, ym, funcConfig::AbstractFuncConfig)
    push!(values, ym_id)
    return nothing
end

# We force Wrapping to be identity mapping. The unwrapped values are in ym
function save(values,  ym_id, ym, funcConfig::WrappingFuncConfig)
    push!(values, ym)
    return nothing
end

@generated function tuple_map(func, args...)
    vec = []
    N = 0
    for arg in args
        if arg <: Tuple
            N = length(arg.types)
        end
    end
    for i in 1:N
        ex = []
        for j in 1:length(args)
            if args[j] <: Tuple
                push!(ex, :(args[$j][$i]))
            else
                push!(ex, :(args[$j]))
            end
        end
        push!(vec, :(func($(ex...))))
    end
    :(Core.tuple($(vec...),))
end

end


text:
export truncated_svd, TSVD

struct TSVD end

truncated_svd(A::AbstractMatrix, alg = SVDFact(); tol = sqrt(eps(eltype(A))), rmin = 1, rmax = minimum(size(A))) = truncated_svd(A, alg, tol, rmin, rmax)
function truncated_svd(A::AbstractMatrix, rank::Int, alg = SVDFact())
    @assert rank <= minimum(size(A)) "rank ≤ min(n,m)"
    truncated_svd(A, alg; rmin = rank, rmax = rank)
end

function truncated_svd(A::AbstractMatrix, ::SVDFact, tol, rmin::Int, rmax::Int)
    U, S, V = svd(A)
    r = (rmax == rmin) ? rmax : max(min(truncate_to_tolerance(S, tol), rmax), rmin)
    return SVDLikeRepresentation(U[:,1:r], diagm(S[1:r]), V[:,1:r])
end

function truncated_svd(A::AbstractMatrix, ::TSVD, tol, rmin::Int, rmax::Int)
    U, S, V = tsvd(A, rmax)
    return SVDLikeRepresentation(U, diagm(S), V)
end

function truncate_to_tolerance(S, tol)
    s = 0
    r = length(S)
    for σ in reverse(S)
        s += σ^2
        if s > tol^2
            break
        end
        r -= 1
    end
    return r
end

function blockdiagonal(A::T1, B::T2) where {T1, T2 <: Union{AbstractMatrix, AbstractVector}}
    n1,m1 = size(A,1), size(A,2)
    n2,m2 = size(B,1), size(B,2)
    C = zeros(eltype(A), n1+n2, m1+m2)
    C[1:n1, 1:m1] .= A
    C[n1+1:end, m1+1:end] .= B
    return C
end


text:
@testset "Partitioning" begin
  setify(lists) = Set(Set.(lists))

  d = RegularGrid(10,10)
  p = partition(d, UniformPartitioner(100))
  @test sprint(show, p) == "100 SpatialPartition"
  @test sprint(show, MIME"text/plain"(), p) == "100 SpatialPartition\n  N° points\n  └─1\n  └─1\n  └─1\n  └─1\n  └─1\n  ⋮\n  └─1\n  └─1\n  └─1\n  └─1\n  └─1"

  @testset "UniformPartitioner" begin
    grid = RegularGrid(3,3)

    Random.seed!(123)
    p = partition(grid, UniformPartitioner(3, false))
    @test setify(subsets(p)) == setify([[1,2,3], [4,5,6], [7,8,9]])
    p = partition(grid, UniformPartitioner(3))
    @test setify(subsets(p)) == setify([[8,6,9], [4,1,7], [2,3,5]])

    grid = RegularGrid(2,3)
    p = partition(grid, UniformPartitioner(3, false))
    @test setify(subsets(p)) == setify([[1,2], [3,4], [5,6]])
  end

  @testset "DirectionPartitioner" begin
    grid = RegularGrid(3,3)

    # basic checks on small regular grid data
    p = partition(grid, DirectionPartitioner((1.,0.)))
    @test setify(subsets(p)) == setify([[1,2,3], [4,5,6], [7,8,9]])

    p = partition(grid, DirectionPartitioner((0.,1.)))
    @test setify(subsets(p)) == setify([[1,4,7], [2,5,8], [3,6,9]])

    p = partition(grid, DirectionPartitioner((1.,1.)))
    @test setify(subsets(p)) == setify([[1,5,9], [2,6], [3], [4,8], [7]])

    p = partition(grid, DirectionPartitioner((1.,-1.)))
    @test setify(subsets(p)) == setify([[1], [2,4], [3,5,7], [6,8], [9]])

    # opposite directions produce same partition
    dir1 = (rand(), rand()); dir2 = .-dir1
    p1 = partition(grid, DirectionPartitioner(dir1))
    p2 = partition(grid, DirectionPartitioner(dir2))
    @test setify(subsets(p1)) == setify(subsets(p2))

    # partition of arbitrarily large regular grid always
    # returns the "lines" and "columns" of the grid
    for n in [10,100,200]
      grid = RegularGrid(n,n)

      p = partition(grid, DirectionPartitioner((1.,0.)))
      @test setify(subsets(p)) == setify([collect((i-1)*n+1:i*n) for i in 1:n])
      ns = [nelms(d) for d in p]
      @test all(ns .== n)

      p = partition(grid, DirectionPartitioner((0.,1.)))
      @test setify(subsets(p)) == setify([collect(i:n:n*n) for i in 1:n])
      ns = [nelms(d) for d in p]
      @test all(ns .== n)
    end
  end

  @testset "FractionPartitioner" begin
    grid = RegularGrid(10,10)

    p = partition(grid, FractionPartitioner(0.5))
    @test nelms(p[1]) == nelms(p[2]) == 50
    @test length(p) == 2

    p = partition(grid, FractionPartitioner(0.7))
    @test nelms(p[1]) == 70
    @test nelms(p[2]) == 30

    p = partition(grid, FractionPartitioner(0.3))
    @test nelms(p[1]) == 30
    @test nelms(p[2]) == 70
  end

  @testset "SLICPartitioner" begin
    𝒮 = georef((z=[ones(10,10) 2ones(10,10); 3ones(10,10) 4ones(10,10)],))
    p = partition(𝒮, SLICPartitioner(4, 1.0))
    @test length(p) == 4
    @test all(nelms.(p) .== 100)
    @test mean(coordinates(p[1]), dims=2) == [ 4.5, 4.5][:,:]
    @test mean(coordinates(p[2]), dims=2) == [14.5, 4.5][:,:]
    @test mean(coordinates(p[3]), dims=2) == [ 4.5,14.5][:,:]
    @test mean(coordinates(p[4]), dims=2) == [14.5,14.5][:,:]

    𝒮 = georef((z=[√(i^2+j^2) for i in 1:100, j in 1:100],))
    p = partition(𝒮, SLICPartitioner(50, 1.0))
    @test length(p) == 49

    if visualtests
      @plottest plot(p) joinpath(datadir,"slic.png") !istravis
    end
  end

  @testset "BlockPartitioner" begin
    grid = RegularGrid(10,10)

    p = partition(grid, BlockPartitioner(5.,5.))
    @test length(p) == 4
    @test all(nelms.(p) .== 25)

    p = partition(grid, BlockPartitioner(5.,2.))
    @test length(p) == 12
    @test Set(nelms.(p)) == Set([5,10])
  end

  @testset "BisectPointPartitioner" begin
    grid = RegularGrid(10,10)

    p = partition(grid, BisectPointPartitioner((0.,1.), (5.,5.1)))
    @test nelms(p[1]) == 60
    @test nelms(p[2]) == 40

    # all points in X₁ are below those in X₂
    X₁ = coordinates(p[1])
    X₂ = coordinates(p[2])
    M₁ = maximum(X₁, dims=2)
    m₂ = minimum(X₂, dims=2)
    @test all(X₁[2,j] < m₂[2] for j in 1:size(X₁,2))
    @test all(X₂[2,j] > M₁[2] for j in 1:size(X₂,2))

    # flipping normal direction is equivalent to swapping subsets
    p₁ = partition(grid, BisectPointPartitioner(( 1.,0.), (5.1,5.)))
    p₂ = partition(grid, BisectPointPartitioner((-1.,0.), (5.1,5.)))
    @test nelms(p₁[1]) == nelms(p₂[2]) == 60
    @test nelms(p₁[2]) == nelms(p₂[1]) == 40
  end

  @testset "BisectFractionPartitioner" begin
    grid = RegularGrid(10,10)

    p = partition(grid, BisectFractionPartitioner((1.,0.), 0.2))
    @test nelms(p[1]) == 20
    @test nelms(p[2]) == 80

    # all points in X₁ are to the left of X₂
    X₁ = coordinates(p[1])
    X₂ = coordinates(p[2])
    M₁ = maximum(X₁, dims=2)
    m₂ = minimum(X₂, dims=2)
    @test all(X₁[1,j] < m₂[1] for j in 1:size(X₁,2))
    @test all(X₂[1,j] > M₁[1] for j in 1:size(X₂,2))

    # flipping normal direction is equivalent to swapping subsets
    p₁ = partition(grid, BisectFractionPartitioner(( 1.,0.), 0.2))
    p₂ = partition(grid, BisectFractionPartitioner((-1.,0.), 0.8))
    @test nelms(p₁[1]) == nelms(p₂[2]) == 20
    @test nelms(p₁[2]) == nelms(p₂[1]) == 80
  end

  @testset "BallPartitioner" begin
    pset = PointSet([
      0 1 1 0 0.2
      0 0 1 1 0.2
    ])

    # 3 balls with 1 point, and 1 ball with 2 points
    p = partition(pset, BallPartitioner(0.5))
    n = nelms.(p)
    @test length(p) == 4
    @test count(i->i==1, n) == 3
    @test count(i->i==2, n) == 1
    @test setify(subsets(p)) == setify([[1,5],[2],[3],[4]])

    # 5 balls with 1 point each
    p = partition(pset, BallPartitioner(0.2))
    @test length(p) == 5
    @test all(nelms.(p) .== 1)
    @test setify(subsets(p)) == setify([[1],[2],[3],[4],[5]])
  end

  @testset "PlanePartitioner" begin
    grid = RegularGrid(3,3)
    p = partition(grid, PlanePartitioner((0.,1.)))
    @test setify(subsets(p)) == setify([[1,2,3],[4,5,6],[7,8,9]])

    grid = RegularGrid(4,4)
    p = partition(grid, PlanePartitioner((0.,1.)))
    @test setify(subsets(p)) == setify([1:4,5:8,9:12,13:16])
  end

  @testset "VariablePartitioner" begin
    z = vec([1 1 1; 2 2 2; 3 3 3])
    sdata = georef(DataFrame(z=z), RegularGrid(3,3))
    p = partition(sdata, VariablePartitioner(:z))
    @test setify(subsets(p)) == setify([[1,4,7],[2,5,8],[3,6,9]])

    # partition with missing values
    z = vec([missing 1 1; 2 missing 2; 3 3 missing])
    sdata = georef(DataFrame(z=z), RegularGrid(3,3))
    p = partition(sdata, VariablePartitioner(:z))
    @test setify(subsets(p)) == setify([[4,7],[2,8],[3,6],[1,5,9]])
  end

  @testset "PredicatePartitioner" begin
    grid = RegularGrid(3,3)

    # partition even from odd locations
    pred(i,j) = iseven(i+j)
    p = partition(grid, PredicatePartitioner(pred))
    @test setify(subsets(p)) == setify([1:2:9,2:2:8])
  end

  @testset "SpatialPredicatePartitioner" begin
    g = RegularGrid(10,10)

    # check if there are 100 partitions, each one having only 1 point
    sp = SpatialPredicatePartitioner((x,y) -> norm(x-y) < 1.0)
    s = subsets(partition(g, sp))
    @test length(s) == 100
    for d in partition(g, sp)
      @test nelms(d) == 1
    end
    # defining a predicate to check if points x and y belong to the square [0.,5.]x[0.,5.]
    pred(x, y) = all([0.,0.] .<= x .<=[5.,5.]) && all([0.,0.] .<= y .<= [5.,5.])
    sp = SpatialPredicatePartitioner(pred)
    p = partition(g, sp)
    s = subsets(p)
    n = nelms.(p)

    # There will be 65 partitions:
    # 1 partition with 36 points (inside square [0.,5.]x[0.,5.])
    # 64 partitions with only 1 point inside each of them
    @test length(s) == 65
    @test maximum(length.(s)) == 36
    @test count(i->i==1, n) == 64
    @test count(i->i==36, n) == 1
  end

  @testset "ProductPartitioner" begin
    g = RegularGrid(100,100)
    bm = BlockPartitioner(10.,10.)
    bn = BlockPartitioner(5.,5.)

    # Bm x Bn = Bn with m > n
    s1 = subsets(partition(g, bm*bn))
    s2 = subsets(partition(g, bn))
    @test setify(s1) == setify(s2)

    # pXp=p (for deterministic p)
    for p in [BlockPartitioner(10.,10.),
              BisectFractionPartitioner((0.1,0.1))]
      s1 = subsets(partition(g, p*p))
      s2 = subsets(partition(g, p))
      @test setify(s1) == setify(s2)
    end
  end

  @testset "HierarchicalPartitioner" begin
    g = RegularGrid(100,100)
    bm = BlockPartitioner(10.,10.)
    bn = BlockPartitioner(5.,5.)

    # Bn -> Bm = Bm with m > n
    s1 = subsets(partition(g, bm → bn))
    s2 = subsets(partition(g, bn))
    @test setify(s1) == setify(s2)
  end

  @testset "Mixed Tests" begin
    g = RegularGrid(100,100)
    bm = BlockPartitioner(10.,10.)
    bn = BlockPartitioner(5.,5.)

    # Bm*Bn = Bm->Bn
    s1 = subsets(partition(g, bm * bn))
    s2 = subsets(partition(g, bm → bn))
    @test setify(s1) == setify(s2)
  end
end


text:
using ChainRulesCore
using ChainRulesTestUtils
using LinearAlgebra
using Random
using Test

@testset "ChainRulesTestUtils.jl" begin
    include("generate_tangent.jl")
    include("to_vec.jl")
    include("isapprox.jl")
    include("iterator.jl")
    include("testers.jl")
    include("data_generation.jl")
end


text:
seating = open("seating.txt") do f
    collect.(readlines(f))
end

function print_map(arr)
    for r in arr
        println(join(r))
    end
    println()
end

function count_neighborhood(arr, x::Int64, y::Int64, c::Char)::Int64
    n = 0

    if y > 1
        if arr[y-1][x] == c
            n += 1
        end
        if x > 1 && arr[y-1][x-1] == c
            n += 1
        end
        if x < length(arr[1]) && arr[y-1][x+1] == c
            n += 1
        end
    end

    if y < length(arr)
        if arr[y+1][x] == c
            n += 1
        end
        if x > 1 && arr[y+1][x-1] == c
            n += 1
        end
        if x < length(arr[1]) && arr[y+1][x+1] == c
            n += 1
        end
    end

    if x > 1 && arr[y][x-1] == c
        n += 1
    end

    if x < length(arr[1]) && arr[y][x+1] == c
        n += 1
    end

    return n
end

function count_visible(arr, x::Int64, y::Int64, c::Char)::Int64
    n = 0
    unit_vectors = map(x -> [round(Int,sin(x)),round(Int,cos(x))], pi*(-1:1/4:0.9))
    for uv in unit_vectors
        l = 1
        while true
            ox, oy = l * uv
            if y+oy < 1 || y+oy > length(arr) || x+ox < 1 || x+ox > length(arr[1])
                break
            end
            if arr[y+oy][x+ox] == c
                n += 1
            end
            if arr[y+oy][x+ox] != '.'
                break
            end
            l += 1
        end
    end
    return n
end

function simulate(arr, fn, n)
    step = deepcopy(arr)

    for y in 1:length(arr)
        for x in 1:length(arr[1])
            # print("at $(y),$(x): $(arr[y][x]): ")
            if arr[y][x] == '.'
                # println("do nothing")
            elseif arr[y][x] == 'L'
                occupied = fn(arr, x, y, '#')
                if occupied == 0
                    step[y][x] = '#'
                end
            elseif arr[y][x] == '#'
                occupied = fn(arr, x, y, '#')
                if occupied >= n
                    step[y][x] = 'L'
                end
            end
        end
    end

    return step
end

function run(state, fn, n)
    i = 0
    while true
        i += 1
        next_state = simulate(state, fn, n)
        if state == next_state
            break
        end
        state = deepcopy(next_state)
    end

    return state
end

println(sum(map(r -> count(x -> x == '#', r), run(seating, count_neighborhood, 4))))
println(sum(map(r -> count(x -> x == '#', r), run(seating, count_visible, 5))))


text:
using Plots
@recipe plot(sd::SimData...; vars=nothing, names=nothing) = begin
    if vars === nothing
        error("Must specify variables to plot")
    elseif vars == :all
        vars = StateSpaceEcon._names(sd[1])
    end
    if length(vars) > 10
        error("Too many variables. Split into pages.")
    end
    if names === nothing
        names = ["data$i" for i = 1:length(sd)]
    end
    layout --> length(vars)
    title --> reshape(map(string, [vars...]), 1, :)
    titlefont --> font(10, :bold)
    label --> repeat(reshape([names...], 1, :), inner=(1, length(vars)))
    linewidth --> 1.5
    left_margin --> 4 * Plots.mm
    for s in sd
        for v in vars
            @series begin
                s[v]
            end
        end
    end
end


text:
include("utils.jl")

build_so("005")

immutable Foo
   i::Cint
end

x = 13

y = ccall(
    (:bar, "binary/lib005.so"),
    Foo,
    (Cint, ),
    x
)

@assert y.i == 13


text:
# Remove the reads with the length of sequence which is different from the length of quality
function remove_bad_read!(records::Array{FASTX.FASTQ.Record, 1})
	seq_lengths = FASTQ.seqlen.(records)
	qual_lengths = map(i -> length(FASTQ.quality(i)), records)
	return records[seq_lengths .==  qual_lengths]
end

# Methods for parsing H5Record

function parse_f5read_record end

function parse_f5read_record(record::HDF5Group, basecallGroup::String)::FASTX.FASTQ.Record
	fastqpath = string("Analyses/$basecallGroup/BaseCalled_template/Fastq");
	fastq = FASTQ.Record(read(record, fastqpath))
end

function parse_f5read_record(record::HDF5File, basecallGroup::String)::FASTX.FASTQ.Record
	fastqpath = string("Analyses/$basecallGroup/BaseCalled_template/Fastq");
	fastq = FASTQ.Record(read(record, fastqpath))
end

# Methods for Record-parsing
function get_info end

# Parse HDF5.Group
function get_info(recordAttrs::HDF5.Attributes, record::HDF5Group, basecallGroup::String)
	readqual = read(recordAttrs, "mean_qscore")
	readlen = read(recordAttrs, "sequence_length")
	readqc = calculate_GC_content(record, basecallGroup)
	FastqInfo(readqual, readlen, readqc)
end

# Parse HDF5.File
function get_info(recordAttrs::HDF5.Attributes, record::HDF5File, h5ver::Float64, basecallGroup::String)
	readqual = read(recordAttrs, "mean_qscore")
	readlen = read(recordAttrs, "sequence_length")
	readqc = calculate_GC_content(record, h5ver, basecallGroup)
	FastqInfo(readqual, readlen, readqc)
end

# Parse BAM.Record
function get_info(record::XAM.BAM.Record)
	flag_check = BAM.flag(record) |> Int
	qual_check = BAM.quality(record) |> length
	if flag_check == 4
		return nothing
	end

	if qual_check == 0
		return nothing
	end

	nm = record["NM"] |> Int
	matches, gaps, cgaps = BAM.cigar_rle(record) |> parse_cigar
	readident = calculate_identity(nm, matches, gaps, cgaps)
	readlen = BAM.seqlength(record)
	readqual = BAM.quality(record) |> calculate_mean_readqual
	readgc = calculate_GC_content(record)
	return BAMInfo(readqual, readlen, readident, readgc)
end

# Parse FASTQ.Record
function get_info(record::FASTX.FASTQ.Record)
	readqual = FASTQ.quality(record) |> calculate_mean_readqual
	readlen = FASTQ.seqlen(record)
	readgc = calculate_GC_content(record)
	FastqInfo(readqual, readlen, readgc)
end


text:
using Test
using Quiqbox
using Quiqbox: BasisFuncMix, hasBoolRelation
using Suppressor: @capture_out

@testset "Overload.jl" begin

    # function show
    pb1 = ParamBox(-1, canDiff=false)
    @test (@capture_out show(pb1)) == string(typeof(pb1))*"(-1.0)[∂][undef]"

    pb2 = ParamBox(-1, :a, index=1)
    @test (@capture_out show(pb2)) == string(typeof(pb2))*"(-1.0)[∂][a₁]"

    pb3 = ParamBox(-1, :x, abs)
    @test (@capture_out show(pb3)) == string(typeof(pb3))*"(-1.0)[∂][x_x]"

    bf1 = genBasisFunc([1,2,1], (2,1))
    gf1 = bf1.gauss[1]
    @test (@capture_out show(gf1)) == string(typeof(gf1))*"(xpn="*
                                      string(typeof(gf1.param[1]))*"(2.0)[∂][α], con="*
                                      string(typeof(gf1.param[2]))*"(1.0)[∂][d])"
    @test (@capture_out show(bf1)) == string(typeof(bf1))*"(gauss, subshell, center)"*
                                      "[X⁰Y⁰Z⁰][1.0, 2.0, 1.0]"

    bf2 = genBasisFunc(missing, "STO-3G")[]
    @test (@capture_out show(bf2)) == string(typeof(bf2))*"(gauss, subshell, center)"*
                                      "[X⁰Y⁰Z⁰]"*"[NaN, NaN, NaN]"

    bfs1 = genBasisFunc([0,0,0], (2,1), "P")
    @test (@capture_out show(bfs1)) == string(typeof(bfs1))*"(gauss, subshell, center)"*
                                       "[3/3]"*"[0.0, 0.0, 0.0]"

    bfs2 = genBasisFunc([0,0,0], (2,1), [[2,0,0]])
    @test (@capture_out show(bfs2)) == string(typeof(bfs2))*"(gauss, subshell, center)"*
                                       "[X²Y⁰Z⁰]"*"[0.0, 0.0, 0.0]"

    bfm1 = BasisFuncMix([bf1, bf2])
    @test (@capture_out show(bfm1)) == string(typeof(bfm1))*"(BasisFunc, param)"

    GTb1 = GTBasis([bf1, bfs2])
    @test (@capture_out show(GTb1)) == string(typeof(GTb1))*"(basis, S, Te, eeI, getVne, "*
                                       "getHcore)"

    box1 = GridBox(2, 1.5)
    @test (@capture_out show(box1)) == string(typeof(box1))*"(num, len, coord)"

    fVar1 = runHF(GTb1, ["H", "H"], [[0, 0, 0], [1,2,1]], printInfo=false, initialC=:Hcore)

    @test (@capture_out show(fVar1.temp)) == string(typeof(fVar1.temp))*"(shared.Etots="*
                                             "[2.263712269, … , 2.262890932], shared.Dtots"*
                                             ", Cs, Es, Ds, Fs)"

    @test (@capture_out show(fVar1)) == string(typeof(fVar1))*"(E0HF=2.262890932, C, F, D,"*
                                        " Emo, occu, temp, isConverged)"

    # function hasBoolRelation
    pb4 = deepcopy(pb1)
    @test true  == hasBoolRelation(==, pb1, pb4)
    @test false == hasBoolRelation(===, pb1, pb4)
    @test false == hasBoolRelation(==, pb2, pb3)
    @test true  == hasBoolRelation(==, pb2, pb3, ignoreFunction = true)
    @test false == hasBoolRelation(==, pb1, pb2)
    @test true  == hasBoolRelation(==, pb1, pb2, ignoreContainer=true)

    # function +
    testAdd = function (a1, a2)
        hasEqual(a2 +  a1 + a1,
                 a1 +  a2 + a1,
                 a1 +  a1 + a2,
                 a1 + (a1 + a2),
                 add(add(a1, a2), a1),
                 add(a1, add(a2, a1)))
    end

    bf3 = genBasisFunc([1,2,1], ([2,1], [0.2, 0.4]))
    bf4 = genBasisFunc([1,1,1], (3,0.5))
    bfm2 = BasisFuncMix([bf1, bf4, bf1])
    bfm3 = BasisFuncMix([bf4, bf4, bf3])

    @test hasEqual(bf1 + bf3 + bf4,
                   bf1 + bf4 + bf3,
                   bf3 + bf1 + bf4,
                   bf3 + bf4 + bf1,
                   bf4 + bf3 + bf1)
    @test testAdd(bf1, bf3)
    @test testAdd(bf1, bf4)
    @test testAdd(bf3, bf4)
    @test testAdd(bf1, bfm2)
    @test testAdd(bf3, bfm2)
    @test testAdd(bf4, bfm2)
    @test testAdd(bf1, bfm3)
    @test testAdd(bf3, bfm3)
    @test testAdd(bf4, bfm3)
    @test testAdd(bfm2, bfm3)

    # function *
    @test hasEqual(gf1 * GaussFunc(0.2, 1.5), GaussFunc(2.2, 1.5))
    @test hasEqual(gf1 * π, GaussFunc(2, π))
    boolFunc(x, y) = (==)(x,y)
    boolFunc(x::Array{<:Real, 0}, y::Array{<:Real, 0}) = isapprox(x, y, atol=1e-12)
    boolFunc(x::Real, y::Real) = isapprox(x, y, atol=1e-12)
    testMul = function (a1, a2, ignoreContainer=false)
        r1 = hasBoolRelation(boolFunc,
                             a1*(a1 + a2),
                             (a1 + a2)*a1,
                             a1*a1 + a1*a2,
                             a1*a1 + a2*a1,
                             mul(a1, add(a1, a2)),
                             mul(add(a1, a2), a1); ignoreContainer)

        r2 = hasBoolRelation(boolFunc, a1 *  a2 * a1,
                             a1 * (a2 * a1),
                             mul(mul(a1, a2), a1),
                             mul(a1, mul(a2, a1)); ignoreContainer)
        r1 * r2
    end

    bf5 = genBasisFunc([1,1,1], GaussFunc(genExponent(3), genContraction(0.2, x->5x)))
    bfm4 = BasisFuncMix([bf4, bf5, bf4])

    @test testMul(bf1,  bf3)
    @test testMul(bf1,  bf4)
    @test testMul(bf1,  bf5)
    @test testMul(bf3,  bf4)
    @test testMul(bf3,  bf5)
    @test testMul(bf4,  bf5)
    @test testMul(bf1,  bfm2)
    @test testMul(bf3,  bfm2)
    @test testMul(bf4,  bfm2)
    @test testMul(bf1,  bfm3)
    @test testMul(bf3,  bfm3)
    @test testMul(bf4,  bfm3)
    @test testMul(bf1,  bfm4)
    @test testMul(bf3,  bfm4)
    @test testMul(bf4,  bfm4)
    @test testMul(bfm2, bfm3)
    @test testMul(bfm2, bfm4)
    @test testMul(bfm3, bfm4)

    testMul2 = function (a, a1, a2)
        r1 = hasBoolRelation(boolFunc,
                             a*(a1 + a2),
                             (a1 + a2)*a,
                             a*a1 + a*a2,
                             a1*a + a2*a, ignoreContainer=true)

        r2 = hasBoolRelation(boolFunc,
                             a  * a1 * a2,
                             a  * a2 * a1,
                             a1 * a  * a2,
                             a1 * a2 * a ,
                             a2 * a1 * a ,
                             a2 * a  * a1,
                             a  * (a1* a2),
                             a1 * (a2* a ),
                             a2 * (a * a1), ignoreContainer=true)
        r1 * r2
    end

    c = 2.0

    @test testMul2(c, bf1 , bf3 )
    @test testMul2(c, bf1 , bf4 )
    @test testMul2(c, bf1 , bf5 )
    @test testMul2(c, bf3 , bf4 )
    @test testMul2(c, bf3 , bf5 )
    @test testMul2(c, bf4 , bf5 )
    @test testMul2(c, bf1 , bfm2)
    @test testMul2(c, bf3 , bfm2)
    @test testMul2(c, bf4 , bfm2)
    @test testMul2(c, bf1 , bfm3)
    @test testMul2(c, bf3 , bfm3)
    @test testMul2(c, bf4 , bfm3)
    @test testMul2(c, bf1 , bfm4)
    @test testMul2(c, bf3 , bfm4)
    @test testMul2(c, bf4 , bfm4)
    @test testMul2(c, bfm2, bfm3)
    @test testMul2(c, bfm2, bfm4)
    @test testMul2(c, bfm3, bfm4)

    c1 = 1.5
    c2 = 2.5

    @test testMul2.([bf1, bf3, bf4, bf5,
                     bfm2.BasisFunc |> sum, bfm3.BasisFunc |> sum, bfm4.BasisFunc |> sum],
                     c1, c2) |> prod

    @test .!testMul2.([bfm2, bfm3, bfm4], c1, c2) |> prod

    # function iterate, size, length, ndims
    @test iterate(pb1) == (-1, nothing)
    @test iterate(gf1) == (gf1, nothing)
    @test iterate(bf1) == (bf1, nothing)

    cs = [pb1, gf1, bf1]
    @test (iterate.(cs, 1) .=== nothing) |> prod
    @test (size.(cs) .== Ref(())) |> prod
    @test (size.(cs, 1) .== 1) |> prod
    @test (length.(cs) .== 1) |> prod
    @test ndims(pb1) == 0

    @test iterate(bfm1) == (bfm1, nothing)
    @test iterate(bfm1, rand()) === nothing
    @test hasEqual(iterate(bfs1), (bfs1[1], 2))
    @test hasEqual(iterate(bfs1, 1), (bfs1[1], 2))
    @test hasEqual(iterate(bfs1, length(bfs1.ijk)+1), nothing)

    @test size(bfm1) == ()
    @test size(bfm1, 1) == 1
    @test size(bfs1) == (length(bfs1.ijk),)
    @test size(bfs1, 1) == length(bfs1.ijk)

    @test length(bfm1) == 1
    @test length(bfs1) == length(bfs1.ijk)

    # function getindex, setindex!, firstindex, lastindex, axes
    @test getindex(pb1) == pb1[] == pb1[begin] == pb1[end] == -1
    @test (pb1[] = -2; res = (pb1[] == -2); pb1[] = -1; res)
    @test axes(pb1) == ()

    @test getindex(gf1) == gf1[] == gf1[begin] == gf1[end] == (gf1.param |> collect)
    @test getindex(bf1) == bf1[] == bf1[begin] == bf1[end] == (bf1.gauss |> collect)
    bfm11 = Quiqbox.BasisFuncMix([bf1, bf2, bf1])
    gs = getindex(bfm11)
    @test gs == bfm11[] == bfm11[begin] == bfm11[end]
    @test getUnique!(gs) == [bf1.gauss[1], bf2.gauss...]
    bfs1_alter = genBasisFunc.(Ref([0,0,0]), Ref((2,1)), [[1,0,0], [0,1,0], [0,0,1]])
    for i in eachindex(bfs1)
        @test hasEqual(getindex(bfs1, i), bfs1[i], bfs1_alter[i])
    end
    @test hasEqual(bfs1[begin], bfs1[1])
    @test hasEqual(bfs1[end], bfs1[3])
    @test getindex(bfs1) == fill(bfs1.gauss[1], 3)

    # function broadcastable
    @test getfield.(pb1, [:data, :map]) == [pb1.data, pb1.map]
    @test getfield.(gf1, [:xpn, :con]) == [gf1.xpn, gf1.con]
    @test hasEqual(bf1 .* [1.0, 3.0], [bf1*1.0, bf1*3.0])
    bfm12 = bf1 + genBasisFunc([1,2,2], (2,3))
    @test hasEqual(bfm12 .* [1.0, 3.0], [bfm12*1.0, bfm12*3.0])
    bfs11 = genBasisFunc([1,1,1], (2,3), "P")
    @test hasEqual(bfs11 .* [1.0, 2.0, 3.0], [i*j for (i,j) in zip(bfs11, 1:3)])
end


text:
using RunMyNotes
using Test

@test RunMyNotes.package(RunMyNotes)

@test RunMyNotes.folder(joinpath(@__DIR__, "..", "notes")) # the same!


text:
import JuLIP
using JuLIP: JVec, JMat, Atoms
using JuLIP.MLIPs: IPBasis
using LinearAlgebra: norm

using JuLIP.Potentials: ZList, SZList, z2i, i2z
using StaticArrays: SMatrix

export PolyPairBasis

struct PolyPairBasis{TJ, NZ} <: IPBasis
   J::TJ
   zlist::SZList{NZ}
   bidx0::SMatrix{NZ,NZ,Int16}
end

PolyPairBasis(maxdeg::Integer, trans::DistanceTransform, fcut::PolyCutoff) =
   PolyPairBasis(:X, maxdeg, trans, fcut)

PolyPairBasis(species, maxdeg::Integer,
              trans::DistanceTransform, fcut::PolyCutoff) =
   PolyPairBasis( TransformedJacobi(maxdeg, trans, fcut),
                  ZList(species; static=true) )

PolyPairBasis(J::TransformedJacobi, zlist::SZList) =
      PolyPairBasis(J, zlist, get_bidx0(J, zlist))

function get_bidx0(J, zlist::SZList{NZ}) where {NZ}
   NJ = length(J)
   bidx0 = fill(zero(Int16), (NZ, NZ))
   i0 = 0
   for i = 1:NZ, j = i:NZ
      bidx0[i,j] = i0
      bidx0[j,i] = i0
      i0 += NJ
   end
   return SMatrix{NZ, NZ, Int16}(bidx0...)
end

==(B1::PolyPairBasis, B2::PolyPairBasis) =
      (B1.J == B2.J) && (B1.zlist == B2.zlist)

nz(B::PolyPairBasis) = length(B.zlist)

Base.length(pB::PolyPairBasis) = length(pB.J) * (nz(pB) * (nz(pB) + 1)) ÷ 2

JuLIP.cutoff(pB::PolyPairBasis) = cutoff(pB.J)

Dict(pB::PolyPairBasis) = Dict(
      "__id__" => "PolyPairPots_PolyPairBasis",
      "J" => Dict(pB.J),
      "zlist" => pB.zlist.list )

PolyPairBasis(D::Dict) = PolyPairBasis( TransformedJacobi(D["J"]),
                                        ZList(D["zlist"]; static=true) )

convert(::Val{:PolyPairPots_PolyPairBasis}, D::Dict) = PolyPairBasis(D)

alloc_temp(pB::PolyPairBasis) = (J = alloc_B(pB.J),)
alloc_temp_d(pB::PolyPairBasis, args...) = ( J = alloc_B( pB.J),
                                            dJ = alloc_dB(pB.J) )

"""
compute the zeroth index of the basis corresponding to the potential between
two species zi, zj; as precomputed in `PolyPairBasis.bidx0`
"""
_Bidx0(pB, zi, zj) = pB.bidx0[ z2i(pB, zi), z2i(pB, zj) ]

function energy(pB::PolyPairBasis, at::Atoms{T}) where {T}
   E = zeros(T, length(pB))
   stor = alloc_temp(pB)
   for (i, j, R) in pairs(at, cutoff(pB))
      r = norm(R)
      evaluate!(stor.J, nothing, pB.J, r)
      idx0 = _Bidx0(pB, at.Z[i], at.Z[j])
      for n = 1:length(pB.J)
         E[idx0 + n] += 0.5 * stor.J[n]
      end
   end
   return E
end

function forces(pB::PolyPairBasis, at::Atoms{T}) where {T}
   F = zeros(JVec{T}, length(at), length(pB))
   stor = alloc_temp_d(pB)
   for (i, j, R) in pairs(at, cutoff(pB))
      r = norm(R)
      evaluate_d!(stor.J, stor.dJ, nothing, pB.J, r)
      idx0 = _Bidx0(pB, at.Z[i], at.Z[j])
      for n = 1:length(pB.J)
         F[i, idx0 + n] += 0.5 * stor.dJ[n] * (R/r)
         F[j, idx0 + n] -= 0.5 * stor.dJ[n] * (R/r)
      end
   end
   return [ F[:, iB] for iB = 1:length(pB) ]
end

function virial(pB::PolyPairBasis, at::Atoms{T}) where {T}
   V = zeros(JMat{T}, length(pB))
   stor = alloc_temp_d(pB)
   for (i, j, R) in pairs(at, cutoff(pB))
      r = norm(R)
      evaluate_d!(stor.J, stor.dJ, nothing, pB.J, r)
      idx0 = _Bidx0(pB, at.Z[i], at.Z[j])
      for n = 1:length(pB.J)
         V[idx0 + n] -= 0.5 * (stor.dJ[n]/r) * R * R'
      end
   end
   return V
end


text:
"""
```julia
nanmean(A; dims)
```
Compute the mean of all non-`NaN` elements in `A`, optionally over dimensions
specified by `dims`. As `Statistics.mean`, but ignoring `NaN`s.

As an alternative to `dims`, `nanmean` also supports the `dim` keyword, which
behaves identically to `dims`, but also drops any singleton dimensions that have
been reduced over (as is the convention in some other languages).

## Examples
```julia
julia> using NaNStatistics

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> nanmean(A, dims=1)
1×2 Matrix{Float64}:
 2.0  3.0

julia> nanmean(A, dims=2)
2×1 Matrix{Float64}:
 1.5
 3.5
```
"""
nanmean(A; dims=:, dim=:) = __nanmean(A, dims, dim)
__nanmean(A, ::Colon, ::Colon) = _nanmean(A, :)
__nanmean(A, region, ::Colon) = _nanmean(A, region)
__nanmean(A, ::Colon, region) = _reducedims(_nanmean(A, region), region)
export nanmean

# Reduce one dim
_nanmean(A, dims::Int) = _nanmean(A, (dims,))

# Reduce some dims
function _nanmean(A::AbstractArray{T,N}, dims::Tuple) where {T,N}
    sᵢ = size(A)
    sₒ = ntuple(Val(N)) do d
        ifelse(d ∈ dims, 1, sᵢ[d])
    end
    Tₒ = Base.promote_op(/, T, Int)
    B = similar(A, Tₒ, sₒ)
    _nanmean!(B, A, dims)
end

# Reduce all the dims!
function _nanmean(A, ::Colon)
    Tₒ = Base.promote_op(/, eltype(A), Int)
    n = 0
    Σ = ∅ = zero(Tₒ)
    @turbo for i ∈ eachindex(A)
        Aᵢ = A[i]
        notnan = Aᵢ==Aᵢ
        n += notnan
        Σ += ifelse(notnan, A[i], ∅)
    end
    return Σ / n
end
function _nanmean(A::AbstractArray{<:Integer}, ::Colon)
    Tₒ = Base.promote_op(/, eltype(A), Int)
    Σ = zero(Tₒ)
    @turbo for i ∈ eachindex(A)
        Σ += A[i]
    end
    return Σ / length(A)
end

# Metaprogramming magic adapted from Chris Elrod example:
# Generate customized set of loops for a given ndims and a vector
# `static_dims` of dimensions to reduce over
function staticdim_nanmean_quote(static_dims::Vector{Int}, N::Int)
  M = length(static_dims)
  # `static_dims` now contains every dim we're taking the mean over.
  Bᵥ = Expr(:call, :view, :B)
  reduct_inds = Int[]
  nonreduct_inds = Int[]
  # Firstly, build our expressions for indexing each array
  Aind = :(A[])
  Bind = :(Bᵥ[])
  inds = Vector{Symbol}(undef, N)
  for n ∈ 1:N
    ind = Symbol(:i_,n)
    inds[n] = ind
    push!(Aind.args, ind)
    if n ∈ static_dims
      push!(reduct_inds, n)
      push!(Bᵥ.args, :(firstindex(B,$n)))
    else
      push!(nonreduct_inds, n)
      push!(Bᵥ.args, :)
      push!(Bind.args, ind)
    end
  end
  firstn = first(nonreduct_inds)
  # Secondly, build up our set of loops
  block = Expr(:block)
  loops = Expr(:for, :($(inds[firstn]) = indices((A,B),$firstn)), block)
  if length(nonreduct_inds) > 1
    for n ∈ @view(nonreduct_inds[2:end])
      newblock = Expr(:block)
      push!(block.args, Expr(:for, :($(inds[n]) = indices((A,B),$n)), newblock))
      block = newblock
    end
  end
  rblock = block
  # Push more things here if you want them at the beginning of the reduction loop
  push!(rblock.args, :(n = 0))
  push!(rblock.args, :(Σ = ∅))
  # Build the reduction loop
  for n ∈ reduct_inds
    newblock = Expr(:block)
    push!(block.args, Expr(:for, :($(inds[n]) = axes(A,$n)), newblock))
    block = newblock
  end
  # Push more things here if you want them in the innermost loop
  push!(block.args, :(Aᵢ = $Aind))
  push!(block.args, :(notnan = Aᵢ==Aᵢ))
  push!(block.args, :(n += notnan))
  push!(block.args, :(Σ += ifelse(notnan, Aᵢ, ∅)))
  # Push more things here if you want them at the end of the reduction loop
  push!(rblock.args, :($Bind = Σ * inv(n)))
  # Put it all together
  quote
    ∅ = zero(eltype(B))
    Bᵥ = $Bᵥ
    @turbo $loops
    return B
  end
end

# Chris Elrod metaprogramming magic:
# Turn non-static integers in `dims` tuple into `StaticInt`s
# so we can construct `static_dims` vector within @generated code
function branches_nanmean_quote(N::Int, M::Int, D)
  static_dims = Int[]
  for m ∈ 1:M
    param = D.parameters[m]
    if param <: StaticInt
      new_dim = _dim(param)::Int
      @assert new_dim ∉ static_dims
      push!(static_dims, new_dim)
    else
      t = Expr(:tuple)
      for n ∈ static_dims
        push!(t.args, :(StaticInt{$n}()))
      end
      q = Expr(:block, :(dimm = dims[$m]))
      qold = q
      ifsym = :if
      for n ∈ 1:N
        n ∈ static_dims && continue
        tc = copy(t)
        push!(tc.args, :(StaticInt{$n}()))
        qnew = Expr(ifsym, :(dimm == $n), :(return _nanmean!(B, A, $tc)))
        for r ∈ m+1:M
          push!(tc.args, :(dims[$r]))
        end
        push!(qold.args, qnew)
        qold = qnew
        ifsym = :elseif
      end
      push!(qold.args, Expr(:block, :(throw("Dimension `$dimm` not found."))))
      return q
    end
  end
  staticdim_nanmean_quote(static_dims, N)
end

# Efficient @generated in-place mean
@generated function _nanmean!(B::AbstractArray{Tₒ,N}, A::AbstractArray{T,N}, dims::D) where {Tₒ,T,N,M,D<:Tuple{Vararg{Integer,M}}}
  N == M && return :(B[1] = _nanmean(A, :); B)
  branches_nanmean_quote(N, M, D)
end

## ---