text:
# This is the core functionality for manipulating dimension names as express as a tuple
# It is central to the performance of this named dimensions as a zero cost abstraction
# For constant input all functions should constant propagate to constant output.
# Each function in this file is annotated with a comment and a benchmark
# running that benchmark should return 0 allocations,
# and you can swap its `@btime` for `@code_typed` to check if it is constant propergating
# This constant propagate lets us write faily naive code elsewhere in the package
# knowing the creation and destruction of NamedDimsArrays will be optimised away.

"""
    compile_time_return_hack(x)

This is a cunning hack to strongly encourage the compiler to use constant propagation,
when returning a `Tuple` of `Symbols`.
Under normal circumstances, constant propergation will not work to fully compute
returned tuples of non-bits types.
But because `Symbol`s are a allowed in type-parameters (unlike other nonbits types),
we can put them into a type-parameter, which triggers (??different??) constant propergation,
(for ??Reasons??)
and once something is a type-parameter, then extracting it and returning it is non-allocating.

See https://discourse.julialang.org/t/zero-allocation-tuple-subsetting/23122/8

The short version of this is:
if a function returns a tuple of symbols and is allocating when it looks like it shouldn't
then try wrapping its returned value in `compile_time_return_hack`.
and then look at the `@code_lowered` again.
"""
compile_time_return_hack(x::Tuple{Vararg{Symbol}}) = _compile_time_return_hack(Val{x}())
_compile_time_return_hack(::Val{X}) where X = X

"""
    dim(dimnames, name)

For `dimnames` being a tuple of names (symbols) for the dimensions.
and `name` being a one of those names
This returns the dimension corresponding to that `name`,
e.g. `dim((:a, :b), :b) == 2`
If that `name` is not one of the given `dimnames` then an error is thrown.
"""
function dim(dimnames::Tuple, name::Symbol)::Int
    # 0-Allocations see: `@btime  (()->dim((:a, :b), :b))()`
    dimnum = dim_noerror(dimnames, name)
    if dimnum === 0
        throw(ArgumentError(
            "Specified name ($(repr(name))) does not match any dimension name ($dimnames)"
        ))
    end
    return dimnum
end

function dim(dimnames::Tuple, names)
    # 0-Allocations see: `@btime (()->dim((:a,:b), (:a,:b)))()`
    return map(name->dim(dimnames, name), names)
end

function dim(dimnames::Tuple, d::Union{Integer, Colon})
    # This is the fallback that allows `NamedDimsArray`'s to be have dimensions
    # referred to by number. This is required to allow functions on `AbstractArray`s
    # and that use function like `sum(xs; dims=2)` to continue to work without changes
    # `:` is the default for most methods that take `dims`
    return d
end

Base.@pure function dim_noerror(dimnames::Tuple{Vararg{Symbol, N}}, name::Symbol) where N
    # 0-Allocations see: @btime  (()->dim_noerror((:a, :b, :c), :c))()
    for ii in 1:N
        getfield(dimnames, ii) === name && return ii
    end
    return 0
end

"""
    expand_dimnames(dimnames, name)

For `dimnames` being a tuple of names (symbols) for the dimensions.
and `name` being a name.
This expands the `dimnames` if `name` is not in `dimnames`.
e.g. `expand_dimnames((:a, :b), :c) == (:a, :b, :c)`
If `name` is already in `dimnames` then `dimnames` is returned.
"""
function expand_dimnames(dimnames::Tuple, name::Symbol)
    if dim_noerror(dimnames, name) > 0  # name in dimnames, but optimised
        return dimnames
    else
        return (dimnames..., name)
    end
end

function expand_dimnames(dimnames::Tuple, name::Union{Colon, Tuple{}})
    return dimnames
end

function expand_dimnames(dimnames::Tuple, name::Integer)
    if name <= length(dimnames)
        return dimnames
    else
        extra_length = name - length(dimnames)
        new_dimnames = ntuple(i->:_, extra_length)
        return (dimnames..., new_dimnames...)
    end
end

function expand_dimnames(dimnames::Tuple, names)
    expand_dimnames(expand_dimnames(dimnames, names[1]), names[2:end])
end

"""
    permute_dimnames(dimnames, perm)

Reorder `dimnames` according `perm`.
`perm` should be ordered set of numerical indexs for the new position of the name.
Note: this does not throw errors if you give it a permutation that skips some positions
and duplicates others.
"""
function permute_dimnames(dimnames::NTuple{N, Symbol}, perm) where N
    # 0-Allocations, but does not seem to fully calculate at compile time
    # even with the `compile_time_return_hack`, though that is still required to
    # prevent allocations. `@code_typed permute_dimnames((:a,:b,:c), (1,3,2))`

    new_dimnames = ntuple(length(perm)) do ii
        ind = perm[ii]
        dimnames[ind]
    end
    return compile_time_return_hack(new_dimnames)
end

"""
    tuple_issubset
A version of `is_subset` sepecifically for `Tuple`s of `Symbol`s, that is `@pure`.
This helps it get optimised out of existance. It is less of an abuse of `@pure` than
most of the stuff for making `NamedTuples` work.
"""
Base.@pure function tuple_issubset(lhs::Tuple{Vararg{Symbol,N}}, rhs::Tuple{Vararg{Symbol,M}}) where {N,M}
    N <= M || return false
    for a in lhs
        found = false
        for b in rhs
            found |= a === b
        end
        found || return false
    end
    return true
end

"""
    order_named_inds(Val(names); kw...)
    order_named_inds(Val(names), namedtuple)

Returns the tuple of index values for an array with `names`, when indexed by keywords.
Any dimensions not fixed are given as `:`, to make a slice.
An error is thrown if any keywords are used which do not occur in `nda`'s names.
"""
order_named_inds(val::Val{L}; kw...) where {L} = order_named_inds(val, kw.data)

@generated function order_named_inds(val::Val{L}, ni::NamedTuple{K}) where {L,K}
    tuple_issubset(K, L) || throw(DimensionMismatch("Expected subset of $L, got $K"))
    exs = map(L) do n
        if Base.sym_in(n, K)
            qn = QuoteNode(n)
            :(getfield(ni, $qn))
        else
            :(Colon())
        end
    end
    return Expr(:tuple, exs...)
end

"""
    incompatible_dimension_error(names_a, names_b)
Throws a `DimensionMismatch` explaining that these dimension names are not compatible.
"""
function incompatible_dimension_error(names_a, names_b)
    throw(DimensionMismatch("Incompatible dimension names: $names_a ≠ $names_b"))
end

"""
    unify_names(a, b)
    unify_names(a, b, cs...)

Produces the merged set of names for tuples of names `a` and `b`,
or an error if it is not possibly to unify them.
Then continues with further names `cs`, if any.

Two tuples of names can be unified they are the same length
and if for each position the names are either the same, or one is a wildcard (`:_`).
When combining wildcard with non-wildcard the resulting name is the non-wildcard.
(This is somewhat like the very simplest case of unification in e.g prolog).

For example:
 - `(:a, :b)` and `(:a, :b)` can be unified to give `(:a, :b)`
 - similarly: `(:a, :_)` and `(:a, :b)` can also be unified to give `(:a, :b)`
 - `(:a, :b)` and `(:b, :c)` cannot be unified as the names at each position to not match.
 - `(:a, :b)` and `(:a, :b, :c)` cannot be unified as they have different lengths

This is the type of name combination used for binary array operations.
Where the dimensions of both arrays must be the same.

See also `names_are_unifiable(a, b)`, which returns `true` instead of the merged names,
and `false` instead of an error.
"""
function unify_names(names_a, names_b)
    # @btime (()->unify_names((:a, :b), (:a, :_)))()
    ret = try_unify_names(names_a, names_b)
    if ret === nothing
        incompatible_dimension_error(names_a, names_b)
    else
        return compile_time_return_hack(ret)
    end
end
unify_names(a) = a
unify_names(a, b, cs...) = unify_names(unify_names(a, b), cs...)
# @btime (()->unify_names((:a, :b), (:a, :_), (:_, :b)))()

names_are_unifiable(names_a, names_b) = try_unify_names(names_a, names_b) !== nothing

function try_unify_names(names_a, names_b)
    if names_a === names_b
        return names_a
    elseif length(names_a) !== length(names_b)
        return nothing
    end

    ret = ntuple(length(names_a)) do ii  # remove :_ wildcards
        a = getfield(names_a, ii)
        b = names_b[ii]
        a === :_ && return b
        b === :_ && return a
        a === b && return a
        return false  # mismatch occured, we mark this with a non-Symbol result
    end

    if ret isa Tuple{Vararg{Symbol}}
        return compile_time_return_hack(ret)
    else
        return nothing
    end
end

"""
    unify_names_longest(a, b)

This is the same as [`unify_names`](@ref), but with the equal length requirement removed.
It unifies the names up to the length of the shortest, and takes the named from the longest
for the remainder.
It can also be considered as padding the shorter of the two given tuples of names
with trailing wildcards (`:_`).
This is the type of name combination used for broadcating array operations.
Where the smaller (dimensionally) array is broadcast against the longer, repeating it for
all entries in the missing trailing dimensions.
"""
unify_names_longest(names, ::Tuple{}) = names
unify_names_longest(::Tuple{}, names) = names
unify_names_longest(::Tuple{}, ::Tuple{}) = tuple()
function unify_names_longest(names_a, names_b)
    # 0 Allocations: @btime (()-> unify_names_longest((:a,:b), (:a,)))()

    length(names_a) === length(names_b) && return unify_names(names_a, names_b)
    long, short = length(names_a) > length(names_b) ? (names_a, names_b) : (names_b, names_a)
    ret = ntuple(length(long)) do ii
        a = getfield(long, ii)
        b = ii <= length(short) ? short[ii] : :_
        a === :_ && return b
        b === :_ && return a
        a === b && return a
        return false  # mismatch occured, we mark this with a nonSymbol result
    end
    ret isa Tuple{Vararg{Symbol}} || incompatible_dimension_error(names_a, names_b)
    return compile_time_return_hack(ret)
end

function unify_names_longest(names_a, names_b, names_c::Vararg{<:NTuple{Nd, Symbol} where Nd, N}) where N
    return unify_names_longest(names_a, unify_names_longest(names_b, names_c...))
end

unify_names_shortest(names, ::Tuple{}) = ()
unify_names_shortest(::Tuple{}, names) = ()
unify_names_shortest(::Tuple{}, ::Tuple{}) = ()
function unify_names_shortest(names_a, names_b)
    # 0 Allocations: @btime (()-> unify_names_shortest((:a,:b), (:a,)))()

    length(names_a) === length(names_b) && return unify_names(names_a, names_b)
    long, short = length(names_a) > length(names_b) ? (names_a, names_b) : (names_b, names_a)
    ret = ntuple(length(short)) do ii
        a = getfield(long, ii)
        b = getfield(short, ii)
        a === :_ && return b
        b === :_ && return a
        a === b && return a
        return false  # mismatch occured, we mark this with a nonSymbol result
    end
    ret isa Tuple{Vararg{Symbol}} || incompatible_dimension_error(names_a, names_b)
    return compile_time_return_hack(ret)
end

function unify_names_shortest(names_a, names_b, names_c::Vararg{<:NTuple{Nd, Symbol} where Nd, N}) where N
    return unify_names_shortest(names_a, unify_names_shortest(names_b, names_c...))
end

# The following are helpers for remaining_dimnames_from_indexing
# as a generated function it can get unhappy if asked to use anon functions
# and it can only call function declared before it. So we declare them explictly here.
is_noninteger_type(::Type{<:Integer}) = false
is_noninteger_type(::Any) = true

"""
    remaining_dimnames_from_indexing(dimnames::Tuple, inds)

Given a tuple of dimension names,
and a tuple of indices e.g `(1, :, 1:3, [true, false])`,
this drops those indexed with scalars or `CartesianIndex`,
inserts `:_` for `newaxis = [CartesianIndex{0}()]`,
and returns another tuple of names.

It will return an empty tuple to indicate that all names should be dropped.
This happend for scalar indexing by integers, or one `CartesianIndex`.
It also happens e.g. when indexing a matrix by a `BitArray{2}` such as `mat[mat .> 0.5]`:
this returns a vector, the same as vec(mat)[vec(mat .> 0.5)], whose dimension isn't any
of the original dimensions, hence has no name.
"""
@generated function remaining_dimnames_from_indexing(dimnames::Tuple, inds::Tuple)
    # 0-Allocation see:
    # `@btime (()->remaining_dimnames_from_indexing((:a, :b, :c), (:,390,:)))()``
    keep_names = []
    dim_num = 0
    for type in inds.parameters
        if type <: Integer
            dim_num += 1
        elseif type <: CartesianIndex
            dim_num += type.parameters[1]
        elseif type <: AbstractVector{CartesianIndex{0}}
            push!(keep_names, QuoteNode(:_))
        elseif type <: AbstractArray{<:Integer} && ndims(type) > 1
            dim_num += 1
            for _ in 1:ndims(type)
                push!(keep_names, QuoteNode(:_))
            end
        else
            dim_num += 1
            push!(keep_names, :(getfield(dimnames, $dim_num)) )
        end
    end
    return Expr(:call, :compile_time_return_hack, Expr(:tuple, keep_names...))
end

remaining_dimnames_from_indexing(dn::Tuple, inds::Tuple{Vararg{<:Integer}}) = ()
remaining_dimnames_from_indexing(dn::Tuple, ci::Tuple{CartesianIndex}) = ()

function remaining_dimnames_from_indexing(
    dimnames::Tuple{<:Any, <:Any, Vararg}, inds::Tuple{T}
) where T <: Union{Base.LogicalIndex, AbstractVector{<:CartesianIndex}}
    return ()
end

"""
    remaining_dimnames_after_dropping(dimnames::Tuple, dropped_dims)
Given a tuple of dimension names, and either a collection of dimensions,
or a single dimension, expressed as a number,
Returns the dimension names with those dimensions dropped.
"""
function remaining_dimnames_after_dropping(dimnames::Tuple, dropped_dim::Int)
    # 0 allocations. See `@btime remaining_dimnames_after_dropping((:a,:b,:c,:d,:e), 4)`
    return _remaining_dimnames_after_dropping(dimnames, Tuple{dropped_dim})
end

function remaining_dimnames_after_dropping(dimnames::NTuple{N, Symbol}, dropped_dims::Tuple{Vararg{Int}}) where N
    # 0-Allocations see:
    # `@code_typed (()->remaining_dimnames_after_dropping((:a,:b,:c,:d,:e), (1,3)))()`
    return _remaining_dimnames_after_dropping(dimnames, Tuple{dropped_dims...})
end

# dropped_dims must be a Tuple type where the values are Int literal for the dimensions being dropped
@generated function _remaining_dimnames_after_dropping(dimnames::NTuple{N,Symbol}, dropped_dims::Type) where N
    dropped_dims_vals = dropped_dims.parameters[1].parameters
    keep_names = [:(getfield(dimnames, $ii)) for ii in 1:N if ii ∉ dropped_dims_vals]
    return Expr(:call, :compile_time_return_hack, Expr(:tuple, keep_names...))
end

"""
    tuple_cat(x, y, zs...)

This is like `vcat` for tuples, it splats everything into one long tuple.
"""
tuple_cat(x::Tuple, ys::Tuple...) = (x..., tuple_cat(ys...)...)
tuple_cat() = ()
# @btime tuple_cat((1, 2), (3, 4, 5), (6,)) # 0 allocations


text:
## IMPORTS ##
using Pkg
Pkg.activate(joinpath(@__DIR__,"../.."))
using Dates
include(joinpath(@__DIR__,"model.jl"))
include(joinpath(@__DIR__,"utilities.jl"))

## SAVING INFO ##
DO_SAVE = false
INTERIM_SAVE = false

## PARAMETERS ##

# Model
Nx = 4
Ny = 1
t = 1
U = 8
ω = 0.5*t
g0 = 0.1*t
g1 = 0.05*t
λ = 0*t
doping = 0
max_phonons = 5

# Save
date_stamp = Dates.format(now(), "y-m-d_HH:MM:SS")
param_stamp = "_$(Nx)Nx_$(Ny)Ny_$(t)t_$(U)U_$(ω)ω_$(g0)g0_$(g1)g1_$(λ)λ_$(doping)doping_$(max_phonons)phonons"
save_path = joinpath(@__DIR__,"outputs",date_stamp*param_stamp*".h5")

# Simulation
DMRG_numsweeps = 20
DMRG_maxdim = 1000
DMRG_cutoff = 1E-10
DMRG_LBO = true
T = 1#80
τ = 0.05
TEBD_maxdim = 1000
TEBD_cutoff = 1E-10
TEBD_LBO = true

# Specify spectral function operators
A_t0 = "Cup"
A_t = "Cdagup"

## CODE ##

# Initialize
println("Initializing...")
params = parameters(Nx=Nx, Ny=Ny, t=t, U=U, ω=ω, g0=g0, g1=g1, λ=λ, doping=doping,
                    max_phonons=max_phonons,
                    DMRG_numsweeps=DMRG_numsweeps,
                    DMRG_maxdim=DMRG_maxdim, DMRG_cutoff=DMRG_cutoff,DMRG_LBO=DMRG_LBO,
                    T=T, τ=τ, TEBD_cutoff=TEBD_cutoff,TEBD_LBO=TEBD_LBO)
hubbholst = HubbardHolsteinModel(params)
if DO_SAVE
    save_structs(params, save_path)
    save_structs(hubbholst, save_path)
end

# Run DMRG
println("Finding ground state...")
dmrg_results = run_DMRG(hubbholst, params, alg="divide_and_conquer")
if DO_SAVE
    save_structs(dmrg_results, save_path)
end

# Equilibrium correlations
println("Computing equilibrium correlations...")
eq_corr = compute_all_equilibrium_correlations(dmrg_results, hubbholst)
if DO_SAVE
    save_structs(eq_corr, save_path)
end

# Compute correlation functions
println("Computing correlation functions...")
tebd_results = compute_correlations(dmrg_results, A_t0, A_t, hubbholst, params,
                                    interim_save=INTERIM_SAVE, savepath=save_path)
if DO_SAVE
    save_structs(tebd_results, save_path)
end

ext: jl


text:
using PyCall

@pyimport keras.activations as activations
@pyimport keras.callbacks as callbacks
@pyimport keras.objectives as objectives
@pyimport keras.preprocessing as preprocessing
@pyimport keras.utils as utils
@pyimport keras.backend as backend
@pyimport keras.constraints as constraints
@pyimport keras.initializations as initializations
@pyimport keras.models as models
@pyimport keras.optimizers as optimizers
@pyimport keras.regularizers as regularizers
@pyimport keras.wrappers as wrappers
@pyimport keras.layers as layers

# Potential mode for how to do this in the future.
# Perhaps the structure can be recovered automatically so we don't have to hand each module
module datasets
using PyCall
for dataset in [:cifar, :cifar10, :cifar100, :imdb, :mnist, :reuters]
    @eval @pyimport keras.datasets.$dataset as $dataset
end
end


text:
immutable NoncentralT{T<:Real} <: ContinuousUnivariateDistribution
    ν::T
    λ::T

    function (::Type{NoncentralT{T}}){T}(ν::T, λ::T)
        @check_args(NoncentralT, ν > zero(ν))
        @check_args(NoncentralT, λ >= zero(λ))
        new{T}(ν, λ)
    end
end

NoncentralT{T<:Real}(ν::T, λ::T) = NoncentralT{T}(ν, λ)
NoncentralT(ν::Real, λ::Real) = NoncentralT(promote(ν, λ)...)
NoncentralT(ν::Integer, λ::Integer) = NoncentralT(Float64(ν), Float64(λ))

@distr_support NoncentralT -Inf Inf

### Conversions
convert{T <: Real, S <: Real}(::Type{NoncentralT{T}}, ν::S, λ::S) = NoncentralT(T(ν), T(λ))
convert{T <: Real, S <: Real}(::Type{NoncentralT{T}}, d::NoncentralT{S}) = NoncentralT(T(d.ν), T(d.λ))

### Parameters

params(d::NoncentralT) = (d.ν, d.λ)
@inline partype{T<:Real}(d::NoncentralT{T}) = T

### Statistics

function mean{T<:Real}(d::NoncentralT{T})
    if d.ν > 1
        isinf(d.ν) ? d.λ :
        sqrt(d.ν/2) * d.λ * gamma((d.ν - 1)/2) / gamma(d.ν/2)
    else
        T(NaN)
    end
end

function var{T<:Real}(d::NoncentralT{T})
    d.ν > 2 ? d.ν*(1 + d.λ^2) / (d.ν - 2) - mean(d)^2 : T(NaN)
end

### Evaluation & Sampling

@_delegate_statsfuns NoncentralT ntdist ν λ

function rand(d::NoncentralT)
    z = randn()
    v = rand(Chisq(d.ν))
    (z+d.λ)/sqrt(v/d.ν)
end


text:
# eq
concept_eq(x) = x[1] == x[2]

"""
    eq
Equality between two variables.
"""
const eq = Constraint(
    args_length = 2,
    concept = concept_eq,
    error = make_error(:eq),
)

# _all_equal
concept_all_equal(x) = all(y -> y == x[1], x)

function error_all_equal(x::V; param = nothing, dom_size = 0
) where {T <: Number, V <: AbstractVector{T}}
    acc = Dictionary{T,Int}()
    foreach(y -> insert_or_inc(acc, y), x)
    return Float64(length(x) - maximum(acc))
end

"""
    all_equal
Global constraint ensuring that all the values of `x` are all equal.
"""
const all_equal = Constraint(
    concept = concept_all_equal,
    error = make_error(:all_equal),
)

# all_equal_param
concept_all_equal_param(x; param) = all(y -> y == param, x)
error_all_equal_param(x; param, dom_size = 0) = count(y -> y != param, x)

"""
    all_equal_param
Global constraint ensuring that all the values of `x` are all equal to a given parameter `param`.
"""
const all_equal_param = Constraint(
    concept = concept_all_equal_param,
    error = make_error(:all_equal_param),
    param = 1,
)

# _sum_equal_param
concept_sum_equal_param(x; param) = sum(x) == param
error_sum_equal_param(x; param, dom_size) = abs(sum(x) - param) / dom_size

"""
    sum_equal_param
Global constraint ensuring that the sum of the values of `x` is equal to a given parameter `param`.
"""
const sum_equal_param = Constraint(
    concept = concept_sum_equal_param,
    error = make_error(:sum_equal_param),
    param = 1,
)

# _minus_equal_param
concept_minus_equal_param(x; param) = x[1] == param + x[2]
error_minus_equal_param(x; param, dom_size) = abs(x[1] - x[2] - param) / dom_size

"""
    minus_equal_param
Constraint ensuring that the difference between `x[1]` and `x[2]` is equal to a given parameter `param`.
"""
const minus_equal_param = Constraint(
    args_length = 2,
    concept = concept_minus_equal_param,
    error = make_error(:minus_equal_param),
    param = 1,
)


text:
using LogDensityProblems, Test, Distributions, TransformVariables
import LogDensityProblems: capabilities, dimension, logdensity
using LogDensityProblems: logdensity_and_gradient, LogDensityOrder

import ForwardDiff, Flux, Zygote, TransformVariables, Random
using Parameters: @unpack

Random.seed!(1)

####
#### transformed Bayesian problem
####

@testset "transformed Bayesian problem" begin
    t = as((y = asℝ₊, ))
    #d = LogNormal(1.0, 2.0)
    d = Normal(1.0, 2.0)
    logposterior = ((x, ), ) -> logpdf(d, x)

    # a Bayesian problem
    p = TransformedLogDensity(t, logposterior)
    @test repr(p) == "TransformedLogDensity of dimension 1"
    @test dimension(p) == 1
    @test p.transformation ≡ t
    @test capabilities(p) == LogDensityOrder(0)

    # gradient of a problem
    ∇p = ADgradient(:ForwardDiff, p)
    @test dimension(∇p) == 1
    @test parent(∇p).transformation ≡ t

    for _ in 1:100
        x = random_arg(t)
        θ, lj = transform_and_logjac(t, x)
        px = logdensity(p, x)
        @test logpdf(d, θ.y) + lj ≈ (px::Real)
        px2, ∇px = logdensity_and_gradient(∇p, x)
        @test px2 == px
        @test ∇px ≈ [ForwardDiff.derivative(x -> logpdf(d, exp(x)) + x, x[1])]
    end
end


text:
println("Use LM to fit the parameter")
@test_broken begin
    fit = lm_fit(prob1,t,vec(data),[1.0],Tsit5(),show_trace=false,lambda=10000.0)
    param = fit.param
    @test param[1] ≈ 1.5 atol=1e-3

    fit = lm_fit(prob2,t,vec(data),[1.3,2.6],Tsit5(),show_trace=false,lambda=10000.0)
    param = fit.param
    @test param ≈ [1.5; 3.0] atol=0.002

    fit = lm_fit(prob3,t,vec(data),[1.3,0.8,2.6,1.2],Tsit5(),show_trace=false,lambda=10000.0)
    param = fit.param
    @test param ≈ [1.5;1.0;3.0;1.0] atol=1e-2
end


text:
include("../required.jl")

include("../07_comparison_data.jl")

nspecies = size(raw.Y, 2)
nsites_raw = format(length(raw.richness); commas=true)
nsites_sdm_obs = format(length(sdm.richness); commas=true)
pres_df = CSV.read(
        joinpath("data", "proc", "bart_predictions_pres.csv"), DataFrame; missingstrings=["NA"]
    );
nsites_sdm_tot = nrow(dropmissing(pres_df))

lcbdraw_min, lcbdraw_max = sprintf1.("%.3e", extrema(raw.lcbd))
lcbdsdm_min, lcbdsdm_max = sprintf1.("%.3e", extrema(sdm.lcbd))

bdtot_raw = sprintf1("%.3f", raw.beta_total)
bdtot_sdm = sprintf1("%.3f", sdm.beta_total)

##

R"""
library(SpatialPack)
load(file.path("data", "rdata", "cor-dutilleul.Rdata"))

z_richness_values = list(
    corr = z_richness$corr,
    Fstat = z_richness$Fstat * z_richness$dof, # as reported in summary
    pvalue = z_richness$p.value
)
z_lcbd_values = list(
    corr = z_lcbd$corr,
    Fstat = z_lcbd$Fstat * z_lcbd$dof, # as reported in summary
    pvalue = z_lcbd$p.value
)
"""

@rget z_richness_values z_lcbd_values

sprintf1("%.3f", z_richness_values[:corr])
sprintf1("%.3f", z_richness_values[:Fstat])
sprintf1("%.3e", z_richness_values[:pvalue])
sprintf1("%.3f", z_lcbd_values[:corr])
sprintf1("%.3f", z_lcbd_values[:Fstat])
sprintf1("%.3e", z_lcbd_values[:pvalue])

include("../09_comparison_plots.jl")

richdiff_min, richdiff_max = sprintf1.("%i", extrema(richness_diff))

lcbdcomraw_min, lcbdcomraw_max = sprintf1.("%.3e", extrema(raw.lcbdcom))
lcbdcomsdm_min, lcbdcomsdm_max = sprintf1.("%.3e", extrema(sdm.lcbdcom))

lcbddiff_min, lcbddiff_max = sprintf1.("%.3e", extrema(lcbd_diff))

richres_min, richres_max = sprintf1.("%.3f", extrema(richres_nb_layer))
lcbdres_min, lcbdres_max = sprintf1.("%.3f", extrema(lcbdres_br_layer))

##

include("../05_subareas.jl")

coords_NE = (left=-80.0, right=-60.0, bottom=40.0, top=50.0)
coords_SW = (left=-120.0, right=-100.0, bottom=30.0, top=40.0)

bdtot_NE = sprintf1("%.3f", beta_NE)
bdtot_SW = sprintf1("%.3f", beta_SW)

lcbdne_min, lcbdne_max = sprintf1.("%.3e", extrema(lcbd_NE))
lcbdsw_min, lcbdsw_max = sprintf1.("%.3e", extrema(lcbd_SW))

bdtot_s1, bdtot_s2, bdtot_s3 = sprintf1.("%.3f", beta_values[[1, mid_ind, end]])

lcbdmin_s1 = sprintf1("%.3e", lcbd_mins[1])
lcbdmax_s1 = sprintf1("%.3e", lcbd_maxs[1])
lcbdmin_s2 = sprintf1("%.3e", lcbd_mins[mid_ind])
lcbdmax_s2 = sprintf1("%.3e", lcbd_maxs[mid_ind])
lcbdmin_s3 = sprintf1("%.3e", lcbd_mins[end])
lcbdmax_s3 = sprintf1("%.3e", lcbd_maxs[end])

ascending_threshold_NE = minimum(lcbd_NE)
sprintf1("%.3e", ascending_threshold_NE)
min_indx_NE = findall(x -> x == ascending_threshold_NE, lcbd_NE.grid)
abs_min_NE = median(richness_NE.grid[min_indx_NE])

ascending_threshold_SW = minimum(lcbd_SW)
sprintf1("%.3e", ascending_threshold_SW)
min_indx_SW = findall(x -> x == ascending_threshold_SW, lcbd_SW.grid)
abs_min_SW = median(richness_SW.grid[min_indx_SW])

# medrichmin
# medrichmax
# medlcbdmin
# medlcbdmax
# medgammamin
# medgammamax
# medbdtotmin
# medbttotmax

##

cols = [
    "SCIENTIFIC NAME",
    "LATITUDE",
    "LONGITUDE",
    "OBSERVATION DATE",
    "GROUP IDENTIFIER",
    "APPROVED",
    "SAMPLING EVENT IDENTIFIER"
]

# Load data from CSV files (from file cut with terminal)
@time df = CSV.read(joinpath("data", "raw", "ebd_warblers_cut.csv"), DataFrame; select=cols)

# Prepare data (arrange values & columns)
@time prepare_ebd_data!(df)

# Select subset with specific columns
select!(df, [:species, :latitude, :longitude, :groupIdentifier, :samplingEventIdentifier])

# Remove 1 Aleutian Islands observation with positive longitude
filter!(:longitude => <(0.0), df)

# Remove duplicates (BUT not ok for counts, so-so for dates, see group-observation.jl in src/test/)
n_obs = nrow(df)
df_groups = filter(:groupIdentifier => !ismissing, df)
unique!(df_groups, [:species, :groupIdentifier])
filter!(:groupIdentifier => ismissing, df)
append!(df, df_groups)

n_obs_tot = format(n_obs; commas=true)
n_obs_u = format(nrow(df), commas=true)
n_checklists = format(length(unique(df.samplingEventIdentifier)), commas=true)

# nobs_tot = 27_821_881
# nobs_u = 22_974_330
# nchecklists = 9_103_750

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
using Robotlib, Robotlib.Frames, MAT, Test

@testset "Frames" begin

R = toOrthoNormal(randn(3,3))
t = randn(3)
f = Frame(R,t)
f = Frame([R t;0 0 0 1])
F = f

path = joinpath(dirname(@__FILE__),"..","src","applications","frames/")
add_frame_name!("SEAM","Weld seam frame")
add_frame_name!("TAB","Table frame")

T_RB_Tm = MAT.matread(path*"T_RB_T.mat")["T_RB_T"]
T_TF_TCPm = MAT.matread(path*"T_TF_TCP.mat")["T_TF_TCP"]
T_T_TABm = MAT.matread(path*"T_T_Table.mat")["T_T_Table"]

T_RB_T = Frame(T_RB_Tm,"RB","T")
T_S_D = Frame(T_TF_TCPm,"S","D")
T_T_TAB = Frame(T_T_TABm,"T","TAB")

cloud_seam = readcloud(path*"CloudSeam_edge.txt")
plane_seam = readplane(path*"PlaneSeam_edge.txt")
cloud_seam_projected = project(plane_seam,cloud_seam)
line_seam = fitline(cloud_seam_projected)

T_T_SEAM = framefromfeatures(("z+",line_seam),("y-",plane_seam),cloud_seam_projected[1],"SEAM")
T_RB_SEAM = T_RB_T*T_T_SEAM
T_RB_TAB = T_RB_T*T_T_TAB
T_TAB_SEAM = inv(T_T_TAB)*T_T_SEAM

cloud_seam_RB = T_RB_T*cloud_seam
cloud_seam_projected_RB = T_RB_T*cloud_seam_projected
plane_seam_RB = T_RB_T*plane_seam
line_seam_RB = T_RB_T*line_seam

# plot(Frame(I4,"RB","U"),200)
# plot!(cloud_seam_RB,label="Cloud seam_RB")
# plot!(cloud_seam_projected_RB,label="Cloud seam_RB projected")
# plot!(line_seam_RB,500,label="Line seam")
# plot!(plane_seam_RB,200,label="Plane seam")
# plot!(T_RB_SEAM,200)
# plot!(T_RB_TAB,200)
# gui()

T = convert(Matrix, f)
@test convert(Array,f) == T
@test promote(T,f) == (T,T)
@test promote_type(Matrix{Float64}, Frame{Float64}) == Matrix{Float64}

end


text:
# Requires the create_testfuncs function defined in test_fft.jl

# -----------------------------------------------------------------------------
# IFFT's TEST FUNCTIONS

rtol = 1e-12

function test_ifft_cosmx(g::OneDGrid)
    # test fft for cos(mx+φ)
    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th = create_testfuncs(g)
    f1b = real(ifft(f1h))
    isapprox(f1, f1b, rtol=rtol)
end

function test_irfft_cosmx(g::OneDGrid)
    # test irfft for cos(mx+φ)
    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th = create_testfuncs(g)
    f1hr_c = deepcopy(f1hr)
    # deepcopy is needed because FFTW irfft messes up with input!
    f1b = irfft(f1hr_c, nx)
    isapprox(f1, f1b, rtol=rtol)
end

function test_irfft_mul_cosmx(g::OneDGrid)
    # test irfft with ldiv! for cos(mx+φ)
    f1, f1h, f1hr, f1hr_mul, f1h_th, f1hr_th = create_testfuncs(g)
    f1hr_c = deepcopy(f1hr)
    # deepcopy is needed because FFTW irfft messes up with input!
    f1b = Array{Float64}(undef, g.nx)
    ldiv!( f1b, g.rfftplan, f1hr_c )
    isapprox(f1, f1b, rtol=rtol)
end

function test_ifft_cosmxcosny(g::TwoDGrid)
    # test fft for cos(mx)cos(ny)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    f1b = real(ifft(f1h))
    isapprox(f1, f1b, rtol=rtol)
end

function test_irfft_cosmxcosny(g::TwoDGrid)
    # test irfft for cos(mx)cos(ny)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    f1hr_c = deepcopy(f1hr)
    # deepcopy is needed because FFTW irfft messes up with input!
    f1b = irfft(f1hr_c, nx)
    isapprox(f1, f1b, rtol=rtol)
end

        function test_irfft_mul_cosmxcosny(g::TwoDGrid)
    # test irfft with ldiv! for cos(mx)cos(ny)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    f1hr_c = deepcopy(f1hr)
    # deepcopy is needed because FFTW irfft messes up with input!
    f1b = Array{Float64}(undef, g.nx, g.ny)
    ldiv!( f1b, g.rfftplan, f1hr_c )
    isapprox(f1, f1b, rtol=rtol)
end

function test_ifft_sinmxny(g::TwoDGrid)
    # test fft for sin(mx+ny)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    f2hr_c = deepcopy(f2hr)
    # deepcopy is needed because FFTW irfft messes up with input!
    f2b3 = irfft(f2hr_c, g.nx)
    isapprox(f2, f2b3, rtol=rtol)
end

function test_irfft_sinmxny(g::TwoDGrid)
    # test irfft for sin(mx+ny)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    f2b4 = real(ifft(f2h))
    isapprox(f2, f2b4, rtol=rtol)
end

function test_irfft_mul_sinmxny(g::TwoDGrid)
    # test irfft with ldiv! for sin(mx+ny)
    f1, f2, f1h, f2h, f1hr, f2hr, f1hr_mul, f2hr_mul,
            f1h_th, f1hr_th, f2h_th, f2hr_th = create_testfuncs(g)
    f2hr_c = deepcopy(f2hr)
    # deepcopy is needed because FFTW irfft messes up with input!
    f2b = Array{Float64}(undef, g.nx, g.ny)
    ldiv!(f2b, g.rfftplan, f2hr_c)
    isapprox(f2, f2b, rtol=rtol)
end

# -----------------------------------------------------------------------------
# Running the tests

# Test 1D grid
nx = 32             # number of points
Lx = 2π             # Domain width
g1 = OneDGrid(nx, Lx)

# Test 2D rectangular grid
nx, ny = 32, 64     # number of points
Lx, Ly = 2π, 3π     # Domain width
g2 = TwoDGrid(nx, Lx, ny, Ly)

@test test_ifft_cosmx(g1)
@test test_irfft_cosmx(g1)
@test test_irfft_mul_cosmx(g1)

@test test_ifft_cosmxcosny(g2)
@test test_irfft_cosmxcosny(g2)
@test test_irfft_mul_cosmxcosny(g2)
@test test_ifft_sinmxny(g2)
@test test_irfft_sinmxny(g2)
@test test_irfft_mul_sinmxny(g2)


text:
using Printf
using Parameters
using StaticArrays
using BenchmarkTools
using InteractiveUtils
import Base: display

###
function string_to_index(str::String)
    return parse(Int, reverse(str); base=2)
end

function index_to_string(index::Int)
    return [parse(Int, ss) for ss in reverse(bitstring(index))]
end

"""
Type to organize all the configuration DeterminantString
"""
#struct DeterminantString
mutable struct DeterminantString
    no::Int
    ne::Int
    sign::Int
    lin_index::Int
    config::Array{Int,1}
    #ca_lookup::Array{Int, 3}
    max::Int
end

function DeterminantString(no::Int, ne::Int)
    return DeterminantString(no, ne, 1, 1, Vector(1:ne), get_nchk(no,ne))
    #return DeterminantString(no, ne, 1, 1, Vector(1:ne), zeros(get_nchk(no,ne),no,no), get_nchk(no,ne))
end

function display(d::DeterminantString)
    print(d)
end

"""
    get_unoccupied(c::DeterminantString)

get list of orbitals that are unoccupied
"""
@inline function get_unoccupied(c::DeterminantString)
    unoccupied = zeros(Int, c.no-c.ne,1)
    j = 1
    for i in 1:c.no
        if i in c.config
            continue
        end
        unoccupied[j] = i
        j+=1
    end
    return unoccupied
end

"""
    get_unoccupied(c::DeterminantString)

get list of orbitals that are unoccupied
"""
function get_unoccupied!(unoccupied::Array{Int,1}, c::DeterminantString)
    j = 1
    for i in 1:c.no
        if i in c.config
            continue
        end
        unoccupied[j] = i
        j+=1
    end
end

function Base.length(c::DeterminantString)
    #=
    return number of DeterminantStrings
    =#
    return c.max
end

"""
    Base.copy!(c1::DeterminantString,c2::DeterminantString)

copy c2 into c1
"""
function Base.copy!(c1::DeterminantString,c2::DeterminantString)
    c1.config = copy(c2.config)
    c1.no = c2.no
    c1.ne = c2.ne
    c1.sign = c2.sign
    c1.lin_index = c2.lin_index
    c1.max = c2.max
end

function Base.print(c::DeterminantString)
    #=
    Pretty print of an determinant DeterminantString
    =#
    @printf("Index: %-10i (%3i,%-3i) Dim: %-10i Sign: %2i ",c.lin_index, c.no, c.ne, c.max, c.sign)
    print(" Config:")
    [@printf("%3i",i) for i in c.config]
    print('\n')
end

function incr!(c::DeterminantString)
    #=
    Increment determinant DeterminantString
    =#
    #={{{=#
    if c.max == nothing
        calc_max!(c)
    end
    if c.lin_index == c.max
        return
    end
    c.lin_index += 1
    incr_comb!(c.config, c.no)
end
#=}}}=#

function calc_max(c::DeterminantString)
    #=
    Calculate dimension of space accessible to a DeterminantString
    =#
    #={{{=#
    return get_nchk(c.no,c.ne)
end
#=}}}=#

function calc_max!(c::DeterminantString)
    #=
    Calculate dimension of space accessible to a DeterminantString
    =#
    #={{{=#
    c.max = get_nchk(c.no,c.ne)
end
#=}}}=#

function incr_comb!(comb::Array{Int,1}, Mend::Int)
    #=
    For a given combination, form the next combination
    =#
    #={{{=#
    N = length(comb)
    for i in N:-1:1
        if comb[i] < Mend - N + i
            comb[i] += 1
            for j in i+1:N
                comb[j]=comb[j-1]+1
            end
            return
        end
    end
    return
end
#=}}}=#

"""
    calc_linear_index!(c::DeterminantString)

Calculate the linear index
"""
function calc_linear_index!(c::DeterminantString)
    #={{{=#
    c.lin_index = 1
    v_prev::Int = 0

    for i::Int in 1:c.ne
        v = c.config[i]
        for j::Int in v_prev+1:v-1
            c.lin_index += StringCI.binom_coeff[c.no-j+1,c.ne-i+1]
            #@btime $c.lin_index += $binom_coeff[$c.no-$j+1,$c.ne-$i+1]
            #c.lin_index += get_nchk(c.no-j,c.ne-i)
        end
        v_prev = v
    end
end
#=}}}=#

"""
    calc_linear_index!(c::DeterminantString, binomcoeff::Array{Int,2})

Calculate the linear index, passing in binomial coefficient matrix makes it much faster
"""
function calc_linear_index!(c::DeterminantString, binomcoeff::Array{Int,2})
    #={{{=#
    c.lin_index = 1
    v_prev::Int = 0

    for i::Int in 1:c.ne
        v = c.config[i]
        for j::Int in v_prev+1:v-1
            c.lin_index += binomcoeff[c.no-j+1,c.ne-i+1]
        end
        v_prev = v
    end
end
#=}}}=#

"""
    calc_linear_index(c::DeterminantString)

Return linear index for lexically ordered __config DeterminantString
"""
function calc_linear_index(c::DeterminantString)
    #={{{=#
    v_prev::Int = 0
    lin_index = 1
    for i in 1:c.ne
        v = c.config[i]
        #todo: change mchn from function call to data lookup!
        @inbounds @simd for j in v_prev+1:v-1
            lin_index += get_nchk(c.no-j,c.ne-i)
        end
        v_prev = v
    end
    return lin_index
end
#=}}}=#

"""
    fill_ca_lookup(c::DeterminantString)

Create an index table relating each DeterminantString with all ia substitutions
i.e., ca_lookup[Ka][c(p) + a(p)*n_p] = La
"""
function fill_ca_lookup(c::DeterminantString)
    #={{{=#

    ket = DeterminantString(c.no, c.ne)
    bra = DeterminantString(c.no, c.ne)

    max = calc_max(ket)

    tbl = []
    for K in 1:max
        Kv::Array{Tuple{Int,Int},1} = []
        for p in 1:ket.no
            for q in 1:ket.no
                bra = deepcopy(ket)
                apply_annihilation!(bra,p)
                apply_creation!(bra,q)
                @assert(issorted(bra.config))
                #print("--\n")
                #print(p,'\n')
                #print(q,'\n')
                #print(ket)
                #print(bra)
                if bra.sign == 0
                    push!(Kv,(1,0))
                    continue
                else
                    #calc_max!(bra)
                    #print(bra)
                    calc_linear_index!(bra)
                    push!(Kv,(bra.sign , bra.lin_index))
                end
            end
        end
        push!(tbl,Kv)
        incr!(ket)
    end
    return tbl
end
#=}}}=#

"""
    fill_ca_lookup2(c::DeterminantString)

Create an index table relating each DeterminantString with all ia substitutions
i.e., ca_lookup[Ka,p,q] = sign*La

<L|p'q|K> = sign
"""
function fill_ca_lookup2(c::DeterminantString)
    #={{{=#

    ket = DeterminantString(c.no, c.ne)
    bra = DeterminantString(c.no, c.ne)

    max = calc_max(ket)

    tbl = zeros(Int,ket.no, ket.no, max)
    for K in 1:max
        for q in 1:ket.no
            for p in ket.config
                #bra = deepcopy(ket)
                copy!(bra,ket)

                apply_annihilation!(bra,p)
                if bra.sign == 0
                    continue
                end
                apply_creation!(bra,q)
                issorted(bra.config) || throw(Exception)
                if bra.sign == 0
                    continue
                else
                    calc_linear_index!(bra,binom_coeff)
                    #@code_warntype calc_linear_index!(bra,binom_coeff)
                    tbl[q, p, K] = bra.sign*bra.lin_index
                end
            end
        end
        incr!(ket)
    end
    return tbl
end
#=}}}=#

"""
    fill_ca_lookup3(c::DeterminantString)

Create an index table relating each DeterminantString with all ia substitutions
i.e., ca_lookup[Ka,p,q] = (sign,La)

<L|p'q|K> = sign
"""
function fill_ca_lookup3(c::DeterminantString)
    #={{{=#

    ket = DeterminantString(c.no, c.ne)
    bra = DeterminantString(c.no, c.ne)

    max = calc_max(ket)

    tbl = Array{Tuple,3}(undef,ket.no,ket.no,max)
    #tbl = zeros(Int,ket.no, ket.no, max)
    for K in 1:max
        for p in 1:ket.no
            for q in 1:ket.no
                bra = deepcopy(ket)
                apply_annihilation!(bra,p)
                apply_creation!(bra,q)
                #@assert(issorted(bra.config))
                if bra.sign == 0
                    tbl[q, p, K] = (0.0,0)
                else
                    calc_linear_index!(bra)
                    tbl[q, p, K] = (1.0*bra.sign,bra.lin_index)
                end
            end
        end
        incr!(ket)
    end
    return tbl
end
#=}}}=#

#"""
#    fill_vo_lookup(c::DeterminantString)
#
#Create an index table relating each DeterminantString with all ia substitutions
#i.e., ca_lookup[v,o,Ka] = sign*La
#where o and v are occupied and virtual indices
#
#<L|v'o|K> = sign
#"""
#function fill_vo_lookup(c::DeterminantString)
#    #={{{=#
#
#    ket = DeterminantString(c.no, c.ne)
#    bra = DeterminantString(c.no, c.ne)
#
#    max = calc_max(ket)
#
#    no = ket.no
#    ne = ket.ne
#    nv = no-ne
#
#    tbl = zeros(Int,nv, ne, max)
#    #tbl = Array{Tuple,3}(undef,nv,ne,max)
#    #tbl = Array{SVector{4,Int},nv, ne, max}
#    #println("max:",max)
#    for K in 1:max
#        virt = get_unoccupied(ket)
#        for pp in 1:ne
#            p = ket.config[pp]
#            bra1 = deepcopy(ket)
#            apply_annihilation!(bra1,p)
#            for qq in 1:nv
#                q = virt[qq]
#                bra = deepcopy(bra1)
#                apply_creation!(bra,q)
#                @assert(issorted(bra.config))
#                if bra.sign == 0
#                    continue
#                else
#                    calc_linear_index!(bra)
#                    tbl[qq, pp, K] = bra.sign*bra.lin_index
#                    #tbl[qq, pp, K] = (bra.sign, bra.lin_index, q, p)
#                    #println()
#                    #display(bra)
#                    #println("qq=", qq,", pp=", pp,", K=", K)
#                    #println("sgn=", bra.sign, ", lin=", bra.lin_index,", q=", q, ", p=", p)
#                end
#            end
#        end
#        incr!(ket)
#    end
#    #println(maximum(tbl))
#    #println(minimum(tbl))
#    #throw(Exception)
#    return tbl
#end
##=}}}=#

"""
    reset!(c::DeterminantString)

Reset the DeterminantString to the first config
"""
function reset!(c::DeterminantString)
    #={{{=#
    c.config = Vector(1:c.ne)
    c.sign = 1
    c.lin_index = 1
end
#=}}}=#

function destroy_config!(c::DeterminantString)
    #={{{=#
    c.config = []
    c.sign = 0
    c.lin_index = 0
    c.max = 0
    c.ne = 0
end
#=}}}=#

"""
    apply_annihilation!(c::DeterminantString, orb_index::Int)

Apply an annihilation operator to `c` corresponding to orbital `orb_index`
"""
function apply_annihilation!(c::DeterminantString, orb_index::Int)
    #=
    apply annihilation operator a_i to current DeterminantString
    where orb_index is i
    =#
    #={{{=#
    @assert(orb_index <= c.no)
    if c.sign == 0
        return
    end
    found = -1
    for i in 1:c.ne
        if c.config[i] == orb_index
            found = i
            break
        end
    end

    if found == -1
        destroy_config!(c)
        return
    end

    if found % 2 != 1
        c.sign *= -1
    end

    deleteat!(c.config,found)

    c.ne -= 1

    #unset data that need to be recomputed
    c.max = 0
    c.lin_index = 0
end
#=}}}=#

"""
    apply_creation!(c::DeterminantString, orb_index::Int)

Apply a creation operator to `c` corresponding to orbital `orb_index`
"""
function apply_creation!(c::DeterminantString, orb_index::Int)
    #=
    apply creation operator a_i to current DeterminantString
    where orb_index is i
    =#
    #={{{=#
    @assert(orb_index <= c.no)
    if c.sign == 0
        return
    end

    insert_here = 1
    for i in 1:c.ne
        if c.config[i] > orb_index
            insert_here = i
            break
        elseif c.config[i] == orb_index
            destroy_config!(c)
            return
        else
            insert_here += 1
        end
    end

    if insert_here % 2 != 1
        c.sign *= -1
    end

    #print("insert_here ", insert_here, ' ', c.ne, "\n")
    insert!(c.config,insert_here,orb_index)

    c.ne += 1
    #unset data that need to be recomputed
    c.max = 0
    c.lin_index = 0
end
#=}}}=#

function slater_det_energy(h0, h1, h2, deta::DeterminantString, detb::DeterminantString)
    #Compute the energy of a Slater Det specified by
    #   deta is alpha determinant
    #   detb is beta  determinant
    #={{{=#
    E0 = h0
    E1 = 0
    E2 = 0

    na = deta.config
    nb = detb.config
    for ii=1:deta.ne
        i = na[ii]
        E1 += h1[i,i]
        for jj=ii:deta.ne
            j = na[jj]
            E2 += h2[i,i,j,j]
            E2 -= h2[i,j,j,i]
        end
        for jj=1:detb.ne
            j = nb[jj]
            E2 += h2[i,i,j,j]
        end
    end
    for ii=1:detb.ne
        i = nb[ii]
        E1 += h1[i,i]
        for jj=i:detb.ne
            j = nb[jj]
            E2 += h2[i,i,j,j]
            E2 -= h2[i,j,j,i]
        end
    end
    return E0+E1+E2
end
#=}}}=#

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
# # Munk Gyre analytical solution equations
# things between #### blocks are the code for checking exact solution in DGmodel test

using Pkg
Pkg.add("Latexify")
using Latexify
Pkg.add("Plots")
using Plots
Pkg.add("PyPlot")
pyplot()

# Problem parameters
# These will come from DG model/problem settings
###
τ₀=0.1      # Stress
ρ=1000      # density
f=1.e-4     # Coriolis
β=1.e-11    # Rate of change of Coriolis
H=1000      # Depth
L_x=1000.e3 # Zonal extent
L_y=1000.e3 # Meridional extent
Ah=1000.    # Viscosity
g=10.       # Gravity

δ_m=(Ah/β)^(1/3)
####

# Set up Munk formula per https://mitgcm.readthedocs.io/en/latest/examples/examples.html#model-solution
# Build up using expressions to catch typos
t1=:( cos((3^0.5*x)/(2*δ_m))+(1)/(3^0.5)*sin((3^0.5*x)/(2*δ_m)) )
display(latexify(t1))
t2=:( 1-exp((-x)/(2*δ_m) ) * $t1 )
display(latexify(t2))
t3=:( π*sin(π*y/L_y) )
display(latexify(t3))
t4=:( τ₀/(ρ*g*H)*(f/β)*(1-x/L_x) )
display(latexify(t4))
tfull=:( ($t4) * ($t3) * ($t2) )
display(latexify(tfull))
display(tfull)

# Actual function (a copy paste from display(tfull) line)
#### This is formula needed for checking numerical against exact
ηfun(x,y)=(((τ₀ / (ρ * g * H)) * (f / β) * (1 - x / L_x)) * (π * sin((π * y) / L_y)) * (1 - exp(-x / (2δ_m)) * (cos((3 ^ 0.5 * x) / (2δ_m)) + (1 / 3 ^ 0.5) * sin((3 ^ 0.5 * x) / (2δ_m)))))
####

# Now make a plot to check
xvals=collect(0:L_x/100:L_x)
yvals=collect(0:L_y/100:L_y)
nx=length(xvals)
ny=length(yvals)
η=zeros(nx,ny)
for j=1:ny
 for i=1:nx
  η[i,j]=ηfun(xvals[i],yvals[j])
 end
end

p2 = contour(xvals, yvals, η')
plot(p2)

max_stars_repo_name: christophernhill/climahacks

max_issues_repo_name: christophernhill/climahacks

max_forks_repo_name: christophernhill/climahacks


text:
module _rXMSEn
export rXMSEn
using DSP.Filters: filtfilt, Butterworth, Lowpass, digitalfilter
using Statistics: std, mean, median, var
using Plots
    """
        MSx, CI = rXMSEn(Sig, Mobj)

    Returns a vector of refined multiscale cross-entropy values (`MSx`) and
    the complexity index (`CI`) between the data sequences contained in `Sig`
    using the parameters specified by the multiscale object (`Mobj`) and the
    following default parameters:   Scales = 3, Butterworth LPF Order = 6,
    Butterworth LPF cutoff frequency at scale (T): Fc = 0.5/T.
    If the entropy function specified by `Mobj` is `XSampEn` or `XApEn`,
    `rMSEn` updates the threshold radius of the data sequences (Xt) at each scale
    to 0.2std(Xt) when no `r` value is provided by `Mobj`, or r.std(Xt) if
    `r` is specified.

        MSx, CI = rXMSEn(Sig::AbstractArray{T,2} where T<:Real, Mobj::NamedTuple; Scales::Int=3,
                                F_Order::Int=6, F_Num::Float64=0.5, RadNew::Int=0, Plotx::Bool=false)

    Returns a vector of refined multiscale cross-entropy values (`MSx`) and
    the complexity index (`CI`) between the data sequences contained in `Sig`
    using the parameters specified by the multiscale object (`Mobj`) and the
    following keyword arguments:

    # Arguments:
    `Scales`   - Number of temporal scales, an integer > 1 (default: 3)  \n
    `F_Order`  - Butterworth low-pass filter order, a positive integer (default: 6)  \n
    `F_Num`    - Numerator of Butterworth low-pass filter cutoff frequency,
                 a scalar value in range [0 < `F_Num` < 1]. The cutoff frequency
                 at each scale (T) becomes: Fc = `F_Num`/T.  (default: 0.5) \n
    `RadNew`   - Radius rescaling method, an integer in the range [1 4].
                 When the entropy specified by `Mobj` is `XSampEn` or `XApEn`,
                 `RadNew` allows the radius threshold to be updated at each
                 time scale (Xt). If a radius value is specified by `Mobj` (`r`),
                 this becomes the rescaling coefficient, otherwise it is set
                 to 0.2 (default). The value of `RadNew` specifies one of the
                 following methods:  \n
                 [1]    Standard Deviation          - r*std(Xt)  (default)  \n
                 [2]    Variance                    - r*var(Xt)  \n
                 [3]    Mean Absolute Deviation     - r*mean_ad(Xt)  \n
                 [4]    Median Absolute Deviation   - r*med_ad(Xt,1) \n
    `Plotx`    - When `Plotx` == true, returns a plot of the entropy value at
                 each time scale (i.e. the multiscale entropy curve)
                 [default = false] \n

    # See also `MSobject`, `XMSEn`, `cXMSEn`, `hXMSEn`, `XSampEn`, `XApEn`, `MSEn`

    # References:
        [1]   Matthew W. Flood,
            "rXMSEn - EntropyHub Project"
            2021, https://github.com/MattWillFlood/EntropyHub

        [2]   Rui Yan, Zhuo Yang, and Tao Zhang,
            "Multiscale cross entropy: a novel algorithm for analyzing two
            time series."
            5th International Conference on Natural Computation.
            Vol. 1, pp: 411-413 IEEE, 2009.

        [3] José Fernando Valencia, et al.,
            "Refined multiscale entropy: Application to 24-h holter
            recordings of heart period variability in healthy and aortic
            stenosis subjects."
            IEEE Transactions on Biomedical Engineering
            56.9 (2009): 2202-2213.

        [4] Puneeta Marwaha and Ramesh Kumar Sunkaria,
            "Optimal selection of threshold value ‘r’for refined multiscale
            entropy."
            Cardiovascular engineering and technology
            6.4 (2015): 557-576.

        [5] Yi Yin, Pengjian Shang, and Guochen Feng,
            "Modified multiscale cross-sample entropy for complex time
            series."
            Applied Mathematics and Computation
            289 (2016): 98-110.

        [6] Antoine Jamin, et al,
            "A novel multiscale cross-entropy method applied to navigation
            data acquired with a bike simulator."
            41st annual international conference of the IEEE EMBC
            IEEE, 2019.

        [7] Antoine Jamin and Anne Humeau-Heurtier.
            "(Multiscale) Cross-Entropy Methods: A Review."
            Entropy
            22.1 (2020): 45.

    """
    function rXMSEn(Sig::AbstractArray{T,2} where T<:Real, Mobj::NamedTuple; Scales::Int=3,
        F_Order::Int=6, F_Num::Float64=0.5, RadNew::Int=0, Plotx::Bool=false)

    size(Sig,1) == 2  ?  Sig = Sig' : nothing

    (size(Sig,1)>10) ? nothing : error("Sig:   must be a numeric vector" )
    (length(Mobj) >= 1) ? nothing :  error("Mobj:    must be a multiscale entropy object created
         with the function EntropyHub.MSobject")
    (Scales>1) ? nothing : error("Scales:     must be an integer > 1")
    (F_Order>1) ? nothing :  error("F_Order:     must be an integer > 1")
    (0 < F_Num < 1) ? nothing : error("F_Num:     must be a scalar in range 0 < F_Num < 1")
    (RadNew==0 || (RadNew in 1:4 && String(Symbol(Mobj.Func)) in ("XSampEn","XApEn"))) ? nothing :
    error("RadNew:  must be 0, or an integer in range [1 4] with entropy function `XSampEn` or `XApEn`")

    MSx = zeros(Scales)
    Args = NamedTuple{keys(Mobj)[2:end]}(Mobj)
    (RadNew==0 && String(Symbol(Mobj.Func)) in ("XSampEn","XApEn")) ? RadNew=1 : nothing
    if RadNew > 0
        if RadNew == 1
            Rnew = x -> std(x, corrected=false)
        elseif RadNew == 2
            Rnew = x -> var(x, corrected=false)
        elseif RadNew == 3
            Rnew = x -> mean(abs.(x .- mean(x)))
        elseif RadNew == 4
            Rnew = x -> median(abs.(x .- median(x)))
        end

        if haskey(Mobj,:r)
            Cx = Mobj.r
        else
            Cy = ("Standard Deviation","Variance","Mean Abs Deviation",
                    "Median Abs Deviation")
            @warn("No radius value provided in Mobj.
                Default set to 0.2*$(Cy[RadNew]) of each new time-series.")
            Cx = .2
        end
    end

    for T = 1:Scales
        print(" .")
        Temp = refined(Sig,T,F_Order,F_Num)
        RadNew > 0 ? Args = (Args..., r=Cx*Rnew(Temp[:])) : nothing
        Tempx = Mobj.Func(Temp; Args...)
        typeof(Tempx)<:Tuple ? MSx[T] = Tempx[1][end] : MSx[T] = Tempx[end]
        #MSx[T] = Tempx[1][end]
    end
    CI = sum(MSx)
    print("\n")
    if any(isnan.(MSx))
        println("Some entropy values may be undefined.")
    end

    if Plotx
        p1 = plot(1:Scales, MSx, c=RGB(8/255, 63/255, 77/255), lw=3)
        scatter!(1:Scales, MSx, markersize=6, c=RGB(1, 0, 1),
        xlabel = "Scale Factor", ylabel = "Entropy Value",
        guidefont = font(12, "arial", RGB(7/255, 54/255, 66/255)),
        tickfontsize = 10, tickfontfamily="arial", legend=false,
        title = "Refined Multiscale $(Mobj.Func)",
        plot_titlefontsize=16, plot_titlefontcolor=RGB(7/255, 54/255, 66/255))
        display(p1)
    end

    return MSx, CI
    end

    function refined(Z,sx,P1,P2)
        Yt = filtfilt(digitalfilter(Lowpass(P2/sx), Butterworth(P1)), Z)
        return Yt[1:sx:end,:]
    end
end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
mutable struct Parser

    alloc::Bool
    ptr::Ptr{Void}

    function Parser(a, ptr)
        this = new(a, ptr)
        finalizer(this, dealloc)
        this
    end
end

Parser(ext::String) = Parser(true, ccall((:PARSER_new_for_file, libpfnet), Ptr{Void}, (Cstring,), ext))

parse_case(parser::Parser, filename::String; num_periods::Int=1) = Network(ccall((:PARSER_parse, libpfnet),
                                                                                Ptr{Void},
                                                                                (Ptr{Void}, Cstring, Int,),
                                                                                parser.ptr,
                                                                                filename,
                                                                                num_periods))

function dealloc(parser::Parser)
    if parser.alloc
        ccall((:PARSER_del, libpfnet), Void, (Ptr{Void},), parser.ptr)
    end
    parser.alloc = false
    parser.ptr = C_NULL
end


text:
function build_robot!(vis::Visualizer, model::PointFootQuadruped;
	body_height=0.025,
	body_length=0.17,
	body_width=0.15,
	foot_radius=0.025,
	color_opacity=1.0,
	name::Symbol=:robot)

	body_mat = MeshPhongMaterial(color=RGBA(0.0, 0.0, 0.0, color_opacity))
	foot_mat = MeshPhongMaterial(color=RGBA(1.0, 165.0 / 255.0, 0.0, color_opacity))

	default_background!(vis)

	setobject!(vis[name][:body],
    	GeometryBasics.Rect(Vec(-body_length, -body_width, -body_height),Vec(2.0 * body_length, 2.0 * body_width, body_height)),
    	body_mat)

	foot1 = setobject!(vis[name][:foot1], GeometryBasics.Sphere(Point3f0(0),
        convert(Float32, foot_radius)),
        foot_mat)

	foot2 = setobject!(vis[name][:foot2], GeometryBasics.Sphere(Point3f0(0),
		convert(Float32, foot_radius)),
		foot_mat)

	foot3 = setobject!(vis[name][:foot3], GeometryBasics.Sphere(Point3f0(0),
		convert(Float32, foot_radius)),
		foot_mat)

	foot4 = setobject!(vis[name][:foot4], GeometryBasics.Sphere(Point3f0(0),
		convert(Float32, foot_radius)),
		foot_mat)
end

function set_robot!(vis::Visualizer, model::PointFootQuadruped, q::AbstractVector;
	p_shift=[0.0; 0.0; 0.025], name::Symbol=:robot)

	R = mrp_rotation_matrix(q[4:6])

	p_body  = q[0  .+ (1:3)] + p_shift
	p_foot1 = q[6  .+ (1:3)] + p_shift
	p_foot2 = q[9  .+ (1:3)] + p_shift
	p_foot3 = q[12 .+ (1:3)] + p_shift
	p_foot4 = q[15 .+ (1:3)] + p_shift

	settransform!(vis[name][:body], compose(Translation(p_body), LinearMap(mrp_rotation_matrix(q[4:6]))))
	settransform!(vis[name][:foot1], Translation(p_foot1))
	settransform!(vis[name][:foot2], Translation(p_foot2))
	settransform!(vis[name][:foot3], Translation(p_foot3))
	settransform!(vis[name][:foot4], Translation(p_foot4))

end

function visualize!(vis, model::PointFootQuadruped, q;
	Δt=0.1,
	body_height=0.025,
	body_length=0.17,
	body_width=0.15,
	foot_radius=0.025,
	color_opacity=1.0,
	fixed_camera=false,
	anim=MeshCat.Animation(convert(Int, floor(1.0 / Δt))),
	name::Symbol=:robot)

	build_robot!(vis, model,
		body_height=body_height,
		body_length=body_length,
		body_width=body_width,
		foot_radius=foot_radius,
		color_opacity=color_opacity,
		name=name)

	for (t, qt) in enumerate(q)
		MeshCat.atframe(anim, t) do
			set_robot!(vis, model, qt, name=name)
		end
	end

	MeshCat.setanimation!(vis, anim)

	return anim
end


text:
using OptimizeQCQP
using JuMP, Ipopt   # for comparison purposes
using LinearAlgebra
using Test

function checkposdef(M::SymTridiagonal)
    dv, ev = M.dv, M.ev
    n = length(dv)
    α = 4*(cos(π/(n+1)))^2
    for i = 1:n-1
        @test α*ev[i]^2 <= (1+1e-8)*dv[i]*dv[i+1]
    end
    return nothing
end
function checkposdef(M)
    @assert isposdef(M)
    return nothing
end

function checksolution(Q0, g0, Q1, cu; print_level=0)
    solver = IpoptSolver(print_level=print_level,
                         sb="yes")
                         #  tol=tol,
                         #  max_iter=max_iter)
    model = Model(solver=solver)
    n = length(g0)
    @variable(model, x[1:n])
    @constraint(model, x'*Q1*x <= 2*cu)
    @objective(model, Min, (x'*Q0*x)/2 + g0'*x)
    status = solve(model)
    status == :Optimal || error("checksolution solution not optimal")
    y = getvalue(x)
    # Ipopt sometimes has relatively low accuracy, so make sure we enforce the constraint
    ymag = (y'*Q1*y)/2
    if ymag > cu
       y .*= sqrt(cu/ymag)
    end
    return y, (y'*Q0*y)/2 + g0'*y
end

function symtridiag(n, isposdef)
    if isposdef
        dv = rand(n) + 2
        ev = -rand(n-1)
        return SymTridiagonal(dv, ev)
    else
        dv = randn(n)
        ev = randn(n-1)
        return SymTridiagonal(dv, ev)
    end
end

function symdense(eigvals)
    n = length(eigvals)
    U, _ = qr(randn(n, n))
    Q = U*Diagonal(eigvals)*U'
    return (Q + Q')/2
end

@testset "posdefλ" begin
    A = symdense([-1.4, 0.2, 3])
    λ, _ = OptimizeQCQP.posdefλ(A, I)
    @test λ >= 1.4 && λ < 2
    λ, _ = OptimizeQCQP.posdefλ(A, 2*I)
    @test λ >= 0.7 && λ < 1
    λ, _ = OptimizeQCQP.posdefλ(A, 0.5*I)
    @test λ >= 2.8 && λ < 4
    B = Diagonal(fill(0.5, 3))
    λ, _ = OptimizeQCQP.posdefλ(A, B)
    @test λ >= 2.8 && λ < 4
    λ, _ = OptimizeQCQP.posdefλ(A, Matrix(B))
    @test λ >= 2.8 && λ < 4

    A = SymTridiagonal([-1, 0.1, 2], [0.0, 0.0])
    λ, _ = OptimizeQCQP.posdefλ(A, I; growauto=1+1e-9)
    @test λ ≈ 1
    λ, _ = OptimizeQCQP.posdefλ(A, 2*I; growauto=1+1e-9)
    @test λ ≈ 0.5
    A = SymTridiagonal([-1, 0.1, 2], [1.0, 1.0])
    λ, _ = OptimizeQCQP.posdefλ(A, I)
    @test λ >= -minimum(eigvals(Matrix(A))[1])
end

@testset "Q0/g0/Q1" begin
    n = 4
    # With positive-definite Q0
    for gen in (()->symdense([0.1, 1, 2, 3]),
                ()->symtridiag(4, true),)
        for i = 1:5
            Q0 = gen()
            g0 = randn(n)
            xNewton = Q0 \ (-g0)
            Q1 = 1.0*I
            @test OptimizeQCQP.posdefλ(Q0, Q1)[1] == 0
            cu = 1.1*(xNewton'*Q1*xNewton)/2
            x, λ = minimize(Q0, g0, nothing, nothing, Q1, nothing, nothing, cu; tol=1e-8)
            @test (x'*Q1*x)/2 <= 1.01*cu
            @test λ == 0
            x_ipopt, val_ipopt = checksolution(Q0, g0, Q1, cu)
            @test (x'*Q0*x)/2 + g0'*x <= val_ipopt + 1e-5*abs(val_ipopt)
            cu = 0.3*(xNewton'*Q1*xNewton)/2
            x, λ = minimize(Q0, g0, nothing, nothing, Q1, nothing, nothing, cu; tol=1e-8)
            @test (x'*Q1*x)/2 <= 1.01*cu
            @test λ > 0
            x_ipopt, val_ipopt = checksolution(Q0, g0, Q1, cu)
            @test (x'*Q0*x)/2 + g0'*x <= val_ipopt + 1e-5*abs(val_ipopt)
            Q1 = Diagonal(0.25+rand(n))
            @test OptimizeQCQP.posdefλ(Q0, Q1)[1] == 0
            cu = 1.1*(xNewton'*Q1*xNewton)/2
            x, λ = minimize(Q0, g0, nothing, nothing, Q1, nothing, nothing, cu; tol=1e-8)
            @test (x'*Q1*x)/2 <= 1.01*cu
            @test λ == 0
            x_ipopt, val_ipopt = checksolution(Q0, g0, Q1, cu)
            @test (x'*Q0*x)/2 + g0'*x <= val_ipopt + 1e-5*abs(val_ipopt)
            cu = 0.3*(xNewton'*Q1*xNewton)/2
            x, λ = minimize(Q0, g0, nothing, nothing, Q1, nothing, nothing, cu; tol=1e-8)
            @test (x'*Q1*x)/2 <= 1.01*cu
            @test λ > 0
            x_ipopt, val_ipopt = checksolution(Q0, g0, Q1, cu)
            @test (x'*Q0*x)/2 + g0'*x <= val_ipopt + 1e-5*abs(val_ipopt)
        end
    end
    # With indefinite Q0
    for gen in (()->symdense([-1, -0.1, 2, 3]),
                ()->symtridiag(4, false),)
        for i = 1:5
            Q0 = gen()
            g0 = randn(n)
            for cu in (1.3, 0.5, 0.1, 0.01, 1e-6)
                Q1 = 1.0*I
                λ, _ = OptimizeQCQP.posdefλ(Q0, Q1)
                checkposdef(Q0 + λ*Q1)
                x, λ = minimize(Q0, g0, nothing, nothing, Q1, nothing, nothing, cu; tol=1e-8)
                @test (x'*Q1*x)/2 <= 1.01*cu
                x_ipopt, val_ipopt = checksolution(Q0, g0, Q1, cu)
                @test (x'*Q0*x)/2 + g0'*x <= val_ipopt + 1e-5*abs(val_ipopt)
                Q1 = Diagonal(0.25+rand(n))
                λ, _ = OptimizeQCQP.posdefλ(Q0, Q1)
                checkposdef(Q0 + λ*Q1)
                x, λ = minimize(Q0, g0, nothing, nothing, Q1, nothing, nothing, cu; tol=1e-8)
                @test (x'*Q1*x)/2 <= 1.01*cu
                x_ipopt, val_ipopt = checksolution(Q0, g0, Q1, cu)
                @test (x'*Q0*x)/2 + g0'*x <= val_ipopt + 1e-5*abs(val_ipopt)
            end
        end
    end
end


text:
for (elty, ext) in ((:Float32, :s),
                    (:Float64, :d),
                    (:ComplexF32, :c),
                    (:ComplexF64, :z))

    # Distributed sparse gemv
    @eval begin
        function LinearAlgebra.mul!(y::DistMultiVec{$elty}, A::DistSparseMatrix{$elty}, x::DistMultiVec{$elty}, α::$elty, β::$elty)
            ElError(ccall(($(string("ElMultiplyDist_", ext)), libEl), Cuint,
                (Cint, $elty, Ptr{Cvoid}, Ptr{Cvoid}, $elty, Ptr{Cvoid}),
                NORMAL, α, A.obj, x.obj, β, y.obj))
            return y
        end
        function LinearAlgebra.mul!(y::DistMultiVec{$elty}, adjA::Adjoint{<:Any,DistSparseMatrix{$elty}}, x::DistMultiVec{$elty}, α::$elty, β::$elty)
            ElError(ccall(($(string("ElMultiplyDist_", ext)), libEl), Cuint,
                (Cint, $elty, Ptr{Cvoid}, Ptr{Cvoid}, $elty, Ptr{Cvoid}),
                ADJOINT, α, parent(adjA).obj, x.obj, β, y.obj))
            return y
        end
        function LinearAlgebra.mul!(y::DistMultiVec{$elty}, trA::Transpose{<:Any,DistSparseMatrix{$elty}}, x::DistMultiVec{$elty}, α::$elty, β::$elty)
            ElError(ccall(($(string("ElMultiplyDist_", ext)), libEl), Cuint,
                (Cint, $elty, Ptr{Cvoid}, Ptr{Cvoid}, $elty, Ptr{Cvoid}),
                TRANSPOSE, α, parent(trA).obj, x.obj, β, y.obj))
            return y
        end
    end
end


text:
# # Subspace angles
# This example considers the angles between neighbouring Laguerre polynomials with a perturbed measure:
# ```math
# \cos\theta_n = \frac{\langle L_n, L_{n+k}\rangle}{\|L_n\|_2 \|L_{n+k}\|_2},\quad{\rm for}\quad 0\le n < N-k,
# ```
# where the inner product is defined by $\langle f, g\rangle = \int_0^\infty f(x) g(x) x^\beta e^{-x}{\rm\,d}x$.
#
# We do so by connecting Laguerre polynomials to the normalized generalized Laguerre polynomials associated with the perturbed measure. It follows by the inner product of the connection coefficients that:
# ```math
# \cos\theta_n = \frac{(V^\top V)_{n, n+k}}{\sqrt{(V^\top V)_{n, n}(V^\top V)_{n+k, n+k}}}.
# ```
#
using FastTransforms, LinearAlgebra

# The neighbouring index `k` and the maximum degree `N-1`:
k, N = 1, 11

# The Laguerre connection parameters:
α, β = 0.0, 0.125

# We precompute a Laguerre--Laguerre plan:
P = plan_lag2lag(Float64, N, α, β; norm2=true)

# We apply the plan to the identity, followed by the adjoint plan:
VtV = P*I
lmul!(P', VtV)

# From this matrix, the angles are recovered from:
θ = [acos(VtV[n, n+k]/sqrt(VtV[n, n]*VtV[n+k, n+k])) for n in 1:N-k]


text:
abstract type AbstractPooledArray{T,N} <: AbstractArray{T,N}
end

mutable struct PooledAbstractArray{T,N} <: AbstractPooledArray{T,N}
    array::AbstractArray{T,N}
    ref_count::Int
    dispose::Function
    function PooledAbstractArray(a::AbstractArray{T,N}, dispose::Function) where {T,N}
        new{T,N}(a, 1, dispose)
    end
end

resource(a::AbstractPooledArray{T,N}) where {T,N} = a.array

Base.size(a::AbstractPooledArray{T,N}) where {T,N} = size(a.array)

Base.IndexStyle(::Type{<:AbstractPooledArray{T,N}}) where {T,N} = IndexCartesian()

Base.getindex(a::AbstractPooledArray{T,N}, indices::Vararg{Int, N}) where {T,N} = a.array[indices...]

Base.setindex!(a::AbstractPooledArray{T,N}, v, indices::Vararg{Int, N}) where {T,N} = a.array[indices...] = v


text:
# perturb an array m by percentage of m0, when mask is true
function perturb!(m, m0, perc, mask::Function)
    for i in CartesianIndices(m)
        ii = Tuple(i)
        if (mask(ii))
            m[i] = m[i] + perc * 1e-2 * m0
        end
    end
end
function perturb!(m, m0, perc, idx::CartesianIndices)
    for i in idx
        m[i] = m[i] + perc * 1e-2 * m0
    end
end

"""
In-place method to add a rectangular box or cuboid in the medium. Its position is specified by its corners
and perc denotes the percentage of the perturbation.

```julia
update!(medium, field; rectangle=[[0,10],[10,10]], perc=10) # for 2D medium
update!(medium, field; rectangle=[[0,10,15],[10,10,15]], perc=10) # for 3D medium
```

In-place method to add random noise to medium.

```julia
update!(medium, fields; randn_perc=10)
```
"""
function update!(
    medium::Medium,
    fields::Vector{Symbol};
    rectangle = nothing,
    perc = 0.0,
    randn_perc = 0.0,
)

    mgrid = medium.mgrid
    if (!(rectangle === nothing))
        @assert length(rectangle)==2
        @assert length.(rectangle)==fill(ndims(medium),2)
        rectangle_indices = CartesianIndices(
            Tuple([
                argmin(abs.(mg .- rectangle[1][im])):argmin(abs.(mg .- rectangle[2][im])) for (im, mg) in enumerate(mgrid)
            ]),
        )
        if (!iszero(perc))
            for field in fields
                m = medium.m[field]
                m0 = medium.ref[field]
                perturb!(m, m0, perc, rectangle_indices)
            end
        end
    end
    if (!iszero(randn_perc))
        for field in fields
            m = medium.m[field]
            m0 = medium.ref[field]
            for i in eachindex(m)
                m[i] = χ(Float64(χ(m[i], m0, 1) + randn() * randn_perc * 1e-2), m0, -1)
            end

        end
    end

    return medium
end

# # Keyword Arguments

# * `point_loc::Vector{Float64}=[0., 0.,]` : approx location of point pert
# * `point_pert::Float64=0.0` : perturbation at the point scatterer
# * `ellip_loc::Vector{Float64}=nothing` : location of center of perturbation, [z, x]
# * `ellip_rad::Float64=0.0` : size of elliptic perturbation
# * `ellip_pert::Float64=0.1` : perturbation inside the ellipse
# * `ellip_α=0.0` : rotate the ellipse
# * `rect_loc::Array{Float64}=nothing` : rectangle location, [zmin, xmin, zmax, xmax]
# * `rect_pert::Float64=0.1` : perturbation in a rectangle
# * `constant_pert::Float64=0.0` : constant perturbation
# * `randn_pert::Float64=0.0` : percentage of reference values for additive random noise
# * `fields::Vector{Symbol}=[:χvp,:χrho,:χvs]` : which fields are to be modified?
# * `onlyin` : `medium` is modified only when field values are in these ranges
# """
# function update!(
#     medium::Medium,
#     fields::Vector{Symbol};
#     point_loc = fill(0.0, ndims(medium)),
#     point_pert::Float64 = 0.0,
#     ellip_loc = fill(0.0, ndims(medium)),
#     ellip_rad = 0.0,
#     ellip_pert::Real = 0.0,
#     ellip_α = 0.0,
#     rect_loc = fill(0.0, ndims(medium)),
#     rect_pert::Float64 = 0.0,
#     constant_pert::Float64 = 0.0,
#     randn_perc::Real = 0.0,
#     layerlocations = nothing,
# )

#     # only editing basic fields is allowed
#     for field in fields
#         @assert field ∈ [:vp, :vs, :rho, :Q]
#     end
#     rect_loc = convert.(Float64, rect_loc)
#     ellip_loc = convert.(Float64, ellip_loc)

#     temp = zeros(length.(medium.mgrid)...)

#     ipointloc = [
#         Interpolation.indminn(medium.mgrid[i], Float64(point_loc[i]), 1)[1] for
#         i = 1:length(medium.mgrid)
#     ]
#     temp[ipointloc...] += point_pert

#     if (!(ellip_pert == 0.0))
#         α = ellip_α * pi / 180.0
#         # circle or ellipse
#         rads =
#             (length(ellip_rad) == 1) ? [ellip_rad[1], ellip_rad[1]] :
#             [ellip_rad[1], ellip_rad[2]]

#         temp += [
#             (
#                 (
#                     (
#                         (
#                             (medium.mgrid[2][ix] - ellip_loc[2]) * cos(α) +
#                             (medium.mgrid[1][iz] - ellip_loc[1]) * sin(α)
#                         )^2 * inv(rads[1]^2) +
#                         (
#                             (-medium.mgrid[1][iz] + ellip_loc[1]) * cos(α) +
#                             (medium.mgrid[2][ix] - ellip_loc[2]) * sin(α)
#                         )^2 * inv(rads[2]^2)
#                     ) <= 1.0
#                 ) ? Float64(ellip_pert) : 0.0
#             ) for iz = 1:length(medium.mgrid[1]), ix = 1:length(medium.mgrid[2])
#         ]
#     end
#     if (!(rect_pert == 0.0))
#         temp += [
#             (
#                 (
#                     (medium.mgrid[2][ix] - rect_loc[4]) * (medium.mgrid[2][ix] - rect_loc[2]) < 0.0
#                 ) & (
#                     (medium.mgrid[1][iz] - rect_loc[3]) * (medium.mgrid[1][iz] - rect_loc[1]) < 0.0
#                 )
#             ) ? rect_pert : 0.0 for iz = 1:length(medium.mgrid[1]),
#             ix = 1:length(medium.mgrid[2])
#         ]
#     end
#     if (!(constant_pert == 0.0))
#         temp .+= constant_pert
#     end

#     for (iff, field) in enumerate(fields)
#         m = medium.m[field]
#         if (!(layerlocations === nothing))
#             onlyatvalues = []
#             for ipos in layerlocations
#                 ipx = Interpolation.indminn(medium.mgrid[2], Float64(ipos[2]), 1)[1]
#                 ipz = Interpolation.indminn(medium.mgrid[1], Float64(ipos[1]), 1)[1]
#                 push!(onlyatvalues, m[ipz, ipx])
#             end
#         end
#         for i in eachindex(m)
#             # if(((m[i]-onlyin[iff][1])*(m[i]-onlyin[iff][2]))<0.0)
#             if (layerlocations === nothing)
#                 m[i] += temp[i]
#             else
#                 if (m[i] ∈ onlyatvalues)
#                     m[i] += temp[i]
#                 end
#             end
#         end
#     end

#     return nothing
# end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ a6261f66-40ca-11ec-3021-ed2173550635
begin
	using Pkg
	Pkg.activate("../../")
	using Revise
end

# ╔═╡ 77e6389e-3491-45b6-9f76-459b15a8922d
using TestImages, ImgProcMic, ImageShow, Random, PlutoTest, Plots, FourierTools, Statistics

# ╔═╡ ee4b6b2d-119e-4ac3-a1f5-500ff045bfc5
using FFTW, PlutoUI

# ╔═╡ c7150c2c-afeb-4158-b381-3341bb678fd5
using Noise

# ╔═╡ 7148e705-8b90-46ea-8763-778897daf9a7
md"# Homework 03

In this homework we cover several topics like image correction, Fourier transforms, Fourier Shifting and Sampling.
"

# ╔═╡ 83d9c252-2a8d-47a3-b247-a62b6a175201
md"## 1. Flat Fielding
As discussed in the lecture, raw images might be affected by dirty optics or sensor errors.
In this exercise we want to apply such a correction to simulated images.
"

# ╔═╡ 39c623c7-6063-4360-a2b7-95ea42c5262d
begin
	img = Float32.(testimage("mandril_gray.tif"))
	gray_show(img)
end

# ╔═╡ eacd8fdf-971c-46ca-914f-f73f65e28ca8
md"### Simulate Deteriorated Image
In the first step, we degrade our image with some artifacts.
We add an offset to our image which is caused by sensor artifacts.
Second, we add some dirty spots to the image which are for example caused by dust on the sensor.
"

# ╔═╡ 25dde5cb-9ed6-498e-b511-50ed0640d0c0
"""
	simulate_bad_microscope(img)

The ideal `img` is degraded by artifacts.
"""
function simulate_bad_microscope(img::AbstractArray{T}) where T
	# generate same random numbers!
	rng = MersenneTwister(2021);
	# some weird dark field slope offset
	dark = zero.(img)
	dark .+= range(0,0.3, length=size(img, 1))

	xs = 1:size(img, 1)

	ff = one.(img)
	for i = 1:20 # emulates some dust spots as modulations to the flatfield ff
		x, y = rand(rng, xs), rand(rng, xs)
		sca = abs(randn(rng) .* T(0.2))
		r = sqrt.((xs .- x).^2 .+ (xs' .- y).^2)
		ff .-= T(0.2) .* sinc.(r .* sca)
	end

	return abs.(dark .+ ff .* img)
end

# ╔═╡ a6c024a4-7cb2-4f8d-be24-8727d75679da
img_dirty = simulate_bad_microscope(img);

# ╔═╡ d5c127ee-1ecf-4ed0-9cd9-6fcdefa95847
gray_show(img_dirty)

# ╔═╡ 6e95cbcb-0d89-4eca-8599-b266076708a7
md"## 1.1 Task
Extract both a flat field and a dark field image with `simulate_bad_microscope`.
Use the function `simulate_bad_microscope` to obtain a measured flat field and measured dark field.
"

# ╔═╡ 910cb7fa-b1fc-4106-acc4-3bf56edfb83d
flat_field = simulate_bad_microscope(ones(eltype(img), size(img)));

# ╔═╡ cbebb311-c521-4fa6-a909-1a4df8da4d2d
gray_show(flat_field)

# ╔═╡ 8bb737cb-2d59-49ba-9b89-afef6ca7bbf8
dark_field = simulate_bad_microscope(zeros(eltype(img), size(img)));

# ╔═╡ 52bc8b3d-542e-4a35-b38c-35840f4075d4
gray_show(dark_field)

# ╔═╡ e3448ffe-5fbf-4bed-bf03-6eb9031e6c12
"""
	correct_image(img_dirty, dark_field, flat_field)

Correct a dirty image `img_dirty` by providing the `dark_field` and `flat_field`.

Conceptually we divide by the `flat_field` since the measured `flat_field`
describes the deviations from an ideal white image.
The `dark_field` has to be subtracted before since that only contains sensor artifacts.
"""
function correct_image(img_dirty, dark_field, flat_field)
	return (img_dirty .- dark_field) ./ (flat_field .- dark_field)
end

# ╔═╡ d9308666-efaa-4fcc-a083-f3d12dbab1e2
img_fixed = correct_image(img_dirty, dark_field, flat_field);

# ╔═╡ 493b6bb4-1dfc-4a82-9360-40944b29826c
gray_show(img_fixed)

# ╔═╡ c80a2ce8-6620-4de1-bb78-5af97caf466d
md"### 1.1 Test"

# ╔═╡ e2e470e4-ff73-4dd5-973c-14e1535b01ea
PlutoTest.@test img_fixed ≈ img

# ╔═╡ 0a957f4a-a975-4087-b458-3236940eda43
PlutoTest.@test img ≈ correct_image(img, zero.(img), one.(img))

# ╔═╡ af9866f4-2b03-4dc4-9f8a-e02ae46e2f9c
md"# 2 Sampling
If a signal is not properly sampled, we can see a lot of artifacts which are called aliasing.
Try to determine experimentally (via the slider) which is the correct value for Nyquist sampling.

Afterwards, calculate the correct sampling in the variable $N$!

Don't state the pure number but rather a short equation which Julia calculates (something like `1 * 2 * 3`).
Do you see a difference between your expected and your calculated value?
"

# ╔═╡ 7a481af3-4ec8-4a58-bb84-e9d3d5381570
md"N=$(@bind N Slider(2:300, show_value=true))"

# ╔═╡ 4006ccf7-df8e-440c-8b3c-af525bb50ed4
begin
	xs = range(0, 8, length=N+1)[begin:end-1]
	f(x) =  sin(x * 2π * 5)
end

# ╔═╡ c871cdbd-b560-46f2-918c-45d24d2f9b75
plot(xs, f.(xs))

# ╔═╡ f6b8cd92-c1ab-445b-8b4f-04f521b72cad
#           for Nyquist    because of the        because of the interval length
N_correct = 1 +    2 *         5                    *   8

# ╔═╡ 51b2d7be-fe94-44e9-8e67-73aafb636ef1
begin
	xs_correct = range(0, 8, length=round(Int, N_correct + 1))[begin:end-1]
	plot(xs_correct, f.(xs_correct))
end

# ╔═╡ a19c8850-9a7d-4e26-908b-a067a9006e84
md"No tests, since they would reveal the answer."

# ╔═╡ b606e4d0-8cac-4028-8e38-113050531720
plot(resample(f.(xs_correct), (400 ,)))

# ╔═╡ b0c67bd6-b5ef-4359-8727-1a3bfefcf759
md"# 3 FFTs

Another very important topics are FFTs!
The standard library in many languages is based on FFTW.
In Julia FFTW.jl is the recommended library.

The discrete Fourier transform is given by

$$X_k = \sum_{n=0}^{N-1} x_n \exp\left(-i2 \pi \frac{kn}{N} \right)$$
and the inverse by

$$x_n = \frac1{N}\sum_{k=0}^{N-1} X_k \exp\left(i2 \pi \frac{kn}{N} \right)$$.

The FFT algorithm evaluates those equations very efficiently ($\mathcal O(N \cdot \log N)$) by avoiding the explicit sum ($\mathcal O(N^2)$).

Further, the FFT calculates the frequencies in an (maybe) unexpected format
* Data: $[x_1, x_2, ..., x_N]$
* FFT of Data: $[f_0, f_1, ..., f_{N/2}, f_{-N/2}, f_{-N/2+1}, ...,f_{-1}]$
* (for real data negative and positive frequencies are the complex conjugate)

Sometimes we want to have the frequencies centered, which means a `fftshift` operations is needed
* FFT of Data followed by fftshift: $[f_{-N/2}, f_{-N/2+1},..., f_{-1}, f_0, f_1,..., f_{N/2-1}, f_{N/2}]$
"

# ╔═╡ 7276f9e0-80b8-4794-9125-fca43246f818
"""
    my_dft(x::Vector{T})::Vector{Complex{T}}

Calculate the DFT of `data` according to the equation given above.
Note, the output should be a complex Vector!
"""
function my_dft(x::Vector{T})::Vector{Complex{T}} where T
	X = zeros(Complex{T}, size(x))

	# length of x
	N = size(x, 1)

	# loop over output elements
	for k ∈ 0:N-1
		# loop over inputs
		for n ∈ 0:N-1
			X[k + 1] += x[n + 1] * exp(-im * T(2π) * k * n / N)
		end
	end

	return X
end

# ╔═╡ f054e266-e16c-4311-9e5e-1dbb9d024a63
typeof(π)

# ╔═╡ e895cff6-4018-495f-9a0e-1e56ddcfd17d
arr = [1f0,2,3,4]

# ╔═╡ 14d68904-4de5-4f1d-81f6-dae8bf50b749
my_dft(arr)

# ╔═╡ 7307dd33-cb79-48c0-9e73-78ed70e9c628
fft(arr)

# ╔═╡ 52061097-5be4-4b1a-a4f6-fd030a5a52ba
arr123 = randn((1_000))

# ╔═╡ 66cdd72e-ec6f-49ec-8338-c32e744cc359
@time fft(arr123)

# ╔═╡ 3c77082d-20e9-4541-8426-f8457cc728ec
@time my_dft(arr123)

# ╔═╡ 31b60200-4d70-4a3c-995f-949b19800549
# From worker 4:       0.000315 seconds (32 allocations: 33.844 KiB)
# From worker 4:       0.052112 seconds (1 allocation: 15.750 KiB)

# ╔═╡ 393d03e2-461d-480a-9719-f33a4eef3a6f
md"### 3.1 Test"

# ╔═╡ df12fd86-f29d-4926-b360-daf4cb40af86
begin
	arr_r = randn((13,))
	PlutoTest.@test fft(arr_r) ≈ my_dft(arr_r)
end

# ╔═╡ 37eedc2e-27bd-4d1c-a83b-6ec14aba55ab
begin
	arr_r2 = randn((14,))
	PlutoTest.@test fft(arr_r2) ≈ my_dft(arr_r2)
end

# ╔═╡ 5313dbed-bfd7-4bb8-b39c-89c55e50d1bb
md"## 3.2 Task
Calculate the FFT of the following dataset, such that center frequency is in the center!

Use `fft` and `fftshift`.
"

# ╔═╡ e1eac3cd-e6ce-4743-9119-5ee0c69b8e00
begin
	data_odd = [1.7142857142857142, 1.643544850921574, 1.0471849694228637, 1.383330948397344, 1.8300917955219322, 1.2951185346328968, 1.086443186817675]
	data_even = [-2.5, 1.5, -2.5, 3.5]
end

# ╔═╡ 481961ac-ad7c-4768-a7aa-eed34143226b
ffts(x) = fftshift(fft(x))

# ╔═╡ 25a2de2a-612b-4ed3-8aad-d5caf5d59c30
[1,0,0,0,0]

# ╔═╡ 2a74167b-7056-47e1-861b-bf7c1810a1f0
fftshift([1,0,0,0,0])

# ╔═╡ c95f9463-de4d-47d4-abfd-f707641d2021
[1,0,0,0]

# ╔═╡ 6fc3a300-7b96-4b21-8487-fdf4b5eb7521
fftshift([1,0,0,0])

# ╔═╡ c24b2457-4396-48c6-a21f-6d982d80776e
length([1,2,3,4,5]) ÷ 2 + 1

# ╔═╡ 27aa0d27-1cb1-4827-80c7-48741cc0312e
length([1,2,3,4]) ÷ 2 + 1  # \div or div(N, 2)

# ╔═╡ fc3bb52b-f9a2-42c0-b34b-347d0a137e2d
ffts(data_even)

# ╔═╡ a55851b6-3f8d-4af3-8f09-a3777c56df1b
ffts(data_odd)

# ╔═╡ e46d67be-d974-4f07-8b3e-02fff3e749d0
md"## 3.3 Task

Now do the inverse operation such that `data ≈ iffts(ffts(x))`

Check out `ifftshift`.
"

# ╔═╡ deb9eaf0-ee89-49c5-90b8-570d1a88f8fc
iffts(x) = ifft(ifftshift(x))

# ╔═╡ 90dd14ae-fb42-400e-86b4-9ae9e63716fc

# ╔═╡ 34e1c198-abe0-470e-a4bc-348ec09a2032
ifftshift([0,0,1,0])

# ╔═╡ ffb5665d-ecac-4474-b2a4-9685cc40ad7e
ifftshift([0,0,1,0,0])

# ╔═╡ 6ca02109-6114-411f-bfb3-bd638ab93218

# ╔═╡ 2875999a-a61d-491a-8a42-0552a87d3c6a
md"""
### 3.2 and 3.3 Test
"""

# ╔═╡ 333a253e-ce8b-477a-9db2-f99c034ae048
PlutoTest.@test real(iffts(ffts(data_even))) ≈ data_even

# ╔═╡ e81990c8-2194-472c-a102-b0d010f3a266
 PlutoTest.@test real(ffts(iffts(data_odd))) ≈ data_odd

# ╔═╡ f3fc2ea5-033c-4789-8ce6-0c6b15b607b7
 PlutoTest.@test real(iffts(ffts(data_odd))) ≈ data_odd

# ╔═╡ 90596b8a-53d4-457f-9e93-522e6c5dc35d
 PlutoTest.@test real(iffts(ffts(data_even))) ≈ data_even

# ╔═╡ 4e54ddc6-b9cd-472f-b496-18055ee5b667
md"## 3.4 Task
Try to reproduce this $\cos$ curve.
The reason it looks so _sharp_ is that we only sample slightly above the Nyquist limit.

However, you are only allowed to use `fft`, `fftshift`, `ifftshift`, `ifft` and array manipulation (like `x[34] = 1723.12`).
And you can also take the `real` part at the end.

Please explain with some #comments what you did.
"

# ╔═╡ 9844ec03-b169-40a0-ad56-1c07b11bb886
begin
	x = range(0, 4, length=21)[begin:end-1]
	y = cos.(x .* 2π)
	plot(x, y, mark="-*")
end

# ╔═╡ 49c39cea-895e-429a-909c-d1438d0dfd00
function reproduce_cos(y)
	out = zeros(size(y))
	# we go to Fourier space with center frequncy in the center
	out_ft = ffts(out)

	out_ft[15] = 10
	out_ft[7] = 10

	# we go real space again
	out = real(iffts(out_ft))
	return out
end

# ╔═╡ fcdb96af-30bd-42c2-bf45-ecd13f27b57d
collect(fftshift(FFTW.fftfreq(length(y))))

# ╔═╡ 5c4f40df-cde3-4140-ace2-6c0e1b60ebcf
begin
	plot(x, y, mark="-*")
	plot!(x, reproduce_cos(y) .+  0.1, mark="-*")
end

# ╔═╡ 79578e4d-7527-459b-8915-ee5cb9ccad24
md"### 3.4 Test"

# ╔═╡ a74c4034-fa19-4360-9f5e-191d5ba4001a
PlutoTest.@test real(reproduce_cos(copy(y))) ≈ y

# ╔═╡ 8ce249a3-0881-49e5-8ffd-5df21bda50f4
md"## 3.5 Fourier Shift Theorem

The discrete Fourier transform is given by

$$X_k = \sum_{n=0}^{N-1} x_n \exp\left(-i2 \pi \frac{kn}{N} \right)$$
and the inverse by

$$x_n = \frac1{N}\sum_{k=0}^{N-1} X_k \exp\left(i2 \pi \frac{kn}{N} \right)$$.
"

# ╔═╡ 13131381-11b3-4db4-9351-fb7264c94ea4
md"## 3.4 Task
Proof the Fourier shift theorem (with LaTeX) which states

$x_{n+\Delta n} = \mathcal{F}^{-1}\left[ X_{k} \cdot f(k, \Delta n)\right]$

Derive how the function $f(k, \Delta n)$ looks like!

hint: by clicking on the eye symbol next to this task, you can see how to embed LaTeX code in markdown comments.
"

# ╔═╡ bc3ae07f-3877-480f-bd0b-633ce21172e3
md"### Derivation
$x_{n+\Delta n} =  \frac1{N}\sum_{k=0}^{N-1} X_k \exp\left(i2 \pi \frac{k(n+ \Delta n)}{N} \right)$

$=\frac1{N}\sum_{k=0}^{N-1} X_k \exp\left(i2 \pi \frac{kn}{N} \right)\exp\left(i2 \pi \frac{k\Delta n}{N} \right)$

$=\frac1{N}\sum_{k=0}^{N-1} \left( X_k \exp\left(i2 \pi \frac{k\Delta n}{N} \right) \right) \exp\left(i2 \pi \frac{kn}{N} \right)$

$=\mathcal{F}^{-1} \left[ X_k \exp\left(i2 \pi \frac{k\Delta n}{N} \right) \right]$

$f(k, \Delta n) = \exp\left(i2 \pi \frac{k\Delta n}{N} \right)$

"

# ╔═╡ 8bf6ba5b-1d06-41de-b22a-b58c1eaf2c83
"""
	shift(x::AbstractArray{T, N}, Δn::Number)

Shifts a (real) array via the Fourier shift theorem by the amount of
`Δn`.

If `Δn ∈ ℕ` then `circshift(x, -Δn) ≈ shift(x, Δn).`
Otherwise a sub-pixel shift is obtained which is equivalent
to a sinc interpolation.

## Examples
```julia
julia> shift([1.0, 2.0, 3.0, 4.0], 1)
4-element Vector{Float64}:
 2.0
 3.0
 4.0
 1.0

julia> shift([1.0, 2.0, 3.0, 4.0], -1)
4-element Vector{Float64}:
 4.0
 1.0
 2.0
 3.0

julia> shift([1.0, 2.0, 3.0, 4.0], -0.5)
4-element Vector{Float64}:
 2.5
 1.085786437626905
 2.5
 3.914213562373095
```
"""
function shift(x::AbstractArray{T, N3}, Δn::Number) where {T, N3}
	x_ft = ffts(x)

	N = size(x, 1)
	k = (-N÷2) : 1 : (- N÷2 + N - 1)

	x_ft.*= exp.(1im .* 2π .* k .* Δn ./ N)

	return real(iffts(x_ft))

	# alternative solution
	#x_ft = fft(x)
	#f = fftfreq(size(x, 1))
	#phase = exp.(1im * T(2π) .* f .* Δn)
	#return real(ifft(x_ft .* phase))
end

# ╔═╡ 64507581-5751-496e-87b6-88eb320ae44c
N2 = 5

# ╔═╡ 47c0516a-4a37-4823-a447-1207d5d4b8b4
k = (-N2÷2) : 1 : (- N2÷2 + N2 - 1)

# ╔═╡ 250559fb-1098-47ab-a55c-55a89497873d
length(k)

# ╔═╡ a21d8b13-20db-4df1-a2c7-451f8942466b
begin
	x2 = 0:0.1:3
	y2 = sin.(x2.*3) .+ cos.(x2 .* 4).^2
end;

# ╔═╡ 420d1b15-8dc9-4c9a-b91f-9cc4b8101b92
md"
 $\Delta n$=$(@bind Δn Slider(0:0.01:100, show_value=true))"

# ╔═╡ 54ebf4c4-6692-4a40-a83e-c548440b3b85
begin
	plot(shift(y2, Δn), label="FFT based shift", mark="-*")
	plot!(circshift(y2, round(Int, .-Δn)), label="rounded to integer with circshift", mark="-*")
end

# ╔═╡ bbe52a98-9811-4c3a-9b4e-b90a1f5c163f
md"### 3.5 Test"

# ╔═╡ 73fd7360-ffdf-4870-a301-e39a982e0517
arr4 = randn((37,));

# ╔═╡ 37a815a0-c6cd-4ee8-9f46-ee4aa5259690
PlutoTest.@test circshift(arr4, -12) ≈ shift(arr4, 12)

# ╔═╡ 3c09615a-e9c4-462b-9f23-5d5676aa29b7
PlutoTest.@test FourierTools.shift(arr4, -13.1) ≈ shift(arr4, 13.1)

# ╔═╡ d3f2648a-a736-4812-b721-70298151a1a7
md"# 4 Noise Beyond Bandlimit

As seen in the lecture, noise can move information even beyond
the bandlimit of the optical system.
In this task, we want to demonstrate this behaviour with a simulation.
For that we need a few helper functions, try to fill the missing lines.

## 4.1 rr2 Function
In the first step, we create a rr2 function.
"

# ╔═╡ 9ef7592b-92f9-48bb-aa3e-9710bf251379
begin
	"""
		rr2(T, s, center=s .÷ 2 + 1)

	Returns a 2D array which stores the squared distance to the center.

	## Examples
	```julia
	julia> rr2(Float32, (5,5))
	5×5 Matrix{Float32}:
	 8.0  5.0  4.0  5.0  8.0
	 5.0  2.0  1.0  2.0  5.0
	 4.0  1.0  0.0  1.0  4.0
	 5.0  2.0  1.0  2.0  5.0
	 8.0  5.0  4.0  5.0  8.0

	julia> rr2((4,4))
	4×4 Matrix{Float64}:
	 8.0  5.0  4.0  5.0
	 5.0  2.0  1.0  2.0
	 4.0  1.0  0.0  1.0
	 5.0  2.0  1.0  2.0

	julia> rr2((3,3), (3,3))
	3×3 Matrix{Float64}:
	 8.0  5.0  4.0
	 5.0  2.0  1.0
	 4.0  1.0  0.0
	```
	"""
	function rr2(T::DataType, s, center=s .÷ 2 .+ 1)
		is = T.((1:s[1]) .- center[1]).^2
		js = T.(((1:s[2]).- center[2]).^2)'
		return is .+ js
	end

	function rr2(s, center=s .÷ 2 .+ 1)
		rr2(Float64, s, center)
	end
end

# ╔═╡ 9e63206c-82a8-4744-a5a4-a42281077c4f
(1:5).^2 .+ ((1:5).^2)'

# ╔═╡ 6899ad5c-3d31-40d7-94f4-042286b6c1c4
((1:5).^2)'

# ╔═╡ b62519f5-26dd-49ca-afa4-2c72de7ef8a6
rr2((4,4))

# ╔═╡ fc25378b-11b1-4383-b144-75c7dc6e9a1f
rr2((5,5))

# ╔═╡ bd45654d-15e0-4717-962b-d5a302de4a9b
PlutoTest.@test rr2(Float32, (4, 3), (2.5, 3)) == Float32[6.25 3.25 2.25; 4.25 1.25 0.25; 4.25 1.25 0.25; 6.25 3.25 2.25]

# ╔═╡ c6ed3889-0d44-4180-9f81-3095055c6e54
PlutoTest.@test rr2((3, 3), (3, 3)) == [8.0 5.0 4.0; 5.0 2.0 1.0; 4.0 1.0 0.0]

# ╔═╡ 5a16e94e-b01b-4897-9631-41d345229be2
md"## 4.2 circ Function
Having a function which creates the radial distance, we now create a function which returns an array which is 1 inside an circle of radius `r` and 0 outside of it.
"

# ╔═╡ 06a666f9-c99e-41e5-b703-a5e826721d04
begin
	"""
		circ(size, radius)

	`size` the size of the resulting array
	and `radius` the radius of the circle
	"""
	circ(size, radius) = rr2(size) .< radius^2
end

# ╔═╡ 255a1f5a-fe84-4ae2-85c3-a5a512a040ec
md"
if $|a| < |b|$ then also $a^2 < b^2$
for $a,b \in \mathrm{R}^+$
"

# ╔═╡ 9f24aec3-b93c-4e13-a227-0e4deefda611

# ╔═╡ 592e0ed4-bd97-4bb6-a826-b1e2f5ddbd59
circ((5,5), 2)

# ╔═╡ 205a744b-e8a1-4acc-ba6b-7a7307f9a028
md"## 4.2 Test"

# ╔═╡ daf4e474-0db4-4865-859d-96c31e7947f0
PlutoTest.@test circ((10, 2), 2.76) ≈ Bool[0 0; 0 0; 0 0; 1 1; 1 1; 1 1; 1 1; 1 1; 0 0; 0 0]

# ╔═╡ 7fb8fedf-5684-4874-8f17-d01debaff477
PlutoTest.@test circ((4, 5), 2.1) ≈ Bool[0 0 1 0 0; 0 1 1 1 0; 1 1 1 1 1; 0 1 1 1 0]

# ╔═╡ 9b9fb85b-13d0-43a2-ac5f-713593117b96
md" ## 4.3 Frequency Filter
Now we combine all methods together such that we frequency filter
our image.

Procedure:
* Transform image to Fourier space
* take only the inner part of a circle a `radius` of the Fourier space
* go to real space and take real part
"

# ╔═╡ 2a48072b-fc5a-4044-9aa9-5661483ffbfc
function frequency_filter(img, radius)
	img_ft = ffts(img)

	img_ft .*= circ(size(img), radius)

	return real(iffts(img_ft))
end

# ╔═╡ 643bf329-cce9-4c61-931d-65455b59357d
md"r=$(@bind r Slider(1:512, show_value=true))"

# ╔═╡ 2cafea2d-b4b6-4ede-b899-720a71034fda
gray_show(frequency_filter(img, r))

# ╔═╡ 7467608a-ce6d-4f29-8cd8-d766395f7129
md"## 4.4 Noise Breaks Resolution Limit?
We demonstrate now, that noise beyond the frequency limit still contains information

### Procedure
* First frequency filter our image such that it definetely bandlimited
* Apply Poissonian and Gaussian noise
* Now do the exact opposite of the above low-pass filter and take only the high frequency content!
"

# ╔═╡ 6d7fd369-1f6c-487d-b399-c8cbbfede916
function simulate(img, noise_function, radius)
	img_bandlimited = frequency_filter(img, radius)
	img_noisy = noise_function(img_bandlimited)

	img_noisy_ft = ffts(img_noisy)

	# a high pass filter now!
	img_noisy_ft .*= (1 .- circ(size(img), radius))

	return real(iffts(img_noisy_ft))
end

# ╔═╡ e8bd9cb9-bd20-4619-93da-20e7829a61b4
# anonymous noise function we are going to pass to `simulate`
begin
	p(x) = poisson(x, 100000)
	g(x) = add_gauss(x, 0.1)
end

# ╔═╡ 8a479712-7290-4c00-af7c-44e11ec81d0d
md"r=$(@bind r2 Slider(1:300, show_value=true))"

# ╔═╡ 372e190b-0344-4001-9c72-171610b41e9c
img_noisy = abs.(simulate(img, p, r2));

# ╔═╡ b17cd903-bce6-4291-9d71-ea8a8422295a
gray_show(img_noisy, set_one=true)

# ╔═╡ fc91074a-99fc-4a13-802c-740872bf7502
md"
Mean is $(round(mean(img_noisy), sigdigits=3))

Maximum is $(round(maximum(img_noisy), sigdigits=3))

Minimum is $(round(minimum(img_noisy), sigdigits=3))
"

# ╔═╡ fac6c740-93c1-4e12-9610-f3c621f84043
md"## 5.1 Undersampling
Demonstration of aliasing.

Write a function `undersample` which only takes every `n-th` point of an array, keeping the first pixel of the array.
"

# ╔═╡ 15aa846c-7b81-4094-8b5b-dd6575e82afe
img3 = Float32.(testimage("resolution_test_512"));

# ╔═╡ be3b0811-c05a-4e78-9d34-8c52195dd7d0
"""
	undersample(x, factor)

```julia-repl
julia> undersample([1,2,3,4,5,6,7,8,9,10], 2)
5-element Vector{Int64}:
 1
 3
 5
 7
 9

julia> undersample([1,2,3,4,5,6,7,8,9,10], 3)
4-element Vector{Int64}:
  1
  4
  7
 10

julia> undersample([1,2,3,4,5,6,7,8,9,10], 4)
3-element Vector{Int64}:
 1
 5
 9
```
"""
function undersample(x, factor)
	# more Julia idiomatic
	return x[begin:factor:end]
	#return x[1:factor:size(x, 1)]
end

# ╔═╡ 282acd5a-12e8-43f3-8de9-1b7b913a476c
collect(1:5:100)

# ╔═╡ e5188537-b955-4637-bab5-a76b6d74c189
PlutoTest.@test undersample([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 2) == [1, 3, 5, 7, 9]

# ╔═╡ c2b3cdc7-60b5-4287-8268-f95184d452b7
PlutoTest.@test undersample([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 3) == [1, 4, 7, 10]

# ╔═╡ 4e3a272e-0280-4e27-b014-3f9da5ac8039
PlutoTest.@test undersample([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 4) == [1, 5, 9]

# ╔═╡ a76e4081-7c52-4548-b513-ba03954378db

# ╔═╡ f5fbd0da-fe4c-4162-be5e-f9d9d5b27d87
md"## 5.1 Undersampling Visualization

In this plot you can see a nicely oversampled $\sin$ curve.
"

# ╔═╡ a07cb23e-6fca-4c42-a6b6-a525b2b5e232
xs2 = range(0, 10π, length=500)

# ╔═╡ c373046e-04b3-4159-80c0-8fbec141833f
f2(x) = sin(5 * x)

# ╔═╡ 9c1ed9aa-9d52-425d-892d-af3b8a981663
ys2 = f2.(xs2);

# ╔═╡ 75006409-0338-470d-a41d-fa1aba5c4fca
plot(xs2, ys2, xlabel="x pos", ylabel="amplitude")

# ╔═╡ f2b3d8a0-4bea-4443-9bc3-37258106c3d3
md"
Now we only take each `undersample_factor` data point of the series.

`undersample_factor` = $(@bind undersample_factor Slider(1:40, show_value=true))

Analyze what happens once you're below the Nyquist frequency.
See further what happens if you undersample then even more.
"

# ╔═╡ f6b06d6e-b268-45d6-a7e4-d8614cd0884d
y2_undersampled = undersample(ys2, undersample_factor);

# ╔═╡ eab1fe28-f336-4b33-89a3-4c65ed2a818f
plot(fftfreq(length(y2_undersampled), 1 / 10π * 500 / undersample_factor)[1:length(y2_undersampled)÷2],abs.(rfft(y2_undersampled))[1:end-1], xlabel="Positive Frequencies", ylabel="Absolute Value of Frequency amplitude", mark="-*")

# ╔═╡ Cell order:
# ╠═a6261f66-40ca-11ec-3021-ed2173550635
# ╠═77e6389e-3491-45b6-9f76-459b15a8922d
# ╟─7148e705-8b90-46ea-8763-778897daf9a7
# ╟─83d9c252-2a8d-47a3-b247-a62b6a175201
# ╠═39c623c7-6063-4360-a2b7-95ea42c5262d
# ╟─eacd8fdf-971c-46ca-914f-f73f65e28ca8
# ╠═25dde5cb-9ed6-498e-b511-50ed0640d0c0
# ╠═a6c024a4-7cb2-4f8d-be24-8727d75679da
# ╠═d5c127ee-1ecf-4ed0-9cd9-6fcdefa95847
# ╟─6e95cbcb-0d89-4eca-8599-b266076708a7
# ╠═910cb7fa-b1fc-4106-acc4-3bf56edfb83d
# ╠═cbebb311-c521-4fa6-a909-1a4df8da4d2d
# ╠═8bb737cb-2d59-49ba-9b89-afef6ca7bbf8
# ╠═52bc8b3d-542e-4a35-b38c-35840f4075d4
# ╠═e3448ffe-5fbf-4bed-bf03-6eb9031e6c12
# ╠═d9308666-efaa-4fcc-a083-f3d12dbab1e2
# ╠═493b6bb4-1dfc-4a82-9360-40944b29826c
# ╟─c80a2ce8-6620-4de1-bb78-5af97caf466d
# ╠═e2e470e4-ff73-4dd5-973c-14e1535b01ea
# ╠═0a957f4a-a975-4087-b458-3236940eda43
# ╟─af9866f4-2b03-4dc4-9f8a-e02ae46e2f9c
# ╠═ee4b6b2d-119e-4ac3-a1f5-500ff045bfc5
# ╠═7a481af3-4ec8-4a58-bb84-e9d3d5381570
# ╠═4006ccf7-df8e-440c-8b3c-af525bb50ed4
# ╠═c871cdbd-b560-46f2-918c-45d24d2f9b75
# ╠═f6b8cd92-c1ab-445b-8b4f-04f521b72cad
# ╠═51b2d7be-fe94-44e9-8e67-73aafb636ef1
# ╟─a19c8850-9a7d-4e26-908b-a067a9006e84
# ╠═b606e4d0-8cac-4028-8e38-113050531720
# ╟─b0c67bd6-b5ef-4359-8727-1a3bfefcf759
# ╠═7276f9e0-80b8-4794-9125-fca43246f818
# ╠═f054e266-e16c-4311-9e5e-1dbb9d024a63
# ╠═e895cff6-4018-495f-9a0e-1e56ddcfd17d
# ╠═14d68904-4de5-4f1d-81f6-dae8bf50b749
# ╠═7307dd33-cb79-48c0-9e73-78ed70e9c628
# ╠═52061097-5be4-4b1a-a4f6-fd030a5a52ba
# ╠═66cdd72e-ec6f-49ec-8338-c32e744cc359
# ╠═3c77082d-20e9-4541-8426-f8457cc728ec
# ╠═31b60200-4d70-4a3c-995f-949b19800549
# ╟─393d03e2-461d-480a-9719-f33a4eef3a6f
# ╠═df12fd86-f29d-4926-b360-daf4cb40af86
# ╠═37eedc2e-27bd-4d1c-a83b-6ec14aba55ab
# ╟─5313dbed-bfd7-4bb8-b39c-89c55e50d1bb
# ╠═e1eac3cd-e6ce-4743-9119-5ee0c69b8e00
# ╠═481961ac-ad7c-4768-a7aa-eed34143226b
# ╠═25a2de2a-612b-4ed3-8aad-d5caf5d59c30
# ╠═2a74167b-7056-47e1-861b-bf7c1810a1f0
# ╠═c95f9463-de4d-47d4-abfd-f707641d2021
# ╠═6fc3a300-7b96-4b21-8487-fdf4b5eb7521
# ╠═c24b2457-4396-48c6-a21f-6d982d80776e
# ╠═27aa0d27-1cb1-4827-80c7-48741cc0312e
# ╠═fc3bb52b-f9a2-42c0-b34b-347d0a137e2d
# ╠═a55851b6-3f8d-4af3-8f09-a3777c56df1b
# ╟─e46d67be-d974-4f07-8b3e-02fff3e749d0
# ╠═deb9eaf0-ee89-49c5-90b8-570d1a88f8fc
# ╠═90dd14ae-fb42-400e-86b4-9ae9e63716fc
# ╠═34e1c198-abe0-470e-a4bc-348ec09a2032
# ╠═ffb5665d-ecac-4474-b2a4-9685cc40ad7e
# ╠═6ca02109-6114-411f-bfb3-bd638ab93218
# ╟─2875999a-a61d-491a-8a42-0552a87d3c6a
# ╠═333a253e-ce8b-477a-9db2-f99c034ae048
# ╠═e81990c8-2194-472c-a102-b0d010f3a266
# ╠═f3fc2ea5-033c-4789-8ce6-0c6b15b607b7
# ╠═90596b8a-53d4-457f-9e93-522e6c5dc35d
# ╟─4e54ddc6-b9cd-472f-b496-18055ee5b667
# ╠═9844ec03-b169-40a0-ad56-1c07b11bb886
# ╠═49c39cea-895e-429a-909c-d1438d0dfd00
# ╠═fcdb96af-30bd-42c2-bf45-ecd13f27b57d
# ╠═5c4f40df-cde3-4140-ace2-6c0e1b60ebcf
# ╟─79578e4d-7527-459b-8915-ee5cb9ccad24
# ╠═a74c4034-fa19-4360-9f5e-191d5ba4001a
# ╟─8ce249a3-0881-49e5-8ffd-5df21bda50f4
# ╟─13131381-11b3-4db4-9351-fb7264c94ea4
# ╟─bc3ae07f-3877-480f-bd0b-633ce21172e3
# ╠═8bf6ba5b-1d06-41de-b22a-b58c1eaf2c83
# ╠═64507581-5751-496e-87b6-88eb320ae44c
# ╠═47c0516a-4a37-4823-a447-1207d5d4b8b4
# ╠═250559fb-1098-47ab-a55c-55a89497873d
# ╠═a21d8b13-20db-4df1-a2c7-451f8942466b
# ╟─420d1b15-8dc9-4c9a-b91f-9cc4b8101b92
# ╠═54ebf4c4-6692-4a40-a83e-c548440b3b85
# ╟─bbe52a98-9811-4c3a-9b4e-b90a1f5c163f
# ╠═73fd7360-ffdf-4870-a301-e39a982e0517
# ╠═37a815a0-c6cd-4ee8-9f46-ee4aa5259690
# ╠═3c09615a-e9c4-462b-9f23-5d5676aa29b7
# ╟─d3f2648a-a736-4812-b721-70298151a1a7
# ╠═9ef7592b-92f9-48bb-aa3e-9710bf251379
# ╠═9e63206c-82a8-4744-a5a4-a42281077c4f
# ╠═6899ad5c-3d31-40d7-94f4-042286b6c1c4
# ╠═b62519f5-26dd-49ca-afa4-2c72de7ef8a6
# ╠═fc25378b-11b1-4383-b144-75c7dc6e9a1f
# ╠═bd45654d-15e0-4717-962b-d5a302de4a9b
# ╠═c6ed3889-0d44-4180-9f81-3095055c6e54
# ╟─5a16e94e-b01b-4897-9631-41d345229be2
# ╠═06a666f9-c99e-41e5-b703-a5e826721d04
# ╟─255a1f5a-fe84-4ae2-85c3-a5a512a040ec
# ╠═9f24aec3-b93c-4e13-a227-0e4deefda611
# ╠═592e0ed4-bd97-4bb6-a826-b1e2f5ddbd59
# ╟─205a744b-e8a1-4acc-ba6b-7a7307f9a028
# ╠═daf4e474-0db4-4865-859d-96c31e7947f0
# ╠═7fb8fedf-5684-4874-8f17-d01debaff477
# ╟─9b9fb85b-13d0-43a2-ac5f-713593117b96
# ╠═2a48072b-fc5a-4044-9aa9-5661483ffbfc
# ╠═2cafea2d-b4b6-4ede-b899-720a71034fda
# ╠═643bf329-cce9-4c61-931d-65455b59357d
# ╟─7467608a-ce6d-4f29-8cd8-d766395f7129
# ╠═c7150c2c-afeb-4158-b381-3341bb678fd5
# ╠═6d7fd369-1f6c-487d-b399-c8cbbfede916
# ╠═e8bd9cb9-bd20-4619-93da-20e7829a61b4
# ╠═8a479712-7290-4c00-af7c-44e11ec81d0d
# ╠═372e190b-0344-4001-9c72-171610b41e9c
# ╠═b17cd903-bce6-4291-9d71-ea8a8422295a
# ╟─fc91074a-99fc-4a13-802c-740872bf7502
# ╟─fac6c740-93c1-4e12-9610-f3c621f84043
# ╠═15aa846c-7b81-4094-8b5b-dd6575e82afe
# ╠═be3b0811-c05a-4e78-9d34-8c52195dd7d0
# ╠═282acd5a-12e8-43f3-8de9-1b7b913a476c
# ╠═e5188537-b955-4637-bab5-a76b6d74c189
# ╠═c2b3cdc7-60b5-4287-8268-f95184d452b7
# ╠═4e3a272e-0280-4e27-b014-3f9da5ac8039
# ╠═a76e4081-7c52-4548-b513-ba03954378db
# ╠═f5fbd0da-fe4c-4162-be5e-f9d9d5b27d87
# ╠═a07cb23e-6fca-4c42-a6b6-a525b2b5e232
# ╠═c373046e-04b3-4159-80c0-8fbec141833f
# ╠═9c1ed9aa-9d52-425d-892d-af3b8a981663
# ╠═75006409-0338-470d-a41d-fa1aba5c4fca
# ╠═f2b3d8a0-4bea-4443-9bc3-37258106c3d3
# ╠═f6b06d6e-b268-45d6-a7e4-d8614cd0884d
# ╠═eab1fe28-f336-4b33-89a3-4c65ed2a818f


text:
include("run.jl")
include(joinpath(SLIDE_DIR,"src","build","partition.jl"))

# Read and partition the data.
println("  Reading supply/use data...")
io = read_to_check(joinpath(SLIDE_DIR,"data","readfiles","build","partition.yml"));
partition!(io, set)

# Read benchmarking data and benchmark!
println("  Reading i/o data...")
bio_out = read_to_check(joinpath(SLIDE_DIR,"data","readfiles","benchmark","io","io_out.yml"));

io_bench = Dict()
[benchmark!(io_bench, k, bio_out, io) for k in intersect(keys(io), keys(bio_out))];
io_bench

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE